$Z_{cs}(3985)$ in next-to-leading-order chiral effective field theory -- the first truncation uncertainty analysis

We revisit the $D_s^-D^{*0}$/$D_s^{*-}D^0$ interaction and the $Z_{cs}(3985)$ state in chiral effective field theory(EFT) up to the next-to-leading order. We examine the relative importance of the leading-order contact, one-eta exchange, next-to-leading-order contact, and two-kaon-exchange contributions. We show that the leading-order and next-to-leading-order contact contributions play the most important role such that the $Z_{cs}(3985)$ state qualifies as a $D_s^-D^{*0}$/$D_s^{*-}D^0$ resonance. On the other hand, the weak one-eta-exchange and weakly energy-dependent two-kaon-exchange contributions are less important in dynamically generating the $Z_{cs}(3985)$ state, indicating that chiral EFT is less predictive in the present situation. Furthermore, we apply the Bayesian method to estimate chiral truncation uncertainties and find that they are of similar magnitude as their statistical counterparts. Our study shows that if $Z_c(3900)$ exists, then SU(3)-flavor symmetry also predicts $Z_{cs}(3985)$ with a certain robustness, i.e., both can be accommodated in the chiral EFT with dominant contact contributions.


I. INTRODUCTION
In 2020 the BESIII Collaboration reported on the existence of a hidden-charm state with strangeness, Z cs (3985), with a significance of 5.3σ [1]. The obtained mass and width are: Γ Zcs = 12.8 +5.3 −4.4 ± 3.0 MeV.
In 2021 the LHCb Collaboration observed two more hiddencharm tetraquark states with strangeness in the J/ψK + invariant mass spectrum of the B + → J/ψφK + decay, Z cs (4000) + and Z cs (4200) + [2]. Because the width of Z cs (3985) and that of Z cs (4000) differ by one order of magnitude and their production mechanism is also different, they are quite unlikely the same state [2,3]. As a result, in the present work, we only focus on the Z − cs (3985) state. There were quite a number of studies on the likely existence of a ccsq state in the literature before the experimental discovery [4][5][6][7][8][9][10][11]. In Ref. [9], assuming a compact tetraquark picture and considering the chromomagnetic interaction, the authors found that the lowest ccsq tetraquark is about 100 MeV below the DD s mass threshold. In the relativistic quark model, the lowest ccsq state is found to be 10 MeV below the DD s mass threshold [4]. In Ref. [7] no D ( * )D ( * ) s bound state was found in the one-boson-exchange (OBE) model, where only π, σ, ω, and ρ exchanges were considered. In the chiral unitary approach constrained by the hidden gauge symmetry and broken SU(4) symmetry [5], no resonance or bound state near the DD s mass threshold was found. In Ref. [8] an effective * Electronic address: ljxwohool@buaa.edu.cn † Electronic address: lisheng.geng@buaa.edu.cn field theory study with constraints from heavy quark spin symmetry and SU(3)-flavor symmetry predicted several hadronic molecules near the D ( * )D ( * ) s thresholds. 1 After the BESIII discovery [1], many new studies were performed and some of the earlier studies were updated. As the Z cs (3985) state lies close to the mass thresholds of D − s D * 0 and D * − s D 0 , the molecular picture has gained a lot of attention [13][14][15][16][17][18][19][20][21][22][23][24][25][26][27][28][29]. There exist also competing interpretations. For instance, it has been suggested to be either a compact tetraquark state [24,[30][31][32][33][34][35][36][37][38][39], a hadro-charmonium [38], a virtual state [40], a threshold effect [41], a reflection effect [42], or a recoupling effect [43]. A few remarks are in order regarding some of these studies. It was shown that the OBE potential does not support the existence of a D − s D * 0 /D * − D 0 molecule [15], while the QCD sum rule approach cannot distinguish between a molecule and a compact tetraquark state [18,19,24,32]. There is no candidate of either molecular nature or compact tetraquark nature for the Z cs (3985) state in the chiral quark model of Ref. [44].The existence of the Z cs (3985) state in the OBE model of Ref. [25] depends on a number of less known factors, such as the η exchange and the scalar meson exchange. In addition to its mass and spin-parity, the decay and production mechanisms [45][46][47][48], electromagnetic properties [49], and even medium modifications [31,50] of the Z cs (3985) state have been extensively studied.
In the molecular picture, it is particularly interesting to note that in Ref. [20], using the next-to-leading-order chiral poten- 1 These results should be taken with caution because there X(3915) and Y (4140) were assumed to be D * D * andD * s D * s hadronic molecules with quantum numbers J P C = 0 ++ . However, the J P C of Y (4140) has now been determined to be 1 ++ [12]. tial, with the two low-energy constants (LECs) and the cutoff determined from the Z c (3900) data [51] and SU(3)-flavor symmetry, it is shown that the K + recoil-mass spectrum of the e + e − → K + (D − s D * 0 + D * − s D 0 ) reaction can be well described. With the same three parameters, a pole is found at (m, Γ) = (3982.4 +4.8 −3.4 , 11.8 +5.5 −5.2 ) MeV, which can be associated with the Z cs (3985) state.
An effective field theory is a low-energy realization of the underlying (and more fundamental) theory, which satisfies all the relevant global symmetries and allows for modelindependent descriptions of physical processes in its validity domain. In the present work, the underlying theory of the strong interaction is quantum chromodynamics (QCD), while the effective field theory is chiral effective field theory (ChEFT). It has been successfully applied to study many low-energy phenomena, in particular, the nucleon-nucleon interaction [52][53][54][55], pioneered by Weinberg [56,57]. For a review of the application of ChEFT in the heavy quark sector, see Ref. [58]. A key issue in any EFT is that it allows for a controlled estimate of truncation uncertainties in a statistically meaningful way.
In the past, chiral truncation uncertainties at a given order were usually estimated by the variation of the cutoff needed to regularize the corresponding scattering equation, e.g. in Ref. [59]. However, as pointed out in Ref. [60], the resulting uncertainties are actually the residue effect of cutoff dependence, which essentially originate from the missing contact interactions that appear only at even chiral orders. Meanwhile, the uncertainties obtained in this way are somehow arbitrary since they depend on the employed cutoff range. In Ref. [60], the authors proposed to use the differences between the optimal results obtained at different orders as an estimate of truncation uncertainties, which are referred to as the EKM uncertainties [54]. They are motivated by the expectation that at low energies, truncation uncertainties are dominated by powers of expansion quantities such as momentum or light quark masses. For a short review of the application of the EKM approach in nuclear physics, see Ref. [61]. The EKM uncertainties are by construction smaller order by order (given that the EFT converges). However, since only the known contributions are taken into account, a statistical interpretation of this algorithm is missing. More recently, a new method based on the Bayesian model has been proposed [62][63][64], in which by encoding the expectation of the chiral expansion coefficients in a "prior probability density function", the truncation uncertainties for a certain degree-of-belief (DoB) can be determined by integrating out the coefficients of omitted orders. The new Bayesian approach and its modified version have been applied in the latest studies of pion-nucleon scattering [65], NN and 3N scattering [55,66,67] and other nuclear physics observables [68][69][70][71].
In this work, for the first time, we apply the Bayesian model to the heavy quark sector and to studies of exotic hadrons. In particular, we study the D − s D * 0 /D * − s D 0 scattering and the related Z cs (3985) state in chiral effective field theory (ChEFT) up to the next-to-leading order (NLO). In particular, we focus on the decomposition of the NLO chiral potential and apply the Bayesian model to study the chiral truncation un-certainties. We show that the expansion converges well up to the breakdown scale of the cutoff and even taking into account chiral truncation uncertainties, the state remains as a resonance, supporting its identification as the Z cs (3985) state [20].
This work is organized as follows. We briefly explain the ChEFT approach and the Bayesian model in Sec. II. Results and discussions are given in Sec.III, followed by a short summary in the last section.

II. THEORETICAL FORMALISM
In this section, we introduce the chiral potentials up to the next-to-leading order and explain the Bayesian model adopted in the present work.
A. Chiral potentials up to next-to-leading order In this work, we follow Refs. [20,51] and describe the e + e − → K + (D − s D * 0 − D * − s D 0 ) reaction in two steps as shown in Fig. 1: the decay of a virtual pho- The amplitude of the whole process U(E, p p p) can be obtained by solving the following Lippmann-Schwinger equation where E, p p p and µ are the energy, momentum and reduced mass of the D − s D * 0 /D * − s D 0 system in Fig. 1 respectively and |p p p| = 2µ(E − m th ) with m th the threshold of the system. M(E, p p p) is the photon decay amplitude, which can be described by the following effective Lagrangian where g γ denotes the effective coupling constant, F µν is the field strength tensor of the virtual photon, ) collect the charmed vector/pseudoscalar meson fields withP ( * ) the corresponding anti-meson fields. The u ν = i 2 {ξ † , ∂ µ ξ} is the axial-vector field in which ξ =exp( iφ 2f φ ) with φ the pseudoscalar octet and f φ the decay constant.
The potential V(p, q) consists of four parts: the LO and NLO contact potential, the one-eta-exchange (OEE) potential, and the two-kaon-exchange (TKE) potential. Their explicit expressions [20] are where V 1 is As was demonstrated in Ref. [51], the S-D mixing effect is insignificant and in the present work, we neglect the D-wave interaction as was done in Ref. [20]. In the above chiral potential, f η = 116 MeV, f K = 113 MeV, g = 0.57, p and p ′ represent the momenta of initial and final state in the centerof-mass (c.m) system, q q q = p p p ′ − p p p, m η and m K are the masses of eta and kaon, θ is the scattering angle in the c.
, is multiplied to the chiral potential to remove the ultraviolet divergence in solving the Lippmann-Schwinger equation.
It should be noted that the OEE and TKE potentials are pure predictions, while the two LECsC S and C S need to be determined by fitting to experimental data. If the OEE and TKE potentials are attractive and strong enough to allow for formations of bound states or resonances, then the ChEFT is predictive. Otherwise, it is not. It is one of our purposes to check whether this is the case for the D − s D * 0 /D * − s D 0 system.

B. Bayesian model
In this subsection, we briefly explain how the Bayesian model can be extended to study the likely existence of resonances. For details, please refer to Refs. [62][63][64].
Consider the EFT expansion for an observable X at fixed kinematics, where {c n } are dimensionless expansion coefficients. X ref is a natural-sized X taken as a reference which is suggested [66] to be Q is the expansion parameter, which is assumed to take the form [54,60,62,63] with p the momentum in the c.m. frame, and Λ the breakdown scale. However, in the present work, as we show later, the optimal cutoff obtained is less than 0.4 GeV. Therefore, we could not use the kaon mass or eta mass to estimate the expansion parameter. Instead we use an effective mass of 200 MeV to put a lower limit on the expansion parameter, which is about 0.5. ∆X (2) is the difference between the NLO result and the LO result. If the series is truncated at order k, then the truncation uncertainty is defined as the sum of the contributions of chiral orders higher than k removing the dimensional X ref .
To estimate ∆ k , one has to start from the known {c n }(n ≤ k). That is, one needs to specify the probability distribution function p(∆ k |c 0 , c 1 , . . . , c k ), with which the degree-of-belief (DoB) intervals can be calculated as with d (p) k the integral interval corresponding to the DoB p%. Usually one chooses p% = 68%. The truncation uncertainty for the observable X is then k , which means that the possibility of the value of observable X in (X (k) − ∆X (k) , X (k) + ∆X (k) ) is p%.
Following Ref. [63], we ultilize a Gaussian prior probability distribution function (pdf) with a common hyperparameter c for the expansion coefficients c i as wherec itself follows the log-uniform pdf [72] as wherec < (c > ) is the lower(upper) limit ofc and θ(x) here denotes the step function. The posterior pdf p(∆ k |c 0 , c 1 , . . . , c k ) then takes the following form [66] p(∆ k |c 0 , c 1 , . . . , c k ) = 1 with h the next highest order to be taken into account, c c c 2 k ≡ i∈A c 2 i , A ≡ {n ∈ N 0 |n ≤ k ∧ n = 1 ∧ n = m} with c m the expansion coefficient corresponding to the overall scale X ref defined above. In the present work, we take k = 2, h = 10,c < = 0.5 andc > = 10.0 following Ref. [63].

III. NUMERICAL RESULTS AND DISCUSSIONS
In this section, we study the decomposition of the chiral potential up to NLO and perform two types of truncation uncertainty analyses. In Table I, we tabulate the masses of the relevant particles.

A. Decomposition of the chiral potential
First, we examine the relative importance of the different terms of the chiral potential up to NLO. In the nucleonnucleon system, it is well known that the one-pion exchange(OPE) potential provides the longest-range nuclear force, while the two-pion exchange(TPE) contributions describe the intermediate part. As a result, for higher partial waves, e.g., those with angular momentum larger than 3 or 4, the OPE plus TPE potentials can already describe well the corresponding partial wave phaseshifts [73,74]. On the other hand, for those channels of low angular momenta, the short-range contributions encoded in the low-energy constants (LECs) are important, to such an extent that for low energy regions the pions can be integrated out [75]. The above discussion is relevant because if the OPE and TPE contributions are dominant, the theory is more predictive. Otherwise, one has to rely on either experimental data or lattice QCD data or phenomenology (e.g., resonance saturation) to determine the relevant LECs. This is particularly relevant to the D − s D * 0 /D * − s D 0 interaction and the corresponding Z cs (3985) state. This is because if the LO OEE potential or the NLO TKE potential is attractive enough, then without relying on other inputs, the ChEFT approach is predictive by itself. Otherwise, certain experimental inputs are needed to predict the existence of Z cs (3985).
In Fig. 2 with the two LECs and the cutoff determined in Ref. [51] by fitting to the Z c (3900) data, i.e.,C S = 3.6 × 10 2 GeV −2 , C S = −76.9 × 10 2 GeV −4 , and Λ = 0.33 GeV, we compare the LO contact, LO OEE, NLO contact, and NLO TKE contributions around the D − s D * 0 /D * − s D 0 threshold. It is clear that the OEE contribution is negligible compared to the LO contact term. In addition, in the energy region studied, the TKE contribution is nearly energy independent and is only about one quarter of the LO contact contribution. From this, we conclude that the OEE or TKE contribution alone is not enough to dynamically generate the Z cs (3985) state as a D − s D * 0 /D * − s D 0 resonance. On the other hand, the NLO contact contribution is zero at threshold but increases quickly as one moves away from the threshold. At the pole position, its size is about half that of the LO contact contribution. 50 MeV above the threshold, it already becomes four times stronger. As a matter of fact, the NLO contributions (both contact and TKE) are responsible for the appearance of a resonant D − s D * 0 − D * − s D 0 state, because with the energyindependent LO contribution, only bound or virtual states can emerge.
It is interesting to note that the LO contact potential is repulsive while the OEE potential is attractive but negligible in magnitude. The TKE potential is moderately attractive while the NLO contact potential is attractive above the threshold but repulsive below the threshold.

B. Truncation uncertainty analysis
We perform two types of truncation uncertainty analyses. One is based on the chiral potential, which itself is not a physical observable. The other is based on the K + recoil-mass spectrum, which is a more proper observable.

Uncertainty analysis based on chiral potential
In Fig. 3, we plot the NLO chiral potential of Ref. [20] as a function of the D − s D * 0 /D * − s D 0 invariant mass and the corresponding truncation uncertainty for a DoB of 68%. The chiral truncation uncertainty is relatively small close to the threshold, but becomes larger around the breakdown scale. We need to mention that for the potential below the threshold, we have replaced p with |p| in the calculation of the expansion parameter Q. One can translate the truncation uncertainties of the chiral potential shown in Fig. 3 to those of the two LECs,C S and C S . They now becomeC S = 3.60 +1.2+0.5 −1.2−0.5 GeV −2 and C S = −76.9 +6.2+10.0 −6.2−10.0 GeV −4 , where the first uncertainty is statistical and the second is systematic originating from chiral truncations.
With the NLO chiral potential and the corresponding truncation uncertainty as well as statistical uncertainty, we search for poles on the second Riemann sheet. The resulting positions are shown in Fig. 4. 3 Clearly, even with the truncation uncertainties taken into account, ChEFT still supports the interpretation of the Z cs (3985) state as a D − s D * 0 /D * − s D 0 resonance. More specifically, with both uncertainties, we obtain the position as (m, Γ) = (3982.4 +6.7 −5.2 , 11.8 +8.9 −7.4 ) MeV, which should be compared with that obtained in Ref. [20], (m, Γ) = (3982.4 +4.8 −3.4 , 11.8 +5.5 −5.2 ) MeV. Clearly even compared to the relatively large statistical uncertainty, the systematic uncertainty is sizable.

Uncertainty analysis based on K + -recoil mass spectrum
The leading order potential of Ref. [20] describes the K +recoil spectrum very badly (see Fig. 5), therefore to perform a proper truncation uncertainty analysis, we perform a new fit to the experimental data of Ref. [1]. At LO, we have two LECs to determine, i.e.,C S and Λ. At NLO, the LEC C S also contributes. By fitting to the experimental data shown in Fig. 5  we obtain the LECs given in Table II together with the corresponding χ 2 /d.o.f.. The light blue band is generated using the Bayesian method explained in Sec.IIB with Λ = 0.36. It should be mentioned that in the spirit of the Bayesian model, we have fixed the cutoff at its central value in obtaining the statistical uncertainties of the LECs. For the sake of comparison, we also show the LO and NLO predictions of Ref. [20]. It should be noted that the experimental data we fitted are taken from Fig. 4 of Ref. [1] where the wrong-sign combinations of D − s and K − candidates and the excited D * * + s contributions have been subtracted. Our fits indicate that there is no need to consider K rescattering from either D ( * )− s or D ( * )0 , in agreement with Refs. [13,17,20,28,29]. Once more precise data become available, one may need to examine this part of finalstate interactions in more detail. the following observation. Among the three data points neglected, two of them have negative central values and the one at threshold cannot be simultaneously described together with the others. A few things are noteworthy. First, the NLO description of the experimental data is quantitatively better than its LO counterpart. Nonetheless, both the LO and NLO fits can be viewed as reasonable given the relatively large experimental uncertainties. We note that in the literature, EFT studies have been performed at both LO and NLO, and reasonable descriptions of the experimental data have been claimed. Second, the truncation uncertainties increase quickly as one moves away from the threshold. This is mainly because of the relatively small optimal cutoff of 0.36 GeV 5 obtained at NLO, which implies that the NLO ChEFT can be trusted at most up to √ s D − s D * 0 /D * − s D 0 ≈ 4.041 GeV. With the parameters given in Table II, we search for poles on the complex plane. There is no bound state or virtual state at LO, but we find a resonant state at NLO. The position is (M, Γ) = (3983.5 +9.8 −6.5 , 21.3 +18.5 −12.5 ) MeV, where the uncertainties include both systematic and statistical ones. We note that the pole position is consistent with that of Ref. [20] within the relatively large uncertainty.

IV. SUMMARY AND OUTLOOK
We revisited the e + e − → K + (D − s D * 0 + D * − s D 0 ) reaction in chiral effective field theory up to the next-to-leading order in the single-channel approximation. We first examined the relative importance of various contributions in the chiral expansion and showed explicitly the important role played by the leading order and next-to-leading-order contact contributions in describing the BESIII data. This demonstrates that  [20]. The solid dots associated with errors are experimental data after removing the combinatorial backgrounds from Ref. [1] for the particular case studied, the ChEFT approach is not very predictive by itself. Experimental inputs or other symmetries are needed to predict the existence of the Z cs (3985) state. This was done in Ref. [20] using the Z c (3900) data and SU(3) symmetry. Then we applied the Bayesian method to estimate chiral truncation uncertainties. We performed two types of uncertainty analyses, based on either the chiral potential or the event distribution. For the latter, we have refitted the BESIII data. Our results showed that because of the relatively small cutoff obtained from the fitting, the chiral EFT approach can only be trusted up to 50 MeV above the D − s D * 0 /D * − s D 0 threshold. In either case, our results showed that even with chiral truncation uncertainties taken into account, the Z cs (3985) state remains a robust D − s D * 0 /D * − s D 0 resonance.
To the best of our knowledge, the present study constitutes the first exploration of the Bayesian model in studies of exotic hadrons. It is shown that a proper and statistically meaningful uncertainty analysis can be performed. Such an analysis gives us more confidence in the prediction and/or results of the chiral effective field theory studies. We hope that the present study can stimulate more similar studies in the heavy flavor sector in the future.  II: LECs of the LO and NLO chiral potential and the cutoff needed to regularize the chiral potential obtained by fitting to the BESIII data shown in Fig.5 (the solid points).CS is in units of 10 2 GeV −2 , CS is in units of 10 2 GeV −4 , and Λ is in units of GeV. For the NLO LECs, the first uncertainty is statistical (originating from the uncertainties of the BESIII data) and the second is systematic (originating from chiral truncation uncertainties, see the main text for details).