When is this sum of perfect powers bounded

For any positive integers $n,d$, let

$$ A_d(n)=\frac{\sum_{k=1}^n k^{2d}}{n(n+1)(2n+1)} $$

It is easy to see (and well-known) that for fixed $d$, $A_d(.)$ is a polynomial of degree $2d-2$. Writing $A_d(x)$ then makes sense for any $x\in{\mathbb R}$, not just the positive integers.

Is it known for which real numbers $x$ the sequence $(A_d(x))_{d\geq 1}$ is bounded ?

UPDATE : this question was also asked on MO

Markus Scheuer pointed to the literature which contains everything needed for a much
more comprehensive and clean answer. Still let me leave this one as it was, it deserves a separate answer.

... Hope it is not sheer bounty hunting on my side but here is some evidence that $A_d(x)$ as a whole might have nice limit as $d\to\infty$: real roots of $A_d(x)$ very quickly become very close to half-integers. For example, here are the real roots of $A_{31}(x)$:

$$ \begin{array}{rl} -4.&7392448816946959003622314391353847421996102737947803419494492900007\\ -4.&5027652294851661968345082151000807931898466324185824199087980186141\\ -3.&9999998139020492031175531389561140529018995821878583363385222409911\\ -3.&5000000000022935383500658889864185834363214538575577529515462763107\\ -2.&9999999999999999977505041317882254571396873898232742129535442958828\\ -2.&5000000000000000000000000403537157103973597828990762521212139197237\\ -1.&9999999999999999999999999999999999995122183385375255793622791417059\\ -1.&5000000000000000000000000000000000000000000000000000001057707874112\\ 0.&5000000000000000000000000000000000000000000000000000001057707874112\\ 0.&9999999999999999999999999999999999995122183385375255793622791417059\\ 1.&5000000000000000000000000403537157103973597828990762521212139197237\\ 1.&9999999999999999977505041317882254571396873898232742129535442958828\\ 2.&5000000000022935383500658889864185834363214538575577529515462763107\\ 2.&9999998139020492031175531389561140529018995821878583363385222409911\\ 3.&5027652294851661968345082151000807931898466324185824199087980186141\\ 3.&7392448816946959003622314391353847421996102737947803419494492900007 \end{array} $$ and the pattern continues on

This indicates that after some normailization $A_d(x)$ might well tend to something like $\sin(2\pi x)$.

Here are the superimposed plots of $\frac{32x(x+1)(2x+1)A_d(x)}{-3A_d(-\frac14)}$ for $d$ up to 28

Later addition - let me reproduce here the stuff from an answer on MO:

for $S_{2d}(x):=x(x+1)(2x+1)A_d(x)$ I can prove $$ \lim_{d\to\infty}\frac{4^{2d+1}}{E_{2d}}S_{2d}(x)=\sin(2\pi x) $$ where $E_{2d}$ are the Euler numbers, modulo another limit expression, now purely numerical, involving Euler and Bernoulli numbers (see below).

Before starting the proof let me explain why this particular expression. As I said above, there was some evidence that $S_{2d}(x)$ behaves as $d$ grows more and more like $c_d\sin(2\pi x)$ where $c_d$ are some (rapidly growing) constants. Now if we hope that $\lim_{d\to\infty}S_{2d}(x)/c_d$ is indeed $\sin(2\pi x)$, then we can fix any $x$ where $\sin(2\pi x)$ is nonzero and where we can calculate $S_{2d}(x)$, and take $c_d=S_{2d}(x)/\sin(2\pi x)$. And $x=-1/4$ seems to be reasonable choice, see below.

I will use Faulhaber's formula $$ \sum_{k=1}^nk^d=\frac1{d+1}\sum_{j=0}^d(-1)^j\binom{d+1}jB_jn^{d+1-j} $$ where $B_j$ are the Bernoulli numbers. On one hand it explicitly gives the coefficients of the polynomial $x(x+1)(2x+1)A_d(x)=S_{2d}(x)$: we can rewrite Faulhaber's formula as $$ S_{2d}(x)=\sum_{i=0}^{d-1}\binom{2d}{2i+1}\frac{B_{2(d-i)}}{2(d-i)}x^{2i+1}+\frac{x^{2d}}2+\frac{x^{2d+1}}{2d+1}. $$

On the other hand, changing slightly the proof of Faulhaber's formula as given in Wikipedia (or maybe using some other argument) one can compute $S_{2d}(-1/4)$. First for general $n$, one considers the exponential generating function: $$ \sum_{d=0}^\infty\frac{S_{2d}(n)}{(2d)!}z^{2d} =\sum_{d=0}^\infty\frac1{(2d)!}\sum_{k=1}^nk^{2d}z^{2d} =\sum_{k=1}^n\sum_{d=0}^\infty\frac{(kz)^{2d}}{(2d)!}=\sum_{k=1}^n\cosh(kz); $$
the latter can be further collected as $$ \sum_{k=1}^n\cosh(kz) =\sum_{k=1}^n\frac{e^{kz}+e^{-kz}}2 =\frac12\left(\sum_{k=1}^n(e^z)^k+\sum_{k=1}^n(e^{-z})^k\right) =\frac12\left(e^z\frac{1-e^{nz}}{1-e^z}+e^{-z}\frac{1-e^{-nz}}{1-e^{-z}}\right) =\frac12\left(\cosh(nz)-1+\coth\left(\frac z2\right)\sinh(nz)\right). $$ Substituting $n=-\frac14$ gives $$ \frac12\left(\cosh\left(\frac z4\right)-1-\coth\left(\frac z2\right)\sinh\left(\frac z4\right)\right) =\frac12\operatorname{sech}\left(\frac z4\right)-1 =-\frac12+\sum_{d=1}^\infty\frac{E_{2d}}{2\times4^{2d}(2d)!}z^{2d} $$ Thus $S_{2d}(-1/4)=\frac{E_{2d}}{2\times4^{2d}}$.

Putting these two calculations together we see that we will have $$ \lim_{d\to\infty}\frac{S_{2d}(x)}{-S_{2d}(-1/4)}=\sin(2\pi x) $$ if we can prove $$ \lim_{d\to\infty}\frac{\binom{2d}{2i+1}\frac{B_{2(d-i)}}{2(d-i)}}{-\frac{E_{2d}}{2\times4^{2d}}}=(-1)^i\frac{(2\pi)^{2i+1}}{(2i+1)!} $$ and from here on I don't know how to proceed, except for seemingly slightly more manageable modification $$ 2\lim_{d\to\infty}4^{2d-i}\frac{B_{2(d-i)}}{(2(d-i))!}/\frac{E_{2d}}{(2d)!}=(-1)^{i+1}\pi^{2i+1} $$

I've been told on MO that this follows from known asymptotics of Bernoulli and Euler numbers but I don't know how to estimate precise errors of their asymptotic approximations and without that I do not know how to go about computing the required limit of their ratio.

In a comment above Markus Scheuer gave links to materials which readily contain an answer. I will use highly informative lecture notes by Omran Kouba pointed out by him.

It follows from the formulæ on page 6 of these notes that what we are looking at are actually the Bernoulli polynomials: $$ S_{2d}(x):=x(x+1)(2x+1)A_d(x)=\frac1{2d+1}B_{2d+1}(x+1) $$ which by Corollary 5.2 (page 18) of the same notes gives that $(-1)^{d+1}\frac{(2\pi)^{2d+1}}{2(2d)!}S_{2d}(x)$ converges to $\sin(2\pi x)$ uniformly on every compact subset of $\mathbb C$.

It thus follows that for large $d$, $A_d(x)$ behaves like $(-1)^{d+1}\frac{2(2d)!}{(2\pi)^{2d+1}x(x+1)(2x+1)}\sin(2\pi x)$, i. e. grows quite rapidly - by Stirling's formula, $\log(A_d(x))$ grows like $2d\log\frac d{\pi e}+$ const.

Well in fact as Ewan Delanoy points out in the comment below this only shows unboundedness of $A_d(x)$ when $x$ is not a half-integer. To capture this case too, first note that $$ A_d(0)=\lim_{x\to0}\frac{B_{2d+1}(x+1)}{(2d+1)x(x+1)(2x+1)} =\lim_{x\to0}\frac{B_{2d+1}'(x+1)}{(2d+1)(6x^2+6x+1)}=B_{2d}(0) $$ since $B_n'(t)=nB_{n-1}(t)$ and $B_{2d}(1)=B_{2d}(0)$, and similarly $$ A_d(-\frac12)=-2B_{2d}(\frac12)=2(1-\frac1{2^{2d-1}})B_{2d}(0). $$ For other half-integers use $B_n(x+1)=B_n(x)+nx^{n-1}$ (all the needed equalities are in the mentioned notes) to obtain $$ A_d(x)=\frac{(x-1)(2x-1)A_d(x-1)+x^{2d-1}}{(x+1)(2x+1)} $$ for $x\ne-1,-1/2$ and $$ A_d(x-1)=\frac{(x+1)(2x+1)A_d(x)-x^{2d-1}}{(x-1)(2x-1)} $$ for $x\ne1,1/2$. This gives $A_d(-1)=A_d(0)$, $$ A_d(1/2)=A_d(-3/2)=\frac1{3\cdot2^{2d-1}}, $$ $$ A_d(1)=A_d(-2)=\frac16 $$ and shows (by induction) that at all other half-integers and integers $A_d$ is unbounded.

Thus $A_d(x)$ goes to infinity with $d$ everywhere except that it goes to zero at $x=1/2,-3/2$ and is constantly $1/6$ at $x=1,-2$.

PS all this could be slightly simplified using $A_d(-1-x)=A_d(x)$...