4-component Relativistic Calculations in a Multiwavelet Basis with Improved Convergence

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 non-relativistic one, yet in a four-component framework, where Multiwavelet tools and algorithms developed for the non-relativistic 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 \textit{GRASP} code otherwise.