Asymptotic expansion of Polygamma functions According to Wikipedia, the log-Gamma and Polygamma functions have the following asymptotic behaviour on the real line for $x\to\infty$:
$$\ln\Gamma(x) = (x - \tfrac{1}{2}) \ln(x) - x + \tfrac{1}{2}\ln (2\pi) + \sum_{k=1}^\infty \frac{B_{k+1}}{k(k+1)x^{k}} \\
\psi^{(0)}(x) \sim \ln(x) - \sum_{k=1}^\infty \frac{B_k}{kx^k}\\
\psi^{(m)}(x) \sim (-1)^{m+1} \sum_{k=0}^\infty \frac{(k+m-1)}{k!} \frac{B_k}{x^{k+m}}\quad\text{for}\ m>0$$
The second one can be obtained from the first one by term-wise differentiation, and the third one from the second one by iterated term-wise iteration.
The problem is that term-wise differention of asymptotic power series such as this is not valid in general. The book Asymptotic Expansions by A. Erdelyi mentions that one may perform term-wise differentiation of asymptotic power series like this one the function being expanded is a complex-valued function that is holomorphic on a suitably-shaped set (which $\ln\Gamma$ is). The proof uses Cauchy's integral formula and requires the power series expansion to hold uniformly.
However, I am working on a computer-aided proof and for technical reasons, I would like to avoid talking about the complex-valued $\ln\Gamma$ function and uniform asymptotic expansions if possible.
Is there any other way of obtaining the above results directly for the real $\psi^{(m)}$ functions? I already have the result
$$\ln\Gamma(x) = (x - \tfrac{1}{2})\ln(x) - s + \tfrac{1}{2}\ln(2\pi) + \sum_{k=1}^m \frac{B_{k+1}}{k(k+1)x^k}-\\\hskip-1em\frac{1}{m+1} \int_0^\infty \frac{B_{n+1}([t])}{(t+x)^{n+1}}\,\text{d}t$$
and, since the integral is ${\!}\in O(x^{-m-1})$, the above asymptotic expansion for $\ln\Gamma$, so I can use that – but is there any other way to get the other two expansions than the one described above?
One way that I thought about is somehow directly estimating the growth of the derivatives of the integral, but I did not get very far with that.
 A: Thanks for all the suggestions. Unfortunately, none of them struck me as straightforward and easy enough for my purposes and I eventually ended up using complex analysis after all. I'd like to show the solution I used in the end. It's quite similar to Erdelyi's proof I mentioned before.
First of all, I had already available the following inequality for all complex numbers $s$ with $\Re (s) > 0$ from Bruce Berndt's Rudiments of the Gamma Function (p. 13, Theorem 3):
$$\ln\Gamma(s) = (s - \tfrac{1}{2})\ln s - s + \tfrac{1}{2}\ln(2\pi) + \sum_{k=1}^n \frac{B_{2k}}{(2k-1)\,2k\, s^{2k-1}} - R_n(s)\\[8mm]
R_n(s) = \frac{1}{2n+1} \int_0^\infty\! \frac{B_{2n+1}([x])}{(x+s)^{2n+1}}\,\text{d}x$$
Also following Berndt's reasoning, it is easy to show that there exists some $c\in\mathbb{R}$ such that $|R_n(s)| \leq c \,\Re(s)^{-(2n+1)}$ for $R_n(s)$ for all $s$ with $\Re(s) > 0$.
If we let
$$F(s) := (s-\tfrac{1}{2})\ln s - s + \tfrac{1}{2}\ln(2\pi) + \sum_{k=1}^n\frac{B_{2k}}{(2k-1)\,2k\,s^{2k-1}}$$
then we have 
$$\ln\Gamma(s) - F(s) = R_n(s)$$
and since the left-hand side is clearly holomorphic on the positive real half-space, so is the right-hand side. Therefore, for any $x\in\mathbb{R}$, we can apply Cauchy's integral formula and integrate along a circle centred around $x$ with radius $\tfrac{x}{2}$:
$$\psi^{(m-1)}(x) - F^{(m)}(x) = R_n^{(m)}(x) = 
\frac{m!}{2\pi i} \oint\limits_{|z-x| = \tfrac{x}{2}} \frac{R_n(z)}{(z-x)^{m+1}}\,\text{d}z$$
The standard estimate for the integral gives us:
$$|\psi^{(m-1)}(x) - F^{(m)}(x)| \leq \frac{m!}{2\pi} \frac{2\pi\tfrac{x}{2} }{(\tfrac{x}{2})^{m+1}} \sup\limits_{|z-x|=\tfrac{x}{2}} |R_n(z)| \leq\\ m!\frac{2^m}{x^m}  \sup\limits_{|z-x|=\tfrac{x}{2}} c\,\Re(z)^{-(2n+1)}\leq\\
m!\frac{2^m}{x^m} c \cdot (\tfrac{x}{2})^{-(2n+1)} \in O\left(x^{-(2n+m+1)}\right)$$
which shows that $F^{(m)}$ is a valid asymptotic expansion of $\psi^{(m-1)}$.
A: Just a point of interest:
$$\psi^{(0)}(x)=c+\sum_{n=1}^x\frac1n$$
For a fixed constant $c$ and $x\in\mathbb N$.  More generally,
$$\psi^{(m)}(x)=c_m+(-1)^m\sum_{n=1}^x\frac1{n^{m+1}}$$
And now just tackle with the Euler-Maclaurin formula.
$$\sum_{n=1}^x\frac1{n^{m+1}}=\zeta(m+1)-\frac1m+\int_1^x\frac1{t^{m+1}}\ dt+\frac1{2x^{m+1}}-\frac{m+1}{12x^{m+2}}+\mathcal O\left(\frac1{x^{m+4}}\right)$$
