  We present a new estimator for the self-energy based on a combination of two
equations of motion and discuss its benefits for numerical renormalization
group (NRG) calculations. In challenging regimes, NRG results from the standard
estimator, a ratio of two correlators, often suffer from artifacts: the
imaginary part of the retarded self-energy is not properly normalized and, at
low energies, overshoots to unphysical values and displays wiggles. We show
that the new estimator resolves the artifacts in these properties as they can
be determined directly from the imaginary parts of auxiliary correlators and do
not involve real parts obtained by Kramers-Kronig transform. Furthermore, we
find that the new estimator yields converged results with reduced numerical
effort (for a lower number of kept states) and thus is highly valuable when
applying NRG to multiorbital systems. Our analysis is targeted at NRG
treatments of quantum impurity models, especially those arising within
dynamical mean-field theory, but most results can be straightforwardly
generalized to other impurity or cluster solvers.
