Let two random variables $X$ and $Y$ are jointly normal as follows:
\begin{eqnarray} \left( \begin{array}{c} X \\ Y \end{array} \right) \sim \text{Normal} \left( \left( \begin{array}{c} \mu_x \\ \mu_y \end{array} \right), \left( \begin{array}{cc} 1 & \rho \\ \rho & 1 \end{array} \right) \right), \end{eqnarray}
where $\rho >0$.
The question is how to calculate $$\text{Cov}(\, \, \mid X\mid \, \, , \, \, \mid Y\mid \, \, )\, \, ?$$
For a special case $\mu_x=\mu_y=0$,
$$\text{Cov}(\, \, \mid X\mid \, \, , \, \, \mid Y\mid \, \, )=\frac{2}{\pi}\left( \sqrt{1-\rho^2}+\rho \arcsin (\rho) -1 \right).$$
I know $\mid X \mid$ is a folded normal random variable (Folded_normal_distribution) so $$E(\mid X \mid)=\sqrt{\frac{2}{\pi}}e^{-\frac{\mu_x^2}{2}}+\mu_x \text{erf}\left(\frac{\mu_x}{\sqrt{2}}\right)$$
and
$$E(\mid Y \mid)=\sqrt{\frac{2}{\pi}}e^{-\frac{\mu_y^2}{2}}+\mu_y \text{erf}\left(\frac{\mu_y}{\sqrt{2}}\right).$$
But I stuck here: $$E(\mid X Y\mid)=$$