Integral $\int_0^1 \int_0^1 \cdots \int_0^1\frac{x_{1}^2+x_{2}^2+\cdots+x_{n}^2}{x_{1}+x_{2}+\cdots+x_{n}}dx_{1}\, dx_{2}\cdots \, dx_{n}=?$ How to evaluate this integral:
$$\int_0^1 \int_0^1 \cdots \int_0^1\frac{x_{1}^2+x_{2}^2+\cdots+x_{n}^2}{x_{1}+x_{2}+\cdots+x_{n}}dx_{1}\, dx_{2}\cdots \, dx_{n}=?$$
I'm making use of the integral identity:
$$\int_{0}^{+\infty }e^{-t(x_{1}+x_{2}\cdots +x_{n})}dt=\frac{1}{x_{1}+x_{2}\cdots +x_{n}}$$ and then reversing the order of integration with respect to time and space variables.
But for $n=1$, then such that, $$\int_{0}^{\infty }dt\int_{0}^{1}x^{2}e^{-tx}dx=\int_{0}^{\infty }\frac{2 - e^{-t}(2 + 2t+t^2)}{t^3}dt=\int_{0}^{1}x\,dx=\frac{1}{2},$$ and $$\int_0^1 \int_0^1 \cdots \int_0^1\frac{x_{1}^2+x_{2}^2+\cdots+x_{n}^2}{x_{1}+x_{2}+\cdots+x_{n}}dx_{1}\, dx_{2}\cdots \, dx_{n}\\=n\int_{0}^{+\infty }\frac{2 - e^{-t}\left ( 2 + 2t+t^2 \right )}{t^3}\left ( \frac{1-e^{-t}}{t} \right )^{n-1}dt.$$
 A: A little long as a comment:
Thanks to Maple, some first few terms are as follows:
$n=1:$ $I=\frac{1}{2}$
$n=2:$ $I=\frac{2^2\ln(2)}{3}-\frac{1}{3}$
$n=3:$ $I=-2^3\ln(2)+\frac{3^3\ln(3)}{4}-\frac{5}{4}$
$n=4:$ $I={2^6\ln(2)}-\frac{189\ln(3)}{5}-\frac{11}{5}$
$n=5:$ $I=-\frac{5\times 584\ln(2)}{9}+3^4\ln(3)+\frac{5^4\ln(5)}{36}-\frac{19}{6}$
$n=6:$ $I=\frac{6\times 11696\ln(2)}{63}+\frac{6\times 405\ln(3)}{14}-\frac{6\times 6250\ln(5)}{63}-\frac{29}{7}$
and so on. It seems there is some patterns which may help to obtain closed formula for the integral.
A: Ah, we may simply integrate by parts!
Denote $h(t)=(2 - e^{-t}( 2 + 2t+t^2))(1-e^{-t})^{n-1}$. Integrating by parts $(n+1)$ times we get $$\int_0^\infty h(t)t^{-n-2}dt=\frac1{(n+1)!}\int_0^\infty h^{(n+1)}(t)t^{-1}dt,$$
all off-integral terms vanish since $h(0)=h'(0)=\dots=h^{(n+1)}(0)=0$ and at infinity we have $h(t)=O(1)$, $h^{(i)}(t)=o(1)$ for $i>0$.
We have $h^{(n+1)}(t)=\sum_{k=1}^n a_ke^{-kt}+u(t)$, where $u(t)=t\times \text{polynomial}(e^{-t},t)$. Note that $\sum a_k=h^{(n+1)}(0)=0$, thus $$\int_0^\infty \sum a_ke^{-kt}t^{-1}dt=\int_0^\infty \sum a_k(e^{-kt}-e^{-t})t^{-1}dt=-\sum a_k \log k$$
by Frullani integrals. It is easy to see that $$a_k=(-1)^{k+n+1}\binom{n-1}{k-1}k^{n-1}(n^2+n-2k).$$
It remains to evaluate $\int_0^\infty u(t) t^{-1}dt$. We have $$u(t)=-(2t+t^2)(g(t))^{(n+1)}-2(n+1)t(g(t))^{(n)},\,g(t)=e^{-t}(1-e^{-t})^{n-1}.$$
Therefore $$\int_0^\infty u(t) t^{-1}dt=-2\int g^{(n+1)}(t)dt-2(n+1)\int g^{(n)}(t)dt-\int_0^{\infty}tg^{(n+1)}(t)dt.$$
We get $\int (g^{(n)}+tg^{(n+1)})=tg^{(n)}$, and the definite integral against $(0,\infty)$ equals 0. It remains $\int_0^\infty -2g^{(n+1)}-(2n+1)g^{(n)}=2g^{(n)}(0)+(2n+1)g^{(n-1)}(0)$. We have $$g(t)=(1-t+\dots)t^{n-1}(1-t/2+\dots)^{n-1}=t^{n-1}-\frac{n+1}2t^n+\dots,$$
$g^{(n-1)}(0)=(n-1)!$, $g^{n}(0)=-\frac{(n+1)!}2$, $2g^{(n)}(0)+(2n+1)g^{(n-1)}(0)=-(n^2-n-1)(n-1)!$,
confirming the guess of Sylvain JULIEN. 
PREVIOUS NON COMPLETE VERSION
Still not a complete answer, but a method to be completed or improved.
Denote $f_n(t)=\frac{2 - e^{-t}\left ( 2 + 2t+t^2 \right )}{t^3}\left ( \frac{1-e^{-t}}{t} \right )^{n-1}$. Assume that we have found the coefficients $c_1,c_2,\dots,c_n$ and polynomials $g_1,\dots,g_n$ such that $$f_n(t)=2t^{-n-2}+\sum_{k=1}^n(g_k(t^{-1})e^{-kt})'+\sum_{k=1}^n c_k \frac{e^{-kt}}t.\,\,(*)$$
Then $\sum c_k=0$ (else $f_n$ would have a non-zero residue at 0, which is absurd). We have $\int_0^\infty \sum_{k=1}^n c_k \frac{e^{-kt}}t dt=-\sum c_k\log k$ by Frullani integrals.
The integral of the other part of our sum $(*)$ is minus the limit at zero of the function $-\frac2{n+1}t^{-n-1}+\sum g_k(t^{-1})e^{-kt}$. The limit must exist, since the initial integral converges. 
Now I cheat a bit. Note that if we write $f_n(t)=\sum_{k=0}^n q_k(1/t) e^{-kt}$ for polynomials $q_k$, then $c_k$ equals to the residue of the $q_k(1/t) e^{-kt}$, which is pretty computable. If I am not mistaken, $$q_k(t^{-1})=(-1)^kt^{-n-2}\left(2\binom{n-1}k+\binom{n-1}{k-1}(2+2t+t^2)\right)=\\=(-1)^kt^{-n-2}\left(2\binom{n}k+\binom{n-1}{k-1}(2t+t^2)\right).$$ Thus $$c_k=2(-1)^{k+n+1}\binom{n}k\frac{k^{n+1}}{(n+1)!}+2(-1)^{k+n}\binom{n-1}{k-1}\frac{k^{n}}{n!}+(-1)^{k+n+1}\binom{n-1}{k-1}\frac{k^{n-1}}{(n-1)!}=\\
=(-1)^{k+n+1}\frac{k^{n-1}(n^2+n-2k)}{n(n+1)(n-k)!(k-1)!}.$$
This matches a coefficient of $\log 3$ for $n=4$ from Shahrooz Janbaz's answer, you may check others for be sure. 
It remains to prove Sylvain JULIEN's guess for the rational part. 
