Assume a pair of normal variables $a,b\sim N(\mu, \Sigma)$, with $\rho_{ab}\neq0$. We know their joint distribution in the (shorthand) form: $$f_X(a,b)=\frac{1}{2\pi\sqrt{|\Sigma|}}exp(-\frac{1}{2}(x-\mu)^T\Sigma^{-1}(x-\mu))$$
Now, let $z=g(a)$ some transformation of $a$. We find the $k^{th}$ moment of $z$ by solving $$\int_{-\infty}^\infty\int_{-\infty}^\infty z^k\cdot f(a,b)da db$$
and $Cov(z,b)$ is given by solving $$\int_{-\infty}^\infty\int_{-\infty}^\infty zb\cdot f(a,b)da db$$
How do I solve this when $g(x)=\Phi(x)=0.5erf(x/\sqrt{2}+1)$?
