  We consider the production of $Z$ bosons associated with heavy (charm and
beauty) jets at the LHC energies using two scenarios based on the transverse
momentum dependent (TMD) parton densities in a proton. The first of them
employs the Catani-Ciafaloni-Fiorani-Marchesini gluon evolution and is
implemented in the Monte-Carlo event generator PEGASUS. Here, the heavy quarks
are always produced in the hard partonic scattering. The second scheme is based
on the parton branching approach, currently implemented into the Monte-Carlo
event generator CASCADE. In this scenario, the $Z$ + jets sample is generated
and then events containing the heavy flavor jet in a final state are selected.
We compare the predictions obtained within these two TMD-based approaches to
each other, investigate their sensitivity to the TMD gluon densities in a
proton and estimate the effects coming from parton showers and double parton
scattering mechanism. Additionally, we compare our predictions with the results
of traditional (collinear) pQCD calculations performed at NLO accuracy. It is
shown that the TMD-based results agree with the LHC experimental data collected
at $\sqrt s = 8$ and $13$ TeV. We discuss the sensitivity of observables to the
quark distributions in a proton and present predictions to search for the
intrinsic charm signal in forthcoming analyses of the LHC experimental data.
