Q:

Patterns of the zeros of the Faulhaber polynomials (modified)

Faulhaber polynomial of order $p \in \Bbb{N}$ is defined as the unique polynomial of degree $p+1$ satisfying
$$ S_{p}(n) = \sum_{k=1}^{n} k^p $$
for $n = 1, 2, 3, \cdots$. For example,
\begin{align*}
S_0(x) &= x, \\
S_1(x) &= \frac{x(x+1)}{2}, \\
S_2(x) &= \frac{x(x+1)(2x+1)}{6}, \\
S_3(x) &= \frac{x^2 (x+1)^2}{4}.
\end{align*}
In order to grasp some intuition on the partial decomposition of $1/S_p (x)$, I tried plotting the complex zeros of $S_p (x)$. The following graphics shows the distribution of the zeros of $S_{800}(x)$.

(The precision of the calculated zeros $z_j$ of $S_{800}(z)$ above satisfy $|f(z_j)| \leq 10^{-300}$.)
It turns out that they exhibits a very neat, yet still a strange pattern as seen above.
So far I have never heard of the topic related to this pattern, and I want to know (out of curiosity) if there are some results concerning the pattern of zeros of $S_p(x)$.
Addendum. The previous fuzzy graphics turned out to be a result of numerical error. Now I removed those errors.

A:

(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.

A:

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

