New Identities
For any positive integers $k\geq1$ and $j\geq2$, let $x_j=\frac{j+\sqrt{j^2-4}}{2}$. Let us define $(\text{A}_k)_{k\geq1}$ by the following constants "which are variants of Stieltjes constants" :
$$\text{A}_k = \lim_{n\to\infty}\bigg(\sum_{j=2}^{n}\frac{\ln^k(x_j)}{j}-\frac{\ln^{k+1}(x_n)}{k+1}\bigg)$$
Then following a straightforward calculus we can get the following identity :
$$\int_{0}^{1}\int_{0}^{1}\int_{0}^{1}\bigg\{x\,y\,z+\frac{1}{x\,y\,z}\bigg\}dx\,dy\,dz=\frac{13}{8}-\gamma-\text{A}_1-\frac{\text{A}_2}{2}$$
And more generally we have :
$$\int_{0}^{1}\int_{0}^{1}\dots\int_{0}^{1}\bigg\{\prod_{k=1}^{n}x_k+\frac{1}{\prod_{k=1}^{n}x_k}\bigg\}dx_1\,dx_2\dots\,dx_n=\frac{13}{8}-\gamma-\sum_{j=1}^{n-1}\frac{\text{A}_j}{j!}$$
Where $\gamma$ represents the Euler-Mascheroni constant and $\{\}$ denotes the fractional part.
Is it possible to evaluate $\text{A}_1$ and $\text{A}_2$ in closed-form ? In such a way we can plug them in the above integral.
The best I can come up with is a well-behaved integral for $A_1.$ There appears little hope of it being specified in terms of other well-known constants. The integral's complexity dissuades me for looking for a 'closed form' for $A_2.$ $$A_1 = \lim_{n \to \infty} \Big( \sum_{j=2}^n \frac{1}{j} \log{(\tfrac{1}{2}(j+\sqrt{j^2-4})) } - \frac{1}{2}\log^2{(\tfrac{1}{2}(n+\sqrt{n^2-4})) } \Big)$$ $$ = \lim_{n \to \infty} \Big( \sum_{j=2}^n \frac{1}{j} \log{(\tfrac{1}{2}(j+\sqrt{j^2-4})) } - \frac{1}{2}\log^2{n) } \Big)$$ $$ = \lim_{n \to \infty} \Big( \sum_{j=2}^n \frac{1}{j} \log{(\tfrac{1}{2}(j+\sqrt{j^2-4})) } - \frac{1}{j}\log{j) } \Big) + \gamma_1$$ where the definition of the ordinary Stieltjes constant has been used. Combine the logs and we now must deal with the infinite sum $$ U:= \sum_{j=2}^\infty \frac{1}{j} \log{(\tfrac{1}{2}(1+\sqrt{1-1/(j/2)^2})) } = \frac{1}{\pi} \int_0^\pi dx \,x\sin{x}\cos{x} \sum_{j=2}^\infty \frac{1}{j} \frac{1}{(j/2)^2 - \sin(x)^2}.$$
The last integral used a formula from Gradshteyn & Rhyzhik 3.812.12, though my version of the book has a typographical error and the correct version is $$\frac{1}{\pi} \int_0^\pi dx \,x\frac{\sin{x}\cos{x} }{a^2 - \sin^2{x}} =\log{(\tfrac{1}{2}(1+\sqrt{1-1/a^2})) }\quad,\quad a^2 > 1.$$ In fact, the integral is good for $a=1$ if the integral is interpreted in principal value sense, which is necessary next step after the definition of $U.$ The reason this integral was chosen was because the sum can be done in closed form, $$\sum_{j=2}^\infty \frac{1}{j} \frac{1}{(j/2)^2 - a^2} = \frac{2-2\gamma-\psi(2-2a)-\psi(2+2a)}{2a^2}.$$ Set $x \to x+\pi/2$ in the integrand and use trig id's to get the form, $$U=\frac{-1}{2\pi} \int_{-\pi/2}^{\pi/2} dx\tan{x}\, \big(x+\pi/2) \big(2-2\gamma -\psi(4\sin^2{(x/2)}) -\psi(4\cos^2{(x/2)}).$$ The integrand breaks into and even and odd part, and the odd part will equate to zero because of the symmetric integral limits, and it also rids us of the implicit interpretation as a principal value integral. Collecting, we finally have $$ A_1 = \gamma_1 - \frac{1}{\pi}\int_{0}^{\pi/2} dx\,x\,\tan{x}\, \big(2-2\gamma -\psi(4\sin^2{(x/2)}) -\psi(4\cos^2{(x/2)}). $$ Numerically, a brute force summation with n=50000 gave $A_1 \sim -0.50735$ whereas the integral and Stieltjes constant method easily give $A_1=-0.50746156086...$ Despite the digamma function, the integrand is very smooth.