Juhan Raidal, Juan Urrutia, Ville Vaskonen, Hardi Veermäe

We study the non-Gaussian features of the gravitational wave (GW) background generated by a population of inspiraling supermassive black hole (SMBH) binaries. We show that the SMBH GW amplitude distribution (GWAD) features a universal heavy power-law tail $\propto A^{-4}$, while the low-amplitude tail depends on the SMBH merger rate and the energy-loss mechanisms of the binaries. The distribution of the induced timing residuals inherits this heavy tail. As a result, the ensemble averaged statistical moments of order three and higher diverge, limiting their usefulness as measures of non-Gaussianity, and the GW background from SMBH binaries exhibits the single loud source principle, according to which the strongest signals are more likely to be caused by a small number of loud sources. We confirm that the variance-averaged Gaussian approximation accurately describes the timing residual statistics. This approximation justifies a factored likelihood structure that combines standard Gaussian-process PTA posteriors with the non-Gaussian population prior, enabling consistent incorporation of non-Gaussian effects into SMBH model inference. We provide a fast and flexible Python implementation to compute the distribution of timing residuals from a given SMBH merger rate or GWAD.