The Belle Collaboration

The $Y(4260)$ state, also known as $\psi(4260)$ pdg , was first seen by the BABAR collaboration in 2005 BaBar in the initial-state radiation (ISR) process $e^{+}e^{-}\to\gamma_{\rm ISR}Y(4260),~Y(4260)\to J/\psi\pi^{+}\pi^{-}$ and confirmed by the Belle and CLEO collaborations using the same process belle ; cleo . The world average mass and decay width of the $Y(4260)$ are $(4230\pm 8)$ MeV/$c^{2}$ and $(55\pm 19)$ MeV pdg , respectively. Due to its observation in ISR production, the $J^{PC}$ of the $Y(4260)$ is expected to be $1^{--}$. The decay of $Y(4260)$ to $J/\psi\pi^{+}\pi^{-}$ indicates the presence of a $c\bar{c}$ pair among its quark constituents. However, its mass and properties are not consistent with those expected for any of the $c\bar{c}$ states in the charmonium spectrum, which makes it problematic to assign the $Y(4260)$ to one of the conventional $c\bar{c}$ states with $J^{PC}=1^{--}$.

Attempts have been made to identify $Y(4260)$ as a candidate for a mixed state, which is an admixture of charmonium and tetraquark states admix , a hybrid charmonium state, which is a bound state of charmonium with a gluon t_5 , a tetraquark tetra , a mesonic molecule t_3 ; t_2 ; t_4 , or a charmonium baryonium t_6 . The $Z_{c}(3900)^{\pm}$ state, which as it is charged makes it a natural tetraquark candidate, has been observed by the BESIII and Belle collaborations in the $J/\psi\pi^{\pm}$ invariant mass spectrum of the $e^{+}e^{-}\to Y(4260)\to J/\psi\pi^{+}\pi^{-}$ process besiii ; belle_Zc , provides further evidence of the unconventional nature of the $Y(4260)$.

A mixed-state model, based upon a QCD sum-rule approach theo , suggests the possible interval on the product of the branching fractions of $B^{+}\to Y(4260)K^{+},~Y(4260)\to J/\psi\pi^{+}\pi^{-}$ to be in the range $3.0\times 10^{-8}-1.8\times 10^{-6}$. The BABAR collaboration has measured a signal for the charged $B$ decay with a statistical significance of 3.1 standard deviations ($\sigma$) based on a data sample of 211 $\rm fb^{-1}$ which contains $(232\pm 3)\times 10^{6}B\bar{B}$ pairs babar . They set the upper limit at the 95% confidence interval to be $\mathcal{B}(B^{+}\to Y(4260)K^{+})\times\mathcal{B}(Y(4260)\to J/\psi\pi^{+}\pi^{-})<2.9\times 10^{-5}$. Further improvement is required on the precision of both the theoretical estimate and experimental measurement to elucidate the structure of $Y(4260)$.

Recently, two resonance structures have been observed by the BESIII collaboration in a fit to the cross section of the $e^{+}e^{-}\to J/\psi\pi^{+}\pi^{-}$ process besIII . The resonance structures are interpreted as $Y(4260)$ and $Y(4360)$ with measured masses $(4222.0\pm 3.1\pm 1.4)$ MeV/$c^{2}$ and $(4320.0\pm 10.4\pm 7.0)$ MeV/$c^{2}$, respectively. The measured $Y(4260)$ mass is not significantly lower than world average pdg , from which it deviates merely about 1$\sigma$, and the $Y(4360)$ has not yet been confirmed. We assume the presence of $Y(4260)$ only in the $J/\psi\pi^{+}\pi^{-}$ invariant mass region of interest as in the previous measurements belle ; BaBar ; cleo instead of adopting the search for the improved mass region.

In this paper, we report a search for $B\to Y(4260)K,~Y(4260)\to J/\psi\pi^{+}\pi^{-}$ charge decays based on a data sample corresponding to an integrated luminosity of 711 $\rm fb^{-1}$ which contains $(771.58\pm 10.57)\times 10^{6}~B\bar{B}$ pairs, collected with the Belle detector detector at the KEKB asymmetric-energy $e^{+}e^{-}$ collider kekb operating at the $\Upsilon(4S)$ resonance.

As the well established $B\to\psi(2S)K$ and $B\to X(3872)K$ decays have the same topology as the $B\to Y(4260)K$ decays, these decays are used as control samples to validate and calibrate Monte Carlo (MC) simulations. The signal simulation sample for each decay mode is generated using EvtGen evtgen . Here, the decays of $\psi(2S)$, $X(3872)$ and $J/\psi$ are specified to be $\psi(2S)\to J/\psi\pi^{+}\pi^{-}$, $X(3872)\to J/\psi\pi^{+}\pi^{-}$ and $J/\psi\to\ell^{+}\ell^{-}$, respectively, while $K_{S}^{0}$ decays generically pdg . All radiation effects are taken into account using PHOTOS photos . The detector response is simulated using GEANT3 geant .

The charged tracks used in the analysis are required to originate from the interaction point (IP) and have their point of closest approach to the IP within 3.5 cm along the beam axis and 1.0 cm in the plane transverse to the beam axis. Identification of charged pions and kaons are based on the information from the aerogel Cherenkov counter system, time-of-flight scintillation counter (TOF) and central drift chamber. All of the information is combined to form the pion (kaon) likelihood, $\mathcal{L}_{\pi}~(\mathcal{L}_{K})$, and the selections are made on the basis of the likelihood ratio $\mathcal{R}_{\pi(K)}={\mathcal{L}_{\pi(K)}}/({\mathcal{L}_{\pi}+\mathcal{L}_{K}})$. Charged pions (kaons) are identified requiring $\mathcal{R}_{\pi}~(\mathcal{R}_{K})>0.6$ with an identification efficiency of 94% (86%) and a misidentification rate of 7.5% (4%) for misidentifying kaon (pion) as a pion (kaon), respectively. These efficiencies and misidentification rates are determined using a control sample of $D^{*+}\to D^{0}(K^{-}\pi^{+})\pi^{+}$ decays in the kinematic region of interest.

A $K_{S}^{0}\to\pi^{+}\pi^{-}$ candidate decay is reconstructed from a pair of oppositely charged tracks with a $\pi^{+}\pi^{-}$ invariant mass in the range 488 MeV/$c^{2}<M_{\pi\pi}<$ 508 MeV/$c^{2}$ ($\pm 4\sigma$ around the nominal $K_{S}^{0}$ mass pdg ). The selected candidates are required to satisfy the criteria described in Ref. fangks .

Muon identification muon utilizes the track penetration depth and hit-distribution pattern in the $K_{L}^{0}$ and $\mu$ detector, which are combined to form the muon likelihood, $\mathcal{L}_{\mu}$, and the selection is made on the basis of the likelihood ratio $\mathcal{R}_{\mu}={\mathcal{L}_{\mu}}/({\mathcal{L}_{\mu}+\mathcal{L}_{\pi}+\mathcal{L}_{K}})$. Muons are identified requiring $\mathcal{R}_{\mu}>0.1$ with an identification efficiency of 93% and a misidentification rate of 3% for misidentifying a pion as a muon. Electron identification electron utilizes the electromagnetic shower shape and $E_{\rm ECL}/p$ ratio, where $E_{\rm ECL}$ is the energy deposition in electromagnetic calorimeter and $p$ is the track momentum, as well as the information used in the charged hadron identification, except that from the TOF. All the information is combined to form the electron likelihood ratio, $\mathcal{R}_{e}$. Electrons are identified requiring $\mathcal{R}_{e}>0.01$.

A $J/\psi$ candidate is reconstructed in its decay mode $J/\psi\to\ell^{+}\ell^{-}$, where $\ell$ stands for $e$ or $\mu$. In the $J/\psi\to e^{+}e^{-}$ mode, the energy loss due to bremsstrahlung photons is recovered by including the four-momenta of the photons detected within 0.05 radians around the electron or positron initial direction in the invariant mass calculation; this mode is, hereinafter, referred to as $J/\psi\to e^{+}e^{-}(\gamma)$. An invariant mass of a $J/\psi$ candidate is required to be in the range $3.05~{\rm GeV}/c^{2}\leq M_{ee(\gamma)}\leq 3.13~{\rm GeV}/c^{2}$ or $3.07~{\rm GeV}/c^{2}\leq M_{\mu\mu}\leq 3.13~{\rm GeV}/c^{2}$. The asymmetric interval is taken for $e^{+}e^{-}(\gamma)$ to include the radiative tail due to the imperfect energy loss recovery. A vertex- and mass-constrained fit is performed to the selected $J/\psi$ candidates in order to improve their momentum resolution.

The selected $J/\psi$ candidate is then combined with a $\pi^{+}\pi^{-}$ pair to form $\psi(2S)$, $X(3872)$, and $Y(4260)$ candidates, requiring the $J/\psi\pi^{+}\pi^{-}$ invariant mass, $M_{J/\psi\pi\pi}$, to be in the range $3.67~{\rm GeV}/c^{2}\leq M_{J/\psi\pi\pi}\leq 3.70~{\rm GeV}/c^{2}$, $3.835~{\rm GeV}/c^{2}\leq M_{J/\psi\pi\pi}\leq 3.910~{\rm GeV}/c^{2}$ and $4.0~{\rm GeV}/c^{2}\leq M_{J/\psi\pi\pi}\leq 4.6~{\rm GeV}/c^{2}$, respectively. To reconstruct a $B^{+}~(B^{0})$ candidate, a $K^{+}$ ($K_{S}^{0}$) candidate is combined with a $\psi(2S)$, $X(3872)$ or $Y(4260)$ candidate.

To identify the $B$ meson, two kinematic variables, the beam-constrained mass ($M_{\rm bc}=\sqrt{(E_{\rm beam}/c^{2})^{2}-\sum_{i}(p_{i}^{*}/c)^{2}}$) and the energy difference ($\Delta E=\sum_{i}E_{i}^{*}-E_{\rm beam}$), are used to discriminate the signal from the background. Here, $E_{\rm beam}$ is the beam energy and $p_{i}^{*}$ ($E_{i}^{*}$) is the momentum (energy) of the $i^{\rm th}$ final-state particle of the reconstructed signal candidate, where both are evaluated in the $e^{+}e^{-}$ center-of-mass (CM) frame. The $B$ candidates with $M_{\rm bc}>$ 5.27 GeV/$c^{2}$ are selected for further analysis.

Even after applying all the selection criteria, multiple $B$ candidates can be reconstructed from wrong combinations of the retained particles in an event. The mean number of $B$ candidates per event is found to be 1.6 (1.6), 1.7 (1.6) and 1.4 (1.2) for the charged (neutral) $B\to\psi(2S)(\to J/\psi\pi^{+}\pi^{-})K$, $B\to X(3872)(\to J/\psi\pi^{+}\pi^{-})K$ and $B\to Y(4260)K$ decays, respectively. In an event with multiple $B$ candidates, we select the best candidate that has the smallest value of $\chi^{2}_{\rm BCS}$ = $\chi^{2}_{\rm vtx}+\chi^{2}_{M_{\rm bc}}+\chi^{2}_{J/\psi}(+\chi^{2}_{K_{S}^{0}})$, where $\chi^{2}_{\rm vtx}$ represents the $\chi^{2}$ value obtained from a kinematic fit to the $B$ decay vertex for all the charged daughter particles, and the other $\chi^{2}$ values are evaluated using the reconstructed mass $M_{i}$ and its resolution $\sigma_{i}$ and the nominal mass $m_{i}^{\rm PDG}$ pdg of the reconstructed meson $i$ as $\chi^{2}_{i}=[(M_{i}-m_{i}^{\rm PDG})/\sigma_{i}]^{2}$. Here, beam-constrained $M_{\rm bc}$ is used for the reconstructed mass in $\chi_{M_{\rm bc}}^{2}$, and $\chi_{K_{S}^{0}}^{2}$ is used only for the neutral $B$ decays. The reconstructed mass resolutions $\sigma_{M_{\rm bc}}$, $\sigma_{J/\psi}$, and $\sigma_{K_{S}^{0}}$ are evaluated in the $B\to\psi(2S)K$ decays to be 2.6 ${\rm MeV}/c^{2}$, 9.8 ${\rm MeV}/c^{2}$ and 1.6 ${\rm MeV}/c^{2}$, respectively. According to MC simulations, the best candidate selection identifies the true signal at rates of 76% (72%) for the charged (neutral) $B\to Y(4260)K$ decays. The same best candidate selection criterion are applied in the reconstruction of the control sample decays.

The dominant background comes from $e^{+}e^{-}\to q\bar{q}$ ($q=u,d,s$ or $c$) continuum events. To suppress this background, we utilize the difference in event topology between the isotropic distribution of particles in $B\bar{B}$ events and the jet-like collimation of particles in $q\bar{q}$ events by placing a requirement on the ratio of the second- and zeroth-order Fox-Wolfram moments r2cut to be less than 0.5.

Among the backgrounds from $B\bar{B}$ events, the main contribution is expected to arise from inclusive $B$ decays to $J/\psi$. To understand possible backgrounds, simulated sample of inclusive $B$ decays with a $J/\psi~(\ell^{+}\ell^{-}$) in the final state are studied; the sample corresponds to an integrated luminosity that is two orders of magnitude larger than that of data. No peaking structures are found in the $M_{J/\psi\pi\pi}$ signal regions of $B\to\psi(2S)K$, $B\to X(3872)K$, and $B\to Y(4260)K$ decays. In order to check possible contributions from non-$J/\psi$ sources, the $J/\psi$ mass sidebands (2.54 GeV/$c^{2}<M_{J/\psi}<$ 2.72 GeV/$c^{2}$ and 3.32 GeV/$c^{2}<M_{J/\psi}<$ 3.50 GeV/$c^{2}$) are studied. The contributions are found to be negligible.

An unbinned extended maximum likelihood (UML) fit is performed to the $\Delta E$ distribution of each decay mode. The statistical weight for each candidate to be a signal decay is determined by using the ${}_{s}\mathcal{P}lot$ technique splot . The statistical weights can be used to effectively subtract the combinatorial background from the $M_{J/\psi\pi\pi}$ distribution of each decay mode. The signal yield of the intended resonance, then, can be extracted from the weighted $M_{J/\psi\pi\pi}$ distribution, having a single background component of the non-resonant $B\to J/\psi\pi^{+}\pi^{-}K$ decays.

The $\Delta E$ variable is required to satisfy $\mathbin{-}0.11~{\rm GeV}<\Delta E<0.11~{\rm GeV}$ for the $B\to\psi(2S)K$, $X(3872)K$ and $Y(4260)K$ decay modes. The UML function used here is

where $N$ is the total number of events, $N_{\rm S}$ ($N_{\rm B}$) is the number of signal (background) events, $P_{\rm S}$ ($P_{\rm B}$) is the signal (background) probability density function (PDF) of the variable $x$, and the index $i$ runs over the total number of events. Here, the signal refers to the charged or neutral $B\to J/\psi\pi^{+}\pi^{-}K$ decays, the background refers to the combinatorial background, and $x$ refers to the $\Delta E$ variable. The signal PDF is modeled by a sum of three Gaussians for the $B\to Y(4260)K$ decay modes and by a sum of two Gaussians and a bifurcated Gaussian for the $B\to\psi(2S)K$ and $B\to X(3872)K$ modes. The mean and resolution of the core Gaussian are allowed to vary in the fit while the remaining shape and normalization parameters are fixed to those obtained in the fit to the signal MC. The background PDF is modeled by a first-order polynomial except for the $B\to X(3872)K$ decay mode, in which a second-order polynomial is used. All parameters of the background PDF are allowed to vary in the fit.

The yields of the $B\to\psi(2S)K$, $X(3872)K$, and $Y(4260)K$ decays are extracted using independent UML fits to the weighted $M_{J/\psi\pi\pi}$ distributions. Here, while the functional form of Eq. 1 is used to evaluate the likelihood, the signal refers to the charged or neutral decay of $B\to\psi(2S)K$, $X(3872)K$, or $Y(4260)K$, the background refers to the corresponding non-resonant $B\to J/\psi\pi^{+}\pi^{-}K$ decay, and $x$ refers to the $M_{J/\psi\pi\pi}$ variable. The signal PDF is modeled by a sum of two Gaussians for the $B\to\psi(2S)K$ and $Y(4260)K$ decays while an additional bifurcated Gaussian is used for the $B\to X(3872)K$ decays. The core Gaussian parameters for the $B\to\psi(2S)K$ and $B^{+}\to X(3872)K^{+}$ decays are allowed to vary in the fit, while those for the $B^{0}\to X(3872)K^{0}$ and $B\to Y(4260)K$ decays are fixed to the values obtained in the fit to the signal MC and calibrated with data; the calibration is based on the comparison of the shape parameters between the data and simulation of the $B^{+}\to X(3872)K^{+}$ decay. All the remaining shape and normalization parameters of the signal PDF are fixed to those obtained in the fit to the signal MC. The background PDF is modeled by a first-order polynomial except for the $B\to\psi(2S)K$ decay modes, in which a second-order polynomial is used. All parameters of the background PDF are allowed to vary in the fit. The $\Delta E$ distributions, weighted $M_{J/\psi\pi\pi}$ distributions and projections of their PDFs obtained from the fits are shown in Fig. 1, 2, and 3 for the $B\to\psi(2S)K$, $X(3872)K$, and $Y(4260)K$ decay samples, respectively. The obtained signal yields of the $B\to\psi(2S)K$, and $B\to X(3872)K$, $X(3872)\to J/\psi\pi^{+}\pi^{-}$ decays are listed in Table 1 and for $B\to Y(4260)K$, $Y(4260)\to J/\psi\pi^{+}\pi^{-}$ decays are listed in Table 2.

For the $B\to Y(4260)K$ decays, the statistical significance of the signal yield is evaluated using the likelihood ratio as $\sqrt{-2\rm ln(\mathcal{L}_{0}/\mathcal{L}_{\rm max})}$, where $\mathcal{L}_{\rm max}$ and $\mathcal{L}_{0}$ denote the maximum likelihood of the nominal fit and that of the fit with the null signal hypothesis. The statistical significances are evaluated to be 2.9$\sigma$ and 1.4$\sigma$ for the charged and neutral $B\to Y(4260)K$ decays, respectively. The likelihood ratio is smeared with the systematic uncertainties, discussed later, and listed in Table 3. The signal significances taking into account the systematic uncertainties are determined to be 2.1$\sigma$ and 0.9$\sigma$ for the charged and neutral $B\to Y(4260)K$ decays, respectively.

The branching fractions ($\mathcal{B}$) of the $B\to\psi(2S)K$ decays are obtained as $\mathcal{B}=N_{\rm S}/[N_{B\bar{B}}\times\epsilon\times f_{K}\times\mathcal{B}(\psi(2S)\to J/\psi\pi^{+}\pi^{-})\times\mathcal{B}(J/\psi\to\ell^{+}\ell^{-})]$, where $N_{\rm S}$ is the number of signal decays, $N_{B\bar{B}}$ is the number of $B\bar{B}$ events in the data sample, and the branching fractions of the secondary decays are taken from Ref. pdg . Here, equal production of $B^{+}B^{-}$ and $B^{0}\bar{B}^{0}$ pairs from $\Upsilon(4S)$ decays is assumed. The reconstruction efficiency, $\epsilon$, is estimated from the signal MC simulation, with the application of calibrations to account for discrepancies between the data and signal MC related to particle identifications and $K^{0}_{S}$ reconstruction. These calibrations use dedicated control samples as discussed later. The coefficient $f_{K}$ is introduced to translate the branching fractions for the final states with $K^{0}_{S}$ into those for the ones with $K^{0}$ and set 1 and 0.5 for the charged and neutral $B\to\psi(2S)K$ decays, respectively. For the $B\to X(3872)K$ and $Y(4260)K$ decays, the branching fraction products are obtained in a similar manner as

where $R$ stands for the $X(3872)$ or $Y(4260)$ resonance. The obtained branching fractions of the $B\to\psi(2S)K$ decays and branching fraction products for the $B\to X(3872)K,~X(3872)\to J/\psi\pi^{+}\pi^{-}$ decays are listed in Table 1 with the associated reconstruction efficiencies and signal yields. The obtained values agree well with the world averages pdg and also with the previous Belle measurements data , indicating the validity of the signal extraction procedure. The branching fraction products of the $B\to Y(4260)K,~Y(4260)\to J/\psi\pi^{+}\pi^{-}$ decays as well as the associated reconstruction efficiencies and signal yields, are listed in Table 2.

With the absence of significant signals for the $B\to Y(4260)K$ decays, an upper limit (U.L.) is set on each signal yield at the 90% confidence level (C.L.) using a frequentist approach freq . The upper limits on the signal yields at the 90% C.L. ($N_{\rm S}^{UL}$) are found to be 259 and 84 events for the $B^{+}\to Y(4260)K^{+}$ and $B^{0}\to Y(4260)K^{0}_{S}$ decays, respectively. The upper limits on the branching fraction products are calculated using Eq.  Search for the $B\to Y(4260)K,~Y(4260)\to J/\psi\pi^{+}\pi^{-}$ decays, with $N_{\rm S}$ replaced by $N_{\rm S}^{UL}$ (systematic uncertainties are included in the upper limit calculation, as will be described later in this paper). The resulting upper limits are listed in Table 2.

In order to improve the signal sensitivity, a simultaneous fit to the charged and neutral signal decays is performed keeping the fit procedure same as in the nominal fits for the individual signal decays, except for incorporating the constraint that $\mathcal{B}(B^{+}\to Y(4260)K^{+})/\mathcal{B}(B^{0}\to Y(4260)K^{0})=\mathcal{B}(B^{+}\to\psi(2S)K^{+})/\mathcal{B}(B^{0}\to\psi(2S)K^{0})$ sys . The simultaneous fit for the $B\to Y(4260)K$ decays obtains $218\pm 68$ signal events, where the quoted uncertainty is statistical only. The combined statistical significance of the $B\to Y(4260)K,~Y(4260)\to J/\psi\pi^{+}\pi^{-}$ decays is found to be 3.2$\sigma$, which reduces to 2.2$\sigma$ once systematic uncertainties are taken into account. The simultaneous fit does not increase the significance of the $Y(4260)$ signal.

All the systematic uncertainties are summarized in Table 3. The tracking efficiency in MC simulation is calibrated using a control sample of $D^{*}\to\pi D^{0},~D^{0}\to\pi^{+}\pi^{-}K_{S}^{0},~K_{S}^{0}\to\pi^{+}\pi^{-}$ decays, and the uncertainty on the calibration factor is 0.35% per track. The calibration factor for the $K^{0}_{S}$ reconstruction efficiency is obtained using $D^{*\pm}\to D^{0}(\to K_{S}^{0}\pi^{0})\pi^{\pm}$ decays with an uncertainty of 0.7%. For the particle identification efficiencies, the calibration factors are obtained using the dedicated control samples mentioned earlier, and the resulting systematic uncertainty is 0.9% and 1.3% for kaon and pion identification, respectively. The dominant systematic uncertainties are due to the PDF modeling, and the values of the $Y(4260)$ mass and decay width pdg assumed in the fit. The changes on the signal yield from the nominal one due to the uncertainty in the PDF modeling is estimated by varying each of the fixed parameters independently by $\pm 1\sigma$. The corresponding changes due to the uncertainties on the $Y(4260)$ mass and decay width are estimated by separately applying the variation in the signal PDF based on the alternative signal MC simulations, which are generated varying each of the mass and decay width in the same manner. The resulting changes are added in quadrature. The uncertainty in the PDF modeling for the $B^{0}\to Y(4260)K^{0}_{S}$ decay gives an exceptionally large systematic uncertainty of 77.0%. This is due to the systematic uncertainty associated with the background PDF modeling. The fit procedures are validated in fully simulated MC experiments with ensembles of signal and inclusive $B$ decays involving $J/\psi$. The small biases of 4.3%–4.8% seen in the validation are taken as systematic uncertainties. The uncertainties on $N_{B\bar{B}}$ and $\mathcal{B}(J/\psi\to\ell^{+}\ell^{-}$), 1.4% and 0.4%, respectively, are also included in the systematic uncertainties. The total systematic uncertainties are estimated to be ${}_{-23.0}^{+30.5}$%, ${}_{-79.2}^{+17.5}$% and ${}_{-24.3}^{+26.2}$% on the results for the charged, neutral and combined $B\to Y(4260)K,~Y(4260)\to J/\psi\pi^{+}\pi^{-}$ decays, respectively, by adding all the sources in quadrature.

In summary, a search for the $B\to Y(4260)K$, $Y(4260)\to J/\psi\pi^{+}\pi^{-}$ decays is performed using $B\bar{B}$ pairs collected at the $\Upsilon(4S)$ resonance by the Belle experiment at the KEKB. The observed signal yields are $179\pm 53^{+55}_{-41}$ events and $39\pm 28^{+~7}_{-31}$ events for the charged and neutral $B\to Y(4260)K,~Y(4260)\to J/\psi\pi^{+}\pi^{-}$ decays, respectively, from fits to the individual decay samples; the first and second uncertainties are statistical and systematic, respectively. The signal significances are obtained to be $2.1\sigma$ and $0.9\sigma$ for the charged and neutral decays, respectively, taking into account the systematic uncertainties associated with the signal extraction. In the absence of any significant signals, the upper limits on the branching fraction products at the 90% C.L. are determined to be $1.4\times 10^{-5}$ and $1.7\times 10^{-5}$ for the charged and neutral decays, respectively, taking into account the systematic uncertainties.

The obtained results give the most stringent upper limits, to date, on the branching fraction products of the charged and neutral $B\to Y(4260)K,~Y(4260)\to J/\psi\pi^{+}\pi^{-}$ decays. The upper limits on the branching fraction products at the 95% C.L. are also determined and are $1.56\times 10^{-5}$ and $2.16\times 10^{-5}$ for the charged and neutral decays, respectively. The upper limit for the charged decay is consistent with the 95% confidence interval set by the BABAR collaboration babar and the one for the neutral decay is given for the first time. While an excess of events above background is seen, improved measurements with a larger data sample are demanded to establish signals and to elucidate the nature of the $Y(4260)$ state.

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 under Grant No. P 26794-N20; 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 Grant of the Russian Federation Government, Agreement No. 14.W03.31.0026; 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.