Observation of the decay $B_{s}^{0} \to \eta_{c} \phi$ and evidence for $B_{s}^{0} \to \eta_{c} \pi^{+} \pi^{-} $

TL;DR

This paper reports the first observation of the decay $B_{s}^{0} o ext{eta}_c \, \phi$ and evidence for $B_{s}^{0} o ext{eta}_c \, \pi^{+} \, \pi^{-}$ using LHCb data, measuring their branching fractions.

Contribution

It provides the first observation of $B_{s}^{0} o \eta_c \phi$ decay and evidence for $B_{s}^{0} \to \eta_c \pi^{+} \pi^{-}$, with measured branching fractions using Run 1 LHCb data.

Findings

Observation of $B_{s}^{0} \to \eta_c \phi$ decay.

Evidence for $B_{s}^{0} \to \eta_c \pi^{+} \pi^{-}$ decay.

Measured branching fractions for both decays.

Abstract

A study of $B^{0}_{s} \to \eta_{c} \phi$ and $B^{0}_{s} \to \eta_{c} \pi^{+}\pi^{-}$ decays is performed using $pp$ collision data corresponding to an integrated luminosity of 3.0$\,\rm fb^{-1}$, collected with the LHCb detector in Run~1 of the LHC. The observation of the decay $B^{0}_{s} \to \eta_{c} \phi$ is reported, where the $\eta_{c}$ meson is reconstructed in the $p\bar p$, $K^+K^-\pi^+\pi^-$, $\pi^+\pi^-\pi^+\pi^-$ and $K^+K^-K^+K^-$ decay modes and the $\phi(1020)$ in the $K^+ K^-$ decay mode. The decay $B^{0}_{s} \to J/\psi \phi$ is used as a normalisation channel. Evidence is also reported for the decay $B^{0}_{s} \to \eta_{c} \pi^{+}\pi^{-}$, where the $\eta_{c}$ meson is reconstructed in the $p\bar p$ decay mode, using the decay $B^{0}_{s} \to J/\psi \pi^+ \pi^-$ as a normalisation channel. The measured branching fractions are \begin{eqnarray*} {\mathcal B (B^{0}_{s} \to \eta_{c} \phi)} &=& \left(5.01 \pm 0.53 \pm 0.27 \pm 0.63 \right) \times 10^{-4} \,, \nonumber \\ {\mathcal B (B^{0}_{s} \to \eta_{c} \pi^+ \pi^-)} &=& \left(1.76 \pm 0.59 \pm 0.12 \pm 0.29 \right) \times 10^{-4} \,, \end{eqnarray*} 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.

Tables

Table 1. Table 1: Ratio of efficiencies between the normalisation and signal channels for each final state.

	$B_s^0\to 2K2\pi \varphi$	$B_s^0\to 4\pi \varphi$	$B_s^0\to 4K\varphi$	$B_s^0\to p\overline{p}\varphi$	$B_s^0\to p\overline{p}{\pi }^{+}{\pi }^{-}$
$\frac{\epsilon \left(J/\psi \right)}{\epsilon \left({\eta }_c\right)}$	$1.047\pm 0.011$	$1.068\pm 0.016$	$0.962\pm 0.028$	$1.038\pm 0.009$	$1.004\pm 0.015$

Table 2. Table 2: Yields of the different final states as obtained from the fit to the p ​ p ¯ ​ π + ​ π − 𝑝 ¯ 𝑝 superscript 𝜋 superscript 𝜋 p\bar{p}{{\pi}^{+}}{{\pi}^{-}} invariant-mass distribution and from the 2D fits in the p ​ p ¯ ​ ( 4 ​ h ) ​ K + ​ K − × K + ​ K − 𝑝 ¯ 𝑝 4 ℎ superscript 𝐾 superscript 𝐾 superscript 𝐾 superscript 𝐾 p\bar{p}(4h){{K}^{+}}{{K}^{-}}\times{{K}^{+}}{{K}^{-}} invariant-mass planes. Only statistical uncertainties are reported. The abbreviation “n/a” stands for “not applicable”.

	Yield
Mode	$B_s^0\to \mathrm{Mode}$	$B^0\to \mathrm{Mode}$	Combinatorial	$B^0\to p\overline{p}{K}^{+}{\pi }^{-}$
$p\overline{p}{\pi }^{+}{\pi }^{-}$	$179\pm 32$	$384\pm 43$	$3261\pm 119$	$897\pm 69$
$p\overline{p}\varphi$	$447\pm 24$	$13\pm 7$	$43\pm 17$	$11\pm 14$
$p\overline{p}{K}^{+}{K}^{-}$	$10\pm 11$	$-4\pm 5$	$106\pm 19$
$2K2\pi \varphi$	$586\pm 34$	$7\pm 17$	$419\pm 39$	n/a
$2K2\pi K^{+}{K}^{-}$	$86\pm 21$	$18\pm 16$	$329\pm 33$	n/a
$4\pi \varphi$	$502\pm 33$	$77\pm 23$	$380\pm 43$	n/a
$4\pi K^{+}{K}^{-}$	$111\pm 25$	$67\pm 24$	$599\pm 43$	n/a
$4K\varphi$	$151\pm 15$	$6\pm 5$	$44\pm 13$	n/a
$4K{K}^{+}{K}^{-}$	$-3\pm 4$	$-10\pm 9$	$44\pm 11$	n/a

Table 3. Table 3: Results of the fit to the p ​ p ¯ 𝑝 ¯ 𝑝 p\bar{p} invariant-mass spectra weighted for B s 0 → p ​ p ¯ ​ π + ​ π − → subscript superscript 𝐵 0 𝑠 𝑝 ¯ 𝑝 superscript 𝜋 superscript 𝜋 {{B}^{0}_{s}}\rightarrow p\bar{p}{{\pi}^{+}}{{\pi}^{-}} candidates. Uncertainties are statistical only. The parameter N NR subscript 𝑁 NR N_{\rm NR} corresponds to the yield of B s 0 → ( p ​ p ¯ ) NR ​ π + ​ π − → subscript superscript 𝐵 0 𝑠 subscript 𝑝 ¯ 𝑝 NR superscript 𝜋 superscript 𝜋 B^{0}_{s}\rightarrow(p\bar{p})_{\rm NR}\pi^{+}\pi^{-} candidates. The η c subscript 𝜂 𝑐 \eta_{c} yield does not appear since it is parametrised as a function of ℬ ​ ( B s 0 → η c ​ π + ​ π − ) ℬ → subscript superscript 𝐵 0 𝑠 subscript 𝜂 𝑐 superscript 𝜋 superscript 𝜋 {\mathcal{B}}({{B}^{0}_{s}}\rightarrow{\eta_{c}}{{\pi}^{+}}{{\pi}^{-}}) , the measured value of which is reported in Eq. ( 22 ).

$m_{{\eta }_c}$ ($\mathrm{MeV}/c^2$)	$2973\pm 8$
$m_{J/\psi }$ ($\mathrm{MeV}/c^2$)	$3096.9\pm 1.0$
${\sigma }_{\mathrm{res}}^{c\overline{c}}$ ($\mathrm{MeV}/c^2$)	$4.8\pm 0.8$
$N_{J/\psi }$	$113\pm 48$
$N_{\mathrm{NR}}$	$0.5\pm 8.5$

Table 4. Table 4: Result of the simultaneous fit to the p ​ p ¯ 𝑝 ¯ 𝑝 p\bar{p} and 4 ​ h 4 ℎ 4h invariant-mass spectra. Uncertainties are statistical only. The J / ψ 𝐽 𝜓 {J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu} and η c subscript 𝜂 𝑐 \eta_{c} magnitudes do not appear since they are set to unity as reference and parametrised as a function of ℬ ​ ( B s 0 → η c ​ ϕ ) ℬ → subscript superscript 𝐵 0 𝑠 subscript 𝜂 𝑐 italic-ϕ {\mathcal{B}}({{B}^{0}_{s}}\rightarrow{\eta_{c}}\phi) , respectively. In the simultaneous fit, the m η c subscript 𝑚 subscript 𝜂 𝑐 m_{{\eta_{c}}} and m J / ψ subscript 𝑚 𝐽 𝜓 m_{{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}} parameters are shared across the four modes. The measured value of ℬ ​ ( B s 0 → η c ​ ϕ ) ℬ → subscript superscript 𝐵 0 𝑠 subscript 𝜂 𝑐 italic-ϕ {\mathcal{B}}({{B}^{0}_{s}}\rightarrow{\eta_{c}}\phi) is reported in Eq. ( 23 ). The abbreviation “n/a” stands for “not applicable”.

Parameter	Mode
	$2K2\pi$	$4\pi$	$4K$	$p\overline{p}$
$m_{{\eta }_c}$ ($\mathrm{MeV}/c^2$)	$2980.0\pm 2.3$
$m_{J/\psi }$ ($\mathrm{MeV}/c^2$)	$3097.2\pm 0.2$
${\kappa }_{\text{NR}}$ ($\mathrm{MeV}/c^2$)	$-4\pm 4$	$-1\pm 8$	$5\pm 4$	n/a
${\kappa }_{\text{bkg}}$ ($\mathrm{MeV}/c^2$)	-$13\pm 8$	$-14\pm 18$	$4.9\pm 2.9$	n/a
${\xi }_{\text{NR}}$ ($\mathrm{MeV}/c^2$)	-0$0.62\pm 0.29$	-0$0.42\pm 0.31$	$0.5\pm 0.6$	n/a
${\xi }_{\text{bkg}}$ ($\mathrm{MeV}/c^2$)	-0$0.31\pm 0.25$	-0$0.09\pm 0.11$	$1.1\pm 0.7$	n/a
$\delta \phi$ (rad)	-0$1.73\pm 0.18$	-0$2.9\pm 0.6$	$0.3\pm 0.9$	n/a

Table 5. Table 5: Fit fractions obtained from the parameters of the simultaneous fit to the p ​ p ¯ 𝑝 ¯ 𝑝 p\bar{p} and 4 ​ h 4 ℎ 4h invariant-mass spectra. Uncertainties are statistical only. Due to interference between B s 0 → η c ( → 4 h ) ϕ {{B}^{0}_{s}}\rightarrow{\eta_{c}}(\rightarrow 4h)\phi and B s 0 → ( 4 ​ h ) NR ​ ϕ → subscript superscript 𝐵 0 𝑠 subscript 4 ℎ NR italic-ϕ {{B}^{0}_{s}}\rightarrow(4h)_{\rm NR}\phi amplitudes, for the 4 ​ h 4 ℎ 4h final states the sum of fit fractions, ∑ k FF k subscript 𝑘 subscript FF 𝑘 \sum_{k}{\rm FF}_{k} , may be different from unity. The abbreviation “n/a” stands for “not applicable”.

	$2K2\pi$	$4\pi$	$4K$	$p\overline{p}$
FF${}_{J/\psi }$	$0.39\pm 0.03$	$0.28\pm 0.03$	$0.29\pm 0.05$	$0.76\pm 0.03$
FF${}_{{\eta }_c}$	$0.49\pm 0.05$	$0.63\pm 0.07$	$0.31\pm 0.03$	$0.24\pm 0.03$
FF${}_{\text{bkg}}$	$0.12\pm 0.10$	$0.02\pm 0.03$	$0.32\pm 0.19$	n/a
FF${}_{\text{NR}}$	$0.24\pm 0.11$	$0.12\pm 0.09$	$0.15\pm 0.16$	n/a
${\sum }_k{\mathrm{FF}}_k$	$1.24\pm 0.07$	$1.05\pm 0.11$	$1.08\pm 0.08$	$1.00$

Table 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 % percent \% of the measured branching fractions. The abbreviation “n/a” stands for “not applicable”.

Source	Value [%]
	$ℬ\left(B_s^0\to {\eta }_c{\pi }^{+}{\pi }^{-}\right)$	$ℬ\left(B_s^0\to {\eta }_c\varphi \right)$
Fixed PDF parameters	$5.7$	$1.4$
Efficiencies	$3.4$	$0.8$
Fit bias	$1.7$	$1.4$
Resolution model	$0.6$	$4.4$
$\varphi \left(1020\right)$ barrier radius	n/a	$1.6$
Acceptance $\left(4h\right)$	n/a	$1.6$
Nonresonant $p\overline{p}$	n/a	$1.0$
Sum	$6.8$	$5.4$
External branching fractions	$16.4$	$12.6$

Table 7. Table 7: Statistical correlation matrix for the parameters from the simultaneous fit to the p ​ p ¯ 𝑝 ¯ 𝑝 p\bar{p} and 4 ​ h 4 ℎ 4h invariant-mass spectra for B s 0 → p ​ p ¯ ​ ϕ → subscript superscript 𝐵 0 𝑠 𝑝 ¯ 𝑝 italic-ϕ {{B}^{0}_{s}}\rightarrow p\bar{p}\phi and B s 0 → 4 ​ h ​ ϕ → subscript superscript 𝐵 0 𝑠 4 ℎ italic-ϕ {{B}^{0}_{s}}\rightarrow 4h\phi candidates.

	${\kappa }_{\mathrm{NR}}^{2K2\pi }$	${\kappa }_{\mathrm{NR}}^{4K}$	${\kappa }_{\mathrm{NR}}^{4\pi }$	${\kappa }_{\text{bkg}}^{2K2\pi }$	${\kappa }_{\text{bkg}}^{4K}$	${\kappa }_{\text{bkg}}^{4\pi }$	${\xi }_{\mathrm{NR}}^{2K2\pi }$	${\xi }_{\mathrm{NR}}^{4K}$	${\xi }_{\mathrm{NR}}^{4\pi }$
$ℬ\left(B_s^0\to {\eta }_c\varphi \right)$	$+0.22$	$-0.00$	$-0.05$	$-0.19$	$+0.00$	$-0.00$	$+0.55$	$+0.04$	$+0.13$
${\kappa }_{\mathrm{NR}}^{2K2\pi }$		$-0.00$	$-0.01$	$+0.65$	$-0.00$	$-0.00$	$+0.07$	$+0.01$	$+0.02$
${\kappa }_{\mathrm{NR}}^{4K}$			$-0.01$	$+0.00$	$-0.06$	$-0.00$	$+0.00$	$+0.37$	$+0.02$
${\kappa }_{\mathrm{NR}}^{4\pi }$				$-0.00$	$-0.01$	$+0.25$	$-0.10$	$+0.09$	$-0.72$
${\kappa }_{\text{bkg}}^{2K2\pi }$					$+0.00$	$-0.00$	$-0.47$	$-0.01$	$-0.01$
${\kappa }_{\text{bkg}}^{4K}$						$-0.00$	$+0.01$	$-0.21$	$+0.02$
${\kappa }_{\text{bkg}}^{4\pi }$							$-0.01$	$+0.01$	$-0.07$
${\xi }_{\mathrm{NR}}^{2K2\pi }$								$-0.02$	$+0.17$
${\xi }_{\mathrm{NR}}^{4K}$									$-0.12$
		${\xi }_{\mathrm{bkg}}^{2K2\pi }$	${\xi }_{\mathrm{bkg}}^{4K}$	${\xi }_{\mathrm{bkg}}^{4\pi }$	$m_{{\eta }_c}$	$m_{J/\psi }$	$\delta {\phi }_{2K2\pi }$	$\delta {\phi }_{4K}$	$\delta {\phi }_{4\pi }$
$ℬ\left(B_s^0\to {\eta }_c\varphi \right)$		$-0.40$	$-0.04$	$-0.11$	$-0.02$	$+0.00$	$-0.23$	$+0.28$	$-0.44$
${\kappa }_{\mathrm{NR}}^{2K2\pi }$		$-0.00$	$-0.01$	$-0.02$	$+0.01$	$+0.00$	$-0.02$	$+0.06$	$-0.09$
${\kappa }_{\mathrm{NR}}^{4K}$		$-0.00$	$-0.26$	$-0.01$	$-0.03$	$-0.00$	$+0.01$	$-0.10$	$+0.00$
${\kappa }_{\mathrm{NR}}^{4\pi }$		$+0.05$	$-0.05$	$+0.52$	$+0.45$	$-0.04$	$-0.15$	$+0.02$	$+0.09$
${\kappa }_{\text{bkg}}^{2K2\pi }$		$+0.55$	$+0.01$	$+0.01$	$-0.03$	$+0.00$	$+0.49$	$-0.05$	$+0.08$
${\kappa }_{\text{bkg}}^{4K}$		$-0.00$	$+0.04$	$-0.01$	$-0.03$	$-0.00$	$+0.01$	$-0.17$	$-0.00$
${\kappa }_{\text{bkg}}^{4\pi }$		$+0.00$	$-0.00$	$-0.33$	$+0.04$	$-0.00$	$-0.01$	$+0.00$	$+0.03$
${\xi }_{\mathrm{NR}}^{2K2\pi }$		$-0.82$	$-0.00$	$-0.12$	$-0.18$	$+0.01$	$-0.23$	$+0.14$	$-0.24$
${\xi }_{\mathrm{NR}}^{4K}$		$-0.00$	$-0.52$	$+0.07$	$+0.21$	$-0.01$	$-0.08$	$+0.07$	$-0.01$
${\xi }_{\mathrm{NR}}^{4\pi }$		$-0.09$	$+0.06$	$-0.49$	$-0.61$	$+0.01$	$+0.18$	$-0.01$	$+0.05$
${\xi }_{\mathrm{bkg}}^{2K2\pi }$			$+0.01$	$+0.07$	$+0.07$	$-0.02$	$+0.41$	$-0.11$	$+0.18$
${\xi }_{\mathrm{bkg}}^{4K}$				$-0.03$	$-0.11$	$-0.02$	$+0.05$	$-0.26$	$+0.01$
${\xi }_{\mathrm{bkg}}^{4\pi }$					$+0.34$	$-0.02$	$-0.09$	$-0.01$	$+0.28$
$m_{{\eta }_c}$						$-0.01$	$-0.35$	$+0.07$	$+0.02$
$m_{J/\psi }$							$+0.01$	$+0.01$	$-0.01$
$\delta {\phi }_{2K2\pi }$								$-0.09$	$+0.10$
$\delta {\phi }_{4K}$									$-0.12$

Equations

$$\mathcal{B}(B_s^0\to {\eta }_c\varphi )$$

$$\mathcal{B}(B_s^0\to {\eta }_c{\pi }^{+}{\pi }^{-})$$

$$\frac{\mathcal{B}(B_s^0\to {\eta }_c\varphi )}{\mathcal{B}(B_s^0\to J/\psi \varphi )}=\frac{\mathcal{B}(B^0\to {\eta }_c{K}^0)}{\mathcal{B}(B^0\to J/\psi K^0)}=\frac{\mathcal{B}(B_s^0\to {\eta }_c{\pi }^{+}{\pi }^{-})}{\mathcal{B}(B_s^0\to J/\psi {\pi }^{+}{\pi }^{-})}$$

$$\mathcal{B}(B_s^0\to {\eta }_c\varphi )$$

$$\mathcal{B}(B_s^0\to {\eta }_c{\pi }^{+}{\pi }^{-})$$

$$\mathcal{L}(\mathbf{N},\mathbf{a})=\frac{e^{-{\sum }_j{N}_j}}{n!}\prod _{l=1}^n\left(\sum _j{N}_j{\mathcal{P}}_j(m;a)\right),$$

$$-\ln \mathcal{L}(\mathbf{N},\mathbf{a})=\zeta \sum _j{N}_j-\zeta \sum _l{\omega }_l\ln \left(\sum _j{N}_j{\mathcal{P}}_j(m;a)\right)+\ln (n!),$$

$$\mathcal{B}(B_s^0\to {\eta }_{c}X)=\frac{N_{{\eta }_c,f}}{N_{J/\psi ,f}}\times \mathcal{B}(B_s^0\to J/\psi X)\times \frac{\mathcal{B}(J/\psi \to f)}{\mathcal{B}({\eta }_c\to f)}\times \frac{\epsilon (J/\psi {)}_f}{\epsilon ({\eta }_c{)}_f},$$

$$R_{\mathrm{res}}(m;m_{\mathrm{res}},{\Gamma }_{\mathrm{res}})\propto \frac{1}{m_{\mathrm{res}}^2-m^2-i{m}_{\mathrm{res}}{\Gamma }_{\mathrm{res}}},$$

$$N_{{\eta }_c}=N_{J/\psi }\times \frac{\mathcal{B}(B_s^0\to {\eta }_c{\pi }^{+}{\pi }^{-})}{\mathcal{B}(B_s^0\to J/\psi {\pi }^{+}{\pi }^{-})}\times \frac{\mathcal{B}({\eta }_c\to p\bar{p})}{\mathcal{B}(J/\psi \to p\bar{p})}\times \frac{\epsilon ({\eta }_c{)}_{p\bar{p}}}{\epsilon (J/\psi {)}_{p\bar{p}}},$$

$$\Gamma (m)={\Gamma }_{\varphi }{\left(\frac{q}{q_{\varphi }}\right)}^3\left(\frac{m_{\varphi }}{m}\right)X^2(qr),$$

$$q=\frac{1}{2}\sqrt{m^2-4m_{K^{\pm }}^2}.$$

$$X(z)=\sqrt{\frac{1+z_{\varphi }^2}{1+z^2}},$$

$${\left|A(m_f;c_k^f,\mathbf{a})\right|}^2=\sum _k{\left|c_k^f{R}_k(m_f;\mathbf{a})\right|}^2+2\mathcal{R}e(c_{{\eta }_c}^f{R}_{{\eta }_c}(m_f;\mathbf{a})c_{\mathrm{NR}}^{f\ast }{R}_{\mathrm{NR}}^{\ast }(m_f;\mathbf{a})),$$

$$\begin{array}{rl}{\mathcal{F}}_{\mathrm{tot}}(m_f;c_k^f,\mathbf{a},{\mathbf{a}}^{\prime })& ={\left|A(m_f;c_k^f,\mathbf{a})\right|}^2\otimes \mathcal{R}({\mathbf{a}}^{\prime }(m_f))\\ & ={\xi }_{{\eta }_c}^f\frac{{\mathcal{F}}_{{\eta }_c}(m_f)}{{\int }_{m_f}{\mathcal{F}}_{{\eta }_c}(m_f)\mathrm{d}{m}_f}+{\xi }_{J/\psi }^f\frac{{\mathcal{F}}_{J/\psi }(m_f)}{{\int }_{m_f}{\mathcal{F}}_{J/\psi }(m_f)\mathrm{d}{m}_f}\\ & +{\xi }_{\mathrm{NR}}^f\frac{{\mathcal{F}}_{\mathrm{NR}}(m_f)}{{\int }_{m_f}{\mathcal{F}}_{\mathrm{NR}}(m_f)\mathrm{d}{m}_f}+{\xi }_{\mathrm{bkg}}^f\frac{{\mathcal{F}}_{\mathrm{bkg}}(m_f)}{{\int }_{m_f}{\mathcal{F}}_{\mathrm{bkg}}(m_f)\mathrm{d}{m}_f}\\ & +2\sqrt{{\xi }_{{\eta }_c}^f{\xi }_{\mathrm{NR}}^f}\frac{{\mathcal{F}}_{\mathrm{I}}(m_f)}{{\int }_{m_f}\sqrt{{\mathcal{F}}_{{\eta }_c}(m_f){\mathcal{F}}_{\mathrm{NR}}(m_f)}\mathrm{d}{m}_f},\end{array}$$

$${\mathcal{F}}_{{\eta }_c}(m_f)$$

$${\mathcal{F}}_{J/\psi }(m_f)$$

$${\mathcal{F}}_{\mathrm{NR}}(m_f)$$

$${\mathcal{F}}_{\mathrm{bkg}}(m_f)$$

$${\mathcal{F}}_{\mathrm{I}}(m_f)$$

$${\mathrm{FF}}_k^f$$

$${\mathrm{FF}}_{\mathrm{I}}^f$$

$${\xi }_{{\eta }_c}^f={\xi }_{J/\psi }^f\times \frac{\mathcal{B}(B_s^0\to {\eta }_c\varphi )}{\mathcal{B}(B_s^0\to J/\psi \varphi )}\times \frac{\mathcal{B}({\eta }_c\to f)}{\mathcal{B}(J/\psi \to f)}\times \frac{\epsilon ({\eta }_c{)}_f}{\epsilon (J/\psi {)}_f}.$$

$$\mathcal{B}(B_s^0\to {\eta }_c{\pi }^{+}{\pi }^{-})=(1.76\pm 0.59\pm 0.12\pm 0.29)\times 10^{-4},$$

$$\mathcal{B}(B_s^0\to {\eta }_c\varphi )=\left(5.01\pm 0.53\pm 0.27\pm 0.63\right)\times 10^{-4},$$

$$\mathcal{B}(B_s^0\to {\eta }_c\varphi )$$

$$\mathcal{B}(B_s^0\to {\eta }_c{\pi }^{+}{\pi }^{-})$$

Full text

EUROPEAN ORGANIZATION FOR NUCLEAR RESEARCH (CERN)

​​​

CERN-EP-2017-026

LHCb-PAPER-2016-056

Observation of the decay

${{B}^{0}_{s}}\rightarrow{\eta_{c}}\phi$** and evidence**

for ${{B}^{0}_{s}}\rightarrow{\eta_{c}}{{\pi}^{+}}{{\pi}^{-}}$

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

A study of ${{B}^{0}_{s}}\rightarrow{\eta_{c}}\phi$ and ${{B}^{0}_{s}}\rightarrow{\eta_{c}}{{\pi}^{+}}{{\pi}^{-}}$ decays is performed using $pp$ collision data corresponding to an integrated luminosity of 3.0$\mbox{\,fb}^{-1}$, collected with the LHCb detector in Run 1 of the LHC. The observation of the decay ${{B}^{0}_{s}}\rightarrow{\eta_{c}}\phi$ is reported, where the $\eta_{c}$ meson is reconstructed in the $p\bar{p}$, $K^{+}K^{-}\pi^{+}\pi^{-}$, $\pi^{+}\pi^{-}\pi^{+}\pi^{-}$ and $K^{+}K^{-}K^{+}K^{-}$ decay modes and the $\phi(1020)$ in the $K^{+}K^{-}$ decay mode. The decay ${{B}^{0}_{s}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}\phi$ is used as a normalisation channel. Evidence is also reported for the decay ${{B}^{0}_{s}}\rightarrow{\eta_{c}}{{\pi}^{+}}{{\pi}^{-}}$, where the $\eta_{c}$ meson is reconstructed in the $p\bar{p}$ decay mode, using the decay ${{B}^{0}_{s}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}{{\pi}^{+}}{{\pi}^{-}}$ as a normalisation channel. The measured branching fractions are

[TABLE]

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.

Published in JHEP 07 (2017) 021

© CERN on behalf of the LHCb collaboration, licence CC-BY-4.0.

1 Introduction

When a ${B}^{0}_{s}$ meson decays through the $\bar{b}\rightarrow\bar{c}c\bar{s}$ process, interference between the direct decay amplitude, and the amplitude after ${{B}^{0}_{s}}-{{\kern 1.79993pt\overline{\kern-1.79993ptB}{}}{}^{0}_{s}}$ oscillation, gives rise to a $C\!P$-violating phase, $\phi_{{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 $\phi_{{s}}$ phase is best measured using the “golden” channel111The simplified notation $\phi$ and ${\eta_{c}}$ are used to refer to the $\phi(1020)$ and the ${\eta_{c}}(1S)$ mesons throughout this article. ${{B}^{0}_{s}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}\phi$ [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}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}\phi$, other modes have been used to constrain $\phi_{{s}}$: ${{B}^{0}_{s}}\!\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}{{\pi}^{+}}{{\pi}^{-}}$ [6], ${{B}^{0}_{s}}\!\rightarrow{{D}^{+}_{s}}{{D}^{-}_{s}}$ [11], and ${{B}^{0}_{s}}\!\rightarrow\psi({\rm 2S})\phi$ [12].

In this paper, the first study of ${{B}^{0}_{s}}\rightarrow\eta_{c}\phi$ and ${{B}^{0}_{s}}\rightarrow\eta_{c}{{\pi}^{+}}{{\pi}^{-}}$ decays is presented.222The use of charge-conjugate modes is implied throughout this article. These decays also proceed dominantly through a $\bar{b}\rightarrow\bar{c}c\bar{s}$ tree diagram as shown in Fig. 1. Unlike in ${{B}^{0}_{s}}\!\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}\phi$ decays, the ${\eta_{c}}\phi$ final state is purely $C\!P$-even, so that no angular analysis is required to measure the mixing phase $\phi_{s}$. However, the size of the data sample recorded by the LHCb experiment in LHC Run 1 is not sufficient to perform time-dependent analyses of ${{B}^{0}_{s}}\rightarrow\eta_{c}\phi$ and ${{B}^{0}_{s}}\rightarrow\eta_{c}{{\pi}^{+}}{{\pi}^{-}}$ decays. Instead, the first measurement of their branching fractions is performed. No prediction is available for either ${\mathcal{B}}({{B}^{0}_{s}}\rightarrow{\eta_{c}}\phi)$ or ${\mathcal{B}}({{B}^{0}_{s}}\rightarrow{\eta_{c}}{{\pi}^{+}}{{\pi}^{-}})$. Assuming

[TABLE]

allows ${\mathcal{B}}({{B}^{0}_{s}}\rightarrow{\eta_{c}}\phi)$ and ${\mathcal{B}}({{B}^{0}_{s}}\rightarrow{\eta_{c}}{{\pi}^{+}}{{\pi}^{-}})$ to be estimated. From the known values of ${\mathcal{B}}({{B}^{0}}\rightarrow{\eta_{c}}{{K}^{0}})$, ${\mathcal{B}}({{B}^{0}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}{{K}^{0}})$, ${\mathcal{B}}({{B}^{0}_{s}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}\phi)$ and ${\mathcal{B}}({{B}^{0}_{s}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}{{\pi}^{+}}{{\pi}^{-}})$ [13], one finds

[TABLE]

The measurements presented in this paper are performed using a dataset corresponding to 3$\mbox{\,fb}^{-1}$ of integrated luminosity collected by the LHCb experiment in $pp$ collisions during 2011 and 2012 at centre-of-mass energies of 7$\mathrm{\,Te\kern-1.00006ptV}$ and 8$\mathrm{\,Te\kern-1.00006ptV}$, 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 ${\mathcal{B}}({{B}^{0}_{s}}\rightarrow{\eta_{c}}\phi)$ and ${\mathcal{B}}({{B}^{0}_{s}}\rightarrow{\eta_{c}}{{\pi}^{+}}{{\pi}^{-}})$ is given in Sec. 3; the selection of candidate signal decays is described in Sec. 4; the methods to determine the reconstruction and selection efficiencies are discussed in Sec. 5. Section 6 describes the fit models. The results and associated systematic uncertainties are discussed in Secs. 7 and 8. Finally, conclusions are presented in Sec. 9.

2 Detector and simulation

The LHCb detector [14, 15] is a single-arm forward spectrometer covering the pseudorapidity range $2<\eta<5$, designed for the study of particles containing $b$ or $c$ quarks. The detector includes a high-precision tracking system consisting of a silicon-strip vertex detector surrounding the $pp$ interaction region, a large-area silicon-strip detector located upstream of a dipole magnet with a bending power of about $4{\mathrm{\,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${\mathrm{\,Ge\kern-1.00006ptV\!/}c}$. The minimum distance of a track to a primary vertex (PV), the impact parameter (IP), is measured with a resolution of (15+29/\mbox{p_{\mathrm{T}}}){\,\upmu\mathrm{m}}, where $p_{\mathrm{T}}$ is the component of the momentum transverse to the beam, in ${\mathrm{\,Ge\kern-1.00006ptV\!/}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, *Sjostrand:2007gs] 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, *Agostinelli:2002hh] 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 ${{B}^{0}_{s}}\rightarrow{\eta_{c}}\phi$ decays, the $\phi$ meson is reconstructed in the ${{K}^{+}}{{K}^{-}}$ final state and the $\eta_{c}$ meson is reconstructed in the $p\bar{p}$, $K^{+}K^{-}\pi^{+}\pi^{-}$, $\pi^{+}\pi^{-}\pi^{+}\pi^{-}$ 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}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}\phi$ is used as a normalisation channel, where the ${{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}$ meson is reconstructed in the same decay modes as the $\eta_{c}$ meson. A similar strategy is adopted for the measurement of the branching fraction of ${{B}^{0}_{s}}\rightarrow{\eta_{c}}{{\pi}^{+}}{{\pi}^{-}}$ decays. However, due to the higher expected level of combinatorial background compared to ${{B}^{0}_{s}}\rightarrow{\eta_{c}}\phi$ decays, the $\eta_{c}$ and ${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$ mesons are reconstructed only in the $p\bar{p}$ final state in the measurement of ${\mathcal{B}}({{B}^{0}_{s}}\rightarrow{\eta_{c}}{{\pi}^{+}}{{\pi}^{-}})$.

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}}\rightarrow{\eta_{c}}(\rightarrow p\bar{p}){{\pi}^{+}}{{\pi}^{-}}$ decay the fit is done to the $p\bar{p}{{\pi}^{+}}{{\pi}^{-}}$ mass distribution, while for the decays ${{B}^{0}_{s}}\rightarrow{\eta_{c}}(\rightarrow p\bar{p})\phi(\rightarrow{{K}^{+}}{{K}^{-}})$ and ${{B}^{0}_{s}}\rightarrow{\eta_{c}}(\rightarrow 4h)\phi(\rightarrow{{K}^{+}}{{K}^{-}})$ it is made to the two-dimensional $p\bar{p}{{K}^{+}}{{K}^{-}}$ versus ${{K}^{+}}{{K}^{-}}$ or $4h{{K}^{+}}{{K}^{-}}$ versus ${{K}^{+}}{{K}^{-}}$ mass distributions, respectively. The likelihood function is

[TABLE]

where $j$ stands for the event species, $N_{j}$ is the corresponding yield and ${\bf N}$ is the vector of yields $N_{j}$, ${\bf a}$ is the vector of fitted parameters other than yields, $n$ is the total number of candidates in the sample, and ${\cal 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, $\omega_{l}$, are obtained using the ${}_{s}{\cal P}lot$ technique [28].

In the second stage, for ${{B}^{0}_{s}}\rightarrow p\bar{p}{{\pi}^{+}}{{\pi}^{-}}$ candidates a weighted UML fit is made to the $p\bar{p}$ invariant-mass spectrum, and weighted UML fits of the $p\bar{p}$ and the $4h$ invariant-mass spectra are done for ${{B}^{0}_{s}}\rightarrow p\bar{p}\phi$ and ${{B}^{0}_{s}}\rightarrow 4h\phi$ candidates, respectively, to disentangle ${\eta_{c}}$ and ${{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}$ candidates from nonresonant (NR) and remaining background contributions, as described in Sec. 6. For the weighted fits, the NLL function is given by

[TABLE]

where $\zeta=\sum_{l}\omega_{l}/\sum_{l}\omega_{l}^{2}$ ensures proper uncertainty estimates from the weighted likelihood fit [29]. For the observed numbers of $\eta_{c}$ and ${{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}$ candidates in final state $f$, $N_{{\eta_{c}},f}$ and $N_{{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}},f}$, the measured branching fraction is

[TABLE]

where $X$ refers to either the $\phi$ meson or the ${{\pi}^{+}}{{\pi}^{-}}$ pair. The branching fractions $\mathcal{B}({{B}^{0}_{s}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}\phi)$, $\mathcal{B}({{B}^{0}_{s}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}{{\pi}^{+}}{{\pi}^{-}})$, $\mathcal{B}({{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}\rightarrow f)$ and $\mathcal{B}({\eta_{c}}\rightarrow f)$ are taken from Ref. [13], and the efficiency correction factors, $\varepsilon$, are obtained from simulation. In order to maximise the sensitivity to $\mathcal{B}({{B}^{0}_{s}}\rightarrow{\eta_{c}}\phi)$, a simultaneous fit to the $p\bar{p}$ and $4h$ invariant-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 $p_{\mathrm{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 $\phi\rightarrow 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}}\rightarrow p\bar{p}{{\pi}^{+}}{{\pi}^{-}}$ and ${{B}^{0}_{s}}\rightarrow p\bar{p}{{K}^{+}}{{K}^{-}}$ $({{B}^{0}_{s}}\rightarrow 4h{{K}^{+}}{{K}^{-}})$ decays are required to have four (six) good quality, high-$p_{\mathrm{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 $\chi^{2}_{\text{IP}}$, which is defined as the difference in the vertex-fit $\chi^{2}$ of the considered PV reconstructed with and without the candidate. When forming the ${{B}^{0}_{s}}$ candidates for ${{B}^{0}_{s}}\rightarrow p\bar{p}{{\pi}^{+}}{{\pi}^{-}}$ and ${{B}^{0}_{s}}\rightarrow p\bar{p}{{K}^{+}}{{K}^{-}}$ decays, the $p\bar{p}$ 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 $\chi^{2}_{\text{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}}\rightarrow 4h{{K}^{+}}{{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}}\rightarrow 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 $\chi^{2}_{\rm IP}$ of each reconstructed track, as well as the vertex-fit $\chi^{2}$ per degree of freedom, the $\chi^{2}_{\rm 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}}\rightarrow p\bar{p}{{\pi}^{+}}{{\pi}^{-}}$ and ${{B}^{0}_{s}}\rightarrow p\bar{p}{{K}^{+}}{{K}^{-}}$ decays. For the $p\bar{p}K^{+}K^{-}$ final state, the direction angle, the flight distance significance and the $\chi^{2}_{\rm IP}$ of the reconstructed ${B}^{0}_{s}$ candidate are also used as input to the BDT, while the $p_{\mathrm{T}}$ of the ${B}^{0}_{s}$ candidate is used for the $p\bar{p}\pi^{+}\pi^{-}$ final state. The difference in the choice of input variables for the $p\bar{p}K^{+}K^{-}$ and the $p\bar{p}\pi^{+}\pi^{-}$ 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}}\rightarrow p\bar{p}{{\pi}^{+}}{{\pi}^{-}}$ $({{B}^{0}_{s}}\rightarrow p\bar{p}{{K}^{+}}{{K}^{-}})$ decays retain $\sim 45\%$ ($40\%$) of the signal and reject more than $99\%$ ($99\%$) of the combinatorial background, inside the mass-fit ranges defined in Sec. 6.

Dedicated BDT classifiers are trained to select candidate ${{B}^{0}_{s}}\rightarrow 4h{{K}^{+}}{{K}^{-}}$ decays using the following set of input variables: the $p_{\mathrm{T}}$ and the IP with respect to the SV of all reconstructed tracks; the vertex-fit $\chi^{2}$ of the ${\eta_{c}}$ and $\phi$ candidates; the vertex-fit $\chi^{2}$, the $p_{\mathrm{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 Sec. 6.

From simulation, after all requirements for ${{B}^{0}_{s}}\rightarrow 4h{{K}^{+}}{{K}^{-}}$ decays, a significant contamination is expected from ${{B}^{0}_{s}}\rightarrow{{D}^{+}_{s}}3h$ decays, where the ${{D}^{+}_{s}}$ decays to $\phi{{\pi}^{+}}$ 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 $\phi$ candidate and any signal candidate track consistent with a pion hypothesis is required to be $>17{\mathrm{\,Me\kern-1.00006ptV\!/}c^{2}}$. This requirement is optimised using the significance of ${{B}^{0}_{s}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}{{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}}\rightarrow{{D}^{+}_{s}}3h$ decays, with no significant signal loss.

5 Efficiency correction

The efficiency correction factors appearing in Eq. 6 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}}\rightarrow{\eta_{c}}X$ and ${{B}^{0}_{s}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}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 ${\eta_{c}}\phi$ $({\eta_{c}}{{\pi}^{+}}{{\pi}^{-}})$ final state, being purely $C\!P$-even ($C\!P$-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}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}\phi$ (${{B}^{0}_{s}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}{{\pi}^{+}}{{\pi}^{-}}$) 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 Fit models

In this section the fit models used for the measurement of the branching fractions are described, first the model used for ${{B}^{0}_{s}}\rightarrow{\eta_{c}}{{\pi}^{+}}{{\pi}^{-}}$ decays in Sec. 6.1, then the model used for ${{B}^{0}_{s}}\rightarrow{\eta_{c}}\phi$ decays in Sec. 6.2.

6.1 Model for ${{B}^{0}_{s}}\rightarrow{\eta_{c}}{{\pi}^{+}}{{\pi}^{-}}$ decays

Candidates are fitted in two stages. First, an extended UML fit to the $p\bar{p}\pi^{+}\pi^{-}$ invariant-mass spectrum is performed in the range $5150$–5540$${\mathrm{\,Me\kern-1.00006ptV\!/}c^{2}}, to discriminate ${{B}^{0}_{s}}\rightarrow p\bar{p}{{\pi}^{+}}{{\pi}^{-}}$ events from combinatorial background, $B^{0}\rightarrow p\bar{p}{{\pi}^{+}}{{\pi}^{-}}$ decays, and $B^{0}\rightarrow p\bar{p}K\pi$ decays, where the kaon is misidentified as a pion. The $p\bar{p}{{\pi}^{+}}{{\pi}^{-}}$ mass distribution of ${{B}^{0}_{s}}\rightarrow p\bar{p}{{\pi}^{+}}{{\pi}^{-}}$ and $B^{0}\rightarrow p\bar{p}{{\pi}^{+}}{{\pi}^{-}}$ 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}\rightarrow p\bar{p}K\pi$ 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\bar{p}{{\pi}^{+}}{{\pi}^{-}}$ spectrum, signal weights are then computed and the background components are subtracted using the ${}_{s}{\cal P}lot$ technique [28]. Correlations between the $p\bar{p}$ and $p\bar{p}{{\pi}^{+}}{{\pi}^{-}}$ invariant-mass spectra, for both signal and backgrounds, are found to be negligible.

Second, a UML fit to the weighted $p\bar{p}$ invariant-mass distribution is performed in the mass range $2900$–3200$${\mathrm{\,Me\kern-1.00006ptV\!/}c^{2}}. In this region, three event categories are expected to populate the $p\bar{p}$ spectrum: the $\eta_{c}$ and ${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$ resonances, as well as a possible contribution from nonresonant ${{B}^{0}_{s}}\rightarrow(p\bar{p})_{\rm NR}{{\pi}^{+}}{{\pi}^{-}}$ decays. The $p\bar{p}$ mass distribution of $\eta_{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

[TABLE]

where $m_{\rm res}$ and $\Gamma_{\rm res}$ are the pole mass and the natural width, respectively, of the resonance. From simulation, in the mass range considered, the $p\bar{p}$ invariant-mass resolution is found to be a few ${\mathrm{\,Me\kern-1.00006ptV\!/}c^{2}}$, while $\Gamma_{\eta_{c}}=31.8\pm 0.8{\mathrm{\,Me\kern-1.00006ptV\!/}c^{2}}$ [13]. Thus, the $p\bar{p}$ distribution of $\eta_{c}$ candidates is expected to be dominated by the RBW, with only small effects on the total $\eta_{c}$ lineshape from the resolution. On the other hand, due to the small natural width of the ${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$ resonance [13], the corresponding lineshape is assumed to be described to a very good approximation by the resolution function only. For the $\eta_{c}$ and ${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$ 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, $\sigma_{\rm res}^{c\bar{c}}$, shared between these two functions, is free to vary in the fit to data. The $\eta_{c}$ pole mass and the mode of the Hypatia function describing the ${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$ lineshape, which can be approximated by the pole mass of the resonance, are also free to vary, while the $\eta_{c}$ natural width is constrained to its known value [13]. The possible contribution from ${{B}^{0}_{s}}\rightarrow(p\bar{p})_{\rm NR}{{\pi}^{+}}{{\pi}^{-}}$ 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 ${\eta_{c}}(\rightarrow p\bar{p}){{\pi}^{+}}{{\pi}^{-}}$ and ${{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}(\rightarrow p\bar{p}){{\pi}^{+}}{{\pi}^{-}}$ amplitudes can arise in the $p\bar{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\bar{p}$ component, interference between the ${\eta_{c}}(\rightarrow p\bar{p}){{\pi}^{+}}{{\pi}^{-}}$ and $(p\bar{p})_{\rm NR}{{\pi}^{+}}{{\pi}^{-}}$ amplitudes is neglected.

In order to fully exploit the correlation between the yields of $\eta_{c}$ and ${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$ candidates, the former is parametrised in the fit, rearranging Eq. (6), as

[TABLE]

where ${\mathcal{B}}({{B}^{0}_{s}}\rightarrow{\eta_{c}}{{\pi}^{+}}{{\pi}^{-}})$ and $N_{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}$ are free parameters. The yield of the NR $p\bar{p}$ component is also free to vary.

6.2 Model for ${{B}^{0}_{s}}\rightarrow{\eta_{c}}\phi$ decays

The procedure and the fit model used to measure $\mathcal{B}({{B}^{0}_{s}}\rightarrow{\eta_{c}}\phi)$ is based on that described in Sec. 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\bar{p}{{K}^{+}}{{K}^{-}}$ or $4h{{K}^{+}}{{K}^{-}}$ and ${{K}^{+}}{{K}^{-}}$ invariant mass spectra. This allows the contributions from $\phi\rightarrow{{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\bar{p}{{K}^{+}}{{K}^{-}}$ versus ${{K}^{+}}{{K}^{-}}$ and $4h{{K}^{+}}{{K}^{-}}$ versus ${{K}^{+}}{{K}^{-}}$ invariant-mass spectra in the ranges $5200$–5500$${\mathrm{\,Me\kern-1.00006ptV\!/}c^{2}} and $990$–1050$${\mathrm{\,Me\kern-1.00006ptV\!/}c^{2}}, respectively.333In order to better constrain the combinatorial background shape, the upper limit of the $p\bar{p}{{K}^{+}}{{K}^{-}}$ invariant-mass range is extended to 5550$${\mathrm{\,Me\kern-1.00006ptV\!/}c^{2}}.

Similar 2D fit models are used for each $4h$ mode. The $4h{{K}^{+}}{{K}^{-}}$ distributions of ${{B}^{0}_{s}}\rightarrow 4h\phi$ signal and ${{B}^{0}}\rightarrow 4h\phi$ background contributions, as well as those of ${{B}^{0}_{s}}\rightarrow 4h{{K}^{+}}{{K}^{-}}$ and ${{B}^{0}}\rightarrow 4h{{K}^{+}}{{K}^{-}}$ backgrounds, are described by Hypatia functions. The $4h{{K}^{+}}{{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 $\phi$ meson. The $K^{+}K^{-}$ distribution of each contribution including a $\phi$ 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 including a nonresonant ${{K}^{+}}{{K}^{-}}$ pair are parametrised by linear functions. The expression of the RBW with mass-dependent width describing the $\phi$ resonance is the analogue of Eq. (7), with the mass-dependent width given by

[TABLE]

where m_{\phi}=1019.461\pm 0.019$${\mathrm{\,Me\kern-1.00006ptV\!/}c^{2}}, \Gamma_{\phi}=4.266\pm 0.031$${\mathrm{\,Me\kern-1.00006ptV\!/}c^{2}} [13], and $q$ is the magnitude of the momentum of one of the $\phi$ decay products, evaluated in the resonance rest frame such that

[TABLE]

with $m_{K^{\pm}}=493.677\pm 0.016{\mathrm{\,Me\kern-1.00006ptV\!/}c^{2}}$ [13]. The symbol $q_{\phi}$ denotes the value of $q$ when $m=m_{\phi}$. 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\mathrm{\,(GeV/c)^{-1}}$. Defining the quantity $z=qr$, the Blatt-Weisskopf barrier function for a spin-1 resonance is given by

[TABLE]

where $z_{\phi}$ represents the value of $z$ when $m=m_{\phi}$.

The same 2D fit model is used for the $p\bar{p}$ mode with an additional component accounting for the presence of misidentified $B^{0}\rightarrow p\bar{p}K\pi$ background events. The $p\bar{p}{{K}^{+}}{{K}^{-}}$ and ${{K}^{+}}{{K}^{-}}$ distributions of $B^{0}\rightarrow p\bar{p}K\pi$ 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 $p\bar{p}$ 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 $p\bar{p}$ or $4h$ spectrum are found to be negligible. A simultaneous UML fit is then performed to the weighted $p\bar{p}$ and $4h$ invariant-mass distributions, with identical mass ranges of $2820$–3170$${\mathrm{\,Me\kern-1.00006ptV\!/}c^{2}}. Different models are used to describe the $p\bar{p}$ and $4h$ spectra.

The $p\bar{p}$ invariant-mass spectrum is modelled similarly to the description in Sec. 6.1. However, as shown in Sec. 7, the fit to the $p\bar{p}$ spectrum for ${{B}^{0}_{s}}\rightarrow p\bar{p}{{\pi}^{+}}{{\pi}^{-}}$ decays yields a contribution of NR $p\bar{p}$ 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 Sec. 8.

For the $4h$ modes, in addition to ${{B}^{0}_{s}}\rightarrow{\eta_{c}}\phi$ and ${{B}^{0}_{s}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}\phi$ decays, other contributions are expected in the mass range considered: ${{B}^{0}_{s}}\rightarrow 4h\phi$ 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}}\rightarrow PP^{\prime}\phi$ decays, where $P$ and $P^{\prime}$ could be any resonance such as ${{K}^{*}}(892)$, $\rho(770)$, $\phi(1020)$, $\omega(782)$, $f_{2}(1270)$, $f^{\prime}_{2}(1525)$ and $a_{2}(1320)$. Similarly to ${{B}^{0}_{s}}\rightarrow{{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 merged into one category, denoted $(4h)_{\rm bkg}$. The $4h$ nonresonant contribution is denoted $(4h)_{\rm NR}$. The ${\eta_{c}}$ being a pseudoscalar particle, interference between ${{B}^{0}_{s}}\rightarrow{\eta_{c}}(\rightarrow 4h)\phi$ and ${{B}^{0}_{s}}\rightarrow(4h)_{\rm NR}\phi$ 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)_{\rm bkg}$ event category, the net effect of all interference terms is assumed to cancel. Similarly to the $p\bar{p}$ fit model, terms describing residual effects of the interference between the ${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$ 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

[TABLE]

where $R_{k}(m_{f};\mathbf{a})$ is the line-shape of the component $k$, $\mathbf{a}$ represents the line-shape parameters, $c_{k}^{f}$ are complex numbers such that $c^{f}_{k}=\alpha_{k}^{f}\,e^{i\varphi_{k}^{f}}$ where $\alpha_{k}^{f}$ and $\varphi_{k}^{f}$ are the magnitude and the strong phase of amplitude $k$, and $m_{f}$ is one of the $4h$ invariant masses. The ${\eta_{c}}$ and the ${{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}$ resonances are described similarly to the $p\bar{p}$ mode, and the $(4h)_{\rm NR}$ and $(4h)_{\rm bkg}$ components are described using exponential functions.

Finally, taking into account the detector resolution, the total function, ${\mathcal{F}}_{\rm tot}$, used to describe the invariant-mass spectra $m_{f}$ is given by

[TABLE]

with $\xi^{f}_{k}=(\alpha_{k}^{f})^{2}$ and where the expressions for ${\mathcal{F}}_{k}(m_{f})$ are

[TABLE]

where $\delta\varphi$ is the difference between the strong phases of $(4h)_{\rm NR}\phi$ and ${\eta_{c}}(\rightarrow 4h)\phi$ amplitudes. The integrals in Eq. (13) are calculated over the mass range in which the fit is performed. Only the $\eta_{c}$ and ${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$ components are used in the expression for ${\mathcal{F}}_{\rm tot}(m_{p\bar{p}})$. The fit fractions FFk measured for each component, as well as the interference fit fraction FFI between the ${\eta_{c}}$ and the NR amplitudes for the $4h$ modes, are calculated as:

[TABLE]

The resolution, $\mathcal{R}(\mathbf{a}^{\prime}(m_{f}))$, is described by a Hypatia function, with parameters $\mathbf{a}^{\prime}(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\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$ resonance, and a low-mass region around the $\eta_{c}$ resonance.

As in the model for ${{B}^{0}_{s}}\rightarrow p\bar{p}{{\pi}^{+}}{{\pi}^{-}}$ decays, the branching fraction $\mathcal{B}({{B}^{0}_{s}}\rightarrow{\eta_{c}}\phi)$ is directly determined in the fit. In this configuration, the squared magnitudes of the $\eta_{c}$ amplitudes, $\xi^{f}_{\eta_{c}}$, are parametrised as

[TABLE]

In the simultaneous fit to the $p\bar{p}$ and $4h$ invariant-mass spectra several parameters are allowed to take different values depending on the final state: the intensities $\xi^{f}_{k}$ (free to vary), the slopes $\kappa_{\rm bkg}$ and $\kappa_{\rm NR}$ of the $(4h)_{\rm bkg}$ and $(4h)_{\rm NR}$ exponentials, respectively, (free to vary), the relative strong phase between the $(4h)_{\rm NR}$ and $\eta_{c}$ amplitudes (free to vary) as well as the low and high mass resolution parameters (fixed). The $\eta_{c}$ pole mass, the mode of the Hypatia function describing the ${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$ and the branching fraction $\mathcal{B}({{B}^{0}_{s}}\rightarrow{\eta_{c}}\phi)$ are common parameters across all final states and are free to vary in the fit. The $\eta_{c}$ width is fixed to the world average value taken from Ref. [13]. For each mode, $\xi_{{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}}$ and $\varphi_{\eta_{c}}$ are fixed as reference to 1 and 0, respectively.

7 Results

The yields of the various decay modes determined by the UML fit to the $p\bar{p}{{\pi}^{+}}{{\pi}^{-}}$ invariant mass distribution, and from the 2D fits to the $p\bar{p}(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 $p\bar{p}{{\pi}^{+}}{{\pi}^{-}}$ and 2D fit models are validated using large samples of pseudoexperiments, from which no significant bias is observed.

The $p\bar{p}$ invariant-mass distribution for ${{B}^{0}_{s}}\rightarrow p\bar{p}{{\pi}^{+}}{{\pi}^{-}}$ candidates, and the projection of the fit are shown in Fig. 2. The values of the $\eta_{c}$ and ${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$ shape parameters as well as the yields are given in Table 3. The branching fraction for the ${{B}^{0}_{s}}\rightarrow{\eta_{c}}{{\pi}^{+}}{{\pi}^{-}}$ decay mode is found to be

[TABLE]

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 Sec. 8. The significance of the presence of ${{B}^{0}_{s}}\rightarrow{\eta_{c}}{{\pi}^{+}}{{\pi}^{-}}$ decays in the $p\bar{p}$ invariant-mass spectrum is estimated, as $\sqrt{-2\Delta\ln L}$, from the difference between the log-likelihood ($\ln L$) values for $N_{\eta_{c}}=0$ and the value of $N_{\eta_{c}}$ that minimises $\ln L$. For the estimation of the significance, $N_{\eta_{c}}$ is not parametrised as a function of ${\mathcal{B}}({{B}^{0}_{s}}\rightarrow{\eta_{c}}{{\pi}^{+}}{{\pi}^{-}})$, but is a free parameter in the fit. As shown in Fig. 3, the significance of the $\eta_{c}$ component in the fit to the $p\bar{p}$ invariant-mass distribution is $5.0$ standard deviations ($\sigma$) with statistical uncertainties and $4.6\sigma$ 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}}\rightarrow{\eta_{c}}{{\pi}^{+}}{{\pi}^{-}}$ decays.

The $p\bar{p}$ and $4h$ invariant-mass distributions for ${{B}^{0}_{s}}\rightarrow p\bar{p}\phi$ and ${{B}^{0}_{s}}\rightarrow 4h\phi$ candidates, and the projection of the simultaneous fit are shown in Fig. 4. The values of the shape parameters, of the magnitudes and of the relative strong phases are given in Table 4. The statistical correlation matrix of the simultaneous fit is given in Appendix B. The fit fractions are given in Table 5. The measured branching fraction for the ${{B}^{0}_{s}}\rightarrow{\eta_{c}}\phi$ decay mode is

[TABLE]

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}}\rightarrow{\eta_{c}}\phi$ decays. As a cross-check, individual fits to the $p\bar{p}$ and to each of the $4h$ invariant-mass spectra give compatible values of $\mathcal{B}(B^{0}_{s}\rightarrow\eta_{c}\phi)$ within statistical uncertainties. The precision of the $\mathcal{B}(B^{0}_{s}\rightarrow\eta_{c}\phi)$ measurement obtained using each of the $4h$ modes is limited compared to the $p\bar{p}$ mode. This is expected due to the presence of additional components below the $\eta_{c}$ and ${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$ resonance in the $4h$ invariant-mass spectra, and due to the interference between ${{B}^{0}_{s}}\rightarrow{\eta_{c}}(\rightarrow 4h)\phi$ and ${{B}^{0}_{s}}\rightarrow(4h)_{\rm NR}\phi$ amplitudes. The measurement of $\mathcal{B}(B^{0}_{s}\rightarrow\eta_{c}\phi)$ from the simultaneous fit is largely dominated by the $p\bar{p}$ mode.

8 Systematic uncertainties

As the expressions for $\mathcal{B}({{B}^{0}_{s}}\rightarrow{\eta_{c}}{{\pi}^{+}}{{\pi}^{-}})$ and $\mathcal{B}({{B}^{0}_{s}}\rightarrow{\eta_{c}}\phi)$ are based on the ratios of observed quantities, only sources of systematic uncertainties inducing different biases to the number of observed $\eta_{c}$ and ${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$ 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.

The Blatt-Weisskopf parameter $r$ of the $\phi$ is arbitrarily set to $3\mathrm{\,({\mathrm{\,Ge\kern-1.00006ptV\!/}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\mathrm{\,({\mathrm{\,Ge\kern-1.00006ptV\!/}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$ invariant-mass spectra is constructed from the ratio of the normalised $4h$ invariant-mass distributions taken from simulated samples of ${{B}^{0}_{s}}\rightarrow(4h)\phi$ 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 $p\bar{p}$ contribution in the $p\bar{p}$ spectrum for ${{B}^{0}_{s}}\rightarrow p\bar{p}\phi$ 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 ${{B}^{0}_{s}}\rightarrow{\eta_{c}}\phi$ decays and the first evidence for ${{B}^{0}_{s}}\rightarrow{\eta_{c}}{{\pi}^{+}}{{\pi}^{-}}$ decays. The branching fractions are measured to be

[TABLE]

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}}\rightarrow{\eta_{c}}{{\pi}^{+}}{{\pi}^{-}}$ decay mode, including systematic uncertainties, is $4.6\sigma$. The results for ${\mathcal{B}}({{B}^{0}_{s}}\rightarrow{\eta_{c}}{{\pi}^{+}}{{\pi}^{-}})$ and ${\mathcal{B}}({{B}^{0}_{s}}\rightarrow{\eta_{c}}\phi)$ are in agreement with expectations based on Eqs. (1), (2) and (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 $C\!P$-violating phase $\phi_{{s}}$ from time-dependent analysis of ${{B}^{0}_{s}}\rightarrow{\eta_{c}}\phi$ or ${{B}^{0}_{s}}\rightarrow{\eta_{c}}{{\pi}^{+}}{{\pi}^{-}}$ 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 $\phi_{{s}}$.

Acknowledgements

We express our gratitude to our colleagues in the CERN accelerator departments for the excellent performance of the LHC. We thank the technical and administrative staff at the LHCb institutes. We acknowledge support from CERN and from the national agencies: CAPES, CNPq, FAPERJ and FINEP (Brazil); MOST and NSFC (China); CNRS/IN2P3 (France); BMBF, DFG and MPG (Germany); INFN (Italy); 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 (USA). 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 (USA). We are indebted to the communities behind the multiple open source software packages on which we depend. Individual groups or members have received support from AvH Foundation (Germany), EPLANET, Marie Skłodowska-Curie Actions and ERC (European Union), Conseil Général de Haute-Savoie, Labex ENIGMASS and OCEVU, Région 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).

Appendix

Appendix A Fit projections

The $p\bar{p}{{\pi}^{+}}{{\pi}^{-}}$ invariant mass distribution and the fit projection are shown in Fig. 5. The four $p\bar{p}(4h){{K}^{+}}{{K}^{-}}$ and ${{K}^{+}}{{K}^{-}}$ invariant-mass distributions and the corresponding two-dimensional fit projections are shown in Figs. 6 to 9.

Appendix B Correlation matrix

The statistical correlation matrix for the simultaneous fit to the $p\bar{p}$ and $4h$ invariant-mass distributions for ${{B}^{0}_{s}}\rightarrow p\bar{p}\phi$ and ${{B}^{0}_{s}}\rightarrow 4h\phi$ candidates is given in Table 7.