Finite convolution sum of power function

Question

I am interested in understanding the behavior of the convolution sum $$\sum_{k = 1}^{n} k^r (n - k)^s$$ for integer $n$, $r$, and $s$. I imagine that there must be some larger theory surrounding sums of this form, but the only idea that has occurred to me to simplify this expression is to apply the binomial theorem to the $(n-k)^s$ term, to get

$$\sum_{k = 1}^n k^r \sum_{i = 0}^{s} {s \choose i}n^{s - i} (-k)^i$$ We can rewrite this expression as

$$\sum_{i = 0}^s (-1)^i {s \choose i }n^{s - i} \sum_{k = 1}^{n} k^{r + i}$$ I noted that the inner-most sum is a Faulhaber sum (https://en.wikipedia.org/wiki/Faulhaber%27s_formula), but applying Faulhaber's formula doesn't seem to clarify the expression at all. I just thought I'd ask if there is some obvious technique I am missing here, or if there are any resources that explicitly deal with this expression?

Thank you!

Answer 1

Not very useful but, for a given value of $s$, the result is a combination of powers of $n$, and combinations of Riemann and Hurwitz $\zeta$ functions with negative arguments.

For example, for $s=3$, the sum will be $$\color{red}{1}\Big[\zeta (-r)-\zeta (-r,(n+1))\Big]n^3-\color{red}{3}\Big[\zeta (-(r+1))-\zeta (-(r+1),(n+1))\Big]n^2+$$ $$\color{red}{3}\Big[\zeta (-(r+2))-\zeta (-(r+2),(n+1))\Big]n-\color{red}{1}\Big[\zeta (-(r+3))-\zeta (-(r+3),(n+1))\Big]$$ and we can easily conjecture the general form since the "red" coefficients are those of the expansion of $(n-1)^s$ and the pattern of the remaining terms is quite obvious.

Answer 2

To get a simpler formula we replace $n$ with $n-1$ and start the sum with $k=0$ (with $0^0=1$), defining $$S(n,r,s) = \sum_{k=0}^{n-1}k^r(n-1-k)^s.$$ Then $$ \sum_{r,s=0}^\infty S(n,r,s) \frac{x^r}{r!}\frac{y^s}{s!} = \frac{e^{nx}-e^{ny}}{e^x-e^y}. $$ The coefficients can be expressed in terms of Bernoulli numbers, since $$ \begin{aligned} \frac{e^{nx}-e^{ny}}{e^x-e^y}&= \frac{e^{nx}-e^{ny}}{x-y}\frac{x-y}{e^x-e^y}\\ &=\frac{e^{nx}-e^{ny}}{x-y} \cdot e^{-y}B(x-y), \end{aligned} $$ where $B(z)=z/(e^z-1)$ is the exponential generating function for the Bernoulli numbers, though it's not clear if this formula will really be helpful.