Solution (source code)

= Solution
{c}

Define the dimensionless coupling $C=\pi|V_0|/E_{\rm ZB}$, so that $\epsilon=C|\delta|$. Substitution into the <energy density> from part (b) gives, within the prescribed continuation of the two-wave <nearly-free electron model>,
$$
\boxed{\mathcal U_{\rm el}(\delta)=\frac{E_{\rm ZB}}a
\left[\frac83-\sqrt{4+C^2\delta^2}
-\frac{C^2\delta^2}{2}\ln\frac{2+\sqrt{4+C^2\delta^2}}{C|\delta|}\right].}
$$
Its value at $\delta=0$ is understood by continuity. For $C|\delta|\ll1$,
$$
\boxed{\mathcal U_{\rm el}(\delta)=\frac{2E_{\rm ZB}}{3a}
-\frac{\pi^2V_0^2}{4aE_{\rm ZB}}\delta^2
\left[1-2\ln\frac{\pi|V_0\delta|}{4E_{\rm ZB}}\right]
+O\left(\frac{E_{\rm ZB}}a C^4\delta^4\right).}
$$
Both <electrons> per doubled <unit cell> have been included. The fixed $2G$ harmonic affects the underlying <energy bands> beyond the two-wave approximation; the expression here is the <energy density> requested from part (b), rather than an exact solution for all harmonics of the full <periodic potential>.