  In this paper we present the perturbative computation of the difference
between the renormalization factors of flavor singlet
($\sum_f\bar\psi_f\Gamma\psi_f$, $f$: flavor index) and nonsinglet
($\bar\psi_{f_1} \Gamma \psi_{f_2}, f_1 \neq f_2$) bilinear quark operators
(where $\Gamma =
\mathbb{1},\,\gamma_5,\,\gamma_{\mu},\,\gamma_5\,\gamma_{\mu},\,
\gamma_5\,\sigma_{\mu\,\nu}$) on the lattice. The computation is performed to
two loops and to lowest order in the lattice spacing, using Symanzik improved
gluons and staggered fermions with twice stout-smeared links. The stout
smearing procedure is also applied to the definition of bilinear operators. A
significant part of this work is the development of a method for treating some
new peculiar divergent integrals stemming from the staggered formalism. Our
results can be combined with precise simulation results for the renormalization
factors of the nonsinglet operators, in order to obtain an estimate of the
renormalization factors for the singlet operators. The results have been
published in Physical Review D.
