This finite sum involving roots of unity is bounded, but why? For $N \in \mathbb{N}$, let $\omega := \exp(\frac{2 \pi}{N} \sqrt{-1} )$. For $k,l \in \{1, \dots, N \}$, define
$$d_{k,l}^{(N)} := \begin{cases}
 \frac{\sin\left(\frac{2\pi}{N}(k-1)\right)}{\sin\left(\frac{2\pi}{N}(k-1)\right)^2 + \sin\left(\frac{2\pi}{N}(l-1)\right)^2}&\mbox{ if } \sin\left(\frac{2\pi}{N}(k-1)\right) \neq 0\\
 0 & \mbox{ otherwise }
 \end{cases}$$
For $i,j,m,n \in \{1, \dots, N\}$, let
$$a_{i,j,m,n}^{(N)} :=\frac{1}{N^2}  \sum_{l=1}^N \sum_{k=1}^N (\omega^{i-j})^{k-1} (\omega^{m-n})^{l-1} d^{(N)}_{k,l} $$
Numerical experiments very strongly suggest that there exists $C>0$ independent of $N$ ($C \approx 0.25$) such that $\forall N$,
$$
 \max_{1 \le i,j,m,n \le N} | a^{(N)}_{i,j,m,n} | \le C
$$
How can one formally show that this is true?
 A: I eventually found a (fairly lengthy) proof of this fact. I will not write it in detail, but I sketch the ideas, in hopes somebody might eventually find them helpful -- I would be happy to clarify any question in the comments. Let for convenience $p := i-j$, $q := m-n$.

*

*Using Euler's formula $e^{ix} = \cos(x) + i\sin(x)$, expand the terms $\omega^{p(k-1)}$ and $\omega^{q(l-1)}$ as sines and cosines.


*The sum can then be identified with four Riemann sums, for integrals
$$\int_{y=0}^{2\pi}\int_{x=0}^{2\pi} \frac{\sin(x)}{\sin(x)^2 + \sin(y)^2 }h^x(px)h^y(qy) \,dx\, dy $$
for $h^x, h^y \in \{ \sin, \cos\}$, up to multiplicative constants.


*A simple argument of parity shows that three of these four integrals are zero, and the only nonzero is:
$$I_{p,q} := \int_{y=0}^{2\pi}\int_{x=0}^{2\pi} \frac{\sin(x)}{\sin(x)^2 + \sin(y)^2 }\sin(px)\cos(qy) \,dx\, dy $$


*I refer to user TheSimpliFire's excellent proof here that this integral tends to zero as $q\to\infty$. A completely analogous argument shows that it also tends to zero as $p\to\infty$.


*A careful application of Riemann-Lebesgue's lemma to the integrals in the proof above determines that the rate of convergence to zero is $I_{p,q} = O(p^{-1}q^{-2})$, and $I_{p,0} = O(p^{-1})$. The asymmetry with respect to $p$ and $q$ is caused by the following fact: if one takes the innermost integral to be the one with respect to $x$, as in the aforementioned proof, then eventually one arrives to an integral
of the form $\int_0^{2\pi} g(y)\cos(qy)\,dy $, with $g$ smooth and such that $g(0) = g(2\pi)$. Riemann-Lebesgue thus gives a bound of $O(q^{-2})$. Conversely, integrating first by $y$, one reaches an integrand of the form $\int_0^{2\pi} f(x)\sin(px) \,dx$ with $f$ smooth but $f(0)\neq f(2\pi)$; Riemann-Lebesgue only gives a $O(p^{-1})$ bound.


*The integrands $\frac{\sin(x)}{\sin(x)^2 + \sin(y)^2 }h^x(px)h^y(qy)$ are (smooth) periodic functions, integrated over one period. Periodic functions display accelerated convergence of its Riemann sums (e.g. see Guillemin and Stroock's "Some Riemann sums are better than others"); so the error of the Riemann sums, $\epsilon_N$, decays faster than any polynomial of $N (\ge |p|,|q|)$, and can effectively be ignored.


*In conclusion, if $m\neq n$,
$$a_{i,j,m,n}^{(N)} = \begin{cases}\epsilon_N + O(|i-j|^{-1}|m-n|^{-2}) &\text{ if }m\neq n \\ \epsilon_N + O(|i-j|^{-1}) &\text{ if } m=n,\,i\neq j\\
\epsilon_N &\text{ if } m=n,\,i=j\end{cases}$$
Which is obviously bounded for all $i,j,m,n$.
