Exotic hadron bound state production at hadron colliders

The non-relativistic wave function framework is applied to study the production and decay of exotic hadrons, which can be effectively described as bound states of other hadrons. Employing the factorized formulation, with the help of event generators, we investigate the production of exotic hadrons in multiproduction processes at high energy hadron colliders. 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 to be multiquark states and/or bound states of other ingredient hadrons in many 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 about the inner structure can be deduced. In general, the production processes are more complex.
At the same time, studies of the production processes also provide information on the cross section, rapidity and transverse momentum distributions, etc., of the relevant particles at a specific collider, which can help experimentalists to set proper triggers and cutoffs for their measurements [5]. A good example is the 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 men-tioned above are crucial for studies of exotic hadrons with a specific detector at the LHC.
Furthermore, the direct production of exotic hadrons in the multiproduction process of high energy scattering can set a crucial point for understanding of the hadronization mechanism. Since exotic hadrons always refer to states with more than three constituent quarks (here we do not discuss hybrids or glueballs), one feasible way of 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 out in Refs. [8,9], the produced color-singlet (anti)quark system eventually transitions 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 is beyond the present approach of calculation. This leaves 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, i.e. quark confinement. The introduction of multiquark states sets a challenge for hadronization models dealing with the transition from a color-singlet (anti)quark system to a hadron system. From experiments, the production of general mesons and baryons is dominant, i.e., where 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 the number of constituent quarks in the hadron, in the case of large number of quarks produced, such as in high energy nuclear collisions, the more constituent quarks a hadron contains, the larger production rate one gets. So to regain unitarity, one needs a 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 present knowledge is not enough to judge how many kinds of multiquark states there are and how they 'share' the total probability of ε, one cannot predict the production rate of a specific multiquark state. We can suspect, though, that if there are many 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 a 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., the hydrogen atom). When the number of hadrons produced is 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 the large production rate of X(3872) [10][11][12][13][14][15][16][17][18][19][20][21][22] and X(5568) [23-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 considered as possible 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 respect to the hadron mass (almost of the order of the charm quark), so the hadronic molecule is in principle a NR system. Besides its application in the investigation of positronium, the NR wave function is also used for heavy quarkonium production and decay, generally referred to as the 'color-singlet model'.
Non-relativistic quantum chromodynamics (NRQCD) implies that the quark pair in a 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 the hadron as basic degree of freedom, there is no problem of color confinement because every object is a 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 mesons and baryons in high energy hadronic collisions, it is impossible to construct the effective Lagrangian. In Ref. [30], we suggested using a 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 respect 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 the 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 the ingredient hadrons A and B, are relevant. This means that the details of the structure of the bound state and the production of the ingredient hadrons in the whole phase space are not fully used. So this framework provides a benchmark for the information least sensitive to 083106-2 the models and details of the strong dynamics at a momentum scale of several MeV to several tens of 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 the present approach in Refs. [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 giving further discussion.

Non-relativistic wave function formulation for pp→A+B+X →H(A,B)+X
For the case where 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 the 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: 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 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 states, as illlustrated in the following examples. 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 discussion). 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 u(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 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: where 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: For the S-wave case, L,L z ;S,S z |J,J z =δ J,S δ Jz ,Sz in Eq. (5). The wave function in Eq. (3) is symmetric and only the first part in Eq. (9) contributes, hence where the coordinate space Schrödinger S-wave function . 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 terms do, 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 simplification of the factorized formulation at cross section level. We need 1) the production cross section rather than the amplitude of a 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).
2) the Schrödinger wave function to describe the bound state. It can be obtained from potential models of the relevant bound system.
The way to get the distribution of free pair production can sometimes rely on effective theories based on QCD, to calculate the amplitude then the cross section. For example, in Ref. [30], the processes of J/ψ decay into the baryon, anti-baryon 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) collisions, the effective Lagrangian cannot be constructed. In this case, to employ event generators to give the distribution is the only practical way. We will show that the generators can also be applied to the P -wave case.

Studying free hadron pairs 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 over 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 calculate 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 the Shandong Quark Combination Model [34], etc., for the case that A and B are both on shell. The advantage of the above framework that we employ is that only the on shell case is considered, so that numerical calculation with an event generator is plausible. The value of Eq. (12) describes the correlation of the two hadrons/clusters (A and B) in phase space. For the hadron case, by proper integration over the components of P H and/or q, the resulting correlations can be directly compared with data and serve for tuning the parameters. Because of 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 aŝ 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 k = √ −q 2 is exactly the absolute value of the 3-relative momentum | p A − p B |.
Employing the event generator, one gets then extrapolate 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, Equation (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 the 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 discussed in the next section. Figure 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 extrapolation result, or equivalently use 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 the derivative.
Related to this extrapolation, we would like to address that the formulation we employ here is well factorized. The NR wave function is universal. However, this factorization is not a priori correct, since we do not have an (at least effective) field theory of hadron production. So there is no way to prove the factorization. We can only use the generator as a numerical effective theory, which itself is finite anyway. The validity of the factorization can only be checked for certain particles, 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 can also be seen as a benchmark for the study of higher order corrections. As shown from Refs. [10][11][12][13][14][15][16][17][18][19][20][21][22], there is an intuitive point of view in which 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, detailed investigation of 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.
Analysis of 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 point of view, this energy and momentum resolution is determined by the experiment, especially the two-hadron correlation data, which is yet not available for high energy hadronic collisions. On 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 deeper problem and determines to what scale and energy resolution our factorized formulation has 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 dependence with 'coarse resolution', which should be more consistent with high energy processes.

Numerical results and discussion
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 collisions. Nowadays, a lot of hadron bound states have been studied theoretically and experimentally. Some of the bound states and the literature relevant to their corresponding NR wave functions are listed in Table 1. 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 used to predict the production in other phase space regions, other collision energies and processes.
Recently, the D0 collaboration announced that they found 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 a multiquark state directly produced in the multi-production process of a high energy collision, rather than from hadron decays. The bottom flavor here is decisive. Because of the bottom flavor and the mass, it is hardly possible to be produced from the 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 the 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 it. Furthermore, theoretical analyses, like the kinematic distributions of the signal particles from the decay of X(5568), are useful for vari-083106-6 ous detectors. Though recently searches for X(5568) by LHCb and CMS have given negative results [25,26], one must also consider the impacts of the different cuts and trigger conditions used by those detectors compared with those used by D0.
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 a bound state of two hadrons (see the 14th line of the middle column of Table 1). 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 the 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 no other collaborations have yet reported 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 because the two clusters are required to be close to each other in phase space for combination. The above formulations in Section 2 give vanishing relative momentum, so this framework generally predicts a softer p T spectrum: the behavior looks like a single particle at small p T . However, the production rate is suppressed for large p T . The reason is very clear. If we assume random correlations between two ingredient particles, the larger their single transverse momentum, the smaller the probability of the two particles having almost vanishing relative momentum, and hence they are much more suppressed. This is a very typical property of the NR formulations, in contrast to the fragmentation spectrum, where the spectrum is harder for more massive particles. For high energy and then high p T processes, the question of whether the NR formulations are still valid or not is not settled yet. Quarkonium production has 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 useful information. 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 harder spectra.
Furthermore, both charm and bottom quarks are heavy. Relevant processes can be calculated with perturbative 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, the kaon, 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 detect: 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 pion can eliminate the possibility that the X c is produced from the decay of the bottom. Of course just by keeping or not keeping 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 multiproduction in hard processes and the wave function of a certain bound state are universal. The only difference between 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 in center-of-mass energy, parton distribution functions and rapidity region. But the behavior of X(5568) is quite similar to that of B s . Fatal violation of the factorization will lead to more difference in 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 a hadron bound state as shown in Table 1, and employing 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 momentum of the 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. N Ψ (2S) is the number of Ψ (2S) for the corresponding collision. dN X(3872) dyd 2 p T d 3q is calculated by the event generator. σ exp X is the prompt production cross section of X(3872) from 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, and the universality of the wave function and hence factorization is confirmed. It is considered in some studies [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 the multiproduction process, as has been discussed in Section 1. If it is effectively dealt with as a hadron molecule, its wave function at the origin should be much smaller than that of X(3872). A similar investigation can also be applied to Z b .
The distributions of the rapidity and transverse momentum are shown in Figs. 6 (with p T > 5 GeV/c) and 7 for the corresponding bound states listed in Table 1. All the bound states have similar shapes and properties, like X(5568) and X c [24]. Figures 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 Ref. [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 the produced pions not being energetic, e.g., only around 10% of the signal pions have 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 that of the X(5568) for the same  momentum. This means that, whether Tevatron or LHC, in both central and large rapidity regions, the signal pions are more energetic and can be detected.
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 results for bound states with larger spin. This is a good interplay arena between hadron polarization and exotic hadron state studies.