Solving an integral coming from Perron's formula In analytic number theory, Perron's formula says that 
$$
\sum_{1 \leq k < n} a_k + \frac{1}{2}a_n = \int_{c - i\infty}^{c+i\infty} f(s)\frac{n^s}{s}ds, 
$$
where $f(s) = \sum_{k \geq 1} a_k/k^s$ and $c$ is greater than the abscissa of convergence 
of $f(s)$.  My question is whether Perron's formula can be used as the first step in a derivation of the formula for 
a sum of powers of the first $n$ integers. 
Specifically, for any positive integer $m$, Perron's formula 
for $f(s) = \zeta(s-m) = \sum_{k \geq 1} k^m/k^s$ says 
$$\sum_{k=1}^n k^m = \frac{1}{2}n^m+\frac{1}{2\pi i}\int_{c-i\infty}^{c+i\infty}\zeta(s-m)\frac{n^s}{s}ds$$ 
for any $c>m+1$.  The integrand $\zeta(s-m)n^s/s$ has poles at $s = 0$ and $s = m+1$.
If we formally shift the line of integration to the left across the pole at $s = m+1$, then a formal use of the residue theorem gives us 
$$
\sum_{k=1}^n k^m = \frac{1}{m+1}n^{m+1} + \frac{1}{2}n^m+\frac{1}{2\pi i}\int_{b-i\infty}^{b+i\infty}\zeta(s-m)\frac{n^s}{s}ds
$$ 
for $b < m+1$.  When $m-1/2 < b < m+1$, I can justify this shift (i.e., if I truncate the integral on the top and bottom and make a rectangle 
with sides along $x = c$ and $x = b$, the integrals along the top and bottom go to 0 as the height of the rectangle goes to $\pm \infty$).
The two powers of $n$ on the right side are exactly the first two dominant terms in the standard formula for power sums 
in terms of Bernoulli numbers: 
$$\sum_{k=1}^n k^m =  \sum_{j=0}^{m} \dbinom{m}{j}\frac{B_{m-j}}{j+1}n^{j+1} = 
\frac{1}{m+1}n^{m+1}+\frac{1}{2}n^m + \sum_{j=0}^{m-2} \dbinom{m}{j}\frac{B_{m-j}}{j+1}n^{j+1}. 
$$
What I'd like to know is if anyone sees a way to extract all lower order terms in this standard formula from the integral 
along the line ${\rm Re}(s) = b$ when $m-1/2 < b  < m+1$. That is, could one show
$$
\frac{1}{2\pi i}\int_{b-i\infty}^{b+i\infty}\zeta(s-m)\frac{n^s}{s}ds = \sum_{j=0}^{m-2} \dbinom{m}{j}\frac{B_{m-j}}{j+1}n^{j+1}
$$
by some method that is independent of knowledge of the power sum formula? 
For concreteness, the first four cases of the standard power sum formula are 
\begin{eqnarray*}
\sum_{k=1}^n k &=& \frac{1}{2}n^2+\frac{1}{2}n, \\
\sum_{k=1}^n k^2 &=& \frac{1}{3}n^3+\frac{1}{2}n^2 + \frac{1}{6}n, \\
\sum_{k=1}^n k^3 &=& \frac{1}{4}n^4+\frac{1}{2}n^3 + \frac{1}{4}n^2, \\ 
\sum_{k=1}^n k^4 &=& \frac{1}{5}n^5+\frac{1}{2}n^4 + \frac{1}{3}n^3-\frac{1}{30}n.
\end{eqnarray*}
So when $m=1$ the integral is 0, when $m=2$ the integral is $\frac{1}{6}n$, when 
$m = 3$ the integral is $\frac{1}{4}n^2$, and when $m = 4$ the integral is $\frac{1}{3}n^3 - \frac{1}{30}{n}$.
Is there a way to evaluate the integral directly to recover these computations (thereby leading to an alternate proof of the 
power sum formula)?  Or is there at least a way to see that the integral is $O(n^{m-1})$ as $n \rightarrow \infty$?
The obvious thing to try is to push the contour further to the left, past more poles.  Even if that could be justified (is it valid?), the only pole remaining in $\zeta(s-m)n^s/s$ is at $s = 0$, where the residue is $\zeta(-m) = -B_{m+1}/(m+1)$, which doesn't involve $n$ and is not even a term in the power sum formula.  Without poles to work with, I don't see a way to extract any lower order terms from the integral.
 A: I believe this was done (in a more general context) by P. Flajelet et alin 1995 (Mellin transforms and asymptotics, Harmonic sums). A thesis all about this is Marko Riedel's thesis at UBC (which you can google).
A: Very nice question. You are on the right track when you conjecture that Faulhaber's formula can be proved by Mellin summation. However Mellin-Perron will only give you the first two terms. You need to apply Mellin directly to a certain harmonic sum if you want more terms.
Re-write your sum like this:
$$S_m(n) = \sum_{k=1}^n k^m = n^m \sum_{k=1}^n \left(\frac{k}{n}\right)^m
= n^m + n^m \sum_{k=1}^n \left(1 - \frac{k}{n}\right)^m.$$
Now consider the harmonic sum
$$S(x) = \sum_{k\ge 1} H_m(kx)$$
where $H_m(x)$ is the Heaviside step function defined by
$$H_0(x) = \begin{cases} & 1&\text{if}\quad x\in[0,1] \\
& 0 &\text{otherwise}\end{cases}
\quad\text{and}\quad H_m(x) = (1-x)^m H_0(x)
\quad\text{when}\quad m\in\mathbb{Z}^+.$$
We see that $$S(1/n) = \sum_{k=1}^n \left(1 - \frac{k}{n}\right)^m,$$
i.e. the sum we are trying to compute.
But $S(x)$ is harmonic and may be evaluated by inverting its Mellin transform.
Recall the harmonic sum identity
$$\mathfrak{M}\left(\sum_{k\ge 1} \lambda_k g(\mu_k x);s\right) =
\left(\sum_{k\ge 1} \frac{\lambda_k}{\mu_k^s} \right) g^*(s)$$
where $g^*(s)$ is the Mellin transform of $g(x).$
In the present case we have
$$\lambda_k = 1, \quad \mu_k = k \quad \text{and} \quad
g(x) = H_m(x).$$
We need the Mellin transform $H_m^*(s)$ of $H_m(x)$ which is
$$\int_0^\infty H_m(x) x^{s-1} dx
= \frac{m!}{s(s+1)(s+2)\cdots(s+m)}$$
as is easily seen by induction using integration by parts.
It follows that the Mellin transform $Q(s)$ of $S(x)$ is given by
$$ Q(s) = \frac{m!}{s(s+1)(s+2)\cdots(s+m)} \zeta(s)
\quad\text{because}\quad \sum_{k\ge 1} \frac{\lambda_k}{\mu_k^s} = \zeta(s).$$
We thus obtain the Mellin inversion integral
$$ S(x) = \frac{1}{2\pi i} \int_{3/2-i\infty}^{3/2+i\infty} Q(s)/x^s ds$$
which we evaluate by shifting it to the left for an expansion about zero.
We have
$$\mathrm{Res}(Q(s)/x^s; s=1) = \frac{1}{m+1} \frac{1}{x}.$$
For the pole at $s=0$ we get
$$\mathrm{Res}(Q(s)/x^s; s=0) = -\frac{1}{2}.$$
For the poles at $s=-1, -2, \ldots -m$ we obtain
$$\mathrm{Res}(Q(s)/x^s; s=-q)
= \left.\frac{m!}{s(s+1)\cdots(s+q-1)(s+q+1)\cdots(s+m)}
\frac{\zeta(-q)}{x^s}\right|_{s=-q}
\\= - \frac{m!}{(-1)^q q!\times (m-q)!} \frac{B_{q+1}}{q+1} x^q
= (-1)^{q+1} {m\choose q}  \frac{B_{q+1}}{q+1} x^q.$$
It is readily seen that this formula correctly encapsulates the case of the trivial zeros of the zeta function canceling the pole from the fractional term. Now returning to the sum and using $S(1/n)$ we finally have
$$S_m(n) = n^m +
n^m \left(\frac{1}{m+1} n -\frac{1}{2}
+ \sum_{q=1}^m (-1)^{q+1} {m\choose q}  \frac{B_{q+1}}{q+1} \frac{1}{n^q} \right).$$
This is $$\frac{1}{m+1} n^{m+1} + \frac{1}{2} n^m
+ \sum_{q=1}^m (-1)^{q+1} {m\choose q}  \frac{B_{q+1}}{q+1} n^{m-q}.$$
Here we are confronted with a phenomenon that was noted by the OP namely that there is a spurious contribution from the pole at $s=-m$ which is a constant term that does not appear in any of the power sum formulas. So we should only shift the integral to $\Re(s) = -m+1/2$ and not $\Re(s) = -m-1/2.$ To properly justify this is likely to be a delicate matter since the zeta function outgrows the $O(s^{m+1})$ cancellation effect from the fractional term on vertical lines and the line integral does not vanish as we shift further to the left. I would conjecture that it attains some kind of minimum on  $\Re(s)\in(-m-1/2, -m+1/2).$ 
Without the non-contributing term we obtain
$$\frac{1}{m+1} n^{m+1} + \frac{1}{2} n^m
+ \sum_{q=1}^{m-1} (-1)^{q+1} {m\choose q}  \frac{B_{q+1}}{q+1} n^{m-q}$$
which some authors simplify as follows: we get for the sum term
$$\sum_{q=0}^{m-2} (-1)^{m-q} {m\choose m-q-1}  \frac{B_{m-q}}{m-q} n^{q+1} 
= \sum_{q=0}^{m-2} (-1)^{m-q} {m\choose q+1}  \frac{B_{m-q}}{m-q} n^{q+1}\\
= \sum_{q=0}^{m-2} (-1)^{m-q} {m\choose q}  \frac{B_{m-q}}{q+1} n^{q+1}.$$
Now we have $$(-1)^{m-(m-1)}{m\choose m-1} \frac{B_{m-(m-1)}}{(m-1)+1} n^{(m-1)+1}
= -B_1 n^m = - \left(-\frac{1}{2}\right) n^m = \frac{1}{2} n^m$$
and we can absorb the second term into the sum, getting
$$\sum_{q=0}^{m-1} (-1)^{m-q} {m\choose q}  \frac{B_{m-q}}{q+1} n^{q+1}.$$
Finally note that
$$(-1)^{m-m}{m\choose m} \frac{B_{m-m}}{m+1} n^{m+1} = \frac{1}{m+1} n^{m+1}$$
and we can absorb the first term as well, for the ultimate answer
$$\sum_{q=0}^m (-1)^{m-q} {m\choose q}  \frac{B_{m-q}}{q+1} n^{q+1}.$$
Some additional simplification is possible noting that the $B_{m-q}$ are zero when $m-q$ is not even with the exception of $B_1$, giving
$$\sum_{q=0}^m (-1)^{\delta_{q,m-1}} {m\choose q}  \frac{B_{m-q}}{q+1} n^{q+1}.$$
Here we have used the Bernoulli numbers with exponential generating function
$$\frac{t}{\exp t -1}.$$
