# Find closed form for quadruple integral

I am trying to find a closed form of the following integral $$\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$$ where $$a,b,c,d>0$$ are some constants.

My idea is to change variable to to polar system by letting $$x=\frac{r \cos (\alpha ) \cos (\beta ) \cos (\theta )}{\sqrt{a}}, \quad y=\frac{r \cos (\alpha ) \cos (\beta ) \sin (\theta )}{\sqrt{a}}$$ $$z=\frac{r \sin (\alpha ) \cos (\beta )}{\sqrt{c}}, \quad w=\frac{r \sin (\beta )}{\sqrt{d}}$$

This reduces the original integral into $$\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$$ But then I get stuck here.

PS: I am interested in this because I found that $$\int _0^{\infty }\int _0^x \exp \left(-\frac{a x^2}{2}-\frac{b y^2}{2}\right)\,\mathrm{d}y\,\mathrm{d}x = \frac{\tan ^{-1}\left(\sqrt{\frac{b}{a}}\right)}{\sqrt{a b}}$$ and $$\int _0^{\infty }\int _0^x\int _0^y \exp \left( -\frac{a x^2}{2}-\frac{b y^2}{2}-\frac{c z^2}{2} \right) \,\mathrm{d}z\,\mathrm{d}y\,\mathrm{d}x = \frac{\sqrt{\pi/2 } }{\sqrt{a b c}} \left(\tan ^{-1}\left(\sqrt{\frac{c}{b}}\right)-\tan ^{-1}\left(\sqrt{\frac{a c}{b (a+b+c)}}\right)\right)$$ So I am trying to generalize this. Maybe this is already known?

• At least for all constants equal to $1$ Wolfram Alpha claims that the answer is $\tfrac {\pi^2}{96}$. Apr 28, 2019 at 15:29
• Well from the pattern I deduce that there must be a $\sf{\dfrac{\tan^{-1}\left(\sqrt{\dfrac dc}\right)}{\sqrt{abcd}}}$ term... Apr 29, 2019 at 6:26

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 $$$$F[z,a,b] := \log(z+a) \log\left( \frac{z+b}{b-a}\right) + Li_2\left( \frac{z+a}{a-b}\right)$$$$ 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.

• This is an impressive amount of algebra, but looking at Your Mathematica code at the end we're still left with one integral to be computed. It's unclear to me how is this more efficient than the original integral. What were we going for here? Apr 30, 2019 at 16:03
• @Radost Sorry, the code was outdated. Now the code uses no numerical integrals. You can easily get the result with fifty digits of accuracy as I do above. But going back to your point I have played around with integrals like that many times. They simply reduce to poly-logarithms. there is no way around that. Yet polylogarithms can be computed with arbitrary precision whereas your original integral cannot. Apr 30, 2019 at 16:22

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.

Here we give an answer using a different method. Assume that $$a\ge0$$, $$b\ge 0$$, $$c\ge 0$$ and $$d\ge 0$$. Define: $$$$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$$$$ 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}] )


Define the function $$\mathcal{I}:\mathbb{R}_{>0}^{4}\rightarrow\mathbb{R}$$ via the quadruple integral

$$\mathcal{I}{\left(a,b,c,d\right)}:=\int_{0}^{\infty}\mathrm{d}x\int_{0}^{x}\mathrm{d}y\int_{0}^{y}\mathrm{d}z\int_{0}^{z}\mathrm{d}w\,\exp{\left(-\frac{ax^{2}}{2}-\frac{by^{2}}{2}-\frac{cz^{2}}{2}-\frac{dw^{2}}{2}\right)}.\tag{1}$$

If we look at how this integral transforms under the substitution $$\left(x,y,z,w\right)\mapsto\left(px,py,pz,pw\right)$$ for some fixed but arbitrary positive real $$p$$, we obtain the following scaling relation:

$$\forall\left(a,b,c,d,p\right)\in\mathbb{R}_{>0}^{5}:\mathcal{I}{\left(a,b,c,d\right)}=p^{4}\mathcal{I}{\left(ap^{2},bp^{2},cp^{2},dp^{2}\right)}.\tag{2}$$

As such, in our general evaluation of $$\mathcal{I}{\left(a,b,c,d\right)}$$ it will suffice to consider the $$a=1$$ case.

To begin with, we derive some quick integration formulas that will be helpful below.

For any $$p\in\mathbb{R}_{>0}$$,

\begin{align} \int_{0}^{\infty}\mathrm{d}x\,x^{3}\exp{\left(-\frac{px^{2}}{2}\right)} &=\int_{0}^{\infty}\mathrm{d}y\,\frac{y}{2}\exp{\left(-\frac{py}{2}\right)};~~~\small{\left[x=\sqrt{y}\right]}\\ &=\int_{0}^{\infty}\mathrm{d}z\,\frac{2}{p}\cdot\frac{z}{p}\exp{\left(-z\right)};~~~\small{\left[y=\frac{2z}{p}\right]}\\ &=\frac{2}{p^{2}}\int_{0}^{\infty}\mathrm{d}z\,z\exp{\left(-z\right)}\\ &=\frac{2}{p^{2}}.\tag{3a}\\ \end{align}

Next, given $$p\in\mathbb{R}_{>0}$$ and setting $$q:=\sqrt{p}$$,

\begin{align} \int_{0}^{1}\mathrm{d}t\,\frac{2t^{2}}{\left(1+pt^{2}\right)^{2}} &=\int_{0}^{1}\mathrm{d}t\,\frac{2t^{2}}{\left(1+q^{2}t^{2}\right)^{2}}\\ &=\frac{1}{q^{3}}\int_{0}^{q}\mathrm{d}x\,\frac{2x^{2}}{\left(1+x^{2}\right)^{2}};~~~\small{\left[qt=x\right]}\\ &=\frac{1}{q^{3}}\int_{0}^{q}\mathrm{d}x\,\frac{\mathrm{d}}{\mathrm{d}x}\left[\arctan{\left(x\right)}-\frac{x}{1+x^{2}}\right]\\ &=\frac{1}{q^{3}}\left[\arctan{\left(q\right)}-\frac{q}{1+q^{2}}\right]\\ &=\left[\frac{\arctan{\left(q\right)}}{q^{3}}-\frac{1}{q^{2}\left(1+q^{2}\right)}\right]\\ &=\left[\frac{\arctan{\left(\sqrt{p}\right)}}{p\sqrt{p}}-\frac{1}{p\left(1+p\right)}\right].\tag{3b}\\ \end{align}

Defining the function $$\mathcal{J}:\mathbb{R}_{>0}^{3}\rightarrow\mathbb{R}$$ by

$$\mathcal{J}{\left(p,q,x\right)}:=\int_{0}^{x}\mathrm{d}y\,\left[\frac{\arctan{\left(\sqrt{p+qy^{2}}\right)}}{\left(p+qy^{2}\right)^{3/2}}-\frac{1}{\left(p+qy^{2}\right)\left(1+p+qy^{2}\right)}\right],\tag{3c}$$

we then find that for any $$\left(p,q,x\right)\in\mathbb{R}_{>0}^{3}$$,

\begin{align} \mathcal{J}{\left(p,q,x\right)} &=\int_{0}^{x}\mathrm{d}y\,\left[\frac{\arctan{\left(\sqrt{p+qy^{2}}\right)}}{\left(p+qy^{2}\right)^{3/2}}-\frac{1}{\left(p+qy^{2}\right)\left(1+p+qy^{2}\right)}\right]\\ &=\int_{0}^{x}\mathrm{d}y\,\frac{\arctan{\left(\sqrt{p+qy^{2}}\right)}}{\left(p+qy^{2}\right)^{3/2}}-\int_{0}^{x}\mathrm{d}y\,\frac{1}{\left(p+qy^{2}\right)\left(1+p+qy^{2}\right)}\\ &=\frac{x\arctan{\left(\sqrt{p+qx^{2}}\right)}}{p\sqrt{p+qx^{2}}}-\int_{0}^{x}\mathrm{d}y\,\frac{qy^{2}}{p\left(p+qy^{2}\right)\left(1+p+qy^{2}\right)};~~~\small{I.B.P.}\\ &~~~~~-\int_{0}^{x}\mathrm{d}y\,\frac{1}{\left(p+qy^{2}\right)\left(1+p+qy^{2}\right)}\\ &=\frac{x\arctan{\left(\sqrt{p+qx^{2}}\right)}}{p\sqrt{p+qx^{2}}}-\int_{0}^{x}\mathrm{d}y\,\frac{1}{p\left(1+p+qy^{2}\right)}\\ &=\frac{x\arctan{\left(\sqrt{p+qx^{2}}\right)}}{p\sqrt{p+qx^{2}}}-\frac{1}{p\sqrt{q\left(1+p\right)}}\int_{0}^{x\sqrt{\frac{q}{1+p}}}\mathrm{d}t\,\frac{1}{\left(1+t^{2}\right)};~~~\small{\left[y=t\sqrt{\frac{1+p}{q}}\right]}\\ &=\frac{x\arctan{\left(\sqrt{p+qx^{2}}\right)}}{p\sqrt{p+qx^{2}}}-\frac{\arctan{\left(\frac{x\sqrt{q}}{\sqrt{1+p}}\right)}}{p\sqrt{q\left(1+p\right)}}.\tag{3d}\\ \end{align}

Suppose $$\left(b,c,d\right)\in\mathbb{R}_{>0}^{3}$$.

We begin by reducing $$\mathcal{I}$$ to a single-variable integral in the following way:

\begin{align} \mathcal{I}{\left(1,b,c,d\right)} &=\int_{0}^{\infty}\mathrm{d}x\int_{0}^{x}\mathrm{d}y\int_{0}^{y}\mathrm{d}z\int_{0}^{z}\mathrm{d}w\,\exp{\left(-\frac{x^{2}}{2}-\frac{by^{2}}{2}-\frac{cz^{2}}{2}-\frac{dw^{2}}{2}\right)}\\ &=\int_{0}^{\infty}\mathrm{d}x\int_{0}^{1}\mathrm{d}t\int_{0}^{t}\mathrm{d}u\int_{0}^{u}\mathrm{d}v\,x^{3}\\ &~~~~~\times\exp{\left(-\frac{x^{2}}{2}-\frac{bx^{2}t^{2}}{2}-\frac{cx^{2}u^{2}}{2}-\frac{dx^{2}v^{2}}{2}\right)};~~~\small{\left[\left(y,z,w\right)=\left(xt,xu,xv\right)\right]}\\ &=\int_{0}^{\infty}\mathrm{d}x\int_{0}^{1}\mathrm{d}t\int_{0}^{t}\mathrm{d}u\int_{0}^{u}\mathrm{d}v\,x^{3}e^{-\frac{\left(1+bt^{2}+cu^{2}+dv^{2}\right)x^{2}}{2}}\\ &=\int_{0}^{1}\mathrm{d}t\int_{0}^{t}\mathrm{d}u\int_{0}^{u}\mathrm{d}v\int_{0}^{\infty}\mathrm{d}x\,x^{3}e^{-\frac{\left(1+bt^{2}+cu^{2}+dv^{2}\right)x^{2}}{2}}\\ &=\int_{0}^{1}\mathrm{d}t\int_{0}^{t}\mathrm{d}u\int_{0}^{u}\mathrm{d}v\,\frac{2}{\left(1+bt^{2}+cu^{2}+dv^{2}\right)^{2}},\\ \end{align}

and then,

\begin{align} \mathcal{I}{\left(1,b,c,d\right)} &=\int_{0}^{1}\mathrm{d}t\int_{0}^{t}\mathrm{d}u\int_{0}^{u}\mathrm{d}v\,\frac{2}{\left(1+bt^{2}+cu^{2}+dv^{2}\right)^{2}}\\ &=\int_{0}^{1}\mathrm{d}t\int_{0}^{1}\mathrm{d}x\int_{0}^{x}\mathrm{d}y\,\frac{2t^{2}}{\left(1+bt^{2}+ct^{2}x^{2}+dt^{2}y^{2}\right)^{2}};~~~\small{\left[\left(u,v\right)=\left(tx,ty\right)\right]}\\ &=\int_{0}^{1}\mathrm{d}x\int_{0}^{x}\mathrm{d}y\int_{0}^{1}\mathrm{d}t\,\frac{2t^{2}}{\left(1+bt^{2}+ct^{2}x^{2}+dt^{2}y^{2}\right)^{2}}\\ &=\int_{0}^{1}\mathrm{d}x\int_{0}^{x}\mathrm{d}y\int_{0}^{1}\mathrm{d}t\,\frac{2t^{2}}{\left[1+\left(b+cx^{2}+dy^{2}\right)t^{2}\right]^{2}}\\ &=\int_{0}^{1}\mathrm{d}x\int_{0}^{x}\mathrm{d}y\,\bigg{[}\frac{\arctan{\left(\sqrt{b+cx^{2}+dy^{2}}\right)}}{\left(b+cx^{2}+dy^{2}\right)^{3/2}}\\ &~~~~~-\frac{1}{\left(b+cx^{2}+dy^{2}\right)\left(1+b+cx^{2}+dy^{2}\right)}\bigg{]}\\ &=\int_{0}^{1}\mathrm{d}x\,\mathcal{J}{\left(b+cx^{2},d,x\right)}\\ &=\int_{0}^{1}\mathrm{d}x\,\left[\frac{x\arctan{\left(\sqrt{b+cx^{2}+dx^{2}}\right)}}{\left(b+cx^{2}\right)\sqrt{b+cx^{2}+dx^{2}}}-\frac{\arctan{\left(\frac{x\sqrt{d}}{\sqrt{1+b+cx^{2}}}\right)}}{\left(b+cx^{2}\right)\sqrt{d\left(1+b+cx^{2}\right)}}\right].\tag{4}\\ \end{align}

Define the auxiliary functions $$\mathcal{F}:\mathbb{R}_{>0}^{3}\rightarrow\mathbb{R}$$ and $$\mathcal{G}:\mathbb{R}_{>0}^{3}\rightarrow\mathbb{R}$$ by the respective integrals

\begin{align} \mathcal{F}{\left(b,c,d\right)} &:=\int_{0}^{1}\mathrm{d}x\,\frac{x\arctan{\left(\sqrt{b+cx^{2}+dx^{2}}\right)}}{\left(b+cx^{2}\right)\sqrt{b+cx^{2}+dx^{2}}}\tag{5a}\\ \end{align}

and

\begin{align} \mathcal{G}{\left(b,c,d\right)} &:=\int_{0}^{1}\mathrm{d}x\,\frac{\arctan{\left(\frac{x\sqrt{d}}{\sqrt{1+b+cx^{2}}}\right)}}{\left(b+cx^{2}\right)\sqrt{d\left(1+b+cx^{2}\right)}}.\tag{5b}\\ \end{align}

Suppose $$\left(b,c,d\right)\in\mathbb{R}_{>0}^{3}$$, and set $$p:=\sqrt{\frac{b}{c}}\land q:=\sqrt{\frac{1+b}{c}}\land r:=\sqrt{\frac{d}{c}}$$. Next, noting that $$0, set $$s:=\frac{p}{q}\land z:=q^{-1}$$. We then obtain the following expressions for $$\mathcal{G}$$ and $$\mathcal{F}$$:

\begin{align} \mathcal{G}{\left(b,c,d\right)} &=\int_{0}^{1}\mathrm{d}x\,\frac{\arctan{\left(\frac{x\sqrt{d}}{\sqrt{1+b+cx^{2}}}\right)}}{\left(b+cx^{2}\right)\sqrt{d\left(1+b+cx^{2}\right)}}\\ &=\frac{1}{c\sqrt{cd}}\int_{0}^{1}\mathrm{d}x\,\frac{\arctan{\left(\frac{rx}{\sqrt{q^{2}+x^{2}}}\right)}}{\left(p^{2}+x^{2}\right)\sqrt{q^{2}+x^{2}}}\\ &=\frac{1}{c\sqrt{cd}}\int_{0}^{q^{-1}}\mathrm{d}y\,\frac{q\arctan{\left(\frac{rqy}{\sqrt{q^{2}+q^{2}y^{2}}}\right)}}{\left(p^{2}+q^{2}y^{2}\right)\sqrt{q^{2}+q^{2}y^{2}}};~~~\small{\left[x=qy\right]}\\ &=\frac{1}{c\sqrt{cd}}\int_{0}^{q^{-1}}\mathrm{d}y\,\frac{\arctan{\left(\frac{ry}{\sqrt{1+y^{2}}}\right)}}{\left(p^{2}+q^{2}y^{2}\right)\sqrt{1+y^{2}}}\\ &=\frac{1}{\left(1+b\right)\sqrt{cd}}\int_{0}^{z}\mathrm{d}y\,\frac{\arctan{\left(\frac{ry}{\sqrt{1+y^{2}}}\right)}}{\left(s^{2}+y^{2}\right)\sqrt{1+y^{2}}}\\ &=\frac{1}{\left(1+b\right)\sqrt{cd}}\int_{0}^{\frac{z}{\sqrt{1+z^{2}}}}\mathrm{d}t\,\frac{\arctan{\left(rt\right)}}{s^{2}+\left(1-s^{2}\right)t^{2}};~~~\small{\left[\frac{y}{\sqrt{1+y^{2}}}=t\right]}\\ &=\frac{1}{\left(1+b\right)\sqrt{cd}}\int_{0}^{\frac{rz}{\sqrt{1+z^{2}}}}\mathrm{d}u\,\frac{r\arctan{\left(u\right)}}{r^{2}s^{2}+\left(1-s^{2}\right)u^{2}};~~~\small{\left[rt=u\right]}\\ &=\int_{0}^{\sqrt{\frac{d}{1+b+c}}}\mathrm{d}u\,\frac{\arctan{\left(u\right)}}{bd+cu^{2}},\\ \end{align}

and

\begin{align} \mathcal{F}{\left(b,c,d\right)} &=\int_{0}^{1}\mathrm{d}x\,\frac{x\arctan{\left(\sqrt{b+cx^{2}+dx^{2}}\right)}}{\left(b+cx^{2}\right)\sqrt{b+cx^{2}+dx^{2}}}\\ &=\int_{0}^{1}\mathrm{d}x\,\frac{2x\arctan{\left(\sqrt{b+\left(c+d\right)x^{2}}\right)}}{2\left(b+cx^{2}\right)\sqrt{b+\left(c+d\right)x^{2}}}\\ &=\int_{0}^{1}\mathrm{d}y\,\frac{\arctan{\left(\sqrt{b+\left(c+d\right)y}\right)}}{2\left(b+cy\right)\sqrt{b+\left(c+d\right)y}};~~~\small{\left[x^{2}=y\right]}\\ &=\int_{b}^{b+c+d}\mathrm{d}t\,\frac{\arctan{\left(\sqrt{t}\right)}}{2\left(bd+ct\right)\sqrt{t}};~~~\small{\left[y=\frac{t-b}{c+d}\right]}\\ &=\int_{\sqrt{b}}^{\sqrt{b+c+d}}\mathrm{d}u\,\frac{\arctan{\left(u\right)}}{bd+cu^{2}};~~~\small{\left[\sqrt{t}=u\right]}.\\ \end{align}

Hence, we can express $$\mathcal{I}$$ as

\begin{align} \mathcal{I}{\left(1,b,c,d\right)} &=\int_{0}^{1}\mathrm{d}x\,\left[\frac{x\arctan{\left(\sqrt{b+cx^{2}+dx^{2}}\right)}}{\left(b+cx^{2}\right)\sqrt{b+cx^{2}+dx^{2}}}-\frac{\arctan{\left(\frac{x\sqrt{d}}{\sqrt{1+b+cx^{2}}}\right)}}{\left(b+cx^{2}\right)\sqrt{d\left(1+b+cx^{2}\right)}}\right]\\ &=\mathcal{F}{\left(b,c,d\right)}-\mathcal{G}{\left(b,c,d\right)}\\ &=\int_{\sqrt{b}}^{\sqrt{b+c+d}}\mathrm{d}u\,\frac{\arctan{\left(u\right)}}{bd+cu^{2}}-\int_{0}^{\sqrt{\frac{d}{1+b+c}}}\mathrm{d}u\,\frac{\arctan{\left(u\right)}}{bd+cu^{2}}\\ &=\frac{1}{\sqrt{bcd}}\int_{\sqrt{\frac{c}{d}}}^{\sqrt{\frac{(b+c+d)c}{bd}}}\mathrm{d}x\,\frac{\arctan{\left(x\sqrt{\frac{bd}{c}}\right)}}{1+x^{2}}\\ &~~~~~-\frac{1}{\sqrt{bcd}}\int_{0}^{\sqrt{\frac{c}{b(1+b+c)}}}\mathrm{d}x\,\frac{\arctan{\left(x\sqrt{\frac{bd}{c}}\right)}}{1+x^{2}};~~~\small{\left[u=x\sqrt{\frac{bd}{c}}\right]}\\ &=\frac{1}{\sqrt{bcd}}\int_{0}^{\sqrt{\frac{(b+c+d)c}{bd}}}\mathrm{d}x\,\frac{\arctan{\left(x\sqrt{\frac{bd}{c}}\right)}}{1+x^{2}}-\frac{1}{\sqrt{bcd}}\int_{0}^{\sqrt{\frac{c}{d}}}\mathrm{d}x\,\frac{\arctan{\left(x\sqrt{\frac{bd}{c}}\right)}}{1+x^{2}}\\ &~~~~~-\frac{1}{\sqrt{bcd}}\int_{0}^{\sqrt{\frac{c}{b(1+b+c)}}}\mathrm{d}x\,\frac{\arctan{\left(x\sqrt{\frac{bd}{c}}\right)}}{1+x^{2}}.\tag{6}\\ \end{align}

Finally, it can be shown (see Appendix) that the following integration formula holds for all $$\left(p,z\right)\in\mathbb{R}_{>0}^{2}$$:

$$\int_{0}^{z}\mathrm{d}x\,\frac{2\arctan{\left(px\right)}}{1+x^{2}}=\arctan^{2}{\left(z\right)}+\operatorname{Li}_{2}{\left(-\frac{1-p}{1+p}\right)}-\operatorname{Li}_{2}{\left(\frac{1-p}{1+p},\pi-2\arctan{\left(z\right)}\right)},\tag{7}$$

where the two-variable variant of the dilogarithm is defined by the integral representation

$$\operatorname{Li}_{2}{\left(r,\theta\right)}:=-\frac12\int_{0}^{r}\mathrm{d}x\,\frac{\ln{\left(1-2x\cos{\left(\theta\right)}+x^{2}\right)}}{x};~~~\small{\left(r,\theta\right)\in\mathbb{R}^{2}}.$$

Since each of the three remaining integrals in the last line of $$(6)$$ above can be evaluated with formula $$(7)$$, this in principle completes the derivation. I don't see much point in going through the tedium of actually writing out the explicit expression though.

Appendix.

Define the function $$\mathcal{K}:\mathbb{R}_{>0}^{2}\rightarrow\mathbb{R}$$ via the definite integral

$$\mathcal{K}{\left(a,z\right)}:=\int_{0}^{z}\mathrm{d}x\,\frac{\arctan{\left(ax\right)}}{1+x^{2}}.$$

In the special case $$a=1$$ the integral is elementary, and we have

$$\mathcal{K}{\left(1,z\right)}=\int_{0}^{z}\mathrm{d}x\,\frac{\arctan{\left(x\right)}}{1+x^{2}}=\frac12\arctan^{2}{\left(z\right)};~~~\small{z\in\mathbb{R}_{>0}}.$$

Suppose $$\left(a,z\right)\in\mathbb{R}_{>0}^{2}$$, and assume $$a\neq1$$. Then, $$-1<\frac{1-a}{1+a}<1\land\frac{1-a}{1+a}\neq0$$.

Set $$\omega:=2\arctan{\left(z\right)}$$, and note that $$0<\omega<\pi\land z=\tan{\left(\frac{\omega}{2}\right)}$$.

\begin{align} \mathcal{K}{\left(a,z\right)} &=\int_{0}^{z}\mathrm{d}x\,\frac{\arctan{\left(ax\right)}}{1+x^{2}}\\ &=\int_{0}^{z}\mathrm{d}x\,\frac{1}{1+x^{2}}\int_{0}^{a}\mathrm{d}t\,\frac{x}{1+x^{2}t^{2}}\\ &=\int_{0}^{z}\mathrm{d}x\int_{0}^{a}\mathrm{d}t\,\frac{x}{\left(1+x^{2}\right)\left(1+t^{2}x^{2}\right)}\\ &=\int_{0}^{a}\mathrm{d}t\int_{0}^{z}\mathrm{d}x\,\frac{x}{\left(1+x^{2}\right)\left(1+t^{2}x^{2}\right)}\\ &=\frac12\int_{0}^{a}\mathrm{d}t\int_{0}^{z^{2}}\mathrm{d}y\,\frac{1}{\left(1+y\right)\left(1+t^{2}y\right)};~~~\small{\left[x^{2}=y\right]}\\ &=\frac12\int_{0}^{a}\mathrm{d}t\int_{0}^{z^{2}}\mathrm{d}y\,\frac{d}{dy}\left[\frac{\ln{\left(1+y\right)}-\ln{\left(1+t^{2}y\right)}}{1-t^{2}}\right]\\ &=\frac12\int_{0}^{a}\mathrm{d}t\,\frac{\ln{\left(1+z^{2}\right)}-\ln{\left(1+z^{2}t^{2}\right)}}{1-t^{2}}\\ &=-\frac12\int_{0}^{a}\mathrm{d}t\,\frac{\ln{\left(\frac{1+z^{2}t^{2}}{1+z^{2}}\right)}}{1-t^{2}}\\ &=-\frac12\int_{\frac{1-a}{1+a}}^{1}\mathrm{d}x\,\frac{1}{2x}\ln{\left(\frac{1+z^{2}\left(\frac{1-x}{1+x}\right)^{2}}{1+z^{2}}\right)};~~~\small{\left[t=\frac{1-x}{1+x}\right]}\\ &=-\frac12\int_{\frac{1-a}{1+a}}^{1}\mathrm{d}x\,\frac{1}{2x}\ln{\left(\frac{\left(1+x\right)^{2}+z^{2}\left(1-x\right)^{2}}{\left(1+z^{2}\right)\left(1+x\right)^{2}}\right)}\\ &=-\frac12\int_{\frac{1-a}{1+a}}^{1}\mathrm{d}x\,\frac{1}{2x}\ln{\left(\frac{1+2\left(\frac{1-z^{2}}{1+z^{2}}\right)x+x^{2}}{\left(1+x\right)^{2}}\right)}\\ &=-\frac12\int_{\frac{1-a}{1+a}}^{1}\mathrm{d}x\,\frac{1}{2x}\left[\ln{\left(1+2x\cos{\left(\omega\right)}+x^{2}\right)}-2\ln{\left(1+x\right)}\right]\\ &=-\frac12\int_{\frac{1-a}{1+a}}^{1}\mathrm{d}x\,\frac{\ln{\left(1+2x\cos{\left(\omega\right)}+x^{2}\right)}}{2x}+\frac12\int_{\frac{1-a}{1+a}}^{1}\mathrm{d}x\,\frac{\ln{\left(1+x\right)}}{x}\\ &=-\frac12\int_{0}^{1}\mathrm{d}x\,\frac{\ln{\left(1+2x\cos{\left(\omega\right)}+x^{2}\right)}}{2x}+\frac12\int_{0}^{\frac{1-a}{1+a}}\mathrm{d}x\,\frac{\ln{\left(1+2x\cos{\left(\omega\right)}+x^{2}\right)}}{2x}\\ &~~~~~-\frac12\operatorname{Li}_{2}{\left(-1\right)}+\frac12\operatorname{Li}_{2}{\left(-\frac{1-a}{1+a}\right)}\\ &=+\frac12\operatorname{Li}_{2}{\left(-1\right)}+\frac18\omega^{2}+\frac12\int_{0}^{\frac{1-a}{1+a}}\mathrm{d}x\,\frac{\ln{\left(1+2x\cos{\left(\omega\right)}+x^{2}\right)}}{2x}\\ &~~~~~-\frac12\operatorname{Li}_{2}{\left(-1\right)}+\frac12\operatorname{Li}_{2}{\left(-\frac{1-a}{1+a}\right)}\\ &=\frac12\arctan^{2}{\left(z\right)}+\frac12\operatorname{Li}_{2}{\left(-\frac{1-a}{1+a}\right)}\\ &~~~~~-\frac12\int_{0}^{\frac{1-a}{1+a}}\mathrm{d}x\,\frac{(-1)\ln{\left(1-2x\cos{\left(\pi-\omega\right)}+x^{2}\right)}}{2x}\\ &=\frac12\arctan^{2}{\left(z\right)}+\frac12\operatorname{Li}_{2}{\left(-\frac{1-a}{1+a}\right)}-\frac12\operatorname{Li}_{2}{\left(\frac{1-a}{1+a},\pi-\omega\right)}.\blacksquare\\ \end{align}