I am playing around with some conditional distribution in R and I have run into the integral of the form $$ G(x) := \int_{m}^{\infty} \frac{1}{\sqrt{2 \pi} \sigma} \phi \left( \frac{z - x}{\sigma} \right) \Phi \left( \frac{z}{\sqrt{1 + \sigma^2}} \right) dz,$$ where $m > 0$, $\sigma > 0$, and $x \in \mathbb{R}$. Equivalently, the integrand is the product of the PDF of a $N (0, \sigma^2)$ variable, evaluated at $z - x$, and the CDF of a $N (0, 1 + \sigma^2)$ variable, evaluated at $z$.
In particular, I am interested in this integral as a function of $x$. What I have observed numerically, is that for any $m \in \mathbb{R}$ (particularly $m \in [0, 1]$), and for any $\sigma > 0$, there seems to exist some $\omega = \omega(m, \sigma) \in (0, 1)$ such that $$ G(x) \approx \Phi \left( \frac{x}{\sqrt{1 + 2 \sigma^2}} \right) \Phi \left( \frac{x - m}{\sigma} \right)^{\omega}, \quad \quad \text{for any } x \in \mathbb{R}. $$
What I am looking for is any insight as to why this approximation holds. In particular I am trying to construct an analytic proof of this approximation.
I am aware that the integral has a closed form solution, in terms of $\Phi$ for $m = -\infty$. I am also aware that there is no closed form solution for $G (x)$, outside of using a Bivariate Normal CDF or the Owen T function. I have tried substituting in all sorts of known approximations to $\Phi$, all to no avail. Any guidance would be greatly appreciated!