Is $\sum_{k=1}^{n} \sin(k^2)$ bounded by a constant $M$?

I know $\sum_{k=1}^{n} \sin(k)$ is bounded by a constant. How about $\sum_{k=1}^{n} \sin(k^2)$?

We have: $$\begin{eqnarray*}2\left(\sum_{k=1}^{n}\sin(k^2)\right)^2 &=& \sum_{j,k=1}^{n}\cos(j^2-k^2)-\sum_{j,k=1}^{n}\cos(j^2+k^2)\\&=&n+2\sum_{m=1}^{n^2-1}d_1(m)\cos(m)-2\sum_{m=2}^{2n^2}d_2(m)\cos m\end{eqnarray*}$$ where $d_1(m)$ accounts for the number of ways to write $m$ as $j^2-k^2$ with $1\leq k<j\leq n$ and $d_2(m)$ accounts for the number of ways to write $m$ as $j^2+k^2$ with $1\leq j,k\leq n$. Since both these arithmetic functions do not deviate much from their average order (by Dirichlet's hyperbola method $d_1(m)$ behaves on average like $\log m$ and $d_2(m)$ behaves on average like $\frac{\pi}{4}$), I think it is not difficult to prove that for infinitely many $n$s $$ \left|\sum_{k=1}^{n}\sin(k^2)\right|\geq C\sqrt{n} $$ holds for some absolute constant $C\approx \frac{1}{\sqrt{2}}$ through the Cauchy-Schwarz inequality.

However, this is quite delicate: Weyl bounds work just in the opposite direction. As an alternative approach, we may consider only the set of $n$ that appear in the numerators of the convergents of $2\pi$, in order that our sum behaves like a Gaussian sum (with magnitude $\sqrt{n}$ or $\sqrt{2n}$) plus a small error. In both cases we get that the sequence given by $$S_n=\sum_{k=1}^{n}\sin(k^2)$$ is not bounded.

In fact, Terence Tao already proved that on MO.