Observation of the decay Bs0 → ηcϕ and evidence for Bs0 → ηcπ+π−

A study of Bs0 → ηcϕ and Bs0 → η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 Bs0 → ηcϕ is reported, where the ηc meson is reconstructed in the pp¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ p\overline{p} $$\end{document}, K+K−π+π−, π+π−π+π− and K+K−K+K− decay modes and the ϕ(1020) in the K+K− decay mode. The decay Bs0 → J/ψϕ is used as a normalisation channel. Evidence is also reported for the decay Bs0 → ηcπ+π−, where the ηc meson is reconstructed in the pp¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ p\overline{p} $$\end{document} decay mode, using the decay Bs0 → J/ψπ+π− as a normalisation channel. The measured branching fractions are ℬBs0→ηcϕ=5.01±0.53±0.27±0.63×10−4,ℬBs0→ηcπ+π−=1.76±0.59±0.12±0.29×10−4,\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \begin{array}{l}\mathrm{\mathcal{B}}\left({B}_s^0\to {\eta}_c\phi \right)=\left(5.01\pm 0.53\pm 0.27\pm 0.63\right)\times {10}^{-4},\hfill \\ {}\mathrm{\mathcal{B}}\left({B}_s^0\to {\eta}_c{\pi}^{+}{\pi}^{-}\right)=\left(1.76\pm 0.59\pm 0.12\pm 0.29\right)\times {10}^{-4},\hfill \end{array} $$\end{document} 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.


Introduction
When a B 0 s meson decays through theb →ccs process, interference between the direct decay amplitude, and the amplitude after B 0 s − B 0 s oscillation, 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][3][4][5]. The φ s phase is best measured using the "golden" channel 1 B 0 s → J/ψ φ [6][7][8][9][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 B 0 s → J/ψ φ, other modes have been used to constrain φ s : B 0 s → J/ψ π + π − [6], B 0 s → D + s D − s [11], and B 0 s → ψ(2S)φ [12]. In this paper, the first study of B 0 s → η c φ and B 0 s → η c π + π − decays is presented. 2 These decays also proceed dominantly through ab →ccs tree diagram as shown in figure 1. Unlike in B 0 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 1 The simplified notation φ and ηc are used to refer to the φ(1020) and the ηc(1S) mesons throughout this article. 2 The use of charge-conjugate modes is implied throughout this article.  Figure 1. Leading diagram corresponding to B 0 s → η c φ and B 0 s → η c π + π − decays, where the π + π − pair may arise from the decay of the f 0 (980) resonance.
analyses of B 0 s → η c φ and B 0 s → η c π + π − decays. Instead, the first measurement of their branching fractions is performed. No prediction is available for either B(B 0 s → η c φ) or B(B 0 s → η c π + π − ). Assuming The measurements presented in this paper are performed using a dataset corresponding to 3 fb −1 of 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 2 describes the LHCb detector and the procedure used to generate simulated events; an overview of the strategy for the measurements of B(B 0 s → η c φ) and B(B 0 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 6 describes the fit models. The results and associated systematic uncertainties are discussed in sections 7 and 8. Finally, conclusions are presented in section 9.

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 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 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/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 [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 B 0 s meson 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.

Analysis strategy
In the analysis of B 0 s → η c φ decays, the φ meson is reconstructed in the K + K − final state and the η c meson is reconstructed in the pp, K + K − π + π − , π + π − π + π − and K + K − K + K − final states. For clarity, the three four-body final states are referred to as 4h throughout the paper. In determining the branching fraction, the decay B 0 s → J/ψ φ is used as a normalisation channel, where the J/ψ meson is reconstructed in the same decay modes as the η c meson. A similar strategy is adopted for the measurement of the branching fraction of B 0 s → η c π + π − decays. However, due to the higher expected level of combinatorial background compared to B 0 s → η c φ decays, the η c and J/ψ mesons are reconstructed only in the pp final state in the measurement of B(B 0 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 B 0 s → η c (→ pp)π + π − decay the fit is done to the ppπ + π − mass distribution, while for the decays B 0 where j stands for the event species, N j is the corresponding yield and N is the vector of yields N j , a is the vector of fitted parameters other than yields, n is the total number -3 -

JHEP07(2017)021
of candidates in the sample, and P j (m) is the probability density function (PDF) used to parametrise the set of invariant-mass distributions m considered. 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 s Plot technique [28].
In the second stage, for B 0 s → ppπ + π − candidates a weighted UML fit is made to the pp invariant-mass spectrum, and weighted UML fits of the pp and the 4h invariant-mass spectra are done for B 0 s → ppφ and B 0 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 where ζ = l ω l / l ω 2 l ensures proper uncertainty estimates from the weighted likelihood fit [29]. For the observed numbers of η c and J/ψ candidates in final state f , N ηc,f and N J/ψ ,f , the measured branching fraction is where X refers to either the φ meson or the π + π − pair. The branching fractions [13], and the efficiency correction factors, ε, are obtained from simulation. In order to maximise the sensitivity to B(B 0 s → η c φ), a simultaneous fit to the pp and 4h invariant-mass spectra is performed.

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 p T or 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 b hadron into charged hadrons. In addition, for the 4h final 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 a b hadron.
In the initial stage of the offline selection, candidates for B 0 s → ppπ + π − and B 0 s → ppK + K − (B 0 s → 4hK + K − ) decays are required to have four (six) good quality, high-p T tracks 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 B 0 s candidates, formed by the combination of the final-state candidates, are required to originate from a PV by requiring a small angle between the B 0 s candidate momentum vector and the vector joining this PV and the B 0 s decay vertex, and a small χ 2 IP , which is defined as the difference in the vertex-fit χ 2 of the considered PV reconstructed with and without the candidate. When forming the B 0 s candidates for B 0 s → ppπ + π − and B 0 s → ppK + K − decays, the pp mass resolution is improved by performing a kinematic fit [33] in which the B 0 s candidate is constrained to originate from its associated PV (that with the smallest value of χ 2 IP for the B 0 s ), and its reconstructed invariant mass is constrained to be equal to the known value of the B 0 s mass [13]. No significant improvement of the 4h mass resolution is observed for B 0 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 B 0 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 B 0 s candidates 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 χ 2 per degree of freedom, the χ 2 IP , and the pointing angle of the B 0 s candidates are used as input to the BDT classifiers used to select candidate B 0 s → ppπ + π − and B 0 s → ppK + K − decays. For the ppK + K − final state, the direction angle, the flight distance significance and the χ 2 IP of the reconstructed B 0 s candidate are also used as input to the BDT, while the p T of the B 0 s candidate is used for the ppπ + π − final state. The difference in the choice of input variables for the ppK + K − and the ppπ + π − 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 B 0 s → ppπ + π − (B 0 s → ppK + 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 B 0 s → 4hK + K − decays using the following set of input variables: the p T and the IP with respect to the SV of all reconstructed tracks; the vertex-fit χ 2 of the η c and φ candidates; the vertex-fit χ 2 , the p T , the flight-distance significance with respect to the PV of the B 0 s candidate, and the angle between the momentum and the vector joining the primary to the secondary vertex of the B 0 s candidate. The optimised requirements on the BDT output and PID variables, for each of the 4h modes, retain about 50% of the signal and reject more than 99% of the combinatorial background inside the mass-fit ranges defined in section 6.
From simulation, after all requirements for B 0 s → 4hK + K − decays, a significant contamination is expected from B 0 s → D + s 3h decays, where the D + s decays 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 4h invariant-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 + s mass [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/c 2 . This requirement is optimised using the significance of B 0 s → J/ψ K + K − candidates with respect to background contributions. This significance is stable for cut values in the range [9, 25] MeV/c 2 , with a maximum at 17 MeV/c 2 , which removes about 90% of B 0 s → D + s 3h decays, with no significant signal loss.

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 B 0 s → η c X and B 0 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 B 0 s decays 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) B 0 s state [35]. The effective lifetime of B 0 s → J/ψ φ (B 0 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. In this section the fit models used for the measurement of the branching fractions are described, first the model used for B 0 s → η c π + π − decays in section 6.1, then the model used for B 0 s → η c φ decays in section 6.2.
6.1 Model for B 0 s → η c π + π − decays Candidates are fitted in two stages. First, an extended UML fit to the ppπ + π − invariantmass spectrum is performed in the range 5150-5540 MeV/c 2 , to discriminate B 0 s → ppπ + π − events from combinatorial background, B 0 → ppπ + π − decays, and B 0 → ppKπ decays, where the kaon is misidentified as a pion. The ppπ + π − mass distribution of B 0 s → ppπ + π − and B 0 → ppπ + π − 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 B 0 → ppKπ 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 ppπ + π − spectrum, signal weights are then computed and the background components are subtracted using the s Plot technique [28]. Correlations between the pp and ppπ + π − invariant-mass spectra, for both signal and backgrounds, are found to be negligible.
Second, a UML fit to the weighted pp invariant-mass distribution is performed in the mass range 2900-3200 MeV/c 2 . In this region, three event categories are expected to populate the pp spectrum: the η c and J/ψ resonances, as well as a possible contribution from nonresonant B 0 s → (pp) NR π + π − decays. The pp mass distribution of η c candidates 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 where m res and Γ res are the pole mass and the natural width, respectively, of the resonance. From simulation, in the mass range considered, the pp invariant-mass resolution is found to be a few MeV/c 2 , while Γ ηc = 31.8 ± 0.8 MeV/c 2 [13]. Thus, the pp distribution of η c candidates is expected to be dominated by the RBW, with only small effects on the total η c lineshape 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 η c and 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, σ cc res , shared between these two functions, is free to vary in the fit to data. The η c pole mass and the mode of the Hypatia function describing the J/ψ lineshape, which can be approximated by the pole mass of the resonance, are also free to vary, while the η c natural width is constrained to its known value [13]. The possible contribution from B 0 s → (pp) 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 (→ pp)π + π − and J/ψ (→ pp)π + π − amplitudes can arise in the pp 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 pp component, interference between the η c (→ pp)π + π − and (pp) NR π + π − amplitudes is neglected.
In order to fully exploit the correlation between the yields of η c and J/ψ candidates, the former is parametrised in the fit, rearranging eq. (3.3), as where B(B 0 s → η c π + π − ) and N J/ψ are free parameters. The yield of the NR pp component is also free to vary.

Model for B 0
s → η c φ decays The procedure and the fit model used to measure B(B 0 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 ppK + 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 ppK + K − versus K + K − and 4hK + K − versus K + K − invariant-mass spectra in the ranges 5200-5500 MeV/c 2 and 990-1050 MeV/c 2 , respectively. 3 Similar 2D fit models are used for each 4h mode. The 4hK + K − distributions of B 0 s → 4hφ signal and B 0 → 4hφ background contributions, as well as those of B 0 s → 4hK + K − and B 0 → 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 JHEP07(2017)021 including a nonresonant K + K − pair are parametrised by linear functions. The expression of the RBW with mass-dependent width describing the φ resonance is the analogue of eq. (6.1), with the mass-dependent width given by where m φ = 1019.461 ± 0.019 MeV/c 2 , Γ φ = 4.266 ± 0.031 MeV/c 2 [13], and q is the magnitude of the momentum of one of the φ decay products, evaluated in the resonance rest frame such that with m K ± = 493.677 ± 0.016 MeV/c 2 [13]. The symbol q φ denotes the value of q when m = m φ . The X(qr) function is the Blatt-Weisskopf barrier factor [38,39] with a barrier radius of r. The value of the parameter r is fixed at 3 (GeV/c) −1 . Defining the quantity z = qr, the Blatt-Weisskopf barrier function for a spin-1 resonance is given by where z φ represents the value of z when m = m φ .
The same 2D fit model is used for the pp mode with an additional component accounting for the presence of misidentified B 0 → ppKπ background events. The ppK + K − and K + K − distributions of B 0 → ppKπ candidates are described by a Crystal Ball function and a linear function, respectively.
Using the sets of signal weights computed from the 2D fits, the pp and 4h spectra are obtained after subtraction of background candidates from B 0 decays and B 0 s decays with nonresonant K + K − pairs as well as combinatorial background. Correlations between the invariant-mass spectra used in the 2D fits and the pp or 4h spectrum are found to be negligible. A simultaneous UML fit is then performed to the weighted pp and 4h invariantmass distributions, with identical mass ranges of 2820-3170 MeV/c 2 . Different models are used to describe the pp and 4h spectra.
The pp invariant-mass spectrum is modelled similarly to the description in section 6.1. However, as shown in section 7, the fit to the pp spectrum for B 0 s → ppπ + π − decays yields a contribution of NR pp decays compatible with zero. Thus, here, the contribution of such decays is fixed to zero and only considered as a source of systematic uncertainty, as described in section 8.
For the 4h modes, in addition to B 0 s → η c φ and B 0 s → J/ψ φ decays, other contributions are expected in the mass range considered: B 0 s → 4hφ decays, where the 4h system is in a nonresonant state with a total angular momentum equal to zero, and where B 0 s decays proceed via intermediate resonant states decaying in turn into two or three particles for instance, B 0 s → P P φ decays, where P and P could be any resonance such as K * (892), ρ(770), φ(1020), ω(782), f 2 (1270), f 2 (1525) and a 2 (1320). Similarly to B 0 s → D + s 3h decays, all these decays are expected to have smooth distributions in the 4h invariant-mass spectra. Therefore, lacking information from previous measurements, all these contributions are -9 -JHEP07(2017)021 merged into one category, denoted (4h) bkg . The 4h nonresonant contribution is denoted (4h) NR . The η c being a pseudoscalar particle, interference between B 0 s → η c (→ 4h)φ and B 0 s → (4h) NR φ amplitudes for each 4h final state are accounted for in the model. On the other hand, given the large number of amplitudes contributing to the (4h) bkg event category, the net effect of all interference terms is assumed to cancel. Similarly to the pp fit model, terms describing residual effects of the interference between the J/ψ and the other fit components are neglected. The total amplitude for each of the 4h modes, integrated over the helicity angles, is then given by Finally, taking into account the detector resolution, the total function, F tot , used to describe the invariant-mass spectra m f is given by with ξ f k = (α f k ) 2 and where the expressions for F k (m f ) are Re R ηc (m f ; a)e iδϕ ⊗ R(a (m f )), (6.12) where δϕ is the difference between the strong phases of (4h) NR φ and η c (→ 4h)φ amplitudes. The integrals in eq. (6.7) are calculated over the mass range in which the fit is performed. Only the η c and J/ψ components are used in the expression for F tot (m pp ). The fit fractions FF k measured for each component, as well as the interference fit fraction FF I between the -10 -

JHEP07(2017)021
η c and the NR amplitudes for the 4h modes, are calculated as: 14) The resolution, R(a (m f )), is described by a Hypatia function, with parameters a (m f ) that depend on the final state and the invariant-mass region. They are estimated using dedicated simulation samples in two mass regions: a high-mass region around the J/ψ resonance, and a low-mass region around the η c resonance. As in the model for B 0 s → ppπ + π − decays, the branching fraction B(B 0 s → η c φ) is directly determined in the fit. In this configuration, the squared magnitudes of the η c amplitudes, ξ f ηc , are parametrised as In the simultaneous fit to the pp and 4h invariant-mass spectra several parameters are allowed to take different values depending on the final state: the intensities ξ f k (free to vary), the slopes κ bkg and κ NR of the (4h) bkg and (4h) NR exponentials, respectively, (free to vary), the relative strong phase between the (4h) NR and η c amplitudes (free to vary) as well as the low and high mass resolution parameters (fixed). The η c pole mass, the mode of the Hypatia function describing the J/ψ and the branching fraction B(B 0 s → η c φ) are common parameters across all final states and are free to vary in the fit. The η c width is fixed to the world average value taken from ref. [13]. For each mode, ξ J/ψ and ϕ ηc are fixed as reference to 1 and 0, respectively.

Results
The yields of the various decay modes determined by the UML fit to the ppπ + π − invariant mass distribution, and from the 2D fits to the pp(4h)K + K − versus K + K − invariant mass planes, are summarised in table 2. The mass distributions and the fit projections are shown in appendix A. The ppπ + π − and 2D fit models are validated using large samples of pseudoexperiments, from which no significant bias is observed.
The pp invariant-mass distribution for B 0 s → ppπ + π − candidates, and the projection of the fit are shown in figure 2. The values of the η c and J/ψ shape parameters as well as the yields are given in table 3. The branching fraction for the B 0 s → η c π + π − decay mode is found to be B(B 0 s → η c π + π − ) = (1.76 ± 0.59 ± 0.12 ± 0.29) × 10 −4 , where the two first uncertainties are statistical and systematic, respectively, and the third uncertainty is due to the limited knowledge of the external branching fractions. The systematic uncertainties on the branching fraction are discussed in section 8. The significance -11 -

JHEP07(2017)021
Yield  Table 2. Yields of the different final states as obtained from the fit to the ppπ + π − invariantmass distribution and from the 2D fits in the pp(4h)K + K − × K + K − invariant-mass planes. Only statistical uncertainties are reported. The abbreviation "n/a" stands for "not applicable". Table 3. Results of the fit to the pp invariant-mass spectra weighted for B 0 s → ppπ + π − candidates. Uncertainties are statistical only. The parameter N NR corresponds to the yield of B 0 s → (pp) NR π + π − candidates. The η c yield does not appear since it is parametrised as a function of B(B 0 s → η c π + π − ), the measured value of which is reported in eq. (7.1).
of the presence of B 0 s → η c π + π − decays in the pp invariant-mass spectrum is estimated, as √ −2∆ ln L, from the difference between the log-likelihood (ln L) values for N ηc = 0 and the value of N ηc that minimises ln L. For the estimation of the significance, N ηc is not parametrised as a function of B(B 0 s → η c π + π − ), but is a free parameter in the fit. As shown in figure 3, the significance of the η c component in the fit to the pp invariant-mass distribution is 5.0 standard deviations (σ) with statistical uncertainties and 4.6σ when including systematic uncertainties. The latter is obtained by adding Gaussian constraints to the likelihood function. This result is the first evidence for B 0 s → η c π + π − decays. The pp and 4h invariant-mass distributions for B 0 s → ppφ and B 0 s → 4hφ candidates, and the projection of the simultaneous fit are shown in figure 4. The values of the shape parameters, of the magnitudes and of the relative strong phases are given in  where the two first uncertainties are statistical and systematic, respectively, and the third uncertainty is due to the limited knowledge of the external branching fractions. This measurement corresponds to the first observation of B 0 s → η c φ decays. As a cross-check, individual fits to the pp and to each of the 4h invariant-mass spectra give compatible values of B(B 0 s → η c φ) within statistical uncertainties. The precision of the B(B 0 s → η c φ) measurement obtained using each of the 4h modes is limited compared to the pp mode. This is expected due to the presence of additional components below the η c and J/ψ resonance in the 4h invariant-mass spectra, and due to the interference between B 0 s → η c (→ 4h)φ and B 0 s → (4h) NR φ amplitudes. The measurement of B(B 0 s → η c φ) from the simultaneous fit is largely dominated by the pp mode.    Data  Table 4. Result of the simultaneous fit to the pp and 4h invariant-mass spectra. Uncertainties are statistical only. The J/ψ and η c magnitudes do not appear since they are set to unity as reference and parametrised as a function of B(B 0 s → η c φ), respectively. In the simultaneous fit, the m ηc and m J/ψ parameters are shared across the four modes. The measured value of B(B 0 s → η c φ) is reported in eq. (7.2). The abbreviation "n/a" stands for "not applicable". 1.08 ± 0.08 1.00 Table 5. Fit fractions obtained from the parameters of the simultaneous fit to the pp and 4h invariant-mass spectra. Uncertainties are statistical only. Due to interference between B 0 s → η c (→ 4h)φ and B 0 s → (4h) NR φ amplitudes, for the 4h final states the sum of fit fractions, k FF k , may be different from unity. The abbreviation "n/a" stands for "not applicable".

Systematic uncertainties
As the expressions for B(B 0 s → η c π + π − ) and B(B 0 s → η c φ) are based on the ratios of observed quantities, only sources of systematic uncertainties inducing different biases to the number of observed η c and 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.  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 r of the φ is arbitrarily set to 3 (GeV/c) −1 . To assign a systematic uncertainty due to the fixed value of this r parameter, 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 4h invariant-mass distributions for acceptance effects. A histogram describing the acceptance effects in each of the 4h invariantmass spectra is constructed from the ratio of the normalised 4h invariant-mass distributions taken from simulated samples of B 0 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 4h invariant-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 pp contribution in the pp spectrum for B 0 s → ppφ 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.

Conclusions
This paper reports the observation of B 0 s → η c φ decays and the first evidence for B 0 s → η c π + π − decays. The branching fractions are measured to be 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 B 0 s → η c π + π − decay mode, including systematic uncertainties, is 4.6σ. The results for B(B 0 s → η c π + π − ) and B(B 0 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 φ s from time-dependent analysis of B 0 s → η c φ or B 0 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:

A Fit projections
The ppπ + π − invariant mass distribution and the fit projection are shown in figure 5. The four pp(4h)K + K − and K + K − invariant-mass distributions and the corresponding twodimensional fit projections are shown in figures 6 to 9.

B Correlation matrix
The statistical correlation matrix for the simultaneous fit to the pp and 4h invariant-mass distributions for B 0 s → ppφ and B 0 s → 4hφ candidates is given in table 7.  Figure 5. Distribution of the ppπ + π − 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 B 0 s → ppπ + π − , B 0 → ppπ + π − , B 0 → ppK + π − and combinatorial background contributions, respectively.      Table 7. Statistical correlation matrix for the parameters from the simultaneous fit to the pp and 4h invariant-mass spectra for B 0 s → ppφ and B 0 s → 4hφ candidates.

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. [12] LHCb collaboration, First study of the CP-violating phase and decay-width difference in B 0 s → ψ(2S)φ decays, Phys. -24 -