Came across the integral
$$\int_{x_1, x_2 > 0} x_1 x_2 dX$$
where $X \sim \mathcal{N}(0, \Sigma)$ is a 2D Gaussian. So we could write it as
$$= \int_{x_1, x_2 > 0} x_1x_2\left(\frac{1}{2\pi}\exp\left( - \frac12 x^T\Sigma^{-1}x \right)\right)dx_1dx_2$$
The claim is that if $\Sigma = \begin{bmatrix} a & b \\ b & d\end{bmatrix}$, and we abbreviate $\theta := \cos^{-1}\left( \frac{b}{\sqrt{ad}} \right)$, then the answer is
$$\frac{1}{2\pi} \sqrt{ad} (\sin \theta + (\pi - \theta) \cos \theta).$$
In the case that $\Sigma = I$, it works, because then we get
$$\int_{0}^\infty \int_0^\infty x_1x_2 \left(\frac{1}{2\pi} \exp\left(- \frac{1}{2}(x_1^2 + x_2^2)\right) \right) dx_1dx_2$$
$$= \frac{1}{2\pi} \cdot \left( \int_0^\infty x \exp\left( - \frac12 x^2 \right) \right)^2 = \frac{1}{2\pi}$$
which agrees when $\theta = \cos^{-1}(0) = \pi/2$.
I've been banging my head against the proof for a while with no success. Please help!
(I would link to the solution that I've found, but you don't want to see it. Everything is phrased totally differently, and it would be a slog to see the claim in the weeds. If you really want, it's here on page 12, linking to this on page 2.)
The general approach here is a change of variables. We have using standard notation for the bivariate normal distribution
$$\Sigma = \pmatrix{\sigma_1^2 & \rho\sigma_1 \sigma_2 \\ \rho \sigma_1 \sigma_2 & \sigma_2^2} \\ \Sigma^{-1} = \frac{1}{1 - \rho^2}\pmatrix{\sigma_1^{-2} & -\rho\sigma_1^{-1} \sigma_2^{-1} \\ -\rho \sigma_1^{-1} \sigma_2^{-1} & \sigma_2^{-2}}, $$
which translates to $\sigma_1 = \sqrt{a}$, $\sigma_2 = \sqrt{d}$ and $\rho = b/\sqrt{ad}$.
Since $\Sigma^{-1}$ is a positive definite matrix, there is a lower triangular, nonsingular matrix $C$ such that
$$C'\Sigma^{-1}C = I \\ |C'||\Sigma^{-1}||C| = 1 \implies |C| = |\Sigma|^{1/2} \\ I = I^{-1} = (C'\Sigma^{-1}C)^{-1} = C^{-1} \Sigma (C')^{-1} \implies CC' = \Sigma.$$
The actual form is
$$C = \pmatrix{ \sigma_1 & 0 \\ \rho \sigma_2 & \sigma_2\sqrt{1-\rho^2}} $$
You are trying to find the off-diagonal component of
$$ K =\frac{|\Sigma|^{-1/2}}{2 \pi}\int_D \mathbb{x} \mathbb{x}'\exp\left(-\frac{1}{2}\mathbb{x}'\Sigma^{-1}\mathbb{x}\right) \, d \mathbb{x},$$
where $D = [0,\infty) \times [0,\infty)$.
Make the change of variables $\mathbb{x} = C\mathbb{z}.$ The integral transforms to
$$K = \frac{|\Sigma|^{-1/2}}{2 \pi}C\int_{D'} \mathbb{z}' \mathbb{z}\exp\left(-\frac{1}{2}\mathbb{z}'\mathbb{z}\right) \, |C|d \mathbb{z}C' = C\frac{1}{2 \pi}\int_{D'} \mathbb{z}' \mathbb{z}\exp\left(-\frac{1}{2}\mathbb{z}'\mathbb{z}\right) \, d \mathbb{z}C'.$$
The integrand now factors into the form $f(z_1)g(z_2)$. You just have to (1) determine how the region $D$ transforms to $D'$ under the change of variables, (2) carry out the integration of all four components, and (3) pre- and post-multiply by $C$ and the transpose $C'$.
In this case, we have $D' = \{(z_1,z_2) \, |\, 0 \leqslant z_1 < \infty, \, -(\rho/\sqrt{1- \rho^2})z_1 \leqslant z_2 < \infty \}$ and the problem reduces to computing the integrals
$$\int_0^ \infty e^{-z_1^2/2} \left(\int_{-\left(\rho/\sqrt{1- \rho^2}\right)z_1}^ \infty z_i \, z_j e^{-z_2^2/2}\, dz_2\right) \, dz_1$$
for $i, j = 1, 2$. Using polar coordinates we get integrals of the form
$$\int_0^\infty r^3 e^{-r^2/2} \, dr \int_{- \tan^{-1}\left(\rho / \sqrt{1-\rho^2}\right)}^{\pi/2}\cos^\alpha \theta \sin^\beta \theta \, d\theta$$