Recursive computation of matrix elements in the numerical renormalization group

The numerical renormalization group is an efficient method to diagonalize model Hamiltonians describing correlated orbitals coupled to conduction states. While only the resulting eigenvalues are needed to calculate the thermodynamical properties for such models, matrix elements of Fermi operators must be evaluated before excitation and transport properties can be computed. The traditional procedure to calculate matrix elements is typically as expensive as the diagonalization of the model Hamiltonian. Here, we present a substantially faster alternative that demands much less memory, yields equally accurate matrix elements and is easier to code.