  We present a novel analysis of the gluon gap equation, where its full
nonlinear structure is duly taken into account. In particular, while in
previous treatments the linearization of this homogeneous integral equation
introduced an indeterminacy in the scale of the corresponding mass, the current
approach determines it uniquely, once the value of the gauge coupling at a
given renormalization point is used as input. A crucial ingredient for this
construction is the "kinetic term" of the gluon propagator, whose form is not
obtained from the complicated equation governing its evolution, but is rather
approximated by suitable initial {\it Ans\"atze}, which are subsequently
improved by means of a systematic iterative procedure. The multiplicative
renormalization of the central equation is carried out following an approximate
method, which is extensively employed in the studies of the standard quark gap
equation. This approach amounts to the effective substitution of the vertex
renormalization constants by kinematically simplified form factors of the
three- and four-gluon vertices. The resulting numerical interplay, exemplified
by the infrared suppression of the three-gluon vertex and the mild enhancement
of the four-gluon vertex, is instrumental for obtaining positive-definite and
monotonically decreasing running gluon masses. The resulting gluon propagators,
put together from the gluon masses and kinetic terms obtained with this method,
match rather accurately the data obtained from large-volume lattice
simulations.
