  We study in detail the impact of dynamical quarks on the gluon mass
generation mechanism, in the Landau gauge, for the case of a small number of
quark families. As in earlier considerations, we assume that the main bulk of
the unquenching corrections to the gluon propagator originates from the fully
dressed quark-loop diagram. The nonperturbative evaluation of this diagram
provides the key relation that expresses the unquenched gluon propagator as a
deviation from its quenched counterpart. This relation is subsequently coupled
to the integral equation that controls the momentum evolution of the effective
gluon mass, which contains a single adjustable parameter; this constitutes a
major improvement compared to the analysis presented in Phys. Rev. D86 (2012)
014032, where the behaviour of the gluon propagator in the deep infrared was
estimated through numerical extrapolation. The resulting nonlinear system is
then treated numerically, yielding unique solutions for the modified gluon mass
and the quenched gluon propagator, which fully confirm the picture put forth
recently in several continuum and lattice studies. In particular, an infrared
finite gluon propagator emerges, whose saturation point is considerably
suppressed, due to a corresponding increase in the value of the gluon mass.
This characteristic feature becomes more pronounced as the number of active
quark families increases, and can be deduced from the infrared structure of the
kernel entering in the gluon mass equation.
