How can we show that $\int_{0}^{1}(\ln{x})^n\ln(-\ln{x})\mathrm dx=(-1)^n(a-b\gamma)\cdot{n!\over b}?$ Consider the integral $(1)$

$$\int_{0}^{1}(\ln{x})^n\ln(-\ln{x})\mathrm dx=I\tag1$$
  $H_n={a\over b}$, it is the nth-harmonic number

$H_n:=1,{3\over 2},{11\over 6},{25\over 12},...$ for $n:=1,2,3,4,... $ respectively.
We have
$$I=(-1)^n(a-b\gamma)\cdot{n!\over b}$$
[We think the closed form is sort of interesting one, because we haven't seem Euler's constant is mixing with the numerator and denumerator of harmonic number]
An attempt:
$u=-\ln{x}$ then $-xdu=dx$ 
simplified to
$$\int_{0}^{\infty}\color{red}{u^n\ln{u}}e^{-u}\mathrm du\tag2$$
We can apply the Laplace transform to $(2)$, but I can't find the Laplace transform of $F(u)=u^n\ln{u}$
or else we can use $$e^{-x}=1-x+{x^2\over 2!}-{x^3\over 3!}+\cdots$$
$$\int_{0}^{\infty}\left(u^n-u^{n+1}+{u^{n+2}\over 2!}-{x^{n+3}\over 3!}+\cdots\right)\ln{u}\mathrm du\tag3$$
$(3)$ does not converge; let $n=1$ for example, the first term of the integral
$$\int_{0}^{\infty}u\ln{u}\mathrm du={1\over 2}u^2\ln{u}-{u^2\over 2}|_{0}^{\infty}$$
How else can we go about tackling $(1)?$

Hint. One may recall the standard integral representation of the Euler gamma function,
$$
\Gamma(s+1)=\int_{0}^{\infty}u^{s} e^{-u}\mathrm du,\quad s>-1,
$$ then by differentiating under the integral sign, one gets
$$
\Gamma'(s+1)=\int_{0}^{\infty}u^{s}\cdot\ln u \cdot e^{-u}\mathrm du,\quad s>-1,
$$ or
$$
\frac{\Gamma'(s+1)}{\Gamma(s+1)}=\frac1{\Gamma(s+1)}\int_{0}^{\infty}u^{s}\cdot\ln u \cdot e^{-u}\mathrm du,\quad s>-1,
$$ then one may recall that the digamma function $\psi$ is such that 
$$
\psi(n+1):=\frac{\Gamma'(n+1)}{\Gamma(n+1)}=H_{n}-\gamma, \quad n=1,2,\cdots,
$$ where $\gamma$ is the Euler-Mascheroni constant and we conclude with the change of variable $u=-\ln x$ using $\Gamma(n+1)=n!$.
