Close form for a sequence of integrals involving the Gamma function I'm trying to find a close form for these integrals (with $n\in\mathbb{N}$) :
$$\int_{0}^{1}\ln(\Gamma(x+1))\cdot x^n \, \text{d}x$$
Wolfram is able to give a close form for every specific value of n, and the result leads to believe there is indeed a pattern ; however it fails to find the general formulae for $n\in\mathbb{N}$.
An integration by part turns it into similar-looking integrals involving the digamma function.
Do you think there is any way to find a close form, or is it a lost battle ?
 A: Feynman's trick does the job. We have
$$\begin{eqnarray*} \int_{0}^{1}\psi(s+x)\,dx &=& \log(s),\\
\int_{0}^{1}x\,\psi(s+x)\,dx &=& \log\Gamma(s+1)+\psi^{(-2)}(s)-\psi^{(-2)}(s+1),\\\int_{0}^{1}x^2\,\psi(s+x)\,dx &=&\log\Gamma(s+1)-2\psi^{(-3)}(s)+2\psi^{(-3)}(s+1)-2\psi^{(-2)}(s+1),\\
\int_{0}^{1}x^3\,\psi(s+x)\,dx &=&\log\Gamma(s+1)+6\psi^{(-4)}(s)-6\psi^{(-4)}(s+1)+6\psi^{(-3)}(s+1)-3\psi^{(-2)}(s+1)
\end{eqnarray*}$$
etcetera, where $\int_{0}^{1}\log\Gamma(s+1)\,ds =-1+\frac{1}{2}\log(2\pi)$ follows from the reflection formula and
$$ \int_{0}^{1}\psi^{(-2)}(s)\,ds = \log A+\tfrac{1}{4}\log(2\pi), $$
$$ \int_{0}^{1}\psi^{(-2)}(s+1)\,ds = \log A-\tfrac{3}{4}+\tfrac{3}{4}\log(2\pi), $$
$$ \int_{0}^{1}\psi^{(-3)}(s)\,ds = \tfrac{1}{2}\log A+\tfrac{1}{48}\log(2\pi)+\tfrac{1}{12}\cdot\frac{\zeta(3)}{\zeta(2)}, $$
$$ \int_{0}^{1}\psi^{(-3)}(s+1)\,ds = \tfrac{3}{2}\log A-\tfrac{11}{36}+\tfrac{7}{12}\log(2\pi)+\tfrac{1}{48}\cdot\frac{\zeta(3)}{\zeta(2)} $$
etcetera, where $A$ is the Glaisher-Kinkelin constant. Summarizing the given integrals can be expressed in terms of $\zeta$ and $\zeta'$ evaluated at $\mathbb{N}\setminus\{1\}$, i.e. in terms of the integrals $\int_{0}^{+\infty}\frac{x^s}{e^x-1}\,dx$ and $\int_{0}^{+\infty}\frac{x^s\log(x)}{e^x-1}\,dx$. This is also a consequence of Binet's first $\log\Gamma$ formula and Fubini's theorem.
Yet another way is to express $x^n$ in terms of Bernoulli polynomials, then exploit the Malmsten-Kummer Fourier series of $\log\Gamma$. This is probably the most efficient way.
No battle is lost with enough weapons in the arsenal :D
A: $$S_n=\int_{0}^{1}x^n\ln(\Gamma(x+1))\text{d}x$$
As already mentioned the integral is a finite sum of polygamma functions, which is not very exciting.
Infinite series are of interest, but not a closed form :
$$\ln(\Gamma(x+1))=\ln(x\Gamma(x))=ln(x)+\ln(\Gamma(x))$$
$$\ln(\Gamma(x))=-\ln(x)-\gamma x+\sum_{k=0}^\infty \frac{(-1)^k \zeta(k+2)x^{k+2}}{(k+2)}$$
http://functions.wolfram.com/GammaBetaErf/LogGamma/06/01/02/01/
$$S_n=\int_{0}^{1} x^n\left(-\gamma x+\sum_{k=0}^\infty \frac{(-1)^k \zeta(k+2)x^{k+2}}{(k+2)}\right)\text{d}x$$
$$S_n=-\frac{\gamma}{n+1} +\sum_{k=0}^\infty \frac{(-1)^k \zeta(k+2)}{(k+2)(k+3+n)}$$
$\gamma$ is the Euler-Macsheroni constant. $\zeta(n)$ is the Riemann zeta function.
To answer to the question, I doubt that a close form made of presently available special functions exists. 
