 Solution

The proposed probability density function is nonnegative. To check its integral, put $I={\int }_{\mathbb{R}}{e}^{-x^2/2}dx$. The Gaussian integral can be evaluated by squaring and using polar coordinates:Since $I>0$, $I=\sqrt{2\pi }$ and $\int \varphi =1$. Thus $X$ has the standard normal distribution. Completing the square in its moment-generating function gives, for every real $t$,Differentiation under the integral is justified near every finite $t$ by Gaussian decay. Expanding $e^{t^2/2}={\sum }_{m\ge 0}{t}^{2m}/(2^{m}m!)$ gives all moments:

To obtain the two-term bounds for Mills ratio, use ${\varphi }^{\prime }(t)=-t\varphi (t)$. Integration by parts gives, for $x>0$,The positive remainder proves $r(x)<1/x$. A second integration by parts givesSubstitution and division by the strictly positive $\varphi (x)$ showThe last integral is strictly positive for every finite $x>0$, proving $\boxed{1/x-1/x^3<r(x)<1/x}$. Both remainders are finite for every fixed positive $x$; the lower bound remains valid even when it is negative.