  We calculate the quarkonic $O(\alpha_s^2)$ massive operator matrix elements
$\Delta A_{Qg}(N), \Delta A_{Qq}^{\rm PS}(N)$ and $\Delta A_{qq,Q}^{\rm NS}(N)$
for the twist--2 operators and the associated heavy flavor Wilson coefficients
in polarized deeply inelastic scattering in the region $Q^2 \gg m^2$ to
$O(\varepsilon)$ in the case of the inclusive heavy flavor contributions. The
evaluation is performed in Mellin space, without applying the
integration-by-parts method. The result is given in terms of harmonic sums.
This leads to a significant compactification of the operator matrix elements
and massive Wilson coefficients in the region $Q^2 \gg m^2$ derived previously
in \cite{BUZA2}, which we partly confirm, and also partly correct. The results
allow to determine the heavy flavor Wilson coefficients for $g_1(x,Q^2)$ to
$O(\alpha_s^2)$ for all but the power suppressed terms $\propto (m^2/Q^2)^k, k
\geq 1$. The results in momentum fraction $z$-space are also presented. We also
discuss the small $x$ effects in the polarized case. Numerical results are
presented. We also compute the gluonic matching coefficients in the two--mass
variable flavor number scheme to $O(\varepsilon)$.
