Search for a doubly-charged $DDK$ bound state in $\Upsilon(1S,2S)$ inclusive decays and via direct production in $e^+e^-$ collisions at $\sqrt{s}$ = 10.520, 10.580, and 10.867 GeV

We report the results of a first search for a doubly-charged $DDK$ bound state, denoted the $R^{++}$, in $\Upsilon(1S)$ and $\Upsilon(2S)$ inclusive decays and via direct production in $e^+e^-$ collisions at $\sqrt{s}$ = 10.520, 10.580, and 10.867 GeV. The search uses data accumulated with the Belle detector at the KEKB asymmetric-energy $e^+e^-$ collider. No significant signals are observed in the $D^{+}D_{s}^{*+}$ invariant-mass spectra of all studied modes. The 90\% credibility level upper limits on their product branching fractions in $\Upsilon(1S)$ and $\Upsilon(2S)$ inclusive decays (${\cal B}(\Upsilon(1S,2S) \to R^{++} + anything) \times {\cal B}(R^{++} \to D^{+} D_{s}^{*+})$), and the product values of Born cross section and branching fraction in $e^+e^-$ collisions ($\sigma(e^+e^- \to R^{++} + anything) \times {\cal B}(R^{++} \to D^{+} D_{s}^{*+})$) at $\sqrt{s}$ = 10.520, 10.580, and 10.867 GeV under different assumptions of $R^{++}$ masses varying from 4.13 to 4.17 GeV/$c^2$, and widths varying from 0 to 5 MeV are obtained.

We report the results of a first search for a doubly-charged DDK bound state, denoted the R ++ , in Υ(1S) and Υ(2S) inclusive decays and via direct production in e + e − collisions at √ s = 10.520, 10.580, and 10.867 GeV. The search uses data accumulated with the Belle detector at the KEKB asymmetric-energy e + e − collider. No significant signals are observed in the D + D * + s invariant-mass spectra of all studied modes. The 90% credibility level upper limits on their product branching fractions in Υ(1S) and Υ(2S) inclusive decays (B(Υ(1S, 2S) → R ++ + anything) × B(R ++ → D + D * + s )), and the product values of Born cross section and branching fraction in e + e − collisions (σ(e + e − → R ++ + anything) × B(R ++ → D + D * + s )) at √ s = 10.520, 10.580, and 10.867 GeV under different assumptions of R ++ masses varying from 4.13 to 4.17 GeV/c 2 , and widths varying from 0 to 5 MeV are obtained.
By exchanging a kaon, a D + D * s0 (2317) + molecular state can be formed with a binding energy of 5 − 15 MeV, regardless of whether the D * s0 (2317) + is treated as a cs state or a DK molecule [35]. In Ref. [36], the authors studied the DDK system in a coupledchannel approach, where an isospin 1/2 state, denoted the R ++ , is formed at 4140 MeV/c 2 when the D * s0 (2317) + is generated from the DK subsystem. The R ++ can be interpreted as a D + D * s0 (2317) + molecule-like state with exotic properties: doubly-charged and doubly-charmed. Hereinafter we also refer to this predicted state as R ++ .
An R ++ with the properties described above would be able to decay via R ++ → D + D * s0 (2317) + , where D * s0 (2317) + → D + s π 0 is an isospin-violating process. The alternative processes are via triangle diagrams into R ++ → D + D * + s and R ++ → D + s D * + [36,37]. The mass of R ++ is predicted to be in the range of 4.13 to 4.17 GeV/c 2 [37].
A D + D * + s final state could also arise from the decay of a doubly-charged tetraquark [38]. The question whether QQqq tetraquarks with two heavy quarks Q and two light antiquarksq are stable or unstable against decay * now at University of Hiroshima into two Qq mesons has a long history. It has been largely undecided, mainly due to a lack of experimental information about the strength of the interaction between two heavy quarks. The discovery of the doubly-charmed baryon Ξ ++ cc by LHCb [39] has provided the crucial experimental input. With the same method based on the quark model to calculate the Ξ ++ cc mass, a doublycharged tetraquark with a mass of 4156 MeV/c 2 , a J P of 1 + and a final state of D + D * + s has been predicted [38]. Thus the D + D * + s final state is a good channel to search for such a tetraquark state.
In this paper, we search for a doubly-charged DDK bound state in the D + D * + s final state in Υ(1S) and Υ(2S) inclusive decays, and via direct production in e + e − collisions at √ s = 10.520, 10.580, and 10.867 GeV. Inclusion of charge-conjugate decays is implicitly assumed throughout this analysis. We report a search for the R ++ with masses varying from 4.13 to 4.17 GeV/c 2 and widths varying from 0 to 5 MeV.

II. THE DATA SAMPLE AND THE BELLE DETECTOR
This analysis utilizes (5.74 ± 0.09) fb −1 of data collected at the Υ(1S) peak [(102 ± 3) million Υ(1S) events], (24.91 ± 0.35) fb −1 of data collected at the Υ(2S) peak [(158 ± 4) million Υ(2S) events], a data sample of (89.5 ± 1.3) fb −1 collected at √ s = 10.520 GeV, a data sample of (711.0 ± 10.0) fb −1 collected at √ s = 10.580 GeV [Υ(4S) peak], and a data sample of (121.4 ± 1.7) fb −1 collected at √ s = 10.867 GeV [Υ(5S) peak]. All the data were collected with the Belle detector [40] operating at the KEKB asymmetric-energy e + e − collider [41]. The Belle detector is described in detail in Ref. [40]. It is a large-solid-angle magnetic spectrometer consisting of a silicon vertex detector, a 50-layer central drift chamber (CDC), an array of aerogel threshold Cherenkov counters (ACC), a barrellike arrangement of time-of-flight scintillation counters (TOF), and an electromagnetic calorimeter comprising CsI(TI) crystals (ECL) located inside a superconducting solenoid coil that provides a 1.5 T magnetic field. An iron flux return placed outside the coil is instrumented to detect K 0 L mesons and to identify muons. Monte Carlo (MC) signal samples are generated with EvtGen [42] to determine signal shapes and efficiencies. Initial-state radiation (ISR) is taken into account by assuming that the cross sections follow a 1/s dependence in e + e − → R ++ + anything reactions, where s is the center-of-mass energy squared. The mass of R ++ is varied from 4.13 to 4.17 GeV/c 2 in steps of 0.01 GeV/c 2 , and with a width varying from 0 to 5 MeV in steps of 1 MeV. These events are processed by a detector simulation based on geant3 [43].
s , and e + e − → qq (q = u, d, s, c) at √ s = 10.520, 10.580, and 10.867 GeV corresponding to four times the integrated luminosity of data are used to study possible peaking backgrounds.

III. COMMON EVENT SELECTION CRITERIA
For well-reconstructed charged tracks, except those from K 0 S → π + π − decays, the impact parameters perpendicular to and along the beam direction with respect to the nominal interaction point (IP) are required to be less than 0.5 cm and 2 cm, respectively, and the transverse momentum in the laboratory frame is required to be larger than 0.1 GeV/c. For the particle identification (PID) of a well-reconstructed charged track, information from different detector subsystems, including specific ionization in the CDC, time measurement in the TOF, and the response of the ACC, is combined to form a likelihood L i [44] for particle species i, where i = π or K. Tracks with R K = L K /(L K +L π ) < 0.4 are identified as pions with an efficiency of 96%, while 5% of kaons are misidentified as pions; tracks with R K > 0.6 are identified as kaons with an efficiency of 95%, while 4% of pions are misidentified as kaons. Except for tracks from K 0 S decays, all charged tracks are required to be positively identified by the above procedures.
An ECL cluster is taken as a photon candidate if it does not match the extrapolation of any charged track. The energy of the photon is required to be greater than 50 MeV.
The K 0 S candidates are first reconstructed from pairs of oppositely charged tracks, which are treated as pions, with a production vertex significantly separated from the average IP, then selected using an artificial neural network [45] based on two sets of input variables [46]. The φ andK * (892) 0 candidates are reconstructed using K + K − and K − π + decay modes, respectively. The invariant masses of the K 0 S and φ candidates are required to be within 7 MeV/c 2 of the corresponding nominal masses (> 90% signal events are retained).
We reconstruct D + mesons in the K − π + π + and K 0 S (→ π + π − )π + decay channels, and D + s mesons in the φπ + andK * (892) 0 K + decay channels. We perform vertexand mass-constrained fits for D + and D + s candidates, and require χ 2 vertex /n.d.f. < 20 (> 97% selection efficiency according to MC simulation). The selected D + s candidate is combined with a photon to form a D * + s candidate, and a mass-constrained fit is performed to improve its momentum resolution.
The signal mass windows forK * (892) 0 , D + , D + s , and D * + s candidates have been optimized by maximizing the Punzi parameter S/(3/2 + √ B) [47]. Here, S is the number of R ++ signal events in the MC-simulated Υ(2S) → R ++ + anything sample with the mass and width of R ++ fixed at 4.13 GeV/c 2 and 2 MeV assuming B(Υ(2S) → R ++ + anything) × B(R ++ → D + D * + s ) = 10 −4 , and B is the number of background events in the R ++ signal window. The number of background events is obtained from the normalized M D + and M D * + s sidebands in the data requiring 4.12 GeV/c 2 ¡M D + D * + s ¡4.14 GeV/c 2 as the R ++ signal region (about 3σ according to signal MC simulations). The optimized signal regions are Finally, when the D + and D * + s candidates are combined to form R ++ candidates, all the combinations are preserved for further analysis. The fraction of events where multiple combinations are selected as R ++ candidates is 14% in data, which is consistent with the MC simulation. In the fits, the D + s and D + signal shapes are described by double-Gaussian functions, and the D * + s signal shape is described by a Novosibirsk function [49], where the values of parameters are fixed to those obtained from the fits to the corresponding signal MC distributions. The backgrounds are parametrized by first-order polynomial functions for D + s and D + , and a second-order polynomial function for D * + s . Figure [50], so first-order polynomial functions with free parameters are taken as background shapes. The fitted results with the R ++ mass fixed at 4.14 GeV/c 2 and width fixed at 2 MeV are shown in Fig. 4 as an example. The local R ++ significance is calculated using −2 ln(L 0 /L max ), where L 0 and L max are the likelihoods of the fits without and with a signal component, respectively. The fitted R ++ signal yields with the M R ++ fixed at values ranging from 4.13 to 4.17 GeV/c 2 in steps of 0.01 GeV/c 2 and Γ R ++ fixed at values ranging from 0 to 5 MeV in steps of 1 MeV and the corresponding statistical significances are listed in Table I.
where N fit is the fitted number of signal events, N Υ(1S) = 1.02×10 8 and N Υ(2S) = 1.58×10 8 are the total numbers of Υ(1S) and Υ(2S) events, the index i runs for all finalstate modes with ε i being the corresponding efficiency and B i the product of all secondary branching fractions of the mode    Table I.
Since the statistical significance in each case is less than 3σ, Bayesian upper limits at the 90% credibility level (C.L.) on the numbers of signal events (N UL ) assuming it follows a Poisson distribution with a uniform prior probability density function are determined by solving the equation x is the number of fitted signal events and L(x) is the likelihood function in the fit to data. Taking into account the systematic uncertainties discussed below, the likelihood curve is convolved with a Gaussian function whose width equals the corresponding total multiplicative systematic uncertainty. The calculated        Fig. 10 as an example. The local R ++ significance is calculated using the same method as described in Sec. IV. The fitted R ++ signal yields with M R ++ fixed at values ranging from 4.13 to 4.17 GeV/c 2 in steps of 0.01 GeV/c 2 and Γ R ++ fixed at values ranging from 0 to 5 MeV in steps of 1 MeV and the corresponding statistical significances are listed in Table II. The product of Born cross section and branching fraction σ(e + e − → R ++ + anything) × B(R ++ → D + D * + s ) is calculated from the following formula: where N fit is the number of fitted signal yields in data, |1 − | 2 is the vacuum polarization factor, L is the integrated luminosity, the index i runs for all final-state modes with ε i being the corresponding efficiency and  Table II. Since the statistical significance in each case is less than 3σ, Bayesian upper limits at the 90% C.L. on N UL are obtained using the same method as described in Sec. IV. The results for N UL and product values of Born cross section and branching fraction (σ UL (e + e − → R ++ + anything) × B(R ++ → D + D * + s )) in e + e − collisions at √ s = 10.520, 10.580, and 10.867 GeV with M R ++ fixed at values ranging from 4.13 to 4.17 GeV/c 2 and Γ R ++ fixed at values ranging from 0 to 5 MeV are listed in Table II

VI. SYSTEMATIC UNCERTAINTIES
The systematic uncertainties in the branching fraction and Born cross section measurements can be divided into two categories: multiplicative systematic uncertainties and additive systematic uncertainties.
The sources of multiplicative systematic uncertainties include detection-efficiency-related uncertainties, the statistical uncertainty of the MC efficiency, the modeling of MC event generation, branching fractions of intermediate states, energy dependence of the cross sections, the total numbers of Υ(1S, 2S) events as well as the integrated luminosity.
The detection-efficiency-related uncertainties include those for tracking efficiency (0.35% per track), particle identification efficiency (1.8% per kaon, 1.0% per pion), as well as momentum-weighted K 0 S selection efficiency (2.2%) [53]. The photon reconstruction contributes 2.0% per photon, as determined from radiative Bhabha events. The above individual uncertainties from different reconstructed modes are added linearly, weighted by the product of the detection efficiency and the product of all secondary branching fractions (ǫ i × B i ). Assuming these uncertainties are independent and adding them in quadrature, the final uncertainty related to the reconstruction efficiency is 6.6%.
The MC statistical uncertainties are estimated using the yields of selected and generated events; these are 1.0% or less.
We use the EvtGen generator to generate the signal MC samples. By changing the recoil mass of the R ++ , the efficiencies are changed by (1 − 3)%. To be conservative, we take 1% and 3% as the systematical uncertainties related to signal MC statistics and generation.
The relative uncertainties of branching fractions for D + → K − π + π + , D + → K 0 S π + , K 0 S → π + π − ,  under different assumptions of R ++ mass (M R ++ in GeV/c 2 ) and width (Γ R ++ in MeV), where N fit is the number of fitted signal events, N UL is the 90% C.L. upper limit on the number of signal events taking into account systematic uncertainties, Σ(σ) is the local R ++ significance, Σi(ǫiBi) is the sum of product of the detection efficiency and the product of all secondary branching fractions for each reconstruction mode, σ multi is the total multiplicative systematic uncertainty, σ add is the additive systematic uncertainty, σ × B (σ(e + e − → R ++ + anything) × B(R ++ → D + D * + s )) is the product value of Born cross section and branching fraction, and σ UL × B (σ UL (e + e − → R ++ + anything) × B(R ++ → D + D * + s )) is the 90% C.L. upper limit on the product value of Born cross section and branching fraction with systematic uncertainties included.
e + e − → R ++ + anything at √ s = 10.520/10.580/10.867 GeV,   D * + s → γD + s , D + s → φ(→ K + K − )π + , and D + s → K * (892) 0 (→ K − π + )K + are taken from Ref. [48] and summed in quadrature to obtain the total uncertainty of the branching fractions of the intermediate states for each reconstructed mode. The above individual uncertainties from different reconstructed modes are added linearly with a weighting factor of ǫ i × B i to obtain 2.5% as the uncertainty due to the branching fractions of intermediate states.
Changing the s dependence of the cross sections of e + e − → R ++ + anything from 1/s to 1/s 4 , the radiative correction factors (1 + δ) ISR become 0.712, 0.711, and 0.709 for √ s = 10.520, 10.580, and 10.867 GeV, respectively. The differences are less than 0.3%. Thus, the systematic uncertainty related to the radiative correction factors is negligible with respect to the other sources.
The uncertainties on the total numbers of Υ(1S) and Υ(2S) events are 2.0% and 2.3%, respectively, which are mainly due to imperfect simulations of the charged track multiplicity distributions from inclusive hadronic MC events. The total luminosity is determined to 1.4% precision using wide-angle Bhabha scattering events.
Additive systematic uncertainties due to the mass resolution and fit are considered as follows.
The uncertainty due to the mass resolution is studied by using the control sample of B 0 → D − D * + s ; the difference in mass resolution between MC simulation and data is around 10%. Thus, the uncertainty due to the mass resolution is estimated by enlarging the mass resolution by 10% when fitting the D + D * + s invariantmass distributions.
To estimate the uncertainties associated with the fit, the order of the background polynomial is changed from first to second or third and the range of the fit is changed by ±30 MeV/c 2 .
The upper limits on the branching fraction and Born cross section at the 90% C.L. are determined and the systematic uncertainties are taken into account in two steps.
First, when we study the additive systematic uncertainties described above, we take the most conservative upper limit at the 90% C.L. on the number of R ++ signal yields. The differences between the most conservative upper limits and the nominal fits are in the range of 2.3% − 16.4% (see Tables I  and II for detailed vaules), depending on the centerof-mass energy, the mass and width of the R ++ state. Then, to take into account the multiplicative systematic uncertainties, the likelihood with the most conservative upper limit is convolved with a Gaussian function    whose width is the corresponding total multiplicative systematic uncertainty.
The sources of uncertainties are assumed independent, and the total multiplicative systematic uncertainties are obtained by adding all uncertainties in quadrature. The total multiplicative systematic uncertainties are listed in Tables I and II for