Measurement of branching fractions of $\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}K_{S}^{0}$ and $\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}\eta$ at Belle

The weak decays of charmed baryons provide an excellent platform for understanding Quantum Chromodynamics with transitions involving the charm quark. The decay amplitudes consist of factorizable and non-factorizable contributions. The latter may play a non-trivial or essential role and are approached in various ways, including the pole model [1, 2], the covariant confined quark model [3, 4], current algebra [5, 6, 7] and SU(3)_F symmetry [8, 9, 10]. To date, there is no established phenomenological model that consistently describes baryon decays. Precise measurements of branching fractions of charmed baryon weak decays are useful for studying the dynamics of charmed baryons and testing the predictions of theoretical models. In addition, the singly Cabibbo-suppressed (SCS) charm decays are essential probes of $C\!P$ violation in the charm sector [11, 12, 13] and new physics beyond the standard model [14, 15, 16].

Experimentally, the investigation of charmed baryons is more challenging than that of charmed mesons, mainly due to lower production rates. For the lightest state, $\mathchar 28931\relax_{c}^{+}$, hadronic modes have been studied at several experiments, but some have yet to be observed or are measured with low precision [17]. For the Cabibbo-favored (CF) channel ${\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}\eta}$ ¹¹1Throughout this paper charge-conjugate modes are implied., the world average branching fraction, $\mathcal{B}({\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}\eta})\!=\!{(4.15\pm 0.90)\times 10^{-3}}$ [17], still has a large uncertainty ($22\%$). The SCS mode ${\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}K_{S}^{0}}$, for which the predicted branching fraction is $\mathcal{B}({\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}K_{S}^{0}})\!=\!{(1.9\pm 0.4)\times 10^{-3}}$ based on SU(3)_F symmetry [19], has not previously been observed.

In this paper, we present a precise measurement of $\mathcal{B}({\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}K_{S}^{0}})$ and $\mathcal{B}({\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}\eta})$ based on the full Belle data set. For both of these three-body decays, the Dalitz plot is of interest for the study intermediate resonances. Understanding the nature of $N^{*}(1535)$ is very challenging and important for hadronic physics. The mass of $N^{*}(1535)$, with spin parity $J^{P}\!=\!1/2^{-}$, is larger than that of the radial excitation $N^{*}(1440)$, in opposition to predictions of classical constituent quark models [20]. The $N^{*}(1535)$ also couples strongly to channels with strangeness, such as $\eta{}N$ and $K\mathchar 28931\relax$, which is difficult to explain within the naive constituent quark models [21, 22]. The inclusion of five-quark components gives a natural explanation for these properties [23]. The ${\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}\eta}$ decay, in which the final-state $p\eta$ is in a pure isospin $I\!=\!1/2$ state, is an ideal process for studying the $N^{*}(1535)$ resonance, as $N^{*}(1535)$ has a large branching ratio to $p\eta$, in $S$-wave. Other intermediate resonances of interest are the light scalars $a_{0}(980)$ and $f_{0}(980)$, which both couple to $K\overline{K}$ in ${\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}K_{S}^{0}}$. They contribute to the SCS $\mathchar 28931\relax_{c}^{+}$ decays $\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK\overline{K}$ and $\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}p\pi\pi$, as predicted in Ref. [24], and likely contribute to ${\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}K_{S}^{0}}$, based on isospin symmetry. The nature of $f_{0}(980)$ and $a_{0}(980)$ remains poorly understood and continues to be controversial [25, 26, 27]. They are often interpreted as compact tetraquark states [28, 29, 30] or $K\overline{K}$ bound states [31, 32]. Therefore, we reconstruct the Dalitz plots of ${\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}K_{S}^{0}}$ and ${\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}\eta}$ decays to check such interesting intermediate resonances.

This analysis uses the full dataset recorded by the Belle detector [33] operating at the KEKB energy-asymmetric $e^{+}e^{-}$ collider [34]. This data sample corresponds to a total integrated luminosity of $980~{}\mbox{\,fb}^{-1}$ collected at or near the $\Upsilon(nS)$ ($n\!=\!1,\,2,\,3,\,4,\,5$) resonances. The Belle detector is a large-solid-angle magnetic spectrometer consisting of a silicon vertex detector, a 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) consisting of CsI(Tl) crystals. These components are all located inside a superconducting solenoid coil that provides a 1.5 T magnetic field. The iron flux-return of the magnet is instrumented to detect $K^{0}_{L}$ mesons and to identify muons. The detector is described in detail elsewhere [33].

Monte Carlo (MC) simulated events are generated with evtgen [35] and pythia [36], and are subsequently processed through the full detector simulation based on geant3 [37]. Final-state radiation from charged particles is included at the event generation stage using photos [38]. “Generic” MC samples include $B\overline{B}$ events and continuum processes $e^{+}e^{-}\!\rightarrow\!{}q\bar{q}$ ($q\!=\!u,\,d,\,s,\,c)$ corresponding to an integrated luminosity three times that of the data. Samples of MC events of $\mathchar 28931\relax_{c}^{+}$ signal decay modes are produced in the $e^{+}e^{-}\!\rightarrow\!{}c\overline{c}$ process, decayed uniformly in three-body phase space, and used to study the efficiency.

We reconstruct the two signal modes ${\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}K_{S}^{0}}$ and ${\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}\eta}$ and their reference mode ${\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}}$. The event selections are optimized based on a figure of merit (FOM), defined as ${\rm FOM}\!=\!\varepsilon_{\scriptscriptstyle{S}}/\sqrt{N_{B}}$ for ${\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}K_{S}^{0}}$ due to its branching fraction having not yet been measured, and ${\rm FOM}\!=\!N_{S}/\sqrt{N_{S}+N_{B}}$ for ${\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}\eta}$ assuming its current world average branching fraction [17]. Here $\varepsilon_{\scriptscriptstyle{S}}$ is the selection efficiency of signal, $N_{S}$ and $N_{B}$ are the expected yields of signal and background, respectively, based on numbers of candidates in the $M(\mathchar 28931\relax_{c}^{+})$ signal regions, where $M(\mathchar 28931\relax_{c}^{+})$ is the invariant mass of reconstructed $\mathchar 28931\relax_{c}^{+}$ candidates. These signal regions are defined to be within 10, 22, and 18 MeV/$c^{2}$ of the nominal $\mathchar 28931\relax_{c}^{+}$ mass [17] for the ${\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}K_{S}^{0}}$, ${\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}\eta}$, and ${\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}}$ channels, respectively; each signal band includes $\approx$98% of the signal. For the expected background, $N_{B}$, the number found in MC is multiplied by the data/MC yield ratio in the $M(\mathchar 28931\relax_{c}^{+})$ sideband region ($30<|M(\mathchar 28931\relax_{c}^{+})-m_{\mathchar 28931\relax_{c}^{+}}|\!<\!{50~{}{\rm MeV}/c^{2}}$), where $m_{\mathchar 28931\relax_{c}^{+}}$ is the nominal $\mathchar 28931\relax_{c}^{+}$ mass [17].

The particle identification (PID) likelihood for a given particle hypothesis, $\mathcal{L}_{i}$ ($i\!=\!\pi,\,K,\,p$), is calculated from the Cherenkov photon yield in the ACC, energy-loss measurements in the CDC, and time-of-flight information from the TOF [39]. Charged tracks satisfying $\mathcal{R}(p|K)\!=\!\mathcal{L}_{p}/(\mathcal{L}_{p}+\mathcal{L}_{K})\!>\!0.6$ and $\mathcal{R}(p|\pi)\!=\!\mathcal{L}_{p}/(\mathcal{L}_{p}+\mathcal{L}_{\pi})\!>\!0.6$, are identified as protons. These PID requirements have signal efficiencies of 94% for ${\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}K_{S}^{0}}$ and 97% for ${\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}\eta}$.

For proton candidates, the point on the track nearest to the axis defined by the positron beam and in the direction opposite to it (“$z$-axis”) is required to be within 3.0 cm of the interaction point in the $z$-direction and within 1.0 cm on the transverse ($x$-$y$) plane. This requirement rejects tracks not originating at the interaction point (IP) and introduces a negligible signal efficiency loss ($<0.01\%$).

Candidate $K_{S}^{0}$’s are reconstructed from pairs of oppositely-charged tracks, treated as pions, using an artificial neural network (NN) [40]. The NN utilizes the following 13 input variables: the $K_{S}^{0}$ momentum in the laboratory frame; the separation in $z$ between the two $\pi^{\pm}$ tracks at their intersection in the $x$-$y$ plane; for each track, the nearest distance to the IP in the $x$-$y$ plane; the $K_{S}^{0}$ flight length in the $x$-$y$ plane; the angle between the $K_{S}^{0}$ momentum and the vector joining the IP to the $K_{S}^{0}$ decay vertex; in the $K_{S}^{0}$ rest frame, the angle between the $\pi^{+}$ momentum and the laboratory-frame boost direction; and, for each $\pi^{\pm}$ track, the number of CDC hits in both stereo and axial views, and the presence or absence of SVD hits. Detailed information is provided elsewhere [41]. The invariant mass of the reconstructed $K_{S}^{0}\!\rightarrow\!\pi^{+}\pi^{-}$ candidate is required to lie within $10~{}{\rm MeV}/c^{2}$ of the nominal $K_{S}^{0}$ mass [17]; this includes 99.9% of the $K_{S}^{0}$ signal. The two pion tracks from each $K_{S}^{0}$ candidate are refitted to originate from a common vertex and constrained to have invariant mass equal to the nominal $K_{S}^{0}$ mass [17]. The corresponding fit quality $\chi_{\rm mv}^{2}(K_{S}^{0})$ is required to be smaller than 100. The selected $K_{S}^{0}$ sample has a purity of greater than 98%.

Photon candidates are identified as energy clusters in the ECL that are not associated with any charged track. The ratio of the energy deposited in the 3$\times$3 array of crystals centered on the crystal with the highest energy, to the energy deposited in the corresponding 5$\times$5 array of crystals, is required to be greater than 0.8. The photon energy is required to be greater than 50 MeV in the barrel region (covering the polar angle $32^{\circ}\!<\!\theta\!<\!129^{\circ}$), and greater than 100 MeV in the endcap region ($12^{\circ}\!<\!\theta\!<\!31^{\circ}$ or $132^{\circ}\!<\!\theta\!<\!157^{\circ}$).

Candidate $\eta\!\rightarrow\!\gamma\gamma$ decays are reconstructed from photon pairs having an invariant mass satisfying $500~{}{\rm MeV}/c^{2}\!<\!M(\gamma\gamma)\!<\!580~{}{\rm MeV}/c^{2}$ ($3\sigma$ in $M_{\eta}(\gamma\gamma)$ resolution). The invariant mass of each $\eta$ candidate is constrained to the nominal $\eta$ mass [17] at the $\mathchar 28931\relax_{c}^{+}$ decay vertex (described below). The fit quality of this mass constraint is required to satisfy $\chi_{m}^{2}(\eta)\!<\!8$, and the resulting $\eta$ momentum in the laboratory frame is required to be greater than 0.4 GeV/$c$. To further suppress the background, $\eta$ candidates are vetoed if either of daughters can be paired with another photon such that the $\gamma\gamma$ pair has an invariant mass within $2.5\sigma$ of the nominal $\pi^{0}$ mass ($\sigma\!=\!5$ MeV/$c^{2}$). This $\pi^{0}$-veto results in a signal loss of 28% and removes 72% of background.

The $\mathchar 28931\relax_{c}^{+}$ candidates are assembled by forming combinations of the final-state particles for each mode. The $p$ and $K_{S}^{0}$ are required to originate from a common vertex (denoted the $\mathchar 28931\relax_{c}^{+}$ decay vertex and the $K_{S}^{0}$ production vertex) with a fit quality $\chi_{\rm vtx}^{2}\!<\!24$. To reduce combinatorial background, the scaled momentum of the $\mathchar 28931\relax_{c}^{+}$ candidate, defined as $x_{p}\!=\!p^{*}c/\sqrt{s/4-M^{2}(\mathchar 28931\relax_{c}^{+})\cdot c^{4}}$, is required to be greater than 0.48, where $s$ is the square of the center-of-mass energy and $p^{*}$ is the momentum of reconstructed $\mathchar 28931\relax_{c}^{+}$ candidates in the $e^{+}e^{-}$ center-of-mass frame.

For the SCS decay ${\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}K_{S}^{0}}$, a non-$K_{S}^{0}$ peaking background from the CF decay $\mathchar 28931\relax_{c}^{+}\!\rightarrow\!pK_{S}^{0}\pi^{+}\pi^{-}$ exists, even though it is suppressed by the vertex fit and $K_{S}^{0}$ selection. The $K_{S}^{0}$ decay length $L$ is determined by the projection of the vector joining the $K_{S}^{0}$ production and decay vertices onto the $K_{S}^{0}$ momentum direction, and its corresponding uncertainty $\sigma_{L}$ is calculated by propagating uncertainties in the vertices and the $K_{S}^{0}$ momentum, including their correlations. To suppress the non-$K_{S}^{0}$ peaking CF background, we require the significance of the $K_{S}^{0}$ decay length $L/\sigma_{L}(K_{S}^{0})\!>\!10$ for the slower of the two $K_{S}^{0}$’s in ${\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}K_{S}^{0}}$. This requirement reduces the signal efficiency by 3%, and rejects 80% of non-$K_{S}^{0}$ peaking background. The remaining non-$K_{S}^{0}$ peaking background is ignored in the $M(\mathchar 28931\relax_{c}^{+})$ fits because it has a tiny ratio 0.4% to signal based on the MC studies with the branching fraction $(1.6\pm 0.12)\%$ [17], but considered in the systematic uncertainty.

After applying all selection criteria to the data, we find 1.03, 1.06, and 1.01 candidates per event for ${\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}K_{S}^{0}}$, ${\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}\eta}$, and ${\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}}$, respectively, in candidates selected from the entire $M(\mathchar 28931\relax_{c}^{+})$ fit region ($|M(\mathchar 28931\relax_{c}^{+})-m_{\mathchar 28931\relax_{c}^{+}}|\!<\!0.05~{}{\rm GeV}/c^{2}$). Correspondingly, about 3.1%, 5.7% and 1.2% of events have multiple signal candidates, which do not introduce any peaking background. We retain all candidates for this branching fraction measurement.

The signal yield is extracted by an unbinned extended maximum likelihood fit to the $M(\mathchar 28931\relax_{c}^{+})$ distribution. The signal probability density function (PDF) is a sum of three symmetric Gaussian functions for the ${\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}}K_{S}^{0}$ mode, a sum of one symmetric Gaussian and two asymmetric Gaussians for the ${\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}\eta}$ mode, and a sum of one symmetric Gaussian and three asymmetric Gaussians for the ${\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}}$ mode. The Gaussian functions share a common mean parameter but have different width parameters. The fit is first performed on truth-matched signal MC events.

In fitting data, the mean is allowed a common shift ($\delta_{\mu}$) from the value found in MC, and the widths are those found in MC, multiplied by a common scaling factor ($k_{\sigma}$). The background PDF is a first-order polynomial function for ${\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}K_{S}^{0}}$ and a second-order polynomial function for ${\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}\eta}$ and ${\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}}$. The background parameters are floated to account for differences between the experimental data and MC simulated samples. The results are shown in Fig. 1, along with the pulls $(N_{\rm data}-N_{\rm fit})/\sigma_{\rm data}$ where $\sigma_{\rm data}$ is the error on $N_{\rm data}$. The pull distributions demonstrate that the data are statistically consistent with the fitted shapes. The signal and background yields are listed in Table 1.

For the ${\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}K_{S}^{0}}$ mode, we obtain the difference in the log likelihoods obtained from fits performed with and without a signal PDF, ${\Delta\ln\mathcal{L}\!=\!524}$; as the number of degrees of freedom without a signal component is three less than that with a signal component (parameters $N_{\rm sig}$, $\delta_{\mu}$ and $k_{\sigma}$ are dropped), and this value of ${\Delta\ln\mathcal{L}}$ corresponds to a statistical significance greater than $10\sigma$. This measurement constitutes the first observation of this SCS $\mathchar 28931\relax_{c}^{+}$ decay.

For the three-body decay modes, the Dalitz plots for candidates in the $M(\mathchar 28931\relax_{c}^{+})$ signal region and sideband region are shown in Figs. 2(a, b) for ${\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}K_{S}^{0}}$ and Figs. 2(d, e) for ${\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}\eta}$. For ${\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}K_{S}^{0}}$, Bose symmetry requires invariance under the exchange of the two $K_{S}^{0}$’s, hence the Dalitz plot for two $pK_{S}^{0}$ masses is symmetric. We plot $M^{2}(pK_{S}^{0})_{\rm max}$ versus $M^{2}(pK_{S}^{0})_{\rm min}$ in half of the Dalitz plot, as shown in Figs. 2(a–c), and use it to measure the branching fraction.

For each mode, a large MC sample of signal events, generated uniformly across the decay phase space, is used to determine the reconstruction efficiency. For ${\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}K_{S}^{0}}$ and ${\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}\eta}$, the efficiencies are calculated in bins across the phase space, based on truth-matched signal yield in the $M(\mathchar 28931\relax_{c}^{+})$ signal region. The results are shown in Fig. 2(c) for ${\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}K_{S}^{0}}$ and Fig. 2(f) for ${\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}\eta}$.

In order to calculate the efficiency-corrected yield, properly taking into account the variations in efficiency and uncertainties in signal yield over the Dalitz plot, we make a bin-by-bin correction. The Dalitz plots are divided uniformly into 7$\times$7 bins for ${\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}K_{S}^{0}}$ and 5$\times$5 bins for ${\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}\eta}$, as shown in Figs. 2(c, f) respectively. The efficiency-corrected yields are

where $N_{i}^{\rm{tot}}$ is the raw yield in the $i^{\rm th}$ bin of the Dalitz plot in $M(\mathchar 28931\relax_{c}^{+})$ signal region, $N_{\rm bkg}^{\rm SR}$ is the fitted background yield as listed in Table 1, $f_{i}^{\rm{bkg}}$ is the fraction of background in the $i^{\rm th}$-bin, with $\sum_{i}{}f_{i}\!=\!1$. These fractions are obtained from the Dalitz plot distribution of events in the $M(\mathchar 28931\relax_{c}^{+})$ sideband region, shown in Fig. 2(b) for ${\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}K_{S}^{0}}$ and Fig. 2(e) for ${\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}\eta}$. Using the generic MC sample, we find that the Dalitz plot in the chosen $M(\mathchar 28931\relax_{c}^{+})$ sideband region is consistent with the generic background in the $M(\mathchar 28931\relax_{c}^{+})$ signal region. The uncertainties on each variable in Eq. (1) have been considered and are propagated into the efficiency-corrected yields, $N_{\rm corr}$. We obtain

The relative branching fractions of signal modes to reference mode are determined by Eqs. (4, 5).

Here, $\varepsilon_{0}\!=\!(33.09\pm 0.05)\%$ is the efficiency of the reference mode ${\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}}$ in the $M(\mathchar 28931\relax_{c}^{+})$ signal region. Inserting the efficiency-corrected yields in Eqs. (2, 3), $N^{\rm SR}_{\rm{sig}}({\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}})$ in Table 1, and the world averages $\mathcal{B}(K_{S}^{0}\!\rightarrow\!\pi^{+}\pi^{-})\!=\!(69.20\pm 0.05)\%$ and $\mathcal{B}(\eta\!\rightarrow\!\gamma\gamma)\!=\!(39.41\pm 0.20)\%$ [17], we find

Combining with the world average branching fraction of reference mode $\mathcal{B}({\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}})\!=\!(1.59\pm 0.08)\%$ [17], we have the absolute branching fractions:

where the uncertainties are statistical and from the uncertainty on $\mathcal{B}({\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}})$.

We examine the Dalitz plots for ${\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}K_{S}^{0}}$ and ${\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}\eta}$, after background subtraction and efficiency correction, for intermediate resonances. In ${\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}K_{S}^{0}}$, clear evidence for $f_{0}(980)$ or $a_{0}(980)^{0}$ (labeled as $S_{0}(980)$) near the $K_{S}^{0}K_{S}^{0}$ threshold is seen, as shown in Fig. 3. In ${\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}\eta}$, a significant enhancement consistent with $N^{*}(1535)$ is found near the $p\eta$ threshold, as shown in Fig. 4. In the future, amplitude analyses of these decays can be expected to improve our understanding of the nature of $S_{0}(980)$ and $N^{*}(1535)$.

In measuring the ratio of branching fractions, many systematic uncertainties cancel, as they affect both the signal and reference modes. The remaining systematic uncertainties are summarized in Table 2 and introduced in detail below.

The systematic uncertainty associated with the $K_{S}^{0}$ reconstruction is considered as follows. A table of $K_{S}^{0}$ efficiency ratios of data to MC in eight bins of the $K_{S}^{0}$ momentum distribution, $R{}_{\varepsilon}^{K_{S}^{0}}$, is determined based on a control sample ${D^{*\pm}\!\rightarrow\!(D^{0}\!\rightarrow\!K_{S}^{0}\pi^{0})\pi^{\pm}}$. The unfolded momentum distribution in data of $K_{S}^{0}$ from signal is obtained using the ${}_{s}\mathcal{P}lot$ technique [42]. From one $R{}_{\varepsilon}^{K_{S}^{0}}$ table, we can determine the average ratios: (1) for $\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK^{0}_{S,{\rm fast}}K^{0}_{S,{\rm slow}}$ where the subscript ‘fast’ (‘slow’) indicates the faster (slower) of two $K_{S}^{0}$’s in the final state, $\overline{R}{}_{\varepsilon}^{K_{S}^{0}}\!=\!\sum_{i}^{8}\sum_{j}^{8}N_{ij}(R{}_{\varepsilon,i}^{K^{0}_{S,{\rm fast}}}R{}_{\varepsilon,j}^{K^{0}_{S,{\rm slow}}})/\sum_{i}^{8}\sum_{j}^{8}N_{ij}$ calculated on the two-dimensional $(p_{\scriptscriptstyle{K^{0}_{S,{\rm fast}}}},\,p_{\scriptscriptstyle{K^{0}_{S,{\rm slow}}}})$ distribution due to the correlations between the momenta of two $K_{S}^{0}$’s. Here $N_{ij}$ and $(R{}_{\varepsilon,i}^{K^{0}_{S,{\rm fast}}}R{}_{\varepsilon,j}^{K^{0}_{S,{\rm slow}}})$ are the yield and the averaged $R{}_{\varepsilon}^{K_{S}^{0}}$, respectively, in the bin of $i^{\rm th}$ raw and $j^{\rm th}$ column of such two-dimensional momenta distribution; (2) for ${\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}}(\eta)$, ${\overline{R}{}_{\varepsilon}^{K_{S}^{0}}\!=\!\sum_{i}^{8}N_{i}R{}_{\varepsilon,i}^{K_{S}^{0}}/\sum_{i}^{8}N_{i}}$ calculated on the one-dimensional $p_{\scriptscriptstyle K_{S}^{0}}$ distribution. Here $N_{i}$ and $R{}_{\varepsilon,i}^{K_{S}^{0}}$ are the yield and the averaged $R{}_{\varepsilon}^{K_{S}^{0}}$, respectively, in the $i^{\rm th}$ bin of such one-dimensional distribution. We build 10000 $R{}_{\varepsilon}^{K_{S}^{0}}$ tables by randomly fluctuating $R{}_{\varepsilon,i}^{K_{S}^{0}}$ in each bin according to its uncertainty and calculate $\overline{R}{}_{\varepsilon}^{K_{S}^{0}}$ for each. We take the mean and root-mean-square (RMS) values from the distribution of $\overline{R}{}_{\varepsilon,\,{\rm sig.}}^{K_{S}^{0}}/\overline{R}{}_{\varepsilon,\,{\rm ref.}}^{K_{S}^{0}}-1$, where the subscripts ‘sig.’ and ‘ref.’ refer to the signal and reference modes, respectively, and add in quadrature as the estimate of the systematic uncertainty.

Since the protons in the signal and reference modes have different kinematic distributions, the systematic effects due to PID do not cancel completely. The data/MC ratio of proton PID efficiency depends on the proton momentum and polar angle: $R{}_{\varepsilon}^{p}(p,\cos\theta)$. Such a $R{}_{\varepsilon}^{p}$ map is determined based on an inclusive sample of $\mathchar 28931\relax\!\rightarrow\!{}p\pi^{-}$. Following steps similar to those used above for $K_{S}^{0}$ efficiency, we obtain the unfolded ($p,\,\cos\theta$) two-dimensional distribution for protons using the ${}_{s}\mathcal{P}lot$ technique [42], and plot the $\overline{R}{}_{\varepsilon,\,{\rm sig.}}^{p}/\overline{R}{}_{\varepsilon,\,{\rm ref.}}^{p}\!-\!1$ values based on 10000 maps of $R_{\varepsilon}^{p}(p,\,\cos\theta)$. The systematic uncertainty due to PID is obtained by adding in quadrature the mean and RMS values of the $\overline{R}{}_{\varepsilon,\,{\rm sig.}}^{p}/\overline{R}{}_{\varepsilon,\,{\rm ref.}}^{p}-1$ distribution.

The uncertainty due to $\eta\!\rightarrow\!\gamma\gamma$ reconstruction is estimated to be 4%, considering 2% per photon according to a study of radiative Bhabha events.

The systematic uncertainties from the $M(\mathchar 28931\relax_{c}^{+})$ fits for ${\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}K_{S}^{0}}$ and ${\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}\eta}$ channels are evaluated to be 1.8% and 2.3%, respectively, after considering two sources below. (a) The uncertainty due to fixing the signal parameters in the fits is estimated by randomly varying them via a multiple-dimensional Gaussian function (including these parameters’ uncertainties and their correlation matrix from the $M(\mathchar 28931\relax_{c}^{+})$ fit of truth-matched signals). We produce 1000 sets of such signal parameters and repeat the $M(\mathchar 28931\relax_{c}^{+})$ fits. We take the ratio of RMS to mean value of the distribution of fitted yield as the relative systematic uncertainty: 0.2% for ${\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}K_{S}^{0}}$, 0.4% for ${\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}\eta}$, and 0.2% for ${\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}}$. (b) To evaluate the potential fit bias, we perform a bias check for the fitted signal yield based on 1000 sets of MC samples, of which the signals are randomly sampled from a large signal MC sample and the backgrounds from the generic $B\overline{B}$ and continuum MC samples. Their sampled yields are equal to the fitted yields in Table 1. We perform $M(\mathchar 28931\relax_{c}^{+})$ fits for these samples. The fitted signal yields are plotted and fitted with a Gaussian function. The shifts of the fitted mean values of the Gaussian functions from the corresponding input values are assigned as systematic uncertainties: 1.9% for ${\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}K_{S}^{0}}$, 2.3% for ${\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}\eta}$, and 0.1% for ${\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}}$. The uncertainties for signal modes and reference mode are added in quadrature, as listed in Table 2.

The systematic effects from the efficiency corrections for the ${\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}K_{S}^{0}}$ and ${\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}\eta}$ channels are evaluated to be 0.8% and 0.4%, respectively, which are obtained by taking the quadratic sum of the following sources: (a) Varying bin size: the 7$\times$7 bins are changed to 6$\times$6 and 8$\times$8 bins for ${\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}K_{S}^{0}}$ and the 5$\times$5 bins are changed to 4$\times$4 and 6$\times$6 bins for ${\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}\eta}$. The changes of efficiency-corrected yields, $0.2\%$ for ${\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}K_{S}^{0}}$ and $0.1\%$ for ${\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}\eta}$, are assigned as the systematic uncertainties. (b) To estimate the uncertainties due to the background Dalitz plot, we shift the $M(\mathchar 28931\relax_{c}^{+})$ sideband region by $\pm 5$ MeV, and repeat the efficiency correction. The resulting changes of efficiency-corrected yields, 0.1% for both channels, are assigned as systematic uncertainty. (c) The signal efficiency effects due to the additional requirements in the signal mode with respect to the reference mode, such as $p(\eta)$, $\chi^{2}_{m}(\eta)$, and $L/\sigma_{L}(K_{S}^{0})$, are neglected, as the signal distributions unfolded from data using the ${}_{s}\mathcal{P}lot$ technique [42] and truth-matched signal distributions from MC are consistent. (d) Systematic effects from the $\chi_{\rm vtx}^{2}$ requirement are considered, since the signal and reference modes have different $\chi_{\rm vtx}^{2}$ distributions. We change the requirement to $\chi_{\rm vtx}^{2}\!<\!21$ and repeat our measurement. The resulting changes to the nominal results, 0.6% and 0.3%, are small as expected and assigned as the corresponding systematic uncertainties. (e) The uncertainty due to the $\pi^{0}$ veto for $\eta$ candidates in ${\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}\eta}$ is estimated by enlarging the veto region from $\pm 12.5$ MeV/$c^{2}$ to be $\pm 15$ MeV/$c^{2}$. The resulting change on the branching fraction is 0.2%, and is assigned as a systematic uncertainty. (f) The uncertainty due to possible data/MC differences in $M(\mathchar 28931\relax_{c}^{+})$ resolution is estimated as follows. Defining $R$ as the ratio of the signal yield in the $M(\mathchar 28931\relax_{c}^{+})$ signal region to that in the fit region, we calculate $r\!=\!R_{\rm data}/R_{\rm MC}$ for the signal and reference modes. The fractional difference in $r$ between signal and reference modes and the uncertainty thereon are summed in quadrature and taken as the systematic uncertainty, which we find to be 0.5% for $\mathcal{B}({\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}K_{S}^{0}})/\mathcal{B}({\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}})$ and 0.1% for $\mathcal{B}({\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}\eta})/\mathcal{B}({\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}})$. (g) The uncertainty due to limited MC statistics for the efficiency value is 0.1%.

The uncertainty due to the non-$K_{S}^{0}$ peaking background is estimated based on the generic MC sample aforementioned. As the rate of this background may depend on intermediate processes, we double its size, and take the resulting ratio with the signal yield, 0.8%, as the associated systematic uncertainty. The uncertainties on $\mathcal{B}(K_{S}^{0}\!\rightarrow\!\pi^{+}\pi^{-})\!=\!{(69.20\pm 0.05)\%}$ ($\delta\mathcal{B}/\mathcal{B}=0.1\%$) and $\mathcal{B}(\eta\!\rightarrow\!\gamma\gamma)\!=\!{(39.41\pm 0.20)\%}$ ($\delta\mathcal{B}/\mathcal{B}=0.5\%$) are also considered. All uncertainties above are added in quadrature to give an overall systematic uncertainty, as listed in Table 2. Additionally, the uncertainty from the world average branching fraction of the reference mode ($5.0\%$) is considered.

In summary, based on the entire dataset with integrated luminosity $980~{}\mbox{\,fb}^{-1}$ collected by the Belle detector at the KEKB energy-asymmetric $e^{+}e^{-}$ collider, we present the first observation of the SCS decay ${\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}K_{S}^{0}}$ with a statistical significance of ${>\!10\sigma}$ and measure the branching fractions of ${\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}K_{S}^{0}}$ and ${\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}\eta}$ relative to ${\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}}$:

where the uncertainties are statistical and systematic, respectively, Using the world average $\mathcal{B}({\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}})=(1.59\pm 0.08)\%$ [17], we obtain the absolute branching fractions

where the first uncertainties are statistical, the second systematic, and the third from the uncertainty on $\mathcal{B}({\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}})$. The first of these branching fractions is measured for the first time and found to be much smaller than the theoretical prediction of $(1.9\pm 0.4)\!\times\!10^{-3}$ [19]. The latter is consistent with the world average, $(4.15\pm 0.90)\!\times\!10^{-3}$ [17], with a threefold improvement in precision.

We reconstruct the Dalitz plots for ${\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}K_{S}^{0}}$ and ${\mathchar 28931\relax_{c}^{+}\!\rightarrow\!{}pK_{S}^{0}\eta}$, with background subtractions and efficiency corrections. We note two clear structures that are consistent with $f_{0}(980)\!\rightarrow\!K_{S}^{0}K_{S}^{0}$ or $a_{0}(980)\!\rightarrow\!K_{S}^{0}K_{S}^{0}$ and $N^{*}(1535)\!\rightarrow\!p\eta$, raising the expectation that the nature of these intermediate resonances will be probed in the future with amplitude analyses on the larger data sets anticipated from BESIII [43] and Belle II [44].

This work, based on data collected using the Belle detector, which was operated until June 2010, was supported by 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, DE220100462, DP150103061, FT130100303; Austrian Federal Ministry of Education, Science and Research (FWF) and FWF Austrian Science Fund No. P 31361-N36; the National Natural Science Foundation of China under Contracts No. 11675166, No. 11705209; No. 11975076; No. 12135005; No. 12175041; No. 12161141008; Key Research Program of Frontier Sciences, Chinese Academy of Sciences (CAS), Grant No. QYZDJ-SSW-SLH011; the Ministry of Education, Youth and Sports of the Czech Republic under Contract No. LTT17020; the Czech Science Foundation Grant No. 22-18469S; Horizon 2020 ERC Advanced Grant No. 884719 and ERC Starting Grant No. 947006 “InterLeptons” (European Union); the Carl Zeiss Foundation, the Deutsche Forschungsgemeinschaft, the Excellence Cluster Universe, and the VolkswagenStiftung; the Department of Atomic Energy (Project Identification No. RTI 4002) and the Department of Science and Technology of India; the Istituto Nazionale di Fisica Nucleare of Italy; National Research Foundation (NRF) of Korea Grant Nos. 2016R1D1A1B02012900, 2018R1A2B3003643, 2018R1A6A1A06024970, RS202200197659, 2019R1I1A3A01058933, 2021R1A6A1A03043957, 2021R1F1A1060423, 2021R1F1A1064008, 2022R1A2C1003993; 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 Ministry of Science and Higher Education of the Russian Federation, Agreement 14.W03.31.0026, and the HSE University Basic Research Program, Moscow; University of Tabuk research grants S-1440-0321, S-0256-1438, and S-0280-1439 (Saudi Arabia); the Slovenian Research Agency Grant Nos. J1-9124 and P1-0135; 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. These acknowledgements are not to be interpreted as an endorsement of any statement made by any of our institutes, funding agencies, governments, or their representatives. 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 6 (SINET6) for valuable network support. We thank Li-Sheng Geng and Ju-Jun Xie for helpful discussions on the $N^{*}(1535)$.