Finding the general/closed form of $\sum_{k=1}^n k^a$ I recently noticed the following:
$$\sum_{i=1}^ni=\frac{n(n+1)}{2}$$
$$\sum_{i=1}^{n}i^2 = \dfrac{n(n+1)(2n+1)}{6} = \Bigg(\sum_{i=1}^ni\Bigg) \cdot \frac{2n+1}{3}$$
$$\sum\limits_{i = 1}^n i^3 = \Bigg( \frac{n(n+1)}{2}\Bigg)^2 = \Bigg(\sum_{i=1}^ni\Bigg)^2$$
$$\sum_{i=1}^ni^4=\frac{n(n+1)(2n+1)(3n^2+3n-1)}{30} = \Bigg(\sum_{i=1}^ni^2\Bigg) \cdot \frac{3 n^2 + 3 n - 1}{5}$$
It seems that the sums build on each other in some interesting way, and I'm interested in finding a general form. A cursory investigation found that $$\sum_{k=0}^n k^a = H_n^{(-a)}$$ where $H_n^{(-a)}$ is the $n$th hyperharmonic number of $r$th order. However, I doubt this would help because this (a) is sort of circular given the definition of $H_n$, (b) has no closed form, and (c) doesn't reveal the multiplicative connection among the terms.
Is it possible to derive a closed form for $\sum_{k=0}^n k^a$? Why or why not? And, of course, if one exists, what is it?
 A: Let $$ S_n(p)=\sum_{k=1}^{n} k^p\qquad n, p\in\mathbb N ~~~~~\text{called Cavalieri sum of oder p}$$
Therefore, using the  Binomial formula we get
$$ (k+1)^p = k^p+ \sum_{i=0}^{p-1}\binom{p}{i} k^i$$
where $\binom{p}{i}= \frac{p!}{i!(p-i)!}$.
summing up both side yields, 
$$\sum_{k=1}^{n} (k+1)^p =\sum_{k=1}^{n} k^p+\sum_{i=0}^{p-1}\binom{p}{i} \sum_{k=1}^{n} k^i = S_n(p) +\sum_{i=0}^{p-1}\binom{p}{i} S_n(i) $$
However, $$\sum_{k=1}^{n} (k+1)^p = \sum_{k=2}^{n+1} k^p = S_{n+1}(p) -1 = S_n(p) +(n+1)^p -1$$
Hence finally we get the formula :
$$\color{red}{(n+1)^p -1  =\sum_{i=0}^{p-1}\binom{p}{i} S_n(i)} $$
or $$\color{red}{S_{n+1}(p) -1 =S_{n}(p) =\sum_{i=0}^{p}\binom{p}{i} S_n(i)} $$

From this it is possible to compute the sum for any $p\ge 1 $ in $ \mathbb N $.

A: So apparently, the "multiplicative" behavior I was looking for is almost certainly a coincidence. Perhaps what I'm seeing is commonality between terms, since the sum I'm looking for is defined by a summation published by Jacob Bernoulli in 1713:
$$\sum_{k=1}^n k^p = \frac{1}{p+1} \sum_{j=0}^p \binom{p+1}{j} B_j n^{p+1-j}$$
where $B_j$ is the $j$th Bernoulli number of the second kind.
So, for instance, if we had $\sum k^4$, we'd note that $$B_0 = 1, B_1 = \frac 12, B_2 = \frac 16, B_3 = 0, B_4 = -\frac {1}{30}$$
and then compute
$$\sum_{k=1}^n k^4 = \frac{1}{5} \sum_{j=0}^p \binom{5}{j} B_j n^{5-j}$$
$$=\frac15(B_0 n^5 + 5B_1 n^4 + 10 B_2 n^3 + 10 B_3 n^2 + 5 B_4 n)$$
$$=\frac15n^5 + \frac12 n^4 + \frac 13 n^3 - \frac {1}{30}n$$
A: The Summa Potestatum involves many other relations.
Some of the more interesting are
Recursion
$$ \bbox[lightyellow] {  
S_m (n) = \sum\limits_{0\, \le \,k\, \le \,n - 1} {k^{\,m} }  = \left[ {1 \le n} \right]\left( {\left[ {0 = m} \right] + \sum\limits_k {\left( \matrix{
  m \cr 
  k \cr}  \right)S_k (n - 1)} } \right)\quad \left| {\;0 \le {\rm integer }m,n} \right.
} \tag{1}$$
where $[P]$ denotes the Iverson bracket
Relation with Eulerian, Stirling and Bernoulli Numbers
$$ \bbox[lightyellow] {  
\eqalign{
  & S_m (n) = \sum\limits_{0\, \le \,k\, \le \,n - 1} {k^{\,m} } \quad \left| {\;0 \le {\rm integer }m,n} \right. =   \cr 
  &  = \sum\limits_{\left( {0\, \le } \right)\,j\,\left( { \le \,m} \right)} {\left\langle \matrix{
  m \hfill \cr 
  j \hfill \cr}  \right\rangle \left( \matrix{
  n + j \cr 
  m + 1 \cr}  \right)}  = \sum\limits_{\left( {0\, \le } \right)\,j\,\left( { \le \,m} \right)} {\;j!\;\left\{ \matrix{
  m \cr 
  j \cr}  \right\}\left( \matrix{
  n \cr 
  j + 1 \cr}  \right)}  =   \cr 
  &  = {1 \over {m + 1}}\sum\limits_{0\, \le \,j\, \le \,m} {\left( \matrix{
  m + 1 \cr 
  j \cr}  \right)\;B_{\,j} \;n^{\,m + 1 - j} }  \cr} 
}\tag{2}$$
where 
 - the angle braces denote the Eulerian Numbers (1st kind)
 - the curly braces correspond to the Stirling Numbers of 2nd kind
 - $B_j$ are the Bernoulli Numbers ($B_1=-1/2$).
The first two, at least, use numbers well defined in its own and so avoid the "circularity" you are lamenting.
We can see how for instance it comes that
$$
\eqalign{
  & S_3 (n) = \sum\limits_{\left( {0\, \le } \right)\,j\,\left( { \le \,m} \right)} {\;j!\;\left\{ \matrix{
  3 \cr 
  j \cr}  \right\}\left( \matrix{
  n \cr 
  j + 1 \cr}  \right)}  = \left( {\left( \matrix{
  n \cr 
  2 \cr}  \right) + 6\left( \matrix{
  n \cr 
  3 \cr}  \right) + 6\;\left( \matrix{
  n \cr 
  4 \cr}  \right)} \right) =   \cr 
  &  = {{n\left( {n - 1} \right)} \over 2} + n\left( {n - 1} \right)\left( {n - 2} \right) + {{n\left( {n - 1} \right)\left( {n - 2} \right)\left( {n - 3} \right)} \over 4} =   \cr 
  &  = {{n\left( {n - 1} \right)\left( {n^{\,2}  - n} \right)} \over 4} = \left( {{{n\left( {n - 1} \right)} \over 2}} \right)^{\,2}  = S_1 (n)^{\,2}  \cr} 
$$
The Faulhaber's Formula that Achille hui already indicated, provides much more insight.
Relation with Bernoulli Polynomial and Hurwitz Zeta
Finally it is interesting to note that the Indefinite Sum is
$$ \bbox[lightyellow] {
\eqalign{
  & \sum\nolimits_{\;x\;} {x^{\,m} }  + c\quad \left| {\;0 \le {\rm integer }m} \right. =   \cr 
  &  = \left( {\sum\limits_{\left( {0\, \le } \right)\,j\,\left( { \le \,m} \right)} {j!\left\{ \matrix{
  m \cr 
  j \cr}  \right\}\left( \matrix{
  x \cr 
  j + 1 \cr}  \right)} } \right) =   \cr 
  &  = \sum\limits_{\left( {0\, \le } \right)\,j\,\left( { \le \,m} \right)} {\left\langle \matrix{
  m \cr 
  j \cr}  \right\rangle \left( \matrix{
  x + j \cr 
  m + 1 \cr}  \right)}  = {1 \over {m + 1}}\sum\limits_{\left( {0\, \le } \right)\,j\,\left( { \le \,m + 1} \right)} {\left( \matrix{
  m + 1 \cr 
  j \cr}  \right)\;B_j \;x^{\,m + 1 - j} }  =   \cr 
  &  = {1 \over {m + 1}}B(m + 1,x) =  - \,\zeta ( - m,x) \cr} 
} \tag{3}$$
where $B(n,x)$ are the Bernoully polynomials and $\zeta(n,x)$ the Hurwitz Zeta Function.
And it is also possible to establish a connection with the Generalized Harmonic numbers (as you already found) and many others.
