  We study the point spectrum of the linearization at a solitary wave solution
$\phi_\omega(x)e^{-\mathrm{i}\omega t}$ to the nonlinear Dirac equation in
$\mathbb{R}^n$, $n\ge 1$, with the nonlinear term given by
$f(\psi^*\beta\psi)\beta\psi$ (known as the Soler model). We focus on the
spectral stability, that is, the absence of eigenvalues with nonzero real part,
in the non-relativistic limit $\omega\lesssim m$, in the case when $f\in
C^1(\mathbb{R}\setminus\{0\})$, $f(\tau)=|\tau|^k+O(|\tau|^K)$ for $\tau\to 0$,
with $0<k<K$. For $n\ge 1$, we prove the spectral stability of small amplitude
solitary waves ($\omega\lesssim m$) for the charge-subcritical cases $k\lesssim
2/n$ ($1<k\le 2$ when $n=1$) and for the "charge-critical case" $k=2/n$,
$K>4/n$.
  An important part of the stability analysis is the proof of the absence of
bifurcations of nonzero-real-part eigenvalues from the embedded threshold
points at $\pm 2m\mathrm{i}$. Our approach is based on constructing a new
family of exact bi-frequency solitary wave solutions in the Soler model, using
this family to determine the multiplicity of $\pm 2\omega\mathrm{i}$
eigenvalues of the linearized operator, and the analysis of the behaviour of
"nonlinear eigenvalues" (characteristic roots of holomorphic operator-valued
functions).
