Search for tetraquark states $X_{cc\bar{s}\bar{s}}$ in $D_{s}^{+}D_{s}^{+}~{}(D_{s}^{*+}D_{s}^{*+})$ final states at Belle

The hadron spectrum was successfully categorized based on the quark model as early as the 1960s quarkmodel . For a long time, all known hadrons could be classified as mesons or baryons with components of a quark-antiquark pair ($q\bar{q}$) or three quarks ($qqq$), respectively. However, Quantum Chromodynamics (QCD) also allows the existence of more complex structures, such as the tetraquark, pentaquark, or glueball, which possess properties that are forbidden for conventional hadrons. The states that do not fit into the ordinary $q\bar{q}$ or $qqq$ scheme in the quark model are referred to as exotic states.

The experimental discovery of exotic states began in 2003 with the observation of the $X(3872)$ 3872-1 . This new state did not fit any ordinary $c\bar{c}$ quarkonia in the quark model. After that, the $X(3872)$ was observed in multiple decay modes and confirmed by various experiments 3872-2 ; 3872-3 ; 3872-4 . Many different theoretical interpretations of this state have been proposed, such as meson molecule, tetraquark, and conventional bound state 3872-theory-1 ; 3872-theory-2 ; 3872-theory-3 ; 3872-theory-4 ; 3872-theory-5 . During the past two decades, there has been considerable world-wide activity in exotic state research using various processes, such as $e^{+}e^{-}$ annihilation (e.g., at $\tau$-charm facilities and B-factories), hadron collisions (e.g., at the Tevatron and the LHC), or photo- and leptoproduction (e.g., at the SPS, HERA or at Jefferson Lab), and many exotic state candidates were observed xyz-yuan ; xyz-shen .

In searches for exotic states, a clear feature that helps distinguish exotic from ordinary hadrons would be a nonzero electric charge in a state which contains a heavy quark-antiquark pair of the same flavor. Such a state must contain at least one more quark-antiquark pair, and is thus not a conventional quark-antiquark meson. Furthermore, a state with a pair of two identical heavy flavor quarks (for example, $cc$), has even more pronounced features as an exotic state. Very recently, the LHCb experiment announced observation of an open-double-charm state $T_{cc}^{+}$ in the $D^{0}D^{0}\pi^{+}$ mass spectrum near threshold tcc-lhcb-1 ; tcc-lhcb-2 . It contains two charm quarks and two light quarks, thus it is a clear evidence for an exotic state. On the theoretical side, in addition to tetraquark models based on a heavy quark pair and two light quarks, the double-heavy tetraquark states are studied using QCD sum rules the_qcd , quark models the_qm1 ; the_qm2 , and lattice QCD computations the_lattqcd . Besides, a QCD-inspired chiral quark model gives a prediction on the tetraquark states denoted as $X_{cc\bar{s}\bar{s}}$ with $+2$ electric charge in spin-parity channels $J^{P}=0^{+}$ and $2^{+}$, which are expected to be found in $D_{s}^{+}D_{s}^{+}$ and $D_{s}^{*+}D_{s}^{*+}$ final states the_qqss2 . The predicted masses and widths of those resonances are listed in Table 1. Among the three predicted resonances in $D_{s}^{*+}D_{s}^{*+}$ final state, the narrowest one has the highest observable probability.

In this paper, we present a search for double-heavy tetraquark candidates using the $D_{s}^{+}D_{s}^{+}$ and $D_{s}^{*+}D_{s}^{*+}$ final states in $\Upsilon(1S,2S)$ inclusive decays, and $e^{+}e^{-}\to D_{s}^{+}D_{s}^{+}(D_{s}^{*+}D_{s}^{*+})+anything$ processes at $\sqrt{s}$ = 10.52, 10.58, and 10.867 GeV. The $D_{s}^{*+}$ candidates are reconstructed in decays to $D_{s}^{+}\gamma$, while the $D_{s}^{+}$ candidates are reconstructed in the $D_{s}^{+}\to\phi(\to K^{+}K^{-})\pi^{+}$ and $\bar{K}^{*}(892)^{0}(\to K^{-}\pi^{+})K^{+}$ decays. Inclusion of charged-conjugate modes is implicitly assumed throughout this analysis.

The data samples used in this analysis include: a 5.74 fb^-1 data sample collected at the $\Upsilon(1S)$ peak (102 million $\Upsilon(1S)$ events); a 24.7 fb^-1 data sample collected at the $\Upsilon(2S)$ peak (158 million $\Upsilon(2S)$ events); an 89.5 fb^-1 data sample collected at $\sqrt{s}$ = 10.52 GeV; a 711 fb^-1 data sample collected at $\sqrt{s}$ = 10.58 GeV, and a 121.4 fb^-1 data sample collected at $\sqrt{s}$ = 10.867 GeV, where $s$ is the center-of-mass energy squared. All the data were collected with the Belle detector, which is described in detail in Ref. detector , operating at the KEKB asymmetric-energy $e^{+}e^{-}$ collider collider . 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 barrel-like arrangement of time-of-flight scintillation counters (TOF), and an electromagnetic calorimeter comprising CsI(Tl) crystals (ECL) located inside a superconducting solenoid coil that provides a $1.5~{}\hbox{T}$ magnetic field. An iron flux return comprising resistive plate chambers placed outside the coil was instrumented to detect $K^{0}_{L}$ mesons and to identify muons.

Monte Carlo (MC) signal events are generated with EvtGen evtgen and processed through a full simulation of the Belle detector based on GEANT3 geant . Initial-state radiation (ISR) is taken into account assuming that the cross sections follow a $1/s$ dependence in $e^{+}e^{-}\to X_{cc\bar{s}\bar{s}}+anything$ reactions. The processes $\Upsilon(1S,2S)\to D_{s}^{+}D_{s}^{+}(D_{s}^{*+}D_{s}^{*+})+anything$ and $e^{+}e^{-}\to D_{s}^{+}D_{s}^{+}(D_{s}^{*+}D_{s}^{*+})+anything$ at $\sqrt{s}$ = 10.52 GeV, 10.58 GeV, and 10.867 GeV are taken into account, where the $D_{s}^{*+}$ decays into $D_{s}^{+}\gamma$ using a $P$-wave model, and the $D_{s}^{+}$ decays to $K^{+}K^{-}\pi^{+}$ final states using a Dalitz plot decay model of Ref. cleo-dalitz . The mass of $X_{cc\bar{s}\bar{s}}$ is chosen in the interval from 4882 MeV/$c^{2}$ to 4922 MeV/$c^{2}$ (4801 MeV/$c^{2}$ to 4841 MeV/$c^{2}$) in steps of 5 MeV/$c^{2}$, with a width varying from 0.54 MeV to 6.54 MeV (2.58 MeV to 8.58 MeV) in steps of 1 MeV for $X_{cc\bar{s}\bar{s}}\to D_{s}^{+}D_{s}^{+}$ ($D_{s}^{*+}D_{s}^{*+}$). Inclusive MC samples of $\Upsilon(1S,2S)$ decays, $\Upsilon(4S)\to B^{+}B^{-}/B^{0}\bar{B}^{0}$, $\Upsilon(5S)\to B_{s}^{(*)}\bar{B}_{s}^{(*)}$, and $e^{+}e^{-}\to q\bar{q}$ $(q=u,d,s,c)$ at $\sqrt{s}$ = 10.52 GeV, 10.58 GeV, and 10.867 GeV corresponding to four times the integrated luminosity of data are used to study possible peaking backgrounds.

For reconstructed charged tracks, the impact parameters perpendicular to and along the beam direction with respect to the interaction point (IP) are required to be less than 0.2 cm and 1.5 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 ${\mathcal{L}}_{i}$ pidcode for particle species $i$, where $i$ = $\pi$ or $K$. Tracks with $R_{K}=\mathcal{L}_{\textrm{K}}/(\mathcal{L}_{K}+\mathcal{L}_{\pi})<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.

An ECL cluster is taken as a photon candidate if it does not match the extrapolation of any charged tracks. The energy of the photon candidate from the $D_{s}^{*+}$ decay is required to be greater than 50 MeV. For $D_{s}^{+}$ candidates, vertex and mass-constrained fits are performed, and then $\chi^{2}_{\textrm{vertex}}/n.d.f<20$ is required ($>97\%$ selection efficiency according to MC simulation). For $D_{s}^{*+}$ candidates, a mass-constrained fit is performed to improve its momentum resolution. The best $D_{s}^{*+}$ candidate with $\chi^{2}$ of $D_{s}^{*+}$ mass-constrained fit for each $D_{s}^{+}$ candidate is kept to suppress the combinational background.

The signal mass windows for $\bar{K}^{*}(892)^{0}$, $\phi$, $D_{s}^{+}$, and $D_{s}^{*+}$ candidates have been optimized by maximizing the Punzi parameter $S/(3/2+\sqrt{B})$ fom , where $S$ is the number of selected events in the simulated signal process by fitting the $X_{cc\bar{s}\bar{s}}$ invariant-mass spectrum. $B$ is the number of selected events obtained from the normalized $M_{D_{s}^{+}D_{s}^{+}}$ sidebands in inclusive MC samples. The optimized mass window requirements are $|M_{K^{+}K^{-}}-m_{\phi}|<8$ MeV/$c^{2}$, $|M_{\phi\pi^{+}}-m_{D_{s}^{+}}|<7$ MeV/$c^{2}$, $|M_{K^{-}\pi^{+}}-m_{\bar{K}^{*}(892)^{0}}|<50$ MeV/$c^{2}$, $|M_{\bar{K}^{*}(892)^{0}K^{+}}-m_{D_{s}^{+}}|<7$ MeV/$c^{2}$, and $|M_{\gamma D_{s}^{+}}-m_{D_{s}^{*+}}|<14$ MeV/$c^{2}$, where $m_{\phi}$, $m_{\bar{K}^{*}(892)^{0}}$, $m_{{D}_{s}^{+}}$, and $m_{{D}_{s}^{*+}}$ are the nominal masses of $\phi$, $\bar{K}^{*}(892)^{0}$, ${D}_{s}^{+}$, and ${D}_{s}^{*+}$ PDG . There are no multiple candidates after processing all selections in both $D_{s}^{+}D_{s}^{+}$ and $D_{s}^{*+}D_{s}^{*+}$ cases. Figure 1 shows the scatter plots of $D_{s}^{+}$ versus $D_{s}^{+}$ invariant masses from the selected $e^{+}e^{-}\to X_{cc\bar{s}\bar{s}}(\to D_{s}^{+}D_{s}^{+}(D_{s}^{*+}D_{s}^{*+}))+anything$ candidates from data at $\sqrt{s}$ = 10.58 GeV as an example. Here we define the two-dimensional $D_{s}^{+}D_{s}^{+}$ sidebands, and the normalized contribution from $D_{s}^{+}$ and $D_{s}^{+}$ sidebands is estimated using 25% of the number of events in the blue dashed line boxes and reduced by 6.25% of the number of events in the red dotted line boxes.

The $D_{s}^{+}D_{s}^{+}$ and $D_{s}^{*+}D_{s}^{*+}$ invariant mass distributions of selected events from data samples in the kinematically allowed region are shown in Figs. 2 and  3 together with the backgrounds estimated from the normalized $D_{s}^{+}D_{s}^{+}$ sideband events. No peaking backgrounds are found in the normalized sideband events in either $D_{s}^{+}D_{s}^{+}$ and $D_{s}^{*+}D_{s}^{*+}$ invariant mass distributions from data, nor in the $D_{s}^{+}D_{s}^{+}$ and $D_{s}^{*+}D_{s}^{*+}$ mass spectra from inclusive MC samples topoo . Thus in the following we only focus on the mass spectra from the theoretically predicted regions for $X_{cc\bar{s}\bar{s}}$ the_qqss2 which are shown in Figs. 4 and  5.

Since no clear signals are observed in the invariant-mass spectra, the 90% confidence level (C.L.) upper limits on the numbers of signal events are given. The upper limit is calculated by the frequentist approach pole-1 implemented in the POLE (Poissonian limit estimator) program pole-2 , where the mass window is obtained by giving 95% acceptance to the corresponding simulated signal events, the number of signal candidate events is counted directly, and the number of expected background events is estimated from the normalized mass sidebands. The possible non-resonant contributions in the $D_{s}^{+}D_{s}^{+}$ and $D_{s}^{*+}D_{s}^{*+}$ invariant-mass spectra are not subtracted and taken as potential signals, in order to set more conservative upper limits.

The upper limit calculation is repeated with $M_{X_{cc\bar{s}\bar{s}}}$ varying from 4882 MeV/$c^{2}$ to 4922 MeV/$c^{2}$ in steps of 5 MeV/$c^{2}$ and $\Gamma_{X_{cc\bar{s}\bar{s}}}$ varying from 0.54 MeV to 6.54 MeV in steps of 1.0 MeV for the $M_{D_{s}^{+}D_{s}^{+}}$ distribution, and with $M_{X_{cc\bar{s}\bar{s}}}$ varying from 4801 MeV/$c^{2}$ to 4841 MeV/$c^{2}$ in steps of 5 MeV/$c^{2}$ and $\Gamma_{X_{cc\bar{s}\bar{s}}}$ varying from 2.58 MeV to 8.58 MeV in steps of 1.0 MeV for the $M_{D_{s}^{*+}D_{s}^{*+}}$ distribution.

There are several sources of systematic uncertainties on the branching fraction and Born cross section measurements, which can be divided into multiplicative and additive systematic uncertainties. The multiplicative systematic uncertainties include detection-efficiency-related (DER) sources (tracking efficiency, PID, and photon reconstruction), the statistical uncertainty of the MC efficiency, branching fractions of intermediate states, the total numbers of $\Upsilon(1S)$ and $\Upsilon(2S)$ events, and the integrated luminosities at $\sqrt{s}$ = 10.52 GeV, 10.58 GeV, and 10.867 GeV.

The systematic uncertainties related to detection efficiency ($\sigma_{\textrm{DER}}$) include the tracking efficiency (0.35% per track, estimated using partially reconstructed $D^{\ast}$ decays in $D^{*+}\to\pi^{+}D^{0},D^{0}\to K_{S}^{0}\pi^{+}\pi^{-}$), PID efficiency ($2.2\%$ per kaon and $1.8\%$ per pion, estimated using $D^{*+}\to D^{0}\pi^{+}$, $D^{0}\to K^{-}\pi^{+}$ samples), and photon reconstruction (2.0% per photon, estimated using a radiative Bhabha sample). The statistical uncertainty in the signal MC simulation efficiency can be calculated as $\Delta\varepsilon$ = $\sqrt{\varepsilon(1-\varepsilon)/N}$, where $\varepsilon$ is the reconstruction efficiency after all event selections, and $N$ is the total number of generated events. Its relative uncertainty $\sigma_{\textrm{MC stat.}}=\Delta\varepsilon/\varepsilon$ is at most at the 1.0% level. Changing the $s$ dependence of the cross sections of $e^{+}e^{-}\to X_{cc\bar{s}\bar{s}}(\to D_{s}^{*+}D_{s}^{*+})+anything$ from $1/s$ to $1/s^{4}$, the product of efficiency and radiative correction factor $\epsilon(1+\delta)_{\textrm{ISR}}$ changes by less than 0.3% ($\sigma_{\textrm{ISR}}$).

The relative uncertainties of branching fractions for $D_{s}^{*+}\to\gamma D_{s}^{+}$, $D_{s}^{+}\to\phi(\to K^{+}K^{-})\pi^{+}$, and $D_{s}^{+}\to\bar{K}^{*}(892)^{0}(\to K^{-}\pi^{+})K^{+}$ are 0.75%, 3.52%, and 3.45% PDG , respectively. The total uncertainties are calculated using $\sigma_{{\cal B}}=\frac{\sqrt{\Sigma{(\varepsilon_{i}\times{\cal B}_{i}\times\sigma_{{\cal B}_{i}})^{2}}}}{\Sigma{(\varepsilon_{i}\times{\cal B}_{i}})}$, where $\varepsilon_{i}$ is the efficiency, $\sigma_{{\cal B}_{i}}$ is the relative uncertainty of intermediate states’ branching fractions, and ${\cal B}_{i}$ is the product of branching fractions of the intermediate states for each reconstructed mode $i$.

The total numbers of $\Upsilon(1S)$ and $\Upsilon(2S)$ events are estimated to be ($102\pm 2)\times 10^{6}$ and ($157.8\pm 3.6)\times 10^{6}$, which are determined by counting the numbers of inclusive hadrons. The uncertainties are mainly due to imperfect simulations of the charged multiplicity distributions from inclusive hadronic MC events ($\sigma_{\textrm{N}_{\Upsilon(1S,2S)}}$). Belle measures luminosity with 1.4% precision using wide angle Bhabha events ($\sigma_{{\cal L}}$).

All the multiplicative uncertainties are summarized in Table 2 for the measurements of $\Upsilon(1S,2S)\to X_{cc\bar{s}\bar{s}}+anything$ and $e^{+}e^{-}\to X_{cc\bar{s}\bar{s}}+anything$ at $\sqrt{s}$ = 10.52 GeV, 10.58 GeV, and 10.867 GeV, respectively. The total multiplicative uncertainty is calculated by adding all sources of multiplicative uncertainty in quadrature,

The additive uncertainty due to the number of expected background is considered by counting normalized background distributions directly, fitting the distributions with a constant, and a 1st-order polynominal.

Since no signal traces are observed in the $D_{s}^{+}D_{s}^{+}$ or $D_{s}^{*+}D_{s}^{*+}$ distributions from data at all energy points, the 90% C.L. upper limits on the numbers of signal events ($N^{\textrm{UP}}$) are determined. To take into account the additive and multiplicative uncertainties, we first study the additive systematic uncertainty and take the most conservative case, then use the total multiplicative systematic uncertainty as an input parameter to the POLE program.

Since there are few events observed from data sample at $\sqrt{s}$ = 10.52 GeV, the continuum contributions are neglected for the $\Upsilon(1S,2S)$ decays. The conservative upper limit on the product branching fractions in $\Upsilon(1S,2S)$ decays ${\cal B}^{\textrm{UP}}(\Upsilon(1S,2S)\to X_{cc\bar{s}\bar{s}}+anything)\times{\cal B}(X_{cc\bar{s}\bar{s}}\to D_{s}^{+}D_{s}^{+}(D_{s}^{*+}D_{s}^{*+}))$ are obtained by the following formula:

where $N^{\textrm{UP}}$ is the 90% C.L. upper limit on the number of events from the data signal yields including all systematic uncertainties that are mentioned above from other variables in this expression, $N_{\Upsilon(1S,2S)}$ is the total number of $\Upsilon(1S,2S)$ events, $\varepsilon_{i}$ is the corresponding detection efficiency, and ${\cal B}_{i}$ is the product of all secondary branching fractions for each reconstructed channel.

The conservative upper limit on the product values of Born cross section and branching fraction $\sigma^{\textrm{UP}}(e^{+}e^{-}\to X_{cc\bar{s}\bar{s}}+anything)\times{\cal B}(X_{cc\bar{s}\bar{s}}\to D_{s}^{+}D_{s}^{+}(D_{s}^{*+}D_{s}^{*+}))$ are calculated by the following formula:

where $N^{\textrm{UP}}$ is the 90% C.L. upper limit on the number of events in data signal yields including all systematic uncertainties that are mentioned above from other variables in this expression, $|1-\Pi|^{2}$ is the vacuum polarization factor, ${\cal L}$ is the integrated luminosity, $\varepsilon_{i}$ is the corresponding detection efficiency, ${\cal B}_{\textrm{i}}$ is the product of all secondary branching fractions for each reconstructed channel, and $(1+\delta)_{\textrm{ISR}}$ is the radiative correction factor. The values of $|1-\Pi|^{2}$ are 0.931, 0.930, and 0.929 for $\sqrt{s}$ = 10.52 GeV, 10.58 GeV, and 10.867 GeV vacuum , and the uncertainty is calculated to be less than 0.1%, which is negligible. The radiative correction factors $(1+\delta)_{\textrm{ISR}}$ are 0.686, 0.694, and 0.738, as calculated using the formula given in Ref. ISR for $\sqrt{s}$ = 10.52 GeV, 10.58 GeV, and 10.867 GeV, respectively, where we assume that the dependence of cross sections on $s$ is $1/s$.

The calculated 90% C.L. upper limits on the product branching fractions of $\Upsilon(1S,2S)\to X_{cc\bar{s}\bar{s}}+anything$ and the product values of Born cross section and branching fraction of $e^{+}e^{-}\to X_{cc\bar{s}\bar{s}}+anything$ at $\sqrt{s}$ = 10.52 GeV, 10.58 GeV, and 10.867 GeV for the mode $X_{cc\bar{s}\bar{s}}\to D_{s}^{+}D_{s}^{+}$ ($X_{cc\bar{s}\bar{s}}\to D_{s}^{*+}D_{s}^{*+}$) are displayed in Fig. 6 (7). Numerical values for the mode $X_{cc\bar{s}\bar{s}}\to D_{s}^{+}D_{s}^{+}$ can be found in Tables 3 and 4, while those for the mode $X_{cc\bar{s}\bar{s}}\to D_{s}^{*+}D_{s}^{*+}$ are shown in Tables 5 and 6.

Using the data samples of 102 million $\Upsilon(1S)$ events, 158 million $\Upsilon(2S)$ events, and data samples at $\sqrt{s}$ = 10.52 GeV, 10.58 GeV, and 10.867 GeV corresponding to integrated luminosities 89.5 fb^-1, 711.0 fb^-1, and 121.4 fb^-1, respectively, we search for the double-heavy tetraquark states $X_{cc\bar{s}\bar{s}}$ in the processes of $\Upsilon(1S,2S)\to D_{s}^{+}D_{s}^{+}(D_{s}^{*+}D_{s}^{*+})+anything$ and $e^{+}e^{-}\to D_{s}^{+}D_{s}^{+}(D_{s}^{*+}D_{s}^{*+})+anything$ at $\sqrt{s}$ = 10.52 GeV, 10.58 GeV, and 10.867 GeV. No peaking structures are observed in the $M_{D_{s}^{+}D_{s}^{+}}$ and $M_{D_{s}^{*+}D_{s}^{*+}}$ distributions from data. The 90% C.L. upper limits on the product branching fractions in $\Upsilon(1S,2S)$ inclusive decays [${\cal B}(\Upsilon(1S,2S)\to X_{cc\bar{s}\bar{s}}+anything)\times{\cal B}(X_{cc\bar{s}\bar{s}}\to D_{s}^{+}D_{s}^{+}(D_{s}^{*+}D_{s}^{*+}))$] and the product values of Born cross section and branching fraction for $e^{+}e^{-}\to X_{cc\bar{s}\bar{s}}+anything$ [$\sigma(e^{+}e^{-}\to X_{cc\bar{s}\bar{s}}+anything)\times{\cal B}(X_{cc\bar{s}\bar{s}}\to D_{s}^{+}D_{s}^{+}(D_{s}^{*+}D_{s}^{*+}))$] at $\sqrt{s}$ = 10.52 GeV, 10.58 GeV, and 10.867 GeV as functions of various assumed $X_{cc\bar{s}\bar{s}}$ masses and widths are determined.