$\DeclareMathOperator{\P}{\mathbb{P}}$
Thanks to the fact that $\rho>0$, the following holds \begin{align}p_k(t):&=\P(X_1\le t,\dots,X_k\le t)\\&=\int_{\mathbb{R}}\Phi\left(\frac{t - \sqrt{\rho} z}{\sqrt{1 - \rho}}\right)^k\phi(z)\,dz,\end{align}
where $\Phi$ and $\phi$ are the CDF and the PDF of a $\mathop{N}(0,1)$ - distributed random variable, respectively. Since $\phi$ is symmetric around $0$, we get
$$p_k(t)=\int_{\mathbb{R}}\Phi\left(\frac{t + \sqrt{\rho} z}{\sqrt{1 - \rho}}\right)^k\phi(z)\,dz.$$
For $x,y\in\mathbb{R}^d$ define the notations
\begin{align}
x\wedge y&\equiv\bigl(\min(x_1,y_1),\min(x_2,y_2),\dots,\min(x_d,y_d)\bigr),\\
x\vee y&\equiv\bigl(\max(x_1,y_1),\max(x_2,y_2),\dots,\max(x_d,y_d)\bigr).
\end{align}
Recall some definitions
- A function $f\colon \mathbb{R}^d\to[0,+\infty)$ is $\mathrm{MTP_2}$ (i.e multivariate totally positive of order 2) if $$f(x\wedge y)f(x\vee y)\ge f(x) f(y),\;\;\forall x,y\in \mathbb{R}^d.$$
- A function $f\colon \mathbb{R}^d\to[0,+\infty)$ is less than another function $g\colon \mathbb{R}^d\to[0,+\infty)$ in the $\mathrm{MTP_2}$ sense if
$$f(x\wedge y)g(x\vee y)\ge f(x) g(y),\;\;\forall x,y\in \mathbb{R}^d,$$
and we will write $f\preceq g$.
In the unidimensional case $d=1$, we have
$$f\preceq g\iff \frac{f(x)}{g(x)}\;\text{is non-decreasing},$$ and that is why the $\mathrm{MTP_2}$ ordering is also called the multivariate likelihood ratio ordering. Recall the following result
If $f(z)\preceq g(z)$ and $H(x,z)\preceq L(x,z)$, then
\begin{equation}\int H(x,z)f(z)\,dz\preceq \int L(x,z)g(z)\,dz.\label{eq:0}\tag1.\end{equation}
Note that $r_k(t)$ is non-decreasing iff $p_{k+1}(t)\preceq p_k(t)$.
We want to apply $\eqref{eq:0}$ with $f(z)=g(z)=\phi(z)$, and $H(t,z)=\Phi\left(\frac{t + \sqrt{\rho} z}{\sqrt{1 - \rho}}\right)^{k+1}$ and $L(t,z)=\Phi\left(\frac{t + \sqrt{\rho} z}{\sqrt{1 - \rho}}\right)^k$. Clearly, $f\preceq g$, because all functions of one variable are $\mathrm{MTP_2}$.
To see why $H(t,z)\preceq L(t,z)$, note that $\Phi(\cdot)^{k+1}\preceq\Phi(\cdot)^{k}$ and use the fact that the transformation $(t,z)\mapsto\frac{t + \sqrt{\rho} z}{\sqrt{1 - \rho}}$ is coordinatewise increasing in both $t$ and $z$.
Karlin, Samuel; Rinott, Yosef, Classes of orderings of measures and related correlation inequalities. I. Multivariate totally positive distributions, J. Multivariate Anal. 10, 467-498 (1980). ZBL0469.60006.
