  The many-body wave function of an $N$-electron system within a relativistic
framework can be described by the Dirac equation. Unfortunately, the Dirac
operator $\hat{\mathfrak{D}}$ is unbounded and in case we would describe anions
we will observe the variational collapse of wavefunction of $N+1$ electron.
Thus, it is necessary to avoid it and an alternative approach is based on
applying the square of the Dirac operator $\hat{\mathfrak{D}}^{2}$. This
approach is especially suitable for a multiwavelet framework: its
implementation in an integral equation form is readily available and the
completeness of the multiwavelet basis guarantees that the matrix
representation of $\hat{\mathfrak{D}}^{2}$ in the occupied space is identical
to the square of the matrix representation of $\hat{\mathfrak{D}}$. This
equivalence is not valid in a finite basis set. The use of the
$\hat{\mathfrak{D}}^{2}$ operator within a multiwavelet framework is able to
achieve higher precision compared to the $\hat{\mathfrak{D}}$, consistent with
the quadratic nature of the error estimates in a variational framework.
