given by Wick’s formula $$\begin{gathered}
\Gamma_\gamma^{(2)}(x_1,\sigma_1,x_2,\sigma_2;y_1,\sigma'_1,y_2,\sigma'_2)\\
=\gamma(x_1,\sigma_1;y_1,\sigma'_1)\gamma(x_2,\sigma_2;y_2,\sigma_2')-\gamma(x_1,\sigma_1;x_2,\sigma_2)\gamma(y_1,\sigma'_1;y_2,\sigma'_2).
\label{eq:2PDM_quasi_free_state}\end{gathered}$$ In particular, the corresponding interaction energy with pair potential $w$ is $$\frac12 \iint_{{\mathbb{R}}^d\times{\mathbb{R}}^d}w(x-y)\left(\rho(x)\rho(y)-\sum_{\sigma,\sigma'=1}^q|\gamma(x,\sigma;y,\sigma')|^2\right)\,dx\,dy.$$ From this we immediately obtain the following.

\[cor:upper\_bound\_E\_rho\] We have $$E(\rho){\leqslant}q^{-\frac2d}c_{\rm TF}(1+\kappa_1{\varepsilon})\int_{{\mathbb{R}}^d} \rho(x)^{1+\frac2d}\,dx+\frac{\kappa_2(1+\sqrt{\varepsilon})^2}{{\varepsilon}}\int_{{\mathbb{R}}^d}|\nabla\sqrt{\rho}(x)|^2\,dx
\label{eq:upper_bound_E}$$ for every $\rho{\geqslant}0$ such that $\sqrt{\rho}\in H^1({\mathbb{R}}^d)$.

This is for the grand-canonical version  of the Levy-Lieb functional which is the object of concern in this paper. It was proved in [@Lieb-81] that any fermionic $\gamma$ with integer trace $N={{\rm Tr}\,}(\gamma)$ is also the one-particle density matrix of an $N$-particle mixed state $\Gamma$ on the fermionic space ${\mathfrak{H}}^N$, such that the corresponding two-particle density matrix satisfies $$\Gamma^{(2)}{\leqslant}\Gamma^{(2)}_\gamma$$ in the sense of operators. From the positivity of the Coulomb potential we deduce immediately the following result for mixed canonical states.

Let $\rho\in L^1({\mathbb{R}}^d,{\mathbb{R}}_+)$ be such that $\int_{{\mathbb{R}}^d}\rho=N\in\mathbb{N}$, and $\sqrt{\rho}\in H^1({\mathbb{R}}^d)$. Then there exists a mixed state $\Gamma$ on the fermionic space $\bigwedge_1^NL^2({\mathbb{R}}^d\times\{1,...,q\})$ such that $\rho_\Gamma=\rho$ and $$\begin{gathered}
{{\rm Tr}\,}\left(\sum_{j=1}^N-\Delta_{x_j}+\sum_{1{\leqslant}j<k{\leqslant}N}\frac{1}{|x_j-x_k|}\right)\Gamma\\
{\leqslant}q^{-\frac2d}c_{\rm LT}(1+\kappa_1{\varepsilon})\int_{{\mathbb{R}}^d} \rho(x)^{1+\frac2d}\,dx+\frac{\kappa_2(1+\sqrt{\varepsilon})^2}{{\varepsilon}}\int_{{\mathbb{R}}^d}|\nabla\sqrt{\rho}(x)|^2\,dx.\end{gathered}$$

Proof of Theorems \[thm:UEG\_thermo\_limit\_quantum\] and \[thm:LDA\] {#sec:proof}
=====================================================================

Our proof is divided into several steps. The first is to show that the energy is essentially local, that is, to prove that $$E(\rho)\approx \sum_k  E\big(\rho\chi_k\big)
\label{eq:energy_is_almost_local}$$ where $\{\chi_k\}$ is a smooth partition of unity, $\sum_k\chi_k=1$. The 