Patterns of the zeros of the Faulhaber polynomials (modified)

(This is just a comment, really.)

First it appears that the zeros are symmetric about the line $x=-1/2$, and indeed the polynomials

$$ F_p(z) = S_p(z-1/2) $$

appear to have only even or only odd powers of $z$.

It seems that the zeros not on the real axis grow on the order of $p/(2\pi)$. Numerically they have the same limiting behavior as the zeros of the partial sums of the sine and cosine series (see this paper [PDF]) as well as the partial sums of the Bessel functions (see this preprint).

Below is a plot of the zeros of $F_p\left(\frac{p}{2\pi}z\right)$ for $p=400$, along with the modified Szegő curve

$$ \begin{align} &\left\{z \in \mathbb{C} \,\colon \Im(z) \geq 0,\,\,\, |z| \leq 1, \,\,\,\text{and}\,\,\, \left|ze^{1+iz}\right| = 1 \right\} \\ &\qquad \cup \,\left\{z \in \mathbb{C} \,\colon \Im(z) \leq 0,\,\,\, |z| \leq 1, \,\,\,\text{and}\,\,\, \left|ze^{1-iz}\right| = 1 \right\} \\ &\qquad \cup \,\left\{x \in \mathbb{R} \,\colon -1/e \leq x \leq 1/e \right\} \end{align} $$

in blue.

There is possibly a connection to the zeros of the Bernoulli polynomials $B_p(x)$ as a result of the fact that

$$ S_p(z) = \frac{B_{p+1}(z+1) - B_{p+1}(0)}{p+1}. $$

You may wish to take a look at Karl Dilcher's memoir Zeros of Bernoulli, Generalized Bernoulli, and Euler Polynomials and this paper by John Mangual.

In addition to Antonio's comment/answer: Looking at the real roots (symmetrized by adding +1/2 (!)) of the 1,5,9,13,... polynomial we get the following list, where only the first three real roots are rational numbers. The rate of convergence to the half-integers is impressive... $$\small \begin{matrix} 0 & 1/2 & . & . & . & . & . & . \\ 0 & 1/2 & -1/2 & 0.763762615826 & -0.763762615826 & . & . & . \\ 0 & 1/2 & -1/2 & 0.949106003964 & -0.949106003964 & . & . & . \\ 0 & 1/2 & -1/2 & 0.999056597832 & -0.999056597832 & . & . & . \\ 0 & 1/2 & -1/2 & 0.999997848581 & -0.999997848581 & . & . & . \\ 0 & 1/2 & -1/2 & 0.999999998198 & -0.999999998198 & -1.50196566814 & 1.50196566814 & 1.74815179290 \\ 0 & 1/2 & -1/2 & 0.999999999999 & -0.999999999999 & -1.50001155318 & 1.50001155318 & 1.93305092402 \\ 0 & 1/2 & -1/2 & 1.00000000000 & -1.00000000000 & -1.50000003663 & 1.50000003663 & 1.99704558735 \\ 0 & 1/2 & -1/2 & 1.00000000000 & -1.00000000000 & -1.50000000007 & 1.50000000007 & 1.99997147602 \\ 0 & 1/2 & -1/2 & 1.00000000000 & -1.00000000000 & -1.50000000000 & 1.50000000000 & 1.99999984071 \\ 0 & 1/2 & -1/2 & 1.00000000000 & -1.00000000000 & -1.50000000000 & 1.50000000000 & 1.99999999943 \\ 0 & 1/2 & -1/2 & 1.00000000000 & -1.00000000000 & -1.50000000000 & 1.50000000000 & 2.00000000000 \\ 0 & 1/2 & -1/2 & 1.00000000000 & -1.00000000000 & -1.50000000000 & 1.50000000000 & 2.00000000000 \\ 0 & 1/2 & -1/2 & 1.00000000000 & -1.00000000000 & -1.50000000000 & 1.50000000000 & 2.00000000000 \\ 0 & 1/2 & -1/2 & 1.00000000000 & -1.00000000000 & -1.50000000000 & 1.50000000000 & 2.00000000000 \\ 0 & 1/2 & -1/2 & 1.00000000000 & -1.00000000000 & -1.50000000000 & 1.50000000000 & 2.00000000000 \end{matrix} $$ The complex roots may be fit into that pattern perhaps by their absolute values, but this naive idea is not yet convincing to me

As an addendum to Antonio Vargas's answer, let's prove that the roots of $S_p$ are indeed symmetrically distributed around $-1/2$, or in other words that if $r$ is a root, then so is $-1-r$. A somewhat more precise result is that $S_p(-1-x) = S_p(x)$ when $p$ is odd, and $S_p(-1-x) = -S_p(x)$ when $p$ is even (thus confirming that the polynomial $T(x) := S_p(x-1/2)$ is an odd function (has only odd powers of $x$) if $p$ is even, and is an even function (has only even powers of $x$) if $p$ is odd).

The key principle is that a polynomial is the zero polynomial iff it has infinitely many roots. For example, $f(x) = g(x)$ if $f(n) = g(n)$ for all $n \in \mathbb{N}$, and a polynomial $h$ is constant if $h(n) = h(0)$ for all $n \in \mathbb{N}$. If $h$ is a polynomial such that its first difference $(\Delta h)(x) := h(x) - h(x-1)$ satisfies $(\Delta h)(n) = 0$ for all $n \in \mathbb{N}$, then $h$ is constant with constant value $h(0)$.

Thus we have $(\Delta S_p)(x) = S_p(x) - S_p(x-1) = x^p$ and so $S_p(-x) - S_p(-1-x) = (-1)^p x^p$. For fixed $p$, put $g(x) := S_p(-1-x)$, so that $g(x-1) = S_p(-1-(x-1)) = S_p(-x)$. Thus

$$\Delta g(x) = g(x) - g(x-1) = S_p(-1-x) - S_p(-x) = (-1)^{p+1} x^p$$

which equals $\Delta S_p(x)$ if $p$ is odd, and $-\Delta S_p(x)$ if $p$ is even. Thus in the case where $p$ is odd, $\Delta(g-S_p) = 0$; applying the principle above, $g$ and $S_p$ differ by a constant: $S_p(-1-x) = S_p(x) + c$ for some constant $c$. For $x=0$, we note that $S_p(0) = 0$, and $S_p(0) - S_p(-1) = 0^p = 0$ so $S_p(-1) = 0$. It follows that $c = 0$, and we conclude $S_p(-1-x) = S_p(x)$ if $p$ is odd.

The case where $p$ is even is wholly similar. There we derive $\Delta(g + S_p) = 0$, so $g(x) + S_p(x) = c$ for some constant $c$, and we argue as before that $c = 0$, and so in this case $S_p(-1-x) = -S_p(x)$.