It is not hard to show, by induction on $m\in\mathbb N$, that there exist a sequence $(B_n)_{n\geq0}$ of rational numbers such that
$$\sum_{k=1}^nk^m=\frac1{m+1}\sum_{k=0}^m\binom{m+1}kB_k\,n^{m+1-k},\ \style{font-family:inherit;}{\text{for all}}\ n\geq1\,.$$
These are the Bernoulli numbers of the second kind, that is, $B_1=+1/2$. A corollary of the proof (by induction) of the fact above is a recurrence formula for such numbers $B_n$, which are known as Bernoulli numbers:
$$\sum_{k=0}^{n-1}\binom nkB_k=0\,.\tag{$\ast$}$$
On the other hand, there are a number of explicit formulae for $B_n$, which are obtained using the equality $\frac x{e^x-1}=\sum_{n=0}^\infty B_n\,\frac{x^n}{n!}$; see here and here. For example:
$$B_n=\sum_{k=0}^n\frac1{k+1}\sum_{r=0}^k(-1)^r\binom krr^n\,.\tag{$\ast\ast$}$$
Is there some clever way to manipulate the recurrence $(\ast)$ in order to obtain $(\ast\ast)$?
Best Answer
[20-08] I will add something that sheds light on why the latter works out nicely. In what follows we provide a proof of Faulhaber's Formula using the operation of binomial convolution on sequences of real numbers (in fact, we needn't worry about the whole real numbers, but merely the rational numbers). This operation arises naturally when multiplying (formal) exponential generating series of the form $$\sum_{n\geqslant 0}\gamma_n \frac{z^n}{n!}$$ and collecting powers of $z$. The problem of finding $1^n+2^n+\cdots+m^n$ motivates us to look into each of the sequences $1,j,j^2,j^3,\ldots$ that have the EGF $e^{jz}$. Thus the sum has EGF $$\sum_{j=1}^m e^{jz}=e^z\frac{e^{mz}-1}{e^{z}-1}=\frac{e^{mz}-1}z \frac{(-z)}{e^{-z}-1}$$
It is easy to obtain explicit generating sequences for $\dfrac{e^{mz}-1}z$; indeed $$\frac{e^{mz}-1}z=\sum_{n\geqslant 0}\frac{m^{n+1}}{n+1}\frac{z^n}{n!}$$
In particular when $m=1$ we get the Harmonic numbers generate $$H(z)=\frac{e^{z}-1}z$$
The Bernoulli numbers arise naturally as the inverses of the Harmonic numbers, in the sense that $$B(z)=\frac{z}{e^z-1}=\sum_{n\geqslant 0}B_n\frac{z^n}{n!}$$
The first equation then gives Faulhaber's Formula immediately: we have $$\sum_{j=1}^m e^{jz}=e^z\frac{e^{mz}-1}{e^{z}-1}=B(-z)\frac{e^{mz}-1}z $$
Equating coefficients and using the convolution gives $$\sum_{j=1}^m j^n =\sum_{k=0}^n\binom nk (-1)^k B_k\frac{m^{n-k+1}}{n-k+1}$$
This is the shortest proof I know of. We work things formally, without worrying about convergence. We use the fact that a sequence with non-zero first term $a_0$ admits a unique inverse (as is the case of the Harmonic numbers). This is proven by observing the convolution gives a recurrence for the terms of this inverse.
The above is very similar to what will follow. The first lemma gives $$H(z)e^{jz}=\frac{e^{(j+1)z}-1}z-\frac{e^{jz}-1}z$$ and the second gives $$e^z\frac{e^{mz}-1}z=\frac{e^{(m+1)z}-1}z-\frac{e^{z}-1}z$$ Telescoping gives $$H(z)\sum\limits_{j = 1}^m {{e^{jz}}} = \frac{{{e^{(m + 1)z}} - {e^z}}}{z} = {e^z}\frac{{{e^{mz}} - 1}}{z}$$
Thus $$\sum\limits_{j = 1}^m {{e^{jz}}} = \frac{{{e^{(m + 1)z}} - {e^z}}}{z} = {e^z}B\left( z \right)\frac{{{e^{mz}} - 1}}{z}$$ and using $B(-z)=e^zB(z)$ $-$ this $1\star B=\mu B$ which was quite ugly to prove without EGFs $-$ finishes things to give Faulhaber's Formula $$\sum\limits_{j = 1}^m {{e^{jz}}} = B\left( { - z} \right)\frac{{{e^{mz}} - 1}}{z}$$
We're ready to prove the
P By the above $$\begin{align}{\bf x}\star {\bf 1}&={\bf y}\\({\bf x}\star {\bf 1})\star \mu&={\bf y}\star \mu\\{\bf x}\star ({\bf 1}\star\mu)&={\bf y}\star \mu\\{\bf x}\star {\rm id }&={\bf y}\star \mu\\{\bf x}&={\bf y}\star \mu\end{align}$$
and the other direction is analogous.
[7-28] I have decided to define the Bernoulli numbers as the inverse of the Harmonic numbers, since this puts in evidence their importance when proving Faulhaber's formula. I will keep looking for a proof of the closed form formula.
This is well-defined since inverses are unique. A few values are $$\left(1,-\dfrac 1 2,\dfrac 1 6,0,-\dfrac 1 {30},0,\dfrac1{42},0,-\dfrac{1}{30},\ldots\right)$$
P (Credits to Rob) From ${\bf B}\star {\bf 1}={\bf B}+[n=1]$ we obtain by inversion that ${\bf B} ={\bf B}\star \mu+[n=1]\star \mu$. But $[n=1]\star \mu=-(-1)^nn$; so ${\bf B}\star \mu=B+\mu{\bf N}$. Since ${\bf x}\star \mu {\bf y}=\mu(\mu{\bf x}\star {\bf y})$ we obtain $\mu{\bf B}\star {\bf 1}=\mu{\bf B}+{\bf N}$. Now we write things explicitly $$\sum_{k=0}^n\binom{n+1}{k}B_k=[n=0]\;\;,\;\;\sum_{k=0}^{m-1}\binom mk(-1)^kB_k=m$$
Thus by $m\mapsto n+1$ $$\sum_{k=0}^n\binom{n+1}{k}B_k=[n=0]\;\;,\;\;\sum_{k=0}^{n}\binom {n+1}k(-1)^kB_k=n+1$$
Thus $$n+1-[n=0]=\sum_{k=0}^n\left((-1)^k-1\right)\binom{n+1}kB_k$$ and $n\mapsto 2n+1$ gives $$2n+2=\sum_{k=0}^{2n+1}\left((-1)^k-1\right)\binom{2n+2}kB_k$$ but the $k=1$ term is $-2\binom{2n+2}1\left(-\frac1 2\right)=2n+2$ so, since even terms vanish we get $$0=-2\sum_{k=1}^{n}\binom{2n+2}{2k+1}B_{2k+1}$$ and induction does the rest. $\blacktriangle$
First, we have:
P This follows from direct calculation and use of the binomial theorem.
P This should also be fairly obvious by integrating the binomial expansion and finding the appropriate constant of integration.
Now we can prove
P We have that $$\sum\limits_{j = 1}^m {{{\left( {{\bf{H}} \star {j^k}} \right)}_n}} = \left( {1 \star \frac{{{m^{k + 1}}}}{{k + 1}}} \right)$$
By the distributive property $$\left( {1 \star \frac{{{m^{k + 1}}}}{{k + 1}}} \right) = {\bf{H}} \star \sum\limits_{j = 1}^m {{j^n}} $$ and after multiplication by $\bf B$, this is $$\eqalign{ & {\bf{B}} \star \left( {1 \star \frac{{{m^{k + 1}}}}{{k + 1}}} \right) = \sum\limits_{j = 1}^m {{j^n}} \cr & \left( {\mu {\bf{B}} \star \frac{{{m^{k + 1}}}}{{k + 1}}} \right) = \sum\limits_{j = 1}^m {{j^n}} \cr} $$
as desired. $\blacktriangle$