Study of the lineshape of $X(3872)$ using $B$ decays to $D^{0}\overline{D}{}^{*0}K$

The charmonium-like $X(3872)$ state, also known as $\chi_{c1}(3872)$ [1], was discovered by the Belle experiment as a narrow peak in the vicinity of the $D^{0}\overline{D}{}^{*0}$ threshold in the $J/\psi\pi^{+}\pi^{-}$ invariant mass distribution in exclusive $B^{+}\to J/\psi\pi^{+}\pi^{-}K^{+}$ decays [2]. Its existence has been confirmed by multiple experiments: D0 [3], BABAR [4], CDF [5], LHCb [6], and BESIII [7]. In addition to the $J/\psi\pi^{+}\pi^{-}$ decay, other decays such as $J/\psi\omega$ [8], $J/\psi\gamma$, $\psi(2S)\gamma$ [9], $D^{0}\overline{D}{}^{*0}$ [10, 11], $D^{0}\overline{D}{}^{0}\pi^{0}$ [12], and $\pi^{0}\chi_{c0}$ [13] have been observed. The $X(3872)$ quantum numbers $J^{PC}$ have been determined to be $1^{++}$ [14, 15]. Various interpretations such as a loosely bound state [16, 17, 18, 19], an admixture of a molecular state and a pure charmonium resonance [20], a tetraquark [21], and a cusp at the $D^{0}\overline{D}{}^{*0}$ threshold [22, 23, 24] have been proposed, and the structure of the state remains uncertain. Measurement of the lineshape in various decay modes can help to discriminate among different choices for the structure. In this paper, we examine two models for the lineshape in the decay to $D^{0}\overline{D}{}^{*0}$: a Breit-Wigner, and a Flatté-inspired parametrization.

The $X(3872)$ peak has already been analyzed with the Breit-Wigner lineshape commonly used for resonance states. Based on the analyses of the decays including $J/\psi$, the mass is $3871.65\pm 0.06~{\rm MeV}/c^{2}$ and the width is $1.19\pm 0.21~{\rm MeV}$ [1], with two measurements at the LHCb experiment [25, 26] contributing significantly to these averages. Analyses of the decay to $D^{0}\overline{D}{}^{*0}$ based on the Breit-Wigner lineshape tend to yield a higher mass and a larger width, with the width measurement subject to large uncertainties [10, 11]. Discrepancies in the lineshape between the decays to the $J/\psi\pi^{+}\pi^{-}$ and $D^{0}\overline{D}{}^{*0}$ final states can arise near the threshold due to coupled-channel effects [22]. This may be significant for the $X(3872)$, as the observed mass coincides with the $D^{0}\overline{D}{}^{*0}$ threshold of $3871.69\pm 0.10~{\rm MeV}/c^{2}$, and a $1^{++}$ state can couple to the $D^{0}\overline{D}{}^{*0}$ channel in S-wave. One model to account for coupled-channel effects is the Flatté-inspired parametrization [22, 27], a Breit-Wigner model with an explicit expression for the energy-dependent partial width. At LHCb, an analysis of the $J/\psi\pi^{+}\pi^{-}$ invariant mass distribution was performed using this Flatté-inspired model [26]. It is difficult to determine all of the parameters using only this distribution, due to a scaling behavior in which the lineshape near the threshold does not change under a linear transformation of four of the five parameters [26, 28]. To determine all the parameters, it is important to analyze not only the $J/\psi\pi^{+}\pi^{-}$ decay but also the $D^{0}\overline{D}{}^{*0}$ decay, as proposed in the theoretical analysis [27]. By analyzing the $D^{0}\overline{D}{}^{*0}$ decay, we aim to provide more information on the lineshape, and in particular on the coupling strength of $X(3872)\to D^{0}\bar{D}^{*0}$.

In this paper, we present a study of the $X(3872)$ lineshape using a sample of $X(3872)\to D^{0}\overline{D}{}^{*0}$ candidates produced in the exclusive decay $B\to D^{0}\overline{D}{}^{*0}K$ using the full Belle dataset. There have been three previous studies [12, 10, 11]. Reference [12] is an analysis of the $B\to D^{0}\overline{D}{}^{0}\pi^{0}K$ decay at Belle, and Refs [10, 11] are analyses of the $B\to D^{0}\overline{D}{}^{*0}K$ decays at BABAR and Belle, respectively. The latter two analyses apply a $D^{*0}$ selection and a mass-constrained fit to the $D^{*0}$ candidates. While this has the advantage of improving the signal-to-noise ratio, it has the disadvantage of disallowing entries below the $D^{0}\overline{D}{}^{*0}$ threshold, which is important for studying the structure. Given the limited size of our data sample, we adopt a similar technique to the latter analyses, i.e. subtracting the reconstructed $D^{*0}$ mass and adding the nominal mass. The disadvantage of requiring the $D^{*0}$ is partially compensated for by analyzing the Flatté model, in which we can obtain a lineshape reflecting poles of the scattering amplitude. Compared to Refs [10, 11], additional $D^{0}$ decay modes are included, increasing the efficiency to reconstruct $D^{0}\overline{D}{}^{*0}$ decays. Throughout this paper, charge conjugation is always included. We do not distinguish $D^{0}\overline{D}{}^{*0}$ from $\overline{D}{}^{0}D^{*0}$ unless otherwise indicated.

The paper is organized as follows. In Sec. II, the Belle detector and data set are described. In Sections III and IV, the event selection and the fitted model are presented. In Sec. V, the results of fitting the data with the relativistic Breit-Wigner model and the Flatté model are presented. Section VI contains a discussion of the results, and the conclusions of the paper.

We use a data sample of $772\times 10^{6}$ $B\overline{B}$ pairs, collected at a center-of-mass energy of $\sqrt{s}=10.58~{\rm GeV}$, corresponding to the $\Upsilon(4S)$ resonance, with the Belle detector at the KEKB asymmetric-energy $e^{+}e^{-}$ collider [29, 30]. The Belle detector is a large-solid-angle magnetic spectrometer that consists of a silicon vertex detector (SVD), 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 comprised of CsI(Tl) crystals (ECL) located inside a super-conducting solenoid coil that provides a 1.5 T magnetic field. An iron flux-return located outside of the coil is installed to detect $K_{L}^{0}$ mesons and to identify muons (KLM). The detector is described in detail elsewhere [31, 32].

To determine the event selection and the detector response, we use a sample of Monte Carlo (MC) simulated events generated using the EvtGen event generator [33]. The detector response is simulated using the GEANT3 package [34].

The event selection is determined using the MC samples in two steps. First, the selection criteria for the final-state particles are determined based on our previous studies [12, 11]. Second, the selection criteria for the intermediate-state particles are optimized by maximizing the figure-of-merit $S/\sqrt{S+B}$, where $S$ and $B$ are the estimated numbers of signal and background events, respectively. The resulting selection is described below.

Tracks are selected using vertex information measured by the tracking system. A track candidate is accepted if its distance along the detector axis from the point of closest approach to the interaction point is less than 4.0 cm, and its distance transverse to the detector axis is less than 1.0 cm. These requirements are not imposed for tracks in $K_{S}^{0}\to\pi^{+}\pi^{-}$ candidates. In addition, pion and kaon candidates are selected using likelihoods ${\cal L}_{\pi}$ and ${\cal L}_{K}$ based on the time-of-flight measured by the TOF, the number of Cherenkov photons detected by the ACC, and the ionization loss in the CDC. Tracks with a likelihood ratio ${\cal L}_{\pi}/({\cal L}_{\pi}+{\cal L}_{K})>0.1$ are used as charged pion candidates, and tracks with ${\cal L}_{\pi}/({\cal L}_{\pi}+{\cal L}_{K})<0.9$ are used as charged kaon candidates. The hadron identification efficiency is approximately 97% for both pions and kaons. Tracks satisfying ${\cal L}_{e}/({\cal L}_{e}+{\cal L}_{\tilde{e}})>0.95$ are identified as electrons and eliminated. Here, ${\cal L}_{e}$ and ${\cal L}_{\tilde{e}}$ are distinct likelihoods for the electron and non-electron hypotheses, based on ECL, tracking, and other information. The particle identification is described in detail elsewhere [35].

$K_{S}^{0}$ candidates are reconstructed from charged pion pairs with opposite charges. The $\pi^{+}\pi^{-}$ invariant mass is required to agree with the known $K_{S}^{0}$ mass [1] within $7~{\rm MeV}/c^{2}$ ($\approx 3.6\sigma$ of the resolution). Candidates are selected using a neural network classifier [36] with various kinematic variables as input. To improve the four-momentum resolution, a mass- and vertex-constrained fit is applied.

Photon candidates are reconstructed from ECL clusters with no matching charged tracks. Candidates are selected based on the ratio, $E_{9}/E_{25}$, of the energy deposited in the $3\times 3$ array of crystals centered on the crystal with the highest energy deposition to that in the $5\times 5$ array: we require $E_{9}/E_{25}>0.8$.

Neutral pions are reconstructed from photon pairs. The photons are required to have energy greater than $30~{\rm MeV}$ in the barrel region or $50~{\rm MeV}$ in the endcaps. The $\gamma\gamma$ invariant mass is required to agree with the $\pi^{0}$ nominal mass [1] within $12~{\rm MeV}/c^{2}$. This mass window corresponds to 92% signal efficiency. A mass-constrained fit is applied to improve the momentum resolution.

$D^{0}$ candidates are reconstructed in six decay modes: $K^{-}\pi^{+}$, $K^{-}\pi^{+}\pi^{0}$, $K^{-}\pi^{+}\pi^{-}\pi^{+}$, $K_{S}^{0}\pi^{+}\pi^{-}$, $K_{S}^{0}\pi^{+}\pi^{-}\pi^{0}$, and $K^{+}K^{-}$. The $\pi^{0}$ candidates used in this reconstruction are required to have momentum in the center-of-mass system greater than $100~{\rm MeV}/c$, and energy in the laboratory system greater than $150~{\rm MeV}$. If a $\pi^{0}$ is included, the reconstructed $D^{0}$ invariant mass is required to be within $16~{\rm MeV}/c^{2}$ of the nominal mass [1] corresponding to 85% signal efficiency; otherwise, it is required to be within $8.5~{\rm MeV}/c^{2}$ corresponding to 91% efficiency. To improve the momentum resolution, a mass- and vertex-constrained fit is applied. Candidates where the $\chi^{2}$ probability of the fit is less than 0.0001 are eliminated.

$\overline{D}{}^{*0}$ candidates are reconstructed in two decay modes: $\overline{D}{}^{0}\gamma$ and $\overline{D}{}^{0}\pi^{0}$. For the $\overline{D}{}^{0}\gamma$ mode, only $\gamma$ candidates with an energy greater than $90~{\rm MeV}$ in the laboratory system are used. For the $\overline{D}{}^{0}\pi^{0}$ mode, only $\pi^{0}$ candidates with a momentum in the center-of-mass system of less than $100~{\rm MeV}/c$ and an energy in the laboratory system of less than $200~{\rm MeV}$ are used. The difference in the reconstructed mass between $\overline{D}{}^{*0}$ and $\overline{D}{}^{0}$ is required to agree with the nominal value [1] within $9.0~{\rm MeV}/c^{2}$ and $2.0~{\rm MeV}/c^{2}$ for $\overline{D}{}^{0}\gamma$ and $\overline{D}{}^{0}\pi^{0}$, respectively, corresponding to 90% signal efficiency in each case.

$B$ meson candidates are then reconstructed in the decay modes $D^{0}\overline{D}{}^{*0}K^{+}$ and $D^{0}\overline{D}{}^{*0}K_{S}^{0}$. To reduce wrong combinations, the daughter $K^{+}$ is required to have ${\cal L}_{K}/({\cal L}_{\pi}+{\cal L}_{K})>0.6$, corresponding to an identification efficiency of 89%. The $B$ candidates are selected based on the beam-energy constrained mass, $M_{\rm bc}\equiv\sqrt{(E_{\rm beam}^{\rm cms})^{2}-(p_{B}^{\rm cms})^{2}}$ and the difference of the energy in the center-of-mass system between the $B$ candidate and the beam, $\Delta E\equiv E_{B}^{\rm cms}-E_{\rm beam}^{\rm cms}$, where $E_{\rm beam}^{\rm cms}$ is the beam energy in the center-of-mass system corresponding to half of $\sqrt{s}$, and $p_{B}^{\rm cms}$ and $E_{B}^{\rm cms}$ are the energy and momentum of $B$ candidates in the center-of-mass system, respectively. We retain events with $M_{\rm bc}>5.2~{\rm GeV}/c^{2}$ and $|\Delta E|<50~{\rm MeV}$ for later analysis. The $M_{\rm bc}$ signal region is defined as $|M_{\rm bc}-m_{B}|<4.5~{\rm MeV}/c^{2}$ ($\approx 2\sigma$) for $\overline{D}{}^{*0}\to\overline{D}{}^{0}\gamma$ and $|M_{\rm bc}-m_{B}|<6.0~{\rm MeV}/c^{2}$ ($\approx 2.5\sigma$) for $\overline{D}{}^{*0}\to\overline{D}{}^{0}\pi^{0}$, where $m_{B}$ denotes the nominal $B$ mass [1]. The $\Delta E$ signal region is defined as $|\Delta E|<12~{\rm MeV}$ ($\approx 2\sigma$). For suppression of continuum events, we use a FastBDT classifier [37] trained on the simulation sample with the following event-shape information as input: modified Fox-Wolfram moments [38], the momentum flow in concentric cones around the thrust axis [39], and thrust-related quantities. Events for which the classifier output is less than 0.15 are eliminated. This requirement retains 96% of the signal candidates and rejects 49% of the candidates of continuum events.

After this selection, the average number of $B$ candidates per event is 1.8, because $D^{0}\overline{D}{}^{*0}$ and $D^{*0}\overline{D}{}^{0}$ are often indistinguishable and double-counted. To avoid multiple counting of signal events, we select the candidate that has the highest value of the product of the following likelihood ${\cal L}$ and prior probability ${\cal P}$

where ${\cal L}$ is the product of the likelihoods of the measured $D^{0}$, $\overline{D}{}^{0}$, and $\overline{D}{}^{*0}$ masses, and $\Delta E$; and, for the $\overline{D}{}^{*0}\to\overline{D}{}^{0}\pi^{0}$ mode, the likelihood of the measured $\pi^{0}$ mass. Each likelihood is obtained using probability density functions (PDFs) determined using the MC samples. The probability ${\cal P}$ is obtained from the probability that a signal event can be reconstructed $\varepsilon_{ijk}$, the average number of $B$ candidates per event $\zeta_{ijk}$, and the decay branching fraction, when $D^{0}$, $\overline{D}{}^{0}$, and $\overline{D}{}^{*0}$ are reconstructed in the $i$, $j$, and $k$ modes, respectively. The values of $\varepsilon_{ijk}$ and $\zeta_{ijk}$ are determined using the MC samples. The $(M_{\rm bc},\Delta E)$ distribution of the selected $B$ candidates is shown in Fig. 1. The red solid (blue dashed) rectangle shows the $(M_{\rm bc},\Delta E)$ signal region for $\overline{D}{}^{*0}\to\overline{D}{}^{0}\gamma$ ($\overline{D}{}^{0}\pi^{0}$): $B$ candidates used in the lineshape study are selected from this region.

For all events remaining in the selection, the following $D^{0}\overline{D}{}^{*0}$ invariant mass is calculated instead of applying a mass-constrained fit to improve the mass resolution,

where the reconstructed $\overline{D}{}^{*0}$ invariant mass, $M(\overline{D}{}^{0}\gamma)$ or $M(\overline{D}{}^{0}\pi^{0})$, is subtracted, and the $\overline{D}{}^{*0}$ nominal mass, $m_{\overline{D}{}^{*0}}$, is added. The lineshape and signal yield are determined by fitting the distribution in the region below $4.0~{\rm GeV}/c^{2}$.

In this work, the obtained $M(D^{0}\overline{D}{}^{*0})$ distributions are fitted with two lineshape models: the relativistic Breit-Wigner, and a Flatté-inspired model. The fits with these two models, shown in the next section, are performed with the following procedure.

When signal events are reconstructed correctly, the invariant mass distribution has a peak consisting of the natural lineshape convolved with the mass-dependent detector response. This response, i.e. the mass dependence of the signal efficiency and the mass resolution, is studied and parameterized using a set of $X(3872)\to D^{0}\overline{D}{}^{*0}$ MC samples generated with zero width, and a range of mass values from the $D^{0}\overline{D}{}^{*0}$ threshold to $4.0~{\rm GeV}/c^{2}$. Here, the $X(3872)\to D^{0}\overline{D}{}^{*0}$ decays are generated using a uniform phase space model; the $D^{*0}$ width is assumed to be around 60 keV [23]. Since the signal-to-noise ratio depends on the $D^{*0}$ decay mode, fits are performed separately for $D^{*0}\to D^{0}\gamma$ and $D^{*0}\to D^{0}\pi^{0}$. In addition, fits are performed separately for $B^{0}$ and $B^{+}$ candidates to determine the ratio of branching fractions between $B^{0}\to X(3872)K^{0}$ and $B^{+}\to X(3872)K^{+}$.

The fit function for correctly reconstructed $X(3872)\to D^{0}\overline{D}{}^{*0}$ decays, which we refer to as “signal”, is constructed as follows. The signal efficiency varies depending on the mass by a few tens of percent, especially around the threshold, as shown in Fig. 2 (a). It is parameterized by the threshold function $p_{0}\{1-p_{1}e^{p_{2}(M-m_{D^{0}}-m_{D^{*0}})}+p_{3}(M-m_{D^{0}}-m_{D^{*0}})\}$ with parameters $p_{0}$–$p_{3}$ in the low-mass region, which is continuously connected to a constant value in the high-mass region. The mass resolution for the signal is modeled as the sum of a Gaussian and a reversed Crystal Ball function [40] with a common mean. Figure 2 (b) shows the $M(D^{0}\overline{D}{}^{*0})$ spread due to the resolution, and the resolution function used in this work, for several choices of the $X(3872)$ mass. As noted in the previous Belle study, the resolution degrades with the square root of the difference between the mass and the threshold [11]. The convolution with the mass-dependent resolution function entails longer computation times. The effect of smearing due to the resolution is small at masses away from the peak, since the natural lineshape is broad [10, 11]. For example, the full width at half maximum (FWHM) of the natural lineshape is a few MeV, while the FWHM of the mass resolution near the peak is only about 220 keV. Therefore, instead of convolution with the mass-dependent resolution function, convolution with the specific resolution function at the mass of $3871.9~{\rm MeV}/c^{2}$ (near the peak) is adopted as an approximation. To reproduce the behavior near threshold, the signal function is multiplied by a soft threshold function that rises from zero to one at the threshold using an error function. The procedure is validated on $X(3872)\to D^{0}\overline{D}{}^{*0}$ MC samples generated with a broad lineshape. The effect of the approximation is negligible.

The ratios of signal yields among the decay modes are fixed in the fit using the product of the expected total signal efficiency and the branching fraction of each decay mode. Here the expected signal efficiency depends on the lineshape because of the mass dependence of the signal efficiency. The total signal efficiency is obtained by averaging the signal efficiencies as a function of mass weighted by the values of the lineshape function. It is then corrected by taking account of the signal which may leak out of the fit range, depending on the lineshape, and by taking the ratio of the area of the signal function in the fit range to that from the threshold to $m_{B}-m_{K}$. Here, a mass-dependent resolution function is convolved with the signal function, because smearing due to the resolution is important at higher masses. The calculation of the total signal efficiency is validated on MC samples for a broad range of lineshape parameters.

A separate fit function is used for “broken signal”: cases where the $D^{0}$ is reconstructed incorrectly, a wrong $\pi^{0}$ or $\gamma$ is combined in the $\overline{D}{}^{*0}$ reconstruction, or a ${D}^{*0}\overline{D}{}^{0}$ signal event is misinterpreted as $D^{0}\overline{D}{}^{*0}$ by combining $\pi^{0}$ or $\gamma$ from ${D}^{*0}$ incorrectly with the $\overline{D}{}^{0}$ to make a fake $\overline{D}{}^{*0}$. For the $\overline{D}{}^{*0}\to\overline{D}{}^{0}\pi^{0}$ mode, such events produce a broad peak in the $M(D^{0}\overline{D}{}^{*0})$ signal region and possibly distort the lineshape of the signal. The fit function for the broken signal therefore takes account of the mass dependencies of the resolution and the efficiency as in the case of correctly reconstructed signal events. The mass dependence of the efficiency is parameterized using the same threshold function used for the signal. The resolution is reproduced by a triple Gaussian multiplied by a soft threshold function at the $D^{0}\overline{D}{}^{*0}$ threshold, and its mass dependence is studied and parameterized using zero-width MC samples. Since the resolution for the broken signal is several times worse than that for the signal, we do not use approximated convolution: we instead use a discrete convolution with the mass-dependent resolution function. In the fit, the yield of the broken signal relative to the signal is fixed to the value expected from the lineshape and the ratio of the total efficiency of the broken signal to that of the signal.

The broken signal peak due to the $\overline{D}{}^{*0}\to\overline{D}{}^{0}\gamma$ mode has little sensitivity to the natural lineshape. To reduce the systematic uncertainty due to the shape, we use a histogram PDF depending on the lineshape. This PDF is obtained by plotting the broken-signal histogram for each of the zero-width MC samples, scaling it by the value of the assumed lineshape at the generated mass, and summing up all of the scaled histograms. Here, the bin widths are adjusted to increase as the mass increases to suppress statistical fluctuations.

The background from $e^{+}e^{-}\to q\bar{q}$ ($q=u,d,s,c$) continuum events, and $e^{+}e^{-}\to\Upsilon(4S)\to B\overline{B}$ events other than signal, is studied using the background MC sample. The shape of the invariant mass distribution is reproduced using a threshold function, $\sqrt{M-(m_{D^{0}}+m_{D^{*0}})}$, where $m_{D^{0}}$ and $m_{D^{*0}}$ are the nominal masses of $D^{0}$ and $D^{*0}$ [1], respectively.