Observation of a resonant structure near the $D_s^+ D_s^-$ threshold in the $B^+\to D_s^+ D_s^- K^+$ decay

An amplitude analysis of the $B^+\to D_s^+ D_s^- K^+$ decay is carried out to study for the first time its intermediate resonant contributions, using proton-proton collision data collected with the LHCb detector at centre-of-mass energies of 7, 8 and 13 TeV. A near-threshold peaking structure, referred to as $X(3960)$, is observed in the $D_s^+ D_s^-$ invariant-mass spectrum with significance greater than 12 standard deviations. The mass, width and the quantum numbers of the structure are measured to be $3956\pm5\pm10$ MeV, $43\pm13\pm8$ MeV and $J^{PC}=0^{++}$, respectively, where the first uncertainties are statistical and the second systematic. The properties of the new structure are consistent with recent theoretical predictions for a state composed of $c\bar{c}s\bar{s}$ quarks. Evidence for an additional structure is found around 4140 MeV in the $D_s^+ D_s^-$ invariant mass, which might be caused either by a new resonance with the $0^{++}$ assignment or by a $J/\psi \phi\leftrightarrow D_s^+ D_s^-$ coupled-channel effect.

invariant-mass fit can be found in the accompanying paper [45].
To improve the resolution on the masses of the two-body combinations that are used in the amplitude analysis, the four momentum of each final-state particle is determined from a kinematic fit [48] where the B + mass is constrained to its known value [32]. Figure 1 shows the resulting Dalitz-plot distribution for the B + → D + s D − s K + signal decays, where the non-B + background is subtracted by the sPlot technique [49] with the reconstructed B + mass as the discriminating variable. The most evident feature is the band near the D + s D − s threshold. To validate that this peaking structure is not due to the combinatorial background, the D + s D − s invariant-mass distribution of candidates in the B + mass region from 5360 to 5600 MeV is investigated and no peak is observed.
Employing an unbinned maximum-likelihood method, an amplitude fit with the sFit technique [50] is performed to investigate the intermediate states and determine the quantum numbers J P C of any new particle. Two known 1 −− charmonium states, ψ(4260) and ψ(4660) [32,46,47], and two new 0 ++ X states are needed to fit the structures in the D + s D − s spectrum. One of these scalars, X(3960), describes the D + s D − s threshold enhancement and the other, designated X 0 (4140), is necessary to model the dip around 4140 MeV, as shown in Fig. 2. The subscript 0 is used to distinguish the latter from the 1 ++ X(4140) state seen in the J/ψϕ final state [32]. Additionally, an S-wave three-body phase-space function [32] is employed to model the nonresonant (NR) B + → D + s D − s K + component. Since no significant contribution of any state is observed in either the D − s K + or D + s K + systems, these five contributions constitute the baseline model. The helicity formalism [51] is used to construct the amplitude model of the where M 0 is the mass of the resonance, g j denotes the coupling of this resonance to the j-th channel, ρ j (m) is the phase-space factor [32] for the j-th two-body decay.  value of m is below the threshold of the channel j, i.e. q 2 j < 0, an analytic continuation is applied for q j = i −q 2 j [55, 56]. The total width of the resonance is calculated as Γ 0 = j g j ρ j (M 0 ). In the baseline model, only the D + s D − s channel (j = 1) is included in the Flatté-like parameterisation.
The total probability density function is the squared modulus of the total decay amplitude multiplied by the efficiency, normalised to ensure that the integral over the Dalitz plot is unity. The fit fraction F i expresses the fraction of the total rate due to the component i, and the interference fraction I ij describes the interference between components i and j. They are defined in Eqs. (18) and (19) of Ref.
[53], such that i F i + i<j I ij = 1. As shown in Fig. 2, the two-body mass distributions are well modelled by the baseline amplitude fit. The corresponding numerical results are summarised in Table 1, including the mass, width, fit fraction, and significance (S) of each component. The significance of a given component is evaluated by assuming that the change of twice the negative log-likelihood (−2 ln L) between the baseline fit and the fit without that component obeys a χ 2 distribution, where the number of degrees of freedom (n.d.f.) is given by the Table 1: Summary of the main results obtained with the baseline model, where the first uncertainty is statistical and the second systematic. The last column shows the signal significance with (without) the systematic uncertainty included.
The J P C assignment for the system of a pair of oppositely-charged pseudoscalar mesons must be in the series 0 ++ , 1 −− , 2 ++ , etc. States with higher intrinsic spin are not expected to contribute significantly in the current dataset. To determine the X(3960) quantum numbers, fits with the baseline model are performed under alternative J P C hypotheses, 1 −− , 2 ++ , instead of 0 ++ . The significance to reject a J P C hypothesis is computed as ∆(−2 ln L), where ∆(−2 ln L) = −(2 ln L(0 ++ ) − 2 ln L(J P C )), and indicates the likelihood difference between the fits for the preferred 0 ++ assignment and for each alternative J P C hypothesis. To ensure that for different J P C hypotheses this resonance corresponds to the same particle, the mass and the width are limited to be within a ±3 σ range of the baseline fit results. The 0 ++ assignment is preferred over 1 −− and 2 ++ hypotheses by 9.3 σ and 12.3 σ, respectively. Similarly, replacing the baseline 0 ++ assignment by 1 −− or 2 ++ for the X 0 (4140) state deteriorates the fit quality. The 0 ++ assignment is favoured over 1 −− (2 ++ ) hypothesis at a 3.5 σ (4.2 σ) level. Within the baseline model this 0 ++ state produces the dip around 4140 MeV via destructive interference with the 0 ++ NR and X(3960) components, with the interference fractions of, respectively, (−22.4 ± 6.4)% and (−5.2 ± 3.9)%, where the uncertainties are statistical only.
Systematic uncertainties on the measured resonance properties are evaluated, and are summarised in Table S1 in the supplemental material [61]. Corrections, derived from calibration samples, are applied to account for possible discrepancies between data and simulation in the hardware trigger and particle-identification responses. The uncertainty due to the limited size of the simulation samples is evaluated using the bootstrap method [62]. Additional resonances, not included in the baseline model (states in the D + s D − s system: 0 ++ χ c0 (4500) and χ c0 (4700) [12], 1 −− ψ(4040), ψ(4160) and ψ(4415) [32], and 2 ++ χ c2 (3930) [14]; and in the D − s K + system: 0 [32,63] and D * 1 (2760) 0 [64], and 2 + D * 2 (2460) 0 [32]) are utilised to estimate the uncertainty due to insufficient consideration of possible amplitude components. None of these states significantly improve the baseline model. The ccss candidates χ c0 (4500) and χ c0 (4700) have statistical significances of 0.8 σ and 1.3 σ, respectively, and their fit fractions are (0.6 ± 1.0)% (< 3.5% at 90% confidence level) and (2.4 ± 1.8)% (< 6.7% at 90% confidence level), where the uncertainties are statistical. The Blatt-Weisskopf hadron size is varied between 1.5 and 4.5 GeV −1 . The fixed masses and widths of two baseline ψ states are varied by their corresponding uncertainties. The Flatté-like parameterisation for the X(3960) state is replaced by a constant-width relativistic Breit-Wigner function. The uncertainty due to the possible bias of the sFit method is evaluated using pseudoexperiments. The total systematic uncertainties on mass, width, and fit fraction are obtained by adding all contributions in quadrature, assuming that each source is independent. Regarding the total significance for each component in the baseline model, the smallest significance among these systematic tests is selected.
The measured mass and width of the X(3960) state are consistent with those of the χ c0 (3930) meson [14] within 3 σ. Assuming that the X(3960) in the D + s D − s system and the χ c0 (3930) in the D + D − system are the same state, the baseline model is extended by adding a second channel (D + D − ) in the Flatté-like parameterisation. The corresponding fit projections and numerical results can be found in the supplemental material [61]. The likelihood is essentially unchanged while the n.d.f. is increased by one compared to the baseline fit. The coupling strength of the X(3960) state to D + s D − s (D + D − ) is found to be 0.33 ± 1.18 (0.15 ± 0.33) GeV. The masses and fit fractions of all components are consistent with those in the baseline one-channel Flatté-like model.
In the case that the X(3960) and χ c0 (3930) states are the same particle, the partial width ratio of such an X resonance decaying to D + s D − s and D + D − final states is calculated as where the superscripts (1) and (2) X is the fit fraction of the X(3960) resonance presented in this Letter, and the branching fraction ratio B (1) /B (2) is taken from the accompanying paper [45]. The first uncertainty is statistical, the second systematic, and the third is due to uncertainties in the measured branching fractions, B(D + →K − π + π + ) and B(D + s →K − K + π + ) [32], and the uncertainty on F (1) . This ratio is compatible with that of the couplings mentioned above.
It is well known that the creation of an ss quark pair from the vacuum is suppressed relative to uu or dd pairs. Moreover, the X → D + s D − s decay, occurring near the threshold, has a rather smaller phase-space factor than that of X → D + D − . These two features indicate that Γ(X → D + D − ) should be considerably larger than Γ(X → D + s D − s ) if X does not have any intrinsic ss content. However, the value measured in Eq. (2) contradicts this expectation. This implies that the X(3960) and χ c0 (3930) are either not the same resonance, or they are the same non-conventional charmonium-like state, for instance, a candidate containing the dominant ccss constituents predicted in recent theoretical models [38][39][40][41][42][43]65]. Further studies are needed to gain insights into the nature of the D + s D − s threshold enhancement, in particular the measurement of the relative branching fraction for the D (s) D (s) and ωJ/ψ channels, produced in a different environment, such as from two-photon fusion processes by the Belle II experiment.
To test the possibility that the dip in the D + s D − s invariant mass around 4140 MeV can be produced by the opening of the nearby J/ψϕ threshold, without introducing an additional resonance, we employ a simple K-matrix model that contains the single resonance X(3960) and two coupled channels, D + s D − s and J/ψϕ. The K-matrix reads where K 12 = K 21 , and the subscripts 1 and 2 represent D + s D − s and J/ψϕ final states, respectively. One possible choice for the 2 × 2 K-matrix parameterisation [32] is where M R refers to the bare mass of the resonance R, m is the D + s D − s invariant mass, g R a denotes the bare coupling of the resonance R to the channel a, and the f ba is a real matrix parameterising the non-pole part of the K-matrix. As the X(3960) mass is about 160 MeV lower than the J/ψϕ threshold and its width is less than 50 MeV, the coupling of the X(3960) state to J/ψϕ should be negligible, giving g R 2 = 0. This results in the X(3960) resonance entering the K 11 element only. The production amplitude is expressed in the P -vector formalism [32,68,69], which gives where β R and β b are complex free parameters due to rescattering effects or missing channels [60]. The amplitude M is where ρ = diag{ρ 11 , ρ 22 } is the diagonal matrix composed of phase-space factors, I represents the identity matrix, and a = 1 for the D + s D − s channel under consideration. The fit demonstrates that the dip around the J/ψϕ threshold can also be modelled by the J/ψϕ → D + s D − s rescattering, and results in a −2 ln L that is worse by 6.0, while the n.d.f. is increased by one, compared to the baseline fit. The fit projections and numerical results can be found in the supplemental material [61]. Since the fit quality of the K-matrix parameterisation is close to that of the baseline model, a strong conclusion cannot be drawn whether the dip is due to destructive interference with the X 0 (4140) resonance or caused by the J/ψϕ → D + s D − s rescattering.
In addition, it is found that the fits with the two-channel Flatté-like and K-matrix parameterisations are unstable, due to having too many free parameters for the limited data sample size. Consequently, the statistical uncertainties for some parameters are large. Therefore, neither of these parameterisations are taken as the baseline model.
In conclusion, the first amplitude analysis of the B + → D + s D − s K + decay is performed using pp collision data with an integrated luminosity of 9 fb −1 collected by the LHCb experiment between 2011 and 2018. A peaking structure near the D + s D − s mass threshold, denoted as X(3960), is observed with a significance larger than 12 σ. Its quantum numbers are determined to be J P C = 0 ++ , favoured over 1 −− or 2 ++ with a significance greater than 9 σ. As argued above, the X(3960) and χ c0 (3930) states are unlikely to be the same pure conventional charmonium state. The X(3960) resonance presented in this Letter could be a candidate for an exotic state predominantly consisting of ccss constituents, as suggested in recent theoretical literature [38][39][40][41][42][43]65]. If predominant ccss content is confirmed, this state should be labelled T f ψϕ (3960) in the new naming scheme for exotic hadrons [70]. In addition, a dip around 4140 MeV can be described either by a 0 ++ X 0 (4140) resonance having a significance of 3.5 σ, or the coupled-channel effect of the J/ψϕ ↔ D + s D − s reaction. The data from the forthcoming Run 3 of the LHCb experiment and from the Belle II experiment will be critical to clarify the nature of these phenomena. Table S1: Summary of the systematic uncertainties on the parameters measured in the amplitude analysis. Mass (M 0 ) and width (Γ 0 ) are in units of MeV, while fit fraction (F) is in %. Contribution J P C M 0 (MeV) Γ 0 (MeV) F (%) X(3960) 0 ++ 3951 ± 14 38 ± 104 25.0 ± 7.6 X 0 (4140) 0 ++ 4133 ± 7 67 ± 16 16.7 ± 4.6 ψ(4260) 1 −− 4230 [60] 55 [60] 3.6 ± 0.4 ψ(4660) 1 −− 4633 [32] 64 [32] 2.2 ± 0.2 NR 0 ++ --45.9 ± 10.9 Table S3: Main results found from the K-matrix fit. Uncertainties are statistical only.         [50] Y. Xie, sFit: a method for background subtraction in maximum likelihood fit, arXiv:0905.0724.