Introduction

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 data1,2,3,4. The Fourier spectrum of delays and advances in pulsar pulse arrival times exhibited the expected spectral properties of the background. Most recently, PTAs — with varying levels of statistical significance5,6,7,8 — showed that this stochastic process exhibits Hellings-Downs correlations9 consistent with the isotropic unpolarised stochastic gravitational wave background.

Properties of the gravitational wave background are assessed using the power law model of its characteristic strain spectrum: \({h}_{{{\rm{c}}}}(f)=A{(f/{{\mbox{yr}}}^{-1})}^{-\alpha }\). Here, A is the strain amplitude at the gravitational wave frequency f of fref  = yr−1, and  − α is the power law spectral index. The corresponding spectral index of the power spectral density [s3] 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.

Supermassive black hole binaries at subparsec separations are expected to be a dominant source of the stochastic gravitational wave background at nanohertz frequencies10. However, the expected amplitude of the background is lower than the observations suggest11. This is visible in Figure 7 from ref. 12 and Fig. 5 from ref. 13, where the best-fit strain amplitude is at the edge of the values simulated from supermassive black hole binary population synthesis models. The inferred strain amplitude lies at the theoretical upper limit of the predicted astrophysical range10,14,15. It might also be in tension with the observed black hole mass function16,17, although see18 for a different view. Furthermore, the strain spectral index of the gravitational wave background is in about 2σ tension with the value corresponding to binary inspirals driven by gravitational wave emission alone (γ = 13/3, α = 2/3)19. The tension is visible in Fig. 5 in ref. 13 and Fig. 11 in ref. 5, where the posterior on the spectral index is compared not with a single value, but with the distribution of values expected from realistic backgrounds made up of discrete point sources. Although previous PTA results were consistent with a very broad range of assumptions about binary black hole populations20, they suggested deviations from purely gravitational wave-driven binary evolution. Furthermore, deviations of the strain spectral index from  − 2/3 can point to the early-Universe origin of the signal13,21.

Modelling of pulsar-intrinsic noise is important because it can affect the conclusions about the properties of the gravitational wave background, such as A and γ22. It is relevant for the European Pulsar Timing Array (EPTA), one of the world’s leading PTAs, where the latest 10-year data set showed evidence for the gravitational wave background at 3.5σ6. Despite dedicated noise modelling studies being performed for the EPTA23, there are four indications that some noise is still mismodelled. First, it was pointed out in several studies2,4,24 that the standard PTA models of how noise parameters are distributed across pulsars are incorrect. These models manifest as prior probabilities in PTA data analyses. 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 background25,26,27,28, it may introduce systematic errors in the inferred strain spectrum29. Second, epoch-correlated temporally-uncorrelated (white) noise term30 was neglected in the previous EPTA analysis. Third, van Haasteren31 pointed out that the procedure of removing certain noise terms in the search for gravitational waves based on the results of single-pulsar noise analysis, performed in the original EPTA analysis6, is prone to systematic errors. Fourth, the original EPTA analysis treats the transient noise effect in pulsar PSR J1713 + 0747 as from a sudden change in dispersion measure, although Goncharov et al.32 suggested that it is associated with the pulse shape change.

In this work, we address the above four limitations and revise the properties of the gravitational wave background inferred by the European Pulsar Timing Array (EPTA)6. As a result, we find a steeper characteristic strain spectrum of the background, which is in better agreement with the hypothesis that the background originates from a superposition of adiabatically inspiraling supermassive binary black holes in circular orbits. Based on the revised amplitude and spectral index of the background, we discuss implications for the evolution of supermassive black hole binaries. Finally, we describe the observational impact of our new noise model. In particular, we show that despite our new model results in less evidence for Hellings-Downs correlations, evidence remains strong. Our results highlight the importance of accurate noise models for correctly inferring background properties.

Results

Posterior distributions for A and γ are shown as contours in Fig. 1. Solid blue contours correspond to our improved model. For comparison, dashed red contours correspond to the noise model used in the original analysis of the EPTA data6. 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 visible33. 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 Fig. 1. Maximum-aposteriori (lgA, γ) obtained with our improved model and the respective measurement uncertainties based on 1σ credible levels are reported in the caption of Fig. 1.

Fig. 1: Posterior distributions for the power-law amplitude A and spectral index γ of the putative gravitational wave background in the European Pulsar Timing Array (EPTA) data.
figure 1

The results of a fit to Hellings-Downs correlations are shown in the left panel (a, full 25-year data, \(({\mathrm{lg}}\,A,\gamma )=(-14.5{3}_{-0.27}^{+0.15},4.1{4}_{-0.26}^{+0.56})\)) and the middle panel (b, the 10-year subset of the data described in ref. 6, \(({\mathrm{lg}}\,A,\gamma )=(-14.3{7}_{-0.57}^{+0.17},3.3{8}_{-0.25}^{+1.50})\)). The results of a fit of only temporal correlations to the 10-year data are shown in the right panel (c, \(({\mathrm{lg}}\,A,\gamma )=(-14.5{2}_{-0.80}^{+0.15},4.2{0}_{-0.72}^{+1.32})\)). Dashed red contours correspond to the result using the standard pulsar noise priors, and the blue contours correspond to our improved model. The horizontal dashed line corresponds to the background from supermassive binary black holes inspiralling entirely due to gravitational wave emission (subject to cosmic variance). Our improved model results in a lower median-aposteriori strain amplitude of the background and mitigates tensions with γ = 13/3. In marginalised distributions, shaded areas correspond to 1σ credible levels. In joint distributions, an inner dark area corresponds to the 1σ level, and the outer lighter area corresponds to the 2σ level.

Full size image

A detailed inspection of the contribution of four new components of our noise models, not shown in Fig. 1, reveals the following. Using noise model averaging instead of model selection, as recommended by ref. 31, has mostly affected the measurement uncertainty. Every other component of our improved noise model results in a consistently lower background strain amplitude and a steeper spectral index. Resolving noise prior misspecification using hierarchical inference24, the first component of our improved noise model, has mostly affected the 25-year data (Fig. 1a) and the 10-year data when inter-pulsar correlations are not modelled (Fig. 1c). When modelling Hellings-Downs correlations in the 10-year data (Fig. 1b), other components of our revised noise model have made most of the impact. In particular, we found that the amplitude of the transient noise event in PSR J1713+0747 depends on the radio frequency consistently with ref. 32, not as assumed in the original EPTA analysis6. Modelling a transient noise event in PSR J1713 + 0747 correctly results in a significant shift in the posterior. Finally, epoch-correlated noise has not been considered in previous EPTA analyses because of a low number (1–4) of pulse arrival time measurements per observation epoch compared to other PTAs6,23. Nevertheless, we find evidence for this noise in the data from certain backend-received systems. In Fig. 1b, the choice leads to the consistency of the fully marginalised posterior on γ with 13/3 at the 1σ level. Revising noise priors alone, without considering epoch-correlated noise, yields a weaker consistency of the joint posterior on (lgA, γ) with γ = 13/3 at the 1σ level (not shown in Fig. 1b). Stronger impact of hierarchical inference in Fig. 1a, c compared to Fig. 1b illustrates how inter-pulsar correlations influence the posterior.

Implications for supermassive binary black holes

If the energy loss in binary inspirals is dominated by the adiabatic emission of gravitational waves and if binaries are circular, the characteristic strain spectrum of the gravitational wave background is ref. 19

$${h}_{\,{\mbox{c}}\,}^{2}(f)=\frac{4{G}^{5/3}}{3{\pi }^{1/3}{c}^{2}}{f}^{-4/3}\int\frac{{d}^{2}N}{dVdz}\frac{{{{\mathcal{M}}}}^{5/3}}{{(1+z)}^{1/3}}dz,$$
(1)

where G is the Newton’s constant, c is the speed of light, z is redshift, \({{\mathcal{M}}}\) is the binary chirp mass, and d2N/(dVdz), a function of \(({{\mathcal{M}}},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 hc ∝ 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 derived value of α = 2/3 (γ = 13/3) is confirmed by population synthesis simulations, e.g., Figure 7 in ref. 12, where the theoretical uncertainty is only δγ of approximately 0.1 at 1σ due to cosmic variance34,35.

Tensions of γ of 3 inferred from the previous EPTA analysis with 13/3 = 4.3(3) reported during the announcement of evidence for the gravitational background5,6,7,8 has led to discussions on whether the signal is influenced by certain effects of binary evolution that make the strain spectrum to appear flatter. Mechanisms of flattening hc(f) typically involve the introduction of a more rapid physical mechanism of a reduction in binary separation compared to a gravitational wave emission at  < 0.1 parsec separations. Such a mechanism could be an environmental effect36, such as stellar scattering37,38, the torques of a circumbinary gas disc39. Furthermore, it could be due to the abundance of binaries in eccentric orbits, which leads to a more prominent gravitational wave emission40. Although eccentricity also results in a steeper hc(f > 10−8Hz)41, but PTA sensitivity declines towards high frequencies. In contrast, the results of our improved analysis are consistent with binary evolution driven only by the emission of gravitational waves.

Our improved model also impacts the measurement of the strain amplitude, which is directly related to the number density of supermassive black hole binaries. The new results suggest that the supermassive black hole binaries are not as (over-)abundant as the earlier measurements implied. This is shown in Fig. 2, the bottom panel of which is based on Figure A1 in ref. 12. Grey horizontal bands correspond to theoretical uncertainties on the strain amplitude at the 16th - 84th percentile level in 26 studies10,15,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65. For ref. 50 the model is “HS-nod”, and for ref. 51 the model is “HS_nod_SN_high_accr”. Our improved model reduces tensions of the gravitational wave background strain amplitude with theoretical and observationally-based predictions for supermassive black hole binaries. Furthermore, the revised background strain amplitude corresponds to longer delay times between galaxy mergers and supermassive binary black hole mergers51. The caveat is that the reported amplitude is referenced to f = yr−1, so the covariance between lgA and γ in posterior changes following a rotation of a power law about this frequency. We tested that at f = 10yr−1 our improved model mostly affects the posterior on γ, introducing consistency with 13/3, leaving the posterior on lgA10yr almost unchanged. Precisely, with 10-year EPTA data, we measure (lgA10yr, γ) to be \((-14.0{2}_{-0.11}^{+0.08},4.3{5}_{-1.02}^{+0.42})\) with our new noise model and \((-14.0{0}_{-0.13}^{+0.08},3.0{2}_{-0.25}^{+0.89})\) with the original model. Consistently, one may notice in Figure A1 in ref. 12 that the posterior on lgA10yr referenced to f = 10yr−1 are in less tension with theoretical predictions.

Fig. 2: Consistency of the inferred strain amplitude of the gravitational wave background with theoretical expectations.
figure 2

Bottom panels show the predicted log-10 strain amplitude lgA from 26 studies. The top panels show posteriors on lgA marginalised over γ. Red lines correspond to the original noise model, and blue lines correspond to our revised model. Solid lines correspond to the 10-year data, and dashed lines correspond to the 25-year data. Blue bands correspond to 1σ credible levels for 10-year data and our improved noise model. Panel (a) shows the characteristic strain amplitude at yr−1, lgA. Panel (b) shows the strain amplitude at the frequency of (10yr)−1.

Full size image

As shown in ref. 16, Equation (1) can be presented as a function of black hole number density ρBH, black hole mass scale M*, and the intrinsic scatter ϵ0 in galaxy stellar velocity dispersion, which is a proxy of mass for supermassive black holes,

$${h}_{\,{\mbox{c}}\,}^{2}(f)= 1.2\times 1{0}^{-30}{\left(\frac{f}{{{\mbox{yr}}}^{-1}}\right)}^{-4/3}\times {\left(\frac{{M}_{*}}{5.8\times 1{0}^{7}{M}_{\odot }}\right)}^{2/3} \\ \times \left(\frac{{\rho }_{{{\rm{BH}}}}}{4.5\times 1{0}^{5}{M}_{\odot }{{\mbox{Mpc}}}^{-3}}\right)\times \left(\frac{2}{{e}^{\frac{8}{9}{\epsilon }_{0}^{2}{\ln }^{2}(10)}}\right).$$
(2)

The mass scale M* is expressed as \({\mathrm{lg}}\,{M}_{*}={a}_{\bullet }+{b}_{\bullet }\frac{{\sigma }_{*}}{200{{\mbox{km s}}}^{-1}}\). It is also worth noting that ρBH ∝ M*. Based on the equation above, we perform parameter estimation on (ρBH, a•, b•, σ*) with the 10-year EPTA data and Hellings-Downs correlations. We assume Normal priors for (a•, b•, σ*) as well as the value of ϵ0 = 0.38 from ref. 16. The priors for (a•, b•) for black hole masses are inferred from kinematic observations of local black holes in galaxy catalogues66. The prior for σ* is based on the velocity dispersion measured in the Sloan Digital Sky Survey67. We adopt a Normal prior on ρBH from ref. 18 because a prior from ref. 16 yields the strain amplitude in significant tension with the EPTA observations. The posterior and the prior are shown in Fig. 3. All parameters are degenerate with each other due to representing a single observed quantity, the strain amplitude. Nevertheless, the figure illustrates which properties of the supermassive black hole population are constrained the most by the pulsar timing data, given astrophysical uncertainties. Provided the scaling in Equation (2), (0.6%, 5.7%, 0.9%) uncertainties on (a•, b•, σ*), which contribute to M*, and a larger 44.4% uncertainty on ρBH, the posterior informs solely on the number density ρBH. Our improved model influences the posterior to a lesser extent compared to Fig. 1 because the spectral index is fixed at γ = 13/3. Accordingly, it can be noticed in Fig. 1 that lgA does not change much for our improved noise model at the slices of fixed γ = 13/3.

Fig. 3: Posterior of the population parameters for supermassive black hole binaries with the 10-year data.
figure 3

Effectively, the constraints are provided solely by the inferred characteristic strain amplitude of the gravitational wave background, assuming the strain power law index of  − 2/3. Hellings-Downs correlations are included in the model.

Full size image

Discussion

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 (lgA, γ) spans 0.44 for the 10-year data and 0.39 for the 25-year data before our improved noise model. As shown in Fig. 1c, our improved analysis yields 0.95 for the 10-year data and 0.42 for the 25-year data. A larger reduction in the measurement uncertainty when adding ten extra years of data is in agreement with our expectations. A decrease in 1σ uncertainty levels in Fig. 1c compared to Fig. 1b indicates that inter-pulsar correlations in the 10-year data provide additional constraints on the background amplitude and spectral index. A larger maximum-aposteriori lgA in Fig.1b compared to both Fig. 1a and Fig. 1c may also be driven by inter-pulsar correlations.

The shift of the posteriors in Fig. 1 towards larger spectral indices and smaller amplitudes with our improved analysis suggests that the louder pulsar-intrinsic noise with flatter spectra leaks into our measurement of the background strain spectrum when the four new sources of noise we point out are not modelled. The covariance between lgA and γ in Fig. 1 is along the line of equal noise power. Furthermore, in the early part of the 25-year data, our improved model may ameliorate instrumental noise and a lack of frequency coverage in addition to better modelling of millisecond pulsar spin noise33.

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 (lgA, γ) in the 25-year data differ from those of the 10-year data by (0.29, 0.76). It is visible in red contours across all three panels in Fig. 1. A smaller difference between the best-fit (lgA, γ) in the 25-year data and in the 10-year data of (0.16, 0.76) is achieved with our improved analysis. When not modelling Hellings-Downs correlations (Fig. 1c), the difference between the best-fit (lgA, γ) in the 25-year data and in the 10-year data is only (0.01, 0.06). Because the best-fit (lgA, γ) obtained with only temporal correlations is still expected to match those obtained with temporal and Hellings-Downs correlations, it is also possible that we have not removed all noise model misspecification from the analysis of the EPTA data.

Let us briefly hypothesise about the nature of any other mismodeled noise. 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 noise5, as in the EPTA analysis23. Therefore, one potential source of a systematic error could be the mismodelling of the pulsar-specific noise that depends on the radio frequency. A very nearby binary is another example68,69,70,71. Frequency-wise comparison of the inferred strain spectrum against black hole population synthesis models performed earlier by the EPTA (Fig. 3 in ref. 13) suggests that the deviation from γ = 13/3 may occur due to excess noise in two frequency bins, 1.3 × 10−8 Hz and 2.9 × 10−8 Hz. The rest of the spectrum appears 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.

Finally, we calculate evidence for Hellings-Downs correlations in the EPTA data with the revised noise model. Namely, the Bayes factor \({{\mathcal{B}}}\) in favour of the hypothesis that all pulsars contain Hellings-Downs correlated signal with (A, γ) against the hypothesis of the same signal without Hellings-Downs correlations. For 10-year EPTA data, we find \({{\mathcal{B}}}=38\) (\(\ln {{\mathcal{B}}}=3.63\)). For the original EPTA noise model, we find a higher value of \({{\mathcal{B}}}=59\), consistent with the previously reported value of 60 in Table 5 in ref. 6. Therefore, our improved model has slightly decreased evidence for Hellings-Downs correlations. It is a combination of the increased measurement uncertainty and, potentially, a fraction of the previously-reported evidence due to the propagation of unmodelled noise into evidence. For 25-year EPTA data, we find \({{\mathcal{B}}}=3\) (\(\ln {{\mathcal{B}}}=1.17\)). For the original EPTA noise model, we find a consistent value (\({{\mathcal{B}}}=3\), \(\ln {{\mathcal{B}}}=0.93\)), which is similar to the previously reported value \({{\mathcal{B}}}=4\) from Table 5 in ref. 6.

To evaluate the broad form of inter-pulsar correlations in the 10-year EPTA data, we present Fig. 4. It shows the posterior on inter-pulsar correlations, as a function of angular separation between EPTA pulsar pairs, with our revised noise model. The result shows broad consistency with Hellings-Downs correlations. With this, we should point out a caveat that some of the parameter space in correlation coefficients may be ruled out solely because they do not result in a positive definite covariance matrix of the likelihood. Thus, strictly speaking, flexible correlation models are ill-posed (Rutger van Haasteren, private communication).

Fig. 4: Posterior on inter-pulsar correlations of the common stochastic process in the 10-year data.
figure 4

Here, we assume our improved noise model. Correlations are shown as a function of angular separation between the observed pulsar pairs. Hellings-downs function of the gravitational wave background is shown as the dashed line.

Full size image

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 modelled as a multivariate Gaussian distribution \({{\mathcal{L}}}({{\boldsymbol{\delta t}}}| {{\boldsymbol{\theta }}})\), where θ is a vector of parameters of models that describe the data. From Bayes’ theorem, it follows that the posterior distribution of model parameters is \({{\mathcal{P}}}({{\boldsymbol{\theta }}}| {{\boldsymbol{\delta t}}})={{{\mathcal{Z}}}}^{-1}{{\mathcal{L}}}({{\boldsymbol{\delta t}}}| {{\boldsymbol{\theta }}})\pi ({{\boldsymbol{\theta }}})\), where \({{\mathcal{Z}}}\) is Bayesian evidence, a fully-marginalised likelihood. The term π(θ) is called a prior probability distribution, a model of how likely it is to find a certain value of θ in Nature. Model selection is performed by computing the ratio of \({{\mathcal{Z}}}\) for pairs of models, which 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.

The time-domain likelihood of the observed radio pulsar pulse delays-advances time series δt, at the position and of the Solar System barycenter, is a multivariate Gaussian,

$${{\mathcal{L}}}({{\boldsymbol{\delta t}}}| {{\boldsymbol{\theta }}})=\frac{\exp \left(-\frac{1}{2}{({{\boldsymbol{\delta t}}}-{{\boldsymbol{\mu }}})}^{{{\rm{T}}}}{{{\boldsymbol{C}}}}^{-1}({{\boldsymbol{\delta t}}}-{{\boldsymbol{\mu }}})\right)}{\sqrt{\det \left(2\pi {{\boldsymbol{C}}}\right)}}.$$
(3)

Here, θ are model parameters. Vector μ, a function of θ, represents the model prediction for timing residuals. The covariance matrix is C = N + TTB−1T. The likelihood is marginalised over coefficients b which determine the time series realisation Tb of the corresponding signal or noise. Square matrix N only contains diagonal elements that represent temporally uncorrelated noise, which is referred to as white. The so-called design matrix T = [M, F, U] is made up of three blocks: for the timing model, M, for the Fourier series of time-correlated signals, F, and for the epoch-correlated noise, U. Design matrix maps the model parameter space (columns) to the pulse time of arrival space (rows), such that Tb is the time series, like μ. Precisely, parameters ϵ for M are those of the EPTA timing model72, including pulsar spin frequency, its derivative, pulsar sky position, etc. Parameters for F are power spectral density amplitudes a corresponding to Fourier sine and cosine terms at Fourier frequencies. We model temporal correlations in 30 frequencies for a common time-correlated stochastic signal attributed to the gravitational wave background; the number of frequencies for time-correlated noise is determined based on single-pulsar noise analysis23. Parameters for U are delay time series for every epoch, j; epochs are groups of nearly simultaneous pulse arrival times at different radio frequencies. Coefficients b = [ϵ, a, j], corresponding to the same signals as in T, represent signal amplitudes before a transformation to the time domain. Finally, the covariance matrix of coefficients, B, is such that the coefficients b can be generated from the zero-mean Gaussian distribution described by the covariance matrix B. It is the matrix with diagonal blocks \([{{\boldsymbol{\xi }}},{{\boldsymbol{\varphi }}},{{\boldsymbol{{{\mathcal{J}}}}}}]\). The first component ξ, a diagonal matrix of the values of 1040, corresponds to a wide, uninformative prior on the timing model coefficients ϵ. Component \({{\boldsymbol{{{\mathcal{J}}}}}}\) is the covariance matrix of white epoch-correlated noise. The second component, φ, is such that

$${\varphi }_{(ai),(bj)}={P}_{ai}{\delta }_{ab}{\delta }_{ij}+{\Gamma }_{ab}{P}_{i}{\delta }_{ij},$$
(4)

where (a, b) are pulsar indices, (i, j) are frequency indices, Pai is the power spectral density of noise, Pi is the power spectral density of a common time-correlated signal, and Γab determines the degree of spatial correlations between pulsars a and b. In the case of an isotropic gravitational wave background from black hole binaries in circular orbits, Γab is given by the Hellings-Downs function9:

$${\Gamma }_{ab}{| }_{a\ne b}=\frac{1}{2}-\frac{{x}_{ab}}{4}+\frac{3}{2}{x}_{ab}\ln {x}_{ab},$$
(5)

where ζab is the sky separation angle for a given pair of pulsars and \({x}_{ab}=(1-\cos {\zeta }_{ab})/2\). In addition, Γaa = 1.

Time-correlated noise, which yields stochastic timing delays, can be categorised into achromatic and chromatic. Achromatic noise does not depend on the pulsar’s pulse radio frequency. It is also referred to as spin noise because it represents pulsar rotational irregularities. Chromatic noise73 is characterised by an amplitude that depends on the radio frequency at which it is observed. Most EPTA pulsars have measurable levels of dispersion measure (DM) noise, and many pulsars are shown to have spin noise23. Epoch-correlated noise, associated with pulse jitter (short-term pulse profile variations), can influence the inferred strain spectrum of the gravitational wave background at the highest frequencies because of the “white” nature of this noise. Sudden dips with exponential relaxations in timing delays, associated with pulse shape changes on timescales of days-months32, can also influence parameter estimation for time-correlated signals if not modelled. Noise terms missed in the model, despite evidence, are expected, in many cases, to incorrectly increase the inferred signal amplitude. In particular, the gravitational wave model would attempt to absorb excess noise power.

The priors are our models of how parameters governing pulsar-specific noise are distributed across the pulsars. PTA data provides information about the distribution of θ in Nature, so failing to model this distribution could lead to prior misspecification. To address this, we parametrise prior distributions for amplitudes and spectral indices of pulsar spin and DM noise spectra. Our improved analysis introduces hyperparameters Λ to parametrise priors: π(θ∣Λ)π(Λ). We then perform a numerical marginalisation over Λ. We use a new procedure of marginalisation over hyperparameters, which is described in Section 2.2.2 of the companion paper24.

Because the gravitational wave background is a stochastic time-correlated signal, it is important to hierarchically model other stochastic time-correlated signals in the data. Pulsar spin noise is described by amplitudes and spectral indices (ASN, γSN), which are different for every pulsar. Similarly, DM noise is characterised by (ADM, γDM), where the amplitude is referenced to 1400 MHz. Hierarchical noise model directly impacts posteriors on (ASN, γSN, ADM, γDM). Because the data is only compatible with a certain total amount of timing fluctuations, hierarchical noise inference also impacts posteriors on the gravitational wave background’s amplitude.