Renormalization of quark bilinears I'm looking at the one-loop corrections to the amputated quark two-point functions ($\Gamma_i$) with insertions of quark bilinears (indexed by $i\in\{S,P,V,A,T\}$) with off-shell legs in Euclidean QCD. I approach the renormalization through an OPE, where the renormalized operators are defined through $O_i^B=Z_{ij}O_j^R$. Since the bilinears are dimension 3 (or $d-1$ in dim. reg.), they should only mix with themselves and $m^3$ through vacuum correlation functions in the case of the scalar. Since the two-point correlation function of $m^3$ is disconnected, I expect that this mixing does not appear in the current calculation. Then, for example, the bare, dimensionally regularized  one-loop correction to the scalar should look like
\begin{equation}
\Gamma_S^B=\Gamma_S^{(0)}\left\{1+g_0^2\frac{C_F}{(4\pi)^2}\left[\frac{A}{\epsilon}+B\right]+\mathcal{O}(g_0^4)\right\}
\end{equation}
for some constants $A$ and $B$ ($B$ includes the logarithms). If I retain all soft scales however, there is an extra term that appears:
\begin{equation}
\Gamma_S^B=\Gamma_S^{(0)}\left\{1+g_0^2\frac{C_F}{(4\pi)^2}\left[\frac{A}{\epsilon}+B\right]+\mathcal{O}(g_0^4)\right\}+C g_0^2\frac{C_F}{(4\pi)^2}\frac{i m_0 p\!\!/}{p^2}.
\end{equation}
This appears to be an off-diagonal mixing with some off-shell operator, but I cannot figure out what this operator could be. Since I'm working in dim. reg., the operator should be dimension 3 as well, but the inverse power of the momentum $p$ makes me nervous. How can an operator like this exist? Initially, I had assumed that my disconnectedness argument above was incorrect and that this was the amputated tree-level of the operator $m^3$, since that would in principle look like $m^3S^{-1}=m^3(i p\!\!/+m)$, where $S$ is the propagator. The new term can be recast as
\begin{equation}
\frac{i m_0 p\!\!/}{p^2}=(\alpha-1)\left(\frac{1}{m_0}S^{-1}-1\right),
\end{equation}
where $\alpha=1+\frac{m_0^2}{p^2}$ is a dimensionless quantity that already appears in $A$, $B$, and $C$. Since $C$ contains no UV poles, and there will therefore be no scale dependence in the residue, this seems safe (there would be a shift in $B$ by $-C(\alpha-1)$). Unfortunately, the $1/m_0$ above is equally unsatisfactory, and moreover this does not even have the desired form for the $m^3$ operator. Further, there is an equivalent structure for each other bilinear, including in particular the parity-odd bilinears which cannot mix with vacuum operators.
I'm not overly worried, since the extra term certainly does not renormalize, so extraction of the $Z$ factors should go smoothly, but I can imagine a situation where the new operator is related to this new term by a quantity proportional to the identity (as in the previous equation), which will shift the finite bit of the correlation function.
In the end, I'm back to asking what the structure $\frac{i m_0 p\!\!/}{p^2}$ could even be. Any ideas? Thank you all in advance.
 A: I think just to clear the air, let's be concrete. Let's calculate the amputated diagram in $d$-dimensions:

where the $\otimes$ is the $\overline{\psi} \psi$ insertion. The diagram is:
$$i \mathcal{M} = \mu^{4-d} \int \frac{d^d k}{(2 \pi)^d}(i g \gamma^\mu T^a) \frac{i}{ \displaystyle{\not} k \, + \displaystyle{\not} p \, - m} \frac{i}{\displaystyle{\not} k \, + \displaystyle{\not} p \, - m} (i g \gamma_\mu T^a) \frac{-i}{k^2}$$
$$ = -i g^2 \mu^{4-d} C_F \int \frac{d^d k}{(2 \pi)^d} \frac{d(k+p)^2 + 2 m (2-d) (\displaystyle{\not} k \, +\displaystyle{\not} p \,) + dm^2 }{k^2 ((k+p)^2 - m^2)^2}$$
The part of interest is the nontrivial $\gamma$-structure part (the middle term in the numerator). We can write this piece as:
$$\begin{split}
i\mathcal{M}_{\mathrm{traceless}} &= -2 (2-d) ig^2 \mu^{4-d} C_F m \int \frac{d^d k}{(2 \pi)^d} \frac{\displaystyle{\not} k \, +\displaystyle{\not} p \, }{k^2 ((k+p)^2 - m^2)^2} \\
&= -2 (2-d) ig^2 \mu^{4-d} C_F m \displaystyle{\not} p \, \left( \frac{1}{2 p^2} I(1,1) + \frac{1 + \frac{m^2}{p^2}}{2} I(1,2) - \frac{1}{2 p^2} I(0,2) \right) \\
&= \displaystyle{\not} p \, \left( \frac{C_F \alpha_S}{2 \pi m} + \frac{C_F \alpha_S \epsilon}{4 \pi m}\ln \frac{\tilde{\mu}^2}{m^2} + O(\epsilon^2)\right)
\end{split}$$
In the OPE picture, this seems to correspond to mixing with the dimension-4 operator $\overline{\psi} \displaystyle{\not} \partial \, \psi $, which is allowed due to the presence of the mass scale $m$. At this point, everything seems consistent to me, so maybe I am misunderstanding the question

Details:
The master integral is defined:
$$I(m,n) = \int \frac{d^d k}{(2 \pi)^d} \frac{1}{(k^2)^m ((k+p)^2 - m^2)^n}$$
