Bottomonium production and polarization in the NRQCD with $k_T$-factorization. III: $\Upsilon(1S)$ and $\chi_b(1P)$ mesons

The $\Upsilon(1S)$ meson production and polarization at high energies is studied in the framework of the $k_T$-factorization approach. Our consideration is based on the non-relativistic QCD formalism for a bound states formation and off-shell production amplitudes for hard partonic subprocesses. The direct production mechanism, feed-down contributions from radiative $\chi_b(mP)$ decays and contributions from $\Upsilon(3S)$ and $\Upsilon(2S)$ decays are taken into account. The transverse momentum dependent (TMD) gluon densities in a proton were derived from the Ciafaloni-Catani-Fiorani-Marchesini evolution equation and the Kimber-Martin-Ryskin prescription. Treating the non-perturbative color octet transitions in terms of multipole radiation theory, we extract the corresponding non-perturbative matrix elements for $\Upsilon(1S)$ and $\chi_b(1P)$ mesons from a combined fit to transverse momenta distributions measured at various LHC experiments. Then we apply the extracted values to investigate the polarization parameters $\lambda_\theta$, $\lambda_\phi$ and $\lambda_{\theta\phi}$, which determine the $\Upsilon(1S)$ spin density matrix. Our predictions have a reasonably good agreement with the currently available Tevatron and LHC data within the theoretical and experimental uncertainties.

Since it was first observed, the production of heavy quarkonium states in high energy hadronic collisions remains a subject of considerable theoretical and experimental interest [1,2]. These processes are sensitive to the interaction dynamics both at small and large distances: the production of heavy (c or b) quarks with a high transverse momentum is followed by a bound states formation with a low relative quark momentum. Accordingly, a theoretical description of these processes involves both perturbative and non-perturbative methods, as it was proposed in the non-relativistic QCD (NRQCD) [3][4][5][6]. However, it is known that the NRQCD at the next-to-leading order (NLO) accuracy meets difficulties in a simultaneous description of all the collider data in there entirety (see also discussions [7][8][9][10][11][12]). In particular, it has a long-standing challenge in the J/ψ and ψ(2S) polarization and provides an inadequate description [13][14][15][16][17][18] of the recent η c production data taken by the LHCb Collaboration at the LHC [19]. One of possible solutions of the problems mentioned above, which implies a certain modification of the NRQCD rules, has been proposed recently [20]. As it was shown, the approach [20] allows one to describe well the recent data on the production and polarization of the entire charmonia family. The bottomonium production, namely Υ(nS) and χ b (mP ) mesons, provides an alternative laboratory for understanding the physics of the hadronization of heavy quark pairs. Due to heavier masses and a smaller quark relative velocity v (in a produced quarkonium rest frame), these processes could be even a more suitable case to apply the double NRQCD expansion in QCD coupling α s and v. The NLO NRQCD predictions for the Υ(nS) production at the LHC were presented [21][22][23]. Of course, it is important to apply also the approach [20] to the bottomonium family.
Our present work continues the line started in the previous studies [24,25]. We have considered there the inclusive production of Υ(3S), Υ(2S), χ b (3P ) and χ b (2P ) mesons and now come to Υ(1S) and χ b (1P ) mesons. The motivation for the whole business has been already given [24,25]. Below we present a systematic analysis of the CMS [26][27][28], ATLAS [29] and LHCb [30][31][32][33] data on the Υ(1S) and χ b (1P ) production collected at √ s = 7, 8 and 13 TeV (including the different relative production rates) and we extract from these data non-perturbative matrix elements (NMEs) for the Υ(1S) and χ b (1P ) mesons. Then we make predictions for polarization parameters λ θ , λ φ , λ θφ (and a frame-independent parameterλ), which determine the Υ(1S) spin density matrix and compare them to the currently available data [34,35]. As it is known, the feed-down contributions from χ b (2P ), χ b (3P ), Υ(2S) and Υ(3S) decays give a significant impact on the Υ(1S) production and polarization, so studies [24,25] are important and necessary for our present consideration. Another important issue concerns the relative production rate σ(χ b2 )/σ(χ b1 ) recently measured by the CMS [28] and LHCb [33] Collaborations. This ratio is sensitive to the color singlet (CS) and color octet (CO) production mechanisms and provides information complementary to the study of the S-wave bottomonium states.
In the present note we follow mostly the same steps as in [24,25]. So, to describe the perturbative production of the bb pair in the hard scattering subprocesses we apply the k Tfactorization approach [36,37], which is mainly based on the Balitsky-Fadin-Kuraev-Lipatov (BFKL) [38] or Ciafaloni-Catani-Fiorani-Marchesini (CCFM) [39] gluon evolution equations. A detailed description and discussion of the different aspects of the k T -factorization can be found in the reviews [40]. As usual, we see certain advantages in the ease of including into the calculations a large piece of higher order pQCD corrections taking them into account in the form of transverse momentum dependent (TMD), or unintegrated, gluon densities in a proton. Our consideration is based on the off-shell gluon-gluon fusion subprocesses representing the true leading order (LO) in QCD: where J = 0, 1 or 2 and the four-momenta of all particles are given in the parentheses. The color states taken into account are directly indicated. Both initial gluons are off mass shell, that means that they have non-zero transverse four-momenta k 2 1 = −k 2 1T = 0, k 2 2 = −k 2 2T = 0 and an admixture of longitudinal component in the polarization four-vectors (see [36,37] for more information). The corresponding off-shell (k T -dependent) production amplitudes contain projection operators [41] for spin and color, that guarantee the proper quantum numbers of the final state bottomonium. Following the ideas [20], to describe the nonperturbative transformations of the color-octet bb pairs produced in hard subprocesses into observed final state mesons we employ the classical multipole radiation theory (where the electric dipole E1 transition dominates [42]) under the key physical assumption that the lifetime of intermediate color-octet states is rather long. According to [20], only a single E1 transition is needed to transform a P -wave state into an S-wave state 1 , whereas the transformation of the color-octet S-wave state into the color-singlet S-wave state is treated as two successive E1 transitions 3 S     states with J = 0, 1, 2. An essential consequence of the idea above is the nonconservation of the spin momentum S z during the transformation of color octet 3 S (1) 1 one. In fact, the intermediate P -wave state is a state with definite total momentum J and its projection J z rather than a state with definite L z and S z . To describe the formation of the intermediate state, we have to contract the electric dipole transition amplitude [42] (which by its own conserves S z ) with Clebsch-Gordan coefficients |L, S, J, J z L, L z , S, S z | which are symmetric with respect to L z and S z . Then the resulting expression comprises both: the terms containing * (S z [ 3 S 1 ]) · (S z [ 3 P J ]) and the terms containing Consequently, there is no direct transfer from the initial spin polarization to the final spin polarization (see [20] for more information).
Below we apply the gauge invariant expressions for quarkonia production and decay amplitudes implemented into the Monte-Carlo event generator pegasus [43]. The derivation steps are explained in [24,25] in detail.
According to the k T -factorization prescription, to calculate the cross sections of a considered process one has to convolute the partonic cross sectionσ * (related with an off-shell production amplitude) and TMD gluon densities in a proton f g (x, k 2 T , µ 2 ): where x 1 and x 2 are the longitudinal momentum fractions of initial off-shell gluons, φ 1 and φ 2 are their azimuthal angles and µ is the hard interaction scale. Following [24,25], we have tested several sets of TMD gluon densities in a proton. Two of them (A0 [44] and JH'2013 set 1 [45]) were obtained from the CCFM equation where all input parameters were fitted to the proton structure function F 2 (x, Q 2 ). We have applied the TMD gluon densities obtained within the Kimber-Martin-Ryskin (KMR) prescription [46], which provides a method to construct the TMD quark and gluon distributions from the conventional (collinear) ones. For the input, we have applied the recent LO NNPDF3.1 set [47]. The parton level calculations according to (4) were performed using the Monte-Carlo generator pegasus. Of course, we take into account the feed-down contributions from χ b (3P ), χ b (2P ), χ b (1P ), Υ(3S) and Υ(2S) decays. Numerically, everywhere we set the masses m Υ(1S) = 9.4603 GeV, m Υ(2S) = 10.02326 GeV, [48] and adopt the usual non-relativistic approximation m b = m Q /2 for the beauty quark mass, where m Q is the mass of bottomonium Q. We set the necessary branching ratios as they are given in [48]. Note that there is no experimental data for the branching ratios of χ b (3P ), so we use the results of assumption [22] that the total decay widths of χ b (mP ) are approximately independent on m. So, we have B(χ b1 (3P ) → Υ(1S) + γ) = 0.0381 and B(χ b2 (3P ) → Υ(1S) + γ) = 0.0192 [22]. We use the one-loop formula for the QCD coupling α s with n f = 4(5) quark flavours at Λ QCD = 250(167) MeV for A0 (KMR) gluon density and two-loop expression for α s with n f = 4 and Λ QCD = 200 MeV for JH'2013 set 1. We set the color-singlet NMEs 0 ] = 2.30 GeV 5 as obtained from the potential model calculations [49]. All the NMEs for Υ(2S), Υ(3S), χ b (2P ) and χ b (3P ) mesons were derived in [24,25].
To determine the NMEs for both Υ(1S) and χ b (1P ) mesons we have performed a global fit to the Υ(1S) production data at the LHC. We have included in the fitting procedure the Υ(1S) transverse momentum distributions measured by the CMS [26,27] and ATLAS [29] Collaborations at √ s = 7 and 13 TeV. Similar to the NRQCD analyses [21][22][23], we have excluded from our fit the low p T region and considered only data at p T > p cut T = 10 GeV. We note that at low transverse momenta a more accurate treatment of large logarithms ∼ ln m 2 Υ /p 2 T and other nonperturbative effects become necessary. To determine NMEs for χ b (1P ) mesons, we also included into the fit the recent LHCb data [32] on the radiative χ b (1P ) → Υ(1S) + γ decays collected at √ s = 7 and 8 TeV and the recent CMS [28] and LHCb data [33] on the ratio σ(χ b2 (1P ))/σ(χ b1 (1P )) collected at √ s = 8 TeV. Our analysis strategy is the following. First, we found that the p T shape of the direct 1 ] contributions to the Υ(1S) production is almost the same in all kinematical regions probed at the LHC. Thus, the ratio can be well approximated by a constant for a wide Υ(1S) transverse momentum p T and rapidity y range at different energies. For example, we estimate the mean-square average r = 1.743 ± 0.010 for the A0 set, which is practically the same for all other TMD gluon densities in a proton. So, we construct a linear combination which can be only extracted from the measured Υ(1S) transverse momentum distributions. Note that here we considered the color singlet wave functions of χ b (1P ) mesons as independent (not necessarily identical) free parameters, as it was proposed [50] to describe the LHC data on relative σ(χ c2 )/σ(χ c1 ) production rate. Of course, we understand that doing so is at odds with the Heavy Quark Effective Theory (HQET) and Heavy Quark Spin Symmetry (HQSS). However, it was argued [50] that the HQSS predictions must not be taken for granted 2 . Thus, here we try two alternative scenarios for the mesons. We assume HQSS violation either solely for the color singlet states ("fit A") or for both color singlet and color octet states ("fit B"). In the latter case, the color octet NMEs for χ b (1P ) mesons are also treated as independent parameters not related to each other through the (2J + 1) factor. Thus, we introduce the ratios: and obtain the mean-square average values r CO 0 = 0.01737 ± 0.00006, r CO 1 = 0.304 ± 0.002 and r CO 2 = 0.1626 ± 0.0010 (for the A0 gluon density). Then, instead of (6), we have a modified linear combination for the color octet NMEs: Next, we found that the p T shapes of the direct Υ[ 3 P 2 ] contributions to the Υ(1S) production are also the same in all kinematical regions. So, the ratios 1 ]/dp T can be approximated by constants for a wide Υ(1S) transverse momentum p T and rapidity y range at different energies. For example, we estimate the mean-square average r 1 = 0.91 ± 0.02 and r 2 = 104 ± 2 for the A0 set. Then we construct a linear combination which can be extracted from the measured Υ(1S) transverse momentum distributions. As 2 The possible reason may be seen in the spin-orbilal interactions or in radiative corrections which can be large (see more discussion [51]). the next step, we use the recent LHCb data [32] on the ratio of Υ(1S) mesons originating from the χ b (1P ) radiative decays measured at √ s = 7 and 8 TeV: In the "fit A" scenario, from the known M r , M r 1 r 2 and R 0 ] and the linear combination In the case of "fit B", we can determine 1 ] . Finally, we use recent CMS [28] and LHCb data [33] measured at √ s = 8 TeV on the ratio From the known 2 ] values for the first fit. For the second one we use only the CMS data [28], because the LHCb data [33] are very few and only increase the total error of the fitted quantities. So, from the known M CS , M CO and R 1 ] values. Therefore, we have reconstructed the full map of the NMEs for both Υ(1S) and χ b (1P ) mesons.
The fitting procedure described above was separately done in each of the rapidity subdivisions (using the fitting algorithm as implemented in the commonly used gnuplot package [52]) under the requirement that all the NMEs are strictly positive. Then, the meansquare average of the fitted values was taken. The corresponding uncertainties are estimated in the conventional way using Student's t-distribution at the confidence level P = 80%. The results of our fits are collected in Tables 1 and 2. For comparison, we also presented there the NMEs obtained in the conventional NLO NRQCD by other authors [21]. Note that the results [21] were obtained from the fit on the same data set as in our analysis. The corresponding χ 2 /d.o.f. are listed in Table 3, where we additionally show their dependence on the minimal Υ(1S) transverse momenta involved into the fit p cut T . As one can see, the χ 2 /d.o.f. tends to stay the same or slightly increase when p cut T grows up and the best fit of the LHC data is achieved with the A0 and KMR gluon, although other gluon densities also return reliable χ 2 /d.o.f. values. We note that including into the fit the latest CMS data [27] taken at √ s = 13 TeV leads to 2 3 times higher values of χ 2 /d.o.f., as it was with the data on Υ(2S) [25]. We have checked that this is true for both the k T -factorization and collinear approaches 3 and, therefore, it could be a sign of some inconsistency between these CMS data and other measurements.
Both fit scenarios result in unequal values for χ b1 (1P ) and χ b2 (1P ) color singlet wave functions. So, for "fit A" we achieved the ratio O χ b2 (1P ) [ 3 P 0 ] ∼ 2.6 : 4.8 : 1 for the JH'2013 set 1, ∼ 2.6 : 3.9 : 1 for the KMR and ∼ 1 : 3 : 1 for the A0 gluon densities, respectively. This is an obvious contradiction with naive expectations based on the number of spin degrees of freedom, ∼ 5 : 3 : 1. The difference between the predictions for this ratio obtained with the considered TMD gluon densities could be a sign of a sensitivity of the relative production rate R χ b2 (1P ) χ b1 (1P ) to the gluon distributions and/or due to lack of the experimental data. If we assume the HQSS violation in the color octet sector as well (the "fit B" scenario), the fitted values of the color singlet NMEs of the χ b1 (1P ) and χ b2 (1P ) mesons also differ from each other (see Table2). The latter qualitatively agrees with the observations [50,51] done in the case of the χ c mesons.
All the data used in the fits above are compared with our predictions in Figs. 1 4. The shaded areas represent the theoretical uncertainties of our calculations, which include the uncertainties coming from the NME fitting procedure and the scale uncertainties. To estimate the latter, the standard variations in default renormalization scale (which is set to be equal to µ 2 R = m 2 Q +p 2 T ), namely, µ R → 2µ R or µ R → µ R /2 were introduced with replacing the A0 and JH'2013 set 1 gluon densities by the A0+ and JH'2013 set 1+, or by the A0− and JH'2013 set 1− ones. This was done to preserve the intrinsic correspondence between the TMD gluon set and the factorization scale taken as µ 2 F =ŝ + Q 2 T (where Q T is the net transverse momentum of incoming off-shell gluon pair) according to the TMD gluon fits (see [44,45] for more information). Of course, in the case of KMR gluons both factorization and renormalization scales have been varied to estimate the scale uncertainties. One can see that we have achieved a reasonably good description of the CMS [26,27] and ATLAS [29] data for the Υ(1S) transverse momentum distributions in the whole p T range within the experimental and theoretical uncertainties. The relative production rates R χ b2 (1P ) χ b1 (1P ) measured by the CMS [28] and LHCb [33] Collaborations and the R χ b (1P ) Υ(1S) ratios measured by the LHCb Collaboration [32] at √ s = 7 and 8 TeV are also reproduced well. However, our predictions for the R Υ(1S) rates tend to overestimate a bit the LHCb data [32], although they are rather close to the measurements within the uncertainties bands (see Fig. 3). The same situation is observed in the conventional NRQCD scenario, where the NLO NRQCD calculations [21][22][23] also overestimate the experimental data for the R Evaluation of these observables involves the NMEs for Υ(2S), Υ(3S), χ b (2P ) and χ b (3P ) mesons determined previously [24,25]. Of course, both scenarios, "fit A" and "fit B", lead to exactly the same results for the Y (1S) transverse momentum distributions due to full correspondence between (6) and (10). The corresponding predictions differ to each other for the relative production rates R χ b (mP ) Υ(nS) and/or R χ b2 (1P ) χ b1 (1P ) only (see Figs. 3 and 4). In addition, we have checked our results with the data, not included into the fit procedure: namely, the rather old CDF data [53] taken at √ s = 1.8 TeV and the LHCb data [30,31] taken in the forward rapidity region 2 < y < 4.5 at √ s = 7, 8 and 13 TeV (see Fig. 5). As one can see, we acceptably describe all the data above. Moreover, we find that the KMR gluon density does a much better job here than the JH'2013 set 1 or A0 distributions. Remarkably, the KMR gluon is only one TMD gluon density which is able to reproduce well the measurements in the low p T region.
Based on that we have investigated the sensitivity of our fit to the low p T region using the KMR gluon density function. To do this we include the low p T region into our fit ("fit C"). Our results can be find in Table 2 and Figs. 6 10, where we used the KMR gluon density using the NMEs from the "fit A" and "fit C". One can see that both fit scenarios give overall almost the same results, and only for the ratios R gives slightly better results. However, the uncertainties for the "fit C" are larger than for the "fit A" scenario. Additionally, we present the corresponding χ 2 /d.o.f. using the NMEs from the "fit A" for A0, JH'2013 set 1, KMR distributions and from the "fit C" for KMR only, that are listed in Table 4. We should note that using the ATLAS [29] and CMS [26,27] data give us lower values of χ 2 /d.o.f. than using only the CMS data. This is due to large uncertainties of the ATLAS data in the low p T region. One can see that χ 2 /d.o.f. for the "fit C" are larger not only compared with the results using the NMEs from the "fit A" scenario in the low p T region, but also compared with NMEs from Table 1 in the p T > 10 GeV region. These results justify our exclusion of the low p T region from our fit. Now we turn to the polarization of Υ(1S) mesons at the LHC conditions. It is well known that the polarization of any vector meson can be described with three parameters λ θ , λ φ and λ θφ , which determine the spin density matrix of a meson decaying into a lepton pair and can be measured experimentally. The double differential angular distribution of the decay leptons can be written as [54]: where θ * and φ * are the polar and azimuthal angles of the decay lepton measured in the meson rest frame. The case of (λ θ , λ φ , λ θφ ) = (0, 0, 0) corresponds to an unpolarized state, while (λ θ , λ φ , λ θφ ) = (1, 0, 0) and (λ θ , λ φ , λ θφ ) = (−1, 0, 0) refer to fully transverse and fully longitudinal polarizations. The CMS Collaboration has measured all of these polarization parameters for the Υ(1S) mesons as functions of their transverse momentum in three complementary frames: the Collins-Soper, helicity and perpendicular helicity ones at √ s = 7 TeV [34]. The CDF Collaboration has measured the polarization parameters in the helicity frame at √ s = 1.96 TeV [35]. The frame-independent parameterλ = (λ θ + 3λ φ )/(1 − λ φ ) has been additionally studied. As it was done previously [24,25], to estimate λ θ , λ φ , λ θφ andλ we generally follow the experimental procedure. We collect the simulated events in the kinematical region defined by the experimental setup, generate the decay lepton angular distributions according to the production and decay matrix elements and then apply a three-parametric fit based on (16).
Our predictions are shown in Figs. 11 14. The calculations were performed using the A0 gluon density which provides the best description of the measured Υ(1S) transverse momenta distributions at the LHC conditions. The NMEs from Table1 (the "fit A" scenario) were applied. As one can see, we find only a weak or zero polarization in the all kinematic regions, that perfectly agrees with the CMS and CDF measurements. This agreement shows no fundamental problems in describing the Υ(1S) polarization data. Moreover, the calculated polarization parameters λ θ , λ φ , λ θφ andλ are stable with respect to variations in the model parameters. In fact, there is no dependence on the strong coupling constant and/or TMD gluon densities in a proton. As it was already pointed out above, our results for λ θ , λ φ , λ θφ andλ are based on the key assumption [20] that the intermediate color octet states are states with a definite total angular momentum J and its projection J z , rather than states with definite projections of a spin S z and orbital angular momentum L z . Given that, the transition amplitudes only involve the polarization vector associated with J z and not with L z . As a result, we have no conservation of S z in the electric dipole transitions. Under this assumption, we have achieved a reasonable simultaneous description for all of the available data for the Υ(1S) and χ b (1P ) mesons (the transverse momentum distributions, relative production rates and polarization observables). We have obtained earlier similar results for charmonia (J/ψ, ψ , χ c ), Υ(2S) and Υ(3S) polarizations [24,25,51,55,56]. Thus, keeping in mind the remarkable absence of tension with the η c production data (see [55,56]), one can conclude that the approach [20] results in the self-consistent and simultaneous description of charmonium and bottomonium data and therefore can be considered as providing an easy and natural solution to the long-standing quarkonia production and polarization puzzle.
phenomenology of the LHC processes and TMD parton densities.
1 ] /GeV 3 0.030 ± 0.018 0.009 ± 0.011 0.008 ± 0.009 0.007 Table 2: The NMEs for the Υ(1S) and χ b (1P ) mesons as determined from our fit at p cut T = 10 GeV (the "fit B" scenario) and only for the KMR density at p cut T = 0 GeV (the "fit C" scenario).   Table 4: The dependence of χ 2 /d.o.f. achieved in the fit procedure (using NMEs from the "fit A" and "fit C" scenarios) on the choice of p cut T for only the ATLAS data and the ATLAS and CMS data combined.  TeV (lower histograms, divided by 100) in the different rapidity regions. The notation of all histograms is the same as in Fig. 1. The experimental data are from CMS [26,27]. calculated as a function of the Υ(1S) transverse momentum calculated at √ s = 7 and 8 TeV using the NMEs from Tables 1 and 2 (for the A0 gluon density). The notation of all histograms is the same as in Fig. 1. The experimental data are from LHCb [32].  Tables 1 and 2 (for the A0 gluon density). The notation of all histograms is the same as in Fig. 1. The experimental data are from CMS [28] and LHCb [33]. TeV in the different rapidity regions. The notation of all histograms is the same as in Fig. 1. The experimental data are from CDF [53] and LHCb [30,31]. Figure 6: Transverse momentum distribution of the inclusive Υ(1S) production calculated at √ s = 7 TeV in the different rapidity regions. The green and purple histograms correspond to the "fit A" and "fit C" predictions obtained with the KMR gluon density. Shaded bands represent the NME uncertainties of our calculations, as it is described in the text. The experimental data are from ATLAS [29]. calculated as a function of the Υ(1S) transverse momentum calculated at √ s = 7 and 8 TeV. The notation of all histograms is the same as in Fig. 6. The experimental data are from LHCb [32]. TeV in the different rapidity regions. The notation of all histograms is the same as in Fig. 6. The experimental data are from CDF [53] and LHCb [30,31]. Figure 11: The polarization parameters λ θ , λ φ , λ θφ andλ of the Υ(1S) mesons calculated in the CS frame as functions of its transverse momentum at √ s = 7 TeV. The A0 gluon density is used. The blue and red histograms correspond to the predictions obtained at |y| < 0.6 and 0.6 < |y| < 1.2, respectively. The experimental data are from CMS [34]. Figure 12: The polarization parameters λ θ , λ φ , λ θφ andλ of the Υ(1S) mesons calculated in the helicity frame as functions of its transverse momentum at √ s = 7 TeV. The notation of all histograms is the same as in Fig. 11. The experimental data are from CMS [34]. Figure 13: The polarization parameters λ θ , λ φ , λ θφ andλ of the Υ(1S) mesons calculated in the perpendicular helicity frame as functions of its transverse momentum at √ s = 7 TeV. The notation of all histograms is the same as in Fig. 11. The experimental data are from CMS [34]. Figure 14: The polarization parameters λ θ andλ of the Υ(1S) mesons calculated in the helicity frame as functions of its transverse momentum at √ s = 1.96 TeV. The notation of all histograms is the same as in Fig. 11. The experimental data are from CDF [35].