We revive an approach to solve the Dirac equation originally proposed by Kutzelnigg which makes use of the squared Dirac operator \(\hat{\mathfrak {D}}^{2}\) . This approach holds the promise to avoid the negative energy solution because the negative energy spectrum is now “folded” on the positive energy side and at the same time provides a convex equation, which is amenable to a minimization process and increased precision in the final result. The \(\hat{\mathfrak {D}}^{2}\) yields an equation similar to the nonrelativistic one, yet in a four-component framework, where Multiwavelet tools and algorithms developed for the nonrelativistic case can be employed with minor modifications. On the other hand, the use of Multiwavelets is here essential to achieve the full potential of the approach. We implemented and validated this approach for one- and two-electron systems with increasing nuclear charge. Numerical tests were performed to gauge the actual precision of the approach with respect to either analytical reference values when possible or numerical results obtained with the GRASP code otherwise.
