[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*}
