Unfortunately, this is not new, though I would like to offer another derivation; \begin{align*} \int_{0}^{t}H_ydy &= \int_{0}^{t}\int_{0}^{1}\frac{1-x^y}{1-x}dxdy\\ &= \int_{0}^{1}\int_{0}^{t}\frac{1-x^y}{1-x}dydx\\ &= \int_{0}^{1}\frac{t}{1-x}+\frac{1-x^t}{(1-x)\ln(x)}dx\\ (1)&= \int_{0}^{1}\frac{t}{1-x}+\sum_{j=0}^{t-1}\frac{x^j}{\ln(x)}dx\\ &=\lim_{x\rightarrow 1^{-}}\left( -t\ln(1-x)+\sum_{j=0}^{t-1}\text{li}(x^{j+1})\right)\\ (2)&=\lim_{x\rightarrow 1^{-}}\sum_{j=0}^{t-1}\left(\text{li}(x^{j+1})-\ln(1-x)\right)\\ &= \gamma t + \sum_{j=0}^{t-1}\ln(j+1)\\ &= \gamma t + \ln(t!)\\ &= \gamma t +\ln\Gamma(t+1) \end{align*}

Now writing $H_y$ as $\sum_{n=1}^{\infty}\frac{y}{n(n+y)}$ gives \begin{align*} \int_{0}^{t}H_ydy &= \gamma t + \ln\Gamma(t+1)\\ &= \int_{0}^{t}\sum_{n=1}^{\infty}\frac{y}{n(n+y)}dy\\ &= \sum_{n=1}^{\infty}\frac{t}{n}-\ln \left(1+\frac{t}{n}\right) \end{align*} solving for $\Gamma(t+1)$ and using $\Gamma(t+1) = t\Gamma(t)$ gives $$\Gamma(t) = \frac{e^{-\gamma t}}{t}\prod_{n=1}^{\infty}\left(1+\frac{t}{n}\right)^{-1}e^{\frac{t}{n}}$$

EDIT First, I should say that I am just starting to study analysis - I just try stuff and see if it works, so I apologize for the lack of rigour.

I'm not too sure which steps I should fill in, so I will try to fill in any apparent holes.

assuming we only know that $\gamma$ shows up as $$\lim_{N \rightarrow \infty}(H_N - \ln(N)) = \gamma,$$ one can show that $\int_{0}^{1}H_ydy=\gamma$ using $H_y=\sum_{k=2}^{\infty}(-1)^k\zeta(k)y^{k-1}$ for $|y|<1$ (I have a crude derivation here).

Knowing that, we get \begin{align*} \gamma = \int_{0}^{1}H_ydy &= \int_{0}^{1}\int_{0}^{1}\frac{1-x^y}{1-x}dxdy\\ &= \int_{0}^{1}\int_{0}^{1} \frac{1}{1-x}-\frac{x^y}{1-x}dydx\\ &= \int_{0}^{1}\frac{1}{1-x}+\frac{1}{\ln(x)}dx\\ (*)&= \lim_{x\rightarrow 1^{-}}(\text{li}(x)-\ln(1-x)),\\ \end{align*} where $\text{li}(x) = \int_{0}^{x}\frac{1}{\ln(t)}dt$ is the logarithmic integral.

Now we need to evaluate $\int_{0}^{t}\frac{x^k}{\ln(x)}dx$ for (1); Make the substitution $\;x=e^u$ to get $$\int_{0}^{t}\frac{x^k}{\ln(x)}dx = \int_{-\infty}^{\ln(t)}\frac{e^{(k+1)u}}{u}du,$$ now let $v=(k+1)u$ and see that \begin{align*} \int_{0}^{t}\frac{x^k}{\ln(x)}dx &= \int_{-\infty}^{(k+1)\ln(t)}\frac{e^x}{x}dx\\ &= \int_{0}^{t^{k+1}}\frac{1}{\ln(x)}dx\\ &=\text{li}(t^{k+1}) \end{align*} after the substitution $x=\ln(u)$.

Next we need to evaluate $\lim_{x\rightarrow 1^{-}}(\text{li}(x^k)-\ln(1-x))$. so by using $(*)$ - \begin{align*} \lim_{x\rightarrow 1^{-}}(\text{li}(x^k)-\ln(1-x)) &= \lim_{x\rightarrow 1^{-}}(\text{li}(x) - \ln(1-\sqrt[k]{x})\\ &= \lim_{x\rightarrow 1^{-}}\left(\text{li}(x) - \ln \left(\frac{1-x}{1+\sqrt[k]{x} + \cdots +\sqrt[k]{x^{k-1}}}\right)\right)\\ &= \gamma + \ln(k) \end{align*}

I just have to say that I really like how you derived it.


Integral of the Logarithmic Derivative of Gamma

The logarithmic derivative of the Gamma function is the digamma function: $$ \frac{\Gamma'(x)}{\Gamma(x)}=\psi(x)=-\gamma+H_{x-1}\tag{1} $$ Therefore, $$ \begin{align} \log(\Gamma(x)) &=\int_1^x\left(-\gamma+H_{\phi-1}\right)\mathrm{d}\phi\\ &=\int_1^x\left(-\gamma+\int_0^1\frac{1-t^{\phi-1}}{1-t}\,\mathrm{d}t\right)\mathrm{d}\phi\\ &=\int_0^{x-1}\left(-\gamma+\int_0^1\frac{1-t^\phi}{1-t}\,\mathrm{d}t\right)\mathrm{d}\phi\tag{2} \end{align} $$


Verification of the Gamma Function

The Bohr-Mollerup Theorem says that the Gamma function is uniquely determined as the log-convex function so that $\Gamma(1)=1$ and $\Gamma(x+1)=x\,\Gamma(x)$.

We can verify these assuming only $H_x-H_{x-1}=\frac1x$ and $H'_x\ge0$.

$\boldsymbol{\Gamma(1)=1}$: $$ \begin{align} \log(\Gamma(1)) &=\int_1^1\left(-\gamma+H_{\phi-1}\right)\mathrm{d}\phi\\ &=0\tag{3} \end{align} $$ $\boldsymbol{\Gamma(x+1)=x\,\Gamma(x)}$: Since $\lim\limits_{n\to\infty}\left(H_n-\log(n)\right)=\gamma$ and $H_n-\frac1n\le\int_{n-1}^nH_\phi\,\mathrm{d}\phi\le H_n$, $$ \begin{align} \log(\Gamma(x+1))-\log(\Gamma(x)) &=\int_x^{x+1}\left(-\gamma+H_{\phi-1}\right)\mathrm{d}\phi\\ &=-\gamma+\lim_{n\to\infty}\left(\int_x^nH_{\phi-1}\,\mathrm{d}\phi-\int_{x+1}^nH_{\phi-1}\,\mathrm{d}\phi\right)\\ &=-\gamma+\lim_{n\to\infty}\left(\int_x^nH_{\phi-1}\,\mathrm{d}\phi-\int_x^{n-1}H_\phi\,\mathrm{d}\phi\right)\\ &=-\gamma+\lim_{n\to\infty}\left(-\int_x^n\frac1\phi\,\mathrm{d}\phi+\int_{n-1}^nH_\phi\,\mathrm{d}\phi\right)\\[6pt] &=-\gamma+\log(x)+\lim_{n\to\infty}\left(-\log(n)+H_n\right)\\[8pt] &=\log(x)\tag{4} \end{align} $$ $\boldsymbol{\Gamma}$ is log-convex: $$ \begin{align} \frac{\mathrm{d}^2}{\mathrm{d}x^2}\log(\Gamma(x)) &=H'_{x-1}\\ &\ge0\tag{5} \end{align} $$


Verifying the Necessary Properties of the Extension $$ \begin{align} H_x-H_{x-1} &=\int_0^1\frac{1-t^x}{1-t}\,\mathrm{d}t-\int_0^1\frac{1-t^{x-1}}{1-t}\,\mathrm{d}t\\ &=\int_0^1t^{x-1}\,\mathrm{d}t\\ &=\frac1x\tag{6} \end{align} $$ $$ \begin{align} H'_x &=\int_0^1\frac{-\log(t)t^x}{1-t}\,\mathrm{d}t\\ &\ge0\tag{7} \end{align} $$


Limitations of the Extension

One limitation of the extension $$ H_x=\int_0^1\frac{1-t^x}{1-t}\,\mathrm{d}t\tag{8} $$ is that it doesn't converge for $\mathrm{Re}(x)\le-1$.

An extension of the Harmonic Numbers that works for all $x\in\mathbb{C}$ is $$ H_x=\sum_{k=1}^\infty\left(\frac1k-\frac1{k+x}\right)\tag{9} $$