Given some fixed integer $n$ (but one may take $n = 0$ for simplicity), I would like to compute the following limit:
$$\lim_{N \to \infty}\sum_{\substack{|u_1|, |u_2|, |u_3|, |u_4| \leq N\\|u_1 + u_2 + u_3 + u_4 - n| \leq N}}\frac{1}{\left(u_1 - \frac{1}{2}\right)\left(u_2 - \frac{1}{2}\right)\left(u_3 - \frac{1}{2}\right)\left(u_4 - \frac{1}{2}\right)}$$
I conjecture numerically and I have sound analytic (yet somehow indirect) arguments that the result is $\frac{\pi^4}{5}$. I have managed to treat the easier case with two variables instead of four, that is to find
$$\lim_{N \to \infty}\sum_{\substack{|u_1|, |u_2| \leq N\\|u_1 + u_2 - n| \leq N}}\frac{1}{\left(u_1 - \frac{1}{2}\right)\left(u_2 - \frac{1}{2}\right)} = -\frac{\pi^2}{3}$$
In the latter case, I used that $\frac{1}{u_1 - \frac{1}{2}}$ is antisymmetric with respect to the change of variable $u_1 \to 1 - u_1$ and the resulting cancellation allowed me to restrict the domain of summation and to eventually regard the sum as a double Riemann sum converging (all changes of variables done) to the integral:
$$2\int_{\substack{0 \leq x \leq 1\\-1 \leq y \leq x - 1}}\!\mathrm{d}x\mathrm{d}y\,\frac{1}{xy} = -\frac{\pi^2}{3}$$
The four-dimensional case (and a fortiori the higher ones) seem more intricate to me though... Ideally, I would like to get an error term besides the actual limit. I am pretty sure this has something to do with the Riemann $\zeta$ functions and Bernoulli numbers (by the way, I wondered whether using Fourier analysis -notably the Fourier expansion of the Bernoulli polynomials- could help, but I have got nothing convincing in this direction up to now).
Which technique(s) do you suggest I used to achieve this?
EDIT:
It seems that in the two and four variable cases (denote by $k$ the number of variables), the sought limit coincides with: $$(2\pi i)^k\int_0^1\!\mathrm{d}x\,B_1(x)^k$$ I think this is the general result, but do you see any rationale for that?