Asymptotic expansion double integral

$\newcommand{\sinc}{\operatorname{sinc}}$ $\newcommand{\phi}{\varphi}$ We are interested in the integral $$\frac{l^2}{(2\pi)^2}\int_{[0, 2\pi]^2}\sinc^2((k-q)l/2)\phi(k, q)dkdq = \frac{l^2}{(2\pi)^2}\int_{[0, 2\pi]^2}\frac{\sinc^2((k-q)l/2)}{\sinc^2((k-q)/2)}\psi(k, q)dkdq,$$ where $\psi(k, q) = n(k)(1-n(q))\frac{e^{\beta(\epsilon(k)-\epsilon(q))}-1}{\beta(\epsilon(k)-\epsilon(q))}$ is some concrete smooth function in the closed square precise form of which would be completely irrelevant to us, and $\phi(k, q) = \frac{\psi(k, q)}{\sinc^2((k-q)/2)}$ which is smooth in the closed square except the points $(2\pi, 0)$ and $(0, 2\pi)$, where it has singularities.

I will show that $I(l) = a_1 l + a_2 \log l + O(1)$, where $a_1 = \frac{1}{2\pi}\int_0^{2\pi}\phi(k, k)dk$ and $a_2 = \frac{1}{32\pi^4}(\psi(2\pi, 0)+\psi(0, 2\pi))$, and then explain a possible approach to getting further terms of its asymptotics.

We will split the integral into four regions: I: $|k-q| < \frac{1}{100}$, II: $2\pi - q + k < \frac{1}{100}$, III: $2\pi - k+q < \frac{1}{100}$, IV: the rest. In the region IV we will use the first formula for the integral above, using that $\phi(k, q)$ is bounded in this region, while $\sinc^2((k-q)l/2) \lesssim \frac{1}{l^2}$, so the integral over this region is $O(1)$.

In the region I function $\phi$ is smooth but the factor $\sinc^2((k-q)l/2)$ is not particularly small. There, the business is happening on the scale $\sim \frac{1}{l}$, so the function $\phi$ does not change much and so the area of that region is $\sim \frac{1}{l}$, thus we get something linear and explicit in $l$, which is exactly what you did in your post.

In the regions II and III (clearly, they are morally identical) the factor $\sinc^2((k-q)l/2)$ is small, but we are near the singularity of $\phi$. In particular, to even ensure that the integral converges we have to assume that $l\in \mathbb{N}$. Specifically, we will additionally split it into the subregions IIa: $2\pi - q + k \le \frac{1}{l}$, IIb: $\frac{1}{l}< 2\pi - q +k < \frac{1}{100}$ and the same for III. In IIa (and IIIa) we have $$\frac{\sinc^2((k-q)l/2)}{\sinc^2((k-q)/2}\lesssim \frac{1}{l^2} \frac{l^2}{1} = 1,$$ where the first term is the ratio of $\frac{1}{x^2}$ in the $\sinc^2$, while the second term is the ratio of squared sines (via a crude bound with the derivative, the derivative of $\sin(t)$ around $\pi$ is proportional to $1$, while the derivative of $\sin(lt)$ is proportional to $l$). Since the measures of the regions IIa and IIIa are $\sim \frac{1}{l^2}$, in total they give us $O(1)$ integral. Finally, in IIb and IIIb we have $\sinc^2((k-q)l/2)$ is close to $\frac{\sin^2((k-q)l/2)}{(\pi l)^2}$, the denominator here is constant, while the average of the numerator is $\frac{1}{2}$ since it is fast oscillating. Finally, $\sinc^2((k-q)/2)$ is about $\frac{1}{(2\pi)^2 (2\pi -q + k)^2}$ in IIb and $\frac{1}{(2\pi)^2 (2\pi - k + q)^2}$ in IIIb, because $\frac{1}{2\pi}$ is the derivative of $\sinc(x/2)$ at $x = 2\pi$.

Doing the change of variables $r = 2\pi - q + k, s = k + q - 2\pi, \frac{1}{l}<r < \frac{1}{100}, -r < s < r$, we get the following value (for IIb): $$\frac{l^2}{(2\pi)^2}\frac{1}{(\pi l)^2}\frac{1}{2}\frac{1}{2}\int_{1/l}^{1/100} \int_{-r}^r \frac{1}{4\pi^2 r^2}dsdr = \frac{1}{32\pi^4}(\log l - \log(100)) = \frac{1}{32\pi^4}\log l + O(1),$$ where the first factor is our prefactor for the integral, the second factor is from the denominator of $\sinc^2((k-q)l/2)$, the third term is the average of sine squared, and the fourth term is the Jacobian of our change of variables. Finally, we have to multiply by $\psi(0, 2\pi)$. By annoying but not very difficult estimates with Taylor series with integral reminder or something similar we can show that all our approximations result in errors which are $O(1)$, and we get the claimed value for $a_2$.

Now, if we want to go further than $O(1)$, we have to get a better control of all the $O(1)$ errors that we've accumulated so far, which I think (although I am not completely sure and I did not do the relevant computations) might be doable with more Taylor expansions. Note that once we care about $a_3$ region IV starts to play role, so while $a_1$ only depends on $\phi(k, k)$, $a_2$ depends on $\psi(2\pi, 0)$ and $\psi(0, 2\pi)$, $a_3$ must depend on the whole function $\phi$.