Gaussian Measure of a Simplex

Question

I'm interested in evaluating the following integral $$ I(n) =\int_{-\infty}^0\int_{-\infty}^{z_2 \sqrt{\frac{3-2}{3}}}\cdots \int_{-\infty}^{z_{n-1}\sqrt{\frac{n-2}{n}}}\exp\left(-\sum_{k=2}^n z_k^2\right)dz_n\cdots dz_3dz_2 $$ where $n\geq 2$. Clearly for $n=2$, this is just (half of) a standard Gaussian integral, and gives $$ I(2) = \frac{\sqrt{\pi}}{2} $$ In addition, Mathematica gives $$ I(3) = \frac{\pi}{3\cdot 2} $$ Comparison with numerical results seems to suggest that in general $$ I(n) = \frac{\pi^{\frac{n-1}{2}}}{n!} $$ I.e. that this simplex captures a fraction (1/n!) of the entire Gaussian weight ($\pi^{\frac{n-1}{2}})$. This seems like a very clean result, but I'm absolutely at a loss for how to prove this. A proof by induction doesn't pan out (unless there's some neat integration trick that simplifes the resulting integral of a Gaussian and an error function...), and I'm not aware of any symmetry arguments that could simplify this. Any ideas are welcome.

Edit: If you'd like to test this in Mathematica, you can use the following function

(*Define the function to compute the integral*)
SimplexIntegral[n_Integer] := 
 Module[{vars, integrand, limits, I, sum, 
   i},(*Create a list of variables z2,z3,...,zn*)
  vars = Table[Symbol["z" <> ToString[k]], {k, 2, n}];
  (*Define the integrand as the exponential of the sum of squares*)
  integrand = Exp[-Sum[vars[[k - 1]]^2, {k, 2, n}]];
  (*Set up the integration limits:first from-\[Infinity] to 0,
  then based on the previous variable with the sqrt factor*)
  limits = 
   Join[{{vars[[1]], -Infinity, 0}}, 
    Table[{vars[[i]], -Infinity, 
      vars[[i - 1]] Sqrt[(i - 1)/(i + 1)]}, {i, 2, Length[vars]}]];
  (*Perform the numerical integration*)
  I = NIntegrate[integrand, Evaluate[Sequence @@ limits], 
    Method -> "GlobalAdaptive", MaxRecursion -> 15, 
    PrecisionGoal -> 5, WorkingPrecision -> 25];
  (*Return the result*)I]

and you can plot the comparison via

DiscretePlot[(Factorial[n] SimplexIntegral[n])/\[Pi]^((n - 1)/2), {n, 2, 5}]

An Idea: Due to the spherical symmetry of the Gaussian kernel, we can view $I(n)$ as the total Gaussian weight times the ratio $$ \frac{\Omega_{simplex}}{\Omega_{n-1}}\equiv\tilde{\Omega}_{simplex} $$ where $\Omega_{simplex}$ refers to the solid angle of the vertex of the simplex (defined by $-\infty \leq z_2 \leq 0,-\infty \leq z_3\leq z_2\sqrt{\frac{3-2}{3}},\cdots,-\infty \leq z_n\leq z_{n-1}\sqrt{\frac{n-2}{n}}$) at the origin, and $\Omega_{n-1}$ is the solid angle of $\mathbb{R}^{n-1}$. This answer suggests that the integral for $\tilde{\Omega}_{simplex}$ may be written as a multivariate Taylor series $$ \begin{align*} \tilde{\Omega}_{simplex}&= \frac{n^{n/2}}{n!(4\pi)^{\frac{n-1}{2}}}\sum_{\boldsymbol{a}\in \mathbb{N}^{\binom{n-1}{2}}}\left[\prod_{i=1}^{n-1}\Gamma\left(\frac{1 + \sum_{m>i}a_{im} + \sum_{m<i}a_{mi}}{2}\right)\right]\prod_{i<j}\frac{s_{ij}^{a_{ij}}}{a_{ij}!}\\ &\equiv \frac{n^{n/2}}{n!(4\pi)^{\frac{n-1}{2}}}S(n) \end{align*} $$ where $$ s_{ij} = -2 \sqrt{\frac{i}{j}\frac{n-j}{n-i}}\hspace{5mm}\text{for }i<j $$ and note that the vector $\boldsymbol{a}$ consists of only the superdiagonal entries of a $(n-1)\times (n-1)$ matrix, hence it is defined as $\boldsymbol{a} = \{a_{ij}\}_{i,j=1, i<j}^{n-1}$ (note that $\binom{n-1}{2}$ is the number of such superdiagonal entries). This result comes from the fact that the normalized vectors $\overline{\boldsymbol{v}}_1,\cdots, \overline{\boldsymbol{v}}_n $ that define the simplex are given by $$ \overline{\boldsymbol{v}}_i = \sum_{k=i}^{n-1}\sqrt{\frac{i n}{k(k+1)(n-i)}}\textbf{e}_k $$ which implies that, if we let $\boldsymbol{V}\in \mathbb{R}^{(n-1)\times(n-1)}$ be the matrix with the $\overline{\boldsymbol{v}}_i$ as columns $$ \det(\boldsymbol{V}) = \frac{n^{n/2}}{n!} $$ The question then becomes how to evaluate this multivariable Taylor series, which admittedly doesn't seem feasible. To obtain the numerically observed result, we'd have to have $$ S(n) = \frac{(4\pi)^{\frac{n-1}{2}}}{n^{n/2}} $$ For $n=3$, we have $$ S(3) = \sum_{a_{12}=0}^\infty \Gamma\left(\frac{1 + a_{12}}{2}\right)^2\frac{(-1)^{a_{12}}}{a_{12}!} $$ and Mathematica gives $S(3)=\frac{4\pi}{3^{3/2}}$. For $n=4$, the sum is much more complicated: $$ S(4) = \sum_{a_{12}\geq 0}\sum_{a_{13}\geq 0}\sum_{a_{23}\geq 0}\Gamma\left(\frac{1 + a_{12} + a_{13}}{2}\right)\Gamma\left(\frac{1 + a_{23} + a_{12}}{2}\right)\Gamma\left(\frac{1 + a_{13} + a_{23}}{2}\right)\frac{\left(-\frac{2}{\sqrt{3}}\right)^{a_{12} + a_{23}}}{a_{12}!a_{23}!}\frac{\left(-\frac{2}{3}\right)^{a_{13}}}{a_{13}!} $$

In Terms of Joint Probability: I don't immediately see how this is useful, but the integral can be neatly written as $$ I(n)/\pi^{\frac{n-1}{2}} = P\left(Z_2>0, Z_3>Z_2 \sqrt{\frac{3-2}{3}}, \cdots, Z_n > Z_{n-1}\sqrt{\frac{n-2}{n}}\right) $$ where $Z_2,\cdots, Z_n =\mathcal{N}(0,1)$ are independent and identically distributed. And somehow we have, coincidently $$ I(n)/\pi^{\frac{n-1}{2}} = \frac{1}{n}P\left(Z_2>0, Z_3>Z_2, \cdots, Z_n > Z_{n-1}\right) $$ For $n=2$, we have $$ I(2)/\sqrt{\pi} = P(Z_2>0)=\frac{1}{2}\equiv \frac{1}{2!} $$ for $n=3$ $$ I(3)/\pi = P\left(Z_2>0, Z_3>Z_2\sqrt{\frac{3-2}{3}}\right) = \frac{1}{4} - T\left(0, \sqrt{\frac{3-2}{3}}\right) = \frac{1}{6}\equiv \frac{1}{3!} $$ where $T(h, a)$ is Owen's T function

Answer 1

Here is a thought (too long for a comment):

Let $(Z_k)_k$ IID standard normal and $c_k:=\sqrt{\frac{k-2}{k}}$. Fix $n \geq 3$. Let $\|Z\|:=\sqrt{Z_2^2+...+Z_n^2}$. Recall $(\frac{Z_2}{\|Z\|},...,\frac{Z_n}{\|Z\|})$ is uniformly distributed on the unit $n-1$-sphere $\partial B^{n-1}$. We have for example: $$P(Z_2<0,Z_3<Z_2c_3)=P\bigg(\frac{Z_2}{\|Z\|}<0,\frac{Z_3}{\|Z\|}<\frac{Z_2}{\|Z\|}\frac{1}{\sqrt{3}}\bigg)=\frac{1}{6}=\frac{1}{3!}$$ because it is the ratio between an arc length of 60 degrees and the circumference of the unit circle.

One would be done if it were true that the ratio between the surface area of $A_n:=\{x \in \partial B^{n-1}:x_1<0,...,x_{n-1}<c_{n-1}x_{n-2}\}$ and the total surface area of $\partial B^{n-1}$ is $\frac{1}{n!}$.