Taylor Series of Gamma Function I have always wondered whether the Taylor Series of Gamma Function exists or not. I tried to find it, but in vain. I googled for it, but couldn't find it. Has anyone ever found its Taylor Series?
 A: If you want the Taylor series, you basically need the $n^{th}$ derivative of $\Gamma(x)$. These express in terms of the polygamma function. Considering 
$$d_n=\frac{\left[\Gamma(x)\right]^{(n)}}{\Gamma(x)}$$ the first terms are
$$d_1=\psi ^{(0)}(x)$$
$$d_2=\psi ^{(0)}(x)^2+\psi ^{(1)}(x)$$
$$d_3=\psi ^{(0)}(x)^3+3 \psi ^{(1)}(x) \psi ^{(0)}(x)+\psi ^{(2)}(x)$$
$$d_4=\psi ^{(0)}(x)^4+6 \psi ^{(1)}(x) \psi ^{(0)}(x)^2+4 \psi ^{(2)}(x) \psi
   ^{(0)}(x)+3 \psi ^{(1)}(x)^2+\psi ^{(3)}(x)$$ which "simplify" (a little !) when you perform the expansion around $x=a$, $a$ being a positive integer. 
A: See the Euler Mascheroni Integrals $\text I_n$:
$$(-1)^n\text I_n=\int_0^\infty e^{-z}\ln^n(z)dz=\Gamma^{(n)}(1)$$
Therefore:
$$\Gamma(x)=\sum_{n=0}^\infty \frac{\Gamma^{(n)}(1)x^n}{n!}= \sum_{n=0}^\infty \frac{\int_0^\infty e^{-z}\ln^n(z)dz\,x^n}{n!}  =\sum_{n=0}^\infty \frac{\text I_n\ (-x)^n}{n!} $$
where $(n)$ is the $n$ derivative. From the linked article:
$$\Gamma^{(0)}(1)=\text I_0=1$$
$$\Gamma^{(1)}(1)=-\text I_1=-\text{eulergamma}=-\gamma$$
$$\Gamma^{(2)}(1)=\text I_2=\gamma^2+\frac{\pi^2}6$$
$$\Gamma^{(3)}(1)=-\text I_3=-\gamma^3-\frac{\gamma\pi^2}2-2\zeta(3)$$
$$…$$
where $\zeta$ is the Riemann Zeta function,$\zeta(3)$ is Apéry’s constant, and $n\in\Bbb N$ for this problem. Please correct me and give feedback!
Here is an evaluation of the “incomplete $\text I_n$” integral with numerical proof for $n=3$,$u=4$:
$$\int_0^u e^{-t} \ln(t)^n dt=\sum_{m=0}^\infty \frac{(-1)^m}{m!}\int_0^u \ln(t)^n t^m dt=-\ln(u)^{n+1}\sum_{m=0}^\infty \frac{(-1)^m}{m!}\text E_{-n}(-(m+1)\ln(u)$$
where appears the Exponential Integral E function. See @Weijun’s comment’s CSIE link for a general integral representation of $\Gamma^{(n)}(x)$. Therefore:
$$(-1)^{n+1}\text I_n=\lim_{u\to\infty}\ln(u)^{n+1}\sum_{m=0}^\infty \frac{(-1)^m}{m!}\text E_{-n}(-(m+1)\ln(u)) $$
A series about $x=1$ is simplified as an $x!$ series. If  one estimates by truncating the series and estimating $t$, then the $n,m$ upper bounds must be much greater than $t$:
$$x! = \lim_{t\to\infty} \sum_{n=0}^\infty\sum_{m=0}^\infty\frac{(-1)^{m+1}t^{n+1}x^n\text E_{-n}(-t(m+1))}{m!n!}$$
shown here with interchangeable sums. Similarly using the CSIE link, a series about $x=a$:
$$\boxed{\Gamma(x)=\lim_{t\to\infty}\sum_{m,n=0}^\infty\frac{(-1)^{m+1}t^{n+1}\text E_{-n}(-t(m+a))(x-a)^n}{m!n!}}$$
shown here.
A: Let $n\ge 1$ be an integer. Then the following holds true:
\begin{equation}
\lim\limits_{\xi \rightarrow 0}\frac{d^n}{d \xi^n} \left( \Gamma[\xi] - \frac{1}{\xi} \right) = 
\frac{1}{n+1} \cdot
\sum\limits_{1\cdot m_1+2\cdot m_2+\cdots+ (n+1) \cdot m_{n+1} = n+1}
\frac{(n+1)!}{\prod\limits_{j=1}^{n+1} m_j!} \cdot
\prod\limits_{j=1}^{n+1} \left( \frac{\psi^{(j-1)}(1)}{j!}\right)^{m_j}
\end{equation}
The Mathematica code below verifies the result for $n \le 10$.
In[3082]:= Clear[f]; xi =.; n = 9;
(*From definition*)
ex0 = Normal[Series[D[ Gamma[xi] - 1/xi, {xi, n}], {xi, 0, 0}]];


(*Using d'Hospital rule*)
ex = (D[ Exp[f[xi]], {xi, n}] + 
      xi/(n + 1) D[Exp[f[xi]], {xi, n + 1}]) /. 
    Exp[f[xi]] :> Gamma[xi] /. 
   Derivative[j_][f][xi] :> mPolyGamma[j - 1, xi];

ex = ex /. 
   mPolyGamma[m_, xi_] :> mPolyGamma[m, xi + 1] - (-1)^m m!/xi^(m + 1);

ex = Simplify[ex] /. xi Gamma[xi] :> Gamma[1 + xi] /. xi :> 0;

(*conjecture for generic n*)
ex1 = 1/(n + 1) Total[
    Flatten[DeleteCases[
      Flatten[Table[
        With[{mn = (n + 1 - Sum[eta m[eta], {eta, 1, n}])/(n + 1)}, 
         If[Element[mn, 
            Integers], {(n + 1)!/
             Product[ 
              m[xi]!, {xi, 1, 
               n + 1}] Product[(mPolyGamma[j - 1, 1]/j!)^m[j], {j, 1, 
               n + 1}]}, {}] /. m[n + 1] :> mn], 
        Evaluate[
         Sequence @@ 
          Table[{m[xi], 0, 
            Floor[(n + 1 - Sum[eta m[eta], {eta, 1, xi - 1}])/
             xi]}, {xi, 1, n}]]], n - 1], {}]]];
(ex - ex1) // Simplify

(ex0 - (ex /. mPolyGamma[j_, 1] :> PolyGamma[j, 1])) // Simplify

Out[3088]= 0

Out[3089]= 0

A: https://dlmf.nist.gov/5.7
has this:
$\dfrac1{\Gamma(z)}
=\sum_{k=1}^{\infty} c_kz^k
$
where
$c_1=1,
c_2=\gamma$,
and
$(k-1)c_k
=\gamma c_{k-1}-\sum_{j=2}^{k-1}\zeta(j)c_{k-j}
$.
It also has
$\ln(\Gamma(1+z))
=-\ln(1+z)+z(1-\gamma)
+\sum_{k=2}^{\infty}(-1)^k(\zeta(k)-1)\dfrac{z^k}{k}
$.
