\[\begin{split} U&=\half\!\int\!\!d^d\!r\>\big[\rho^{av}(\Br)+\delta\rho(\Br)\big]\cdot\big[\phi^{av}(\Br) + \delta\phi(\Br)\big]\\ &=\ \stackrel{\equiv\ U\ns_0}{\overbrace{-\half\!\int\!\!d^d\!r\,\rh...\[\begin{split} U&=\half\!\int\!\!d^d\!r\>\big[\rho^{av}(\Br)+\delta\rho(\Br)\big]\cdot\big[\phi^{av}(\Br) + \delta\phi(\Br)\big]\\ &=\ \stackrel{\equiv\ U\ns_0}{\overbrace{-\half\!\int\!\!d^d\!r\,\rho^{av}(\Br)\,\phi^{av}(\Br) }}+ \int\!\!d^d\!r\>\phi^{av}(\Br)\,\rho(\Br) + \stackrel{\scriptstyle{ignore\ fluctuation\ term}}{\overbrace{ \half\!\int\!\!d^d\!r\,\delta\rho(\Br) \,\delta\phi(\Br)}}\ . \end{split}\] We apply the mean field approximation in each region of space, which leads to
