Exact energies of spherical harmonic oscillator in Dirac equation The potential is given by:
$$
V(r) = {1\over 2} \omega^2 r^2
$$
and we are solving the radial Dirac equation (in atomic units):
$$
c{d P(r)\over d r} + c {\kappa\over r} P(r) + Q(r) (V(r)-2mc^2) = E Q(r)
$$
$$
-c{d Q(r)\over d r} + c {\kappa\over r} Q(r) + P(r) V(r) = E P(r)
$$
What is the analytic expression for the eigenvalues $E$ in atomic units?
It is ok to provide source code (any language) to obtain it if one needs to solve some simple analytic equation. Here are the (I think correct) energies from my numerical code that I would like to compare against analytic solution (for $c = 137.03599907$ and $\omega=1$):
  n  l  k kappa        E

  1  0  0 -1       1.49999501
  2  0  0 -1       3.49989517
  2  1  0 -2       2.49997504
  2  1  1  1       2.49993510
  3  0  0 -1       5.49971547
  3  1  0 -2       4.49983527
  3  1  1  1       4.49979534
  3  2  0 -3       3.49994176
  3  2  1  2       3.49987520
  4  0  0 -1       7.49945592
  4  1  0 -2       6.49961564
  4  1  1  1       6.49957571
  4  2  0 -3       5.49976206

If one only solves the radial Schroedinger equation, then the analytic formula is $$E_{nl} = \omega (2n - l - {1\over 2})$$
I am looking for the relativistic version.
I found for example the paper: Qiang Wen-Chao: Bound states of the Klein-Gordon and
Dirac equations for scalar and vector
harmonic oscillator potentials. Vol 11, No 8, 2002, Chin. Phys. Soc., but it only shows a formula for nonzero scalar and vector potentials in Dirac equation (above we only have the scalar potential, the vector potential is zero).
 A: Have you tried to search "Dirac Oscillator"? there's Moshinsky and Szczepaniak article ... 
A: I found the spectrum in the end of the Appendix in the article http://arxiv.org/ftp/arxiv/papers/1203/1203.2458.pdf (the case of $\epsilon=0$, in order to get external field-free solution) Unfortunately I feel a bit lost about what are exact meaning of the scalar and vector potential in the terms of the hamiltonian there. I mean, the hamiltonian will be in the form
$$H = c \cdot \alpha \cdot p + \beta \cdot M c^2 +  \text{WW}$$
And if I will express the potential WW as block-matrix (distinguishing large and small components in the Standard representation), the WW would look like?
$$\begin{pmatrix}
       S(r) & ?V(r)\\
       ?V(r) & -S(r)
\end{pmatrix}$$
Where $S(r) = 1/2 M \omega_s r^2,  V(r) = 1/2 M \omega_v r^2$ ? I suppose the scalar potential to give bilinear form    $\bar \psi S(r) \psi $
transforming as scalar (of course, the $r^2$ should have to be completed by $-c^2 t^2$ to be Lorentz invariant, but for stacionary states it doesn't matter). I am not sure about the vector potential though.
