Simplifying a rational function in terms of Bernoulli numbers and polynomials

Question

Faulhaber's formula expresses a sum over some finite number of naturals to the $m^{th}$ power in terms of the Bernoulli numbers $B_{j}$ (using the $B_{1} = 1/2$ convention) or polynomials $\hat{B}_{j}$ as follows

$$\sum_{k=1}^{n} k^{m} = \frac{1}{m+1} \sum_{j=0}^{m} \begin{pmatrix} m+1 \\ j \end{pmatrix} B_{j} n^{m+1-j} = \frac{1}{m+1} \left( \hat{B}_{m+1}(n+1) - \hat{B}_{m+1}(1) \right) \tag 1$$

I recently needed to derive an expression similar to $(1)$ but for a double summation over squares $k_{1}^{2} + k_{2}^{2}$ rather than $k$, summing to arbitrary $n_{1}, n_{2}$ and with $m \in \mathbb{N} \setminus \{1\}$. To this end, using the binomial theorem and then Faulhaber's formula on the resultant product gives

\begin{align} \sum_{k_{1},k_{2}=1}^{n_{1},n_{2}} (k_{1}^{2} + k_{2}^{2})^{m} &= \sum_{k_{1},k_{2}=1}^{n_{1},n_{2}} \sum_{i=0}^{m} \begin{pmatrix} m \\ i \end{pmatrix} k_{1}^{2i} k_{2}^{2m-2i} \\ &= \sum_{i=0}^{m} \begin{pmatrix} m \\ i \end{pmatrix} \sum_{k_{1}=1}^{n_{1}} k_{1}^{2i} \sum_{k_{2}=1}^{n_{2}} k_{2}^{2m-2i} \\ &= \sum_{i=0}^{m} \begin{pmatrix} m \\ i \end{pmatrix} \frac{1}{2i+1} \frac{1}{2m-2i+1} \left[\sum_{j_{1}=0}^{2i} \begin{pmatrix} 2i+1 \\ j_{1} \end{pmatrix} B_{j_{1}} n_{1}^{2i+1-j_{1}} \sum_{j_{2}=0}^{2m-2i} \begin{pmatrix} 2m-2i+1 \\ j_{2} \end{pmatrix} B_{j_{2}} n_{2}^{2m-2i+1-j_{2}} \right] \\ &= \sum_{i=0}^{m} \begin{pmatrix} m \\ i \end{pmatrix} \frac{1}{2i+1} \frac{1}{2m-2i+1} \bigg[\left( \hat{B}_{2i+1}(n_{1}+1) - \hat{B}_{2i+1}(1) \right) \left( \hat{B}_{2m-2i+1}(n_{2}+1) - \hat{B}_{2m-2i+1}(1) \right) \bigg] \end{align}

What I then wanted to do was to use this to find a polynomial $G(l, m, n_{1}, n_{2})$ such that

$$G(l, m, n_{1}, n_{2}) = \frac{\sum_{k_{1},k_{2}=1}^{n_{1},n_{2}} (k_{1}^{2} + k_{2}^{2})^{m+l}}{\sum_{k_{1},k_{2}=1}^{n_{1},n_{2}} (k_{1}^{2} + k_{2}^{2})^{m}} \tag 2$$

for arbitrary values of $(l, m, n_{1}, n_{2})$ with $l \in \mathbb{N}$. I was hoping to be able to find $G$ in terms of the Bernoulli polynomials and numbers themselves, though unfortunately I haven't been able to make any headway.

I was wondering if it is possible to construct a $G$ that satisfies $(2)$ for arbitrary naturals $(l, m, n_{1}, n_{2})$? I don't mind a construction for the special case $n_{1} = n_{2} = n$ either. Or does anyone know any references they could point me to? Any help is appreciated.

Answer 1

For the special case $n_1 = n_2 = n$ and $1 \le m \le 9$ it turns out that the double sum yields a polynomial of the form $P(n)(2n+1)(n+1)n^2$ where $P(n)$ is irreducible of degree $2(m-1)$. Obviously the ratio of two irreducible polynomials of different degrees cannot be a polynomial. It may be that for higher values of $m$ the polynomial factors further, but I wouldn't be optimistic about finding polynomial ratios.

---

My Sage code (online demo):

n = var('n')


def extrapolate_polynomial(ys, x):
    """
    Extrapolates a polynomial using Neville's algorithm.

    ys: the values of the polynomial at 0, 1, 2, ...
    x: the value at which to evaluate the polynomial
    """
    n = len(ys)
    for delta in range(1, n):
        ys = [(x - i) * ys[i+1] - (x - i - delta) * ys[i] for i in range(n - delta)]
    return ys[0] / factorial(n - 1)


def faulhaber(k, n):
    ys = [0] * (k + 2)
    for x in range(1, k + 2):
        ys[x] = ys[x-1] + x**k

    return extrapolate_polynomial(ys, n)


def double_sum(m):
    return expand(sum(binomial(m,i) * faulhaber(2*i, n) * faulhaber(2*m-2*i,n) for i in range(m+1)))


print(factor(double_sum(3) / double_sum(2)))
print("")
for m in range(1, 10):
    print(m, "\t", factor(double_sum(m)))