Consider a system enclosed in a cubic domain, of edge $L$ and volume $V = L^3$. Let $\mathbf{k}_i = \frac{2\pi}{L}(n_{xi},n_{yi},n_{zi})$, where $i = 1,2,3,4$, and $n_{xi},n_{yi},n_{zi} = 0,\pm 1, \pm 2,....$ Let the positive number $N$ be given, and let the number $k_F$ be determined by $$\sum_{0\leq |\mathbf{k}|\leq k_{F}} 1 = N,$$ where the notation $$\sum_{0\leq |\mathbf{k}|\leq k_{F}}$$ means $$\sum_{\mathbf{k}}$$ where $\mathbf{k}$ takes on all possible values, subject to the restriction $0\leq |\mathbf{k}|\leq k_{F}$.
Question: How to calculate the sum $$\sum_{0\leq |\mathbf{k}_{1}|\leq k_{F}} \sum_{0\leq |\mathbf{k}_{2}|\leq k_{F}} \sum_{0\leq |\mathbf{k}_{3}|\leq k_{F}} \sum_{0\leq |\mathbf{k}_{4}|\leq k_{F}} \delta (\mathbf{k}_{1} + \mathbf{k}_{2}, \mathbf{k}_{3} + \mathbf{k}_{4})$$ as a function of $N$, in the thermodynamic limit (i.e., when $N\rightarrow\infty$, $V\rightarrow\infty$, such that $\frac{N}{V} = constant$), where $\delta$ is the Kronecker symbol, i.e., $\delta (\mathbf{a},\mathbf{b}) = 1$ if $\mathbf{a} = \mathbf{b}$, and $\delta (\mathbf{a},\mathbf{b}) = 0$ if $\mathbf{a}\neq\mathbf{b}$.
Edit: In the thermodynamic limit one can replace the sum over $\mathbf{k}$ by an integral over $\mathbf{k}$: $$\sum_{\mathbf{k}} f(\mathbf{k}) \longrightarrow \frac{V}{(2\pi )^3}\int f(\mathbf{k})d^3 k$$ for any integrable function $f$.
Note that the required sum for large $V$ is approximately $$I:=\left(\prod_{i=1}^4\,\frac{V}{(2\pi)^3}\iiint\limits_{\left\|\textbf{k}_i\right\|\leq k_F}\,\textbf{d}^3\textbf{k}_i\right)\,\left(\frac{1}{V}\,\iiint\limits_{\mathbb{R}^3}\,\text{d}^3\textbf{x}\,\exp\big(\text{i}\,\left(\textbf{k}_1+\textbf{k}_2-\textbf{k}_3-\textbf{k}_4\right)\cdot\textbf{x}\big)\right)\,.$$ That is, with $\textbf{y}:=k_F\,\textbf{x}$ and $\textbf{w}_i:=s_i\,\frac{\textbf{k}_i}{k_F}$ for $i=1,2,3,4$ with $s_1=s_1=+1$ and $s_2=s_3=-1$, we get $$I=\frac{1}{(2\pi)^3}\left(\frac{k_F^3V}{(2\pi)^3}\right)^3\,\iiint\limits_{\mathbb{R}^3}\,\text{d}^3\textbf{y}\,\left(\prod_{i=1}^4\,\iiint\limits_{\left\|\textbf{w}_i\right\|\leq 1}\,\text{d}^3\textbf{w}_i\,\exp\big(\text{i}\,\textbf{w}_i\cdot\textbf{y}\big)\right)\,.$$ By symmetry, $$I=\frac{1}{(2\pi)^3}\left(\frac{k_F^3V}{(2\pi)^3}\right)^3\,\iiint\limits_{\mathbb{R}^3}\,\text{d}^3\textbf{y}\,\left(\iiint\limits_{\left\|\textbf{w}\right\|\leq 1}\,\text{d}^3\textbf{w}\,\exp\big(\text{i}\,\textbf{w}\cdot\textbf{y}\big)\right)^4\,.$$ Let $w:=\|\mathbf{w}\|$ and $y:=\|\mathbf{y}\|$. For fixed $\textbf{x},\textbf{w}\neq\boldsymbol{0}$, $\theta\in[0,\pi]$ denote the angle between $\mathbf{w}$ and $\textbf{y}$. Therefore, $$I=\frac{1}{(2\pi)^3}\left(\frac{k_F^3V}{(2\pi)^3}\right)^3\,\int\limits_{0}^\infty\,4\pi y^2\,\text{d}y\,\left(\int\limits_{0}^1\,\text{d}w\,\int\limits_0^\pi\,\text{d}\theta\,\exp\big(\text{i}\,wy\cos(\theta)\big)\,2\pi w^2\sin(\theta)\right)^4\,.$$ Let $t:=\cos(\theta)$. Then, $$\begin{align}I&=\frac{1}{(2\pi)^3}\left(\frac{k_F^3V}{(2\pi)^3}\right)^3\,\int\limits_{0}^\infty\,4\pi y^2\,\text{d}y\,\left(\int\limits_{0}^{1}\,\text{d}w\,\int\limits_{-1}^{+1}\,\text{d}t\,\exp\big(\text{i}\,wyt\big)\,2\pi w^2\right)^4 \\&=\frac{1}{(2\pi)^3}\left(\frac{k_F^3V}{(2\pi)^3}\right)^3\,\int\limits_{0}^\infty\,4\pi y^2\,\text{d}y\,\left(\int\limits_{0}^{1}\,\text{d}w\,\left(\frac{2\sin(wy)}{wy}\right)\,2\pi w^2\right)^4 \\&=\frac{1}{(2\pi)^3}\left(\frac{k_F^3V}{(2\pi)^3}\right)^3\,\int\limits_{0}^\infty\,4\pi y^2\,\text{d}y\,\left(\frac{4\pi}{y}\int\limits_{0}^{1}\,w\sin(wy)\,\text{d}w\right)^4 \\ &=\frac{1}{(2\pi)^3}\left(\frac{k_F^3V}{(2\pi)^3}\right)^3\,\int\limits_{0}^\infty\,4\pi y^2\,\text{d}y\,\left(\frac{4\pi}{y^3}\,\big(\sin(y)-y\cos(y)\big)\right)^4 \\ &=128\pi^2\left(\frac{k_F^3V}{(2\pi)^3}\right)^3\,\int\limits_{0}^\infty\, y^2\,\left(\frac{\sin(y)-y\cos(y)}{y^3}\right)^4\,\text{d}y\,. \end{align}$$ Using Mathematica, we obtain $$\int\limits_0^\infty\, y^2\,\left(\frac{\sin(y)-y\cos(y)}{y^3}\right)^4\,\text{d}y=\frac{17\pi}{2835}\,,\tag{*}$$ so that $$I=128\pi^2\left(\frac{k_F^3V}{(2\pi)^3}\right)^3\left(\frac{17\pi}{2835}\right)=\frac{2176\pi^2}{2835}\left(\frac{k_F^3V}{(2\pi)^3}\right)^3\,.$$
Since the OP wants the spatial integral to be over $[0,L]^3$, then the required sum is well approximated by $$J(L):=\frac{64\pi}{8}\left(\frac{k_F^3V}{(2\pi)^3}\right)^3\,\iiint\limits_{\left[0,k_FL\right]^3}\,\text{d}\phi\,\text{d}\vartheta\,\text{d}y\,\sin(\vartheta)\,y^2\,\left(\frac{\sin(y)-y\cos(y)}{y^3}\right)^4\,.$$ Thus, $$J(L)\geq 8\pi\left(\frac{k_F^3V}{(2\pi)^3}\right)^3\,\left(4\pi\,\int\limits_0^{k_FL}\,y^2\,\left(\frac{\sin(y)-y\cos(y)}{y^3}\right)^4\,\text{d}y\right)$$ and $$J(L)\leq 8\pi\left(\frac{k_F^3V}{(2\pi)^3}\right)^3\,\left(4\pi\,\int\limits_0^{\sqrt{3}\,k_FL}\,y^2\,\left(\frac{\sin(y)-y\cos(y)}{y^3}\right)^4\,\text{d}y\right)\,.$$ Hence, $$J(L)\approx \frac{272\pi^2}{2835}\left(\frac{k_F^3V}{(2\pi)^3}\right)^3+\mathcal{O}\left(\frac{1}{\left(k_FL\right)^5}\right)\,.$$ Using Mathematica, the error part $\mathcal{O}\left(\frac{1}{\left(k_FL\right)^5}\right)$ lies between $-\frac{272\pi^2}{2835}\left(\frac{k_F^3V}{(2\pi)^3}\right)^3\left(\frac{1}{360\sqrt{3}\,\left(k_FL\right)^5}\right)$ and $-\frac{272\pi^2}{2835}\left(\frac{k_F^3V}{(2\pi)^3}\right)^3\left(\frac{1}{40 \left(k_FL\right)^5}\right)$ (for sufficiently large $L$).
If the OP mean to do the integral with the spatial part on $\left[-\frac{L}{2},+\frac{L}{2}\right]^3$, then the answer is $\tilde{J}(L)=8\,J\left(\frac{L}{2}\right)$. We have $$\tilde{J}(L)\approx \frac{2176\pi^2}{2835}\left(\frac{k_F^3V}{(2\pi)^3}\right)^3+\mathcal{O}\left(\frac{1}{\left(k_FL\right)^5}\right)\,.$$ The error $\mathcal{O}\left(\frac{1}{\left(k_FL\right)^5}\right)$ lies between $-\frac{2176\pi^2}{2835}\left(\frac{k_F^3V}{(2\pi)^3}\right)^3\left(\frac{1}{45\sqrt{3}\,L^5}\right)$ and $-\frac{2176\pi^2}{2835}\left(\frac{k_F^3V}{(2\pi)^3}\right)^3\left(\frac{1}{5 \left(k_FL\right)^5}\right)$ (for sufficiently large $L$).
The notation $\mathcal{O}\left(\frac{1}{\left(k_FL\right)^5}\right)$ is a bit deceptive. Since $V=L^3$, we have $J(L)$ and $\tilde{J}(L)$ go like $\left(k_FL\right)^9$ for large $L$. The error terms in fact grow like $\frac{\left(k_FL\right)^9}{\left(k_FL\right)^5}=\left(k_FL\right)^4$.
If we include the approximation of the Kronecke delta by an integral (see the first equation in my answer), then there is a larger error term for the required sum with behavior $\frac{(k_FL)^9}{(k_FL)^2}=\left(k_FL\right)^7$ for large $L$. I have lost interest to make a proper analysis when this Kronecke delta approximation is included. However, the point is, with this much larger error at hand, it is not reasonable to demand precision for the value of $J(L)$ or $\tilde{J}(L)$.
P.S. Maybe some analyst can show how to obtain (*) manually.