Evidence for $B^+ \rightarrow h_c K^+$ and observation of $\eta_c(2S) \to p \bar{p} \pi^+ \pi^-$

TL;DR

This paper reports evidence for the decay $B^+ ightarrow h_c K^+$, measures its branching fraction, sets an upper limit for $B^0 ightarrow h_c K_S^0$, and observes the decay $ ext{eta}_c(2S) o p ar{p} \pi^+ \pi^-$ with high significance.

Contribution

The study provides the first evidence for $B^+ ightarrow h_c K^+$ decay and the first observation of $ ext{eta}_c(2S) o p ar{p} \pi^+ \pi^-$ decay, expanding knowledge of charmonium-related B decays.

Findings

Evidence for $B^+ ightarrow h_c K^+$ with 4.8σ significance.

Measured branching fraction of $(3.7^{+1.0}_{-0.9}{}^{+0.8}_{-0.8}) imes 10^{-5}$.

First observation of $ ext{eta}_c(2S) o p ar{p} \pi^+ \pi^-$ with 12.1σ significance.

Abstract

A search for the decays $B^+ \rightarrow h_c K^+$ and $B^0 \rightarrow h_c K_S^0$ is performed. Evidence for the decay $B^+ \rightarrow h_c K^+$ is found; its significance is $4.8\sigma$. No evidence is found for $B^0 \rightarrow h_c K_S^0$. The branching fraction for $B^+ \rightarrow h_c K^+$ is measured to be $(3.7^{+1.0}_{-0.9}{}^{+0.8}_{-0.8}) \times 10^{-5}$; the upper limit for the $B^0 \rightarrow h_c K_S^0$ branching fraction is $1.4 \times 10^{-5}$ at $90\%$ C.L. In addition, a study of the $p \bar{p} \pi^+ \pi^-$ invariant mass distribution in the channel $B^+ \to (p \bar{p} \pi^+ \pi^-) K^+$ results in the first observation of the decay $\eta_c(2S) \to p \bar{p} \pi^+ \pi^-$ with $12.1\sigma$ significance. The analysis is based on the 711 $\mathrm{fb}^{-1}$ data sample collected by the Belle detector at the asymmetric-energy $e^+ e^-$ collider KEKB at the $\Upsilon(4S)$ resonance.

Tables

Table 1. Table 1: Results of the optimization of the selection requirements. The signal-region half-axes R Δ ​ E ( i ) superscript subscript 𝑅 Δ 𝐸 𝑖 R_{\Delta E}^{(i)} ( R M bc ( i ) superscript subscript 𝑅 subscript 𝑀 bc 𝑖 R_{M_{\text{bc}}}^{(i)} ) are in MeV ( MeV / c 2 MeV superscript 𝑐 2 \text{MeV}/c^{2} ); all other values are dimensionless.

Channel	$B^{+}\to h_c{K}^{+}$						$B^0\to h_c{K}_S^0$
	Parameters			Efficiency			Parameters			Efficiency
	$R_{\Delta E}^{\left(i\right)}$	$R_{M_{\text{bc}}}^{\left(i\right)}$	$v_0^{\left(i\right)}$	${ϵ}_{\text{SR}}^{\left(i\right)}$	${ϵ}_S^{\left(i\right)}\left(v_0^{\left(i\right)}\right)$	${ϵ}_B^{\left(i\right)}\left(v_0^{\left(i\right)}\right)$	$R_{\Delta E}^{\left(i\right)}$	$R_{M_{\text{bc}}}^{\left(i\right)}$	$v_0^{\left(i\right)}$	${ϵ}_{\text{SR}}^{\left(i\right)}$	${ϵ}_S^{\left(i\right)}\left(v_0^{\left(i\right)}\right)$	${ϵ}_B^{\left(i\right)}\left(v_0^{\left(i\right)}\right)$
Channel group 1: $h_c\to {\eta }_c\gamma$
${\eta }_c(\to K^{+}{K}_S^0{\pi }^{-})\gamma$	32.7	4.82	0.804	6.27%	59.5%	5.08%	34.2	4.97	0.702	4.31%	69.8%	8.55%
${\eta }_c(\to K^{+}{K}^{-}{\pi }^0)\gamma$	36.2	3.90	0.958	4.27%	31.7%	0.56%	43.5	4.54	0.942	3.38%	39.3%	0.85%
${\eta }_c(\to K_S^0{K}_S^0{\pi }^0)\gamma$	42.3	4.49	0.976	1.79%	17.8%	0.18%	35.9	4.07	0.954	1.05%	35.9%	0.64%
${\eta }_c(\to K^{+}{K}^{-}{\eta }_{2\gamma })\gamma$	34.4	4.16	0.977	4.21%	20.2%	0.22%	37.8	4.49	0.967	3.10%	28.0%	0.41%
${\eta }_c(\to K^{+}{K}^{-}{\eta }_{3\pi })\gamma$	24.9	3.59	0.978	1.75%	29.7%	0.23%	33.2	4.55	0.986	1.50%	23.7%	0.17%
${\eta }_c(\to K^{+}{K}^{-}{K}^{+}{K}^{-})\gamma$	25.3	4.13	0.770	4.89%	53.2%	6.69%	29.9	4.80	0.734	3.71%	56.9%	8.49%
${\eta }_c(\to {\eta }^{\prime }(\to {\eta }_{2\gamma }{\pi }^{+}{\pi }^{-}){\pi }^{+}{\pi }^{-})\gamma$	30.5	4.21	0.958	2.69%	40.5%	0.63%	32.0	4.50	0.946	1.87%	45.6%	0.96%
${\eta }_c(\to {\eta }^{\prime }(\to {\eta }_{3\pi }{\pi }^{+}{\pi }^{-}){\pi }^{+}{\pi }^{-})\gamma$	26.8	4.16	0.990	1.01%	29.2%	0.13%	24.6	3.87	0.986	0.59%	32.3%	0.26%
${\eta }_c(\to p\overline{p})\gamma$	38.9	5.48	0.654	17.70%	75.7%	10.66%	42.5	5.85	0.513	12.47%	82.8%	15.43%
${\eta }_c(\to p\overline{p}{\pi }^0)\gamma$	30.2	3.75	0.954	4.65%	30.3%	0.50%	31.8	4.01	0.934	3.33%	39.4%	0.93%
${\eta }_c(\to p\overline{p}{\pi }^{+}{\pi }^{-})\gamma$	24.1	4.03	0.912	6.31%	30.0%	1.53%	24.3	4.09	0.860	4.23%	41.8%	3.26%
${\eta }_c(\to \Lambda \overline{\Lambda })\gamma$	40.4	5.66	0.727	4.04%	70.6%	6.79%	41.5	5.19	0.586	2.65%	76.3%	11.24%
Channel group 2: $h_c\to p\overline{p}{\pi }^{+}{\pi }^{-}$
$p\overline{p}{\pi }^{+}{\pi }^{-}$	13.5	4.36	0.598	14.81%	64.6%	18.40%	13.8	4.56	0.519	10.30%	71.2%	24.20%

Table 2. Table 2: The resulting signal yields and phases in the default model. The errors are statistical only.

Parameter	$B^{+}\to h_c{K}^{+}$	$B^0\to h_c{K}_S^0$
$N_{{\eta }_c}$	$229\pm 18$	$96\pm 11$
${\phi }_{{\eta }_c}$	$-1.12\pm 0.09$	$-1.47\pm 0.14$
$N_{J/\psi }$	$345\pm 19$	$128\pm 12$
$N_{{\chi }_{c0}}$	$25.5\pm 7.1$	$0.9\pm 1.6$
${\phi }_{{\chi }_{c0}}$	$-1.51\pm 0.23$	$-2.05\pm 1.20$
$N_{{\chi }_{c1}}$	$34.5\pm 8.8$	$21.2\pm 6.1$
$N_{h_c}$	$32.6\pm 8.0$	$3.1\pm 3.8$
$N_{{\chi }_{c2}}$	$-1.6\pm 6.3$	$11.6\pm 5.5$
$N_{{\eta }_c\left(2S\right)}$	$86.1\pm 11.9$	$24.0\pm 6.8$
${\phi }_{{\eta }_c\left(2S\right)}$	$-1.41\pm 0.14$	$-1.90\pm 0.23$
$N_{\psi \left(2S\right)}$	$36.9\pm 8.8$	$13.0\pm 5.6$

Table 3. Table 3: Significances of the charmonium states decaying to p ​ p ¯ ​ π + ​ π − 𝑝 ¯ 𝑝 superscript 𝜋 superscript 𝜋 p\bar{p}\pi^{+}\pi^{-} except the h c subscript ℎ 𝑐 h_{c} in the default model.

State	$B^{+}\to \left(c\overline{c}\right)K^{+}$	$B^0\to \left(c\overline{c}\right)K_S^0$
${\eta }_c$	$20.1\sigma$	$12.5\sigma$
$J/\psi$	$33.9\sigma$	$20.8\sigma$
${\chi }_{c0}$	$6.0\sigma$	$0.6\sigma$
${\chi }_{c1}$	$4.9\sigma$	$4.5\sigma$
${\chi }_{c2}$	$0.3\sigma$	$2.5\sigma$
${\eta }_c\left(2S\right)$	$12.3\sigma$	$5.9\sigma$
$\psi \left(2S\right)$	$5.0\sigma$	$2.8\sigma$

Table 4. Table 4: Model dependence of the h c subscript ℎ 𝑐 h_{c} and η c ​ ( 2 ​ S ) → p ​ p ¯ ​ π + ​ π − → subscript 𝜂 𝑐 2 𝑆 𝑝 ¯ 𝑝 superscript 𝜋 superscript 𝜋 \eta_{c}(2S)\to p\bar{p}\pi^{+}\pi^{-} significance.

Model	$h_c$ significance		${\eta }_c\left(2S\right)\to p\overline{p}{\pi }^{+}{\pi }^{-}$ significance
	$B^{+}\to h_c{K}^{+}$	$B^0\to h_c{K}_S^0$	$B^{+}\to h_c{K}^{+}$	$B^0\to h_c{K}_S^0$
Default	$5.0\sigma$	$0.8\sigma$	$12.3\sigma$	$5.9\sigma$
Free masses and widths	$5.0\sigma$	$0.8\sigma$	$12.3\sigma$	$6.0\sigma$
Polynomial order ($h_c\to {\eta }_c\gamma$)	$4.8\sigma$	$0.8\sigma$	$12.3\sigma$	$5.9\sigma$
Polynomial order ($h_c\to p\overline{p}{\pi }^{+}{\pi }^{-}$ background)	$5.0\sigma$	$0.8\sigma$	$12.2\sigma$	$5.9\sigma$
Polynomial order ($h_c\to p\overline{p}{\pi }^{+}{\pi }^{-}$ signal)	$5.0\sigma$	$0.9\sigma$	$12.2\sigma$	$5.9\sigma$
Fitting range variation ($h_c\to {\eta }_c\gamma$)	$5.0\sigma$	$0.9\sigma$	$12.3\sigma$	$5.9\sigma$
Fitting range variation ($h_c\to p\overline{p}{\pi }^{+}{\pi }^{-}$)	$5.0\sigma$	$0.8\sigma$	$12.1\sigma$	$5.8\sigma$
Scaled resolution	$5.0\sigma$	$0.8\sigma$	$12.3\sigma$	$6.0\sigma$
Fraction of $h_c\to {\eta }_c\gamma$ and $h_c\to p\overline{p}{\pi }^{+}{\pi }^{-}$	$4.9\sigma$	$0.7\sigma$	$12.2\sigma$	$5.9\sigma$

Table 5. Table 5: Relative systematic uncertainties of the branching fractions for the channel B + → ( c ​ c ¯ ) ​ K + → superscript 𝐵 𝑐 ¯ 𝑐 superscript 𝐾 B^{+}\to(c\bar{c})K^{+} ( B + → h c ​ K + → superscript 𝐵 subscript ℎ 𝑐 superscript 𝐾 B^{+}\to h_{c}K^{+} for the h c subscript ℎ 𝑐 h_{c} and B + → ( c c ¯ ) ( → p p ¯ π + π − ) K + B^{+}\to(c\bar{c})(\to p\bar{p}\pi^{+}\pi^{-})K^{+} for all other charmonium states).

Error source	$h_c$	${\eta }_c$	$J/\psi$	${\chi }_{c0}$	${\chi }_{c1}$	${\chi }_{c2}$	${\eta }_c\left(2S\right)$	$\psi \left(2S\right)$
Model dependence	$\left(\stackrel{⃡}{\ast }{[}_{-9.16}^{+9.09}]\right)\%$	$\left(\stackrel{⃡}{\ast }{[}_{-1.47}^{+3.45}]\right)\%$	$\left(\stackrel{⃡}{\ast }{[}_{-0.04}^{+1.96}]\right)\%$	$\left(\stackrel{⃡}{\ast }{[}_{-5.32}^{+4.59}]\right)\%$	$\left(\stackrel{⃡}{\ast }{[}_{-2.98}^{+7.10}]\right)\%$	$\left(\stackrel{⃡}{\ast }{[}_{-65.93}^{+21.13}]\right)\%$	$\left(\stackrel{⃡}{\ast }{[}_{-4.62}^{+1.75}]\right)\%$	$\left(\stackrel{⃡}{\ast }{[}_{-7.07}^{+1.74}]\right)\%$
PID	3.99%	3.64%	3.62%	3.50%	3.50%	3.51%	3.52%	3.53%
Overtraining	0.41%	0.14%
Tracking	1.60%	1.75%
MLP efficiency	12.73%	0.25%
Number of ${\pi }^0$ candidates	11.60%	—
${\eta }_c$ mass and width	0.99%	—
$h_c$ branching fraction	10.22%	—
$ℬ\left({\rm Y}\left(4S\right)\to B^{+}{B}^{-}\right)$	1.17%
Number of ${\rm Y}\left(4S\right)$ events	1.37%
Total	$\left(\stackrel{⃡}{\ast }{[}_{-22.54}^{+22.51}]\right)\%$	$\left(\stackrel{⃡}{\ast }{[}_{-4.67}^{+5.62}]\right)\%$	$\left(\stackrel{⃡}{\ast }{[}_{-4.41}^{+4.83}]\right)\%$	$\left(\stackrel{⃡}{\ast }{[}_{-6.86}^{+6.30}]\right)\%$	$\left(\stackrel{⃡}{\ast }{[}_{-5.25}^{+8.31}]\right)\%$	$\left(\stackrel{⃡}{\ast }{[}_{-66.07}^{+21.57}]\right)\%$	$\left(\stackrel{⃡}{\ast }{[}_{-6.34}^{+4.68}]\right)\%$	$\left(\stackrel{⃡}{\ast }{[}_{-8.29}^{+4.68}]\right)\%$

Table 6. Table 6: Relative systematic uncertainties of the branching fractions for the channel B 0 → ( c ​ c ¯ ) ​ K S 0 → superscript 𝐵 0 𝑐 ¯ 𝑐 subscript superscript 𝐾 0 𝑆 B^{0}\to(c\bar{c})K^{0}_{S} ( B 0 → h c ​ K S 0 → superscript 𝐵 0 subscript ℎ 𝑐 subscript superscript 𝐾 0 𝑆 B^{0}\to h_{c}K^{0}_{S} for the h c subscript ℎ 𝑐 h_{c} and B 0 → ( c c ¯ ) ( → p p ¯ π + π − ) K S 0 B^{0}\to(c\bar{c})(\to p\bar{p}\pi^{+}\pi^{-})K^{0}_{S} for all other charmonium states).

Error source	$h_c$	${\eta }_c$	$J/\psi$	${\chi }_{c0}$	${\chi }_{c1}$	${\chi }_{c2}$	${\eta }_c\left(2S\right)$	$\psi \left(2S\right)$
Model dependence	$\left(\stackrel{⃡}{\ast }{[}_{-24.05}^{+30.94}]\right)\%$	$\left(\stackrel{⃡}{\ast }{[}_{-24.37}^{+4.91}]\right)\%$	$\left(\stackrel{⃡}{\ast }{[}_{-1.80}^{+1.45}]\right)\%$	$\left(\stackrel{⃡}{\ast }{[}_{-24.56}^{+134.01}]\right)\%$	$\left(\stackrel{⃡}{\ast }{[}_{-2.28}^{+6.79}]\right)\%$	$\left(\stackrel{⃡}{\ast }{[}_{-4.78}^{+8.63}]\right)\%$	$\left(\stackrel{⃡}{\ast }{[}_{-5.36}^{+4.51}]\right)\%$	$\left(\stackrel{⃡}{\ast }{[}_{-9.01}^{+10.89}]\right)\%$
PID	4.86%	3.93%	3.93%	3.87%	3.87%	3.86%	3.83%	3.81%
Overtraining	0.15%	0.19%
Tracking	1.95%	2.10%
MLP efficiency	12.79%	0.25%
Number of ${\pi }^0$ candidates	11.66%	—
${\eta }_c$ mass and width	0.96%	—
$h_c$ branching fraction	10.27%	—
$ℬ\left({\rm Y}\left(4S\right)\to B^0{\overline{B}}^0\right)$	1.23%
Number of ${\rm Y}\left(4S\right)$ events	1.37%
Total	$\left(\stackrel{⃡}{\ast }{[}_{-31.86}^{+37.33}]\right)\%$	$\left(\stackrel{⃡}{\ast }{[}_{-24.85}^{+6.89}]\right)\%$	$\left(\stackrel{⃡}{\ast }{[}_{-5.16}^{+5.04}]\right)\%$	$\left(\stackrel{⃡}{\ast }{[}_{-25.02}^{+134.10}]\right)\%$	$\left(\stackrel{⃡}{\ast }{[}_{-5.30}^{+8.31}]\right)\%$	$\left(\stackrel{⃡}{\ast }{[}_{-6.76}^{+9.87}]\right)\%$	$\left(\stackrel{⃡}{\ast }{[}_{-7.17}^{+6.56}]\right)\%$	$\left(\stackrel{⃡}{\ast }{[}_{-10.18}^{+11.88}]\right)\%$

Table 7. Table 7: Measured branching fractions and branching fraction products and their comparison with the current world-average values  Tanabashi:2018oca .

Branching fraction	Value or confidence interval (90 % C. L.)	World-average value
$ℬ\left(B^{+}\to h_c{K}^{+}\right)$	$\left(3.7\stackrel{⃡}{\ast }{[}_{-0.9}^{+1.0}]\stackrel{⃡}{\ast }{[}_{-0.8}^{+0.8}]\right)\times {10}^{-5}$
$ℬ\left(B^{+}\to {\eta }_c{K}^{+}\right)\times ℬ\left({\eta }_c\to p\overline{p}{\pi }^{+}{\pi }^{-}\right)$	$\left(39.4\stackrel{⃡}{\ast }{[}_{-3.9}^{+4.1}]\stackrel{⃡}{\ast }{[}_{-1.8}^{+2.2}]\right)\times {10}^{-7}$	$\left(57.8\pm 20.2\right)\times {10}^{-7}$
$ℬ\left(B^{+}\to J/\psi K^{+}\right)\times ℬ\left(J/\psi \to p\overline{p}{\pi }^{+}{\pi }^{-}\right)$	$\left(56.4\stackrel{⃡}{\ast }{[}_{-3.2}^{+3.3}]\stackrel{⃡}{\ast }{[}_{-2.5}^{+2.7}]\right)\times {10}^{-7}$	$\left(60.6\pm 5.3\right)\times {10}^{-7}$
$ℬ\left(B^{+}\to {\chi }_{c0}{K}^{+}\right)\times ℬ\left({\chi }_{c0}\to p\overline{p}{\pi }^{+}{\pi }^{-}\right)$	$\left(3.7\stackrel{⃡}{\ast }{[}_{-1.0}^{+1.2}]\stackrel{⃡}{\ast }{[}_{-0.3}^{+0.2}]\right)\times {10}^{-7}$	$\left(3.1\pm 1.1\right)\times {10}^{-7}$
$ℬ\left(B^{+}\to {\chi }_{c1}{K}^{+}\right)\times ℬ\left({\chi }_{c1}\to p\overline{p}{\pi }^{+}{\pi }^{-}\right)$	$\left(4.7\stackrel{⃡}{\ast }{[}_{-1.2}^{+1.3}]\stackrel{⃡}{\ast }{[}_{-0.2}^{+0.4}]\right)\times {10}^{-7}$	$\left(2.4\pm 0.9\right)\times {10}^{-7}$
$ℬ\left(B^{+}\to {\chi }_{c2}{K}^{+}\right)\times ℬ\left({\chi }_{c2}\to p\overline{p}{\pi }^{+}{\pi }^{-}\right)$		$\left(0.15\pm 0.06\right)\times {10}^{-7}$
$ℬ\left(B^{+}\to {\eta }_c\left(2S\right)K^{+}\right)\times ℬ\left({\eta }_c\left(2S\right)\to p\overline{p}{\pi }^{+}{\pi }^{-}\right)$	$\left(11.2\stackrel{⃡}{\ast }{[}_{-1.6}^{+1.8}]\stackrel{⃡}{\ast }{[}_{-0.7}^{+0.5}]\right)\times {10}^{-7}$	not seen
$ℬ\left(B^{+}\to \psi \left(2S\right)K^{+}\right)\times ℬ\left(\psi \left(2S\right)\to p\overline{p}{\pi }^{+}{\pi }^{-}\right)$	$[0.5,3.5]\times {10}^{-7}$	$\left(3.7\pm 0.3\right)\times {10}^{-7}$
$ℬ\left(B^0\to h_c{K}_S^0\right)$		not seen
$ℬ\left(B^0\to {\eta }_c{K}_S^0\right)\times ℬ\left({\eta }_c\to p\overline{p}{\pi }^{+}{\pi }^{-}\right)$	$\left(19.0\stackrel{⃡}{\ast }{[}_{-2.9}^{+3.2}]\stackrel{⃡}{\ast }{[}_{-4.7}^{+1.3}]\right)\times {10}^{-7}$	$\left(20.9\pm 7.8\right)\times {10}^{-7}$
$ℬ\left(B^0\to J/\psi K_S^0\right)\times ℬ\left(J/\psi \to p\overline{p}{\pi }^{+}{\pi }^{-}\right)$	$\left(24.3\stackrel{⃡}{\ast }{[}_{-2.2}^{+2.3}]\stackrel{⃡}{\ast }{[}_{-1.3}^{+1.2}]\right)\times {10}^{-7}$	$\left(26.2\pm 2.4\right)\times {10}^{-7}$
$ℬ\left(B^0\to {\chi }_{c0}{K}_S^0\right)\times ℬ\left({\chi }_{c0}\to p\overline{p}{\pi }^{+}{\pi }^{-}\right)$		$\left(1.5\pm 0.6\right)\times {10}^{-7}$
$ℬ\left(B^0\to {\chi }_{c1}{K}_S^0\right)\times ℬ\left({\chi }_{c1}\to p\overline{p}{\pi }^{+}{\pi }^{-}\right)$	$\left(3.7\stackrel{⃡}{\ast }{[}_{-1.0}^{+1.2}]\stackrel{⃡}{\ast }{[}_{-0.2}^{+0.3}]\right)\times {10}^{-7}$	$\left(1.0\pm 0.4\right)\times {10}^{-7}$
$ℬ\left(B^0\to {\chi }_{c2}{K}_S^0\right)\times ℬ\left({\chi }_{c2}\to p\overline{p}{\pi }^{+}{\pi }^{-}\right)$	$[0.7,3.8]\times {10}^{-7}$	not seen
$ℬ\left(B^0\to {\eta }_c\left(2S\right)K_S^0\right)\times ℬ\left({\eta }_c\left(2S\right)\to p\overline{p}{\pi }^{+}{\pi }^{-}\right)$	$\left(4.2\stackrel{⃡}{\ast }{[}_{-1.2}^{+1.4}]\stackrel{⃡}{\ast }{[}_{-0.3}^{+0.3}]\right)\times {10}^{-7}$	not seen
$ℬ\left(B^0\to \psi \left(2S\right)K_S^0\right)\times ℬ\left(\psi \left(2S\right)\to p\overline{p}{\pi }^{+}{\pi }^{-}\right)$		$\left(1.7\pm 0.2\right)\times {10}^{-7}$

Equations

$$\sqrt{{\left(\frac{\Delta E}{{\sigma }_{\Delta E}}\right)}^2+{\left(\frac{M_{\text{bc}}-m_B}{{\sigma }_{M_{\text{bc}}}}\right)}^2}<3,$$

$$3.50<M_{h_c}<3.55\text{GeV}/c^2$$

$$3.515<M_{h_c}<3.535\text{GeV}/c^2$$

$$S(\Delta E,M_{\text{bc}})=$$

$$+N_{\text{G1}}{G}_a^{(21)}(x_2)G_a^{(22)}(y_2)$$

$$+N_{\text{G2}}{G}_a^{(31)}(x_3)G_a^{(32)}(y_3),$$

$$\left(\begin{array}{c}{x}_i\\ y_i\end{array}\right)=\left(\begin{array}{cc}\cos {\alpha }_i& \sin {\alpha }_i\\ -\sin {\alpha }_i& \cos {\alpha }_i\end{array}\right)\left(\begin{array}{c}\Delta E-(\Delta E{)}_0\\ M_{\text{bc}}-(M_{\text{bc}}{)}_0\end{array}\right).$$

$$N_{S}S(\Delta E,M_{\text{bc}})+B(\Delta E,M_{\text{bc}}),$$

$$B(\Delta E,M_{\text{bc}})=$$

$$\times P_3(\Delta E,M_{\text{bc}}),$$

$$3<\sqrt{{\left(\frac{\Delta E}{{\sigma }_{\Delta E}}\right)}^2+{\left(\frac{M_{\text{bc}}-m_B}{{\sigma }_{M_{\text{bc}}}}\right)}^2}<8.$$

$$ϵ(v_0)={ϵ}_{\text{MLP}}(v_0)\times {ϵ}_{\text{multiple}}(v_0),$$

$$F_{\text{opt}}=\frac{\sum _i{N}_{\text{sig}}^{(i)}}{\frac{a}{2}+\sqrt{\sum _i{N}_{\text{bg}}^{(i)}}},$$

$${\left(\frac{\Delta E}{R_{\Delta E}^{(i)}}\right)}^2+{\left(\frac{M_{\text{bc}}-m_B}{R_{M_{\text{bc}}}^{(i)}}\right)}^2<1.$$

$$N_{\text{sig}}^{(i)}=$$

$$\times \mathcal{B}(h_c\to i){ϵ}_{\text{SR}}^{(i)}{ϵ}_S^{(i)}(v_0^{(i)}),$$

$${ϵ}_{\text{SR}}^{(i)}={ϵ}_R^{(i)}\underset{\text{SR}}{\int }{S}_i(\Delta E,M_{\text{bc}})d\Delta Ed{M}_{\text{bc}},$$

$$N_{\text{bg}}^{(i)}={ϵ}_B^{(i)}(v_0^{(i)})\frac{N_{h_c\text{ region}}}{N_{\text{full}}}\underset{\text{SR}}{\int }{B}_i(\Delta E,M_{\text{bc}})d\Delta Ed{M}_{\text{bc}},$$

$$R_{h_c}(\Delta M)=N[F_{\text{CB}}(\Delta M)f_{\text{CB}}+G_a(\Delta M)(1-f_{\text{CB}})],$$

$$A_R(M_R)=\frac{1}{M_R^2-m_R^2+i{M}_R{\Gamma }_R},$$

$$S_{{\eta }_c\gamma }(M)=\left(N_{h_c}|A_{h_c}(M){|}^2\right)\otimes R_{h_c}^{({\eta }_c\gamma )}(\Delta M)+P_2(M),$$

$$S_{p\bar{p}{\pi }^{+}{\pi }^{-}}(M)=$$

$$\displaystyle\qquad\Big{(}|P_{3}(M)+\sum\limits_{R=\eta_{c},\chi_{c0},\eta_{c}(2S)}\sqrt{N_{R}}e^{i\varphi_{R}}A_{R}(M)|^{2}$$

$$\displaystyle\qquad+\sum\limits_{R=J/\psi,\chi_{c1},h_{c},\chi_{c2},\psi(2S)}N_{R}|A_{R}(M)|^{2}\Big{)}$$

$$\otimes R_{h_c}^{(p\bar{p}{\pi }^{+}{\pi }^{-})}(\Delta M),$$

$$S_{p\bar{p}{\pi }^{+}{\pi }^{-}}(M)+w{B}_{p\bar{p}{\pi }^{+}{\pi }^{-}}(M),$$

$$R_{h_c}(\Delta M)\to \frac{1}{\mathcal{S}}{R}_{h_c}\left(\frac{\Delta M}{\mathcal{S}}\right),$$

$$r_{\text{MLP}}=\frac{N(\text{all ηc channels, with MLP cut})}{N({\eta }_c\to K^{+}{K}_S^0{\pi }^{-},\text{no MLP cut})}.$$

$$2N_{{\rm Y}(4S)}\left(\sum _i{ϵ}_i{\mathcal{B}}_i\right)=2N_{{\rm Y}(4S)}{\mathcal{B}}_{h_c\to {\eta }_c\gamma }\left(\sum _i{ϵ}_i\frac{{\mathcal{B}}_i}{{\mathcal{B}}_{h_c\to {\eta }_c\gamma }}\right),$$

Full text

I Introduction

The decays $B^{+}\to\chi_{c0}K^{+}$, $B^{+}\to\chi_{c2}K^{+}$ and $B^{+}\to h_{c}K^{+}$ are suppressed by factorization Bauer:1986bm ; Suzuki:2002sq . The decays $B^{+}\to\chi_{cJ}K^{+}$ have been observed; the current world-average branching fractions are $\mathcal{B}(B^{+}\to\chi_{c0}K^{+})=(1.49\tensor{*}[^{+0.15}_{-0.14}]{}{})\times 10^{-4}$ and $\mathcal{B}(B^{+}\to\chi_{c2}K^{+})=(1.1\pm 0.4)\times 10^{-5}$ Tanabashi:2018oca . While $\mathcal{B}(B^{+}\to\chi_{c0}K^{+})$ is smaller than the branching fraction of the factorization-allowed process $\mathcal{B}(B^{+}\to\chi_{c1}K^{+})=(4.84\pm 0.23)\times 10^{-4}$, it is not strongly suppressed. Before the first experimental searches, this resulted in an assumption that the process $B^{+}\to h_{c}K^{+}$ may also have a large branching fraction $\mathcal{B}(B^{+}\to h_{c}K^{+})\approx\mathcal{B}(B^{+}\to\chi_{c0}K^{+})$ Suzuki:2002sq .

However, the decay $B^{+}\to h_{c}K^{+}$ has not been observed experimentally yet. Neither the Belle Fang:2006bz nor BaBar Aubert:2008kp Collaborations have found a statistically significant signal of $B^{+}\to h_{c}K^{+}$ using the decay mode $h_{c}\to\eta_{c}\gamma$. The current branching-fraction upper limit of $\mathcal{B}(B^{+}\to h_{c}K^{+})<3.8\times 10^{-5}$ at $90\%$ confidence level (C. L.) Tanabashi:2018oca was obtained in the $h_{c}$ search by Belle Fang:2006bz .

Also, the LHCb Collaboration searched for the process $B^{+}\to h_{c}(\to p\bar{p})K^{+}$ Aaij:2013rha and set the upper limit on the branching fraction product $\mathcal{B}(B^{+}\to h_{c}K^{+})\times\mathcal{B}(h_{c}\to p\bar{p})<6.4\times 10^{-8}$ (95% C. L.). However, this measurement does not result in a stronger restriction on $\mathcal{B}(B^{+}\to h_{c}K^{+})$, because the decay $h_{c}\to p\bar{p}$ has never been observed and the upper limit on its branching fraction is $\mathcal{B}(h_{c}\to p\bar{p})<1.5\times 10^{-4}$ (90% C. L.) Tanabashi:2018oca . Note that a newer LHCb analysis of the same channel performed in Ref. Aaij:2016kxn does not update the upper limit on $\mathcal{B}(B^{+}\to h_{c}K^{+})\times\mathcal{B}(h_{c}\to p\bar{p})$.

Several new theoretical predictions of $\mathcal{B}(B^{+}\to h_{c}K^{+})$ were made after the experimental upper limit was set. The branching fraction has been calculated in the QCD factorization approach to be $2.7\times 10^{-5}$ Meng:2006mi . A calculation using perturbative QCD was performed in Ref. Li:2006vj ; the result is $\mathcal{B}(B^{+}\to h_{c}K^{+})=3.6\times 10^{-5}$. Another calculation performed in Ref. Beneke:2008pi results in $\mathcal{B}(B^{+}\to h_{c}K^{+})$ [$\mathcal{B}(B^{0}\to h_{c}K^{0})$] in the interval from $3.1\times 10^{-5}$ to $5.7\times 10^{-5}$ [from $2.9\times 10^{-5}$ to $5.3\times 10^{-5}$], depending on the assumed value of the $c$ quark mass. All the results mentioned above are close to each other, and the theoretical values of $\mathcal{B}(B^{+}\to h_{c}K^{+})$ are slightly below the current experimental upper limit. This motivates an updated study of the decays $B^{+}\to h_{c}K^{+}$, which may be able to find them.

Here we present such an updated search for the decays $B^{+}\to h_{c}K^{+}$, and also include a search for the decays $B^{0}\to h_{c}K^{0}_{S}$. The analysis is performed using the $711\ \text{fb}^{-1}$ data sample collected by the Belle detector at the asymmetric-energy $e^{+}e^{-}$ collider KEKB kekb . The data sample was collected at the $\Upsilon(4S)$ resonance and contains $772\times 10^{6}$ $B\bar{B}$ pairs. The integrated luminosity is 2.8 times greater than the luminosity used in the previous analysis Fang:2006bz . For further improvement of the sensitivity, the new analysis uses ten $\eta_{c}$ decay channels to reconstruct the decay $h_{c}\to\eta_{c}\gamma$; only two channels were used in the old one. The new $h_{c}$ decay channel $h_{c}\to p\bar{p}\pi^{+}\pi^{-}$ observed recently by the BESIII Collaboration Ablikim:2018ewr is also used for its reconstruction; in addition, we study the decays of other charmonium states to $p\bar{p}\pi^{+}\pi^{-}$. The discrimination of the signal and background events is improved by performing a multivariate analysis.

II The Belle Detector

The Belle detector is a large-solid-angle magnetic spectrometer that consists of a silicon vertex detector (SVD), a 50-layer central drift chamber (CDC), an array of aerogel threshold Cherenkov counters (ACC), a barrel-like arrangement of time-of-flight scintillation counters (TOF), and an electromagnetic calorimeter (ECL) comprised of CsI(Tl) crystals located inside a superconducting solenoid coil that provides a 1.5 T magnetic field. An iron flux-return located outside of the coil is instrumented to detect $K_{L}^{0}$ mesons and to identify muons. The detector is described in detail elsewhere Belle . Two inner detector configurations were used. A 2.0 cm radius beampipe and a 3-layer silicon vertex detector were used for the first sample of 140 $\text{fb}^{-1}$, while a 1.5 cm radius beampipe, a 4-layer silicon detector and a small-cell inner drift chamber were used to record the remaining data svd2 .

We use a geant-based Monte Carlo (MC) simulation geant to model the response of the detector, identify potential backgrounds and determine the acceptance. The MC simulation includes run-dependent detector performance variations and background conditions. Signal MC events are generated with EvtGen evtgen in proportion to the relative luminosities of the different running periods.

III Event selection

We select events of the type $B^{+}\to h_{c}K^{+}$ and $B^{0}\to h_{c}K^{0}_{S}$. Inclusion of charge-conjugate modes is implied hereinafter. The reconstruction is performed with a conversion from Belle to Belle II data format Gelb:2018agf .

All tracks are required to originate from the interaction point region: we require $dr<0.2\ \text{cm}$ and $|dz|<2\ \text{cm}$, where $dr$ and $dz$ are the cylindrical coordinates of the point of the closest approach of the track to the beam axis (the $z$ axis of the laboratory reference frame coincides with the positron-beam axis).

Charged $\pi$, $K$ mesons and protons are identified using likelihood ratios $R_{h_{1}/h_{2}}=\mathcal{L}_{h_{1}}/(\mathcal{L}_{h_{1}}+\mathcal{L}_{h_{2}})$, where $h_{1}$ and $h_{2}$ are the particle-identification hypotheses ($\pi$, $K$, or $p$) and $\mathcal{L}_{h_{i}}$ are their corresponding likelihoods. The likelihoods are calculated from the combined time-of-flight information from the TOF, the number of photoelectrons from the ACC and $dE/dx$ measurements in the CDC. We require $R_{K/\pi}>0.6$ for $K$ candidates, $R_{\pi/K}>0.6$ for $\pi$ candidates, and $R_{p/\pi}>0.6$, $R_{p/K}>0.6$ for $p$ candidates. The identification efficiency of the above requirements varies in the ranges (94 – 99)%, (84 – 93)%, and (90 – 98)% for $\pi$, $K$, and $p$, respectively, depending on the $h_{c}$ or $\eta_{c}$ decay channel. The misidentification probability for the background particles that are not $\pi$, $K$, and $p$, varies in the ranges (25 – 49)%, (4.9 – 11.3)%, and (0.5 – 1.9)%, respectively. Without the electron background, which is rejected as described below, the $\pi$ fake rate drops to (20 – 35)%.

Electron candidates are identified as CDC charged tracks that are matched to electromagnetic showers in the ECL. The track and ECL cluster matching quality, the ratio of the electromagnetic shower energy to the track momentum, the transverse shape of the shower, the ACC light yield, and the track $dE/dx$ ionization are used in our electron identification criteria. A similar likelihood ratio is constructed: $R_{e}=\mathcal{L}_{e}/(\mathcal{L}_{e}+\mathcal{L}_{h})$, where $\mathcal{L}_{e}$ and $\mathcal{L}_{h}$ are the likelihoods for electrons and charged hadrons ($\pi$, $K$ and $p$), respectively Hanagaki:2001fz . An electron veto ($R_{e}<0.9$) is imposed on $\pi$, $K$, and $p$ candidates. It is not applied for the $K^{0}_{S}$ and $\Lambda$ daughter tracks, because they have independent selection criteria. For the $h_{c}$ or $\eta_{c}$ decay channels other than $\eta_{c}\to K^{0}_{S}K^{0}_{S}\pi^{0}$ and $\eta_{c}\to\Lambda\bar{\Lambda}$, the electron veto rejects from 3.5 to 15% of the background events, while its signal efficiency is not less than 97.5%.

Photons are identified as ECL electromagnetic showers that have no associated charged tracks detected in the CDC. The shower shape is required to be consistent with that of a photon.

The $\pi^{0}$ candidates are reconstructed via their decay to two photons. The photon energies in the laboratory frame are required to be greater than $30\ \text{MeV}$. The $\pi^{0}$ invariant mass is required to satisfy $|M_{\pi^{0}}-m_{\pi^{0}}|<15\ \text{MeV}/c^{2}$. Here and elsewhere, $M_{\text{particle}}$ denotes the reconstructed invariant mass of the specified particle and $m_{\text{particle}}$ stands for the nominal mass of this particle Tanabashi:2018oca . This requirement corresponds approximately to a $3\sigma$ mass window around the nominal mass.

The $V^{0}$-particle ($K^{0}_{S}$ and $\Lambda$) candidates are reconstructed from pairs of oppositely charged tracks that are assumed to be $\pi^{+}\pi^{-}$ and $p\pi^{-}$ for $K^{0}_{S}$ and $\Lambda$, respectively. We require $|M_{K^{0}_{S}}-m_{K^{0}_{S}}|<20\ \text{MeV}/c^{2}$ and $|M_{\Lambda}-m_{\Lambda}|<10\ \text{MeV}/c^{2}$, corresponding approximately to $5.5\sigma$ mass windows in both cases. The $V^{0}$ candidates are selected by a neural network using the following input variables: the $V^{0}$ candidate momentum, decay angle, flight distance in the $xy$ plane, the angle between the $V^{0}$ momentum and the direction from the interaction point to the $V^{0}$ vertex, the shortest $z$ distance between the two daughter tracks, their radial impact parameters, and numbers of hits in the SVD and CDC. The separation of the $K^{0}_{S}$ and $\Lambda$ candidates is performed by another neural network. The input variables of this network are the momenta and polar angles of the daughter tracks in the laboratory frame, their likelihood ratios $R_{\pi/p}$ and the $V^{0}$ candidate invariant mass for the $\Lambda$ hypothesis.

The $\eta$ candidates are reconstructed in $\gamma\gamma$ and $\pi^{+}\pi^{-}\pi^{0}$ channels. The reconstructed $\eta$ candidates are denoted by the $\eta$ decay channel as $\eta_{2\gamma}$ and $\eta_{3\pi}$. The $\eta$ invariant mass is required to satisfy $|M_{\eta_{2\gamma}}-m_{\eta}|<30\ \text{MeV}/c^{2}$ and $|M_{\eta_{3\pi}}-m_{\eta}|<15\ \text{MeV}/c^{2}$; these requirements correspond to $2.5\sigma$ and $4\sigma$ mass windows, respectively.

The $\eta^{\prime}$ candidates are reconstructed in the $\eta\pi^{+}\pi^{-}$ decay mode. The invariant mass window is $|M_{\eta^{\prime}}-m_{\eta^{\prime}}|<15\ \text{MeV}/c^{2}$, corresponding to a $4\sigma$ mass window.

The $\eta_{c}$ candidates are reconstructed in ten decay channels: $K^{+}K^{0}_{S}\pi^{-}$, $K^{+}K^{-}\pi^{0}$, $K^{0}_{S}K^{0}_{S}\pi^{0}$, $K^{+}K^{-}\eta$, $K^{+}K^{-}K^{+}K^{-}$, $\eta^{\prime}(\to\eta\pi^{+}\pi^{-})\pi^{+}\pi^{-}$, $p\bar{p}$, $p\bar{p}\pi^{0}$, $p\bar{p}\pi^{+}\pi^{-}$, and $\Lambda\bar{\Lambda}$. The selected $\eta_{c}$ candidates are required to satisfy $|M_{\eta_{c}}-m_{\eta_{c}}|<50\ \text{MeV}/c^{2}$; the mass-window width is about $1.6$ widths of the $\eta_{c}$.

The $h_{c}$ candidates are reconstructed in the $h_{c}\to\eta_{c}\gamma$ and $h_{c}\to p\bar{p}\pi^{+}\pi^{-}$ decay channels. The invariant mass of the $h_{c}$ candidates is not restricted for the channel $\eta_{c}\gamma$; for the channel $p\bar{p}\pi^{+}\pi^{-}$, it is required to be greater than $2.7\ \text{GeV}/c^{2}$. The lower mass limit is selected to be very low to study other charmonium states decaying to the same final state.

The $B$-meson candidates are reconstructed via the decay modes $B^{+}\to h_{c}K^{+}$ and $B^{0}\to h_{c}K^{0}_{S}$. The $B$ candidates are selected by their energy and the beam-energy-constrained mass. The difference of the $B$-meson and beam energies is defined as $\Delta E=\sum_{i}E_{i}-E_{\text{beam}}$, where $E_{i}$ are the energies of the $B$ decay products in the center-of-mass frame and $E_{\text{beam}}$ is the beam energy in the same frame. The beam-energy-constrained mass is defined as $M_{\text{bc}}=\sqrt{E_{\text{beam}}^{2}-(\sum_{i}\vec{p}_{i})^{2}}$, where $\vec{p}_{i}$ are the momenta of the $B$ decay products in the center-of-mass frame. We retain $B$ candidates satisfying the conditions $5.2<M_{\text{bc}}<5.3\ \text{GeV}/c^{2}$ and $|\Delta E|<0.2\ \text{GeV}$. A mass-constrained fit is applied to the selected $B$-meson candidates.

In addition, for the channel $h_{c}\to\eta_{c}\gamma$, the $h_{c}$ daughter $\gamma$ energy is required to be greater than $200\ \text{MeV}$ in the $B$ rest frame. This requirement removes the background from low-energy photons, including the peaking backgrounds from $B$ decays to the same final state without the photon. The signal efficiency of this requirement is 100%, because the $\eta_{c}\gamma$ invariant mass of all excluded events is smaller than the $h_{c}$ mass.

Also, the ratio of the Fox-Wolfram moments Fox:1978vu $F_{2}/F_{0}$ is required to be less than 0.3. This requirement reduces the continuum background, rejecting from 18% to 53% of background events, depending on the $h_{c}$ or $\eta_{c}$ decay channel. Its signal efficiency is from 93.3% to 96.3%.

IV Multivariate analysis and optimization of the selection requirements

IV.1 General analysis strategy and data samples

To improve the separation of the signal and background events, we perform a multivariate analysis followed by an optimization of selection requirements. The first stages of the analysis are performed individually for $h_{c}\to p\bar{p}\pi^{+}\pi^{-}$ and each $\eta_{c}$ decay channel for the $h_{c}$ candidates reconstructed in the $\eta_{c}\gamma$ mode (the channels $\eta_{c}\to K^{+}K^{-}\eta$ and $\eta_{c}\to\eta^{\prime}(\to\eta\pi^{+}\pi^{-})\pi^{+}\pi^{-}$ are optimized separately for $\eta_{2\gamma}$ and $\eta_{3\pi}$). They include the determination of two-dimensional $(\Delta E,M_{\text{bc}})$ resolution and the distribution of the background in $(\Delta E,M_{\text{bc}})$, and the multivariate-analysis stage. The optimization of the selection requirements uses the results of all initial stages as its input. The resolution is used to determine the expected number of the signal events and the distribution of the background in $(\Delta E,M_{\text{bc}})$ is used to determine the expected number of the background events in the signal region. The optimization is performed individually for the channel $h_{c}\to p\bar{p}\pi^{+}\pi^{-}$ and globally for all $\eta_{c}$ decay channels for the channel $h_{c}\to\eta_{c}\gamma$. The data selected using the resulting channel-dependent criteria are merged into a single sample for the $h_{c}\to\eta_{c}\gamma$ channel. The final fit is performed simultaneously to the $h_{c}\to\eta_{c}\gamma$ and $h_{c}\to p\bar{p}\pi^{+}\pi^{-}$ samples.

The experimental data are used for determination of the $(\Delta E,M_{\text{bc}})$ distribution, selection of the background samples for the neural network, and final fit to the selected events. During the development of the analysis procedure, the $h_{c}$ signal region was excluded to avoid bias of the $h_{c}$ significance. The final fit described in Sec. V was performed on MC pseudoexperiments generated in accordance with the fit result without the $h_{c}$ mixed with the $B^{+}\to h_{c}K^{+}$ ($B^{0}\to h_{c}K^{0}_{S}$) signal MC. The $h_{c}$ signal region is defined by

[TABLE]

where $\sigma_{\Delta E}=18\ \text{MeV}$ and $\sigma_{M_{\text{bc}}}=2.5\ \text{MeV}/c^{2}$ are the approximate resolutions in $\Delta E$ and $M_{\text{bc}}$, respectively, and

[TABLE]

After completion of the analysis procedure development, this requirement is no longer used.

The signal MC is used for the determination of the resolution and the selection of the signal samples for the neural network. The signal MC is generated using the known information about the angular or invariant-mass distributions of the decay products if it is possible; otherwise, uniform distributions are assumed. The angular distribution is known for the channel $h_{c}\to\eta_{c}\gamma$. It does not have any free parameters and is proportional to $\sin^{2}\theta_{h_{c}}$, where $\theta_{h_{c}}$ is the $h_{c}$ helicity angle that is defined as the angle between $-\vec{p}_{B}$ and $\vec{p}_{\eta_{c}}$, where $\vec{p}_{B}$ and $\vec{p}_{\eta_{c}}$ are the momenta of the $B$ and $\eta_{c}$ in the $h_{c}$ rest frame, respectively. In addition, the $\eta_{c}$ decay resonant structure is taken into account if it is known. The distributions for the channels $K^{+}K^{0}_{S}\pi^{-}$, $K^{+}K^{-}\pi^{0}$, $K^{0}_{S}K^{0}_{S}\pi^{0}$, $K^{+}K^{-}\eta_{2\gamma}$, and $K^{+}K^{-}\eta_{3\pi}$ are based on the results of a Dalitz plot analysis performed in Ref. Lees:2014iua . The contributions of intermediate $\phi$ resonances are taken into account for the channel $K^{+}K^{-}K^{+}K^{-}$ based on the world-average branching fractions from Ref. Tanabashi:2018oca .

IV.2 Resolution

The resolution is parameterized by the function

[TABLE]

where $F_{\text{CB}}$ is an asymmetric Crystal Ball function skwarnicki , $G_{a}^{(ij)}$ are asymmetric Gaussian functions, $N_{\text{CB}}$, $N_{\text{G1}}$ and $N_{\text{G2}}$ are normalizations and $x_{i}$ and $y_{i}$ ($i$ = 1, 2, 3) are rotated variables that are given by

[TABLE]

Here, ($(\Delta E)_{0}$, $(M_{\text{bc}})_{0}$) is the central point and $\alpha_{i}$ is the rotation angle. The central point is the same for all three components. The resolution is determined from a binned maximum likelihood fit to signal MC events. Example resolution fit results (for the channel $B^{+}\to h_{c}K^{+}$ with $h_{c}\to\eta_{c}(\to K^{+}K^{-}\pi^{0})\gamma$) are shown in Fig. 1.

IV.3 Fit to the $(\Delta E,M_{\text{bc}})$ distribution

The $(\Delta E,M_{\text{bc}})$ distribution is fitted in order to estimate the expected number of the background events in the signal region. The distribution is fitted to the function

[TABLE]

where $N_{S}$ is the number of signal events and $B$ is the background density function that is given by

[TABLE]

where $m_{0}$ is the threshold mass, $a$ is a rate parameter and $P_{3}$ is a two-dimensional third-order polynomial. The region with $\Delta E<-0.12\ \text{GeV}$ is excluded for the channel $h_{c}\to p\bar{p}\pi^{+}\pi^{-}$ because of the presence of peaking backgrounds from partially reconstructed B decays with an additional $\pi$ meson.

Example $(\Delta E,M_{\text{bc}})$ fit results (for the channel $B^{+}\to h_{c}K^{+}$ with $h_{c}\to\eta_{c}(\to K^{+}K^{-}\pi^{0})\gamma$) are shown in Fig. 2.

IV.4 Multivariate analysis

To improve the separation of signal and background events, we perform a multivariate analysis for each individual channel. The algorithm used for the multivariate analysis is the multilayer perceptron (MLP) neural network implemented in the tmva library tmva . The following variables are always included in the neural network: the angle between the thrust axes of the $B$ candidate and the remaining particles in the event, the angle between the thrust axes of all tracks and all photons in the event, the ratio of the Fox-Wolfram moments $F_{2}/F_{0}$, the $B$ production angle, and the vertex fit quality. For the $h_{c}$ candidates reconstructed in the $\eta_{c}\gamma$ channel, the MLP also includes the $h_{c}$ helicity angle, the $\eta_{c}$ mass, and the number of $\pi^{0}$ candidates that include the $h_{c}$ daughter photon as one of their daughters (separately for two groups of $\pi^{0}$ candidates with the energy of another photon less and greater than $100\ \text{MeV}$).

For the channels $\eta_{c}\to K^{+}K^{0}_{S}\pi^{-}$, $\eta_{c}\to K^{+}K^{-}\pi^{0}$, and $\eta_{c}\to K^{0}_{S}K^{0}_{S}\pi^{0}$, two invariant masses of the $\eta_{c}$ daughter particle pairs (both $K\pi$ combinations) are added to the neural network.

The following particle identification variables are included into the neural network if there are corresponding charged particles in the final state: the minimum likelihood ratio $R_{K/\pi}$ of the $h_{c}$ daughter kaons, the minimum of the two likelihood ratios $R_{p/K}$, $R_{p/\pi}$ of the $h_{c}$ daughter protons, and $R_{K/\pi}$ for the $B$ daughter $K^{+}$ (for the channel $B^{+}\to h_{c}K^{+}$). Here, the $h_{c}$ daughters may be either direct (from the decay $h_{c}\to p\bar{p}\pi^{+}\pi^{-}$) or indirect (the $\eta_{c}$ daughters for the $h_{c}$ candidates reconstructed in the $\eta_{c}\gamma$ mode).

If there is a $\pi^{0}$ or $\eta$ decaying to $\gamma\gamma$ in the final state, four additional variables are added: the $\pi^{0}$ ($\eta$) mass, the minimal energy of the $\pi^{0}$ ($\eta$) daughter photons in the laboratory frame, and the number of $\pi^{0}$ candidates that include a $\pi^{0}$ ($\eta$) daughter photon as one of their daughters (for each of the $\pi^{0}$ ($\eta$) daughter photons). If there is an $\eta$ reconstructed in the $\pi^{+}\pi^{-}\pi^{0}$ decay mode, then only its mass is added to the MLP. If the $\eta_{c}$ has a daughter $\eta^{\prime}$, then the mass of the $\eta^{\prime}$ candidate is also included to the neural network.

The training and testing signal samples are taken from the signal MC. The background sample is taken from a two-dimensional $(\Delta E,M_{\text{bc}})$ sideband. For the channel $B^{+}\to h_{c}K^{+}$, the sideband is defined as

[TABLE]

The background sample is divided into training and testing samples of equal size.

The channel $B^{0}\to h_{c}K^{0}_{S}$ has a small number of background events. In order to avoid overtraining, the background region for this channel is redefined. It includes all selected events except the central region defined by Eq. (1). In addition, the MLP internal architecture is changed. Instead of the default tmva neural network with two hidden layers, only one hidden layer is used.

The resulting efficiency of the requirement ($v>v_{0}$) on the MLP output variable $v$ for the training sample is shown in Fig. 3 for the channel $B^{+}\to h_{c}K^{+}$ with $h_{c}\to\eta_{c}(\to K^{+}K^{-}\pi^{0})\gamma$. Note that the efficiency is given by

[TABLE]

where $\epsilon$ is the full efficiency, $\epsilon_{\text{MLP}}$ is the raw MLP output requirement efficiency and $\epsilon_{\text{multiple}}$ is the best-candidate selection efficiency. Because of the correction by $\epsilon_{\text{multiple}}$, the efficiency at the minimal MLP output value $v_{\text{min}}$ is not 1 but rather $\epsilon_{\text{multiple}}(v_{\text{min}})$.

The best-candidate selection is performed for each of the multivariate-analysis channels separately in the following way. The selected $(\Delta E,M_{\text{bc}})$ region is subdivided into three bins in both $\Delta E$ and $M_{\text{bc}}$. The selection is performed for each of the bins separately. The candidate with the largest MLP output is selected. One of the bins ($-67<\Delta E<67\ \text{MeV}$, $5.267<M_{\text{bc}}<5.3\ \text{GeV}/c^{2}$) always contains the entire signal region selected by the optimization procedure as described in Sec. IV.5. Thus, the signal region of the final data sample does not contain multiple candidates that originate from the same $h_{c}$ channel. However, multiple candidates from different $h_{c}$ channels are possible.

The best-candidate selection efficiency increases for larger values of the MLP output cutoff value $v_{0}$. For the $v_{0}$ values obtained as the result of the optimization of the selection requirements as described in Sec. IV.5, the selection procedure removes from $3$ to 15% of data events, depending on the multivariate-analysis channel.

IV.5 Optimization of the selection requirements

Optimization of the selection requirements is performed by maximizing the value

[TABLE]

where $i$ is the channel index, $N_{\text{sig}}^{(i)}$ is the expected number of the signal events for the $i$-th channel, $N_{\text{bg}}^{(i)}$ is the expected number of the background events in the signal region, and $a=3$ is the target significance. This optimization method is based on Ref. Punzi:2003bu .

The signal region is defined as

[TABLE]

where $R_{\Delta E}^{(i)}$ and $R_{M_{\text{bc}}}^{(i)}$ are the half-axes of the signal region ellipse. The parameters determined by the optimization are $R_{\Delta E}^{(i)}$, $R_{M_{\text{bc}}}^{(i)}$, and the minimal value of the MLP output ($v_{0}^{(i)}$) for each channel.

The expected number of signal events for $B^{+}\to h_{c}K^{+}$ is calculated as

[TABLE]

where $N_{\Upsilon(4S)}$ is the number of $\Upsilon(4S)$ events, $\mathcal{B}(h_{c}\to i)$ is the branching fraction of the $h_{c}$ to its $i$-th decay channel, $\epsilon_{\text{SR}}^{(i)}$ is the reconstruction efficiency for the specific signal region SR, and $\epsilon_{S}^{(i)}(v_{0}^{(i)})$ is the efficiency of the requirement ($v>v_{0}$) on the MLP output variable $v$ for the signal events. The number of $\Upsilon(4S)$ events is assumed to be equal to the number of $B\bar{B}$ pairs; the branching fraction $\mathcal{B}(\Upsilon(4S)\to B^{+}B^{-})$ is calculated under the same assumption Tanabashi:2018oca . The signal-region-dependent reconstruction efficiency is calculated as

[TABLE]

where $\epsilon_{R}^{(i)}$ is the reconstruction efficiency, and $S_{i}$ is the signal PDF for $i$-th $h_{c}$ decay channel (the integral of $S_{i}$ over the signal region is the efficiency of the signal region selection). The unknown branching fraction $\mathcal{B}(B^{+}\to h_{c}K^{+})$ can be set to an arbitrary value because the maximum of $F_{\text{opt}}$ does not depend on it. The expected number of signal events for $B^{0}\to h_{c}K^{0}_{S}$ is calculated similarly.

The expected number of background events is calculated as

[TABLE]

where $\epsilon_{B}^{(i)}(v_{0}^{(i)})$ is the efficiency of the MLP output requirement for the background events, $N_{h_{c}\text{ region}}$ is the number of background events in the $h_{c}$ region defined by Eq. (2), $N_{\text{full}}$ is the full number of the background events, and $B_{i}$ is the background density function defined in Eq. (6) for $i$-th $h_{c}$ decay channel.

The optimization is performed separately for two channel groups. The first group includes the multivariate-analysis channels corresponding to the decay $h_{c}\to\eta_{c}\gamma$; the index $i$ runs over all $\eta_{c}$ decay channels. The second group consists of the single channel $h_{c}\to p\bar{p}\pi^{+}\pi^{-}$. The separate optimization is required by the difference of further data processing: the data from the first group of channels are combined into a single $h_{c}\to\eta_{c}\gamma$ data sample, while the $h_{c}\to p\bar{p}\pi^{+}\pi^{-}$ data are fitted with another function, as described below in Sec. V. The optimization results are shown in Table 1. We also check the improvement achieved by MLP usage by changing the selection method to rectangular cuts. The values of $F_{\text{opt}}$ are found to be about 30% and 10% smaller for the channels $h_{c}\to\eta_{c}\gamma$ and $h_{c}\to p\bar{p}\pi^{+}\pi^{-}$, respectively.

After the optimization, the resulting selection criteria are applied. The selected events for the channel $h_{c}\to\eta_{c}\gamma$ are merged. The resolution and distribution in $(\Delta E,M_{\text{bc}})$ are determined again for the $h_{c}\to p\bar{p}\pi^{+}\pi^{-}$ sample, since the knowledge of the background distribution in $(\Delta E,M_{\text{bc}})$ is necessary for the final fit described in Sec. V. The $(\Delta E,M_{\text{bc}})$ fit results are shown in Fig. 4.

IV.6 Resolution in $M_{h_{c}}$

The resolution in $M_{h_{c}}$ is determined from a fit to the combined $h_{c}\to\eta_{c}\gamma$ or $h_{c}\to p\bar{p}\pi^{+}\pi^{-}$ signal MC samples with $\eta_{c}$ decaying to the reconstructed channels only. All final selection criteria are applied. The distribution of the difference of the reconstructed and true masses is fitted to a sum of an asymmetric Gaussian and asymmetric double-sided Crystal Ball functions:

[TABLE]

where $\Delta M$ is the difference of the reconstructed and true $h_{c}$ masses, $N$ is the common normalization and $f_{\text{CB}}$ is the Crystal Ball fraction. Example resolution fit results (for the channel $B^{+}\to h_{c}K^{+}$ with $h_{c}\to\eta_{c}\gamma$) are shown in Fig. 5.

V Fit to the data

V.1 Default model

For the $h_{c}\to\eta_{c}\gamma$ final sample, the distribution in the $h_{c}$ mass in the $(\Delta E,M_{\text{bc}})$ sideband cannot be used to constrain the background level in the signal region because of the presence of peaking backgrounds, such as the background from $B$ decays to a similar final state with a $\pi^{0}$ instead of the $h_{c}$ daughter $\gamma$. If the second photon from this $\pi^{0}$ has a small energy, then the $\Delta E$ and $M_{\text{bc}}$ values are close to 0 and the $B$ mass, respectively. Thus, the fit is based on the signal distribution only for the channel $h_{c}\to\eta_{c}\gamma$.

For the channel $h_{c}\to p\bar{p}\pi^{+}\pi^{-}$, there is a signal from $B$ decays to the same final state ($B^{+}\to p\bar{p}\pi^{+}\pi^{-}K^{+}$ or $B^{0}\to p\bar{p}\pi^{+}\pi^{-}K^{0}_{S}$) that do not proceed via any charmonium state, called the noncharmonium signal hereinafter. Because of the possible interference of the charmonium and noncharmonium signals, the distribution of the noncharmonium signal in the $p\bar{p}\pi^{+}\pi^{-}$ invariant mass needs to be determined by the fit. Thus, both signal and background distributions are included into the fit for the channel $h_{c}\to p\bar{p}\pi^{+}\pi^{-}$.

We perform a simultaneous extended unbinned maximum likelihood fit to the $h_{c}\to\eta_{c}\gamma$ signal, $h_{c}\to p\bar{p}\pi^{+}\pi^{-}$ background, and $h_{c}\to p\bar{p}\pi^{+}\pi^{-}$ signal distributions. The charmonium states are represented by the Breit-Wigner amplitude:

[TABLE]

where $M_{R}$ is the invariant mass, $m_{R}$ is the nominal mass, and $\Gamma_{R}$ is the width of the resonance $R$. The signal-region density function for the channel $h_{c}\to\eta_{c}\gamma$ is given by

[TABLE]

where $N_{h_{c}}$ is the number of signal events, $R_{h_{c}}^{(\eta_{c}\gamma)}$ is the $h_{c}$ mass resolution for the channel $\eta_{c}\gamma$, and $P_{2}$ is a second-order polynomial. The background density function $B_{p\bar{p}\pi^{+}\pi^{-}}(M)$ for the channel $h_{c}\to p\bar{p}\pi^{+}\pi^{-}$ is a third-order polynomial. The signal density function for the channel $h_{c}\to p\bar{p}\pi^{+}\pi^{-}$ is given by

[TABLE]

where $P_{3}$ is a third-order polynomial representing the noncharmonium signal. The wide states are added coherently to the signal density function, while the states that are narrower than the resolution are added incoherently. The amplitudes are normalized in such a way that all the parameters $N_{R}$ represent the yields of the corresponding states. The signal distribution is fitted to the function

[TABLE]

where $w$ is the weight of the background events in the signal region that is calculated as the ratio of integrals of the background distribution in $(\Delta E,M_{\text{bc}})$ over the signal and background regions. The model described above is the default one; additional models are considered to study systematic uncertainties. In the default model, the masses and widths of all resonances are fixed to their world-average values Tanabashi:2018oca ; all other parameters are free.

The best-candidate selection procedure described in Sec. IV.4 guarantees that there are no multiple candidates in the $h_{c}\to p\bar{p}\pi^{+}\pi^{-}$ signal sample, but multiple candidates in the $h_{c}\to\eta_{c}\gamma$ signal sample are possible if they originate from different $\eta_{c}$ decay channels. However, the fraction of the events with multiple candidates is found to be negligibly small. No events with multiple candidates (for the $h_{c}$ masses within the default fitting regions) are observed for both $B^{+}\to h_{c}K^{+}$ and $B^{0}\to h_{c}K^{0}_{S}$ channels.

The fit results are shown in Fig. 6 for the channel $B^{+}\to h_{c}K^{+}$ and in Fig. 7 for the channel $B^{0}\to h_{c}K^{0}_{S}$. The signal yields and phases are listed in Table 2. The statistical significance of the decays $B^{+}\to h_{c}K^{+}$ and $B^{0}\to h_{c}K^{0}_{S}$, as well as the significances of other charmonium states in the channel $p\bar{p}\pi^{+}\pi^{-}$ are calculated from the difference of $(-2\ln L)$, where $L$ is the maximum likelihood, between the models with and without these states taking the number of degrees of freedom into account. The significance of the decays $B^{+}\to h_{c}K^{+}$ and $B^{0}\to h_{c}K^{0}_{S}$ in the default model is found to be $5.0\sigma$ and $0.8\sigma$, respectively. The significance of the decays $B^{+}\to h_{c}K^{+}$ and $B^{0}\to h_{c}K^{0}_{S}$ with the systematic error taken into account is $4.8\sigma$ and $0.7\sigma$, respectively; the procedure of the calculation of the systematic uncertainty is described in Sec. V.2. Thus, we find evidence for the decay $B^{+}\to h_{c}K^{+}$, but do not find evidence for $B^{0}\to h_{c}K^{0}_{S}$. The significances of charmonium states in the channel $p\bar{p}\pi^{+}\pi^{-}$ (except the $h_{c}$, which is reconstructed in two decay channels) are shown in Table 3. The significance of the $\eta_{c}(2S)$ in the default model is $12.3\sigma$ and $5.9\sigma$ for the processes $B^{+}\to(c\bar{c})(\to p\bar{p}\pi^{+}\pi^{-})K^{+}$ and $B^{0}\to(c\bar{c})(\to p\bar{p}\pi^{+}\pi^{-})K^{0}_{S}$, respectively. The significance including the systematic error is $12.1\sigma$ and $5.8\sigma$, respectively. Consequently, the decay $\eta_{c}(2S)\to p\bar{p}\pi^{+}\pi^{-}$ is observed for the first time in both $B^{+}\to(c\bar{c})(\to p\bar{p}\pi^{+}\pi^{-})K^{+}$ and $B^{0}\to(c\bar{c})(\to p\bar{p}\pi^{+}\pi^{-})K^{0}_{S}$ processes.

V.2 Systematic uncertainty: model dependence

For a systematic-uncertainty study, we consider additional models. They include the models with free masses and widths of the $h_{c}$ and all other charmonium states (with Gaussian constraints in accordance with the errors of their current world average values), with increased order (3) of the background PDF polynomial, different fitting ranges, scaled resolution, and variation of the relative fraction of the channels $h_{c}\to\eta_{c}\gamma$ and $h_{c}\to p\bar{p}\pi^{+}\pi^{-}$. For the model with scaled resolution, the resolution function $R_{h_{c}}(\Delta M)$ is changed to

[TABLE]

where $\mathcal{S}$ is the resolution scaling parameter. The variation of the relative fraction of the channels $h_{c}\to\eta_{c}\gamma$ and $h_{c}\to p\bar{p}\pi^{+}\pi^{-}$ is performed by changing the expected yields in these channels by $\pm 1\sigma$, where the error is due to the error of the corresponding branching fractions. The results are listed in Table 4.

V.3 Branching fraction

Using the number of reconstructed events, we calculate the branching fractions of the decays $B^{+}\to h_{c}K^{+}$ and $B^{0}\to h_{c}K^{0}_{S}$, as well as the branching fraction products $\mathcal{B}(B^{+}\to(c\bar{c})K^{+})\times\mathcal{B}((c\bar{c})\to p\bar{p}\pi^{+}\pi^{-})$ and $\mathcal{B}(B^{0}\to(c\bar{c})K^{0}_{S})\times\mathcal{B}((c\bar{c})\to p\bar{p}\pi^{+}\pi^{-})$. The decay $\psi(2S)\to p\bar{p}\pi^{+}\pi^{-}$ can proceed via the $J/\psi$: $\psi(2S)\to J/\psi(\to p\bar{p})\pi^{+}\pi^{-}$. To remove the events with a $J/\psi$, the $\psi(2S)$ yield is taken from an alternative fit with an additional $J/\psi$ veto defined as $|M_{p\bar{p}}-m_{J/\psi}|>50\ \text{MeV}/c^{2}$.

The sources of the systematic uncertainty of the branching fractions include the model dependence (the same set of alternative fit models is used as in Sec. V.2), overtraining (the difference between the efficiency in the training and testing samples), the error of the difference of the particle identification requirements efficiency between the data and MC, the difference of the MLP efficiency between the data and MC, tracking efficiency, number of $\Upsilon(4S)$ events, the $\eta_{c}\to K^{+}K^{0}_{S}\pi^{-}$ and $\Upsilon(4S)\to B^{+}B^{-}$ or $\Upsilon(4S)\to B^{0}\bar{B}^{0}$ branching fractions. All systematic error sources are listed in Table 5 for the channel $B^{+}\to(c\bar{c})K^{+}$ ($B^{+}\to h_{c}K^{+}$ for the $h_{c}$ and $B^{+}\to(c\bar{c})(\to p\bar{p}\pi^{+}\pi^{-})K^{+}$ for all other charmonium states) and in Table 6 for the channel $B^{0}\to(c\bar{c})K^{0}_{S}$ ($B^{0}\to h_{c}K^{0}_{S}$ for the $h_{c}$ and $B^{0}\to(c\bar{c})(\to p\bar{p}\pi^{+}\pi^{-})K^{0}_{S}$ for all other charmonium states). The errors of the tracking efficiency and the difference of the efficiency of the particle identification requirements depend on the multivariate-analysis channel in case of the calculation of $\mathcal{B}(B^{+}\to h_{c}K^{+})$ and $\mathcal{B}(B^{0}\to h_{c}K^{0}_{S})$; the errors related to the MLP efficiency and overtraining are estimated separately for $h_{c}\to\eta_{c}\gamma$ and $h_{c}\to p\bar{p}\pi^{+}\pi^{-}$. The values presented in Tables 5 and 6 are weighted averages.

The difference of the particle identification requirements efficiency between the data and MC is estimated from several control samples, such as $D^{*+}\to D^{0}(\to K^{-}\pi^{+})\pi^{+}$ for $K$ and $\pi$, $\Lambda\to p\pi^{+}$ for $p$, $\Lambda_{c}^{+}\to\Lambda\pi^{+}$ for $\Lambda$, $D^{0}\to K^{0}_{S}\pi^{+}\pi^{-}$ for $K^{0}_{S}$, and $\tau^{-}\to\pi^{-}\pi^{0}\nu_{\tau}$ for $\pi^{0}$. The resulting overall efficiency ratio depends on the final state and the momenta of the decay products; thus, it is different for all branching fractions. For example, for the $h_{c}$ it is found to be $(95.0\pm 3.8)\%$ for the channel $B^{+}\to h_{c}K^{+}$ and $(93.0\pm 4.5)\%$ for the channel $B^{0}\to h_{c}K^{0}_{S}$.

The error caused by the difference of the MLP efficiency between the data and MC is estimated for the channel $h_{c}\to\eta_{c}\gamma$ using the decay mode $B^{0}\to\eta_{c}\pi^{-}K^{+}$. This decay is reconstructed using selection criteria that are as similar as possible to the signal mode $B^{+}\to h_{c}K^{+}$. The same MLP optimized for $B^{+}\to h_{c}K^{+}$ is applied to the control channel. Some MLP input variables used for the signal channel are undefined for $B^{0}\to\eta_{c}\pi^{-}K^{+}$, for example, the number of $\pi^{0}$ candidates that include the $h_{c}$ daughter $\gamma$ as one of their daughters. Such variables are held constant. The ratio of the number of signal candidates before and after the application of the MLP selection requirements is used to measure the difference of the MLP efficiency between the data and MC:

[TABLE]

Only the channel $\eta_{c}\to K^{+}K^{0}_{S}\pi^{-}$ is used before the MLP selection, because only this channel is sufficiently clean for the determination of the number of the signal events without the MLP selection. The ratio $r_{\text{MLP}}$ is extracted from a simultaneous fit to the $\eta_{c}$ mass distribution before and after the application of the MLP selection. The relative difference between the values of $r_{\text{MLP}}$ in data and MC is found to be $(14.4\pm 9.0)\%$. For conservative treatment, the statistical error is added in quadrature to the central value of the difference. The resulting systematic uncertainty caused by the MLP selection efficiency difference in data and MC is 17.0%.

The estimation of the MLP efficiency error for the hadronic channel $p\bar{p}\pi^{+}\pi^{-}$ is done by performing the fit using the $p\bar{p}\pi^{+}\pi^{-}$ data without the MLP selection and comparing the resulting branching fraction products $\mathcal{B}(B^{+}\to J/\psi K^{+})\times\mathcal{B}(J/\psi\to p\bar{p}\pi^{+}\pi^{-})$ with the results of the default procedure. Their relative difference is found to be $0.2\%$.

The same estimates of the MLP efficiency uncertainty for the channels $h_{c}\to\eta_{c}\gamma$ and $h_{c}\to p\bar{p}\pi^{+}\pi^{-}$ are used for both $B^{+}\to h_{c}K^{+}$ and $B^{0}\to h_{c}K^{0}_{S}$. The final value of the MLP efficiency uncertainty is calculated as a weighted average of the errors for the two $h_{c}$ decay channels. The result is slightly different for the channels $B^{+}\to h_{c}K^{+}$ and $B^{0}\to h_{c}K^{0}_{S}$ because of the difference in the relative number of the expected $h_{c}\to\eta_{c}\gamma$ and $h_{c}\to p\bar{p}\pi^{+}\pi^{-}$ events. The MLP efficiency error for the branching fraction products for the channels $B^{+}\to(c\bar{c})(\to p\bar{p}\pi^{+}\pi^{-})K^{+}$ and $B^{0}\to(c\bar{c})(\to p\bar{p}\pi^{+}\pi^{-})K^{0}_{S}$ is equal to the error for the channel $p\bar{p}\pi^{+}\pi^{-}$, since all charmonium states other than the $h_{c}$ are reconstructed in this channel only.

The MLP efficiency uncertainty for the channel $h_{c}\to\eta_{c}\gamma$ does not include the uncertainty caused by the difference between the data and MC in the distributions of the variables that are not defined for the channel $B^{0}\to\eta_{c}\pi^{-}K^{+}$. There are four such variables: the $\eta_{c}$ mass, the $h_{c}$ helicity angle, and two numbers of $\pi^{0}$ candidates that include the $h_{c}$ daughter photon as one of their daughters. The distribution of the $h_{c}$ helicity angle for the signal events is known precisely; thus, there is no additional uncertainty caused by the difference of its distribution in data and MC. The difference of the numbers of $\pi^{0}$ candidates is taken into account by removing these variables from the neural network for the channel $B^{+}\to h_{c}K^{+}$, performing an alternative optimization, and comparing the resulting $h_{c}$ branching fractions in the channel $h_{c}\to\eta_{c}\gamma$. The relative difference is found to be $15.6\%$. The error due to the $\eta_{c}$ mass distribution uncertainty is estimated by varying the $\eta_{c}$ mass and width by $\pm 1\sigma$ and reweighting the selected MC events in accordance with the relative difference between the modified and default $\eta_{c}$ mass distributions. The largest resulting efficiency difference is considered as the systematic uncertainty related to the $\eta_{c}$ mass distribution. This uncertainty is estimated to be $1.3\%$ for both $B^{+}\to h_{c}K^{+}$ and $B^{0}\to h_{c}K^{0}_{S}$ channels.

The ratio $r_{\text{MLP}}$ used for determination of the MLP efficiency uncertainty includes the number of reconstructed events relatively to the number of events in $h_{c}\to\eta_{c}(\to K^{+}K^{0}_{S}\pi^{-})\gamma$. The total expected number of events can be calculated as

[TABLE]

where $\epsilon_{i}$ and $\mathcal{B}_{i}$ are the efficiency and branching fraction for $i$-th channel, respectively. The last term in Eq. (21) is proportional to $r_{\text{MLP}}$. Consequently, for the channel $h_{c}\to\eta_{c}\gamma$ one needs to take into account only the error of $\mathcal{B}(h_{c}\to\eta_{c}(\to K^{+}K^{0}_{S}\pi^{-})\gamma)$. Errors of the branching fractions of all other channels relatively to $h_{c}\to\eta_{c}(\to K^{+}K^{0}_{S}\pi^{-})\gamma$ enter the MLP efficiency error.

The final $h_{c}$ branching fraction error is calculated as a weighted average of the errors for the channels $h_{c}\to\eta_{c}\gamma$ and $h_{c}\to p\bar{p}\pi^{+}\pi^{-}$. Since branching fraction products are measured for all other charmonium states, they do not have a similar systematic error source.

The resulting branching fractions with both statistical and systematic errors are listed in Table 7. For insignificant decays or decay chains, the confidence intervals are calculated in the frequentist approach Feldman:1997qc using an asymmetric Gaussian as the branching-fraction PDF and the measured central value and errors as its parameters. This is the first measurement of $\mathcal{B}(B^{0}\to h_{c}K^{0}_{S})$. Also, the branching fraction products $\mathcal{B}(B^{+}\to(c\bar{c})K^{+})\times\mathcal{B}((c\bar{c})\to p\bar{p}\pi^{+}\pi^{-})$ and $\mathcal{B}(B^{0}\to(c\bar{c})K^{0}_{S})\times\mathcal{B}((c\bar{c})\to p\bar{p}\pi^{+}\pi^{-})$ are measured directly in $B$ decays for the first time. The current world-average values of the same branching fractions Tanabashi:2018oca are also presented in Table 7 for comparison if they are known. The values of the branching fraction products are calculated by multiplying the individual branching fractions listed in Ref. Tanabashi:2018oca assuming uncorrelated errors. The measured branching fractions are consistent with the world averages; the largest deviation is observed for the branching-fraction product $\mathcal{B}(B^{0}\to\chi_{c1}K^{0}_{S})\times\mathcal{B}(\chi_{c1}\to p\bar{p}\pi^{+}\pi^{-})$, which differs from the world-average value by $2.4\sigma$ taking its error into account. The results for $\mathcal{B}(B^{+}\to\eta_{c}K^{+})\times\mathcal{B}(\eta_{c}\to p\bar{p}\pi^{+}\pi^{-})$ and $\mathcal{B}(B^{0}\to\eta_{c}K^{0}_{S})\times\mathcal{B}(\eta_{c}\to p\bar{p}\pi^{+}\pi^{-})$ have a better precision than the world-average values.

VI Conclusions

A search for the decays $B^{+}\to h_{c}K^{+}$ and $B^{0}\to h_{c}K^{0}_{S}$ has been performed. Evidence for the decay $B^{+}\to h_{c}K^{+}$ is found; its significance is $4.8\sigma$. No evidence is found for $B^{0}\to h_{c}K^{0}_{S}$. The branching fraction of $B^{+}\to h_{c}K^{+}$ is measured to be $(3.7\tensor{*}[^{+1.0}_{-0.9}]{}{}\tensor{*}[^{+0.8}_{-0.8}]{}{})\times 10^{-5}$; the upper limit for the $B^{0}\to h_{c}K^{0}_{S}$ branching fraction is $1.4\times 10^{-5}$ at $90\%$ C. L. The measured value of $\mathcal{B}(B^{+}\to h_{c}K^{+})$ is consistent with the existing upper limit of $3.8\times 10^{-5}$ ($90\%$ C. L.) obtained in the previous Belle analysis Fang:2006bz and supersedes it. The resulting branching fraction $\mathcal{B}(B^{+}\to h_{c}K^{+})$ agrees with the existing theoretical predictions Meng:2006mi ; Li:2006vj ; Beneke:2008pi . In addition, a study of the $p\bar{p}\pi^{+}\pi^{-}$ invariant mass distribution in the channel $B^{+}\to(p\bar{p}\pi^{+}\pi^{-})K^{+}$ results in the first observation of the decay $\eta_{c}(2S)\to p\bar{p}\pi^{+}\pi^{-}$ with $12.1\sigma$ significance.

Acknowledgement

We thank the KEKB group for the excellent operation of the accelerator; the KEK cryogenics group for the efficient operation of the solenoid; and the KEK computer group, and the Pacific Northwest National Laboratory (PNNL) Environmental Molecular Sciences Laboratory (EMSL) computing group for strong computing support; and the National Institute of Informatics, and Science Information NETwork 5 (SINET5) for valuable network support. We acknowledge support from the Ministry of Education, Culture, Sports, Science, and Technology (MEXT) of Japan, the Japan Society for the Promotion of Science (JSPS), and the Tau-Lepton Physics Research Center of Nagoya University; the Australian Research Council including grants DP180102629, DP170102389, DP170102204, DP150103061, FT130100303; Austrian Science Fund (FWF); the National Natural Science Foundation of China under Contracts No. 11435013, No. 11475187, No. 11521505, No. 11575017, No. 11675166, No. 11705209; Key Research Program of Frontier Sciences, Chinese Academy of Sciences (CAS), Grant No. QYZDJ-SSW-SLH011; the CAS Center for Excellence in Particle Physics (CCEPP); the Shanghai Pujiang Program under Grant No. 18PJ1401000; the Ministry of Education, Youth and Sports of the Czech Republic under Contract No. LTT17020; the Carl Zeiss Foundation, the Deutsche Forschungsgemeinschaft, the Excellence Cluster Universe, and the VolkswagenStiftung; the Department of Science and Technology of India; the Istituto Nazionale di Fisica Nucleare of Italy; National Research Foundation (NRF) of Korea Grants No. 2015H1A2A1033649, No. 2016R1D1A1B01010135, No. 2016K1A3A7A09005 603, No. 2016R1D1A1B02012900, No. 2018R1A2B3003 643, No. 2018R1A6A1A06024970, No. 2018R1D1 A1B07047294; Radiation Science Research Institute, Foreign Large-size Research Facility Application Supporting project, the Global Science Experimental Data Hub Center of the Korea Institute of Science and Technology Information and KREONET/GLORIAD; the Polish Ministry of Science and Higher Education and the National Science Center; the Russian Foundation for Basic Research Grant No. 18-32-00277; the Slovenian Research Agency; Ikerbasque, Basque Foundation for Science, Spain; the Swiss National Science Foundation; the Ministry of Education and the Ministry of Science and Technology of Taiwan; and the United States Department of Energy and the National Science Foundation.