From experimentation I am pretty certain that $$ G_3(n) = \int_0^\infty\int_0^\infty \int_0^\infty e^{-x_1-x_2-x_3}\max(|x_1-x_2|,|x_2-x_3|,|x_3-x_1|)^n\; dx_1 \,dx_2 \,dx_3= \Gamma(n+1)(2-2^{-n}) $$ The lower-dimensional version seems to be $$ G_2(n) = \int_0^\infty\int_0^\infty e^{-x_1-x_2}\max(|x_1-x_2|)^n\; dx_1 \,dx_2 = \int_0^\infty x^n e^{-x} \; dx= \Gamma(n+1) $$ so it is not unreasonable. Is there a way of proving/showing: $$ G_3(n) = \Gamma(n+1)(2-2^{-n})\text{?} $$ I'm interested in generalisations to higher dimensions.
Edit:
Thanks to Jack's answer I have also found $$ G_4(n) = \Gamma(1+n)(3-3 \cdot 2^{-n} + 3^{-n}) $$ assuming the max of all $\binom{4}{2}$ differences of variables. It appears that for general $k$ we will have $$ G_k(n) = (k-1) \int_0^\infty x^{n} e^{-(k-1)x}(e^{x}-1)^{k-2} \; dx $$
By symmetry $$\begin{eqnarray*} G_3(n)&=&6\iiint_{0\leq x_1\leq x_2 \leq x_3} e^{-(x_1+x_2+x_3)}(x_3-x_1)^n\,dx_1\,dx_2\,dx_3\\&=&6\iint_{0\leq x_1\leq x_3}e^{-(x_1+x_3)}\left(e^{-x_1}-e^{-x_3}\right)(x_3-x_1)^n\,dx_1\,dx_3\\&\stackrel{x_3\mapsto x_1+t}{=}&6\iint_{(0,+\infty)^2}e^{-(3x_1+t)}\left(1-e^{-t}\right)t^n\,dx_1\,dt\\&=&2\int_{0}^{+\infty}\left(e^{-t}-e^{-2t}\right)t^n\,dt\\&=&2\left(n!-\frac{1}{2^{n+1}}n!\right)=n!\left(2-\frac{1}{2^n}\right).\qquad\square\end{eqnarray*} $$ The generalization to $G_k(n)$ is straightforward.