Despite more than 40 years of research in condensed-matter physics, state-of-the-art approaches for simulating the radial distribution function (RDF) g(r) still rely on binning pair-separations into a histogram. Such methods suffer from undesirable properties, including subjectivity, high uncertainty, and slow rates of convergence. Moreover, such problems go undetected by the metrics often used to assess RDFs. To address these issues, we propose (I) a spectral Monte Carlo (SMC) method that yields g(r) as an analytical series expansion; and (II) a Sobolev norm that assesses the quality of RDFs by quantifying their fluctuations. Using the latter, we show that, relative to histogram-based approaches, SMC reduces by orders of magnitude both the noise in g(r) and the number of pair separations needed for acceptable convergence. Moreover, SMC reduces subjectivity and yields simple, differentiable formulas for the RDF, which are useful for tasks such as coarse-grained force-field calibration via iterative Boltzmann inversion. In simulations of condensed matter systems, one can barely overstate the importance of the radial distribution function (RDF) g(r). To name only a few applications, g(r) is used to (i) link thermodynamic properties to microscopic details [1-3]; (ii) compute structure factors for comparison with X-ray diffraction [4,5]; and more recently, (iii) calibrate interparticle forces for coarse-grained (CG) molecular dynamics (MD) [6][7][8][9][10][11]. Indeed, the RDF is such a key property that in the past few years, much work has been devoted to estimating g(r) via parallel processing on GPUs [12]. Given these observations, it is thus surprising that state-of-the-art techniques still construct g(r) by binning simulated pair-separations into histograms, with little thought given to developing more efficient methods [3,13].In this letter, we address this issue by proposing a spectral Monte Carlo (SMC) method for computing simulated RDFs. The key idea behind our approach is to express g(r) in an appropriate basis set and determine the mode coefficients via Monte Carlo estimates. Relative to binning, we show that this approach decreases subjectivity of the analysis, thereby reducing both the noise in g(r) and the number of pair separations needed to generate useful RDFs. To support these claims, we also discuss how traditional L 2 (or sum-of-squares) metrics are insufficient for assessing convergence of g(r) and propose a Sobolev norm [14] as an appropriate alternative.The motivation for this work stems from the fact that g(r) is increasingly being used in settings in which the details of its functional form play a critical role. For example, scientists now routinely simulate untested materials in an effort to tailor their structural properties without the need for expensive experiments [15,16]; in such applications, objectively computing RDFs is a key task. Along * paul.patrone@nist.gov related lines, structural properties are increasingly being used to calibrate coarse-grained force-fields [6][7][8][9][10][...