  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.
