Analytical calculation of heavy quarkonia production processes in computer

This report is devoted to the analytical calculation of heavy quarkonia production processes in modern experiments such as LHC, B-factories and superB-factories in computer. Theoretical description of heavy quarkonia is based on the factorization theorem. This theorem leads to special structure of the production amplitudes which can be used to develop computer algorithm which calculates these amplitudes automatically. This report is devoted to the description of this algorithm. As an example of its application we present the results of the calculation of double charmonia production in bottomonia decays and inclusive the $\chi_{cJ}$ mesons production in pp-collisions.


Introduction
It is well known that the study of heavy quarkonia lead to a great boost in the understanding of QCD [1].This is connected to the fact that heavy quarkonia contains all manifestations of QCD.At the same time heavy quarkonia are nonrelativistic systems what makes the physics of heavy quarkonia much simpler that that for other mesons.For this reason so far the study of heavy quarkonia in different experiments is very important and intensive.
The study of heavy quarkonia production processes is based on the fact that the mass of heavy quark inside quarkonia is much larger than characteristic QCD energy scale ( m Q ≫ Λ QCD ).This fact leads to the appearance of small parameter -the relative velocity of quark-antiquark pair inside quarkonia v ≪ 1.For the charmonia mesons this parameter is v ∼ 0.3, for bottomonia mesons v ∼ 0.1.The appearance of the small parameter separates different energy scales involved in quarkonia physics and allows to build a rigorous formalism for studying quarkonia properties which is called nonrelativistic QCD (NRQCD) [2].

NRQCD formalism
NRQCD is based on the separation of different energy scales through the factorization theorem where T is the amplitude of the process e + e − → η c γ1 , C n are the Wilson coefficients which parametrize the physics of small distances, η c | Ôn |0 are different NRQCD operators which contribute to the process.The sum in ( 1) is infinite and it runs over all possible operators with correct quantum numbers.In particular, for the production of the η c meson the following operators contribute Ô ∼ χ + ψ, χ + D 2 ψ, χ + H σψ, ...
However, if we restrict the study of the process under consideration by some power of relative velocity O(v n ) then only finite number of the NRQCD operators contibute to (1).In particular, at the leading order approximation in relative velocity only the η c |χ + ψ|0 operator contributes.If we omit the operators which assumes that quarkonia consist of more than quark-antiquark pair ( for instance, χ + H σψ) then factorization theorem (1) can be rewritten in the following form where the matrix elements q n which control relativistic corrections can be written as follows The first few matrix elements v n for the 1S, 2S, 1P charmonia mesons where calculated within QCD sum rules in papers [3,4,5,6].The Wilson coefficients C n which parametrize short distance contribution to the amplitude can be expanded in the strong coupling constants α s The first aim of the present report is the description of the algorithm of analytical calculation of the coefficient c n for any process.

Light cone expansion formalism
If a heavy quarkonium production process contains energy scale E h which is much larger than the masses of final mesons, one can apply light cone expansion formalism (LCF) to study such process [7,8].For instance, LCF can be applied if the energy of e + e − beam is much larger than the masses of charmonia mesons, what takes place at B-factories.The presence of high energy scale E h , which is of order of characteristic energy of the hard exclusive process, allows one to apply factorization theorem for the amplitude of the process (1).As in the case of NRQCD, the coefficients C n describe partons production at small distances, the matrix elements η c |O n |0 describe hadronization of the partons which takes place at large distances.The sum is taken over all possible operators O n .For instance, the operators are few examples of the operator O n .Actually, there are infinite number of the operators O n that contribute to the pseudoscalar meson production.
The cross section of hard exclusive process can be expanded in the inverse powers of the high energy scale To determine if some operator contributes to a given term in 1/E h expansion one uses the concept of the twist of this operator [9].Thus only the leading twist -the twist-2 operators Qẑγ 5 Q, Qẑγ 5 (zD)Q, Qẑγ 5 (zD) 2 Q, ... 2 contribute to the leading term in expansion (6).From this Feynman diagrams for gg → χ cJ g subprocess one sees that already at the leading order approximation infinite number of operators contribute to the cross section.Nevertheless, it is possible to cope with infinite number of contributions if one parametrizes all the twist-2 operators by the moments of some function φ(x) where q is the momentum of the pseudoscalar meson η c , f η is the constant which is defined as η c (q)| Qγ µ γ 5 Q|0 = if η q µ , x is the fraction of momentum of the whole meson η c carried by quark.The function φ(x) is called the leading twist distribution amplitude (DA).One can think of it as about the Fourier transform of the wave function of the meson η c with lightlike distance between quarks.Using the definition of this function, factorization theorem (1) can be rewritten as where H(x) is the hard part of the amplitude, which describes small distance effects.This part of the amplitude can be calculated within perturbative QCD.As it was noted, the leading twist DA parametrizes infinite set of the twist-2 operators.This part of the amplitude describes hadronization of quark-antiquark pair at large distances and parametrizes nonperturbative effects in the amplitude.DA of charmonia mesons were studied in papers [3,4,5,6].
It should be noted that formula (8) resums the contributions of all twist-2 operators.If the meson η c is a nonrelativistic meson, formula (8) resums relativistic corrections to the amplitude T .It should be also noted that if one takes into account the renormalization group running of the DA φ(x), formula (8) resums all leading logarithmic corrections to the amplitude of the process under considerations.
The second aim of the present report is the description of the algorithm of analytical calculation of the hard part of the amplitude H(ξ) for any process.

Hadronic production of P -wave charmonium states
Another example is hadronic production of P -wave charmonium states at high energies.In NRQCD formalism χ cJ mesons are considered as nonrelativistic quark-antiquark pairs in color singlet (CS) or color octet (CO) states The hadronization probability of each component are described by long distance matrix elements |R ′ (0)| 2 , O S , O P , that are usually determined from fit of experimental data.
In high energy limit gluon-gluon interactions give the main contributions to hadronic production of χ c -mesons (see diagrams shown in Fig. 1).To obtain the analytic expressions for the amplitudes of these processes one can write down the amplitude for cc pair production in the reaction gg → ccg and project it on the corresponding state.For example, in the case of color-singlet χ c0 meson one have the expression where P µ , q µ , and ǫ Sν are the momentum of final charmonium meson, relative momentum of quarks and spin polarization vector.

Analytical calculation
The main idea of the algorithm for analytical calculation within NRQCD and LCF is based on the factorization theorem (1).This theorem allows us to make factorization in the calculation of the amplitude which consists in the following.The amplitude is devided into two part.The first part parametrizes the small distance contribution.The calculation of this part is reduced to the calculation of some set of Feynman diagrams.The second part of the amplitude describes hadronization of quarks.In the calculation it is reduced to the some projection operators.In order to get total amplitude these two parts are convoluted by Dirac and Color indices.To carry out the calculation we use the Mathematica package FeynCalc [10] for analytical calculation in high energy physics.
As an example let us consider analytical calculation of H(p)γ(k)|J µ em |0 matrix elements, which can be used to build the whole amplitude of the processes e + e − → H(p)γ(k), H = η c , χ c0 , χ c1 , χ c2 .The calculation will be carried out within LCF.Analogously one can repeat the algorithm for NRQCD.
The program which realizes the algorithm starts from the kinematics  [13] [11] Exp.

hk = GS[k]; hA = GS[A];
After the definitions one can calculate the first part of the amplitude -Feynman diagrams of the processes involved Note that thus calculated diagrams are very simple.If the corresponding diagrams are complicated and the number of diagrams is large then it is better to use package FeynArt ( which is a part of FeynCalc [10]), where it is possible to calculate the diagrams automatically.
The second part of the amplitude which parametrizes hardonization of quarks is represented by projection operators • the η c meson: P βα = (pγ 5 ) βα f P 4 • the χ c0 meson: P βα = (p) βα f χ0 4 • the χ c1 meson: P βα = (pγ 5 ) βα f χ1 4 • the χ c2 meson: It is comfortable to write these projection operators in terms of the array Project Project = {fP/4*hp.GA [5], fc0/4*hp, fc1/4*hp.GA [5], fc2/4*hp} Now two parts of the amplitude can be convoluted by Dirac and color indices At the end of the program one gets correct amplitudes of the processes involved. .3. Applications 3.1.Double charmonia production in exclusive bottomonia decays Using the algorithm described above in [11] all leading twist double charmonia production in exclusive bottomonia decays were calculated.The calculation was done both for NRQCD and LCF.It should be stressed that one program calculated the amplitudes and branching ratios for approximately 30 different processes of bottomonia decays.We are not going to present all the results obtained in paper [11] since it will take a lot of space.The only result that is presented in this report (see Table 1) is the branching ratio for the processes measured at Belle collaboration [14].It is seen from Table 1 that our results (fourth column) are in agreement with the results obtained at Belle experiment (fifth column).It will be highly desirable to trace the origin of the disagreement in the results of papers [11,12,13] 3.2.Heavy quarkonia production in hadronic collisions Using the formalism described above we have calculated also the cross sections of partonic reactions gg → χ cJ g with the help of FeynCalc and Redberry [15] packages.On the basis of these expressions one can calculate the cross sections of hadronic reactions where x 1,2 are momentum fractions of initial gluons and f g (x) is the hadronic structure function.
Using the obtained in this way theoretical predictions for the transverse momentum distributions of the color singlet and octet cross sections one can determine from fit to existing experimental data the values of long distance matrix elements |R ′ (0)| 2 , O S , O P .Such analysis was performed by our group in papers [16,17].Presented in these papers numbers are given in Table 2 and the comparison with experiment is shown in Fig. 2 and Fig. 3.One can see, that good agreement is achieved.

Conclusion
This report is devoted to the analytical calculation of heavy quarkonia production processes in modern experiments such as LHC, B-factories and superB-factories in computer.Theoretical description of heavy quarkonia is based on the factorization theorem.This theorem leads to special structure of the production amplitudes which can be used to develop computer algorithm which calculates these amplitudes automatically.In this report we described this algorithm.As an example of its application we presented the results of the calculation of double charmonia production in bottomonia decays and inclusive the χ cJ mesons production in pp-collisions.The authors would like to thank the ACAT organization committee for invitation to the conference and hospitality.The work was financially supported by Russian Foundation for Basic Research (grant #10-00061a), the grant of the president of Russian Federation (#MK-3513.2012.2),and FRRC grants.Transverse momentum distribution of J/ψ production in radiative χ cJ decays at CDF in comparison with experimental data [18].Solid lines correspond to the total pp → χ cJ +X → J/ψ +X cross section in different choices of µ 2 (upper corresponds to µ 2 = m 2 T , and lower three lines to other schemes).Dashed, dotted and dot-dashed lines correspond to CS, S-wave octet and P -wave octet contributions respectively (all at µ 2 = m 2 T ).[20] (•), [21] (▽) and [22] (filled rectangles).Solid lines correspond to parameter values presented in Tab. 2, while dashed lines show predictions with only CS (lower curve) or CO (upper curve) contributions taken into account.

Figure 2 .
Figure2.Transverse momentum distribution of J/ψ production in radiative χ cJ decays at CDF in comparison with experimental data[18].Solid lines correspond to the total pp → χ cJ +X → J/ψ +X cross section in different choices of µ 2 (upper corresponds to µ 2 = m 2 T , and lower three lines to other schemes).Dashed, dotted and dot-dashed lines correspond to CS, S-wave octet and P -wave octet contributions respectively (all at µ 2 = m 2 T ).

Table 2 .
CS and CO model parameters for different values of the scale µ 2 .