TL;DR: "I've made something work but I do not know why it works - please explain it to me (or disprove my finding)"
DISCLAIMER: A lot of things here are done in a non-rigorous way - which is all the more surprising that in the end it "all works".
Let's define $f_m(n) = \sum_{k=0}^nk^m$ . For example, $f_1(n) = \sum_{k=0}^nk$ which is just a triangular numbers summation. What I wanted to find is the "closed" form for $f_m(n)$ meaning a polynomial form. There are two things I used here as a "mathematical intuition":
- The resulting polynomial must be of m+1 degree. That's because we're adding the m-th powers of n exactly n times
- The differentiation "works" for discrete sums (with some "adjustments" - probably the crux of this question)
Assuming these two points I came out with an idea that - well, $f'_{m+1}(n) = (m+1)\sum_{k=0}^nk^m = (m+1)f_m(n)$ so we can "integrate" $f_m(n)$ to get the expression for $f_{m+1}(n)$ . However if used directly, it will not yield correct results. For instance, if we use it on $f_0(n)$ and $f_1(n)$ we run into a problem: $\int{f_0(n)} = \int{n} = \frac{n^2}{2}$ (I intentionally don't write "dn" here as it makes no sense). As we know the triangular sum is $\frac{n(n+1)}{2} = \frac{n^2}{2} + \frac{n}{2}$, i.e. we are missing a term here.
But wait - when we're integrating, we have to also add a "constant term" - which in this case I will denote $L(n)$. The $L(n)$ is a linear function of $n$, i.e. $L(n) = An+B$ - furthermore, since for all m we have $f_m(0) = 0$ we get that $B = 0$ so $L(n) = An$ for some number $A$. To get the number $A$ we just substitute the $n=1$ to a given $f_m(n)$ and work out what it should be. For example, $f_1(n) = \frac{n^2}{2} + L(n)$, so $L(1) = A\cdot 1 = f_1(1) - \frac{1^2}{2}$ and since $f_1(1) = 1$ we get that $A = \frac{1}{2}$ thus a "correct" expression is: $f_1(n) = \frac{n^2}{2} + L(n) = \frac{n^2}{2} + \frac{n}{2}$ which is exactly what we would expect. How convenient!
Now if all of this sounds like borderline nonsense to you - that's okay - because I also think it is. But the interesting thing is - it works for any $f_m(n)$ (at least I think it does, I've checked it for several higher values). For example:
$$f_1(n) = \frac{n^2}{2} + \mathbf{\frac{n}{2}} \\f_2(n) = \frac{n^3}{3} + \frac{n^2}{2} + \mathbf{\frac{n}{6}} \\f_3(n) = \frac{n^4}{4} + \frac{n^3}{2} + \frac{n^2}{4} + \mathbf{0} \\f_4(n) = \frac{n^5}{5} + \frac{n^4}{2} + \frac{n^3}{3} \mathbf{- \frac{n}{30}} \\f_5(n) = \frac{n^6}{6} + \frac{n^5}{2} + \frac{5n^4}{12} - \frac{n^2}{12} + \mathbf{0}$$
So, in general, $f_{m+1}(n) = (m+1)\int{f_m(n)} + L(n)$ and of course, in cases when we find that $A=0$ we obviously get $L(n)=0$
So - why does it work? The method that I used is .. dubious to say the least. There was also no indication that $L(n)$ is a linear function - because for instance, for cases like $f_4(n)$ we have our sum of known members terminating at a 3rd power (i.e. why can't $L(n)$ be a quadratic?)
Or maybe it doesn't work at all and the expressions I found are just a lucky coincidence?