Let $Z\sim \mathcal{N}(\mu, \Sigma)$ be a Gaussian random vector in $\mathbb{R}^d$.
What are some nontrivial bounds on $p=\mathbb{P}(S)$, where $S$ is the event $Z_1=\max_i Z_i$?
This is motivated by the inequality
$$\mathbb{E}(|Z_1|\mathbf{1_S})\ge \max\bigg(0,p\,\mathbb{E}|Z_1|-\sqrt{\mu_i^2+\Sigma_{ii}-(\mathbb{E}|Z_1|)^2}\sqrt{p-p^2}\bigg)$$
Suppose there are $n$ variables. If we could easily integrate probabilities over all $n$ variables at once, we would find the exact value for $p$, and would have no need of bounds. Instead we will provide bounds by a single two-variable integral.
We go from $r=n$ to $r=2$, where \begin{align} B^n=0,\text{ and generally }&B^r \text{ is }n\times 1\ \text{ and constant }\\ M^n=1,\text{ and generally }&M^r \text{ is }n\times r\ \text{ and constant }\\ Y^n=Z,\text{ and generally }&Y^r \text{ is }r\times 1\ \text{ and random }\\ \lambda^n=\mu,\text{ and generally }&\lambda^r \text{ is }r \times 1\ \text{ and the mean of }Y,\\ \Pi^r=\Sigma,\text{ and generally }&\Pi^r \text{ is }r\times r\ \text{ and the covariance of }Y, \end{align}
At each step we define $$p=P[Y^r_1=\max M^rY^r+B^r, Y^r \sim N(\lambda^r,\Pi^r)]$$ (we omit superscripts on scalars). We will get upper [lower] bounds by choosing variables so that the $p$'s increase [decrease] as $r$ decreases.
Let $\theta^r_{ij} := \arccos(\Pi^r_{kij}/\sqrt{\Sigma_{ii}\Sigma_{jj}}) = \arccos(\text{cor}(Y^r_i,Y^r_j))$
Let $\ell$ be such that $\lambda^r_\ell=\min_{i>1}\lambda^r_i$.
Let $k$ be such that $\lambda^r_k=\min_{i>1,i\neq \ell}\lambda^r_i$.
Let $h$ be the number distinct from $k$ and $\ell$ which maximizes $\cos(\theta_{hk}+\theta_{h\ell})$. (For lower bounds, find $h$ which maximizes $\cos(|\theta_{hk}-\theta_{h\ell}|)$ instead.)
Let $T$ be the same as $\Sigma$ except that $T_{k\ell}=\sigma_k \sigma_\ell \cos(\theta_{hk}+\theta_{h\ell}) \le \Sigma_{k\ell}$
$T$ has the lowest possible correlation between $Y_k$ and $Y_\ell$ given their correlations with $Y_h$, so it describes a distribution where $Y_h$, $Y_k$, and $Y_\ell$ are linearly dependent.
Let $a,b,c$ be such that $Y_\ell = a Y_h + b Y_k + c$ under $T$.
Let $B^{r-1} = B^r + cM^r_{\cdot k}$.
Let $M^{r-1}_i$ (the ith row of the matrix) be the same as $M^r_i$ but with $M^{r-1}_{ih}=M^r_{ih}+aM^r_{i\ell}$, $M^{r-1}_{ik}=M^r_{ik}+bM^r_{i\ell}$ and $M^{r-1}_{i\ell}$ omitted.
Let $\lambda^{r-1}$ be $\lambda^r$ but with the entry for $Z_\ell$ omitted.
Let $\Pi^{r-1}$ be $T$ but with the row and column for $Z_\ell$ omitted.
Then the new $p$ will be larger than the old $p$, as desired. $\square$