Fewer supermassive binary black holes in pulsar timing array observations
Boris Goncharov Shubhit Sardana A. Sesana J. Antoniadis A. Chalumeau D. Champion S. Chen E. F. Keane G. Shaifullah L. Speri






frame=tb,
  	language=Matlab,
  	aboveskip=3mm,
  	belowskip=3mm,
  	showstringspaces=false,
  	columns=flexible,
  	basicstyle=,
  	numbers=none,
  	numberstyle=,
 	keywordstyle=,
	commentstyle=,
  	stringstyle=,
  	breaklines=true,
  	breakatwhitespace=true
  	tabsize=3


 

 boris.goncharov@me.com
Max Planck Institute for Gravitational Physics (Albert Einstein Institute), 30167 Hannover, Germany
Leibniz Universität Hannover, 30167 Hannover, Germany

 
Max Planck Institute for Gravitational Physics (Albert Einstein Institute), 30167 Hannover, Germany
Department of Physics, IISER Bhopal, Bhauri Bypass Road, Bhopal, 462066, India








Dipartimento di Fisica “G. Occhialini", Universitá degli Studi di Milano-Bicocca, Piazza della Scienza 3, I-20126 Milano, Italy
INFN, Sezione di Milano-Bicocca, Piazza della Scienza 3, I-20126 Milano, Italy
INAF - Osservatorio Astronomico di Brera, via Brera 20, I-20121 Milano, Italy


FORTH Institute of Astrophysics, N. Plastira 100, 70013, Heraklion, Greece
Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, 53121, Bonn, Germany


ASTRON, Netherlands Institute for Radio Astronomy, Oude Hoogeveensedijk 4, 7991 PD Dwingeloo, The Netherlands


Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, 53121, Bonn, Germany


Shanghai Astronomical Observatory, Chinese Academy of Sciences, Shanghai 200030, P. R. China


School of Physics, Trinity College Dublin, College Green, Dublin 2, D02 PN40, Ireland


Dipartimento di Fisica “G. Occhialini", Universitá degli Studi di Milano-Bicocca, Piazza della Scienza 3, I-20126 Milano, Italy
INFN, Sezione di Milano-Bicocca, Piazza della Scienza 3, I-20126 Milano, Italy
INAF - Osservatorio Astronomico di Cagliari, via della Scienza 5, 09047 Selargius (CA), Italy


Max Planck Institute for Gravitational Physics (Albert Einstein Institute), Am Mühlenberg 1, 14476 Potsdam, Germany







§ ABSTRACT 






















We reanalyse the second data release of the European Pulsar Timing Array (EPTA) using an observationally-driven model for ensemble properties of pulsar noise. 

We show that the revised gravitational wave background properties are in better agreement with theoretical expectations for the strain spectrum. 


Our improved model for ensemble pulsar noise properties reduces a systematic error at 1σ level and increases Bayesian odds of Hellings-Downs correlations by ∼ 10%.



Fewer supermassive binary black holes in pulsar timing array observations
    L. Speri
    September 9, 2024
=========================================================================




Since 2020, Pulsar Timing Arrays (PTAs) have reported growing evidence for the nanohertz-frequency gravitational wave background in their data.
The first tentative evidence came from a temporally-correlated stochastic process in pulsar timing data <cit.>. 
The Fourier spectrum of delays and advances in pulsar pulse arrival times resembled the expected spectral properties of the background. 
Most recently, PTAs — with varying levels of statistical significance <cit.> — showed that this stochastic process exhibits Hellings-Downs correlations <cit.> consistent with the isotropic unpolarised stochastic gravitational wave background.

Supermassive black hole binaries at subparsec separations are expected to be a dominant source of the stochastic gravitational wave background at nanohertz frequencies <cit.>.

However, the expected amplitude of the background is lower than the observations suggest. 
Previous EPTA analyses found the best-fit strain amplitude to be at the edge of the values simulated from supermassive black hole binary population synthesis models.
This is visible in Figure 7 from <cit.>[This is not apparent in the analogous Figure 5 from <cit.> due to the binary population modelling differences.].


The inferred strain amplitude lies at the theoretical upper limit of the predicted astrophysical range <cit.>.
It might also be in tension with the observed black hole mass function <cit.>, although see <cit.> for a different view.
Furthermore, the strain spectral index of the gravitational wave background is in ≈ 2σ tension with the value corresponding to binary inspirals driven by gravitational wave emission alone. 
This is visible in Figure 5 in  <cit.> and Figure 11 in  <cit.>.
Overall, although previous PTA results were consistent with a very broad range of assumptions about supermassive binary black hole populations <cit.>, they suggested deviations from purely gravitational wave-driven binary evolution. 



It was pointed out in several studies <cit.> that the standard PTA models of how noise parameters are distributed across pulsars are incorrect. 
These models manifest as prior probabilities in PTA data analysis.
To be precise, the models are incorrect because they are `static', i.e., the shape of the distribution of noise parameters is not influenced by the data. 
Although imposing such priors is very unlikely to influence our conclusions about evidence for the gravitational background <cit.>, it may introduce systematic errors in our measurement of the strain spectrum <cit.>. 



In this Letter, we address the aforementioned problem.
We employ hierarchical Bayesian inference to account for the uncertainties in our noise priors. 
In particular, we parametrise noise prior distributions which are our models of how parameters governing pulsar-specific noise are distributed across the pulsars <cit.>. 
These newly-introduced hierarchical parameters are called hyperparameters. 
Furthermore, we use a new procedure of marginalisation over hyperparameters, which is described in Section 2.2.2 of the companion paper (Goncharov et al., in prep.). 




Based on the new procedure, we revisit the measurement of properties of the gravitational wave background with the European Pulsar Timing Array (EPTA) <cit.>. 

We assess the properties of the gravitational wave background using the power law model of its characteristic strain spectrum: h_c(f)=A (f yr^-1)^-α.
Here, A is the strain amplitude, and -α is the power law spectral index[One may also find a notation where α'=2/3 and γ=3-2α'.]. 
The corresponding spectral index of the power spectral density [s^3] of delays(-advances) [s] induced by the background in the timing data is -γ. 
These stochastic timing delays resulting from the gravitational-wave redshift of pulsar spin frequency are referred to as temporal correlations. 
The value of γ=13/3 (α=2/3) corresponds to purely gravitational wave driven inspirals of supermassive black hole binaries <cit.>. 

Modelling of pulsar noise is important because it can affect the conclusions about the properties of the gravitational wave background.
Pulsar noise introduces temporal correlations which may be difficult to distinguish from the effect of gravitational background based on single pulsar data. 
However, unlike for the gravitational background, noise-induced timing delays have different Fourier spectra across pulsars. 
Our hierarchical model describes the ensemble properties of these noise spectra. 

Despite similar delay time series induced by both pulsar noise and gravitational background in pulsar data, the presence of the background can still be established. 

The Hellings-Downs function <cit.> determines inter-pulsar correlations of the stochastic timing delays induced by the isotropic stochastic unpolarised gravitational background. 

The statistical significance of the signal is established based on the Bayesian odds between (1) the Hellings-Downs-correlated stochastic process with the same spectrum of delays across pulsars, and (2) the same spectrum of delays across pulsars and no inter-pulsar correlations. 

Inter-pulsar correlations of the gravitational wave background do not contribute as much to a measurement of its strain spectrum as temporal correlations of the signal <cit.>. 
This is visible in, e.g., Figure 6 in ref <cit.> or Figure 3 in ref. <cit.>. 
Therefore, on one hand, even in the lack of evidence of Hellings-Downs correlations per se, the background amplitude and spectral index can still be well-constrained. 
This is the case for the full 25-year EPTA data which lacks evidence for inter-pulsar correlations, and this was the case for earlier studies that have identified a common-spectrum stochastic process in the PTA data <cit.>. 
On the other hand, improving the measurement of pulsar noise could assist in resolving Hellings-Downs correlations. 

Results.–Posterior distributions for A and γ are shown as contours in Figure <ref>.
Solid blue contours correspond to our improved model.
For comparison, dashed red contours correspond to the standard `static' noise priors used in earlier EPTA analyses <cit.>. 
The results are shown for both the 10-year subset of the EPTA data which showed evidence for the Hellings-Downs correlations, and the full 25-year EPTA data where the evidence is not visible[The correlations are thought to be  `scrambled' by unmodelled noise from the older backend-receiver systems of the telescopes.]. 
The value of γ=13/3 (α=2/3) is shown as a dashed straight line. 
Our improved model results in a lower median-aposteriori strain amplitude of the gravitational wave background, as well as in a steeper spectral index, as shown with solid blue contours in Figure <ref>. 


A detailed inspection of Figure <ref> reveals that the impact is most significant for the full 25-year data, where the  maximum-aposteriori amplitude (best fit) also shifts outside of 1σ levels of fully-marginalized distributions of A and γ, peaking exactly at γ=13/3. 
For the 10-year data, the value of γ=13/3 now lies at the edge of the 1σ credible level, but the maximum-aposteriori value remains almost unaffected. 














































Implications.–It is common to assume that the energy loss in binary inspirals is dominated by the emission of gravitational waves. 
In this case, the characteristic strain spectrum of the gravitational wave background is <cit.> 

    h_c^2(f) = 4G^5/3/3π^1/3 c^2 f^-4/3∫d^2N/dV dzℳ^5/3/(1+z)^1/3 dz,

where (G,c) are the universal constants, z is redshift, ℳ is the binary chirp mass, and d^2N/(dVdz), a function of (ℳ,z), is the number density of binaries per unit comoving volume per unit redshift. 
The integral does not depend on a gravitational wave frequency, thus h_c∝ f^-2/3, as stated earlier.

The background amplitude A depends on the mass spectrum and the abundance of supermassive binary black holes in the universe. 
The α=2/3 (γ=13/3) dependence is confirmed by population synthesis simulations, e.g. Figure 7 in <cit.>, where the theoretical uncertainty is only δγ∼ 0.1 at 1σ due to cosmic variance <cit.>.

Tensions of the previously-estimated γ≈ 3 with 13/3 reported during the announcement of evidence for the gravitational background <cit.> has led to discussions on whether the signal is influenced by certain effects of binary evolution that make strain spectrum to appear flatter. 
Mechanisms of flattening h_c(f) typically involve the introduction of a more rapid physical mechanism of binary hardening[A reduction in binary separation.] compared to a gravitational wave emission at <0.1 parsec separations. 
Such a mechanism could be an environmental effect such as stellar scattering <cit.>, the torques of a circumbinary gas disc <cit.>.
It could also be due to the abundance of binaries in eccentric orbits which lead to a more prominent gravitational wave emission <cit.>[Eccentricity also results in a steeper h_c(f>10^-8 Hz) <cit.>, but PTA sensitivity declines towards high frequencies.].

In contrast, the results of our improved analysis maintain consistency with binary evolution driven only by the emission of gravitational waves. 







Our improved model also impacts the measurement of the strain amplitude which can be recast in terms of the number density of supermassive black hole binaries. 
The new results hint that the supermassive black hole binaries are not as (over-)abundant as the earlier measurements suggested. 

It is visible in Figure <ref>, which is a replica of Figure A1 in <cit.>. 
There, green horizontal bands correspond to theoretical uncertainties on the strain amplitude at the 16th - 84th percentile level in 24 studies <cit.>. 



A consideration of ensemble noise properties of pulsars reduces tensions of the gravitational wave background strain amplitude with theoretical and observationally-based predictions for supermassive black hole binaries. 
The caveat is that the reported amplitude is referenced to f=yr^-1, the covariance between A and γ in a posterior changes following a rotation of a power law about this frequency.
Therefore, visible consistency with theoretical predictions may further improve based on a different refrence frequency. 





Observational impact.–
We find that our improved model yields an increase in evidence for Hellings-Downs correlations. 
For the 10-year data, we find an increase in Bayesian odds by 22%, and in the 25-year data - by 6%. 

While the increase in the signal-to-noise ratio (SNR)[Roughly, log Bayesian odds ∝SNR^2 <cit.>.] is almost negligible, consistency of the increase between the 10- and the 25-year data suggests that we have removed a systematic error arising from prior misspecification. 

When the model closely matches reality, one would expect a reduction of the measurement uncertainty when adding extra data. 
This is not visible in the original EPTA analysis where the 1σ range for ( A, γ) is roughly the same between the 10-year and the 25-year data. 
As shown in Figure <ref>, our improved analysis yields a larger measurement uncertainty based on temporal correlations in the 10-year data, in agreement with our expectation. 
This is no longer visible in Figure <ref> suggesting that inter-pulsar correlations present in the 10-year data provide additional constraints. 
The shift of the posteriors towards larger spectral indices and smaller amplitudes with our improved analysis[The covariance between A and γ in Figure <ref> is along the line of equal noise power.] suggests that the louder pulsar-intrinsic noise with flatter spectra leaks into our measurement of the background strain spectrum when ensemble pulsar noise properties are not modelled. 








 



Because the 10-year data and the 25-year data are not independent data sets, a high degree of consistency is expected. 
In the original EPTA analysis, maximum-aposteriori ( A, γ) in the 25-year data differ from those of the 10-year data by around (0.3, 0.8).
It is visible in red contours across all three panels in Figure <ref>. 
A match of the best-fit ( A, γ) between the 25-year data that does not exhibit Hellings-Downs correlations (Figure <ref>) and the 10-year data when modelling only temporal correlations and not the Hellings-Downs correlations (Figure <ref>) is achieved with our improved analysis. 

However, our improved analysis does not strongly impact the 10-year data when modelling inter-pulsar correlations. 
It supports our previous statement that inter-pulsar pulsar correlations in the 10-year data provide additional constraints on ( A, γ). 
Because ( A, γ) obtained with only temporal correlations is still expected to match those obtained with temporal and Hellings-Downs correlations, it is also possible that the EPTA data contains other systematic errors that are yet to be mitigated. 


Hypothesising where other systematic errors may stem from, we note that the North American Nanohertz Observatory for Gravitational Waves (NANOGrav) has mitigated a tension of the background spectral index with 13/3 by adopting the Gaussian process model of the dispersion variation noise <cit.>, as in the EPTA analysis <cit.>. 
Therefore, one potential source of a systematic error could be the mismodelling of the pulsar-specific noise that depends on a radio frequency. 
A very nearby binary is another example <cit.>. 
Frequency-wise comparison of the inferred strain spectrum against black hole population synthesis models performed earlier by the EPTA (Figure 3 in ref. <cit.>) suggests that the deviation from γ=13/3 may occur due to excess noise in two frequency bins, ≈ 10^-8 Hz and ≈ 3 × 10^-8 Hz. 
The rest of the spectrum seems to be consistent with γ=13/3. 
The aforementioned potential sources of systematic errors may require better temporal and inter-pulsar correlation models of the data as part of future work.

Methods.–PTAs perform precision measurements of pulse arrival times from millisecond radio pulsars. 
The likelihood of delays(-advances) δ t for a vector of pulse arrival times t is a multivariate Gaussian distribution ℒ(δ t | θ), where θ is a vector of parameters of models that describe the data. 
From the Bayes theorem, it follows that the posterior distribution of model parameters outlining our measurement is 𝒫(θ|δ t) = 𝒵^-1ℒ(δ t | θ) π(θ), where 𝒵 is Bayesian evidence, a fully-marginalised likelihood. 
The term π(θ) is called prior, a model of how likely it is to find a certain value of θ in Nature. 
Model selection is performed based on computing the ratio of 𝒵 for pairs of models, it is referred to as the Bayes factor. 
The Bayes factor is equal to the Bayesian odds ratio if both models are assumed to have equal prior odds. 
A form of the likelihood and the description of the standard analysis methodology can be found in <cit.>. 
Because the total PTA noise prior is a product of noise priors for every pulsar, PTA data will inform on the distribution of θ in Nature.
Our improved analysis introduces hyperparameters Λ to parametrise priors: π(θ|Λ) π(Λ). 
We then perform a numerical marginalisation over Λ. 
For more details, please refer to the companion paper (Goncharov et al., in prep.).





Data availability.–Second data release of the European Pulsar Timing Array <cit.> is available at https://zenodo.org/record/8091568zenodo.org and https://gitlab.in2p3.fr/epta/epta-dr2gitlab.in2p3.fr. 

Acknowledgements.–We thank Bruce Allen for helpful comments on the results, contributions to the structuring of this manuscript, and valuable advice on scientific paper writing. 
We also thank Rutger van Haasteren for many insightful discussions about hierarchical Bayesian inference. 
Some of our calculations were carried out using the OzSTAR Australian national facility (high-performance computing) at Swinburne University of Technology. 
European Pulsar Timing Array is a member of the International Pulsar Timing Array <cit.>. 

J. A. acknowledges support from the European Commission (ARGOS-CDS; Grant Agreement number: 101094354).
A. C. acknowledges financial support provided under the European Union's Horizon Europe ERC Starting Grant “A Gamma-ray Infrastructure to Advance Gravitational Wave Astrophysics” (GIGA; Grant Agreement: 101116134).



apsrev4-2





