How to compute the connected spectral form factor for integrable models?

Given a spectrum of $N$ real eigenvalues, ${E_m }$ of some Hermitian operator, the connected spectral form factor is defined as follows:
begin{align}
K_c(t) = langle sum_{m,n=1}^N e^{it (E_m-E_n)} rangle – |langle sum_{m=1}^N e^{it (E_m)} rangle|^2.
end{align}
For random matrices, the average $langle .. rangle$ is over the appropriate ensemble (GOE, GUE, GSE etc). However, for integrable models, I am not sure how to take this average. For example, for a particle in a rectangular box of sides $L_x$ and $L_y$, the energy eigenvalues (in some appropriate units) are generated by the formula
begin{align}
E(n_x,n_y) = frac{n_x^2}{L_x^2}+frac{n_y^2}{L_y^2}.
end{align}
But if I wanted to compute the connected SFF, do I average of random choices of $L_x$, $L_y$? Is there a guiding principle for what the right thing to do is? Are there papers where a computation of this kind has been done before? The only one I could find was this paper by Jens Marklof https://people.maths.bris.ac.uk/~majm/bib/spectral.pdf which is a little hard to follow.

Any help would be appreciated.