Measurement of branching fractions of $\Lambda_c^+\to{}pK_S^0K_S^0$ and $\Lambda_c^+\to{}pK_S^0\eta$ at Belle

We present a study of a singly Cabibbo-suppressed decay $\Lambda_c^+\to{}pK_S^0K_S^0$ and a Cabibbo-favored decay $\Lambda_c^+\to{}pK_S^0\eta$ based on 980 $\rm fb^{-1}$ of data collected by the Belle detector, operating at the KEKB energy-asymmetric $e^+e^-$ collider. We measure their branching fractions relative to $\Lambda_c^+\to{}pK_S^0$: $\mathcal{B}(\Lambda_c^+\to{}pK_S^0K_S^0)/\mathcal{B}(\Lambda_c^+\to{}pK_S^0)={(1.48 \pm 0.08 \pm 0.04)\times 10^{-2}}$ and $\mathcal{B}(\Lambda_c^+\to{}pK_S^0\eta)/\mathcal{B}(\Lambda_c^+\to{}pK_S^0)={(2.73\pm 0.06\pm 0.13)\times 10^{-1}}$. Combining with the world average $\mathcal{B}(\Lambda_c^+\to{}pK_S^0)$, we have the absolute branching fractions: $\mathcal{B}(\Lambda_c^+\to{}pK_S^0K_S^0) = {(2.35\pm 0.12\pm 0.07 \pm 0.12 )\times 10^{-4}}$ and $\mathcal{B}(\Lambda_c^+\to{}pK_S^0\eta) = {(4.35\pm 0.10\pm 0.20 \pm 0.22 )\times 10^{-3}}$. The first and second uncertainties are statistical and systematic, respectively, while the third ones arise from the uncertainty on $\mathcal{B}(\Lambda_c^+\to{}pK_S^0)$. The mode $\Lambda_c^+\to{}pK_S^0K_S^0$ is observed for the first time and has a statistical significance of $>\!10\sigma$. The branching fraction of $\Lambda_c^+\to{}pK_S^0\eta$ has been measured with a threefold improvement in precision over previous results and is found to be consistent with the world average.


I. INTRODUCTION
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 nonfactorizable contributions. The latter may play a nontrivial 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 Cabibbosuppressed (SCS) charm decays are essential probes of CP 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, Λ + 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 Λ + c → pK 0 S η [18], the world average branching fraction, B(Λ + c → pK 0 S η) = (4.15 ± 0.90) × 10 −3 [17], still has a large uncertainty (22%). The SCS mode Λ + c → pK 0 S K 0 S , for which the predicted branching fraction is B(Λ + c → pK 0 S K 0 S ) = (1.9 ± 0.4) × 10 −3 based on SU(3) F symmetry [19], has not previously been observed.
In this paper, we present a precise measurement of B(Λ + c → pK 0 S K 0 S ) and B(Λ + c → pK 0 S η) 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 ηN and KΛ, 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 Λ + c → pK 0 S η decay, in which the final-state pη 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η, in S-wave. Other intermediate resonances of interest are the light scalars a 0 (980) and f 0 (980), which both couple to KK in Λ + c → pK 0 S K 0 S . They contribute to the SCS Λ + c decays Λ + c → pKK and Λ + c → pππ, as predicted in Ref. [24], and likely contribute to Λ + c → pK 0 S K 0 S , 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 KK bound states [31,32]. Therefore, we reconstruct the Dalitz plots of Λ + c → pK 0 S K 0 S and Λ + c → pK 0 S η decays to check such interesting intermediate resonances.

II. DETECTOR AND DATA SET
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 fb −1 collected at or near the Υ(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-offlight 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 BB events and continuum processes e + e − → qq (q = u, d, s, c) corresponding to an integrated luminosity three times that of the data. Samples of MC events of Λ + c signal decay modes are produced in the e + e − → cc process, decayed uniformly in three-body phase space, and used to study the efficiency.

III. EVENT SELECTION
We reconstruct the two signal modes Λ + c → pK 0 S K 0 S and Λ + c → pK 0 S η and their reference mode Λ + c → pK 0 S . The event selections are optimized based on a figure of merit (FOM), defined as FOM = ε S / √ N B for Λ + c → pK 0 S K 0 S due to its branching fraction having not yet been measured, and FOM = N S / √ N S + N B for Λ + c → pK 0 S η assuming its current world average branching fraction [17]. Here ε 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 (Λ + c ) signal regions, where M (Λ + c ) is the invariant mass of reconstructed Λ + c candidates. These signal regions are defined to be within 10, 22, and 18 MeV/c 2 of the nominal Λ + c mass [17] for the Λ + c → pK 0 S K 0 S , Λ + c → pK 0 S η, and Λ + c → pK 0 S channels, respectively; each signal band includes ≈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 (Λ + c ) sideband region ( [17]. The particle identification (PID) likelihood for a given particle hypothesis, L i (i = π, 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 R(p|K) = L p /(L p + L K ) > 0.6 and R(p|π) = L p /(L p + L π ) > 0.6, are identified as protons. These PID requirements have signal efficiencies of 94% for Λ + c → pK 0 S K 0 S and 97% for Λ + c → pK 0 S η. 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 0 S '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 0 S momentum in the laboratory frame; the separation in z between the two π ± 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 0 S flight length in the x-y plane; the angle between the K 0 S momentum and the vector joining the IP to the K 0 S decay vertex; in the K 0 S rest frame, the angle between the π + momentum and the laboratory-frame boost direction; and, for each π ± 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 0 S → π + π − candidate is required to lie within 10 MeV/c 2 of the nominal K 0 S mass [17]; this includes 99.9% of the K 0 S signal. The two pion tracks from each K 0 S candidate are refitted to originate from a common vertex and constrained to have invariant mass equal to the nominal K 0 S mass [17]. The corresponding fit quality χ 2 mv (K 0 S ) is required to be smaller than 100. The selected K 0 S 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×3 array of crystals centered on the crystal with the highest energy, to the energy deposited in the corresponding 5×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 • < θ < 129 • ), and greater than 100 MeV in the endcap region (12 • < θ < 31 • or 132 • < θ < 157 • ).
Candidate η → γγ decays are reconstructed from photon pairs having an invariant mass satisfying 500 MeV/c 2 < M (γγ) < 580 MeV/c 2 (3σ in M η (γγ) resolution). The invariant mass of each η candidate is constrained to the nominal η mass [17] at the Λ + c decay vertex (described below). The fit quality of this mass constraint is required to satisfy χ 2 m (η) < 8, and the resulting η momentum in the laboratory frame is required to be greater than 0.4 GeV/c. To further suppress the background, η candidates are vetoed if either of daughters can be paired with another photon such that the γγ pair has an invariant mass within 2.5σ of the nominal π 0 mass (σ = 5 MeV/c 2 ). This π 0 -veto results in a signal loss of 28% and removes 72% of background.
The Λ + c candidates are assembled by forming combinations of the final-state particles for each mode. The p and K 0 S are required to originate from a common vertex (denoted the Λ + c decay vertex and the K 0 S production vertex) with a fit quality χ 2 vtx < 24. To reduce combinatorial background, the scaled momentum of the Λ + c candidate, defined as x p = p * c/ s/4 − M 2 (Λ + c ) · 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 Λ + c candidates in the e + e − center-of-mass frame.
For the SCS decay Λ + c → pK 0 S K 0 S , a non-K 0 S peaking background from the CF decay Λ + c → pK 0 S π + π − exists, even though it is suppressed by the vertex fit and K 0 S selection. The K 0 S decay length L is determined by the projection of the vector joining the K 0 S production and decay vertices onto the K 0 S momentum direction, and its corresponding uncertainty σ L is calculated by propagating uncertainties in the vertices and the K 0 S momentum, including their correlations. To suppress the non-K 0 S peaking CF background, we require the significance of the K 0 S decay length L/σ L (K 0 S ) > 10 for the slower of the two K 0 S 's in Λ + c → pK 0 S K 0 S . This requirement reduces the signal efficiency by 3%, and rejects 80% of non-K 0 S peaking background. The remaining non-K 0 S peaking background is ignored in the M (Λ + c ) fits because it has a tiny ratio 0.4% to signal based on the MC studies with the branching fraction (1.6 ± 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 . 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.

IV. YIELD EXTRACTION
The signal yield is extracted by an unbinned extended maximum likelihood fit to the M (Λ + c ) distribution. The signal probability density function (PDF) is a sum of three symmetric Gaussian functions for the Λ + c → pK 0 S K 0 S mode, a sum of one symmetric Gaussian and two asymmetric Gaussians for the Λ + c → pK 0 S η mode, and a sum of one symmetric Gaussian and three asymmetric Gaussians for the Λ + c → pK 0 S 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 (δ µ ) from the value found in MC, and the widths are those found in MC, multiplied by a common scaling factor (k σ ). The background PDF is a first-order polynomial function for Λ + c → pK 0 S K 0 S and a second-order polynomial function for Λ + c → pK 0 S η and Λ + c → pK 0 S . 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 data −N fit )/σ data where σ data is the error on N data . The pull distributions demonstrate that the data are statistically consistent with the fitted shapes. The signal and background yields are listed in Table I. For the Λ + c → pK 0 S K 0 S mode, we obtain the difference in the log likelihoods obtained from fits performed with and without a signal PDF, ∆ ln L = 524; as the number of degrees of freedom without a signal component is three less than that with a signal component (parameters N sig , δ µ and k σ are dropped), and this value of ∆ ln L corresponds to a statistical significance greater than 10σ. This measurement constitutes the first observation of this SCS Λ + c decay.
TABLE I. The fitted yields of signal and background in the overall fit region (FR) and the signal region (SR) for the Λ + c → pK 0 S K 0 S , Λ + c → pK 0 S η, and Λ + c → pK 0 S modes. For the definition of these regions, see the text. The yields in signal region, N SR bkg of Λ + c → pK 0 S (K 0 S , η) and N SR sig of Λ + c → pK 0 S , are used to measure the branching fractions. For the three-body decay modes, the Dalitz plots for candidates in the M (Λ + c ) signal region and sideband region are shown in Figs. 2(a, b) for Λ + c → pK 0 S K 0 S and Figs. 2(d, e) for Λ + c → pK 0 S η. For Λ + c → pK 0 S K 0 S , Bose symmetry requires invariance under the exchange of the two K 0 S 's, hence the Dalitz plot for two pK 0 S masses is symmetric. We plot M 2 (pK 0 S ) max versus M 2 (pK 0 S ) 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 Λ + c → pK 0 S K 0 S and Λ + c → pK 0 S η, the efficiencies are calculated in bins across the phase space, based on truthmatched signal yield in the M (Λ + c ) signal region. The results are shown in Fig. 2(c) for Λ + c → pK 0 S K 0 S and Fig. 2(f) for Λ + c → pK 0 S η. 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×7 bins for Λ + c → pK 0 S K 0 S and 5×5 bins for Λ + c → pK 0 S η, as shown in Figs. 2(c, f) respectively. The efficiency-corrected yields are where N tot i is the raw yield in the i th bin of the Dalitz plot in M (Λ + c ) signal region, N SR bkg is the fitted back-ground yield as listed in Table I, f bkg i is the fraction of background in the i th -bin, with i f i = 1. These fractions are obtained from the Dalitz plot distribution of events in the M (Λ + c ) sideband region, shown in Fig. 2 Fig. 2(e) for Λ + c → pK 0 S η. Using the generic MC sample, we find that the Dalitz plot in the chosen M (Λ + c ) sideband region is consistent with the generic background in the M (Λ + c ) signal region. The uncertainties on each variable in Eq. (1) have been considered and are propagated into the efficiency-corrected yields, N corr . We obtain The relative branching fractions of signal modes to reference mode are determined by Eqs. (4,5).
We examine the Dalitz plots for Λ + c → pK 0 S K 0 S and Λ + c → pK 0 S η, after background subtraction and efficiency correction, for intermediate resonances.
In Λ + c → pK 0 S K 0 S , clear evidence for f 0 (980) or a 0 (980) 0 (labeled as S 0 (980)) near the K 0 S K 0 S threshold is seen, as shown in Fig. 3. In Λ + c → pK 0 S η, a significant enhancement consistent with N * (1535) is found near the pη 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).

VI. SYSTEMATIC UNCERTAINTY
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 II and introduced in detail below.
The systematic uncertainty associated with the K 0 S reconstruction is considered as follows. A table of K 0 S efficiency ratios of data to MC in eight bins of the K 0 S momentum distribution, R K 0 S ε , is determined based on a control sample D * ± → (D 0 → K 0 S π 0 )π ± . The unfolded momentum distribution in data of K 0 S from signal is obtained using the s Plot technique [42]. From one R K 0 S ε table, we can determine the average ratios: (1) for Λ + c → pK 0 S,fast K 0 S,slow where the subscript 'fast' ('slow') indicates the faster (slower) of two K 0 S 's in the final state, R S,fast , p K 0 S,slow ) distribution due to the correlations between the momenta of two K 0 S 's. Here N ij and (R  c → pK 0 S K 0 S , the Dalitz plot after background subtraction and efficiency correction bin-by-bin and its projections superimposing with signal MC produced by phase space mode (blue histograms). This symmetric Dalitz plot and its projections show two entries per candidate, one for each possible pK 0 S combination. A dominant structure near the K 0 S K 0 S threshold, which we identify with f0(980) or a0(980) 0 , is clearly seen. 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 p ε (p, cos θ). Such a R p ε map is determined based on an inclusive sample of Λ → pπ − . Following steps similar to those used above for K 0 S efficiency, we obtain the unfolded (p, cos θ) two-dimensional distribution for protons using the s Plot technique [42], and plot the R p ε, sig. /R p ε, ref. −1 values based on 10000 maps of R p ε (p, cos θ). The systematic uncertainty due to PID is obtained by adding in quadrature the mean and RMS values of the R p ε, sig. /R p ε, ref. − 1 distribution. The uncertainty due to η → γγ reconstruction is estimated to be 4%, considering 2% per photon according to a study of radiative Bhabha events.
The systematic uncertainties from the M (Λ + c ) fits for Λ + c → pK 0 S K 0 S and Λ + c → pK 0 S η 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 (Λ + c ) fit of truth-matched signals). We produce 1000 sets of such signal parameters and repeat the M (Λ + 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 Λ + c → pK 0 S K 0 S , 0.4% for Λ + c → pK 0 S η, and 0.2% for Λ + c → pK 0 S . (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 BB and continuum MC samples. Their sampled yields are equal to the fitted yields in Table I. We perform M (Λ + 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 sys-tematic uncertainties: 1.9% for Λ + c → pK 0 S K 0 S , 2.3% for Λ + c → pK 0 S η, and 0.1% for Λ + c → pK 0 S . The uncertainties for signal modes and reference mode are added in quadrature, as listed in Table II. The systematic effects from the efficiency corrections for the Λ + c → pK 0 S K 0 S and Λ + c → pK 0 S η 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×7 bins are changed to 6×6 and 8×8 bins for Λ + c → pK 0 S K 0 S and the 5×5 bins are changed to 4×4 and 6×6 bins for Λ + c → pK 0 S η. The changes of efficiency-corrected yields, 0.2% for Λ + c → pK 0 S K 0 S and 0.1% for Λ + c → pK 0 S η, are assigned as the systematic uncertainties. (b) To estimate the uncertainties due to the background Dalitz plot, we shift the M (Λ + c ) sideband region by ±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(η), χ 2 m (η), and L/σ L (K 0 S ), are neglected, as the signal distributions unfolded from data using the s Plot technique [42] and truth-matched signal distributions from MC are consistent. (d) Systematic effects from the χ 2 vtx requirement are considered, since the signal and reference modes have different χ 2 vtx distributions. We change the requirement to χ 2 vtx < 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 π 0 veto for η candidates in Λ + c → pK 0 S η is estimated by enlarging the veto region from ±12.5 MeV/c 2 to be ±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 (Λ + c ) resolution is estimated as follows. Defining R as the ratio of the signal yield in the M (Λ + c ) signal region to that in the fit region, we calculate r = R data /R 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 B(Λ + c → pK 0 S K 0 S )/B(Λ + c → pK 0 S ) and 0.1% for B(Λ + c → pK 0 S η)/B(Λ + c → pK 0 S ). (g) The uncertainty due to limited MC statistics for the efficiency value is 0.1%.
The uncertainty due to the non-K 0 S 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 B(K 0 S → π + π − ) = (69.20 ± 0.05)% (δB/B = 0.1%) and B(η → γγ) = (39.41 ± 0.20)% (δB/B = 0.5%) are also considered. All uncertainties above are added in quadrature to give an overall systematic uncertainty, as listed in Table II. Additionally, the uncertainty from the world average branching fraction of the reference mode (5.0%) is considered.
We reconstruct the Dalitz plots for Λ + c → pK 0 S K 0 S and Λ + c → pK 0 S η, with background subtractions and efficiency corrections. We note two clear structures that are consistent with f 0 (980) → K 0 S K 0 S or a 0 (980) → K 0 S K 0 S and N * (1535) → pη, 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].