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 problem related to stochastic ordering
$\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)](https://www.jstor.org/stable/2332716)) 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$


