Find closed form for quadruple integral

This integral is up to normalization constant integral of multivariate gaussian distribution. Due too lack of cross terms we get that there are 4 independent zero mean gaussian variables involved.

We can begin by rewriting this a probability of an event.

Let $X,Y,Z,W$ independent normally distributed with zero mean and possibly different variances.

Then integral reduces to: $$ \mathbb{P} (0 < W < Z < Y < X) $$

Which in turn is the same as:

$$ \mathbb{P}(W > 0 \land Z-W >0 \land Y-Z >0 \land X-Y >0) $$

and one can observe that $W,Z-W,Y-Z,X-Y$ are correlated jointly normal variables. Probability that all components of a jointly normal vector are positive is called orthant probability and in general doesn't have closed form expression. I suppose this is quite special case with covariance matrix almost diagonal so maybe there are some articles about how to approach this special case.

For the case with 3 or less variables formulas are known (cf this question for example) and I guess they would coincide with what you've found.

Define: \begin{eqnarray} I(a,b,c,d):=\int_0^{\tan ^{-1}\left(\sqrt{\frac{b}{a}}\right)} \frac{\sin (\theta ) \tan ^{-1}\left(\sin (\theta ) \sqrt{\frac{d}{b+c \sin ^2(\theta )}}\right)}{\sqrt{a b d \left(b+c \sin ^2(\theta )\right)}} \, d\theta \quad (i) \end{eqnarray} Then also define:

\begin{eqnarray} {\mathfrak F}^{(A,B)}_{a,b} &:=& \int\limits_A^B \frac{\log(z+a)}{z+b} dz\\ &=& F[B,a,b] - F[A,a,b] + 1_{t^* \in (0,1)} \left( -F[A+(t^*+\epsilon)(B-A),a,b] + F[A+(t^*-\epsilon)(B-A),a,b] \right) \quad (ii) \end{eqnarray} where \begin{eqnarray} t^*:=-\frac{Im[(A+b)(b^*-a^*)]}{Im[(B-A)(b^*-a^*)]} \end{eqnarray} and \begin{equation} F[z,a,b] := \log(z+a) \log\left( \frac{z+b}{b-a}\right) + Li_2\left( \frac{z+a}{a-b}\right) \end{equation} for $a$,$b$,$A$,$B$ being complex.

Then we have: \begin{eqnarray} I(a,b,c,d)&=& \frac{1}{\sqrt{a}} \int\limits_0^{\sqrt{\frac{d}{a+b+c}}} \frac{u \tan ^{-1}(u)}{\left(d-c u^2\right) \sqrt{d-u^2 (b+c)}} du\\ &=& \frac{\sqrt{d}}{\sqrt{a} (b+c)}\int\limits_0^{\sin ^{-1}\left(\sqrt{\frac{b+c}{a+b+c}}\right)} \frac{\sin (\phi ) \tan ^{-1}\left(\sin (\phi ) \sqrt{\frac{d}{b+c}}\right)}{d-\frac{c d \sin ^2(\phi )}{b+c}} d\phi\\ &=&-\frac{2 i}{\sqrt{a} \sqrt{d}} \int\limits_0^{\frac{\sqrt{\frac{b+c}{a+b+c}}}{\sqrt{\frac{a}{a+b+c}}+1}} \frac{t}{b \left(t^2+1\right)^2+c \left(t^2-1\right)^2} \log \left(\frac{2 i t \sqrt{\frac{d}{b+c}}+t^2+1}{-2 i t \sqrt{\frac{d}{b+c}}+t^2+1}\right) dt\\ &=&\frac{1}{4} \frac{1}{\sqrt{a b c d}} \sum\limits_{\xi=1}^4 \sum\limits_{\eta=1}^4 (-1)^{\left\lfloor \frac{\eta -1}{2}\right\rfloor +\left\lfloor \frac{\xi -1}{2}\right\rfloor } \int\limits_0^{\frac{\sqrt{\frac{b+c}{a+b+c}}}{\sqrt{\frac{a}{a+b+c}}+1}} \frac{\log \left(i (-1)^{\left\lfloor \frac{\xi -1}{2}\right\rfloor } \sqrt{\frac{d}{b+c}}+i (-1)^{\xi -1} \sqrt{\frac{b+c+d}{b+c}}+t\right)}{t-i (-1)^{\left\lfloor \frac{\eta -1}{2}\right\rfloor +\eta +1} e^{i (-1)^{\left\lfloor \frac{\eta -1}{2}\right\rfloor } \tan ^{-1}\left(\frac{\sqrt{c}}{\sqrt{b}}\right)}} dt \\ &=&\frac{1}{4} \frac{1}{\sqrt{a b c d}} \sum\limits_{\xi=1}^4 \sum\limits_{\eta=1}^4 (-1)^{\left\lfloor \frac{\eta -1}{2}\right\rfloor +\left\lfloor \frac{\xi -1}{2}\right\rfloor } % {\mathfrak F}^{(0,\frac{\sqrt{\frac{b+c}{a+b+c}}}{\sqrt{\frac{a}{a+b+c}}+1})}_{i (-1)^{\left\lfloor \frac{\xi -1}{2}\right\rfloor } \sqrt{\frac{d}{b+c}}+i (-1)^{\xi -1} \sqrt{\frac{b+c+d}{b+c}},-i (-1)^{\left\lfloor \frac{\eta -1}{2}\right\rfloor +\eta +1} e^{i (-1)^{\left\lfloor \frac{\eta -1}{2}\right\rfloor } \tan ^{-1}\left(\frac{\sqrt{c}}{\sqrt{b}}\right)}} \end{eqnarray}

In the top line we substituted for $u=\sin(\theta) \sqrt{d/(b+c \sin(\theta)^2)}$. In the second line we substituted $u = \sqrt{d/(c+b)} \sin(\phi)$. In the third line we substituted $t=\tan(\phi/2)$. In the forth line we used partial fraction decomposition and properties of the logarithm. Finally in the fifth line we used the anti-derivative defined in $(ii)$.

Clear[F]; Clear[FF];
F[z_, a_, b_] := 
  Log[a + z] Log[(b + z)/(-a + b)] + PolyLog[2, (a + z)/(a - b)];
FF[A_, B_, a_, b_] := 
  Module[{result, ts, zs, zsp, zsm, eps = 10^(-15)},
   (*This is Integrate[Log[z+a]/(z+b),{z,A,B}] where all a,b,A, 
   and B are complex. *)
   result = F[B, a, b] - F[A, a, b];


   ts = - (Im[(A + b) (Conjugate[b] - Conjugate[a])]/
     Im[(B - A) (Conjugate[b] - Conjugate[a])]);
   If[0 <= ts <= 1,
    zsp = A + (ts + eps) (B - A);
    zsm = A + (ts - eps) (B - A);
    result += -F[zsp, a, b] + F[zsm, a, b];
    ];

   result
   ];

{a, b, c, d} = RandomReal[{0, 3}, 4, WorkingPrecision -> 50];
NIntegrate[
 Exp[-a/2 x^2 - b/2 y^2 - c/2 z^2 - d/2 w^2], {x, 0, Infinity}, {y, 0,
   x}, {z, 0, y}, {w, 0, z}]
NIntegrate[
 Sin[th]/Sqrt[a b d (b + c Sin[th]^2)] ArcTan[
   Sin[th] Sqrt[d/(b + c Sin[th]^2)]], {th, 0, ArcTan[Sqrt[b/a]]}]
  1/Sqrt[a ] NIntegrate[
  u ArcTan[u] 1/((d - c u^2) Sqrt[d - (c + b) u^2]), {u, 0, Sqrt[ 
   d/ (a + b + c)]}]
 Sqrt[d]/(Sqrt[a ] (b + c))
  NIntegrate[
  Sin[phi] ArcTan[
    Sqrt[d/(c + b)] Sin[phi]] 1/(d - c (d/(c + b) Sin[phi]^2)) , {phi,
    0, ArcSin[ Sqrt[( c + b)/ (a + b + c)]]}]
- I 2/(Sqrt[a ] Sqrt[d])
  NIntegrate[
   t /(c (-1 + t^2)^2 + b (1 + t^2)^2) Log[(
    1 + t^2 + 2 I Sqrt[d/(b + c)] t)/(
    1 + t^2 - 2 I Sqrt[d/(b + c)] t)], {t, 0, Sqrt[(b + c)/(
   a + b + c)]/(1 + Sqrt[a/(a + b + c)])}]
- I 2/(Sqrt[a ] Sqrt[d])
  NIntegrate[
   t /(c (-1 + t^2)^2 + 
     b (1 + t^2)^2) Log[((1/
        2 (2 I Sqrt[d/(b + c)] - Sqrt[-4 - (4 d)/(b + c)]) + 
       t) (1/2 (2 I Sqrt[d/(b + c)] + Sqrt[-4 - (4 d)/(b + c)]) + 
       t))/((1/2 (-2 I Sqrt[d/(b + c)] - Sqrt[-4 - (4 d)/(b + c)]) + 
       t) (1/2 (-2 I Sqrt[d/(b + c)] + Sqrt[-4 - (4 d)/(b + c)]) + 
       t))], {t, 0, Sqrt[(b + c)/(a + b + c)]/(
   1 + Sqrt[a/(a + b + c)])}]
 1/Sqrt[a b c d]  1/4 NIntegrate[
  Sum[(-1)^(Floor[(eta - 1)/2]) (-1)^
    Floor[(xi - 1)/2] Log[
     t + (-1)^Floor[(xi - 1)/2] I Sqrt[d/(b + c)] + (-1)^(xi - 1)
        I Sqrt[( b + c + d)/(b + c)]]/(
    t - (-1)^(1 + eta + 
        Floor[(eta - 1)/2]) I Exp[(-1)^(Floor[(eta - 1)/2]) I ArcTan[
         Sqrt[c]/Sqrt[b]]]), {xi, 1, 4}, {eta, 1, 4}], {t, 0, Sqrt[(
   b + c)/(a + b + c)]/(1 + Sqrt[a/(a + b + c)])}]
 1/Sqrt[a b c d]  1/4 Sum[(-1)^(Floor[(eta - 1)/2]) (-1)^
   Floor[(xi - 1)/2] FF[0, Sqrt[(b + c)/(a + b + c)]/(
    1 + Sqrt[a/(
     a + b + c)]), (-1)^Floor[(xi - 1)/2] I Sqrt[d/(b + c)] + (-1)^(
      xi - 1) I Sqrt[( b + c + d)/(
      b + c)], -(-1)^(1 + eta + 
        Floor[(eta - 1)/2]) I Exp[(-1)^(Floor[(eta - 1)/2]) I ArcTan[
        Sqrt[c]/Sqrt[b]]]], {xi, 1, 4}, {eta, 1, 4}]

Update: As a sanity check look at the case $a=b=c=d=1$. Then define: \begin{eqnarray} M1&:=&\left( \begin{array}{cccc} -1+\sqrt{3} & \sqrt{2} & \sqrt{2} & -1+\sqrt{3} \\ 1 & \sqrt{2-\sqrt{3}} & \sqrt{2-\sqrt{3}} & 1 \\ \sqrt{2-\sqrt{3}} & 1 & 1 & \sqrt{2-\sqrt{3}} \\ \sqrt{2} & -1+\sqrt{3} & -1+\sqrt{3} & \sqrt{2} \\ \end{array} \right)\\ M2&:=&\left( \begin{array}{cccc} \frac{1}{\sqrt{2}} & \frac{1}{2} \left(1+\sqrt{3}\right) & \frac{1}{2} \left(1+\sqrt{3}\right) & \frac{1}{\sqrt{2}} \\ \frac{1}{\sqrt{2}} & \frac{1}{2} \left(-1+\sqrt{3}\right) & \frac{1}{2} \left(-1+\sqrt{3}\right) & \frac{1}{\sqrt{2}} \\ \frac{1}{2} \left(-1+\sqrt{3}\right) & \frac{1}{\sqrt{2}} & \frac{1}{\sqrt{2}} & \frac{1}{2} \left(-1+\sqrt{3}\right) \\ \frac{1}{2} \left(1+\sqrt{3}\right) & \frac{1}{\sqrt{2}} & \frac{1}{\sqrt{2}} & \frac{1}{2} \left(1+\sqrt{3}\right) \\ \end{array} \right)\\ A1&:=&\left( \begin{array}{cccc} -\frac{\pi }{6} & \frac{\pi }{12} & -\frac{\pi }{4} & 0 \\ \frac{5 \pi }{6} & \frac{\pi }{12} & \frac{5 \pi }{12} & -\frac{\pi }{3} \\ -\frac{5 \pi }{12} & \frac{\pi }{3} & -\frac{5 \pi }{6} & -\frac{\pi }{12} \\ \frac{\pi }{4} & 0 & \frac{\pi }{6} & -\frac{\pi }{12} \\ \end{array} \right)\\ A2&:=&\left( \begin{array}{cccc} -\frac{\pi }{12} & \frac{\pi }{6} & -\frac{\pi }{6} & \frac{\pi }{12} \\ \frac{7 \pi }{12} & -\frac{\pi }{6} & \frac{\pi }{6} & -\frac{7 \pi }{12} \\ -\frac{\pi }{6} & \frac{7 \pi }{12} & -\frac{7 \pi }{12} & \frac{\pi }{6} \\ \frac{\pi }{6} & -\frac{\pi }{12} & \frac{\pi }{12} & -\frac{\pi }{6} \\ \end{array} \right) \end{eqnarray} and we have \begin{eqnarray} I(1,1,1,1)=\frac{1}{4} \sum\limits_{\xi=1}^4 \sum\limits_{\eta=1}^4 (-1)^{\left\lfloor \frac{\eta -1}{2}\right\rfloor +\left\lfloor \frac{\xi -1}{2}\right\rfloor } \left( Li_2(M1_{\xi,\eta}\exp(\imath A1_{\xi,\eta}))- Li_2(M2_{\xi,\eta}\exp(\imath A2_{\xi,\eta})) \right) \end{eqnarray} We have checked numerically that this quantity above coincides with $\pi^2/96$ to one hundred digits. It would be interesting to prove this analytically.

Here we give an answer using a different method. Assume that $a\ge0$, $b\ge 0$, $c\ge 0$ and $d\ge 0$. Define: \begin{equation} I(a,b,c,d):=\int _0^{\infty }\int _0^x\int _0^y\int _0^z \exp \left( -\frac{a x^2}{2}-\frac{b y^2}{2}-\frac{c z^2}{2}-\frac{d w^2}{2} \right) \,\mathrm{d}w\,\mathrm{d}z\,\mathrm{d}y\,\mathrm{d}x \end{equation} Then we have: \begin{eqnarray} &&I(a,b,c,d)=\\ &&\sqrt{\frac{\pi}{2 c d}} \int\limits_0^\infty \int\limits_0^x \int\limits_0^{\sqrt{c} y}\exp(-\frac{a x^2}{2} - \frac{b y^2}{2}-\frac{ z^2}{2})\cdot \text{erf}(\sqrt{\frac{d}{2 c}} z)\; dz \; dy \;dx=\\ &&\frac{2 \pi}{\sqrt{c d}} \int\limits_0^\infty \int\limits_0^x \exp(-\frac{a x^2}{2}-\frac{b y^2}{2}) \left(\frac{1}{2\pi} \arctan[\sqrt{\frac{d}{c}}] - T(\sqrt{c} y,\sqrt{\frac{d}{c}}\right)\; dy \; dx=\\ && \frac{2 \pi^2}{\sqrt{a b c d}}\int\limits_0^\infty \frac{\exp(-\frac{y^2}{2})}{\sqrt{2\pi}} \text{erfc}(\sqrt{\frac{a}{2 b}y})\left(\frac{1}{2\pi} \arctan[\sqrt{\frac{d}{c}}] - T(\sqrt{\frac{c}{b}} y,\sqrt{\frac{d}{c}}\right) \; dy=\\ && \frac{2 \pi^2}{\sqrt{a b c d}} \left( \frac{\left(\pi -2 \tan ^{-1}\left(\sqrt{\frac{a}{b}}\right)\right) \tan ^{-1}\left(\sqrt{\frac{d}{c}}\right)-\pi \sin ^{-1}\left(\frac{\sqrt{b d}}{\sqrt{(b+c) (c+d)}}\right)}{4 \pi ^2}+ \int\limits_0^\infty \frac{\exp(-\frac{y^2}{2})}{\sqrt{2 \pi}} \text{erf}(\sqrt{\frac{a}{2 b}} y) T(\sqrt{\frac{c}{b}} y,\sqrt{\frac{d}{c}} )\; dy \right)=\\ && \frac{2 \pi^2}{\sqrt{a b c d}} \left( \frac{\tan ^{-1}\left(\sqrt{\frac{d}{c}}\right)-\tan ^{-1}\left(\frac{\sqrt{b d}}{\sqrt{c (b+c+d)}}\right)}{4 \pi } +\right.\\ &&\left. \frac{1}{8 \pi^2} \sum\limits_{i=1}^4 \sum\limits_{j=1}^4 (-1)^{\left\lfloor \frac{i-1}{2}\right\rfloor +j} {\mathfrak F}^{(1,\sqrt{\frac{a+b+c}{b+c}}-\sqrt{\frac{a}{b+c}})}_{i \left((-1)^{\left\lfloor \frac{j-1}{2}\right\rfloor } \sqrt{\frac{b+c+d}{b+c}}+(-1)^j \sqrt{\frac{d}{b+c}}\right),(-1)^{i+1} \sqrt{\frac{c}{b+c}}+i \sqrt{\frac{b}{b+c}} (-1)^{\left\lceil \frac{i-1}{2}\right\rceil +1}} \right) \end{eqnarray} In the first line we integrated over $w$ using the definition of the error function. In the second line we integrated over $z$ using the definition of the Owen's T function. In the third line we swapped the order of integration and integrated over $x$ using the definition of the complementary error function. In the forth line we split the integral into doable integrals and a more complicated one and finally in the fifth line we evaluated the remaining integral using An integral involving a Gaussian, error functions and the Owen's T function. .

Clear[F]; Clear[FF];
F[z_, a_, b_] := 
  Log[a + z] Log[(b + z)/(-a + b)] + PolyLog[2, (a + z)/(a - b)];
FF[A_, B_, a_, b_] := 
  Module[{result, ts, zs, zsp, zsm, eps = 10^(-50)},
   (*This is Integrate[Log[z+a]/(z+b),{z,A,B}] where all a,b,A, 
   and B are complex. *)
   result = F[B, a, b] - F[A, a, b];


   ts = - (Im[(A + b) (Conjugate[b] - Conjugate[a])]/
     Im[(B - A) (Conjugate[b] - Conjugate[a])]);
   If[0 <= ts <= 1,
    zsp = A + (ts + eps) (B - A);
    zsm = A + (ts - eps) (B - A);
    result += -F[zsp, a, b] + F[zsm, a, b];
    ];

   result
   ];
J[a_, b_, c_] := 
  1/ Pi^2 (ArcTan[Sqrt[2] a]/2 ArcTan[ c] + 
     1/8  Sum[
       FF[1, ( Sqrt[1 + 2 a^2 + b^2] - Sqrt[2] a)/Sqrt[
         1 + b^2], ((-1)^j I b c + (-1)^Floor[(j - 1)/2] I Sqrt[
           1 + b^2 + b^2 c^2])/Sqrt[
         1 + b^2], -(((-1)^Ceiling[(i - 1)/2] I + (-1)^i b)/Sqrt[
          1 + b^2])] (-1)^(j + Floor[(i - 1)/2]), {i, 1, 4}, {j, 1, 
        4}] );

{a, b, c, d} = RandomReal[{0, 10}, 4, WorkingPrecision -> 50];
NIntegrate[
 Exp[-a/2 x^2 - b/2 y^2 - c/2 z^2 - d/2 w^2], {x, 0, Infinity}, {y, 0,
   x}, {z, 0, y}, {w, 0, z}]
Sqrt[\[Pi]/2]/(Sqrt[c] Sqrt[d])
  NIntegrate[
  Exp[-a/2 x^2 - b/2 y^2 - 1/2 z^2] Erf[Sqrt[d/(2 c)] z], {x, 0, 
   Infinity}, {y, 0, x}, {z, 0, Sqrt[c] y}]
 (2 \[Pi])/Sqrt[c d]
  NIntegrate[
  Exp[-a/2 x^2 - b/2 y^2] (ArcTan[Sqrt[d/c]]/(2 \[Pi]) - 
     OwenT[Sqrt[c] y, Sqrt[d/c]]), {x, 0, Infinity}, {y, 0, x}]
(2  \[Pi]^2)/Sqrt[a b c d]
  NIntegrate[
  Exp[ -1/2 y^2]/Sqrt[2 Pi]
    Erfc[Sqrt[a/(2 b)] y] (ArcTan[Sqrt[d/c]]/(2 \[Pi]) - 
     OwenT[Sqrt[c/b] y, Sqrt[d/c]]), {y, 0, Infinity}]
(2  \[Pi]^2)/Sqrt[
 a b c d] ((-\[Pi] ArcSin[Sqrt[b d]/
      Sqrt[(b + c) (c + d)]] + (\[Pi] - 2 ArcTan[Sqrt[a/b]]) ArcTan[
      Sqrt[d/c]])/(4 \[Pi]^2) + 
   NIntegrate[
    Exp[ -1/2 y^2]/Sqrt[2 Pi]
      Erf[Sqrt[a/(2 b)] y] OwenT[Sqrt[c/b] y, Sqrt[d/c]], {y, 0, 
     Infinity}])
(2  \[Pi]^2)/Sqrt[
 a b c d] ((-ArcTan[Sqrt[b d]/ Sqrt[c (b + c + d)]] + 
    ArcTan[Sqrt[d/c]])/(4 \[Pi]) +  
   1/( 8 Pi^2)
     Sum[FF[1, Sqrt[(a + b + c)/(b + c)] - Sqrt[a/(b + c)], 
       I  ((-1)^j Sqrt[d/(b + c)] + (-1)^Floor[1/2 (-1 + j)] Sqrt[(
           b + c + d)/(b + c)]), (-1)^(i + 1) Sqrt[c/(b + c)] + 
        Sqrt[b/(b + c)] I (-1)^(1 + Ceiling[1/2 (-1 + i)])] (-1)^(
      j + Floor[(i - 1)/2]), {i, 1, 4}, {j, 1, 4}] )