The gamma function can be defined in a few different ways, the most well known possibly being Euler's second integral definition $$\Gamma(x) = \int_0^\infty t^{x-1}e^{-t}dt, \quad x>0$$ another being Euler's infinite product definition $$\Gamma(x) = \frac{1}{x} \prod_{n=0}^\infty \frac{\left(1+\frac{1}{n}\right)^x}{\left(1+\frac{x}{n}\right)},\quad x>0.$$ It is a standard exercise to show that these agree over the positive real numbers. It is an interesting question to ask whether or not this function is analytic and if so over what intervals? i.e for which $x_0>0$ can we write it as an infinite series of the form $\Gamma(x) = \sum_{n=0}^\infty a_n (x-x_0)^n$ for some real constants $a_n$ such that it converges in an open neighbourhood of $x_0?$ The two definitions of $\Gamma(x)$ given above are well defined for complex $x$ as well (with some restrictions). Since analyticity is a stronger property in the complex plane, it is easier to prove that $\Gamma(z), \Re(z)>0,$ is holomorphic everywhere in it's domain then use analytic continuation to extend elsewhere in the plane. We then know that $\Gamma(x)$ is analytic over $x>0$, since we can just restrict to the positive reals.
There does not seem to be a proof of this theorem that doesn't use complex analysis in someway. To avoid the singularities at negative integer arguments, and to allow us to use the integral definition, I am only interested in proving analyticity over the positive reals. How can we prove the gamma function is analytic over the positive real numbers without using the theorems of complex analysis?
My attempt
Here is my attempt with gaps. Consider Prym's representation $$\Gamma(x)=\sum_{n=0}^\infty \frac{(-1)^n}{n!}\frac{1}{x+n}+\int_1^\infty t^{x-1}e^{-t}dt.$$
We'll show that both parts of this sum is analytic, thus their sum is analytic and $\Gamma(x)$ is analytic.
The infinite sum is analytic because [...]
To prove the right hand is analytic, we'll use Bernstein's Theorem from Krantz and Park's book A Primer of Real Analytic Functions
Theorem: Let $f$ be a smooth function $(C^\infty)$ on an open interval $I\subset \mathbb{R}$. If $f$ and all it's derivatives are non-negative on the entire interval $I$ then $f$ is real analytic on $I$.
The smoothness property follows from the fact that $t^{x-1}e^{-t}$ is smooth over $(x,t)\in \left[0,\infty\right) \times \left[1,\infty\right)$. Then an application of Leibniz theorem gives that \begin{align*} \frac{d^n}{dx^n}\int_1^\infty t^{x-1}e^{-t}dt &= \int_1^\infty \frac{d^n}{dz^n} t^{x-1}e^{-t}dt \\ &= \int_1^\infty \ln(t)^n t^{x-1}e^{-t}dt \end{align*}
now since the integrand is positive for $x>0$ and $1<t<\infty$, the integral must also be positive. Thus $\int_1^\infty t^{x-1}e^{-t}dt$ is smooth and has positive derivatives for all $x>0$ and so $\int_1^\infty t^{x-1}e^{-t}dt$ is analytic over $x>0$. Thus $\Gamma(x)$ is analytic over $x>0$. $\blacksquare$
So it remains to show that the infinite sum defines an analytic function (does the infinite sum of analytic functions converge to an analytic function?) and it would be nice to be less handwavy with the smoothness property.
$$\Gamma(x)=\int^\infty_0 t^{x-1}e^{-t}\,dt=\int^1_0t^{x-1}e^{-t}\,dt+\int^\infty_1t^{x-1}e^{-t}\,dt $$
First the easy part: For any $x\in\mathbb{R}$, monotone convergence yields \begin{align} \int^\infty_1t^{x-1}e^{-t}\,dt&=\int^\infty_1\Big(\sum^\infty_{n=0}\frac{1}{n!}\log^n(t)x^n\Big)t^{-1}e^{-t}\,dt\\ &=\sum^\infty_{n=0}\frac{1}{n!}\Big(\int^\infty_1\log^n(t)\,t^{-1}e^{-t}\,dt\Big)x^n \end{align} Thus, $G_2(x)=\int^\infty_1t^{x-1}e^{-t}\,dt$ is analytic on $\mathbb{R}$.
Dominated convergence, yields for $x>0$ \begin{align} \int^1_0 t^{x-1}e^{-t}\,dt&=\int^1_0 t^{x-1}\Big(\sum^\infty_{n=0}\frac{(-t)^n}{n!}\Big)\,dt\\ &=\sum^\infty_{n=0}\frac{(-1)^n}{n!}\int^1_0 t^{x+n-1}\,dt\\ &=\sum^\infty_{n=0}\frac{(-1)^n}{n!}\frac{1}{x+n}\tag{1}\label{one} \end{align} It remains to show that $G_1(x)=\int^1_0 t^{x-1}e^{-t}\,dt$ is analytic on $(0,\infty)$, The idea is to split the sum in \eqref{one} as $\sum^N_{n=0} +\sum^\infty_{n=N+1}$ so that for all $y$ inside a given neighborhood $B(x;R)$, $y<n$ whenever $n>N$. I leave the details to the OP.