Exotic Hadron Bound State Production at Hadronic Colliders

The non-relativistic wave function framework is applied to study the production and decay of the exotic hadrons which can be effectively described as bound states of other hadrons. Employing the factorized formulation, we investigate the production of exotic hadrons in the multiproduction processes at high energy hadronic colliders with the help of event generators. This study provides crucial information for the measurements of the relevant exotic hadrons.


Introduction
Recently, more and more new exotic hadron states, e.g., the XYZ mesons, have been observed. They are assumed as multiquark states and/or as bound states of other ingredient hadrons in lots of theoretical investigations [1][2][3][4] . Such investigations can be done not only via their decay processes, where the branching ratios and distributions of the decay products can be studied, but also via their production processes (from the decay of the heavier particles or directly from the multiproduction processes). In both cases, the more complex the process is, the more information of inner structure can be drawn. In general, the production processes are more complex.
At the same time, the studies on the production processes also provide the information of the cross section, rapidity and transverse momentum distributions, etc., of the relevant particles on a specific collider, which can help the experimentalists to set the proper triggers and cutoffs for the measurements [5] . A good example is Large hadron Collider (LHC) where various detectors cover a large rapidity range.
They can be used to study exotic hadron production in B-decays [6] , as well as direct production in the multiproduction process of high energy hadronic collisions and nuclear collisions. So the distributions mentioned above are crucial for the studies on exotic hadrons with a specific detector at LHC. Furthermore, the direct production of the exotic hadrons in the multiproduction process of high energy scattering can set a crucial point for the understanding of the hadronization mechanism. Since the exotic hadrons always refer to the states with more than three constituent quarks (here we do not discuss hybrids or glueballs), one of the feasible ways for understanding their production mechanism is to employ the combination model [7] to combine the necessary constituent quarks into the relevant hadron. However, in any hadronization process, as pointed in Refs. [8,9] , the produced color-singlet (anti)quark system eventually transits to various hadron states (the mesons, baryons and beyond) with the total probability exactly 1: Here we introduce the unitary time-evolution operator U to describe the hadronization process. For the quark state |q and the corresponding hadron state |h , the matrix element U hq = h|U |q describes the transition amplitude. U hq is determined by (low energy) Chromodynamics (QCD) but beyond the present approach of calculation. This leaves the space for various hadronization models to mimic the transition process. The unitary operator U reflects the fact that there are no free quarks in the final states of any high energy process, e.g., the so-called quark confinement. The introduction of multiquark states sets a challenge for the hadronization models dealing with the transition from color-singlet (anti)quark system to the hadron system. As a matter of fact from experiments, the production of general mesons and baryons is dominant, i.e., here B,B and M denote baryon, antibaryon and meson, respectively. If the exotic hadrons are produced, ε could be a small but non-vanishing value. Since the production rate is proportional to the quark density to the power of constituent quark number in the hadron, in the cases of large number of quarks produced such as those in high energy nuclear collisions, the more constituent quarks a hadron contains, the larger production rate one gets. So to regain the unitarity, one needs the special 'combination function' which reflects the confinement and may be related with the whole system rather than the several quarks to be combined into a specific hadron [8,9] . Since the present knowledge is not enough to judge how many kinds of multiquark states there are and how they 'share' the total probability of ε, one can not predict the production rate of a specific multiquark state. What we can suspect, though, seems that if there are a lot of kinds of multiquark hadrons, each only shares a small part of the small ε. So the production rate of each is almost vanishing. However, if one of the exotic states is the bound state of other hadrons, i.e., its production can be taken as from the combination of mesons and/or (anti)baryons, there is generally no straightforward unitarity constraint as above on its production rate. Unlike quarks, hadrons are not confined. They can be either free, or bound with other hadron(s), even lepton(s) (e.g., hydrogen atom). As a matter of fact, in the cases that the number of produced hadrons are large, such as in high energy ion collisions, the familiar hadron bound states, such as deuterons, α particles, etc., have been observed. So it is natural to investigate the productions of the 'hadronic molecule' relevant to exotic states. The measurements and explanations of large production rate of X(3872) [10][11][12][13][14][15][16][17][18][19][20][21][22] and X(5568) [23][24][25][26] are also good examples. Their large production rates lead to insights on the investigation of their structure and production mechanism, especially the colour and spin structures [11,27] .
For those XYZ states possibly considered as hadronic molecules, we can describe them with the framework of various Non-Relativistic (NR) effective theories, especially the NR wave function method, and concentrate on their inclusive production. The ingredient hadrons in the hadronic molecule are loosely bound, hence the relative momentum between them is fairly small with respective to the hadron mass (almost of the order of charm quark), so the hadronic molecule is in principle a NR system.
Besides the application of the investigation on positronium, the NR wave function is also used for the heavy quarkonium production and decay, generally referred to as the 'color-singlet model'.
The Non-Relativistic Quantum Chromodynamics (NRQCD) implies that the quark pair in color octet could also transfer into a color-singlet hadron, which is referred to as the color-octet model and used to explain the production rate and transverse momentum of prompt quarkonium in hadronic collisions. But for the case of hadron as basic degree of freedom, there is no problem of color confinement because every object is color-singlet. Hence the relevant complexity [28,29] is eliminated. If the system of bound hadrons is properly modeled and the NR wave function is obtained, the NR wave function framework can be used for various decay and production processes of hadronic molecules. This method has been applied to the near threshold enhancements of mass spectra in the J/Ψ → γpp and J/Ψ → pΛK − channels of the production of bound states X(pp) and X(pΛ) [30] . In that note, the decay to the corresponding ingredient hadrons is described by an effective Lagrangian. However, for the direct production of the mesons and baryons at high energy hadronic collisions, it is impossible to construct the effective Lagrangian. In Ref. [30] , we suggested to employ the general event generator to extract the cross section of the ingredient hadrons. In fact, this is one of the advantages of the NR wave function framework. One can expand the amplitude with respective to the relative momentum between the ingredients because it is relatively small. Thus the cross section is factorized. Only the cross section (rather than the amplitude) of the ingredient hadron production is needed, which can be fixed by fitting the correlation data of the relevant hadrons.
In the next section, we will review and list the formulations for the calculation of hadronic molecule production rate in the NR wave function framework, taking pp scattering at LHC, say pp → A + B + X → H(A, B) + X, as an example. Based on these formulations, only the NR wave function (and/or its derivatives) at the origin and correspondingly, only the square of the absolute value of the production amplitude (and/or its derivatives) of ingredient hadron A and B, are relevant. This means that the details of the structure of the bound state and the production of the ingredient hadron in all the phase space, are not fully used. So this framework provide the benchmark of the information least sensitive to the models and details of the strong dynamics in the momentum scale of several MeV to several ten MeV. At the same time, if the universality of the wave function is confirmed, the factorization is well established. In Section 3, the calculation of the cross section of the ingredient hadron pair of A and B will be investigated, extrapolating to the vanishing relative momentum. We make a comparison with present approach in literatures [10][11][12][13][14][15] . Then the rapidity and transverse momentum distributions of a series of bound states like KK, D * N , DD * (X(3872)), Λ cΛc , Σ cD * , X(5568) and X c [23,24] , etc., are taken as illustrations in Section 4, besides further discussions are given.
2 Non-Relativistic wave function for- For the case that a bound state H is well described by the two ingredients A and B, the only difference between the amplitude of the bound state and that of the free particles lies in that the wave function rather than the plane-waves of free particles is adopted to describe the bound state. The process pp → A+B+X → H(A, B)+X is illustrated in Fig. 1, and the corresponding invariant amplitude is: (3) This formulation is valid in the rest frame of H(A, B), the bound state of the ingredient hadrons A and B. In the above equation, k is the relative 3-momentum between A and B in the rest frame of H(A, B). M( k) is the invariant amplitude for the free (unbound) A and B production. The factor 1/ m A m B m A +m B comes from the normalization of the bound state to be 2E H V just as a single particle. The normalization for the phase space wave function Φ( k) is [31] From the above equation, it is obvious that if we know the analytical form of the wave function as well as the free particle invariant amplitude, we can calculate the amplitude of the bound state simply by integrating the relative momentum k. In practice, certain decomposition and simplification will be taken for concrete 2S+1 L J state, as examples in the following. The resulting formulations are covariant.
The free particle invariant amplitude M( k), as is shown by the part after the product sign in Fig. 1, can be obtained from the corresponding Feynman diagrams with the Feynman rules from 'standard' model or other effective theories (see the following discussions). The general form of this invariant amplitude is B j (p B )Ô ji (P H , k)A i (p A ). In the following we take a spin-1/2 fermion-antifermion pair as an example to illustrate some details. In this case it is a Lorentz scalarū(p B )Ô(P H , k)v(p A ). HereÔ(P H , k) represents a collective product of γ−matrices, spinors, etc., i.e., the internal lines, vertices, and other external lines except those of the A and B. Hence in the rest frame of H(A, B), the invariant amplitude of the bound state production is: where M is the mass of the bound state. The above equation has two C-G coefficients (spin 1/2 -1/2 and spin-orbital angular momentum) for the specified angular momentum of H. The decomposition of the two spin-1/2 addition is: with the polarization vector in the rest frame. Since k is small in the rest frame, the above formulations can be expanded around k = 0 up to terms linear in k for the cases of S-wave and P-wave (L = 0, 1) as: (3) is symmetric and only the first part in Eq. (9) has contribution, hence where the coordinate space Schrödinger S-wave func- . The radial wave function R can be obtained when the interaction potential between the ingredients is determined. For L = 1, the wave function is anti-symmetric, the first part in Eq. (9) does not contribute, but the second and third term give, so the amplitude for the 2S+1 L J = 1 P 1 bound state is: (11) The above formulations can all be found from the literature, e.g., [32] , and are widely applied, e.g., in quarkonium production.
As demonstrated above, by the recognition that the relative momentum between the ingredients of the NR bound system is small, the expansion of the amplitude leads to the simplification of factorized formulation at cross section level. We need • the production cross section rather than the amplitude of free particle with vanishing relative momentum, and the ingredient particles on shell, projected onto the special quantum number of the definite bound state (mainly the isospin and/or angular momentum state).
• Schrödinger wave function to describe the bound state, which can be obtained from the potential models of the relevant bound system.
The way to get the distribution of free pair production sometimes can rely on the effective theories based on QCD concept, to calculate the amplitude then the cross section. For example, in [30] , the processes of J/ψ decay into the baryon, antibaryon and one or more pseudo-scalar particles, e.g., kaons or pions, with the pair of fermions forming a bound state is studied. The relevant decay process J/Ψ → pΛK − is described by a simple effective interaction Lagrangian. We calculated the partial width of the process J/Ψ → pK −Λ to determine the coupling constant G. However, for very complex processes such as the production of hadrons in high energy pp(p) collision, the effective Lagrangian can not be constructed. In this case, to employ the event generators to give the distribution is the only practical way. We will show that the generator can also be applicable to the P-wave case.

Studying free hadron pair with event generators
The free pair cross section p(p 1 )p(p 2 ) → A(p A ) + B(p B ) + X can be expressed as: Here the average is on various spin states, and the proper initial flux factor 1/F and phase space integral are needed.Ô is the amplitude of production of two free ingredient particles (with vanishing relative momentum and proper angular momentum state). It is not possible to be calculated directly with some effective quantum field theory/model when the initial state is (anti) protons and A and B are hadrons or clusters [24] . However, it can be obtained with an event generator such as PYTHIA [33] or equivalently Shandong Quark Combination Model [34] , etc. for the case that A and B are both on shell. It is the advantage that in the above framework we employ, only the on shell case is considered, so that the numerical calculation with event generator is plausible. The quantity of Eq. (12) describes the two hadrons/clusters (A and B) correlation in the phase space. For the hadron case, by proper integral on components of P H and/or q, the resulting correlations can be directly compared with data and serve for tuning the parameters.
Since the special physical picture of the nonrelativistic framework, it is only valid in the rest frame of the two ingredient particles. One can define the following covariant space-like relative momentumq It is clear that in the rest frame of A and B (H(A, B)) where p A + p B = 0,q = (0, k) and the k = √ −q 2 is exactly the absolute value of the 3-relative momentum Employing the event generator, one gets then extrapolates to the special case k = 0. Numerically, one can take an average around k = 0 for the above quantity. Then we get, up to the kinematic factors as for the covariant form, Eq. (15) is exactly the differential cross section of the bound state H(A, B) divided by |Ψ(0)| 2 .
There are several very basic facts supporting the extrapolation. First of all, the amplitude and cross section are analytical in phase space. Any practical generator should reproduce this property, and any ultraviolet divergence is not present. Secondly, the study of strong interaction is complex because of the SU(3) non-Abelian interaction, but its simulation has one simplicity: All particles taking part in the strong interactions are massive, which eliminates the infrared singularities.
Here for simplicity, we only consider the S-wave case. We take two examples: KK and D * N which are to be discussed in the next section. Fig. 2 shows the distribution 1 N dN dyd 3q of these two pairs, at rapidity y ∼ 0. These distributions are quite smooth, so one can get a reasonable extrapolating result, or equivalently using the average value around k ∼ 0.
The smooth line also indicates that for the distribution, the derivative distribution around k ∼ 0 is a relatively small value with respect to the distribution itself, which then can be taken as higher order and neglected. This means that for the P-wave production, the calculation can be much simplified from the numerical calculation of derivative.
Related with this extrapolation, we would like to address that, this formulation we employ here is well factorized. The NR wave function is universal. However, this factorization is not a priori correct, since we do no have an (at least effective) field theory of the hadron production. So one has no way to prove the factorization. We can only use the generator as a numerical effective theory, itself is finite anyway. The validity of the factorization can only be checked for certain particle, for various processes and/or energies (see discussion in section 4). This depends on the concrete production mechanism and the structure of the relevant particles. This also can be seen as a benchmark for the study, for higher order corrections. As shown from [10][11][12][13][14][15][16][17][18][19][20][21][22] , there is an intuitive point of view, high energy collisions tend to produce pairs of hadrons with a very high relative momentum. This is indeed one of the naive points against the production of hadron molecules at these facilities. At the same time, the detailed investigation on the dependence of the range of q (orq), though depending on the wave function/model of the particle such as X(3872), can give more information on the whole picture. The analysis on the factorization also leads to a more subtle consideration, i.e., the concrete and detailed wave function depends on the energy and momentum resolution. From the most practical view point, this energy and momentum resolution is deter-mined by the experiment, especially the two hadron correlation data which is yet not available at high energy hadronic collisions. One the other hand, it is also determined by the interplay between the high energy hadron production mechanism and the static structures of the hadron as well as the hadron molecule. The latter is a more deep problem and determines to what scale and energy resolution our factorized formulation get the best validity. Here we only give another extrapolation at a quite small energy resolution for reference (see Fig. 3). In the following analysis we employ the former ones with 'coarse resolution', which should be more consistent with the high energy processes.

Numerical results and discussions
In the above section, we have studied the production of the ingredient particle pair of A and B. In this section, without the explicit calculation of |Ψ(0)| 2 , we investigate the rapidity and transverse momentum distributions (see Fig. 4, Fig. 6 and Fig. 7 respectively) of some possible exotic hadron states at √ s = 8 TeV in pp collision. Nowadays, a lot of hadron bound states have been studied theoretically and experimentally. Some of the bound states and the literatures relevant to their corresponding NR wave functions are listed in Table I. In general, the wave function can be obtained by solving the relevant potential model. Employing a recent example [24] , we show the wave function at the origin can also be obtained by fitting the available data, and then is used to predict the production in other phase space regions, other collision energies and processes.
Recently, the D0 collaboration announced that they found the evidence of a new state X(5568) [23] , with four flavors in this hadron. If it is really a particle, this will be quite remarkable, as the first solid evidence of multi-quark state directly produced in the multi-production process of high energy collision, rather than from hadron decays. The bottom flavor here is very decisive. Because of the bottom flavor and the mass, it is hardly possible produced from decay of a heavier hadron. Another important fact is that the measurement [23] has given the production ratio ρ of the yield of X(5568) to the yield of the B 0 s meson in two kinematic ranges, 10 < p T (B 0 s ) < 15 GeV/c and 15 < p T (B 0 s ) < 30 GeV/c. The results for ρ are (9.1 ± 2.6 ± 1.6)% and (8.2 ± 2.7 ± 1.6)%, respectively, with an average of (8.6 ± 1.9 ± 1.4)%. If we assume that B 0 s π ± is the dominant decay mode of this new state, this large production rate itself first of all excludes the possibility of decay from heavier particles like B c , and is difficult to be understood by various general hadronization models [24] , such as String fragmentation model [45] , cluster model [46] and Combination model [34] . By employing these general hadronization models, one can hardly raise to a larger production ratio compared with the D0 experiment, since there are no other plausible parameters to tune. So it seems that the X(5568) has a curious structure and unique production mechanism different from the other particles produced. This needs a special theoretical framework to deal with. Furthermore, the theoretical analysis, like the kinematic distributions of the signal particles from the decay of X(5568), are useful for various detectors. Though recently the searches of X(5568) by LHCb and CMS give negative results [25,26] , one must also consider the impacts of different cuts and trigger conditions from the D0 adopted by those detectors.
Exploring the inclusive production formulations proposed in Sections 2 and 3, we can calculate the cross section of the new B 0 s π ± state for various collision processes as well as energies, taking it as the bound state of two hadrons (see the 14th line of the middle column of Table I). By fitting the D0 data at the collision energy 1.96 TeV within the phase space region 10 < p T < 30 GeV/c and |η| < 3, we get the |Ψ(0)| 2 , which is then used to predict the production of X(5568) at LHC for the whole rapidity region, since |Ψ(0)| 2 is independent of collision energy and phase space region. As an example, Fig. 4 shows the transverse momentum distributions and pseudo-rapidity η (p T > 5 GeV/c) distributions for proton-proton collisions at √ s = 8 TeV. All these calculations are useful to study this unconfirmed problem since yet no other collaborations report positive results.
The X(5568) transverse spectrum of p T is softer than that of B s , as demonstrated in Fig. 5 and indicated from the experiment [23] . This is from the fact that the two clusters are required to be close to each other in the phase space for combination. Realized in the above formulations in Section 2, is the relative Table 1. Hadron-hadron states and related molecules.

Hadron-hadron state
Related molecules References for Ψ(0) KK f 0 (980), a 0 (980) [35] ΛΛ Y (2175), η(2225) [36] DD * X(3872) [37][38][39] Z c (3900) [39] D * sD * s Y (4140) [40] BB * Z b (10610) [39] D * D * Z c (4020), Z c (4025) [41] B * B * Z b (10650) [41] D * 0 p Λ c (2940) [42] Λ cΛc Y (4260), Y (4360) [43] Λ bΛb Y (10890) [43] Σ cD * P c (4380) [44] Σ * cD * P c (4450) [44] BK X(5568) [24] DK X c [24] momentum vanishing. So that this framework generally predicts a softer p T spectrum: The behavior looks like a single particle in small p T . However, the production rate is suppressed for large p T . The reason is very clear. If we assume randomly correlations between two ingredient particles, the larger the single transverse momentum of them, the smaller the probability of the two particles with almost vanishing relative momentum, and hence much more suppressed. This is a very typical property of the NR formulations, in contrary to the fragmentation spectrum, the more massive, the harder. For high energy and then high p T processes, the question whether the NR formulations are still valid or not is not settled yet. The quarkonium production had given some implications. A possible way is to employ the Bethe-Salpeter wave function [47,48] . In this formulation, all possible relative momenta of the ingredients are taken into consideration. Once the wave function shows that the probability of large relative momentum between the ingredients is vanishing, the Bethe-Salpeter description might go to the NR formulations. However, the rich Dirac structures in the Bethe-Salpeter wave functions can also introduce sounding informations. In Refs. [10][11][12][13]15] , the phase space wave function is directly used, with a certain significant momentum region. In these kind of models, the bound states can have a harder spectra. Furthermore, both charm and bottom quarks are heavy. Relevant processes can be calculated with per-turbative QCD in the exactly same way once the different values of the masses are taken into account. So, when one replaces the bottom quark by the charm quark in the X(5568) and assumes that the structure and the production mechanism do not change, the production of the 'new D ± s π ± state' (X c ) can be also studied. Since the reduced mass is mainly determined by the mass of the light ingredient hadron, say, Kaon here, the wave function at the origin can be taken as the same. Therefore, the cross section of X c for both Tevatron and LHC should be completely determined. We illustrate the rapidity and transverse momentum distributions of X c in Fig. 6 (a) and Fig. 7 (a). If X c exists, it is not difficult to be detected: the D ± s can be determined from the D ± s → φπ channel, by proper 3-charged particle tracks from the vertex displaced from the primary one; then this reconstructed D ± s can be combined with a proper charged particle track considered as π from the primary vertex to give the invariant mass distribution to look for the resonance. If the K 0 s is well measured, the D ± s can also be reconstructed from the 2K channel and then combined with the π from the primary vertex. This kind of pions can eliminate the possibility that the X c is produced from the decay of bottom. Of course just by keeping or not this restriction, one can preliminarily investigate X c from multi-production or from weak decay.
All the above discussions are based on the factorization formulas, i.e., the mechanism of the multi- production in hard process and the wave function of a certain bound state are universal. The only difference between the pp and pp machines lies in the center of mass energy and the parton distribution functions. Comparing Fig. 4 and 5, the transverse momentum distributions can be different because of the difference on center of mass energy, parton distribution functions and rapidity region. But the behavior of X(5568) are quite similar as that of B s . Fatal violation of the factorization will lead to more difference of the X(5568) production between these two kinds of collisions.
To further demonstrate this framework, we make a global analysis of the prompt X(3872) data. Contrary to the case of X(5568), it is a well-established particle [17,20,49,50] , but without definite evidence for its structure [10][11][12][13]15] . However, as an exotic hadron state, its high production rate is also remarkable. Here, taking it as hadron bound state as shown in Table I, employ the available data from CDF and CMS [17,20] , and the event generator, we calculate the following value: Here y and p T are the rapidity and transverse momen-tum of bound state respectively. y 0 = 0.6 and p T 0 = 5 GeV/c for CDF; y 0 = 1.2 and p T 0 = 10 GeV/c for CMS. The N Ψ(2S) is the number of Ψ(2S) for the corresponding collision. The dN X(3872) dyd 2 p T d 3q is calculated by event generator. The σ exp X is the prompt production cross section of X(3872) of experimental data, and here we use σ exp X = 3.1 nb for CDF and σ exp X = 1.06 nb for CMS. We find R CM S /R CDF is of order of 1, the universality of the wave function and hence factorization is approved. It is considered in some literatures [37][38][39] that Z c shares some similarity with X(3872). Because Z c is apparently a tetraquark, its production rate is relatively small in multiproduction process, as has been discussed in Section 1. If it is effectively dealt with as hadron molecule, its wave function at origin should be much smaller than that of X(3872). The similar investigation can also be applied to Z b .
The distributions of the rapidity and the transverse momentum are shown in Fig. 6 (with p T > 5 GeV/c) and 7 for the corresponding bound states listed in Table I. It is obvious that all the bound states have the similar shape and property, like X(5568) and X c [24] . Fig. 6 and 7 indicate separately that the rapidity distribution can extend to |y| > 5, and the transverse momentum distribution still has some detectable value when p T ∼ 40 GeV. The above results are useful for various detectors. For example, based on our calculation, one can go further to estimate the kinematic distributions of the signal particles which are from the decay of X(5568) and can be directly detected. In [24] , we give the transverse momentum -total momentum distribution for the signal pions from the decay process X → B s + π in the LHCb detector range (2 < η < 5). The mass difference between X(5568) and B s + π is small, and the pion mass is small. These facts lead to that the produced pions are not energetic, e.g., only around 10% of the signal pions with transverse momentum larger than 0.5 GeV/c (the requirement of the relevant measurement by the LHCb Collaboration [25] ). However, since the mass of the X c is around half of the X(5568), it has a larger Lorentz boost factor about two times of that of X(5568) for the same momentum. This means that wherever Tevatron or LHC, in both central and large rapidity regions, the signal pions are more energetic to be detectable.
Additionally, if the present generators are well modified for the production of polarized hadrons when the polarization experimental data are sufficient, the method in this paper can also give the result of the bound state with larger spin. This is a good interplay arena between hadron polarization and exotic hadron states studies. ¯ *