Is there a closed form for the integral $\int_0^1 x^n \log^m (1-x) \, {\rm d}x$? Let $n \in \mathbb{N}$. We know that:
$$\int_0^1 x^n \log(1-x) \, {\rm d}x = - \frac{\mathcal{H}_{n+1}}{n+1}$$
Now, let $m , n \in \mathbb{N}$. What can we say about the integral
$$\int_0^1 x^n \log^m (1-x) \, {\rm d}x$$
For starters we know that $\displaystyle \log^m (1-x)=m! \sum_{k=m}^{\infty} (-1)^k \frac{s(k, m)}{k!} x^k$ where $s(k, m)$ are the Stirling numbers of first kind. 
Thus
\begin{align*}
\int_{0}^{1} x^n \log^m (1-x) \, {\rm d}x &=m! \int_{0}^{1}x^n \sum_{k=m}^{\infty} (-1)^k \frac{s(k, m)}{k!} x^k  \\ 
 &= m! \sum_{k=m}^{\infty} (-1)^k \frac{s(k, m)}{m!} \int_{0}^{1}x^{n+m} \, {\rm d}x\\ 
 &= m! \sum_{k=m}^{\infty} (-1)^k \frac{s(k, m)}{m!} \frac{1}{m+n+1}
\end{align*}
Can we simplify? I know that Striling numbers are related to the Harmonic number but I don't remember all identities. 
 A: Change variables $x=1-t$:
$$
I=\int_0^1 x^n \log^m (1-x) \, {\rm d}x=\int_0^1dt (1-t)^n\log^m t=\sum_{k=0}^n {n\choose k}(-1)^{n-k}\int_0^1 dt \ t^{n-k}\log^m t\ .
$$
Now change variable $t=\exp(z)$ and get
$$
I=\sum_{k=0}^n {n\choose k}(-1)^{n-k}\int_{-\infty}^0 dz\ e^{(n-k+1)z}z^m= 
\boxed{\Gamma (m+1)\sum_{k=0}^n {n\choose k}\frac{(-1)^{n-k+m} }{ (n+1-k)^{m+1}}}\ ,
$$
which is a finite sum.
A: Another closed form follows by differentiating the beta function multiple times and applying the Faà di Bruno's formula.

Claim. For positive integers $m$ and $n$,
  $$ \mathcal{J}_{n,m} := \int_0^1 x^{n-1}\log^m (1-x) \, \mathrm{d}x = (-1)^m \frac{m!}{n} \sum_{\alpha\in I_m} \prod_{k=1}^m \frac{1}{\alpha_k!} \bigg(\frac{H_n^{(k)}}{k}\bigg)^{\alpha_k} \tag{1} $$
  where $\alpha$ runs over the set of indices
  $$I_m = \{(\alpha_1,\cdots,\alpha_m)\in\Bbb{N}_0^m : 1\cdot\alpha_1+\cdots+m\cdot\alpha_m=m\}.$$

This formula gives an almost explicit formula for $\mathcal{J}_{n,m}$ in terms of polynomial of $H_n^{(1)}, \cdots, H_n^{(n)}$ at the expense of introducing certain combinatorial object, namely $I_m$.
Proof. Notice that
$$ \int_0^1 x^{n-1}(1-x)^s \, \mathrm{d}x = \frac{(n-1)!}{(s+1)\cdots(s+n)} = (n-1)!\exp\left(-\sum_{j=1}^n \log(s+j) \right). $$
Letting $f(s) = -\sum_{j=1}^n \log(s+j) $ and applying the Faà di Bruno's formula, we have
$$ \mathcal{J}_{n,m} = (n-1)!e^{f(0)} \sum_{\alpha \in I_m} m! \prod_{k=1}^{m} \frac{1}{\alpha_k !} \bigg( \frac{f^{(k)}(0)}{k!} \bigg)^{\alpha_k}. \tag{2}$$
Plugging $f(0) = -\log n!$ and
$$ f^{(k)}(0) = \sum_{j=1}^n (-1)^k (k-1)! (s+j)^{-k} \bigg|_{s=0} = (-1)^k (k-1)! H_n^{(k)} $$
into $\text{(2)}$ and simplifying the resulting expression yields $\text{(1)}$.
A: We have $$(-1)^m n \int_0^1 x^{n-1} \ln^m(1-x)dx=\lim_{z\to 0} \frac{d^m}{dz^m} \,n \int_0^1 x^{n-1} (1-x)^{-z}dx\\ =\lim_{z\to 0} \frac{d^m}{dz^m} \, \frac{n!\, \Gamma(1-z)}{\Gamma(n-z+1)}=\lim_{z\to 0} \frac{d^m}{dz^m} \, \prod_{k=1}^n \frac{1}{1-z/k}. \tag{1}$$
But, in view of the generating function of the Complete homogeneous symmetric polynomials, we have $$\prod_{k=1}^n \frac{1}{1-z/k} = \sum_{k=0}^{\infty} h_k\,z^k\tag{2}$$
Where $$h_k \equiv h_k(1,1/2,\dots ,1/n)=\sum_{1 \leq a_1\leq a_2\leq\dots\leq a_k\leq n} \,\frac1{a_1\,a_2\cdots a_k}\tag{3}$$
Hence $$(-1)^m n \int_0^1 x^{n-1} \ln^m(1-x)dx= m!\,h_m. \tag{4}$$
The  Newton–Girard formulae connect the symmetric polynomials to their corresponding power sums, which in our case are the generalized harmonic numbers.
The first few cases are:
$$-\int_0^1 x^{n-1} \ln(1-x) dx=\frac{H_n}{n}\tag{5}$$
$$\int_0^1 x^{n-1} \ln^2(1-x) dx=\frac{H_n^2+H_n^{(2)}}{n}\tag{6}$$
$$-\int_0^1 x^{n-1} \ln^3(1-x) dx=\frac1{n}(H_n^3+3H_n\,H_n^{(2)}+2H_n^{(3)})\tag{7}$$
$$\int_0^1 x^{n-1} \ln^4(1-x) dx=\frac1{n}(H_n^4+6H_n^2\,H_n^{(2)}+3H_n^{(2)2}+8H_n\,H_n^{(3)}+6H_n^{(4)})\tag{8}$$
With some more algebraic manipulations, we can get the form
$$f(m)=\int_0^1 x^{n-1} \ln^m(1-x) dx=\frac{m!}{n\,(2m)!} \frac{d^{2m}}{dx^{2m}} \exp\left(\sum_{k=1}^{\infty} \frac{(-1)^k x^{2k}}{k} H_n^{(k)}\right)\Bigg{|}_{x=0}\tag{9}$$
Mathematica code:

f[m_] := m!/(n (2 m)!) D[
        Exp[Sum[(-1)^k x^(2 k)/k H[k], {k, 1, 2 m}]], {x, 2 m}] /. 
      x -> 0 /. H[int_] :> HarmonicNumber[n, int] // Simplify

This code becomes less efficient for large $m$'s.
