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.
