  Numerical Stochastic Perturbation Theory (NSPT) allows for perturbative
computations in quantum field theory. We present an implementation of NSPT that
yields results for high orders in the perturbative expansion of lattice gauge
theories coupled to fermions. The zero-momentum mode is removed by imposing
twisted boundary conditions; in turn, twisted boundary conditions require us to
introduce a smell degree of freedom in order to include fermions in the
fundamental representation. As a first application, we compute the critical
mass of two flavours of Wilson fermions up to order $O(\beta^{-7})$ in a
${\mathrm{SU}}(3)$ gauge theory. We also implement, for the first time,
staggered fermions in NSPT. The residual chiral symmetry of staggered fermions
protects the theory from an additive mass renormalisation. We compute the
perturbative expansion of the plaquette with two flavours of massless staggered
fermions up to order $O(\beta^{-35})$ in a ${\mathrm{SU}}(3)$ gauge theory, and
investigate the renormalon behaviour of such series. We are able to subtract
the power divergence in the Operator Product Expansion (OPE) for the plaquette
and estimate the gluon condensate in massless QCD. Our results confirm that
NSPT provides a viable way to probe systematically the asymptotic behaviour of
perturbative series in QCD and, eventually, gauge theories with fermions in
higher representations.
