Search for B decays to final states with the eta_c meson

We report a search for $B$ decays to selected final states with the $\eta_c$ meson: $B^{\pm}\to K^{\pm}\eta_c\pi^+\pi^-$, $B^{\pm}\to K^{\pm}\eta_c\omega$, $B^{\pm}\to K^{\pm}\eta_c\eta$ and $B^{\pm}\to K^{\pm}\eta_c\pi^0$. The analysis is based on $772\times 10^6$ $B\bar{B}$ pairs collected at the $\Upsilon(4S)$ resonance with the Belle detector at the KEKB asymmetric-energy $e^+e^-$ collider. We set 90\% confidence level upper limits on the branching fractions of the studied $B$ decay modes, independent of intermediate resonances, in the range $(0.6-5.3)\times 10^{-4}$. We also search for molecular-state candidates in the $D^0\bar{D}^{*0}-\bar{D}^0D^{*0}$, $D^0\bar{D}^0+\bar{D}^0D^0$ and $D^{*0}\bar{D}^{*0}+\bar{D}^{*0}D^{*0}$ combinations, neutral partners of the $Z(3900)^{\pm}$ and $Z(4020)^{\pm}$, and a poorly understood state $X(3915)$ as possible intermediate states in the decay chain, and set 90\% confidence level upper limits on the product of branching fractions to the mentioned intermediate states and decay branching fractions of these states in the range $(0.6-6.9)\times 10^{-5}$.


INTRODUCTION
Many exotic charmonium-like states are observed in the mass region above the DD threshold. Decays of B mesons provide a fruitful opportunity to study these states and to find new ones. For example, the state X(3872) was first observed by Belle in exclusive B + → K + π + π − J/ψ decays [1] that was later confirmed by CDF [2], DØ [3] and BaBar [4]. It was also observed in the LHCb experiment [5,6] in pp collisions and B decays. The X(3872) mass is close to the m D 0 + mD * 0 thresh-old, which engendered a hypothesis that this state may be a D 0D * 0 molecule [7]. The observation of the decay X(3872) → γJ/ψ by BaBar [8] and Belle [9] established the charged parity of X(3872) to be positive. Angular analysis of the X(3872) → J/ψπ + π − decay by LHCb [6] determined all its quantum numbers: J P C = 1 ++ .
If X(3872) is indeed a D 0D * 0 molecule, there can exist other "X-like" molecular states with different quantum numbers. Some may reveal themselves in the decays to final states containing the η c meson. For example, a D 0D * 0 −D 0 D * 0 combination (denoted hereinafter by X 1 (3872)) with quantum numbers J P C = 1 +− would have a mass around 3.872 GeV/c 2 and would decay to η c ρ and η c ω. Combinations of D 0D0 +D 0 D 0 , denoted by X(3730), and D * 0D * 0 +D * 0 D * 0 , denoted by X(4014), with quantum numbers J P C = 0 ++ would decay to η c η and η c π 0 . The mass of the X(3730) state would be around 2m D 0 = 3.730 GeV/c 2 while that of the X(4014) state would be near 2m D * 0 = 4.014 GeV/c 2 .
Recently, a new charged state Z(3900) ± was found in Y (4260) decays by Belle [10] and BESIII [11]. Since this particle is observed in the decay to π ± J/ψ, it should contain at least four quarks. BESIII [12] reported subsequently an observation of another decay channel of the seemingly same state Z(3885) ± → (DD * ) ± . The Z(3900) ± was confirmed in the decay to π ± J/ψ by an analysis of CLEO-c data [13] that also reported evidence for its neutral isotopic partner Z(3900) 0 . Another exotic charged state Z(4020) ± was observed by BESIII in decays to π ± h c [14] and (D * D * ) ± [15]. There are some indications from these analyses that the spin and parity of the charged states might be J P = 1 + .
The near-threshold enhancement in the ωJ/ψ invariant mass distribution named Y (3940) was first observed by Belle in exclusive B → KωJ/ψ decays [16]. Later, in the same decay mode, BaBar discovered X(3915) [17], which was confirmed by Belle in two-photon production [18] and by other BaBar measurements [19,20]. The parameters of Y (3940) are consistent with those of X(3915), so they are considered to be the same particle. The quantum numbers of X(3915) are claimed to be J P C = 0 ++ , but its nature is still undetermined, and there are several interpretations describing this state [21][22][23][24][25].

EVENT SELECTION
The analysis is based on a data sample that contains 772 × 10 6 BB pairs, collected with the Belle detector at the KEKB asymmetric-energy e + e − collider [26,27] operating at the Υ(4S) resonance.
The Belle detector [28] (also see detector section in [29]) is a large-solid-angle magnetic spectrometer that consists of a silicon vertex detector, a 50-layer central drift chamber (CDC) for charged particle tracking and specific ionization measurement (dE/dx), an array of aerogel threshold Cherenkov counters (ACC), time-offlight scintillation counters (TOF), and an array of 8736 CsI(Tl) crystals for electromagnetic calorimetry located inside a superconducting solenoid coil that provides a 1.5 T magnetic field. An iron flux return yoke located outside the coil is instrumented to detect K 0 L mesons and identify muons. We use a GEANT3-based Monte Carlo (MC) simulation to model the response of the detector and determine its acceptance [30].
Charged tracks are selected with requirements based on the goodness of fit of the tracks and their impact parameters relative to the interaction point (IP). We require that the polar angle of each track be in the angular range [18 • , 152 • ] and that the track momentum perpendicular to the positron beamline be greater than 100 MeV/c.
Pions and kaons are distinguished by combining the responses of the ACC and the TOF with dE/dx measurements in the CDC to form a likelihood ratio L i/j = L i /(L i + L j ). Here, L i is the likelihood that the particle is of type i. K 0 S candidates are reconstructed via the π + π − final state. The π + π − invariant mass must lie in the range 0.486 GeV/c 2 < M (π + π − ) < 0.510 GeV/c 2 . The flight length of the K 0 S is required to lie within the interval [0. 1,20] cm. The angle ϕ between the pion pair momentum and the line joining the π + π − the vertex to the IP must satisfy cos ϕ > 0.95.
B candidates are identified by their center-of-mass (c.m.) energy difference ∆E = ( i E i ) − E b , and the beam-constrained mass M bc = is the beam energy in the Υ(4S) c.m. frame, and p i and E i are the c.m. three-momenta and energies, respectively, of the B candidate decay products. We define the signal region by |∆E| < 0.02 GeV and 5.273 GeV/c 2 < M bc < 5.285 GeV/c 2 ; the |∆E| range extends to 0.04 GeV for some decay channels.
Although the continuum background (e + e − → qq, where q = u, d, s, c) is not dominant, we suppress it using topological criteria. Since the produced B mesons are nearly at rest in the c.m. frame, the signal tends to be isotropic, while the continuum qq background tends to have a two-jet structure. We use the angle Θ thrust between the thrust axis of the B candidate and that of the rest of the event to discriminate between these two cases. The distribution of | cos Θ thrust | is strongly peaked near 1 for qq events but nearly uniform for Υ(4S) → BB events; we require | cos Θ thrust | < 0.8. In addition, we require that the c.m. polar angle θ B of the reconstructed B meson satisfies | cos θ B | < 0.8, where this angle is measured relative to the z axis that is collinear with the positron beam.
The mean number of multiple B candidates per event varies from 1.2 to 2.0, depending on the decay channel. If an event has multiple B candidates, we select the candidate with the minimum value of the following expressions in the order described below: 3. |m η − M (γγ)| (for η → γγ in the (η) mode) or |m π 0 − M (γγ)| (for the (ω) mode, for η → π + π − π 0 in the (η) mode and for the (π 0 ) mode), 4. maximum difference between z coordinates at the closest-distance point to the z axis among each pair of charged particles in the signal-B final state.
According to MC, this procedure selects the correct B candidate with 95% probability.

RECONSTRUCTION OF FINAL STATES
We search for D ( * )0 D ( * )0 molecular-state candidates Z(3900) 0 , Z(4020) 0 , and X(3915) in the following B meson decays: To determine the branching fractions, we perform a binned maximum-likelihood fit of the ∆E distribution that is modelled by a peaking signal and featureless background. For the (π + π − ), (ω) and (π 0 ) decay modes, the signal function is the sum of two Gaussians (G) and the background function is a linear polynomial; the fitting function is (1) For a given mode, we generate signal MC and fix the mean values x 1 and x 2 , the standard deviations σ 1 and σ 2 , and the fraction of the first Gaussian α; we also obtain the detection efficiency for the mode. Here and in the following the detector resolution is taken into account. To account for the difference in resolution between MC and data, we replace the resolution of the dominant second Gaussian in the fit with σ ′ 2 = σ 2 2 + δ 2 , where the resolution degradation δ = (7.1 ± 2.3) MeV/c 2 is taken from the analysis of the decay [31]. The ∆E distribution for the (ω) mode is shown in Fig. 1. In the (π + π − ) and (π 0 ) modes, we observe some significant signal and so perform a two-dimensional fit of the K S Kπ invariant mass (x) and ∆E (y) distributions: for the (π + π − ) and (π 0 ) modes, respectively. The logarithmic Gaussian function is defined as and P is the asymmetry parameter. The function u(x) characterizes the η c resonance and is described by the convolution of a Breit-Wigner (b) function and a Gaussian detector resolution function with σ res = 6.2 MeV/c 2 obtained from Ref. [31]. The parameter N non-res represents the number of events that do not contain an intermediate η c meson, but have the same final state. According to the fit, most of the events in the ∆E peak are of that origin. The M (K S Kπ) and ∆E distributions for (π + π − ) and (π 0 ) modes are shown in Figs. 2 and 3, respectively.
To improve the X invariant mass resolution in the (η) mode, we modify the energy of the η candidate decaying into photons: [32] is the mass and p η is the reconstructed momentum. Since the η candidate is reconstructed in two decay modes, we perform a combined fit of the ∆E distribution corresponding to η → γγ and η → π + π − π 0 , using the following function: where i refers to either η → γγ or η → π + π − π 0 decay. In particular, B 2γ = B(η → γγ) and B 3π = B(η → π + π − π 0 ) × B(π 0 → γγ). Here, the parameter N eff is the effective number of signal events. To obtain the total yield for each decay channel, it should be multiplied by the corresponding efficiency ε 2γ/3π and η decay branching fraction B i . The combined fit projections of the ∆E distributions for the (η) mode are shown in Fig. 4. The fit results are summarized in Table I.

X1(3872), X(3730) and X(4014)
For the (π + π − ) mode, we perform a maximumlikelihood fit of the η c π + π − invariant mass distribution that is, again, modeled by a peaking signal and a smooth background; the fit function is (6) For each intermediate resonance, we generate signal MC that incorporates the corresponding X quantum numbers. From the signal-MC fit, we fix the mean value M , the standard deviations σ 1 and σ 2 , and the Gaussian fraction α; in addition, we obtain the detection efficiency for the mode. From the fit, we obtain the signal yield N s , which is shown in Table II. The η c π + π − invariant mass distribution is shown in Fig. 5 (left).
We validate our procedure by applying it to the decay B ± → K ± ψ(2S), ψ(2S) → J/ψπ + π − . This decay is similar to the (π + π − ) decay except that we reconstruct the ψ(2S) meson in place of the X 1 (3872) and the J/ψ in place of the η c . The J/ψ meson, like the η c , is reconstructed via the K 0 S K ± π ∓ final state. The selection criteria for this decay are the same as for the (π + π − ) decay except for the invariant mass of the K 0 S K ± π ∓ combination: 3.077 GeV/c 2 < M (K 0 S K ± π ∓ ) < 3.117 GeV/c 2 . The mean number of B candidates per event is 1.7. In case of multiple B candidates, the one with the minimum differences for |m K 0 S − M (π + π − )| and |m J/ψ − M (K 0 S K ± π ∓ )| and the best vertex coordinate is chosen. We fit the J/ψπ + π − invariant mass distribution and obtain the number of signal events N s = 20.2 ± 6.5, which corresponds to a significance of 3.5 standard deviations (σ). The significance is estimated using the value of −2 ln L0 Lmax , where L max (L 0 ) denotes the likelihood value when the yield is allowed to vary (is set to zero). The expected number of events estimated using the world averages of the known branching fractions [32] is 22 ± 4, which is consistent with N s .
In the analysis of the (ω) mode, we use the sum of two Gaussians to describe the signal and a threshold squareroot function to describe the background: The η c ω invariant mass distribution is shown in Fig. 5 (right). In the X(3730) mass region of the (η) mode, the fitting function is where N eff is the effective number of signal events and i refers to either η → γγ or η → π + π − π 0 decay. The η c η invariant mass distributions in the X(3730) mass region are shown in Fig. 6 (top).
In the X(4014) mass region of the (η) mode, the fitting function is The η c η invariant mass distributions in the X(4014) mass region are shown in Fig. 6 (bottom).
For the the (π 0 ) mode, we use The η c π 0 invariant mass distributions are shown in Fig. 7.
The fit results are summarized in Table II.

X(3915)
For the η c η invariant mass distribution, we perform a combined fit of two decay modes of the η meson: (12) where N eff is the effective number of signal events and i refers to either η → γγ or η → π + π − π 0 decay. The detector resolution σ i is obtained from signal MC, taking into account the resolution degradation and is equal to 13.6 MeV/c 2 for η → γγ and 12.3 MeV/c 2 for η → π + π − π 0 . The corresponding η c η invariant mass distribution is shown in Fig. 9.
We fit the η c π 0 invariant mass distribution with the function in Eq. (11). The detector resolution is obtained from signal MC, taking into account the resolution degradation, and is equal to 15.7 MeV/c 2 . The corresponding η c π 0 invariant mass distribution is shown in Fig. 10.
The fit results are summarized in Table III.

SYSTEMATIC UNCERTAINTIES
Systematic uncertainties are categorized as follows: 1. Additive systematic uncertainties affect the number of signal events and are estimated by the variation of the fit conditions. They are displayed as the numbers of events in Tables IV and V and arise from the sources listed below. (i) To obtain the error related to the resolution degradation, we vary the corrected variance of the fitting function within its statistical uncertainty. (ii) We assume that the combinatorial background can be parameterized with a first-order polynomial. To obtain the background shape uncertainty, we describe the background by a second-order polynomial and compare the results. For the X → η c ω channel, we change the square root to the fourth root. (iii) To estimate the systematic error associated with the selection criteria, we relax the criteria on ∆E, M bc , and the invariant masses of η c , ω and η by 50%. (iv) We vary the bin size from 2.5 to 7.5 MeV/c 2 and determine the corresponding systematic error.
2. Multiplicative systematic uncertainties affect the product branching fractions. They are displayed in percent in Table VI and arise from the sources listed below. (i) The number of BB pairs is calculated from the difference between the number of hadronic events on resonance and the scaled number of those off resonance. The systematic error is dominated by the uncertainty in the scale factor and is equal to ∼1.4%. (ii) The uncertainties on the ω, π 0 , η, η c and K 0 S decay branching fractions are taken from Ref. [32]. (iii) The statistical error of the efficiency determined by the signal MC is also taken as a systematic uncertainty. For Z 0 decays, we take into account the efficiency variation in the η c ρ and η c f 0 decay modes. For the decays B ± → K ± η c π + π − and B ± → K ± η c π 0 , we assume alternate decay modes containing intermediate K * 0 , K * (1410) ± and ρ mesons, and take into account the difference of the MC detection efficiency as the corresponding systematic uncertainty. (iv) An analysis of the charged track reconstruction uncertainty as a function of particle momentum gives an estimate of 0.34% per charged track.
(v) To determine the errors due to K and π meson identification, data from the analysis of the process D * + → D 0 π + followed by the decay D 0 → K − π + are used. The uncertainty in K ± identification is 0.8% per K meson and the corresponding value for π ± identification is 0.5% per π meson. (vi) Estimation of the contribution of the η and π 0 reconstruction uncertainty is carried out using the comparison of the number of reconstructed η → 3π 0 and η → γγ events. Such an estimate gives 2% per η and π 0 meson. (vii) The contribution of the K 0 S reconstruction uncertainty is estimated to be 4.4% [33]. (viii) We also take into account the deviation of MC from the data by applying a correction to the efficiency: ε Data /ε MC is 0.9996 for each kaon and 0.9756 for each pion.
In the case of the X → η c η decay, the η meson is reconstructed in two decay modes: η → γγ and η → π + π − π 0 . To estimate the systematic uncertainty on the reconstruction of X, we compare the systematic errors given by each η decay mode and take the maximum as the uncertainty for the reconstruction.
For the Z 0 decays to η c π + π − , the additive systematic uncertainty is assumed to be the same as in X 1 (3872) → η c π + π − . In addition, we vary the resonance width in the [15,65] MeV/c 2 interval for the Z(3900) 0 and [2.5, 27.5] MeV/c 2 interval for the Z(4020) 0 . The intervals are chosen according to the variation of the previously measured Z ± width values. The total additive systematic uncertainty is 28.5 events for the Z(3900) 0 and 25.9 events for the Z(4020) 0 .
For the X(3915) decays, the systematic uncertainty is taken from the similar X(4014) decays (see Tables V  and VI).

RESULTS
Our data sample does not permit us to measure the branching products of production and decay of the states listed above nor the B decay branching fractions so we set upper limits instead, taking into account the statistical and systematic uncertainties. The expression for the required branching value is z = N s /(N B εB i ), where N s is the yield, N B is the number of BB pairs, ε is the detection efficiency, and B i is the product of the branching fractions of the intermediate resonances. This expression can be written as z = xy, where x = N s and y = 1 N BB εBi . The distribution of the numerator N s is assumed to be Gaussian with mean µ and standard deviation σ x = σ 2 stat + σ 2 add.syst . The distribution of the 6.8 6.8 6.8 6.8 B(K 0 S → π + π − ) 0. inverse of the denominator is also assumed to be Gaussian with mean ν = 1/(N B εB i ) and standard deviation σ y = σ mult.syst ν. Thus, the distribution of z can be written in the following way: The 90% confidence level (C.L.) upper limit U on z is defined by Upper limits on the branching fractions and products for all the studied decay modes are shown in Tables VII and VIII.  Figure 11 shows the dependence on the mass bin of the upper limit for the branching product of the Z(3900) 0 and Z(4020) 0 production and decay: no significant signal is seen in any of the invariant mass bins. If we assume that the Z 0 mass is close to that of its charged partner Z ± , we obtain the upper limits on the product branching fractions shown in Table VIII. Similarly, for the study of X(3872)-like particles, we performed a mass scan inside the fitting region, i.e., a sequence of fits similar to the ones described above but with an X mass floating in 20 MeV/c 2 -wide mass bins. No significant signal is found in any of the studied bins.
There are no theoretical predictions for the decay branching fractions of D 0D * 0 molecular states similar to X(3872). In this paper, we obtain an upper limit on the product branching fraction B(B ± → K ± X 1 (3872))× B(X 1 (3872) → η c π + π − ), which is of the same order as the product branching fraction B(B ± → K ± X(3872)) × B(X(3872) → J/ψπ + π − ) measured in Ref. [35]. A similar situation is observed with the upper limit on the B(B ± → K ± X(3915)) × B(X(3915) → η c ω) and the value of B(B ± → K ± X(3915)) × B(X(3915) → J/ψω) obtained in Ref. [19]. A more copious data set expected from the upcoming Belle II experiment [36] can provide an opportunity to determine the ratios of the decay branching fractions of the exotic states described above.