The prescription is actually more or less the same as if one uses dimensional regularization, based on the observation that the power divergences can be absorbed by counterterms at all scales is made. In the example you give, $a_n$ is a polynomial of external momenta, such as $p^2$ where $p$ is an external momentum, then just adding a counterterm $p^2\frac{\Lambda}{m}$ is capable of absorbing this divergence at all scales.
Now let me show a general statement: UV divergences from loop integrals will always be proportional to some polynomials of the external momenta. The proof can be found in Weinberg. To demonstrate it, let us see the following concrete example:
\begin{equation}
I(p)=\int_0^\infty\frac{k^mdk}{(k+p)^n}
\end{equation}
For this integral to have power divergences, we require that $m>n$. If we differentiate $I(p)$ $m+2-n$ times with respect to $p$, we get
\begin{equation}
\frac{d^{m+2-n}I(p)}{dp^{m+2-n}}\propto\int_0^\infty\frac{k^mdk}{(k+p)^{m+2}}
\end{equation}
which is UV finite, and by dimensional analysis, it must be proportional to $\frac{1}{p}$. Now if we integrate the above equation over $p$ for $m-n+2$ times, we get
\begin{equation}
I(p)=\sum_{a=0}^{m-n+1}c_ap^a\Lambda^{m-n-1-a}+d\cdot p^{m-n-1}\ln\frac{\Lambda}{p}
\end{equation}
where $c_a$, $d$ and $\Lambda$ are constants. As we can see in this example, the divergences are always multiplied by polynomials of the external momentum.
The general lesson to take from this example is that, after taking enough derivatives of the external momenta, the loop integrals will be finite and a Laurent polynomial of external momenta, and when we integrate over the external momenta to get the original integrals, we find the divergences are always proportional to some polynomials of external momenta.
Having this in mind, we can always choose local counterterms to cancel power divergences at all scales thus we do not need to worry about them. But logarithmic divergences can only be cancelled at a particular scale. It is this fact that gives the running of the parameters or amplitudes.
In $\overline{MS}$ the $\beta$ function does not depend on the gauge parameter. This means that the dependence on $\mu$ in $Z$ only comes from the coupling constant $a \propto \mu^{-2 \epsilon}$.
For general $Z(a_s,\xi)$, the relation is as follows (eq. 21 in the Chetyrkin paper):
$$-\gamma =\left(-\epsilon + \beta(a_s) \right) a_s \frac{\partial \log Z}{\partial a_s}+\gamma_3(a_s,\xi)\xi \frac{\partial \log Z}{\partial \xi}$$
Best Answer
The anomalous dimension for the field strength is defined as (eqn 12.63 Peskin):
$\gamma = \frac{1}{2} \frac{M}{Z} \frac{\partial Z}{\partial M} = \frac{1}{2} \frac{\partial \log Z}{\partial \log M} $.
This definition always holds. What you actually calculate for the right-hand side of the above equation once you have a Z within a particular scheme will be in general scheme-dependent.
Sorry, I can't help you with the $O(N)$ vector model...