We show how spectral functions for quantum impurity models can be calculated very accurately using a complete set of "discarded" numerical renormalization group eigenstates, recently introduced by Anders and Schiller. The only approximation is to judiciously exploit energy scale separation. Our derivation avoids both the overcounting ambiguities and the single-shell approximation for the equilibrium density matrix prevalent in current methods, ensuring that relevant sum rules hold rigorously and spectral features at energies below the temperature can be described accurately. PACS numbers: 71.27.+a, 75.20.Hr, 73.21.La Quantum impurity models describe a quantum system with a small number of discrete states, the "impurity", coupled to a continuous bath of fermionic or bosonic excitations. Such models are relevant for describing transport through quantum dots, for the treatment of correlated lattice models using dynamical mean field theory, or for the modelling of the decoherence of qubits.The impurity's dynamics in thermal equilibrium can be characterized by spectral functions of the type A BC (ω) = dt 2π e iωt B (t)Ĉ T . Their Lehmann representation readswith Z = a e −βEa and E ba = E b − E a , which implies the sum rule dωA BC (ω) = BĈ T . In this Letter, we describe a strategy for numerically calculating A BC (ω) that, in contrast to previous methods, rigorously satisfies this sum rule and accurately describes both high and low frequencies, including ω T , which we test by checking our results against exact Fermi liquid relations. Our work builds on Wilson's numerical renormalization group (NRG) method [1]. Wilson discretized the environmental spectrum on a logarithmic grid of energies Λ −n , (with Λ > 1, 1 ≤ n ≤ N → ∞), with exponentially high resolution of low-energy excitations, and mapped the impurity model onto a "Wilson tight-binding chain", with hopping matrix elements that decrease exponentially as Λ −n/2 with site index n. Because of this separation of energy scales, the Hamiltonian can be diagonalized iteratively: adding one site at a time, a new "shell" of eigenstates is constructed from the new site's states and the M K lowest-lying eigenstates of the previous shell (the so-called "kept" states), while "discarding" the rest.Subsequent authors [2,3,4,5,6,7,8,9,10] have shown that spectral functions such as A BC (ω) can be calculated via the Lehmann sum, using NRG states (kept and discarded) of those shells n for which ω ∼ Λ −n/2 . Though plausible on heuristic grounds, this strategy entails double-counting ambiguities [5] about how to combine data from successive shells. Patching schemes [9] for addressing such ambiguities involve arbitrariness. As a result, the relevant sum rule is not satisfied rigorously, with typical errors of a few percent. Also, the density matrix (DM)ρ = e −βĤ /Z has hitherto been represented rather crudely using only the single N T -th shell for which, with a chain of length N = N T , resulting in inaccurate spectral information for ω T . In this Letter we avoid these probl...