A problem related to stochastic ordering Let $\boldsymbol{X} = (X_1,X_2)^{\rm T}\sim \mathcal{N}_2(\boldsymbol{\mu}, \mathrm{\Sigma})$, where
\begin{eqnarray*}
\boldsymbol{\mu} = (\mu_1, \mu_2)^{\rm T}& = &(\sqrt{\xi_1\xi_2/(\xi_1+\xi_2)}, 0)^{\rm T}\\
\mathrm{\Sigma} & = &\begin{pmatrix} 1 & -\rho\\
-\rho & 1\end{pmatrix}\\
\rho & = & \sqrt{\xi_1\xi_3/(\xi_1+\xi_2)(\xi_2+\xi_3)}.
\end{eqnarray*}
It is given that $\xi_1\leq\xi_2\leq\xi_3$, $\xi_i\geq 0$ and $\sum_{i=1}^3\xi_i = 1$. I have a function
\begin{equation*}
    \pi(\boldsymbol{\mu};\boldsymbol{\xi}) = 1-\mathbb{P}(\boldsymbol{X}\leq \boldsymbol{t}) = 1-\mathbb{P}(X_1\leq t \cap X_2\leq t),
\end{equation*}
where $t>0$. Under the constraints $\xi_1\leq\xi_2\leq\xi_3$, $\xi_i\geq 0$ and $\sum_{i=1}^3\xi_i = 1$, numerically we are getting the maximum of $\pi(\boldsymbol{\mu};\boldsymbol{\xi})$ where the non-zero mean $\mu_1$ is maximized. This occurs when $\xi_1 = \xi_2 = \xi_3 = 1/3$.
Can we use some kind of stochastic ordering arguments to prove the result theoretically?
 A: $\newcommand{\der}{\mathrm{der}}\newcommand{\tdert}{\mathrm{dert}}\newcommand{\erf}{\operatorname{erf}}\newcommand{\eqs}{\overset{\text{sign}}=}\newcommand{\tder}{\widetilde\der}$The answer is no.
Indeed, let $x:=\xi_1$ and $y:=\xi_2$, so that $\xi_3=1-x-y$, $0\le x\le y\le1-x-y$, whence $x\in[0,1/3]$.
Consider further the case $y=x\in(0,1/3]$, so that
$$\mu_1=M(x):=\sqrt{x/2},$$
$$\rho=-R(x),\quad R(x):=\sqrt{\frac{1-2x}{2(1-x)}}.$$
Let also
$$p_t(x):=1-P(X_1\le t,X_2\le t).$$
Your desired result would imply that
\begin{equation*}
    p_t(x)\le p_t(1/3) \tag{0}
\end{equation*}
for all $x\in(0,1/3]$ and all real $t>0$.
Note that
\begin{equation*}
    p_t(x)=1-P(X_1\le t,X_2\le t)=P(M(x),-R(x)), \tag{1}
\end{equation*}
where
\begin{equation*}
    P(m,r):=1-\int_{-\infty}^t du \int_{-\infty}^t dv\, f_{m,r}(u,v)  
\end{equation*}
and $f_{m,r}$ is the density function of the bivariate normal distribution with means $m,0$, variances $1,1$, and correlation $r$.
The key is Plackett's observation (formula (3)) that
\begin{equation*}
    D_r f_{m,r}(u,v)=D_v D_u f_{m,r}(u,v),
\end{equation*}
where $D_w$ denote the partial derivative with respect to a variable $w$. It follows that
\begin{equation*}
\begin{aligned}
    D_r P(m,r)&:=-\int_{-\infty}^t du \int_{-\infty}^t dv \,    D_r f_{m,r}(u,v) \\ 
    &=-\int_{-\infty}^t du \int_{-\infty}^t dv \,   D_v D_u f_{m,r}(u,v) \\ 
    &=-f_{m,r}(t,t).  
\end{aligned} 
\tag{2}
\end{equation*}
Next,
\begin{equation*}
\begin{aligned}
P(m,r)&=1-\int_{-\infty}^t du \int_{-\infty}^t dv\, f_{0,r}(u-m,v)  \\ 
&=1-\int_{-\infty}^{t-m} dw \int_{-\infty}^t dv\, f_{0,r}(w,v) 
\end{aligned} 
\end{equation*}
and hence
\begin{equation*}
\begin{aligned}
    D_m P(m,r)=\int_{-\infty}^t dv\, f_{0,r}(t-m,v)=\int_{-\infty}^t dv\, f_{m,r}(t,v);  
\end{aligned} 
\tag{3}
\end{equation*}
the latter integral can be easily expressed in terms of the error function $\erf$ and elementary functions.
By (1) and a chain rule of differentiation,
\begin{equation*}
    D_x p_t(x)=D_m P(M(x),-R(x))M'(x)-D_r P(M(x),-R(x))R'(x). 
\end{equation*}
In particular,
\begin{equation}
    D_x p_{1/10}(x)\big|_{x=1/3}=\erf\left(\frac{1}{60} \left(3 \sqrt{6}-10\right)\right)+1-\frac{3
   }{\sqrt{\pi }}\,e^{(30 \sqrt{6}-77)/1800}
   =-0.739\ldots<0. 
\end{equation}
So, inequality (0) fails to hold for $t=1/10$ and all $x$ in a left neighborhood of $1/3$. $\quad\Box$
