Closed form for $1^k + ... + n^k$ (generalized Harmonic number)

This question must have been asked, it's just very hard to search for such questions.

I'm looking for the cleanest method I can find for getting a closed form formula for $\sum_{i=1}^n i^k$

Wikipedia provides https://en.wikipedia.org/wiki/Faulhaber%27s_formula which has a tantalising "There is also a similar (but somehow simpler) expression:..." paragraph but fails to follow up.

http://www.maa.org/press/periodicals/convergence/sums-of-powers-of-positive-integers-conclusion contains a thorough "through the ages" expose, the last two entries being Pascal and Bernoulli.

However, https://www.youtube.com/watch?v=8nUZaVCLgqA seems to contain a fresh approach that I can't find documented anywhere else. However, I find the video hard to follow.

And maybe there is some technique that is cleaner still...

I understand that "cleanest" is maybe subjective and hence it is an imperfect question, but it would be interesting to see the various alternatives (and surely there cannot be many) battle it out.

[2016-08-06]: Connection with Stirling Numbers of the second kind added.

Here we derive a method based upon formal power series to obtain a closed formula for \begin{align*} S_k(n):=\sum_{j=1}^n j^k\qquad\qquad n,k\geq 1 \end{align*}

In order to do so we encode sequences $(a_j)_{j\geq 0}$ by generating functions $A(z)=\sum_{j=0}^\infty a_j z^j$.

Constant sequence $(1)_{j\geq 0}$

We start with the constant sequence $(1)_{j\geq 0}$

\begin{align*} (1)_{j\geq 0}=(1,1,1,\ldots) \quad \rightarrow \quad \frac{1}{1-z}&=\sum_{j=0}^{\infty}z^j\tag{1}\\ &=1+z+z^2+\cdots \end{align*}

We see the constant sequence is encoded by the geometric power series.

$$ $$

Getting $k$-th powers $(j^k)_{j\geq 0}$

Differentiating a power series and multiplication with $z$ results in \begin{align*} zD_z \sum_{j=0}^\infty a_j z^j = \sum_{j=0}^\infty ja_jz^j \end{align*} Here we denote with $D_z:=\frac{d}{dz}$ the differential operator. If we apply in (1) the operator $zD_z$ successively $k$ times, we get $k$-th powers of $j$.

\begin{align*} (j^k)_{j\geq 0}=(0^k,1^k,2^k,\ldots)\quad\rightarrow\quad (zD_z)^k\frac{1}{1-z}&=\sum_{j=0}^\infty j^kz^j\tag{2}\\ &=0^k+1^kz+2^kz^2+\cdots \end{align*}

$$ $$

Summing up $k$-th powers $\left(\sum_{j=0}^n j^k\right)_{j\geq 0}$

A nice fact is that summing up elements is encoded in formal power series by multiplication with $\frac{1}{1-z}$. This is due to the Cauchy product formula \begin{align*} \frac{1}{1-z}\sum_{j=0}^\infty a_jz^j&=\left(\sum_{l=0}^\infty z^l\right)\left(\sum_{j=0}^\infty a_jz^j\right)=\sum_{m=0}^\infty\left(\sum_{l=0}^m a_l\right)z^m \end{align*}

So, multiplying (2) with the operator $\frac{1}{1-z}$ provides a generating function for the sum of the $k$-th powers of the first $n$ numbers.

\begin{align*} \left(\sum_{j=0}^n j^k\right)_{n\geq 0}\quad \rightarrow \quad \frac{1}{1-z}(zD_z)^k\frac{1}{1-z}&=\sum_{n=0}^\infty\left(\sum_{j=0}^n j^k\right)z^n\\ &=0^k+\left(0^k+1^k\right)z+\left(0^k+1^k+2^k\right)z^2+\cdots \end{align*}

It's convenient to use the coefficient of operator $[z^n]$ to denote the coefficient of $z^n$. We summarize this as

Method 1:

The sum of the $k$-th powers of numbers $1$ to $n$ is given by \begin{align*} S_k(n)=\sum_{j=1}^nj^k=[z^n]\frac{1}{1-z}(zD_z)^k\frac{1}{1-z}\tag{3} \end{align*}

A small example with $k=2$.

Example: $S_2(n)$

\begin{align*} S_2(n)=\sum_{j=0}^n j^2&=[z^n]\frac{1}{1-z}(zD_z)^2\frac{1}{1-z}\tag{4}\\ &=[z^n]\frac{z(1+z)}{(1-z)^4}\\ &=[z^n](z+z^2)\sum_{j=0}^{\infty}\binom{-4}{j}(-z)^{j}\tag{5}\\ &=\left([z^{n-1}]+[z^{n-2}]\right)\sum_{j=0}^{\infty}\binom{j+3}{3}z^j\tag{6}\\ &=\binom{n+2}{3}+\binom{n+1}{7}\\ &=\frac{1}{6}n(n+1)(2n+1) \end{align*}

Comment:

• In (4) we apply the operator $\frac{1}{1-z}(zD_z)^2$ to $\frac{1}{1-z}$

• In (5) we use the binomial series expansion

• In (6) we use the linearity of the coefficient of operator, apply the formula

\begin{align*} [z^{p-q}]A(z)=[z^p]z^qA(z) \end{align*} and use the binomial identity \begin{align*} \binom{-p}{q}=\binom{p+q-1}{p-1}(-1)^q \end{align*}

Connection with Stirling Numbers of the second kind:

The operator $(zD_z)^k$ can be expressed with Stirling Numbers of the second kind ${n\brace k}$: (see e.g. this paper) \begin{align*} \left(zD_z\right)^k=\sum_{j=1}^k{k\brace j}z^jD_z^j\tag{7} \end{align*}

Applying this formula to (3) we obtain \begin{align*} [z^n]&\frac{1}{1-z}\left(zD_z\right)^k\frac{1}{1-z}\\ &=[z^n]\frac{1}{1-z}\sum_{j=1}^k{k\brace j}z^jD_z^j\frac{1}{1-z}\\ &=[z^n]\frac{1}{1-z}\sum_{j=1}^k{k\brace j}j!\frac{z^j}{(1-z)^{j+1}}\tag{8}\\ &=[z^n]\frac{1}{1-z}\sum_{j=1}^k{k\brace j}j!\sum_{l=0}^\infty\binom{-(j+1)}{l}(-1)^lz^{j+l}\tag{9}\\ &=[z^n]\frac{1}{1-z}\sum_{j=1}^k{k\brace j}j!\sum_{l=0}^\infty\binom{j+l}{l}z^{j+l}\tag{10}\\ &=\sum_{j=1}^k{k\brace j}j!\sum_{l=0}^{n-j}\binom{j+l}{l}[z^{n-j-l}]\frac{1}{1-z}\tag{11}\\ &=\sum_{j=1}^k{k\brace j}j!\sum_{l=0}^{n-j}\binom{j+l}{l}\tag{12}\\ &=\sum_{j=1}^k{k\brace j}j!\binom{n+1}{j+1}\tag{13}\\ &=\sum_{j=1}^k{k\brace j}\frac{(n+1)^{\underline{j+1}}}{j+1}\tag{14}\\ \end{align*}

Comment:

• In (8) we differentiate $j$ times and get $D_z^j(1-z)^{-1}=j!(1-z)^{j+1}$

• In (9) we use the binomial series expansion

• In (10) we use $\binom{-p}{q}=\binom{p+q-1}{p-1}(-1)^q$ again

• In (11) we do some rearrangements and use $[z^{p-q}]A(z)=[z^p]z^qA(z)$ again

• In (12) we observe the contribution of the coefficient of the geometric series is always one.

• In (13) we use the binomial identity $\sum_{l=0}^{n-j}\binom{j+l}{l}=\sum_{l=j}^{n}\binom{l}{j}=\binom{n+1}{j+1}$

• In (14) we do some simplifications and use Don Knuths falling factorial power notation: $$n^{\underline{j}}=\frac{n!}{(n-j)!}$$

We summarize this as

Method 2:

The sum of the $k$-th powers of numbers $1$ to $n$ is given by \begin{align*} S_k(n)=\sum_{j=1}^nj^k=\sum_{j=1}^k{k\brace j}\frac{(n+1)^{\underline{j+1}}}{j+1} \end{align*}

A small example with $k=2$.

Example $S_2(n)$

\begin{align*} S_2(n)=\sum_{j=1}^nj^2&=\sum_{j=1}^2{2\brace j}\frac{(n+1)^{\underline{j+1}}}{j+1}\\ &={2\brace 1}\frac{(n+1)^{\underline{2}}}{2}+{2\brace 2}\frac{(n+1)^{\underline{3}}}{3}\\ &=\frac{1}{2}(n+1)n+\frac{1}{3}(n+1)n(n-1)\\ &=\frac{1}{6}n(n+1)(2n+1) \end{align*}

Hint: Note that generalised harmonic numbers mean the sum of reciprocal values of $j^k$ \begin{align*} H_{n,k}=\sum_{j=1}^nj^{-k} \end{align*}

Put $$ S_\ell(n):=\sum_{i=1}^{n}i^{\ell} $$ and consider the identity $$ (n+1)^{k+1}-1=\sum_{i=1}^{n}[(i+1)^{k+1}-i^{k+1}]=(k+1)S_k(n)+\binom{k+1}{2}S_{k-1}(n)+\ldots+\\+\binom{k+1}{k}S_1(n)+n. $$ Then, if you know the sums for $\ell=1,\ldots,k-1$, you can calculate $S_k(n)$ with the above recursive relation. For example, $$ (n+1)^2-1=2S_1(n)+n\Rightarrow S_1(n)=\frac12[(n+1)^2-1-n]=\frac12n(n+1),\\ (n+1)^3-1=3 S_2(n)+3S_1(n)+n\Rightarrow S_2(n)=\frac13[n^3+3n^2+3n-n-3S_1(n)]\Rightarrow\\ \Rightarrow S_2(n)=\frac{n^3}{3}+\frac{n^2}{2}+\frac{n}{6} $$

Having posted a much similar question recently, I'm inclined to provide you with my somewhat unique answer.

$$a_p=1-p\int_0^1f(t,p-1)dt,\quad f(x,0)=x$$

$$f(x,p)=a_px+p\int_0^xf(t,p-1)dt$$

$$f(x,p)=\sum_{k=1}^xk^p$$

Easy enough to see that

$$x=\sum_{k=1}^xk^0$$

Also relatively easy to find that

$$a_1=1-\int_0^1t\ dt=\frac12\implies f(x,1)=\frac12x+\int_0^xt\ dt=\frac12x+\frac12x^2=\sum_{k=1}^xk^1$$

And you can keep going with this. Requires very little calculus skill to do.