we can construct a grand canonical trial state $\Gamma_\tau$ having the density $$\rho_{\Gamma_\tau}(x)=\sum_{z\in{\mathbb{Z}}^3}\sum_{j=1}^{24} \chi_{\ell,\delta,j}(x-\ell z-\tau)\,\rho(x).$$ For this we pick $\Gamma_\tau=\bigotimes_{j=1}^{24}\bigotimes_{z\in{\mathbb{Z}}^3} \Gamma_{\tau,z,j}$ where each $\Gamma_{\tau,z,j}$ has the density $$\rho_{\Gamma_{\tau,z,j}}(x)=\chi_{\ell,\delta,j}(x-\ell z-\tau)\rho(x)$$ and minimizes the corresponding energy $E(\chi_{\ell,\delta,j}(\cdot-\ell z-\tau)\rho)$. Since the quantum states $\Gamma_{\tau,z,j}$ have disjoint supports, we can anti-symmetrize the state $\Gamma_\tau$ in the standard manner. We denote by $\Gamma_{\tau,a}$ the anti-symmetrized state. The energy of $\Gamma_{\tau,a}$ is equal to that of $\Gamma_\tau$ and so is its density $\rho_{\Gamma_{\tau,a}}=\rho_{\Gamma_\tau}$. Finally we take as trial state $$\Gamma=\frac{1}{\ell^3}\int_{C_\ell}\Gamma_{\tau,a}\,d\tau$$ which satisfies by construction that $\rho_\Gamma=\rho$. We find the upper bound $$\begin{aligned}
E(\rho)&{\leqslant}\frac{1}{\ell^3}\int_{C_{\ell}}E(\rho_{\Gamma_\tau})\,d\tau+\frac{1}{\ell^3}\int_{C_\ell}D(\rho_{\Gamma_\tau})\,d\tau-D(\rho){\nonumber}\\
&=\frac{1}{\ell^3}\int_{C_\ell}\sum_{z\in{\mathbb{Z}}^3}\sum_{j=1}^{24} E\big(\chi_{\ell,\delta,j}(\cdot-\ell z-\tau)\rho\big)\,d\tau+\frac{1}{\ell^3}\int_{C_\ell}D(\rho_{\Gamma_\tau}-\rho)\,d\tau.
\label{eq:almost-final_upper_bound}\end{aligned}$$ Here we employ the usual notation $$D(\rho):=\frac12\int_{{\mathbb{R}}^3}\int_{{\mathbb{R}}^3}\frac{\rho(x)\rho(y)}{|x-y|}dx\,dy.
\label{eq:def_D_rho_rho}$$ In the second line of  we have used that the energy of a tensor product of states of disjoint supports is the sum of the energies of the pieces [@LewLieSei-18]. This is because the cross terms in the direct energy exactly cancel with the many-particle interactions of different states in the tensor product.

The error term in  is solely due to the nonlinearity of the direct term and it may be rewritten as $$\frac{1}{\ell^3}\int_{C_\ell}D(\rho_{\Gamma_\tau}-\rho)\,d\tau\\
=\frac{1}{\ell^3}\int_{C_\ell}D\left(\sum_{z\in{\mathbb{Z}}^3}\sum_{j=1}^{24}({{\ensuremath {\mathds 1} }}_{\ell{\mathbf{\Delta}}_j}-\chi_{\ell,\delta,j})(\cdot-\ell z-\tau)\rho\right)\,d\tau.$$ For every real-valued $(\ell{\mathbb{Z}}^3)$–periodic function $f$, we have $$\begin{aligned}
&\frac{1}{\ell^3}\int_{C_\ell}D\Big(f(\cdot-\tau)\rho\Big)\,d\tau\\
&\qquad =\frac12\sum_{k\in 