Associated production of $Z$ bosons and $b$-jets at the LHC in the combined $k_T$ + collinear QCD factorization approach

We consider the production of $Z$ bosons associated with beauty quarks at the LHC using a combined $k_T$ + collinear QCD factorization approach, that interpolates between small $x$ and large $x$ physics. Our consideration is based on the off-shell gluon-gluon fusion subprocess $g^* g^* \to Z Q\bar Q$ at the leading order ${\cal O}(\alpha\alpha_s^2)$ (where the $Z$ boson further decays into a lepton pair), calculated in the $k_T$-factorization approach, and several subleading ${\cal O}(\alpha \alpha_s^2 )$ and ${\cal O}(\alpha \alpha_s^3 )$ subprocesses involving quark-antiquark and quark-gluon interactions, taken into account in conventional (collinear) QCD factorization. The contributions from double parton scattering are discussed as well. The transverse momentum dependent (or unintegrated) gluon densities in a proton are derived from Catani-Ciafaloni-Fiorani-Marchesini (CCFM) evolution equation. We achieve reasonably good agreement with the latest data taken by CMS and ATLAS Collaborations. The comparison of our results with next-to-leading-order pQCD predictions, obtained in the collinear QCD factorization, is presented. We discuss the uncertainties of our calculations and demonstrate the importance of subleading quark involving contributions in describing the LHC data in the whole kinematic region.


Introduction
With the LHC in operation, one can access a number of "rare" processes which could have never been systematically studied at previous accelerators. In this article we draw attention to the associated production of Z bosons and b-jets. This process involves both weak and strong interactions and therefore serves as a complex test of Standard Model, perturbative QCD (pQCD) and our knowledge of parton distribution functions in a proton. Similarly to the W + c and W + b processes considered earlier [1,2] it probably provides an arena for double parton scattering (DPS), now widely discussed in the literature. We wish to clarify this point in our paper. Besides that, this process constitutes a substantial background in studying the associated production of Higgs and Z bosons, where the Higgs boson is identified via its decay into a bb pair [3][4][5]. A number of physics scenarios beyond Standard Model also refer to final states containing Z bosons and beauty quarks [6][7][8], thus making the related studies important and topical.
Our present study is greatly stimulated by the recent ATLAS measurements [9] of the total and differential production cross sections of Z bosons associated with beauty quark jets at √ s = 7 TeV accompanied by the CMS measurements [10] of kinematic correlations between Z bosons and b-hadrons at √ s = 7 TeV. We investigate these processes in the framework of a combined QCD approach, based on the k T -factorization formalism [11,12] in the small-x domain and conventional (collinear) QCD factorization at large Bjorken x. Doing so, we employ the k T -factorization approach to calculate the leading contributions from the off-shell gluon-gluon fusion g * g * → ZQQ and, to extend the consideration to the whole kinematic range, take into account several subleading quark-involved subprocesses using collinear QCD factorization. The k T -factorization approach has certain technical advantages in the ease of including higher-order radiative corrections that can be taken into account in the form of transverse momentum dependent (TMD) parton distributions 1 . This approach has become a widely exploited tool and it is of interest and importance to test it in as many cases as possible. Closely related to this is the selection of TMD parton densities best suited to describe the data. These tasks form the major goal of our article. The outline of the paper is the following. In Section 2 we describe our approach and parameter setting. In Section 3 we present the results of our calculations and confront them with the available data. Our conclusions are summarised in Section 4.

The model
Let us start from a short review of calculation steps. The leading contribution comes from the O(αα 2 s ) off-shell gluon-gluon fusion subprocess: where the four-momenta of all particles are given in the parentheses. The corresponding gauge-invariant off-shell amplitude was calculated earlier [14,15] and implemented into the Monte-Carlo event generator cascade [16]. All the details of these calculations have been explained [14,15], we only mention here that the standard QCD Feynman rules were employed with the only exception that the initial off-shell gluon spin density matrix was determined according to the k T -factorization prescription [11,12]: with k i T being the component of the gluon momentum k i (with i = 1 or 2) perpendicular to the beam axis (k 2 i = −k 2 iT = 0). In the collinear limit k 2 iT → 0 this expression converges to the ordinary one after averaging on the azimuthal angle.
In order to fully reproduce the experimental setup [9,10], we simulate the subsequent decay Z → l + l − incorporated with the production step at the amplitude level. Then, the Z boson propagator is parametrised in Breit-Wigner form with mass m Z = 91.1876 GeV and total decay width Γ Z = 2.4952 GeV [17]. The role of virtual photons in the Z boson resonance region is found to be small: it makes not more than a 2% or 3% correction (including the Z/γ * interference effects). This is much less than the scale uncertainty of the main subprocess (see Section 3), and, therefore, is neglected in our analysis.
In addition to off-shell gluon-gluon fusion, we take into account several subprocesses involving quarks in the initial state. These are the flavor excitation at O(αα 2 s ): the quark-antiquark annihilation at O(αα 2 s ): and the quark-gluon scattering at O(αα 3 s ): Quark densities are typically much lower than the gluon density at the LHC conditions; however, these processes may become important at very large transverse momenta (or, respectively, at large parton longitudinal momentum fraction x, which is needed to produce large p T events) where the quarks are less suppressed or can even dominate over the gluon density. Here we find it reasonable to rely upon collinear Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) factorization scheme [18], which provides better theoretical grounds in the large-x region. So, we consider a combination of two techniques with each of them being used at the kinematic conditions where it is best suitable (gluoninduced subprocess (1) at small x and quark-induced subprocesses (3) -(5) at large x values). For the flavor excitation and the quark-antiquark annihilation we apply the on-shell limit of formulas obtained earlier [19] supplementing them by the Z boson decays. The amplitude of quark-gluon scattering subprocess can be easily derived from the gluon-gluon fusion one. As usual, to calculate the contributions of quark-induced subprocesses (3) -(5) one has to convolute the corresponding partonic cross sections dσ ab with the conventional parton distribution functions f a (x, µ 2 ) in a proton: where indices a and b denote quark and/or gluon, x 1 and x 2 are the fractions of longitudinal momenta of colliding protons and µ 2 is the hard scale. In the case of off-shell gluon-gluon fusion we employ the k T -factorization formula: where f g (x, k 2 T , µ 2 ) is the TMD gluon density in a proton. To obtain the latter we use a numerical solution of the CCFM equation [20]. It provides a suitable tool as it smoothly interpolates between the small-x Balitsky-Fadin-Kuraev-Lipatov (BFKL) [21] gluon dynamics and large-x DGLAP one. We adopt the latest JH'2013 parametrization [22], taking JH'2013 set 2 as the default choice. The corresponding TMD gluon density has been fitted to high-precision DIS data on the proton structure functions F 2 (x, Q 2 ) and F c 2 (x, Q 2 ). The fit was based on TMD matrix elements and involves two-loop strong coupling constant, kinematic consistency constraint [23,24] and non-singular terms in the CCFM gluon splitting function [25]. For the conventional quark and gluon densities we use the MSTW'2008 (LO) set [26].
Throughout this paper, all calculations are based on the following parameter setting. In the collinear QCD factorization case we use one-loop running strong and electroweak coupling constants with n f = 4 massless quark flavors and Λ QCD = 200 MeV; the factorization and renormalization scales are both set equal to the Z boson transverse mass, so that we have α s (m 2 Z ) = 0.1232 and α(m 2 Z ) = 1/128. In the k T -factorization case we use a two-loop expression for the strong coupling constant (as it was originally done in the fit [22]) and define the factorization scale as µ 2 F =ŝ + Q 2 T withŝ and Q 2 T being the subprocess invariant energy and the net transverse momentum of the initial off-shell gluon pair, respectively. The latter definition of µ F is unusual and is dictated by the CCFM evolution algorithm [22]. The b-quark mass and Weinberg mixing angle were set to m b = 4.75 GeV and sin 2 θ W = 0.2312 [17]. When necessary, b-quarks were converted into b-hadrons using Peterson fragmentation function [27] with ǫ b = 0.006.
We close our consideration with DPS contributions where we apply a simple factorization formula (for details see the reviews [28][29][30] and references therein): where σ eff is a normalization constant which incorporates all "DPS unknowns" into a single phenomenological parameter. A numerical value of σ eff ≃ 15 mb was earlier obtained from fits to pp and pp data. This will be taken as the default value throughout the paper. The calculation of inclusive cross sections σ(b +b) and σ(Z) is straightforward and needs no special explanations. Here we strictly follow the approach described earlier [31][32][33]. The multidimensional phase space integration was performed by means of the Monte Carlo technique, using the routine vegas [34]. In the next section we confront our predictions with the latest LHC data.

Numerical results
This section presents a detailed comparison between theoretical calculations and recent LHC data. The essential measurements have been carried out by the ATLAS [9] and CMS [10] Collaborations and refer to the following categories: Z bosons produced in association with one beauty jet, Z bosons produced in association with two beauty jets and Z bosons produced in association with explicitly reconstructed b-hadrons. In addition to the above, the ATLAS Collaboration has presented [9] inclusive cross sections for Z bosons associated with any number of b-jets. We do not analyse events of this kind in the present study and only concentrate on the production of Z bosons with one or two b-jets.

Production of Z bosons in association with one b-jet
The ATLAS Collaboration has collected the data [9] at √ s = 7 TeV. Both leptons originating from the Z boson decay are required to have p l T > 20 GeV and |η l | < 2.4, the lepton pair invariant mass lies in the interval 76 < M ll < 106 GeV, the beauty jets are required to have p b T > 20 GeV and |η b | < 2.4. We confront our predictions with the available data in Figs. 1 and 2. To estimate the theoretical uncertainties in the quark-involving subprocesses (3) -(5), calculated using the collinear QCD factorization, we have varied the scales µ R and µ F by a factor of 2 around their default values. In the k T -factorization approach, employed for off-shell gluon-gluon fusion subprocess (1), the scale uncertainties were estimated by using the gluon densities JH'2013 set 2+ and JH'2013 set 2− instead of default density JH'2013 set 2. These two sets refer to the varied hard scales in the strong coupling constant α s in the off-shell amplitude: JH'2013 set 2+ stands for 2µ R , while JH'2013 set 2− refers to µ R /2 (see [22] for more information). The estimated scale uncertainties are shown as shaded bands. As one can see, we achieve reasonably good agreement with the ATLAS data [9] within the experimental and theoretical uncertainties, although we observe some underestimation of these data at high p Z T and a slight overestimation at small transverse momenta. The slight overestimation of the data at low p Z T can probably be attributed to the TMD gluon density used, since the region p Z T < 100 GeV is fully dominated by off-shell gluon-gluon fusion, as it is demonstrated in Fig. 2. The rapidity distribution is well described practically everywhere. The NLO pQCD calculations 2 , performed using mcfm routine [35], tend to slightly overestimate our predictions and better decribe the data at large transverse momenta.
To investigate the importance of k T -factorization, we have repeated the calculation using collinear QCD factorization for all considered subprocesses (dash-dotted histograms in Fig. 1). We find that these effects are significant at low and moderate p Z T (up to p Z T ∼ 100 GeV), where the off-shell gluon-gluon fusion dominates. The effect of using k Tfactorization for gluon-dominated processes is clearly demonstrated in Fig. 2. The quarkinitiated subprocesses (3) -(5) become important only at high transverse momenta, where the typical x values are large, and that supports using of the DGLAP quark and gluon dynamics for these subprocesses (see Fig. 2). The subprocesses (3) -(5) are important to achieve an adequate description of the data in the whole p Z T region. The estimated DPS contributions are found to be small in the considered kinematic region. Some reasonable variations in σ eff ≃ 15 ± 5 mb would affect DPS predictions, though without changing our basic conclusion. We note also that scale uncertainties of the CCFM-based predictions are comparable with the ones of NLO pQCD calculations.

Production of Z bosons in association with two b-jets
The data provided by the ATLAS [9] Collaboration refer to the same energies and kinematic restrictions as in the previous subsection. The observables shown by the ATLAS Collaboration are the Z boson transverse momentum p Z T and rapidity y Z , invariant mass of the b-jet pair M bb and angular separation in η − φ plane between the jets ∆R bb . The latter is useful to identify the contributions where scattering amplitudes are dominated by terms involving gluon splitting g → Q +Q.
The results of our calculations are shown in Fig. 3 in comparison with the ATLAS data [9]. As one can see, our results describe the data reasonably well within the experimental and theoretical uncertainties, although some tendency to slightly underestimate the data at high transverse momentum p Z T and large M bb can be seen. The role of offshell gluon-gluon fusion subprocess is a bit enhanced here compared to the case of Z + b production because the quark-antiquark annihilation subprocess (4) gives a negligible contribution and gluon splitting subprocess (5) populates mainly at low η − φ distances ∆R bb . This subprocess is complementary to the one [36] where quark-gluon scattering q + g * was dominant. The estimated DPS contribution is small and can play a role at low p Z T only. The NLO pQCD calculations, performed using mcfm program 3 , tend to slightly underestimate the ATLAS data at low ∆R bb and M bb , although provide better description of the data at large transverse momentum p Z T and invariant mass M bb .

Production of Z bosons in association with two b-hadrons
In the measurements reported by CMS Collaboration [10], both b-hadrons have been identified explicitly by their full decay reconstruction. This data sample allows to study the production properties of a Zbb system even in the region of small angular seperation between the b quarks (where the usual jet analysis is not possible as the jets would overlap). In a specific subsample, an additional cut on the Z boson transverse momentum is applied, p Z T > 50 GeV. The CMS Collaboration described the angular configuration of the Zbb system in terms of spatial (in η − φ plane) and azimuthal separation between the b-hadrons ∆R bb and ∆φ bb , spatial separation min ∆R Zb between the Z boson and closest b-hadron and the asymmetry in the Zbb system defined as where max ∆R Zb is the distance between the Z boson and remote b-hadron. The correlation observables are useful to identify the different production mechanisms (or specific higher-order corrections). For example, low min ∆R Zb identifies Z bosons in the vicinity of one of the b-hadron (Z bosons promptly radiated from b-quarks), small ∆φ bb indicates gluon to quark splitting g → Q+Q. Moreover, while the configurations where the two b-hadrons are emitted symmetrically with respect to the Z directions leads to a zero value of A Zbb assymetry, the additional final-state gluon radiation results in a non-zero one, that provides us with the possibility to test the high-order pQCD corrections. Our predictions are shown in Figs. 4 and 5 in comparison with the CMS data [10]. As one can see, our results with default b-quark fragmentation parameters reasonably well describe the data within the theoretical and experimental uncertainties. To estimate an additional uncertainty coming from the b-quark fragmentation, we repeated our calculations with varied shape parameter ǫ b = 0.003 (not shown), which is often used in NLO pQCD calculations. We find that the predicted cross sections (in the considered p T region) are larger for smaller ǫ b values. However, the typical dependence of numerical predictions on the fragmentation scheme is much smaller than the scale uncertainties of our calculations. The NLO pQCD predictions, obtained using the amc@nlo [37] event generator 4 , are rather close to our results.

Conclusions
We have considered the associated Z boson and beauty quark production at the LHC conditions. The calculations were done in a "combined" scheme employing both the k Tfactorization and collinear factorization in QCD, with each of them used in the kinematic conditions of its best reliability. The dominant contribution is represented by the gluongluon fusion subprocess g * g * → Zbb with Z boson further decaying into a lepton pair. This subprocess is entirely (for the first time) calculated in the k T -factorization approach. A number of subleading subprocesses contributing at O(αα 2 s ) and O(αα 3 s ) have been considered in the conventional collinear scheme.
Using the TMD gluon densities derived from the CCFM evolution equation, we have achieved reasonably good agreement between our theoretical predictions and latest CMS and ATLAS experimental data collected at √ s = 7. We find that the (formally subleading) quark-involving subprocesses become especially important at high transverse momenta and are necessary to describe the data in the whole kinematic range. Our estimations of the double parton scattering show that the latter is unimportant. This conclusion is also confirmed by the fact that our single parton scattering calculations show no room for additional contributions when compared to the ATLAS and CMS data. The dash-dotted histograms correspond to the collinear limit of our calculations. The estimated DPS contributions and mcfm [35] predictions (taken from [9]) are shown additionally. The data are from ATLAS [9].  calculated as a function of the Z boson transverse momentum, rapidity, invariant mass of the b-jet pair and angular separation between the jets. Notation of the histograms is the same as in Fig. 1. The data are from ATLAS [9]. The mcfm [35] predictions are taken from [9]. Figure 4: Associated production of a Z boson and two b-hadrons at √ s = 7 TeV. Notation of the histograms is the same as in Fig. 1. The data are from CMS [10]. The amc@nlo [37] predictions are taken from [10]. additional kinematical cut on the Z boson transverse momentum p Z T > 50 GeV. Notation of the histograms is the same as in Fig. 1. The data are from CMS [10]. The amc@nlo [37] predictions are taken from [10].