\[\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...U=\half∫ddr[ρav(\Br)+δρ(\Br)]⋅[ϕav(\Br)+δϕ(\Br)]=≡U\ns0⏞−\half∫ddrρav(\Br)ϕav(\Br)+∫ddrϕav(\Br)ρ(\Br)+ignorefluctuationterm⏞\half∫ddrδρ(\Br)δϕ(\Br). We apply the mean field approximation in each region of space, which leads to
