How to derive integrals with error function? How to derive this integral
$\int_{-\infty}^{\infty}erf(\lambda x)\mathcal{N}(\mu, \sigma ^2)dx$
and this 
$\int_{-\infty}^{\infty}(erf(\lambda x)-const)^2\mathcal{N}(\mu, \sigma ^2)dx$
where
$erf(x)=\frac{2}{\sqrt\pi}\int_{0}^{x}{e}^{-t^2}dx$ - is the error function
$\mathcal{N}(\mu, \sigma ^2)=\frac{1}{\sigma\sqrt{2\pi}}{e}^{-\frac{(x-\mu)^2}{2\sigma^2}}$ - is pdf of the normal distribution
From this paper: "...it can be shown by parameter differentiation with respect to mu and then integrating with respect to mu..." they have got quite simple result but I really don't understand how this trick works.
Edit:
So in order to summarize and make everything clear:
My original goal was deriving aforesaid integrals. It is much easier to do if one knows $\int_{-\infty}^{\infty}\Phi(\lambda x)\mathcal{N}(\mu, \sigma ^2)dx$ and $\int_{-\infty}^{\infty}\Phi^2(\lambda x)\mathcal{N}(\mu, \sigma ^2)dx$
Based on what @ir7 has shown $\int_{-\infty}^{\infty}\Phi(\lambda x)\mathcal{N}(\mu, \sigma ^2)dx =\Phi\left(\frac{\lambda \mu}{\sqrt{1+\lambda^2\sigma^2}}\right)$ and from definition  $erf(x)=2\Phi(\sqrt{2}x)-1$ and $\Phi(x)=\frac{1}{2}erf(\frac{x}{\sqrt{2}})+\frac{1}{2}$ we obtain:
\begin{equation}\int_{-\infty}^{\infty}erf(\lambda x)\mathcal{N}(\mu, \sigma ^2)dx = \int_{-\infty}^{\infty}(2\Phi(\lambda\sqrt{2}x)-1)\mathcal{N}(\mu, \sigma ^2)dx=2\int_{-\infty}^{\infty}\Phi(\lambda\sqrt{2}x)\mathcal{N}(\mu, \sigma ^2)dx-1=2\Phi\left(\frac{\lambda\sqrt{2}\mu}{\sqrt{1+2\lambda^2\sigma^2}}\right)-1=erf\left(\frac{\lambda\mu}{\sqrt{1+2\lambda^2\sigma^2}}\right)\end{equation}
\begin{equation}\int_{-\infty}^{\infty}erf(\lambda x)\mathcal{N}(\mu, \sigma ^2)dx = erf\left(\frac{\lambda\mu}{\sqrt{1+2\lambda^2\sigma^2}}\right)\end{equation}
In the case of $\int_{-\infty}^{\infty}\Phi^2(\lambda x)\mathcal{N}(\mu, \sigma ^2)dx$ ...
 A: Hint: The "trick" they are refering to is differentiation under integral. We have:
$$\frac{d}{dy} \int_a^b f(x,y)dx =\int_a^b  \frac{\partial}{\partial y} f(x,y)dx,$$
if $f(x,y)$ and $\frac{\partial}{\partial y} f(x,y)$ are continuous in both variables and are, respectively, bounded in absolute value by two functions $g(x)$ and $h(x)$, both independent of $y$, such that $\int_a^b g(x)dx$ and $\int_a^b h(x)dx$ exist and are finite.
Note that ($\sigma >0$, $\lambda >0$):
$$   F(\sigma,\mu) \triangleq \int_{-\infty}^{\infty} \Phi(\lambda x)\frac{1}{\sigma} {\cal N}\left(\frac{x-\mu}{\sigma}\right) dx = \int_{-\infty}^{\infty} \Phi( \lambda\sigma y+\lambda\mu) {\cal N}\left(y\right) dy .$$
where $\cal N$ is standard normal probability density function and $\Phi$ is standard normal cumulative distribution function. 
Using the diferentiation under integral, its derivative wrt $\mu$ is
$$  \frac{\partial}{\partial \mu} F(\sigma,\mu) = \int_{-\infty}^{\infty}  \lambda{\cal N}(\lambda\sigma y+\lambda\mu) {\cal N}\left(y\right) dy =
\frac{\lambda}{\sqrt{1+\lambda^2\sigma^2}} {\cal N}\left(\frac{\lambda\mu}{\sqrt{1+\lambda^2\sigma^2}}\right). $$
Finally, integrating wrt to $\mu$ from $-\infty$ to $\tilde{\mu}$ and noting that $\lim\limits_{\tilde{\mu}\to -\infty}F(\sigma, \tilde{\mu}) = 0$, we get:
$$ F(\sigma, \tilde{\mu}) = \Phi\left(\frac{ \lambda\tilde{\mu}}{\sqrt{1+\lambda^2\sigma^2}}\right),$$
as per paper's statement.
You can check by yourself that the conditions of the theorem are met, the equalities I stated without full details are true, and adapt the result to the error function. 
EDIT1: The following general formulae should help for the second integral (using the same method):
$$\int_{-\infty}^{\infty}  \exp(Ay) \Phi(By+C) {\cal N}_{d,e^2}\left(y\right) dy  = 
\exp\left(Ad+0.5e^2A^2\right)\Phi\left(\frac{ABe^2+Bd+C}{\sqrt{1+B^2e^2}} \right),$$
$$\int_{-\infty}^{x}  \exp(Ay) \Phi(By+C) {\cal N}_{d,e^2}\left(y\right) dy  = 
\exp\left(Ad+0.5e^2A^2\right)\Phi_2\left(\frac{x-d}{e}-Ae,\frac{ABe^2+Bd+C}{\sqrt{1+B^2e^2}}; \frac{-Be}{\sqrt{1+B^2e^2}}\right),$$
$$\int_{-\infty}^x \Phi(y) dy = x\Phi(x) +  {\cal N}(x), $$
where $\Phi_2(\cdot,\cdot;\rho)$ is the cdf of bivariate standard normal with correlation $\rho$:
$$ \Phi_2(x,y;\rho) =\frac{1}{2\pi\sqrt{1-\rho^2}}\int_{-\infty}^x \int_{-\infty}^y \mathrm e^{-.5\frac{u^2-2\rho uv+v^2}{1-\rho^2}}dudv.$$ 
So, denote $$G\left(\sigma, \mu\right)\triangleq \int_{-\infty}^{\infty}  \Phi(\lambda\sigma y+\lambda\mu)^2 {\cal N}\left(y\right) dy.$$
Then $$\frac{\partial}{\partial \mu} G(\sigma,\mu) = \int_{-\infty}^{\infty}  2\lambda \Phi(\lambda\sigma y+\lambda\mu) {\cal N}(\lambda\sigma y+\lambda\mu){\cal N}\left(y\right) dy $$ 
$$ = \int_{-\infty}^{\infty}  2\lambda \Phi(\lambda\sigma y+\lambda\mu){\cal N}\left(\sqrt{1+\lambda^2\sigma^2} y+\frac{\lambda^2\sigma\mu}{\sqrt{1+\lambda^2\sigma^2}}\right) {\cal N}\left(\frac{\lambda\mu}{\sqrt{1+\lambda^2\sigma^2}}\right)dy$$
$$ = 2\lambda {\cal N}\left(\frac{\lambda\mu}{\sqrt{1+\lambda^2\sigma^2}}\right)
\int_{-\infty}^{\infty}   \Phi(\lambda\sigma y+\lambda\mu){\cal N}\left(\sqrt{1+\lambda^2\sigma^2} y+\frac{\lambda^2\sigma\mu}{\sqrt{1+\lambda^2\sigma^2}}\right) dy $$
$$ =\frac{2\lambda}{\sqrt{1+\lambda^2\sigma^2} } {\cal N}\left(\frac{\lambda\mu}{\sqrt{1+\lambda^2\sigma^2}}\right)
\int_{-\infty}^{\infty}   \Phi\left(\lambda\sigma \frac{z-\frac{\lambda^2\sigma\mu}{\sqrt{1+\lambda^2\sigma^2}}}{\sqrt{1+\lambda^2\sigma^2} }+\lambda\mu\right){\cal N}\left(z\right) dz $$
$$=\frac{2\lambda}{\sqrt{1+\lambda^2\sigma^2} } {\cal N}\left(\frac{\lambda\mu}{\sqrt{1+\lambda^2\sigma^2}}\right)\Phi\left(\frac{\lambda\mu}{\sqrt{1+\lambda^2\sigma^2}}\frac{1}{\sqrt{1+2\lambda^2\sigma^2}} \right). $$
It follows:
$$G(\sigma, \tilde{\mu}) = \frac{2\lambda}{\sqrt{1+\lambda^2\sigma^2} } \int_{-\infty}^{\tilde{\mu}}{\cal N}\left(\frac{\lambda\mu}{\sqrt{1+\lambda^2\sigma^2}}\right)\Phi\left(\frac{\lambda\mu}{\sqrt{1+\lambda^2\sigma^2}}\frac{1}{\sqrt{1+2\lambda^2\sigma^2}} \right) d\mu$$
$$ =2 \int_{-\infty}^{\frac{\lambda\tilde{\mu}}{\sqrt{1+\lambda^2\sigma^2}}}{\cal N}\left(z\right)\Phi\left(\frac{z}{\sqrt{1+2\lambda^2\sigma^2}} \right) dz $$
$$ = 2\Phi_2\left(\frac{\lambda\tilde{\mu}}{\sqrt{1+\lambda^2\sigma^2}},0; -\frac{1}{\sqrt{2+2\lambda^2\sigma^2}}\right). $$
These calculations need to be verified.
Alternatively, the probabilistic interpretation of $G(\sigma,\mu)$ is
$$\mathbb{E}\left[\Phi\left(\lambda\sigma X+\lambda\mu\right)^2\right],$$ and has been partially studied here Closed forms for various expectations involving the standard normal CDF (note Stein's Lemma referenced therein).
EDIT2: 
$$\int_{-\infty}^{\infty}  \lambda{\cal N}(\lambda\sigma y+\lambda\mu) {\cal N}\left(y\right) dy= \lambda \int_{-\infty}^{\infty} \frac{1}{\sqrt{2\pi}}\mathrm e^{-.5(\lambda\sigma y+\lambda\mu)^2}\frac{1}{\sqrt{2\pi}}\mathrm e^{-.5y^2} dy$$
$$ = \lambda \int_{-\infty}^{\infty} \frac{1}{\sqrt{2\pi}}\mathrm e^{-.5\left(\sqrt{1+\lambda^2\sigma^2} y+\frac{\lambda^2\sigma\mu}{\sqrt{1+\lambda^2\sigma^2}}\right)^2}\frac{1}{\sqrt{2\pi}}\mathrm e^{-.5\left(\frac{\lambda\mu}{\sqrt{1+\lambda^2\sigma^2}}\right)^2}  dy$$
$$= \lambda {\cal N}\left(\frac{\lambda\mu}{\sqrt{1+\lambda^2\sigma^2}}\right) \int_{-\infty}^{\infty} \frac{1}{\sqrt{2\pi}}\mathrm e^{-.5\left(\sqrt{1+\lambda^2\sigma^2} y+\frac{\lambda^2\sigma\mu}{\sqrt{1+\lambda^2\sigma^2}}\right)^2}dy$$
$$ =\lambda {\cal N}\left(\frac{\lambda\mu}{\sqrt{1+\lambda^2\sigma^2}}\right) \frac{1}{\sqrt{1+\lambda^2\sigma^2}}, $$
with last equality coming from a variable change.
