Distribution of the Error term in GH Hardy's "curious result" $\sum_{\nu \leq n } \{ \nu \theta \}^2 = \tfrac{1}{12} n + O(1)$

There are several puzzling things about the question: Firstly of course $\theta$ must be irrational, and it is intended for $\{ x\}$ to denote the Bernoulli polynomial $x-[x]-1/2$ rather than the more usual fractional part. Secondly, where is the result of Hardy from? I did find this statement in the Cambridge ICM paper of Hardy and Littlewood where they write

"While engaged on the attempt to elucidate these questions we have found a curious result which seems of sufficient interest to be mentioned separately. It is that $$ \sum_{\nu =1}^{n} \{ \nu \theta\}^2 = \frac{n}{12} +O(1) $$ for all irrational values of $\theta$. When we consider the great irregularity and obscurity of $\ldots$, it is not a little surprising that [this] (and presumably the corresponding sums with higher even powers) should behave with such marked regularity."

Note that Hardy and Littlewood also use $\{x\}$ to denote the Bernoulli polynomial. This is puzzling since it seems completely false if for example $\theta$ is a Liouville number. Indeed then I found a follow up paper of Hardy and Littlewood where they note (see page 36 there)

"We may take this opportunity of correcting a misstatement in our communication to the Cambridge congress $\ldots$. It was stated there that $$ \sum_{\nu =1}^{n} \{ \nu \theta\}^2 = \frac{n}{12} +O(1) $$ for every irrational $\theta$. This is untrue; but the equation holds for very general classes of values of $\theta$, and in particular for any $\theta$ whose partial quotients are bounded."

What Hardy and Littlewood had in mind is presumably to write out the Fourier expansion of the Bernoulli polynomial $(x-[x]-1/2)^2 -1/12$ which is $$ \frac{1}{2\pi^2} \sum_{k\neq 0} \frac{e^{2\pi i kx}}{k^2}, $$ and then sum this over $x = \nu \theta$. From here one can see that what is needed for their result is that if $\Vert k\theta \Vert$ isn't smaller than $k^{-2+\epsilon}$ then the series will converge nicely, but there will be problems for very well approximable numbers. Using the Fourier expansion, for good irrationals (e.g. $\sqrt{7}$) the Fourier expansion shows that the remainder term will have an almost periodic structure.

As indicated in Lucia's answer, we can look at the Fourier transform

$$\{ \nu \theta \} ^2 = \frac{1}{12} + \frac{1}{2 \pi^2} \sum_{k \neq 0} \frac{\exp(2 \pi i k \nu \theta)}{k^2} $$

Interchanging the sums gives (barring errors on my part)

$$ \sum_{\nu=1}^n \{ \nu \theta \}^2 = \frac{n}{12} + \frac{1}{\pi^2} \sum_{k =1}^{\infty} \frac{1}{k^2} \csc(\pi k \theta) \sin(\pi k n \theta) \cos(\pi k (n+1) \theta) $$

Now, forget that $n$ is a discrete variable, and instead varies continuously. The error term is just a linear combination of offset sine waves, each of which has a period exactly dividing $1 / \theta$.

Consequently, the error term is periodic, and the observed errors should extremely closely resemble the distribution we get from assuming $n$ is a uniformly distributed random variable across one period.

After adding by $0.25$ (since you seem to be taking the sum starting from $\nu=0$), an approximate plot of the first few terms for $\theta = \sqrt{7}$ is:

(obtained from wolframalpha)

From which we should expect the errors to have especially high density near $0.18$, and high density near $0.16$, $0.20$, and near the extremes $0.10$ and $0.26$, and that the errors should have relatively low density in the in-between regions.

This is pretty much exactly what you see from your plot, except for an overall shift, which is probably explained from adding up the average of the tail.

With more precision and taking more terms of the series, we should be able to observe all features of the data.

This is an approach on the error term. The idea is to use Koksma's inequality and Erdos-Turan inequality. I used this on a problem in MSE: $\sum \frac{(-1)^n|\sin n|}n$. To recall,

Theorem [Koksma]

Let $f$ be a function on $I=[0,1]$ of bounded variation $V(f)$, and suppose we are given $N$ points $x_1, \ldots , x_N$ in $I$ with discrepancy $$ D_N:=\sup_{0\leq a\leq b\leq 1} \left|\frac1N \#\{1\leq n\leq N: x_n \in (a,b) \} -(b-a)\right|. $$ Then $$ \left|\frac1N \sum_{n\leq N} f(x_n) - \int_I f(x)dx \right|\leq V(f)D_N. $$

Theorem [Erdos-Turan]

Let $x_1, \ldots, x_N$ be $N$ points in $I=[0,1]$. Then there is an absolute constant $C>0$ such that for any positive integer $m$, $$ D_N\leq C \left( \frac1m+ \sum_{h=1}^m \frac1h \left| \frac1N\sum_{n=1}^N e^{2\pi i h x_n}\right|\right). $$

We are using these on the sequence $x_n = n \theta \ \mathrm{mod} \ 1$ and $f(x) = (x -\frac12)^2$. By Erdos-Turan and Lemma 3.3 in Kuipers & Neiderrater 'Uniform Distribution of Sequences', we obtain the following bound of discrepancy:

Lemma

Let $\theta$ be an irrational real number. If $\mu$ is the irrationality measure of $\theta$, then the discrepancy $D_N$ for the sequence $n\theta \ \mathrm{mod} \ 1$ satisfies $$ D_N \ll N^{-\frac1{\mu-1}+\epsilon}. $$

Therefore, we obtain the following estimate:

Let $\theta$ be an irrational real number with irrationality measure $\mu$, then we have $$ \sum_{\nu\leq n} \{ \nu\theta \}^2 = \frac n{12} + O\left( n^{1-\frac1{\mu-1} + \epsilon } \right). $$

Note that $\mu=2$ for almost all $\theta\in\mathbb{R}$ yielding the error term bound $O(n^{\epsilon})$.

Kuipers & Neiderrater's book also contains the following result:

Theorem If $\theta\in\mathbb{R}\backslash \{0\}$ has a bounded partial quotients in the continued fraction, then $$ D_N\ll \frac{\log N}N. $$

For $\theta=\sqrt 7$, it has periodic partial quotients, so we have

$$ \sum_{\nu\leq n} \{ \nu\theta \}^2 = \frac n{12} + O\left( \log n \right) $$ which is better than $O(n^{\epsilon})$.