Not yet reviewed by Pith; the record is open.

The recent detection of gravitational waves (GWs) by pulsar timing array (PTA) collaborations spurred a variety of questions regarding the origin of the signal and the properties of its sources. The amplitude of a GW background produced by inspiralling supermassive black holes (SMBHs) can be predicted in a relatively robust manner from the present-day merged remnants, observed as single SMBHs at the centers of galaxies, but falls short of the signal measured by PTAs by a significant amount, requiring equal mass mergers, extremely short delay times, and no accretion in order to achieve a modest consistency. In this work, we revisit NANOGrav's 15-yr data set and reassess the aforementioned discrepancy using the full spectral information captured by PTA data. As previously noted in the literature, the discrete number of point sources contributing to the background may lead to deviations in the observed spectrum relative to the average ($h^2_c \propto f^{-4/3}$) due to Poisson fluctuations, providing additional information about the source population beyond the background amplitude. We derive a simple expression for the characteristic strain distribution given a SMBH model, which is generally applicable regardless of the method used to model the black hole population. We then refit the NANOGrav free spectrum using a minimal model based on the local mass function, showing that the current GW measurement requires roughly $\sim 10$ times more black holes than suggested by local observations and disfavors mass functions dominated by few very heavy sources, with the typical mass that contributes to the background $\lesssim 10^{10}M_{\odot}$. Given the range of SMBH models found to be consistent with the isotropic background, we address what is the typical number sources that would be individually detectable, given the current sensitivity.