*(Previously I posted a similar question on math.SE, hoping that this question would have an easy answer. As the question appears hard, I am hoping I can perhaps get more feedback here.)*

Let $\mathbf{X} \sim N(\mathbf{0}, \mathbf{\Sigma})$ be a $k$-dimensional Gaussian vector with non-trivial covariance matrix $\mathbf{\Sigma}$. I am interested in $p = \Pr(\mathbf{X} \geq \mathbf{0})$, i.e. the probability that all coordinates are non-negative. This is also known as the *orthant probability* for $\mathbf{X}$, and explicit formulas for orthant probabilities for arbitrary covariance matrices $\mathbf{\Sigma}$ are known in $2$ and $3$ dimensions, with some special cases having been studied in $4$ dimensions, and very few explicit formulas apparently known for higher-dimensional cases. I am interested in these orthant probabilities for various different matrices $\mathbf{\Sigma}$, perhaps the simplest of which is the following:
$$\mathbf{\Sigma} = \begin{pmatrix} 1 & u \\ u & 1 \end{pmatrix} \otimes \begin{pmatrix} 1 & 1/2 & 1/2 \\ 1/2 & 1 & 1/2 \\ 1/2 & 1/2 & 1 \end{pmatrix} = \begin{pmatrix}
1 & 1/2 & 1/2 & u & u/2 & u/2 \\
1/2 & 1 & 1/2 & u/2 & u & u/2 \\
1/2 & 1/2 & 1 & u/2 & u/2 & u \\
u & u/2 & u/2 & 1 & 1/2 & 1/2 \\
u/2 & u & u/2 & 1/2 & 1 & 1/2 \\
u/2 & u/2 & u & 1/2 & 1/2 & 1
\end{pmatrix}$$
Here $u$ is some constant between $0$ and $1$.

Most of the literature I found on this topic (for actually deriving closed-form expressions for these probabilities) dates back to long ago, e.g. works by Cheng, Childs, David, Plackett, Steck in the 1950s and 1960s (mostly in the journal *Biometrika*). I am not so hopeful, but it would be great if

- someone could find a closed form for the above case;
- someone could point me to literature I might have missed on finding closed form expressions for such high-dimensional cases;
- someone could explain which techniques/strategies may generally be useful for finding closed-form expressions for even dimensions.

As for a slightly *fishier* approach, which nonetheless might be rewarding: I know that orthant probabilities for a few low-dimensional cases can be written as follows:
\begin{align}
\mathbf{\Sigma} = \begin{pmatrix} 1 & a \\ a & 1 \end{pmatrix} &\implies p = \frac{\arccos (-a)}{2 \pi}, \\
\mathbf{\Sigma} = \begin{pmatrix} 1 & a & b \\ a & 1 & c \\ b & c & 1 \end{pmatrix} &\implies p = \frac{\arccos (-a)}{4 \pi} + \frac{\arccos (-c)}{4 \pi} - \frac{\arccos (b)}{4 \pi}, \\
\mathbf{\Sigma} = \begin{pmatrix} 1 & a & b & ab \\ a & 1 & ab & b \\ b & ab & 1 & a \\ ab & b & a & 1 \end{pmatrix} &\implies p = \left(\frac{\arccos (-a)}{2 \pi}\right)^2 + \left(\frac{\arccos (-b)}{2 \pi}\right)^2 - \left(\frac{\arccos (a b)}{2 \pi}\right)^2.
\end{align}
Note that the latter case corresponds to:
\begin{align}
\mathbf{\Sigma} = \begin{pmatrix} 1 & a \\ a & 1 \end{pmatrix} \otimes \begin{pmatrix} 1 & b \\ b & 1 \end{pmatrix} = \begin{pmatrix} 1 & a & b & ab \\ a & 1 & ab & b \\ b & ab & 1 & a \\ ab & b & a & 1 \end{pmatrix}
\end{align}
In these cases at least, the probability can be written as a combination of arccosines of the off-diagonal entries, sometimes with minus signs. Alternatively, for the cases of two and four dimensions, we can equivalently express the probabilities in the off-diagonal entries (without changing signs) by considering the inverse matrix $\mathbf{\Sigma}^{-1}$ instead:
\begin{align}
\mathbf{\Sigma}^{-1} = \begin{pmatrix} 1 & -a \\ -a & 1 \end{pmatrix} &\implies p = \frac{\arccos (-a)}{2 \pi}, \\
\mathbf{\Sigma}^{-1} = \begin{pmatrix} 1 & -a & -b & ab \\ -a & 1 & ab & -b \\ -b & ab & 1 & -a \\ ab & -b & -a & 1 \end{pmatrix} &\implies p = \left(\frac{\arccos (-a)}{2 \pi}\right)^2 + \left(\frac{\arccos (-b)}{2 \pi}\right)^2 - \left(\frac{\arccos (a b)}{2 \pi}\right)^2.
\end{align}
Can someone perhaps establish/conjecture a pattern that might allow to extrapolate to higher-dimensional cases?

**Edit 1**: I fixed a mistake in the formula for three dimensions. For references: the two- and three-dimensional cases can, among others, be found in the one-page paper Dav53 or the paper Chi67, while the four-dimensional case stated above appears in the appendix of Che68.

**Edit 2**: For the six-dimensional case, substituting $u = 1/2$ leads to $p \approx 0.115024$. The six decimals seem quite precise; independently, both $10^8$ Monte Carlo experiments in C, and using numerical integration in R, I get these six decimals. So the function $p(u)$ should satisfy $p(0) = 0.25$, $p(1) = 0.0625$, and $p(1/2) \approx 0.115024$.