Production of charm-strange hadronic molecules at the LHC

We explore the inclusive hadroproduction of the D_{s0}^*(2317), D_{s1}(2460) and D_{sJ}(2860) states at the Large Hadron Collider under the assumption that these hadrons are S-wave meson-meson molecules. In addition, the D_{s2}(2910), a predicted bound state of the D_2(2460)K system, is also discussed. We first derive a factorisation formula for the production rates based on effective field theory. Then we make use of two MC event generators, Herwig and Pythia, to simulate the production of pairs of charmed mesons and kaons. Using effective field theory to handle the final state interaction among the meson pairs, we give an estimate of the inclusive production rates for these particles at the order-of-magnitude accuracy. Our results show that the cross sections for the pp -->D_{s0}^*(2317) and pp -->D_{s1}(2460) at the LHC are at O(1 \mu b) level, while the ones for the pp -->D_{sJ}(2860) and the pp -->D_{s2}(2910) are smaller by about one order of magnitude. Such estimates suggest that these four exotic states could be copiously produced at the LHC. A study of these states at the LHC will thus provide valuable information on hadron spectroscopy as well as hadron interactions.


Introduction
Heavy-light mesons can provide a unique window to study heavy quark dynamics and non-perturbative QCD in the presence of a heavy quark. In recent years, a number of the charm-strange mesons have been discovered at various experimental facilities. These observations have not only enriched the particle zoo of the ordinary cs family, but also have provided interesting candidates for exotic states beyond thē qq scenario.
One of the most interesting exotic candidates is the D * s0 (2317) first discovered in 2003 by the BABAR Collaboration in the inclusive D s π 0 mass distribution from e + e − annihilation data [1]. Later, the CLEO Collaboration discovered the D s1 (2460) in the D * s π 0 final state [2], followed by the discovery of the D sJ (2860) by the BABAR Collaboration in the DK channel [3].
The heavy-light hadrons are conveniently classified in spin multiplets, the members of which are degenerate in the heavy quark limit m Q → ∞, labeled by the value of the total angular momentum of the light degrees of freedom s . For instance, the pseudoscalar charmed meson D s and the vector charmed meson D * s belong to the multiplet with s P = 1 2 − . Given that the mass difference between D s1 (2460) and D * s0 (2317) is not large, it is likely that the D * s0 (2317) and D s1 (2460) are spin partners. An important feature of the above mentioned two charm-strange mesons is that their masses are far below the potential model predictions for the lowest states with the same quantum numbers, e.g. in the Godfrey-Isgur quark model [4]. Various explanations were proposed for these states, for reviews, see Refs. [5,6]. Of particular interest is the hadronic molecule assumption under which the D * s0 (2317) and D s1 (2460) are DK and D * K bound states, respectively [7][8][9][10][11]. In this picture, the approximate equality of the mass splittings M D s1 (2460) − M D * s0 (2317) 142 MeV and M D * − M D 143 MeV, as can be seen in Table 1, is a natural consequence of heavy quark spin symmetry [12]. One can also appeal to the parity doublet assumption [13][14][15] to explain this equality. Being below the DK and D * K thresholds, the D * s0 (2317) and D s1 (2460) both have narrow widths, less than 3.8 MeV and 3.5 MeV, respectively [16]. The small widths arise from radiative decays and isospin-violating hadronic decay modes, due to either the η − π 0 mixing [17] or the mass differences between the neutral and charged charmed meson and kaons [18][19][20][21]. These decay modes have been intensively explored under different interpretations of the internal structure, and it was suggested that the isospin breaking decays can be used to distinguish the hadronic molecule picture from the others [18][19][20][21][22]. So far the experimental resolution is not high enough to measure the width of the hadronic width of the order of 100 keV predicted in the hadronic molecule picture. Yet, the planned PANDA experiment is expected to be able to perform the measurement [23]. However, there is support of the hadronic molecule picture from some lattice QCD calculations of the charmed meson-light meson scattering lengths [22,24].
It was proposed in Ref. [25] that the leading order interaction between an excited heavy meson with a small decay width 1 and the kaon should be the same as that for the ground state charmed mesons. This is a consequence of chiral symmetry when the kaons are regarded as pseudo-Goldstone bosons of the spontaneous breaking of chiral symmetry. Analogous to the D * s0 (2317) and D s1 (2460) as the DK and D * K bound states, more kaonic bound states are predicted. Especially, the properties of the predicted D 1 (2420)K bound state agree with the experimental data for the D sJ (2860) very well, including the ratio of branching fractions B(D sJ (2860) → DK)/B(D sJ (2860) → D * K) which is difficult to be accounted for in other models. Its spin partner, a D 2 (2460)K bound state, was predicted to have a mass and a width around 2910 MeV and 10 MeV, respectively. This state has quantum numbers J P = 2 − and can decay into the D * K and D * s η. It will be called D s2 (2910) throughout the paper.
So far, most of the available experimental and theoretical investigations of the exotic candidates in the charm-strange sector have been focused on the spectrum and decays or the production in the e + e − collisions including B meson decays. In order to decipher the nature of the D sJ states, more accurate data and new processes involving these states would be equally important. These processes include heavy ion  Table 1. Masses of the D * s0 (2317), D s1 (2460) and D sJ (2860) taken from the Particle Data Group [16], and the relative energy E 0 of the exotic states to the corresponding thresholds (in units of MeV).
collisions [26,27] and hadroproduction. At low energies, the forthcoming PANDA experiment has a particular focus on the hadroproduction of D * s0 (2317) state [28], as it is supposed to be able to measure the decay width of the D * s0 (2317) at the 100 keV level. At hadron colliders with a very high collision energy like the LHC, the charm quark will be abundantly produced via the QCD processes, which may serve as an ideal platform for the study of the D sJ states. For instance, the LHCb Collaboration has reported a measurement of the decays of the D * s1 (2710), which is the first radially excited state of the D * s with J P = 1 − , and the D sJ (2860) into the D + K 0 S and D 0 K + final states with small statistical uncertainties [29].
In this paper, we will provide a theoretical exploration of the hadroproduction of the D * s0 (2317), D s1 (2460), D sJ (2860) and the D s2 (2910) at the LHC. Our approach is similar to that used in studying the production of the X(3872) [30][31][32] and other heavy quarkonium-like states [33,34], while for different treatments on heavy quarkonium production, see Refs. [35][36][37]. Our calculation is based on the assumption that the four states are S-wave hadron molecules consisting of a kaon and a nonstrange charmed meson. Under this assumption, the charmed meson-kaon pairs will be produced first and these D sJ states are then formed through the final state interaction between them. We will make use of two Monte Carlo (MC) event generators, Herwig and Pythia, to simulate the production of pairs of the charmed meson and the kaon. Then, using effective field theory (EFT) to handle the final state interaction between the kaon and charmed meson, we will derive an estimate of the production rates for these particles at the order-of-magnitude accuracy. As we will show in the following, the prospect to observe these states at the LHC from the accumulated data is rather promising and thus an experimental analysis is urgently called for. This paper is organised as follows. In Sec. 2, an overview of the effective field theory description of the kaonic bound states will be presented. Based on the MC event generators, we will derive the production rates of the molecular states in Sec. 3, while the numerical results are presented subsequently in Sec. 4. We will summarize this work in Sec. 5. Sketch of the resummation of the amplitudes. Here, T denotes the total scattering amplitude, V is the leading order amplitude and G is the loop integral. The heavy meson H and kaon are represented by solid and dashed lines, respectively.

Kaonic bound states in effective field theory
As the isospin of these exotic states is I = 0, the assumed molecular structures can be explicitly written as where H denotes the D, D * , D 1 and D 2 charmed mesons. All the considered bound states contain the light pseudoscalar kaon, which can be considered a Goldstone boson of the spontaneous chiral symmetry breaking from SU(3) L ×SU(3) R to SU(3) V . The interaction between the Goldstone bosons and the heavy mesons can be described by the heavy flavor chiral perturbation theory [38][39][40], which combines the chiral and heavy quark symmetries in the light and heavy quark sectors, respectively. The leading order contact interaction between a Goldstone boson and a heavy meson comes from the kinetic term in the Lagrangian for the heavy mesons with the derivatives being chirally gauged where H denotes the heavy meson field, and the chiral gauge derivative is given by and F is the pion decay constant in the chiral limit. Since we will only consider the leading order interaction, F = 92.2 MeV will be used. From the Lagrangian, one can derive the leading order scattering amplitude. 2 One finds that among all the channels for the light meson-heavy meson S-wave scattering, the kaon-charmed meson one with isospin I = 0 is the most attractive. The tree level scattering amplitude of the process HK → HK reads where s, t and u are the Mandelstam variables with the constraint s + t + u = 2(m 2 H + m 2 K ) for on-shell particles. The S-wave amplitude is projected out for onshell particles by using is the energy of the heavy meson, k is the momentum in the center-of-mass frame, and we have used the fact that the heavy meson is highly nonrelativistic in the energy region of interest. Since the interaction is strong, and a perturbative expansion cannot account for the subthreshold states such as the D * s0 (2317), we need to sum up the amplitudes up to infinite orders, as shown in Fig. 1. This can be done by using the Bethe-Salpeter equation (BSE), which is an integral equation where p 2 = s, and we keep the dependence of the amplitudes on the center-ofmass momenta in the initial and final states explicitly. We will follow the approach developed in Ref. [44] and take the on-shell S-wave amplitude V 0 (s) as the kernel of the BSE. This is an approximation, and its justification can be found by using dispersion relations [45,46] or taking the off-shell contribution into account [47]. 3 Then, the integral equation is simplified to an algebra equation. The resummed S-wave amplitude reads Notice that both here and in Eq. (2.4) we have neglected the polarizations of the charmed mesons with spin 1 or 2. The justification comes from the fact that at the energy region of interest these charmed mesons are highly nonrelativistic (for details we refer to Appendix A). In Eq. (2.7), G(s) is the two-meson scalar loop function The integral is divergent, and the divergence can be regularized using a subtraction constant. Treating the heavy meson nonrelativistically, we get an analytic expression Figure 2. The mechanism considered in the paper for the inclusive production of the D sJ as a HK bound state in proton-proton collisions. Here, X denotes all the produced particles other than the H and K in the collision.
for the renormalised loop function [41] where µ is the regularization scale, a(µ) is the subtraction constant which is scaledependent rendering the loop function scale-independent, A bound state of two mesons is associated with a pole of the resummed scattering amplitude below threshold on the real axis in the first Riemann sheet of the energy plane. Following Ref. [25], we fix the value of a(µ) to reproduce the mass of the D * s0 (2317) at 2317.8 MeV, which gives a(µ = 1 GeV) = −3.84, and use it in all other channels. Certainly, different values of the subtraction constant correspond to different pole positions in a given channel. This allows us to vary the constant to get different binding energies so as to investigate the impact of binding energy on the cross sections later on.

Hadroproduction of kaonic molecules
In principle, the D sJ states can couple to other components such as the cs or a [cq][sq] tetraquark if they have the required quantum numbers even the dominant component of their wave function is a kaon-charmed meson bound state. All the components contribute to the production of the D sJ states. In general, it is a process-dependent question to conclude which one is more important. Here we will assume that the production of the D sJ will occur through generating the charmed meson kaon pairs first which will form the D sJ states afterward. The mechanism is shown in Fig. 2. This mechanism is valid only when the binding energy of the bound state is small so that its constituents H and K are only slightly off-shell.

Factorization of the near-threshold production rate
In Ref. [31], a factorization formula for the DD * production rate in the near-threshold region was used, which reads where k is the center-of-mass momentum of the constituents, µ is the reduced mass, M[D * 0D0 (k) + X] is the amplitude for the D * 0D0 inclusive production, T is the D * 0D0 scattering amplitude which accounts for the final state interaction (FSI) and X means all the other particles in the inclusive process. The phase space integration runs over all the rest particles and the one for the total momentum of the D * 0D0 pair.
In this FSI method, the cross-channel rescattering with other comoving particles is neglected. In fact, as pointed out in Refs. [42,43], the presence of such particles could largely affect the so-calculated cross section. The factorization formula in Eq. (3.1) is based on the an assumption that the momentum dependence of the production amplitude is mainly given by T (k) in the near-threshold region, and M/T is insensitive to k. In general, however, it is not clear how M/T depends on the momentum at all. Especially, if the produced particles contain one or more of the lowest-lying pseudoscalar mesons (pions, kaons and the η-meson) which are the pseudo-Goldstone bosons of the spontaneous breaking of chiral symmetry of QCD, the near-threshold production amplitude without the FSI should have a momentum dependence dictated by chiral symmetry. At first sight, it is not clear whether there exists a factorization formula analogous to Eq. (3.1) or not. Here, we will show that such a formula indeed exists as a consequence of chiral symmetry and if we use the same on-shell approximation as that in Eq. (2.7).
The general differential cross section formula for the inclusive HK production reads where k is the three-momentum in the center-of-mass frame of the HK pair, µ is the reduced mass, and M[HK(k) + X] is the production amplitude. As illustrated in Fig. 3, we separate the inclusive production amplitude of the HK pair into two parts: one is due to the direct production to be denoted by Γ, and the other includes the FSI described by the scattering amplitude T . We thus have an integral equation .  Inclusive production of the HK pair in proton-proton collisions. Here Γ denotes the direct production vertex, T is the resummed HK scattering amplitude, and X denotes all the produced particles other than the HK pair.
whereĜ is the operator for the Green's function of the HK system. Similarly, the BSE given in Eq. (2.6) can be written aŝ Since we are only considering the production of the HK pair in the near-threshold region, the production amplitude should be dominated by the S-wave. Taking the same on-shell approximation as that in the resummed scattering amplitude given in Eq. (2.7), in the near-threshold region we have where we have used the center-of-mass momentum as the variable of all the functions. From Eq. (2.5) we know that V 0 (k) ∝ E K = m 2 K + k 2 . But what is the momentum dependence of the direct production amplitude Γ 0 (k)? The answer comes from a chiral symmetry analysis which is applicable in the near-threshold region.
Since we are only interested in the HK pair and do not care about the details of the other particles in the process, we may parameterize all the other particles involved in the inclusive production process by an external source field S which has the same quantum numbers as the HK pair. This means that we neglect the cross-channel rescattering of the charmed meson and kaon with other particles in the final state. The spirit of this kind of treatment was already used to investigate the ππ(KK) system in the decays J/ψ → φππ(KK) [49]. In the chiral limit when the up, down and strange quarks are massless, the coupling of the Goldstone bosons (pions, kaons and the eta) to any other fields has to be in a derivative form which ensures that the interaction strength vanishes at threshold as required by the Goldstone theorem. Thus, at leading order in the momentum expansion we have the Lagrangian for the near-threshold HK production 4 L prod = c S ∂ µ H ∂ µ K. (3.8) All the short-distance physics has been parameterized into the coefficient c and the source field S. Here, we have made an implicit approximation that in the final states of the inclusive production, there is no other soft chiral particles other than the kaon in the HK pair so that we can neglect the interaction of the kaon with them.
Although there are two derivatives in each term, this Lagrangian is of order O(k), with k being a momentum much smaller than the typical hadron scale Λ χ ∼ 1 GeV. The reason is that the dominant part of the heavy meson four-momentum is the mass of the heavy meson, and thus ∂ µ H is dominated by its temporal component. An insertion of the light quark masses will give higher order corrections. Therefore, at leading order in the chiral expansion, the direct production amplitude is proportional to the kaon energy Γ 0 (k) ∝ E K . The same energy dependence occurs in V 0 (k) as given in Eq. (2.5). Together with Eq. (3.6), this means that the production amplitude can be factorized into the product of a constant C and the HK scattering amplitude in the near-threshold region M(k) = C T 0 (k). (3.9) In brief, we have derived the factorization formula for the production of the HK pair which is valid in the near-threshold region where the constant C may take a value of M[HK(k) + X]/T 0 (k) for any k provided that k Λ χ . It is related to the coefficient c in Eq. (3.8) through |C| = |c|F 2 /2. This is the analogue of the factorization used in studying the DD * production in Ref. [31], and is applicable to the near-threshold production of a pair of a heavy meson and a pseudo-Goldstone boson. It allows for the separation of the long-distance and shortdistance contributions in the amplitudes for the production of the molecules. The latter is the same for the processes pp → HK and pp → D sJ , while the long-distance factor can be deduced from the scattering amplitude given in Eq. (2.7).

Production of the charm-strange hadronic molecules
The cross section for the D sJ production is (3.11) Since in Eq. (3.2) the integrated phase space dφ HK+X already contains the part of the total momentum of the HK pair, the phase space integration in the above equation is the same as that in Eq. (3.2). The resummed scattering amplitude T 0 (k) contains information about the generated hadronic molecules because these D sJ states are the bound state poles of T 0 (s). Thus, it is straightforward to extend the factorization formula Eq. (3.10) to the case of the inclusive D sJ production provided that these states are produced through intermediate HK pairs. The recipe is to replace the HK → HK scattering amplitude T 0 (s) in Eq. (3.10) by the amplitude for the process HK → D sJ , which is given by the effective coupling constant for the D sJ HK vertex in the vicinity of the D sJ pole. We have where C is the same constant as that in Eq. (3.10). The effective coupling constant is given by the residue of the transition matrix element at the pole We have checked that the same equation will be obtained if we use the approach adopted in Ref. [31] which uses the Migdal-Watson theorem [50,51] and the unitary relation for the scattering amplitude. By varying the subtraction constant in the expression of loop function G, we can get different binding energies. Therefore the cross section for the hadronic molecule in Eq. (3.12) is dependent on the subtraction constant, and indeed on the binding energy. With a smaller binding energy, we found that the production rate of the molecules gets smaller. This conclusion is in agreement with Ref [31], in which the universal scattering amplitude was adopted. From the scattering length a = (2µE B ) −1/2 with E B the binding energy, it means that the bound state with a larger scattering length is more difficult to be produced.

Estimating the cross sections using Monte Carlo event generators
As a phenomenological and successful tool that has been used in many other processes, MC event generators are able to simulate the hadronization of partons produced in QCD processes, and therefore provide an estimate of the pp → HK inclusive cross sections. Two commonly adopted programs are Pythia [52] and Herwig [53]. However in the above event generators, the FSI effect which governs the momentum dependence close to threshold, as depicted in the scattering amplitude T 0 (s) in Eq. (2.7), is not incorporated. Thus, the MC cross section corresponds to the case without near-threshold FSI, and can be approximately expressed in terms of the vertex given in the leading order Lagrangian in Eq. (3.8), where K HK ∼ O(1) is introduced because of the overall difference between MC simulation and the experimental data, while for an order-of-magnitude estimate we can roughly take K HK 1. Therefore, the coupling constant c can be determined using Eq. (3.14). Substituting the result into the Eq. (3.12), we obtain the cross section for the D sJ : Although the expression contains explicitly a factor of 1/(E 2 K k 2 ), this factor is completely cancelled by (dσ[HK(k)]/dk) MC and thus a momentum-independent value is obtained.

Results
To form a molecular state, the constituents must move nearly collinear as a multiquark system and thus have a small relative momentum. Such configurations can be realized in an inclusive 2 → 2 QCD process, while the multiquark final states can be produced by soft parton shower radiations. The dominant partonic process for the production of the D sJ is gg → cc, as the gluon density at the LHC energy is much larger than those for quarks. In addition, the process qq → cc will also be included in this analysis.
Using Pythia and Herwig, we have generated 10 8 events which contain a pair of the charm and anti-charm quarks. These events are then analyzed by the Rivet library [54] in order to pick out the charmed-meson kaon pair with a small invariant mass. Since the kaons in these events are mostly produced from the soft gluon emission by the heavy charm quarks, they tend to move together with the charmed hadrons. This portion of the events will contribute significantly to the formation of the molecules. To match the capability of the detectors, we have also implemented the cuts on the transverse momentum of the meson pair p T > 5 GeV, and the rapidity |y| < 2.5 and 2.0 < y < 4.5 for the ATLAS/CMS (denoted as LHC for simplicity)  In the MC calculation, the kinematic cuts used are |y| < 2.5 and p T > 5 GeV, which lie in the phase-space regions of the ATLAS and CMS detectors. Here, we have averaged the events in different charged channels, e.g. D 0 K + and D + K 0 for the DK case. and LHCb detectors, respectively. In Fig. 4, we show the differential cross sections (histograms in this figure) versus the center-of-mass momentum of the constituent mesons up to 0.35 GeV for the inclusive processes pp → DK, D * K, D 1 K, and D 2 K at the LHC with √ s = 8 TeV. Since the production of charmed mesons D 1 and D 2 is not included in Pythia, we only show the results from Herwig.
As analyzed in Sec. 3.1 and shown in Eq. (3.14), one should be able to approximately describe the production mechanism in the MC calculation using the leading order effective Lagrangian given in Eq. (3.8), and the resulting differential cross section is proportional to k 2 E 2 K . To validate this feature, we fit to these distributions with Eq. (3.14) up to 350 MeV for the center-of-mass momentum. The fitted results are shown as curves in Fig. 4. Note that the shape of the cross section in Eq. (3.14) is completely fixed, and the fitting procedure only results in a normalization constant, which is proportional to c. From this figure, one sees that in most of the cases the MC results can be well described. The only exception is the D * K from Herwig where there exists a clear peak at the bin between 150 MeV and 200 MeV. This peak can be attributed to the resonance D s1 (2536) because it decays into D * + K 0 and D * 0 K + with a center-of-mass momentum of 149 and 167 MeV, respectively. Note that this resonance was included in Herwig but not in Pythia. In principle, this resonating contribution can be built in the scattering amplitudes using a coupled-channel formalism, which will improve the consistency with the MC simulations. However, since we are only interested in an order-of-magnitude estimate, we refrain from doing such an analysis.
With Eq. (3.15) and the differential distributions in Fig. 4, we can obtain the production rates of the hadronic D sJ molecules. These results for the cross sections for the inclusive processes pp → D * s0 (2317), D s1 (2460), D sJ (2860) and D s2 (2910) (in units of µb) at the LHC with √ s = (7,8,14) TeV are shown in Table 2. Results outside (inside) brackets are obtained using Herwig (Pythia if applicable). Here the rapidity range |y| < 2.5 has been assumed for the LHC detectors (ATLAS and CMS), while the rapidity range 2.0 < y < 4.5 is used for the LHCb. The differences between the values from Pythia and Herwig are caused by the different hadronization mechanisms to form the charmed and kaon mesons. From this table, one sees that the cross sections for the pp → D sJ (2860) and the pp → D s2 (2910) at the LHC are at 10 −1 µb level, while the ones for the pp → D * s0 (2317) and the pp → D s1 (2460) are larger by roughly one order of magnitude. Such large production rates suggest a promising perspective for searching for these four exotic states at the LHC.
The D s1 (2460) has a large decay branching fraction into D s γ and thus can be reconstructed from the K + K − π + γ final states. With the available decay branching fractions [16], we find that the cross section for the process pp → D s1 (2460) → D + s γ → K + K − π + γ reaches O(10 7 fb), and this would yield O(10 8 ) events when considering the integrated luminosity of 22 fb −1 from ATLAS and CMS in 2012 [55,56]. The D sJ (2860) mainly decays into DK and D * K, while the D s2 (2910) has a large partial width into D * K. These two hadrons can be reconstructed in hadronic final states. In particular, based on the 1 fb −1 data accumulated in 2011, the LHCb Collaboration has used the DK final state to reconstruct the D sJ (2860) [29], where about 3×10 4 D 0 K + events from the D sJ (2860) were observed. Using the cross section shown in Table 2 and assuming that the D sJ (2860) → D 0 K + is of O(10 −1 ) supported by the Babar data [57], we predict about O(10 6 ) D 0 K + events to be generated, which is about two orders of magnitude larger. However, since the detection efficiency of the experiment was not published, a direct comparison is not possible. Currently, we may only conclude that our prediction is not in conflict with the measurement. A more quantitative comparison is expected when the efficiency corrected data is available in future.
The search for the D sJ states at hadron colliders depends on the non-resonant background contributions. To investigate this issue, we take the D s0 (2317) to be constructed in the D s γ final state as an example. To be conservative, we use the cross section σ(pp → D ± s + anything) as an upper bound for the background. The ATLAS collaboration has provided a measurement of this cross section at √ s = 7 TeV [58]: where the p T (D ± s ) > 3.5 GeV and the pseudo-rapitidity η < 2.1. Our results in Table 2 show that the cross section of the pp → D s0 (2317) at √ s = 7 TeV is about 2 µb. Using the integrated luminosity in 2012, 22 fb −1 [55, 56], we have an estimate for the signal/background ratio where 5% is the current upper bound for the branching fraction of the D s0 (2317) → D s γ [16]. It is worthwhile to point out that our theoretical results have to be modified somewhat due to the mismatch in kinematics, however, we believe that the above estimate indicates a great potential for observing the discussed molecular states at the LHC.

Summary
Our understanding of hadron spectroscopy will be largely improved by studies of exotic states that may defy the conventional models of cs meson spectroscopy. Great progress has been made in the past decades, see Ref. [5,6] and references therein. One of the most important aspects is the discrimination of a compact multiquark configuration and a loosely bound hadronic molecule. Many theoretical predictions exist on the decays of the some of these states, and they can be used to distinguish different scenarios including the cs, compact tetraquark states and hadronic molecules. However, only few of the decay branching fractions have been measured. Experimental data with high statistics are desired at this moment. In our work, we have explored the hadroproduction of D * s0 (2317), D s1 (2460), D sJ (2860) and the predicted D s2 (2910) states at the LHC under the assumption that these hadrons are S-wave hadron molecules. We have made use of two MC event generators, Herwig and Pythia, to simulate the production of the charmedmeson kaon pairs. Together with effective field theory to handle the final state interaction among the meson pairs and neglect their interactions with other particles, we have derived an estimate of the production rates for these particles at the order-ofmagnitude accuracy. Our results show that the cross sections for the pp → D sJ (2860) and the pp → D s2 (2910) at the LHC are at 10 −1 µb level, while the ones for the pp → D * s0 (2317) and the pp → D s1 (2460) are larger by roughly one order of magnitude. Thus, these states can be copiously produced at the LHC, and measurements in future would be able to test the molecular description of the above states and also the production mechanism. Such measurements are very important to gain deeper insights into the hadron interactions, in particular the interactions between heavy and light mesons.

Acknowledgments
This work is supported in part by the DFG and the NSFC through funds provided to the Sino-German CRC 110 "Symmetries and the Emergence of Structure in QCD", by the NSFC (Grant No. 11165005) and by the EU I3HP HadronPhysics3 Project.

A Polarization vectors/tensors in the scattering amplitudes
In fact, the scattering amplitude in Eq. (2.4) has a factor of the product of polarization vectors (if the scattered heavy meson is the D * or D 1 ) or tensors (for the D 2 ), and the tree-level amplitude should be of the form V (s, t, u) (λ 1 ) · (λ 2 ) , (A.1) where · means i i if the scattered heavy meson is the D * or D 1 and ij ij for the and that ij is symmetric and traceless. As a result of Eq. (A.2), the inner product of the polarization vectors/tensors can be factorized out and becomes an overall factor of the resummed amplitude It looks different from the analogous equation in Refs. [9,59], where the denominator of the resummed amplitude reads 1 − G(s)V (s) [1 + q 2 cm /(3M V )] for the pseudoscalar meson-vector meson scattering, where q cm is the size of the momentum of the vector meson in the center-of-mass frame, and M V is the vector meson mass. In the nonrelativistic limit, the additional factor 1 + q 2 on /(3M V ) is reduced to 1, and one gets the same equation as Eq. (A.4).