Continuum solutions for the Dirac equation in Coulomb potential - numerical codes

Following the representation used in [1, pag. 11] the solution of the Dirac equation in polar coordinates for energy $E$ is of the type:

$$ psi_{Ekappa m}(bf{r})= dfrac{1}{r} Bigg( begin{matrix}
P_{Ekappa} (r) chi_{kappa}^m(theta,phi)
iQ_{Ekappa} (r) chi_{-kappa}^m(theta,phi)
end{matrix}Bigg) ,$$

where $P$ and $Q$ represent the large and small radial components, respectively. While $chi_{kappa}^m(theta,phi)$ is the spherical spinor function.

The book gives the solutions for the radial equation in terms of the solution of the Kummer’s confluent hypergeometric equation, that is:

$$dfrac{d^2Y(rho)}{drho^2}+(b-rho)dfrac{dY(rho)}{drho}-aY(rho)=0$$.
Thus, the large and small component are given as:

$$P_{Ekappa}proptorho^{gamma}e^{-rho/2} [X(rho)+Y(rho)]$$

$$Q_{Ekappa}proptorho^{gamma}e^{-rho/2} [X(rho)-Y(rho)]$$

where:

$$X(rho)propto biggl(aY(rho)+rhodfrac{dY(rho)}{drho}biggr).$$

and $rho$ is proportional to the radius.

Depending from the type of solution searched (bound or continuum) $P$ and $Q$ assume different shapes.

The code GRASP2K [2] provides the possibility to calculate the Dirac bound wavefunctions, implementing the multiconfiguration Dirac-Hartree-Fock method.

I would like to know if anyone knows a code to calculate the continuum wavefunctions for a given atom or I should try to implement the analytical solutions.

References

[1] Relativistic quantum theory of atoms and molecules, by I. Grant.

[2] https://www-amdis.iaea.org/GRASP2K/