Let
$$S(x,y,z)=\frac14\sqrt{(x+y+z) (-x+y+z) (x-y+z) (x+y-z) }\tag1$$
(note that it's Heron's formula for the area of a triangle with sides of lengths $x,y,z$).
I'm trying to evaluate the following integral in a closed form:
$$\mathcal U=\int_0^1\int_0^x\int_{x - y}^{x + y}\sqrt{S(x,y,z) } dz dy dx.\tag2$$
I wasn't able to evaluate it symbolically (either manually or using Mathematica), but using numerical integration and heuristic methods, I found a plausible form:
$$\mathcal U\stackrel{\color{gray}?}=\frac\pi{12\sqrt{2 }}+\frac{\Gamma\left(\frac14\right)^4\sqrt{10 }}{1440 \pi^2}\left(\left(\small\frac12-\frac1{\sqrt{5 \phi }}\right)\textbf K(\alpha)^2-\textbf K(\alpha) \textbf E(\alpha)\right)\\ \quad\quad\;\approx0.117599420842157114228246644831065494814051852697…,\tag3$$
where $\textbf K(\cdot),\textbf E(\cdot)$ are the complete elliptic integrals of the first and the second kind, and $\alpha=\frac1{\phi \sqrt2}+\frac1{\sqrt{2 \phi}}$.
Can we prove that $(3)$ is indeed the true value of the integral?
Update: You might ask, what heuristic methods could possibly give us the conjectured value? Here is a brief explanation.
First, calculate a sequence of values of $\int\!\int\!\int S^{2n}dz dy dx,n\in\mathbb N;$ their evaluation is straightforward and they are all rational numbers. Then, use FindSequenceFunction to find a possible recurrence relation for that sequence. We need at least $30$ elements of the sequence to get a result:
$$\small{(5+2 n)^2 (7+4 n) (9+4 n) (11+4 n) (27+20 n)\cdot a(n+2)=\\4 (2+n) (7+4 n) \left(6750+16947 n+15764 n^2+6416 n^3+960
n^4\right)\cdot a(n+1)+\\128 (1+n)^2 (2+n) (3+2 n) (3+4 n) (47+20 n)\cdot a(n)}\tag4$$
Of course, there is no guarantee that this relation holds for all larger $n$, but we cross our fingers 🤞 and conjecture that it does. Next, use FunctionExpand to find a possible explicit expression for the general term determined by this recurrence relation; it's complicated and involves generalized hypergeometric functions; we also need to manually adjust a periodic factor to ensure the result remains real even for non-integer $n$. After some simplifications we arrive at this:
$$\small\tfrac{4^n (20 n+27) \Gamma (2 n+2)\sqrt\pi}{(n+1) (16 n+12) \Gamma \left(2 n+\frac{7}{2}\right)}\,{_5F_4}\left(1,n+\tfrac{3}{2},n+\tfrac{3}{2},n+\tfrac{3}{2},n+\tfrac{47}{20};n+\tfrac{27}{20},n+\tfrac{7}{4},n+2,n+\tfrac{9}{4};-4\right)\tag5$$
Then we cross our fingers again 🤞 and conjecture that the same formula holds for all non-integer $n$ as well. Next, substitute $n=1/4$ and use FunctionExpand again to expand the generalized hypergeometric functions in terms of elliptic integrals. The result is quite large and unwieldy ($19$ terms), but it appears to match the integral numerically, which is a sign we might be on the right track. We cross our fingers one more time 🤞 and conjecture that only a few terms of the expression are actually linearly independent over $\mathbb Z,$ and the result can be significantly simplified. Using FindIntegerNullVector and high enough numeric precision, we find a possible basis, and a linear combination of its elements that numerically matches the integral. After some radical denestings (using ResourceFunction["RadicalDenest"]) and other simplifications, this gives us the conjectured value $(3).$
Best Answer
Define the function $S:\{(x,y,z)\in\mathbb{R}^{3}\mid0\le y\le x\land x-y\le z\le x+y\}\rightarrow\mathbb{R}$ bye the expression
$$S{\left(x,y,z\right)}:=\frac14\sqrt{(x+y+z)(-x+y+z)(x-y+z)(x+y-z)},$$
and let $\mathcal{U}$ denote the value of the triple integral
$$\mathcal{U}=\int_{0}^{1}\mathrm{d}x\int_{0}^{x}\mathrm{d}y\int_{x-y}^{x+y}\mathrm{d}z\,\sqrt{S{\left(x,y,z\right)}}\approx0.117599.$$
Observe that $S$ has the following scaling property:
$$\begin{align} S{\left(x,xt,xu\right)} &=\frac14\sqrt{(x+xt+xu)(-x+xt+xu)(x-xt+xu)(x+xt-xu)}\\ &=\frac14\sqrt{x^{4}(1+t+u)(-1+t+u)(1-t+u)(1+t-u)}\\ &=x^{2}\cdot\frac14\sqrt{(1+t+u)(-1+t+u)(1-t+u)(1+t-u)}\\ &=x^{2}S{\left(1,t,u\right)}.\\ \end{align}$$
We can use this property to trivialize one of the integrations of $\mathcal{U}$ and reduce it to a double integral:
$$\begin{align} \mathcal{U} &=\int_{0}^{1}\mathrm{d}x\int_{0}^{x}\mathrm{d}y\int_{x-y}^{x+y}\mathrm{d}z\,\sqrt{S{\left(x,y,z\right)}}\\ &=\int_{0}^{1}\mathrm{d}x\int_{0}^{1}\mathrm{d}t\,x\int_{1-t}^{1+t}\mathrm{d}u\,x\sqrt{S{\left(x,xt,xu\right)}};~~~\small{\left[y=xt, z=xu\right]}\\ &=\int_{0}^{1}\mathrm{d}x\int_{0}^{1}\mathrm{d}t\int_{1-t}^{1+t}\mathrm{d}u\,x^{2}\sqrt{x^{2}S{\left(1,t,u\right)}}\\ &=\int_{0}^{1}\mathrm{d}x\int_{0}^{1}\mathrm{d}t\int_{1-t}^{1+t}\mathrm{d}u\,x^{3}\sqrt{S{\left(1,t,u\right)}}\\ &=\int_{0}^{1}\mathrm{d}t\int_{1-t}^{1+t}\mathrm{d}u\int_{0}^{1}\mathrm{d}x\,x^{3}\sqrt{S{\left(1,t,u\right)}}\\ &=\int_{0}^{1}\mathrm{d}t\int_{1-t}^{1+t}\mathrm{d}u\,\frac14\sqrt{S{\left(1,t,u\right)}}.\\ \end{align}$$
With a couple more simple substitutions, we can rewrite this double integral in a slightly less cumbersome form:
$$\begin{align} \mathcal{U} &=\frac14\int_{0}^{1}\mathrm{d}t\int_{1-t}^{1+t}\mathrm{d}u\,\sqrt{S{\left(1,t,u\right)}}\\ &=\frac14\int_{0}^{1}\mathrm{d}t\int_{1-t}^{1+t}\mathrm{d}u\,\sqrt{\frac14\sqrt{(1+t+u)(-1+t+u)(1-t+u)(1+t-u)}}\\ &=\frac18\int_{0}^{1}\mathrm{d}t\int_{1-t}^{1+t}\mathrm{d}u\,\sqrt[4]{(1+t+u)(-1+t+u)(1-t+u)(1+t-u)}\\ &=\frac18\int_{0}^{1}\mathrm{d}t\int_{1-t}^{1+t}\mathrm{d}u\,\sqrt[4]{\left[(1+t)^{2}-u^{2}\right]\left[u^{2}-(1-t)^{2}\right]}\\ &=\frac18\int_{0}^{1}\mathrm{d}t\int_{\frac{1-t}{1+t}}^{1}\mathrm{d}x\,(1+t)\sqrt[4]{\left[(1+t)^{2}-(1+t)^{2}x^{2}\right]\left[(1+t)^{2}x^{2}-(1-t)^{2}\right]};~~~\small{\left[u=(1+t)x\right]}\\ &=\frac18\int_{0}^{1}\mathrm{d}t\int_{\frac{1-t}{1+t}}^{1}\mathrm{d}x\,(1+t)^{2}\sqrt[4]{\left(1-x^{2}\right)\left[x^{2}-\left(\frac{1-t}{1+t}\right)^{2}\right]}\\ &=\frac18\int_{0}^{1}\mathrm{d}y\,\frac{2}{(1+y)^{2}}\int_{y}^{1}\mathrm{d}x\,\frac{4}{(1+y)^{2}}\sqrt[4]{\left(1-x^{2}\right)\left(x^{2}-y^{2}\right)};~~~\small{\left[t=\frac{1-y}{1+y}\right]}\\ &=\int_{0}^{1}\mathrm{d}y\int_{y}^{1}\mathrm{d}x\,\frac{\sqrt[4]{\left(1-x^{2}\right)\left(x^{2}-y^{2}\right)}}{\left(1+y\right)^{4}}\\ &=\int_{0}^{1}\mathrm{d}x\int_{0}^{x}\mathrm{d}y\,\frac{\sqrt[4]{\left(1-x^{2}\right)\left(x^{2}-y^{2}\right)}}{\left(1+y\right)^{4}}.\\ \end{align}$$
By the general binomial theorem, a binomial with negative integer exponent has the following power series expansion:
$$\left(1-z\right)^{-p}=\sum_{n=0}^{\infty}\binom{p+n-1}{n}z^{n};~~~\small{p\in\mathbb{N}\land|z|<1}.$$
Using the lemma above to expand the $(1+y)^{-4}$ factor of the integrand of $\mathcal{U}$ as a binomial series, we then use technique of switching the order of summation and integration to find
$$\begin{align} \mathcal{U} &=\int_{0}^{1}\mathrm{d}x\int_{0}^{x}\mathrm{d}y\,\left(1+y\right)^{-4}\sqrt[4]{\left(1-x^{2}\right)\left(x^{2}-y^{2}\right)}\\ &=\int_{0}^{1}\mathrm{d}x\int_{0}^{x}\mathrm{d}y\,\sum_{n=0}^{\infty}\binom{n+3}{n}(-y)^{n}\sqrt[4]{\left(1-x^{2}\right)\left(x^{2}-y^{2}\right)}\\ &=\sum_{n=0}^{\infty}(-1)^{n}\binom{n+3}{n}\int_{0}^{1}\mathrm{d}x\int_{0}^{x}\mathrm{d}y\,y^{n}\sqrt[4]{\left(1-x^{2}\right)\left(x^{2}-y^{2}\right)}\\ &=\sum_{n=0}^{\infty}(-1)^{n}\binom{n+3}{n}\int_{0}^{1}\mathrm{d}x\int_{0}^{1}\mathrm{d}v\,x^{n+1}v^{n}\sqrt[4]{\left(1-x^{2}\right)x^{2}\left(1-v^{2}\right)};~~~\small{\left[y=xv\right]}\\ &=\sum_{n=0}^{\infty}(-1)^{n}\binom{n+3}{n}\int_{0}^{1}\mathrm{d}t\int_{0}^{1}\mathrm{d}u\,2^{-2}t^{n/2}u^{(n-1)/2}\sqrt[4]{t\left(1-t\right)\left(1-u\right)};~~~\small{\left[x^{2}=t,\,v^{2}=u\right]}\\ &=\frac14\sum_{n=0}^{\infty}(-1)^{n}\binom{n+3}{n}\int_{0}^{1}\mathrm{d}t\int_{0}^{1}\mathrm{d}u\,t^{n/2+1/4}\left(1-t\right)^{1/4}u^{(n-1)/2}\left(1-u\right)^{1/4}\\ &=\frac14\sum_{n=0}^{\infty}(-1)^{n}\binom{n+3}{n}\int_{0}^{1}\mathrm{d}t\,t^{n/2+1/4}\left(1-t\right)^{1/4}\int_{0}^{1}\mathrm{d}u\,u^{(n-1)/2}\left(1-u\right)^{1/4}\\ &=\frac14\sum_{n=0}^{\infty}(-1)^{n}\binom{n+3}{n}\operatorname{B}{\left(\frac{n}{2}+\frac54,\frac54\right)}\operatorname{B}{\left(\frac{n+1}{2},\frac54\right)}\\ &=\sum_{n=0}^{\infty}\frac{(-1)^{n}}{4}\cdot\frac{(n+3)!}{n!\,3!}\cdot\frac{\Gamma{\left(\frac{n}{2}+\frac54\right)}\,\Gamma{\left(\frac54\right)}}{\Gamma{\left(\frac{n}{2}+\frac52\right)}}\cdot\frac{\Gamma{\left(\frac{n+1}{2}\right)}\,\Gamma{\left(\frac54\right)}}{\Gamma{\left(\frac{n}{2}+\frac74\right)}}\\ &=\frac16\left[\Gamma{\left(\frac54\right)}\right]^{2}\sum_{n=0}^{\infty}(n+2)\,(-1)^{n}\,\frac{\Gamma{\left(\frac{n}{2}+\frac54\right)}}{\Gamma{\left(\frac{n}{2}+\frac74\right)}}\\ &=\frac{\left[\Gamma{\left(\frac54\right)}\right]^{2}}{8\,\Gamma{\left(\frac52\right)}}\sum_{n=0}^{\infty}(n+2)\,(-1)^{n}\,\frac{\Gamma{\left(\frac{n}{2}+\frac54\right)}\,\Gamma{\left(\frac12\right)}}{\Gamma{\left(\frac{n}{2}+\frac74\right)}}\\ &=\frac18\operatorname{B}{\left(\frac54,\frac54\right)}\sum_{n=0}^{\infty}(n+2)\,(-1)^{n}\operatorname{B}{\left(\frac{n}{2}+\frac54,\frac12\right)},\\ \end{align}$$
where the beta function and the gamma function are defined in the usual way by their respective integral representations,
$$\operatorname{B}{\left(x,y\right)}=\int_{0}^{1}\mathrm{d}t\,t^{x-1}\left(1-t\right)^{y-1};~~~\small{x>0\land y>0},$$
$$\Gamma{\left(z\right)}:=\int_{0}^{\infty}\mathrm{d}t\,t^{z-1}\exp{\left(-t\right)};~~~\small{z>0},$$
and are related through the identity
$$\operatorname{B}{\left(x,y\right)}=\frac{\Gamma{\left(x\right)}\,\Gamma{\left(y\right)}}{\Gamma{\left(x+y\right)}};~~~\small{x>0\land y>0}.$$
Define $f$ on $(-1,1)$ by the power series
$$f{\left(z\right)}=\sum_{n=0}^{\infty}z^{n}\operatorname{B}{\left(\frac{n}{2}+\frac54,\frac12\right)};~~~\small{|z|<1}.$$
Differentiating term-by-term, we have
$$f^{\prime}{\left(z\right)}=\sum_{n=0}^{\infty}nz^{n-1}\operatorname{B}{\left(\frac{n}{2}+\frac54,\frac12\right)},$$
and it follows that
$$2f{\left(z\right)}+zf^{\prime}{\left(z\right)}=\sum_{n=0}^{\infty}\left(n+2\right)z^{n}\operatorname{B}{\left(\frac{n}{2}+\frac54,\frac12\right)}.$$
We can convert the series representation for $f$ into an integral by expressing the beta function through its integral representation and summing under the integral sign:
$$\begin{align} f{\left(z\right)} &=\sum_{n=0}^{\infty}z^{n}\operatorname{B}{\left(\frac{n}{2}+\frac54,\frac12\right)}\\ &=\sum_{n=0}^{\infty}z^{n}\int_{0}^{1}\mathrm{d}t\,t^{n/2+1/4}\left(1-t\right)^{-1/2}\\ &=\sum_{n=0}^{\infty}z^{n}\int_{0}^{1}\mathrm{d}t\,\frac{t^{n/2+1/4}}{\sqrt{1-t}}\\ &=\sum_{n=0}^{\infty}z^{n}\int_{0}^{1}\mathrm{d}u\,\frac{2u^{n+1}\sqrt{u}}{\sqrt{1-u^{2}}};~~~\small{\left[t=u^{2}\right]}\\ &=\int_{0}^{1}\mathrm{d}u\,\sum_{n=0}^{\infty}\frac{2z^{n}u^{n+1}\sqrt{u}}{\sqrt{1-u^{2}}}\\ &=\int_{0}^{1}\mathrm{d}u\,\frac{2u\sqrt{u}}{\left(1-zu\right)\sqrt{1-u^{2}}}.\\ \end{align}$$
Then,
$$\begin{align} f^{\prime}{\left(z\right)} &=\frac{d}{dz}\int_{0}^{1}\mathrm{d}u\,\frac{2u\sqrt{u}}{\left(1-zu\right)\sqrt{1-u^{2}}}\\ &=\int_{0}^{1}\mathrm{d}u\,\frac{\partial}{\partial z}\frac{2u\sqrt{u}}{\left(1-zu\right)\sqrt{1-u^{2}}}\\ &=\int_{0}^{1}\mathrm{d}u\,\frac{2u^{2}\sqrt{u}}{\left(1-zu\right)^{2}\sqrt{1-u^{2}}},\\ \end{align}$$
and so
$$\begin{align} \sum_{n=0}^{\infty}\left(n+2\right)z^{n}\operatorname{B}{\left(\frac{n}{2}+\frac54,\frac12\right)} &=2f{\left(z\right)}+zf^{\prime}{\left(z\right)}\\ &=2\int_{0}^{1}\mathrm{d}u\,\frac{2u\sqrt{u}}{\left(1-zu\right)\sqrt{1-u^{2}}}+z\int_{0}^{1}\mathrm{d}u\,\frac{2u^{2}\sqrt{u}}{\left(1-zu\right)^{2}\sqrt{1-u^{2}}}\\ &=\int_{0}^{1}\mathrm{d}u\,\frac{2\left(2-zu\right)u\sqrt{u}}{\left(1-zu\right)^{2}\sqrt{1-u^{2}}}.\\ \end{align}$$
In the limit as $z\to-1^{+}$, this becomes
$$\sum_{n=0}^{\infty}\left(n+2\right)(-1)^{n}\operatorname{B}{\left(\frac{n}{2}+\frac54,\frac12\right)}=2\int_{0}^{1}\mathrm{d}u\,\frac{\left(2+u\right)u\sqrt{u}}{\left(1+u\right)^{2}\sqrt{1-u^{2}}}.$$
For any nonnegative integer $n$, observe how the following class of integrals transform:
$$\begin{align} \int_{0}^{1}\mathrm{d}x\,\frac{1}{\left(1+x\right)^{n}\sqrt{x\left(1-x^{2}\right)}} &=\int_{0}^{1}\mathrm{d}x\,\frac{1}{\left(1+x\right)^{n}\sqrt{x\left(1-x\right)\left(1+x\right)}}\\ &=\int_{0}^{1}\mathrm{d}x\,\frac{1}{\left(1+x\right)^{n+1}\sqrt{x\left(\frac{1-x}{1+x}\right)}}\\ &=\int_{0}^{1}\mathrm{d}y\,\frac{2}{\left(1+y\right)^{2}}\cdot\frac{1}{\left(\frac{2}{1+y}\right)^{n+1}\sqrt{y\left(\frac{1-y}{1+y}\right)}};~~~\small{\left[x=\frac{1-y}{1+y}\right]}\\ &=\int_{0}^{1}\mathrm{d}y\,\frac{\left(1+y\right)^{n}}{2^{n}\sqrt{y\left(1-y^{2}\right)}}.\\ \end{align}$$
Also for any nonnegative integer $n$, we find
$$\begin{align} \int_{0}^{1}\mathrm{d}x\,\frac{x^{n}}{\sqrt{x\left(1-x^{2}\right)}} &=\int_{0}^{1}\mathrm{d}x\,\frac{x^{n-1/2}}{\sqrt{1-x^{2}}}\\ &=\int_{0}^{1}\mathrm{d}t\,\frac{t^{n/2-1/4}}{2\sqrt{t}\sqrt{1-t}};~~~\small{\left[x=\sqrt{t}\right]}\\ &=\frac12\int_{0}^{1}\mathrm{d}t\,t^{n/2-3/4}\left(1-t\right)^{-1/2}\\ &=\frac12\operatorname{B}{\left(\frac{n}{2}+\frac14,\frac12\right)}.\\ \end{align}$$
We then obtain the following result:
$$\begin{align} \int_{0}^{1}\mathrm{d}x\,\frac{\left(2+x\right)x\sqrt{x}}{\left(1+x\right)^{2}\sqrt{1-x^{2}}} &=\int_{0}^{1}\mathrm{d}x\,\frac{\left(2+x\right)x^{2}}{\left(1+x\right)^{2}\sqrt{x\left(1-x^{2}\right)}}\\ &=\int_{0}^{1}\mathrm{d}x\,\left[x-\frac{1}{1+x}+\frac{1}{\left(1+x\right)^{2}}\right]\frac{1}{\sqrt{x\left(1-x^{2}\right)}}\\ &=\int_{0}^{1}\mathrm{d}x\,\frac{x}{\sqrt{x\left(1-x^{2}\right)}}-\int_{0}^{1}\mathrm{d}x\,\frac{1}{\left(1+x\right)\sqrt{x\left(1-x^{2}\right)}}\\ &~~~~~+\int_{0}^{1}\mathrm{d}x\,\frac{1}{\left(1+x\right)^{2}\sqrt{x\left(1-x^{2}\right)}}\\ &=\int_{0}^{1}\mathrm{d}x\,\frac{x}{\sqrt{x\left(1-x^{2}\right)}}-\int_{0}^{1}\mathrm{d}y\,\frac{\left(1+y\right)}{2\sqrt{y\left(1-y^{2}\right)}}\\ &~~~~~+\int_{0}^{1}\mathrm{d}y\,\frac{\left(1+y\right)^{2}}{2^{2}\sqrt{y\left(1-y^{2}\right)}};~~~\small{\left[x=\frac{1-y}{1+y}\right]}\\ &=\int_{0}^{1}\mathrm{d}x\,\frac{x}{\sqrt{x\left(1-x^{2}\right)}}-\frac14\int_{0}^{1}\mathrm{d}y\,\frac{1}{\sqrt{y\left(1-y^{2}\right)}}\\ &~~~~~+\frac14\int_{0}^{1}\mathrm{d}y\,\frac{y^{2}}{\sqrt{y\left(1-y^{2}\right)}}\\ &=\frac12\operatorname{B}{\left(\frac34,\frac12\right)}-\frac18\operatorname{B}{\left(\frac14,\frac12\right)}+\frac18\operatorname{B}{\left(\frac54,\frac12\right)}.\\ \end{align}$$
Returning finally to the evaluation of $\mathcal{U}$, we have
$$\begin{align} \mathcal{U} &=\frac18\operatorname{B}{\left(\frac54,\frac54\right)}\sum_{n=0}^{\infty}\left(n+2\right)(-1)^{n}\operatorname{B}{\left(\frac{n}{2}+\frac54,\frac12\right)}\\ &=\frac14\operatorname{B}{\left(\frac54,\frac54\right)}\int_{0}^{1}\mathrm{d}x\,\frac{\left(2+x\right)x\sqrt{x}}{\left(1+x\right)^{2}\sqrt{1-x^{2}}}\\ &=\frac14\operatorname{B}{\left(\frac54,\frac54\right)}\left[\frac12\operatorname{B}{\left(\frac34,\frac12\right)}-\frac18\operatorname{B}{\left(\frac14,\frac12\right)}+\frac18\operatorname{B}{\left(\frac54,\frac12\right)}\right]\\ &=\frac{\pi}{12\sqrt{2}}-\frac{\left[\Gamma{\left(\frac14\right)}\right]^{4}}{576\sqrt{2}\,\pi}.\blacksquare\\ \end{align}$$
Now, it would be interesting if your conjectured value is also correct, in which case my result provides an exact value for that combination of elliptic integrals in terms of gamma functions. Perhaps that can be explored later. Cheers!