Abstract. We present the calculations of associated production of Z bosons and heavy (charm or beauty) quarks at the LHC energies in the framework of k T -factorization QCD approach. Our consideration is mainly based on the O(αα 2 s ) off-shell gluon-gluon fusion subprocess g * g * → ZQQ, where produced Z boson subsequently decays into a lepton pair. Several subleading contributions from the O(αα 2 s ) and O(αα 3 s ) subprocesses are taken into account. Contributions from the double parton scattering mechanism are estimated. The transverse momentum dependent (or unintegrated) gluon densities in a proton are determined using the Catani-Ciafaloni-Fiorani-Marchesini evolution equation. We achieve reasonably good agreement of our predictions and latest experimental data taken by the CMS and ATLAS Collaborations, discuss the theoretical uncertainties of our calculations and demonstrate the importance of subleading quark contributions in description of the LHC data in the whole kinematical region.Very recently, the CMS Collaboration has presented the measurements of Z boson production in association with charm [1] quarks at the LHC energy √ s = 8 TeV. Also measurements of angular correlations of B-hadrons at √ s = 7 TeV were presented [2] and the ATLAS Collaboration has reported the experimental data on the total and differential cross sections of Z boson production in association with beauty quark jets at √ s = 7 TeV [3]. Such processes involve both strong and weak interactions, so are important as global test of the Standard Model (SM). The Z + b-jets production is an important background for studies of the associated production of Higgs and Z bosons, where the Higgs boson decays into bb pairs [4][5][6]. Many physics scenarios beyond the Standard Model (SM), for example, new generations of heavy quarks (b ′ , t ′ ) decaying into Z bosons and b-quarks [7], supersymmetric Higgs bosons produced in association with beauty quarks [8] and some SM extentions with additional S U(2) doublets with enhanced Zbb coupling [9], predict final states with b-quarks and Z bosons. In addition, Z + b/Z + c cross sections ratio is highly sensitive to the charm content of the proton [10]. Finally, such processes may serve as potential indicators of the Double Parton Scattering (DPS) mechanism [11][12][13].In the present work we analyse recent CMS [1, 2] and ATLAS [3] data using so called k Tfactorization QCD approach [14]. This approach is based on the famous Balitsky-Fadin-KuraevLipatov (BFKL) or Ciafaloni-Catani-Fiorani-Marchesini (CCFM) gluon evolution equations and pro-⋆