Solution (source code)

= Solution

Set $\xi=hk$. For the <discrete Fourier mode>, the displacement difference in the $l$th neighbor pair carries $e^{il\xi}-1$, and
$$
\sum_{l\ne0}\frac{e^{il\xi}-1}{|l|^3}=2\sum_{l=1}^\infty\frac{\cos(l\xi)-1}{l^3}=-F(\xi).
$$
Thus the infinite system reduces to
$$
\boxed{s^2X-2\Omega sY=3\Omega^2X-\frac{Gm}{h^3}F(\xi)X,\qquad s^2Y+2\Omega sX=2\frac{Gm}{h^3}F(\xi)Y.}
$$
The <gravitating lattice Fourier kernel> $F$ is nonnegative, even and $2\pi$-periodic, with $F(0)=0$. The <uniform convergence> and <absolute convergence> of its differentiated series gives
$$
F'(\xi)=2\sum_{l=1}^\infty\frac{\sin(l\xi)}{l^2}.
$$
To prove the maximum, it is not sufficient to maximize each cosine term separately: the even terms do not reach their individual maxima at $\pi$. Instead use the convergent integral representation $l^{-2}=\int_0^\infty t e^{-lt}dt$. Summing the <geometric series> for sines gives
$$
F'(\xi)=2\sin\xi\int_0^\infty\frac{t e^{-t}}{1-2e^{-t}\cos\xi+e^{-2t}}\,dt.
$$
For $0<\xi<\pi$, the denominator is positive and $\sin\xi>0$, so $F'>0$. For $\pi<\xi<2\pi$, $F'<0$. The interchange is justified by $\sum_l\int_0^\infty t e^{-lt}|\sin(l\xi)|dt\le\sum_l l^{-2}<\infty$. Hence
$$
\boxed{F\text{ has its global maxima exactly at }\xi=(2j+1)\pi.}
$$
Their value is
$$
\boxed{F(\pi)=4\sum_{l\text{ odd}}l^{-3}=\frac72\zeta(3),}
$$
where $\zeta$ is the <Riemann zeta function>. This maximizing <Fourier mode> alternates the signs of neighboring particle displacements.