pith. machine review for the scientific record. sign in

arxiv: 2604.08506 · v1 · submitted 2026-04-09 · 🌌 astro-ph.CO · astro-ph.HE· gr-qc

Recognition: unknown

The Heavy Tailed Non-Gaussianity of the Supermassive Black Hole Gravitational Wave Background

Hardi Veerm\"ae, Juan Urrutia, Juhan Raidal, Ville Vaskonen

Authors on Pith no claims yet

Pith reviewed 2026-05-10 16:50 UTC · model grok-4.3

classification 🌌 astro-ph.CO astro-ph.HEgr-qc
keywords supermassive black hole binariesgravitational wave backgroundpulsar timing arraysnon-Gaussianityheavy tailstiming residualssingle loud source
0
0 comments X

The pith

The distribution of gravitational wave amplitudes from supermassive black hole binaries has a universal heavy tail proportional to the inverse fourth power of amplitude.

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The paper shows that the amplitudes of gravitational waves from inspiraling supermassive black hole binaries follow a distribution with a universal heavy power-law tail. This tail carries over to the timing residuals observed in pulsar timing arrays, so that statistical moments above the second order diverge on average. The background therefore obeys the single loud source principle, in which the largest signals arise from a small number of strong sources rather than the combined effect of many weak ones. A variance-averaged Gaussian approximation still holds for the residuals, allowing standard analysis tools to be combined with a non-Gaussian population prior. The authors supply a Python code to compute the residual distribution from any given merger rate model.

Core claim

The SMBH GW amplitude distribution features a universal heavy power-law tail proportional to 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. The variance-averaged Gaussian approximation accurately describes the timing,

What carries the argument

The gravitational wave amplitude distribution (GWAD) of supermassive black hole binaries, which carries a universal heavy tail that sets the non-Gaussian statistics of the background.

If this is right

  • Timing residuals inherit the heavy-tailed distribution of the wave amplitudes.
  • Moments of order three and higher diverge, so they cannot serve as reliable non-Gaussianity diagnostics.
  • The background obeys the single loud source principle, with strongest signals produced by few loud binaries.
  • The variance-averaged Gaussian approximation remains accurate for the residuals.
  • A factored likelihood lets standard Gaussian-process PTA posteriors be combined with the non-Gaussian population prior.

Where Pith is reading between the lines

These are editorial extensions of the paper, not claims the author makes directly.

  • Future PTA observations could search for the predicted power-law tail directly in the distribution of residuals instead of relying on moment statistics.
  • The supplied Python implementation allows rapid checks of how alternate merger-rate histories alter the low-amplitude part of the distribution.
  • Reporting of stochastic-background upper limits in PTA papers may need to account for the possibility that a few loud sources dominate the signal.
  • Similar heavy-tailed behavior could appear in other astrophysical backgrounds whose source populations follow power-law statistics.

Load-bearing premise

The universal heavy tail and its consequences rest on a particular model of the supermassive black hole binary population together with chosen energy-loss mechanisms during inspiral.

What would settle it

Pulsar timing array data that yields finite third- or higher-order moments of the timing residuals, or that shows no power-law tail in the amplitude distribution, would falsify the central claim.

Figures

Figures reproduced from arXiv: 2604.08506 by Hardi Veerm\"ae, Juan Urrutia, Juhan Raidal, Ville Vaskonen.

Figure 1
Figure 1. Figure 1: FIG. 1. The SMBH merger rates of Model I ( [PITH_FULL_IMAGE:figures/full_fig_p004_1.png] view at source ↗
Figure 2
Figure 2. Figure 2: FIG. 2. The solid and dotted curves in the left and middle panels show, respectively, the GWAD for binaries evolving purely [PITH_FULL_IMAGE:figures/full_fig_p005_2.png] view at source ↗
Figure 3
Figure 3. Figure 3: FIG. 3. The distribution of timing residuals generated with [PITH_FULL_IMAGE:figures/full_fig_p007_3.png] view at source ↗
Figure 4
Figure 4. Figure 4: FIG. 4. The distribution of the averaged variances ( [PITH_FULL_IMAGE:figures/full_fig_p008_4.png] view at source ↗
Figure 5
Figure 5. Figure 5: The input of the code is either parameters of the merger rate from Model I or Model II, from which a GWAD is computed, or a GWAD in the form of a broken power law dN dA = Nb(p + q) s h q (A/Ab) p/s + p (A/Ab) q/sis , (35) where Nb is the normalization, p and q are the asymp￾totic powers, Ab is the power-law breaking point, and s determines the smoothness of the break. As discussed in Sec. III B 1, q should… view at source ↗
Figure 6
Figure 6. Figure 6: FIG. 6. Violin plot showing the distributions of variances ( [PITH_FULL_IMAGE:figures/full_fig_p010_6.png] view at source ↗
Figure 7
Figure 7. Figure 7: FIG. 7. Distribution of [PITH_FULL_IMAGE:figures/full_fig_p014_7.png] view at source ↗
Figure 8
Figure 8. Figure 8: FIG. 8. Window function for the mode [PITH_FULL_IMAGE:figures/full_fig_p014_8.png] view at source ↗
Figure 9
Figure 9. Figure 9: FIG. 9. Violin plot showing the distributions of timing residuals obtained by using different window functions. The distri [PITH_FULL_IMAGE:figures/full_fig_p016_9.png] view at source ↗
Figure 10
Figure 10. Figure 10: FIG. 10. Correlations between the Fourier modes, [PITH_FULL_IMAGE:figures/full_fig_p016_10.png] view at source ↗
read the original abstract

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.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, simulated authors' rebuttal, and a circularity audit. Tearing a paper down is the easy half of reading it; the pith above is the substance, this is the friction.

Referee Report

0 major / 2 minor

Summary. The manuscript analyzes the non-Gaussian statistics of the stochastic gravitational wave background from inspiraling supermassive black hole binaries. It shows that the GW amplitude distribution (GWAD) has a universal heavy power-law tail scaling as A^{-4} independent of merger-rate and energy-loss details, while the low-amplitude regime depends on those astrophysical inputs. The induced timing residuals inherit the same heavy tail, causing moments of order three and higher to diverge and enforcing the single-loud-source principle. The authors confirm that a variance-averaged Gaussian approximation accurately captures the residual statistics and supply a fast Python implementation that accepts arbitrary merger-rate or GWAD inputs.

Significance. If the central geometric derivation holds, the result is significant for pulsar-timing-array analyses of the nanohertz gravitational-wave background. It cleanly separates universal Poisson-plus-inverse-distance effects from model-dependent low-amplitude behavior, thereby justifying factored likelihood constructions that combine standard Gaussian-process PTA posteriors with non-Gaussian population priors. The explicit provision of reproducible Python code is a concrete strength that allows immediate community verification and application to specific SMBH population models.

minor comments (2)
  1. The abstract and introduction state that the tail index is exactly -4; a short appendix deriving the index from the Poisson point process and the h ∝ 1/D scaling would make the universality argument fully self-contained without requiring the reader to reconstruct it from the cited literature.
  2. Figure captions and the code documentation should explicitly note the numerical integration limits used to separate the universal tail from the low-amplitude regime so that users can reproduce the plotted distributions with their own merger-rate models.

Simulated Author's Rebuttal

0 responses · 0 unresolved

We thank the referee for their positive and constructive review of our manuscript. Their summary accurately captures the key results, and we are pleased that they recognize the significance for pulsar-timing-array analyses and the value of the provided Python implementation. The recommendation to accept is appreciated.

Circularity Check

0 steps flagged

No significant circularity; universal tail derived from independent geometric statistics

full rationale

The paper's central derivation of the universal heavy tail ∝ A^{-4} in the GW amplitude distribution follows directly from Poisson spatial statistics of sources combined with the amplitude scaling h ∝ 1/D (or equivalent strain scaling). This geometric argument is independent of the specific merger rate or energy-loss mechanisms, which the paper correctly assigns only to the low-amplitude regime. No load-bearing step reduces by construction to a fitted parameter, self-citation, or ansatz imported from the authors' prior work. The inheritance to timing residuals and divergence of higher moments are immediate consequences of the same statistics. The derivation is self-contained against external benchmarks and does not rely on renaming known results or smuggling assumptions via citation.

Axiom & Free-Parameter Ledger

2 free parameters · 2 axioms · 0 invented entities

The central claim rests on standard astrophysical modeling of supermassive black hole binary populations and energy-loss processes; the heavy tail itself is asserted to be universal and therefore independent of the precise values of those inputs.

free parameters (2)
  • SMBH merger rate
    Controls the shape of the low-amplitude tail of the GWAD
  • energy-loss mechanism parameters
    Determine the low-amplitude tail of the amplitude distribution
axioms (2)
  • domain assumption The gravitational wave background is produced by a population of inspiraling supermassive black hole binaries
    Foundational premise stated in the opening sentence of the abstract
  • domain assumption Binaries evolve according to standard energy-loss mechanisms
    Invoked to model the low-amplitude part of the distribution

pith-pipeline@v0.9.0 · 5529 in / 1704 out tokens · 62323 ms · 2026-05-10T16:50:29.812473+00:00 · methodology

discussion (0)

Sign in with ORCID, Apple, or X to comment. Anyone can read and Pith papers without signing in.

Forward citations

Cited by 2 Pith papers

Reviewed papers in the Pith corpus that reference this work. Sorted by Pith novelty score.

  1. A practical theorem on gravitational-wave background statistics

    astro-ph.CO 2026-04 unverdicted novelty 7.0

    For large but finite source counts, the PDF of rescaled GWB characteristic strain squared follows the universal form N^{1/3} times the reflected map-Airy distribution evaluated at N^{1/3}(y-1), fully determined by the...

  2. Are PTA measurements sensitive to gravitational wave non-Gaussianities?

    astro-ph.CO 2026-05 unverdicted novelty 6.0

    PTA statistical tests lose sensitivity to non-Gaussian GW features after decorrelation and cannot distinguish them model-agnostically.

Reference graph

Works this paper leans on

62 extracted references · 46 canonical work pages · cited by 2 Pith papers · 3 internal anchors

  1. [1]

    It originates from the possibility of having nearby sources and can be derived analytically by considering potentialz≪1 bina- ries

    High-amplitude tail The high-Aasymptotic of GWAD is independent of the modeling of the SMBH merger rate. It originates from the possibility of having nearby sources and can be derived analytically by considering potentialz≪1 bina- ries. The probability of finding such a nearby binary is proportional to the area of the shell around the observer, dλ∝D 2 LdD...

  2. [2]

    2, where we also show that the tail varies when the parameter determining the low-Mpower-law of the merger rate is varied

    Low-amplitude tail The emergence of the low-Apower-law can clearly be seen in Fig. 2, where we also show that the tail varies when the parameter determining the low-Mpower-law of the merger rate is varied. This suggests that the low- amplitude tail of the GWAD is largely governed by the low-mass binary population. When ignoring redshift ef- fects, a merge...

  3. [3]

    The width of the band is conditional on the amount of leakage allowed by the window function in Eq

    Strong sources To sample the contributions of the timing residuals, we divide a broad frequency band into a large number Nbins of narrow bins in binary frequency. The width of the band is conditional on the amount of leakage allowed by the window function in Eq. (7). For a given Fourier mode atf k, the amplitudes of theN S strong sources are sampled expli...

  4. [4]

    Weak sources The weak sources are far more numerous but also in- dividually insignificant. Their contribution is included in each binary frequency bin as ˜δt W J,k = NbinsX j=1 h Tjw+ k,j − T ∗ j w− k,j i ,(25) whereT j is a complex random variable whose statistical moments are all finite, and by the central limit theorem, it is expected to follow a compl...

  5. [5]

    (V A)), the statistics of the timing residuals at high| ˜δtJ,k|are dominated by a single loud source, whose amplitude is described by the high-Atail of the GWAD

    Analytic tails Due to the single big jump principle of heavy-tailed distributions (see Sec. (V A)), the statistics of the timing residuals at high| ˜δtJ,k|are dominated by a single loud source, whose amplitude is described by the high-Atail of the GWAD. The distribution in the high-| ˜δtJ,k|tail can be estimated by considering the distribution of timing r...

  6. [6]

    In the K≫1 limit, these intervals correspond tok/K≤1+y < (k+1)/K, wherekis an integer in the range 0≤k <2K

    Formally, it is given by p(|R0|2) = Z 1 −1 dy 2 δ(|R0|2 − |R0|2(y)),(A6) So, we need to invert the mappingy7→ |R 0|by con- sidering intervals of|R 0|where it is one-to-one. In the K≫1 limit, these intervals correspond tok/K≤1+y < (k+1)/K, wherekis an integer in the range 0≤k <2K. For simplicity, we can further assume thatKis an inte- ger, as the contribut...

  7. [7]

    Such terms are removed in the analysis of PTA data [14]

    Low-frequency noise subtraction Deterministic slowly changing components cannot gen- erally be distinguished from low-frequency (LF) back- ground noise. Such terms are removed in the analysis of PTA data [14]. We consider terms with at most quadratic dependence in time, so that δt(t) =δt 0(t)− 2X n=0 An t T n ,(B4) whereδt 0(t) denotes the timing residual...

  8. [8]

    Whitening To mitigate spectral leakage, the time series can first be whitened, i.e., transformed so that the residual noise becomes approximately uncorrelated with unit variance. A whitened time seriesδt W (t) is obtained by convolving the dataδt(t) with a filter kernelW(t): δtW (t) = Z dt′δt(t−t ′)W(t ′).(B9) The Fourier coefficients computed in the whit...

  9. [9]

    9 and 10, respec- tively

    Impact on timing residuals and correlations To illustrate the impact of different window functions on inter-mode leakage, the timing residual PDFs and mode correlations are plotted in Figs. 9 and 10, respec- tively. The correlations, following Eq. (19), depend ex- clusively on the window function, which modulates sig- nal leakage across frequencies, and a...

  10. [10]

    The NANOGrav 15-year Data Set: Evidence for a Gravitational-Wave Background

    G. Agazieet al.(NANOGrav), The NANOGrav 15 yr Data Set: Evidence for a Gravitational-wave Background, Astrophys. J. Lett.951, L8 (2023), arXiv:2306.16213 [astro-ph.HE]

  11. [11]

    Antoniadiset al.(EPTA, InPTA:), Astron

    J. Antoniadiset al.(EPTA, InPTA:), The second data release from the European Pulsar Timing Array - III. Search for gravitational wave signals, Astron. Astrophys. 678, A50 (2023), arXiv:2306.16214 [astro-ph.HE]

  12. [12]

    D. J. Reardonet al., Search for an Isotropic Gravitational-wave Background with the Parkes Pul- sar Timing Array, Astrophys. J. Lett.951, L6 (2023), arXiv:2306.16215 [astro-ph.HE]

  13. [13]

    Searching for the nano-Hertz stochastic gravitational wave background with the Chinese Pulsar Timing Array Data Release I

    H. Xuet al., Searching for the Nano-Hertz Stochastic Gravitational Wave Background with the Chinese Pulsar Timing Array Data Release I, Res. Astron. Astrophys. 23, 075024 (2023), arXiv:2306.16216 [astro-ph.HE]

  14. [14]

    Rajagopal and R

    M. Rajagopal and R. W. Romani, Ultralow frequency gravitational radiation from massive black hole binaries, Astrophys. J.446, 543 (1995), arXiv:astro-ph/9412038

  15. [15]

    J. S. B. Wyithe and A. Loeb, Low - frequency gravita- tional waves from massive black hole binaries: Predic- tions for LISA and pulsar timing arrays, Astrophys. J. 590, 691 (2003), arXiv:astro-ph/0211556

  16. [16]

    Sesana, F

    A. Sesana, F. Haardt, P. Madau, and M. Volonteri, Low 16 −8.8 −8.6 −8.4 −8.2 −8.0 −7.8 −7.6 log10( f [Hz]) −8.5 −8.0 −7.5 −7.0 −6.5 −6.0 −5.5 log10(|˜δt|[s]) GW only −8.8 −8.6 −8.4 −8.2 −8.0 −7.8 −7.6 log10( f [Hz]) GW + env Sinc LF noise subtracted Top-hat Whitened Power law −13/6 RMS mean FIG. 9. Violin plot showing the distributions of timing residuals...

  17. [17]

    Agazieet al.(NANOGrav Collaboration), Astrophys

    G. Agazieet al.(NANOGrav), The NANOGrav 15 yr Data Set: Constraints on Supermassive Black Hole Bina- ries from the Gravitational-wave Background, Astrophys. J. Lett.952, L37 (2023), arXiv:2306.16220 [astro-ph.HE]

  18. [18]

    Antoniadiset al.(EPTA, InPTA), Astron

    J. Antoniadiset al.(EPTA, InPTA), The second data release from the European Pulsar Timing Array - IV. Implications for massive black holes, dark matter, and the early Universe, Astron. Astrophys.685, A94 (2024), arXiv:2306.16227 [astro-ph.CO]

  19. [19]

    Ellis, M

    J. Ellis, M. Fairbairn, G. H¨ utsi, J. Raidal, J. Urru- tia, V. Vaskonen, and H. Veerm¨ ae, Gravitational waves from supermassive black hole binaries in light of the NANOGrav 15-year data, Phys. Rev. D109, L021302 (2024), arXiv:2306.17021 [astro-ph.CO]

  20. [20]

    Afzalet al.(NANOGrav), Astrophys

    A. Afzalet al.(NANOGrav), The NANOGrav 15 yr Data Set: Search for Signals from New Physics, Astrophys. J. Lett.951, L11 (2023), [Erratum: Astrophys.J.Lett. 971, L27 (2024), Erratum: Astrophys.J. 971, L27 (2024)], arXiv:2306.16219 [astro-ph.HE]. 17

  21. [21]

    What is the source of the PTA GW signal?,

    J. Ellis, M. Fairbairn, G. Franciolini, G. H¨ utsi, A. Iovino, M. Lewicki, M. Raidal, J. Urrutia, V. Vaskonen, and H. Veerm¨ ae, What is the source of the PTA GW sig- nal?, Phys. Rev. D109, 023522 (2024), arXiv:2308.08546 [astro-ph.CO]

  22. [22]

    S. R. Taylor, J. Simon, and L. Sampson, Constraints On The Dynamical Environments Of Supermassive Black-hole Binaries Using Pulsar-timing Arrays, Phys. Rev. Lett.118, 181102 (2017), arXiv:1612.02817 [astro- ph.GA]

  23. [23]

    S. R. Taylor, The Nanohertz Gravitational Wave As- tronomer, arXiv preprint (2021), arXiv:2105.13270 [astro-ph.HE]

  24. [24]

    Domcke, G

    V. Domcke, G. Franciolini, and M. Pieroni, Cosmic Vari- ance in Anisotropy Searches at Pulsar Timing Arrays, arXiv preprint (2025), arXiv:2508.21131 [astro-ph.CO]

  25. [25]

    Sesana, A

    A. Sesana, A. Vecchio, and M. Volonteri, Gravitational waves from resolvable massive black hole binary sys- tems and observations with Pulsar Timing Arrays, Mon. Not. Roy. Astron. Soc.394, 2255 (2009), arXiv:0809.3412 [astro-ph]

  26. [26]

    P. A. Rosado, A. Sesana, and J. Gair, Expected proper- ties of the first gravitational wave signal detected with pulsar timing arrays, Mon. Not. Roy. Astron. Soc.451, 2417 (2015), arXiv:1503.04803 [astro-ph.HE]

  27. [27]

    L. Z. Kelley, L. Blecha, L. Hernquist, A. Sesana, and S. R. Taylor, Single Sources in the Low-Frequency Grav- itational Wave Sky: properties and time to detection by pulsar timing arrays, Mon. Not. Roy. Astron. Soc.477, 964 (2018), arXiv:1711.00075 [astro-ph.HE]

  28. [28]

    C. M. F. Mingarelli, B. Larsen, E. Eisenberg, Q. Zheng, and F. Hutchison, Fingerprints of Individual Supermas- sive Black Hole Binaries in Pulsar Timing Arrays, arXiv preprint (2026), arXiv:2603.05722 [astro-ph.HE]

  29. [29]

    S. R. Taylor, R. van Haasteren, and A. Sesana, From Bright Binaries To Bumpy Backgrounds: Mapping Real- istic Gravitational Wave Skies With Pulsar-Timing Ar- rays, Phys. Rev. D102, 084039 (2020), arXiv:2006.04810 [astro-ph.IM]

  30. [30]

    E. C. Gardiner, L. Z. Kelley, A.-M. Lemke, and A. Mitridate, Beyond the Background: Gravitational- wave Anisotropy and Continuous Waves from Supermas- sive Black Hole Binaries, Astrophys. J.965, 164 (2024), arXiv:2309.07227 [astro-ph.HE]

  31. [31]

    Exploring the spectrum of stochastic gravitational-wave anisotropies with pulsar timing arrays,

    G. Sato-Polito and M. Kamionkowski, Exploring the spectrum of stochastic gravitational-wave anisotropies with pulsar timing arrays, Phys. Rev. D109, 123544 (2024), arXiv:2305.05690 [astro-ph.CO]

  32. [32]

    Raidal, J

    J. Raidal, J. Urrutia, V. Vaskonen, and H. Veerm¨ ae, Statistics of supermassive black hole gravitational wave background anisotropy, Astron. Astrophys.706, A159 (2026), arXiv:2411.19692 [astro-ph.CO]

  33. [33]

    Ellis, M

    J. Ellis, M. Fairbairn, G. H¨ utsi, M. Raidal, J. Urru- tia, V. Vaskonen, and H. Veerm¨ ae, Prospects for fu- ture binary black hole gravitational wave studies in light of PTA measurements, Astron. Astrophys.676, A38 (2023), arXiv:2301.13854 [astro-ph.CO]

  34. [34]

    Allen and S

    B. Allen and S. Valtolina, Pulsar timing array source ensembles, Phys. Rev. D109, 083038 (2024), arXiv:2401.14329 [gr-qc]

  35. [35]

    W. G. Lamb and S. R. Taylor, Spectral Variance in a Stochastic Gravitational-wave Background from a Bi- nary Population, Astrophys. J. Lett.971, L10 (2024), arXiv:2407.06270 [gr-qc]

  36. [36]

    Sato-Polito and M

    G. Sato-Polito and M. Zaldarriaga, Distribution of the gravitational-wave background from supermas- sive black holes, Phys. Rev. D111, 023043 (2025), arXiv:2406.17010 [astro-ph.CO]

  37. [37]

    X. Xue, Z. Pan, and L. Dai, Non-Gaussian statistics of nanohertz stochastic gravitational waves, Phys. Rev. D 111, 043022 (2025), arXiv:2409.19516 [astro-ph.CO]

  38. [38]

    W. G. Lamb, J. M. Wachter, A. Mitridate, S. C. Sardesai, B. B´ ecsy, E. L. Hagen, S. R. Taylor, and L. Z. Kelley, Fi- nite Populations & Finite Time: The Non-Gaussianity of a Gravitational Wave Background, arXiv e-print (2025), arXiv:2511.09659 [gr-qc]

  39. [39]

    Lentati, M

    L. Lentati, M. P. Hobson, and P. Alexander, Bayesian Estimation of Non-Gaussianity in Pulsar Timing Anal- ysis, Mon. Not. Roy. Astron. Soc.444, 3863 (2014), arXiv:1405.2460 [astro-ph.IM]

  40. [40]

    R. C. Bernardo, S. Appleby, and K.-W. Ng, Toward a test of Gaussianity of a gravitational wave background, JCAP01, 017, arXiv:2407.17987 [astro-ph.CO]

  41. [41]

    Falxa and A

    M. Falxa and A. Sesana, Modeling non-Gaussianities in pulsar timing array data analysis using Gaussian mixture models, Phys. Rev. D113, 043047 (2026), arXiv:2508.08365 [astro-ph.IM]

  42. [42]

    Kuntz, C

    A. Kuntz, C. Smarra, and M. Vaglio, Looking for non-gaussianity in Pulsar Timing Arrays through the four point correlator, arXiv preprint (2026), arXiv:2603.12311 [gr-qc]

  43. [43]

    Raidal, J

    J. Raidal, J. Urrutia, V. Vaskonen, and H. Veerm¨ ae, GWADpy (2026)

  44. [44]

    Coles, G

    W. Coles, G. Hobbs, D. J. Champion, R. N. Manchester, and J. P. W. Verbiest, Pulsar timing analysis in the pres- ence of correlated noise, Mon. Not. Roy. Astron. Soc. 418, 561 (2011), arXiv:1107.5366 [astro-ph.IM]

  45. [45]

    Raidal, J

    J. Raidal, J. Urrutia, V. Vaskonen, and H. Veerm¨ ae, Ec- centricity effects on the supermassive black hole gravita- tional wave background, Astron. Astrophys.691, A212 (2024), arXiv:2406.05125 [astro-ph.CO]

  46. [46]

    P. J. Armitage and P. Natarajan, Accretion during the merger of supermassive black holes, Astrophys. J. Lett. 567, L9 (2002), arXiv:astro-ph/0201318

  47. [47]

    Merritt, Loss-cone dynamics, Classical and Quan- tum Gravity30, 244005 (2013), arXiv:1307.3268 [astro- ph.GA]

    D. Merritt, Loss-cone dynamics, Classical and Quan- tum Gravity30, 244005 (2013), arXiv:1307.3268 [astro- ph.GA]

  48. [48]

    L. Z. Kelley, L. Blecha, and L. Hernquist, Massive Black Hole Binary Mergers in Dynamical Galactic Environ- ments, Mon. Not. Roy. Astron. Soc.464, 3131 (2017), arXiv:1606.01900 [astro-ph.HE]

  49. [49]

    Y. Tang, A. MacFadyen, and Z. Haiman, On the orbital evolution of supermassive black hole binaries with cir- cumbinary accretion discs, Mon. Not. Roy. Astron. Soc. 469, 4258 (2017), arXiv:1703.03913 [astro-ph.HE]

  50. [50]

    W. H. Press and P. Schechter, Formation of galaxies and clusters of galaxies by selfsimilar gravitational condensa- tion, Astrophys. J.187, 425 (1974)

  51. [51]

    J. R. Bond, S. Cole, G. Efstathiou, and N. Kaiser, Ex- cursion set mass functions for hierarchical Gaussian fluc- tuations, Astrophys. J.379, 440 (1991)

  52. [52]

    Lacey and S

    C. Lacey and S. Cole, Merger rates in hierarchical models of galaxy formation, MNRAS262, 627 (1993)

  53. [53]

    Girelli, L

    G. Girelli, L. Pozzetti, M. Bolzonella, C. Giocoli, F. Marulli, and M. Baldi, The stellar-to-halo mass re- lation over the past 12 Gyr: I. Standard ΛCDM model, Astron. Astrophys.634, A135 (2020), arXiv:2001.02230 [astro-ph.CO]. 18

  54. [54]

    A. E. Reines and M. Volonteri, Relations Between Cen- tral Black Hole Mass and Total Galaxy Stellar Mass in the Local Universe, Astrophys. J.813, 82 (2015), arXiv:1508.06274 [astro-ph.GA]

  55. [55]

    Middleton, W

    H. Middleton, W. Del Pozzo, W. M. Farr, A. Sesana, and A. Vecchio, Astrophysical constraints on massive black hole binary evolution from Pulsar Timing Ar- rays, Mon. Not. Roy. Astron. Soc.455, L72 (2016), arXiv:1507.00992 [astro-ph.CO]

  56. [56]

    Agazieet al., Astrophys

    G. Agazieet al.(NANOGrav), The NANOGrav 15 yr Data Set: Bayesian Limits on Gravitational Waves from Individual Supermassive Black Hole Binaries, Astrophys. J. Lett.951, L50 (2023), arXiv:2306.16222 [astro-ph.HE]

  57. [57]

    Sesana, A

    A. Sesana, A. Vecchio, and C. N. Colacino, The stochastic gravitational-wave background from massive black hole binary systems: implications for observations with Pulsar Timing Arrays, Mon. Not. Roy. Astron. Soc.390, 192 (2008), arXiv:0804.4476 [astro-ph]

  58. [58]

    Agazieet al., Astrophys

    G. Agazieet al., The NANOGrav 15 yr Data Set: Looking for Signs of Discreteness in the Gravitational- wave Background, Astrophys. J.978, 31 (2025), arXiv:2404.07020 [astro-ph.HE]

  59. [59]

    Nanohertz Gravitational Waves,

    A. Sesana and D. G. Figueroa, Nanohertz Gravitational Waves, arXiv preprint (2025), arXiv:2512.18822 [astro- ph.CO]

  60. [60]

    V. P. Chistyakov, A theorem on sums of independent pos- itive random variables and its applications to branching random processes, Theory of Probability & Its Applica- tions9, 640 (1964)

  61. [61]

    S. Foss, D. Korshunov, and S. Zachary,An Introduction to Heavy-Tailed and Subexponential Distributions, 2nd ed., Springer Series in Operations Research and Financial Engineering (Springer, New York, 2013)

  62. [62]

    R. W. Hellings and G. S. Downs, Upper limits on the isotropic gravitational radiation background from pulsar timing analysis, Astrophys. J. Lett.265, L39 (1983)