  In Refs. [1,2] we determined the infinite volume coefficients of the
perturbative expansions of the self-energies of static sources in the
fundamental and adjoint representations in SU(3) gluodynamics to order
$\alpha^{20}$. We used numerical stochastic perturbation theory [3], where we
employed a new second order integrator and twisted boundary conditions. The
expansions were obtained in lattice regularization with the Wilson action and
two different discretizations of the covariant time derivative within the
Polyakov loop. Overall, we obtained four different perturbative series. For all
of them the high order coefficients displayed the factorial growth predicted by
the conjectured renormalon picture, based on the operator product expansion.
This enabled us to determine the normalization constants of the leading
infrared renormalons of heavy quark and heavy gluino pole masses. Here we
present improved determinations of the normalization constants and the
perturbative coefficients by incorporating the four-loop beta-function
coefficient (which we also determine) in the fit function.
