A very hard proof that $\sum_{i=0}^n i = \frac{n(n+1)}2$ (in comparison with the elementary level of the identity) is to compare degrees in the Vandermonde identity, which you prove playing around with the columns and use induction. Using Laplace expansion, the determinant of the Vandermonde matrix is a polynomial of degree $0 + 1+2+ \cdots + n$ and the result is $\prod_{0 \le i < j \le n} (x_j - x_j)$, a polynomial of degree $\binom{n+1}2$.
Knowing that there are formulas for the sum of $k^{\text{th}}$ powers, this begs the following question : is there a neat formula for the "$k$-Vandermonde" polynomials $$ \det \left( \begin{bmatrix} 1 & x_0^{1^k} & x_0^{2^k} & \cdots & x_0^{n^k} \\ 1 & \vdots & \vdots & \ddots & \vdots \\ 1 & \vdots & \vdots & \ddots & \vdots \\ 1 & x_n^{1^k} & x_n^{2^k} & \cdots & x_n^{n^k} \end{bmatrix} \right). $$ I conjecture in the case $k=2$ that $$ \det \left( \begin{bmatrix} 1 & x_0^{1^2} & x_0^{2^2} & \cdots & x_0^{n^2} \\ 1 & \vdots & \vdots & \ddots & \vdots \\ 1 & \vdots & \vdots & \ddots & \vdots \\ 1 & x_n^{1^2} & x_n^{2^2} & \cdots & x_n^{n^2} \end{bmatrix} \right) = \left( \prod_{0 \le i < j \le n} (x_j - x_i) \right) \left( \sum_{0 \le i \le j \le n} x_i x_j \right)( \cdots ) $$ since the formula holds for $n=2$ (checked with Wolfram Alpha, in which case the $(...)$ factor equals $1$ ; the conjecture just suggests that these polynomials divide the rest). I know that the first factor divides the determinant ; this is just because the determinant vanishes when $x_i = x_j$, so in the UFD $\mathbb Z[x_0,\cdots,x_n]$, the distinct irreducibles $x_i - x_j$ divide the determinant. Also note that the $k$-Vandermonde divided by the $1$-Vandermonde is therefore a symmetric polynomial.
EDIT : I checked $n=3$ as the comments suggested, and indeed my conjecture is wrong, but the claims I made during the conjecture (about the product of the linear factors) is true.
Any ideas?
If you divide both sides of your conjecture by the Vandermonde, then the l.h.s. is exactly the Schur polynomial $s_{\lambda}(x_1,\cdots,x_n)$ where $\lambda_i+n-i=(n-i)^k$. Thus the $(\cdots)$ part is quite complicated if you now refer to the fact that $s_\lambda=\sum_T x^T$, the sum over semistandard young tableau of shape $\lambda$. To see that your conjecture doesn't quite work out, pick something like $n=4$ and find a semi-standard young tableau of shape $\lambda$ which has more than two distinct numbers.