Vamos a reescribir la integral un poco. Denotar $X$ una variable aleatoria Gaussiana con cero de la media y la varianza $\sigma_1^2$, e $Y$ independiente variable aleatoria Gaussiana con media de $\mu$ y la varianza $\sigma_2^2$. Entonces su integral lee:
$$
\begin{eqnarray}
\int_\mathbb{R} \left( x f_X(x) F_Y(x) + x f_Y(x) F_X(x) \right) \mathbb{x} &=& \int_\mathbb{R} x f_X(x) F_Y(x) \mathbb{x} + \int_\mathbb{R} y f_Y(y) F_X(y) \mathbb{y}
\\ &=& \mathbb{E}\left( X [ Y \leqslant X] + Y [X < Y] \right) \\
&=& \mathbb{E}\left(\max(X,Y) \right) = \mathbb{E}\left(\max(X-Y, 0)+Y \right) \\ &=& \mathbb{E}\left(\max(X-Y, 0) \right) + \mu
\end{eqnarray}
$$
donde $[X<Y]$ denota Iverson soporte. Pero a diferencia de las variables aleatorias $Z := X-Y$ es también una variable aleatoria Gaussiana, como una combinación lineal de Gaussianas, con la media y la varianza
$$
m:= \mathbb{E}(Z) = \mathbb{E}(X) -\mathbb{E}(Y) = -\mu \qquad
s^2 a := \mathbb{Var}(Z) = \mathbb{Var}(X)+\mathbb{Var}(Y) = \sigma_1^2 + \sigma_2^2
$$
Por lo tanto la integral se reduce a la evaluación de $\mathbb{E}\left( \max(Z, 0)\right)$:
$$\begin{eqnarray}
\mathbb{E}\left( \max(Z, 0)\right) &=& \int_0^\infty z f_Z(z) \mathrm{d} z = \int_0^\infty \frac{z}{s} \frac{1}{\sqrt{2 \pi}} \exp\left(-\frac{1}{2} \frac{(z-m)^2}{s^2} \right) \mathrm{d} z \\
&=& m \int_0^\infty f_Z (z) \mathrm{d} z + s \left. \left(-\frac{1}{\sqrt{2\pi}} \exp\left(-\frac{(z-m)^2}{2s^2} \right)\right) \right|_{z=0}^{z=\infty} \\
&=& m \cdot \Phi\left(\frac{m}{s}\right) + s \cdot \phi\left(\frac{m}{s}\right)
\end{eqnarray}
$$
Por lo tanto:
$$\begin{eqnarray}
\int_\mathbb{R} x \left( \frac{1}{\sigma_1} \phi\left(\frac{x}{\sigma_1}\right) \Phi\left(\frac{x-\mu}{\sigma_2}\right) + \frac{1}{\sigma_2} \phi\left(\frac{x-\mu}{\sigma_2}\right) \Phi\left(\frac{x}{\sigma_1}\right) \right) dx = \\ \mu \Phi\left(\frac{\mu}{\sqrt{\sigma_1^2+\sigma_2^2}}\right) + \sqrt{\sigma_1^2+\sigma_2^2} \phi\left(\frac{\mu}{\sqrt{\sigma_1^2+\sigma_2^2}}\right)
\end{eqnarray}
$$
