Let $\mathcal{E} = \{x\in\Re^n:x^\top P^{-1} x=1\}$ denote an ellipsoid that is centered at origin. The goal is to integrate over the surface of ellipsoid, but only in positive half-spaces $x_1,x_2>0$: \begin{align} I = \int_{x\in \mathcal{E}: x_1,x_2>0} x_1 x_2 \end{align}
Implication: Let $w\sim\mathcal{N}(0,P)$ be a Gaussian random vector, and let $y= w / \|w\|_2$ be its projection to unit-sphere, and finally put $z$ as positive part $z= \frac12 ( y_i + |y_i|)_{i\le n}$. The goal is to now compute the covariance of $z$: $C = E z z^\top $. The answer to the described integral can be used to compute the off-diagonal elements of $C$, and therefore characterize $z$. Here is a summary: \begin{align} &\text{Gaussian random vector }&& w \sim \mathcal{N}(0,P) \\ &\text{project to unit sphere } && y=\frac{w}{\|w\|_2} \\ &\text{take positive part } && z= (y)^{+} \\ &\text{Covariance } && C=\mathbb{E}\ z^\top z \end{align}
A lower bound using convexity: I have tried several approaches to no avail, here I just mention a summary. We can reformulate the integral in the basis of eigenvectors of $P$: $P = U S U^\top$, where $S=\text{diag}(\sigma_1^2,\dots,\sigma_n^2)$, and $u_1,\dots,u_n$ being the eigenvectors: \begin{align} I = \int_{x\in E} x^\top u_1 x^\top u_2 && E:=\{x\in\Re^n: x^\top S^{-1} x=1, x^\top u_1>0, x^\top u_2>0\} \end{align} This does not turn out to be particularly useful, but it at least leads to a rather tight lower bound: if we distribute $M=\sum_i \sigma_i^2$ among the eigenvalues, that will expand the Ellipse and thereby lower the integral. More precisely, $I$ is a convex function in $\sigma_i$s, replacing $\sigma_i^2$ for all $i\neq 1$ by their collective average, its value will be reduced. This will in turn reduce the problem to a 2-dimensional case which is easily integrable. But I wasn't able to establish any meaningful upper bound.