Solution (source code)

= Solution

Set $G=G_1$, $\kappa=\Delta k$, $U=|V_1|>0$ and $E_B=\hbar^2G^2/(8m)$. The two kinetic <energies> satisfy
$$
\frac{\epsilon_1+\epsilon_2}{2}=E_B+\frac{\hbar^2\kappa^2}{2m},\qquad
\frac{\epsilon_1-\epsilon_2}{2}=\frac{\hbar^2G\kappa}{2m}.
$$
Consequently the <nearly-free electron dispersion near a one-dimensional band gap> becomes
$$
E_\pm\left(\frac G2+\kappa\right)=V_0+E_B+\frac{\hbar^2\kappa^2}{2m}
\pm\sqrt{U^2+\left(\frac{\hbar^2G\kappa}{2m}\right)^2}.
$$
The <Taylor expansion> $\sqrt{U^2+s^2}=U+s^2/(2U)+O(s^4/U^3)$ gives
$$
\boxed{E_\pm=V_0+E_B\pm U+
\left(\frac{\hbar^2}{2m}\pm\frac{\hbar^4G^2}{8m^2U}\right)(\Delta k)^2
+O((\Delta k)^4).}
$$
For the paper's <energy> convention, set $V_0=0$. There is no linear or cubic term because each branch is an even function of displacement from the crossing in this two-wave approximation. Each coefficient multiplying $(\Delta k)^2$ has units of <energy> times length squared, as required.

The expansion requires
$$
\left|\frac{\hbar^2G\Delta k}{2m}\right|\ll U,
$$
a more restrictive condition than merely $|\Delta k|\ll G$ when the <band gap> is small. If $U=0$, the square root is proportional to $|\Delta k|$ and this analytic expansion of the energy-ordered branches is unavailable.