Let $R\left ( x\right ) $ be the Mills ratio of the standard Gaussian law. This paper investigates the higher order monotonicity of the function \( x\mapsto Y_{a,b}\left ( x\right ) =\left ( x^{2}+a\right ) R^{2} \left ( x\right ) +bxR\left ( x\right ) \)on $\left ( 0,\infty \right ) $, where $a,b\in \mathbb{R}$. In particular, we prove that the function $x\mapsto Y_{a,a-2}\left ( x\right ) -a+1$ is completely monotonic on $\left ( 0,\infty \right ) $ if and only if $a\geq -1$. As a special case, $x\mapsto Y_{-1,-3}\left ( x\right ) +2$ is completely monotonic on $\left ( 0,\infty \right ) $, which implies the complete monotonicity of the functions $R^{3}\left ( 1/R\right ) ^{\prime \prime }$ and $-R^{3}\left ( \ln R\right ) ^{\prime \prime \prime }$ on the same interval. As applications, several functional inequalities are derived, and some known and new sharp bounds for R are established.