Solution (source code)

= Solution

In the totally symmetric $(r,s)$ basis, the <force constant matrix> reduces to
$$
K_g=\begin{pmatrix}\gamma&-\gamma\\-\gamma&2k+\gamma\end{pmatrix}
=\begin{pmatrix}0&0\\0&2k\end{pmatrix}
+\gamma\begin{pmatrix}1&-1\\-1&1\end{pmatrix}.
$$
The unperturbed eigenvalues are distinct within this symmetry block, so <nondegenerate perturbation theory> applies even though the full unperturbed matrix has degeneracies in different symmetry sectors. The first-order shifts are both $\gamma$. At second order the off-diagonal matrix element gives $-\gamma^2/(2k)$ for the lower mode and $+\gamma^2/(2k)$ for the upper one. Therefore
$$
\boxed{\begin{aligned}
\omega_-^2&=\frac1m\left(\gamma-\frac{\gamma^2}{2k}+\cdots\right),\\
\omega_+^2&=\frac1m\left(2k+\gamma+\frac{\gamma^2}{2k}+\cdots\right).
\end{aligned}}
$$
In frequency form,
$$
\omega_-=\sqrt{\frac\gamma m}\left(1-\frac\gamma{4k}+\cdots\right),\qquad
\omega_+=\sqrt{\frac{2k}m}\left(1+\frac\gamma{4k}+\cdots\right).
$$
The corresponding eigenvectors are $r+(\gamma/2k)s+\cdots$ and $s-(\gamma/2k)r+\cdots$. As a check, the exact squared-frequency eigenvalues are
$$
m\omega_\pm^2=k+\gamma\pm\sqrt{k^2+\gamma^2},
$$
whose expansions reproduce the <perturbation theory> result. The newly lifted zero frequency scales as $\sqrt\gamma$, so expanding $\omega_-$ as an ordinary Taylor series in $\gamma$ would be inappropriate.