Ensemble noise properties of the European Pulsar Timing Array
Boris Goncharov Shubhit Sardana


firstpage–lastpage
Wind turbine condition monitoring based on intra- and inter-farm federated learning
    [
    Accepted XXX. Received YYY; in original form ZZZ
===================================================================================




§ ABSTRACT 
The null hypothesis in Pulsar Timing Array (PTA) analyses includes assumptions about ensemble properties of pulsar time-correlated noise. 
These properties are encoded in prior probabilities for the amplitude and the spectral index of the power-law power spectral density of temporal correlations of the noise.

In this work, we introduce a new procedure for numerical marginalisation over the uncertainties in pulsar noise priors. 
The procedure may be used in searches for nanohertz gravitational waves and other PTA analyses to resolve prior misspecification at negligible computational cost. 
Furthermore, we infer the distribution of amplitudes and spectral indices of the power spectral density of spin noise and dispersion measure variation noise based on the observation of 25 millisecond pulsars by the European Pulsar Timing Array (EPTA). 
Our results may be used for the simulation of realistic noise in PTAs. 






pulsar timing arrays – pulsars – gravitational waves








§  INTRODUCTION

Pulsar Timing Arrays <cit.> are experiments that monitor pulse arrival times from galactic millisecond pulsars with a primary goal of detecting nanohertz-frequency gravitational waves <cit.>. 
The most promising source of such gravitational waves is the stochastic superposition of inspiralling supermassive binary black holes in the nearby universe <cit.>.
Thanks to spacetime metric perturbations from gravitational waves, pulse arrival times experience delays and advances (henceforth, delays). 
The power spectral density (PSD) of delays P(f) induced by stochastic gravitational waves manifests temporal correlations and corresponds to the background's characteristic strain spectrum h_c(f). 
The Fourier frequency of timing delays f is exactly the gravitational wave frequency.
For the isotropic stochastic gravitational wave background from circular binaries where the inspiral is driven by gravitational wave emission alone, h_c∝ f^-2/3 corresponding to P(f) ∝ f^-13/3, highlighting signal prominence towards lower frequencies or, equivalently, longer observational time scales <cit.>.
At low frequencies, PTAs are limited by time-correlated “red” noise which is also modelled using a power law. 
The two most ubiquitous sources of red noise in PTA data are the dispersion measure (DM) variation noise <cit.> and the spin noise <cit.>.
DM noise features an additional dependence P(f) ∝ν^-2, where ν is a radio frequency.
It belongs to a broader class of noise processes called “chromatic”, referring to the dependence of the signal amplitude on radio frequency. 
Achromatic red noise that is independent of ν is also called spin noise because it is associated with irregularities in pulsar rotation.
In data of a single pulsar, spin noise acts as a gravitational wave background with a broader range of possible spectral indices. 
A decisive contribution of the stochastic background to PTA data is ultimately determined through the covariance of the signal across pulsar pairs.
The covariance follows the <cit.> function of pulsar angular separation. 

Contemporary PTA analyses are performed using Bayesian inference, where a uniform prior probability is assumed for pulsar-specific red noise parameters. 
Namely, the log-10 amplitude, A, and the spectral index, γ, of noise PSD. 
Priors are the models of how likely a parameter with a certain value is to be found in a given data realization. 

A prior choice is sufficient as long as it is consistent with an observation. 
Uniform priors on PTA noise parameters are designed to be sufficient for single-pulsar noise analyses. 
However, all pulsars represent different realizations of data with respect to pulsar noise parameters[While remaining a single realization of data with respect to parameters governing h_c(f) of gravitational wave background.].  
So, for full-PTA analyses with multiple pulsars, a mismatch between our prior on ( A,γ) and the observed distribution of these parameters may accumulate across pulsars and lead to biased inference. 
This case is known as prior misspecification. 

Implicitly, noise prior misspecification is shown by <cit.> to bias measurements of the strain amplitude of the gravitational wave background in PTAs. 
The authors have mitigated the bias by allowing the data to choose whether red noise is present or absent in a pulsar. 
Furthermore, in simulations of <cit.>, pulsar spin noise with a wide range of ( A,γ) modelled with the standard uniform priors, not representative of ensemble noise properties, yields evidence for a stochastic signal with the same power-law PSD across pulsars (which is not present in the simulated data). 
Such a signal, under the assumption of uniform noise priors, has been reported in real PTA data as a possible precursor to Hellings-Downs correlations of the gravitational wave background <cit.>.

<cit.> regularised incorrect pulsar noise priors by allowing noise amplitudes of the common-spectrum signal of the putative gravitational wave background to vary across pulsars.
By showing that this variance is consistent with zero, the authors confirmed the consistency of the signal with the gravitational wave background.  







Although the model of <cit.> provides the working solution and has clear use cases, the ultimate solution to the misspecification of noise priors is to parametrise priors for all relevant noise parameters.
Both the statement of the problem and the most general solution are clearly outlined in <cit.>. 
The distribution of pulsar noise parameters is inferred simultaneously with a search for Hellings-Downs correlations and other red processes with pulsar-to-pulsar covariance. 
The technique is shown to mitigate systematic errors introduced by incorrect noise priors. 

In this study, we introduce two new methods for modelling ensemble pulsar noise properties in PTA data. 
However, unlike <cit.>, (1) one of our methods can only be used to infer ensemble pulsar noise properties and (2) another method can only be used to remove a systematic error from incorrect pulsar noise priors without gaining access to ensemble noise properties. 
Next, we perform inference of ensemble noise properties of pulsars from the Second Data Release <cit.> of the European Pulsar Timing Array <cit.>.








The rest of the paper is organized as follows. 
In Section <ref>, we outline the data analysis methodology and an overview of sources of noise in PTAs. 
Our two new methods of modelling ensemble pulsar noise properties are presented in Sections <ref> and <ref>.
In Section <ref>, we report on our main results, where each subsection corresponds to a different model of a distribution of noise parameters ( A,γ). 
In Section <ref>, we test the robustness of our models to circular analysis. 
Finally, in Section <ref>, we draw conclusions. 



§ METHODOLOGY
 






In this study, we analyze the second data release (DR2, <cit.>) of the EPTA. 
In particular, we focus our attention on the “DR2full”, which we will thus refer to as EPTA DR2 or data unless specified otherwise. 
The data is based on pulse arrival times from a set of 25 millisecond radio pulsars observed over time spans varying between 14 and 25 years. 
The data also includes pulsar timing models obtained with the least-squared fitting of pulse arrival times to the models <cit.>. 
Timing models describe how pulse arrival times are affected by deterministic properties of individual pulsars such as pulsar spin frequency and derivatives, pulsar position and proper motion in the sky, dispersion measure (DM) and derivatives, and binary orbital parameters (if applicable). 
Contributions of other signals and noise processes to data are referred to as “residuals”, implying that they yield a difference between pulse arrival times predicted by the timing model and measured arrival times. 
Pulse arrival times are referenced to the position of the Solar System barycenter which is defined based on the DE440 ephemeris <cit.>.



 §.§  Standard PTA data analysis methodology


Contributions to the PTA data beyond the timing model are often determined using Bayesian inference as described below. 
The likelihood of data δ t (a vector of the measured pulse times of arrival, ToA) is a Gaussian distribution which is multivariate with respect to a number of observations, 

    ℒ(δ t | θ) = exp(  -1/2 (δ t - μ)^T C^-1(δ t-μ) )/√(( 2πC)),

where θ is a vector of parameters of models that describe the data.
In other words, ℒ(δ t| θ) is the time-domain likelihood. 
The model prediction for pulse arrival times as a function of time is μ(θ), and C(θ) is a covariance matrix that describes stochastic processes. 
Diagonal elements of C, σ^2, correspond to temporally-uncorrelated “white” noise. 
It is modelled as σ^2(e_f,e_q) = σ_ToA^2 e_f^2 + e_q^2, where σ_ToA is the measurement uncertainty on the data provided by the initial timing model fit and (e_f,e_q) are white noise model parameters[Also referred to as the error factor (EFAC) and the error in quadrature (EQUAD), respectively.] which are chosen to be separate for each telescope backend-receiver combination.
Off-diagonal elements of C describe temporal correlations.
For illustrative purposes, it is convenient to represent temporal correlations as components of μ together with the timing model contributions using a reduced-rank approximation <cit.>:

    μ = Fa + Mϵ + Uj + d(t).

Here, F are the Fourier sine and cosine basis function and a are Fourier amplitudes of red processes.
Timing model contributions are similarly modelled via the design matrix M and the coefficients ϵ represent timing model parameters.
Whereas a type of white noise corresponding to timing delays j in units [s] that are the same for all ToAs obtained in the same observing epoch (“jitter” noise) are modelled using the basis U. 
Vector d(t) corresponds to other deterministic signals as a function of pulsar observation time. 
Let us denote b≡ (a, ϵ, j) ∈θ.
Bayes theorem is used to obtain a posterior distribution of model parameters 𝒫(θ | δ t) from the likelihood and the prior π(θ):

    𝒫(θ | δ t) = ℒ(δ t | θ) π(θ)/𝒵,

where 𝒵 is the integral of the numerator over θ, it is termed Bayesian evidence. 

Hierarchical inference – the approach of parametrising prior distributions – is already a part of the standard PTA data analysis machinery. 
It is employed to reduce the parameter space and to make it physically-motivated. 
Instead of Fourier amplitudes a (two per frequency, one per sine term of F and one per cosine term) it is more convenient to measure power-law parameters (A,γ) of the power spectral density of red processes:

    P(f|A,γ) =  A /12 π ^2(f/f_ yr)^-γ.

Thus, the prior on a joins single pulsar likelihoods into a joint posterior:

    𝒫(θ^', θ^” | δ t) = 𝒵^-1∫ℒ(δ t | b,θ^') π(b|θ^”) π(θ^”) π(θ^') db,

where θ^' are parameters which are not included in the reduced-rank approximation in Equation <ref>, such that θ=(θ^',b).
Values θ^” are referred to as hyperparameters that include (A,γ). 
The integral in Equation <ref> is carried out analytically <cit.>. 
The integral represents marginalisation over b, and it is ultimately performed with (F,M,U) and b modelled in C.
For more details, please refer to <cit.>.

Full PTA analysis with stochastic processes that are correlated between pulsars, such as the gravitational wave background, is performed as follows.
To model the PSD of each pulsar-specific noise term, which may vary from pulsar to pulsar, we have a vector of power-law parameters (A,γ), with a pair of noise (hyper-)parameters (A,γ) per pulsar.
We also have a pair of (A_c,γ_c) to model PSD of each “common” stochastic signal, i.e., applicable to all pulsars.
Using notations from Equation <ref>, an “a” component of the term π(b|θ”) for arbitrary pulsars (a,b) and frequencies (i,j) becomes

    π_(a,i),(b,j)(a|A_a,γ_a,A_c,γ_c) =  P_ai(f_i|A_a,γ_a) δ_abδ_ij + 
     + Γ_ab P_i(f_i|A_c,γ_c) δ_ij,

where Γ_ab is the overlap reduction function such as the Hellings-Downs function.
It encodes pulsar-to-pulsar correlations.







 §.§  Hierarchical inference of ensemble noise properties


As part of the standard PTA analysis routine described in Section <ref>, priors π( A,γ) are assumed to be uniform distributions, 𝒰.
Typically, π( A)=𝒰(-20,-12), π(γ)=𝒰(0,7).
As we pointed out in the Introduction, observational data hints that this model may be incorrect. 
The solution is to propose alternative parametrised models π(θ|Λ) which now depend on hyperparameters Λ.
Therefore, generalising ( A, γ) ∈θ, the posterior becomes

    𝒫(θ,Λ|δ t) = 𝒵^-1ℒ(δ t|θ) π(θ|Λ) π(Λ).

The approach of <cit.> is to directly evaluate 𝒫(θ,Λ|δ t).
In this work, we propose computationally efficient marginalisation of this posterior over (a) θ, or (b) Λ, as described below. 
A summary of the methods is provided in Table <ref>.



  §.§.§  Marginalised posterior of hyperparameters


To obtain 𝒫(Λ) marginalised over θ, we propose a solution based on prior reweighting, the application of importance sampling to priors. 
The calculation is done in two steps, both of which manifest one global fit to data. 
First, one obtains 𝒫(θ',θ”|δ t) for every pulsar assuming a fixed prior on noise parameters π(θ^k_i|∅). 
Here, ∅ means that hyperparameters are fixed as in the standard EPTA analysis <cit.>. 
Second, one uses the resulting posterior samples and evidence values Z_∅, i to construct a likelihood marginalised over all parameters except Λ, ℒ(δ t | Λ).
It is shown in <cit.> that the likelihood from Equation <ref> can be written in the following form:

    ℒ(δ t | Λ) = ∏_i^N_psrs Z_∅, i(δ t_i)/n_i∑_k^n_iπ(θ^k_i|Λ)/π(θ^k_i|∅).

In the above equation, N_psrs is the number of pulsars in a PTA, n_i is the number of posterior samples obtained for a given pulsar, θ^k_i is the k'th posterior sample for i'th pulsar. 
Ensemble noise properties of the EPTA presented in this work are obtained based on Equation <ref>. 
The limitation of this approach is that it is based on expressing the total PTA likelihood as the product of single-pulsar likelihoods. 
Therefore, the approach is blind to Hellings-Downs correlations or other inter-pulsar correlations in the timing data. 
Whereas the presence of temporal correlations with the same ( A, γ) in all pulsars has to be imposed via a prior in individual pulsars. 
To model a common-spectrum process <cit.> associated with a gravitational wave background <cit.> in our hierarchical likelihood defined by Equation <ref>, we impose an additional spin-noise-like term with γ=13/3 in individual pulsar likelihoods. 
Improved methods for separating the effect of gravitational wave background from pulsar-intrinsic work may be explored in future work. 





  §.§.§  Marginalisation over hyperparameters


The posterior from Equation <ref> marginalised over Λ can be written in the following form:

    𝒫(θ|δ t) =  ℒ(δ t|θ) π(θ|∅)/𝒵×1/n_p∑_k^n_pπ(θ|Λ_k)/π(θ|∅),

where Λ_k are n_p samples from the prior. 
The reader may notice that the resulting posterior in Equation <ref> is the product (×) of the standard PTA posterior and the weight factor. 
A derivation of Equation <ref> is provided in the Appendix. 
Unlike Equation <ref>, Equation  <ref> allows to simultaneously model Hellings-Downs and other inter-pulsar correlations in the data while accounting for the uncertainty in pulsar noise priors. 
It is used more extensively in the companion paper (Goncharov et al., in prep.), where the likelihood is multiplied by the weight factor from Equation <ref>, and the rest of the analysis is performed in a standard way. 
Alternatively, the posterior 𝒫(θ|δ t) can be reweighted into a Λ-marginalised posterior using rejection sampling.



§ RESULTS


In this work, we infer distributions for parameters governing pulsar spin noise, π( A_SN,γ_SN), and DM noise, π( A_DM,γ_DM), for EPTA DR2. 
EPTA DR2 contains a common-spectrum stochastic process associated with the gravitational wave background and γ consistent with 13/3 (Goncharov et al., in prep.). 
If we do not account for this noise term, it may appear as a “quasi-common” red noise <cit.>. 
In other words, as a separate cluster in the distribution of spin noise parameters. 
Because we are interested in knowing the distribution of pulsar-intrinsic parameters, we include an additional red noise term with γ_c=13/3 to marginalise over a contribution of the common signal.
This is done during the analysis of single pulsar noise corresponding to the first step in Section <ref>. 
It corresponds to modelling common signal in our data via a prior as per Table <ref>. 
Because we are agnostic about A_c, the approach imposes a suboptimal uncertainty for measuring ( A, γ) in pulsars where γ is consistent with 13/3 ≈ 4. 
Furthermore, to separate contributions of other noise processes such as chromatic noise (including scattering variations) and band- or system-dependent noise, we adopt single-pulsar noise models and optimal numbers of Fourier frequency bins from the EPTA DR2 noise analysis <cit.>.




 §.§  Uniform distribution of noise parameters


We start with the case of uniform distributions for noise parameters and we measure hyperparameters that govern uniform prior range.
For DM noise, these are min(A_DM), max(A_DM), min(γ_DM), max(γ_DM). 
We apply the same principle to spin noise (SN) parameters. 
Pulsar-specific measurements of (A,γ) and the inferred boundaries of the uniform distribution for spin noise and DM noise parameters – with measurement uncertainties – are shown in Figure <ref>. 
The first important observation is that the data rules out, with high credibility, the standard values of uniform prior boundaries. 
The data suggests that pulsar noise parameters are distributed in a more narrow range.
Posterior density 𝒫( A, γ) from pulsars that do not have evidence for the respective noise term according to the <cit.> is shown in colour.
The data from these pulsars has also contributed to the measurement of hyperparameters. 
The inferred noise distribution model extends as little as possible to include as many identified noise parameters in pulsars as possible. 

One may, in principle, extend our model to include selection effects associated with a lack of sensitivity for certain A,γ <cit.>. 
This would enable excess density of π( A|Λ) towards low A.





 §.§  Normal distribution of noise parameters


The uniform prior model may remain a good approximation, but it is not the most natural choice. 
Therefore, we further discuss the possibility that pulsar noise parameters are normally distributed. 
Thus, for spin noise, π( A_SN,γ_SN | Λ_𝒩), and Λ_𝒩=(μ_ A_SN, σ_ A_SN, μ_γ_SN, σ_γ_SN, ρ_SN).
The same model is applied to DM noise. 
Here, μ and σ correspond to the mean and the standard deviation of the normal distribution of either A or γ, as subscripts indicate.
Parameter ρ∈ [-1,1] is the correlation coefficient between A and γ. 
It is important to note that we obtained our posterior samples from individual pulsars on step one from Section <ref> based on uniform priors.
So, for step two, there are no posterior samples to recycle from outside of these boundaries. 
Therefore, we truncate our normal distribution model and apply the same boundaries. 
The inferred boundaries of the truncated normal distribution for spin noise and DM noise parameters are shown in Figure <ref>. 
The shaded area corresponds to an uncertainty in (σ,ρ) and arrows separately show the uncertainty in μ. 
The uncertainty in ρ is seen for DM noise as the difference between the tilt of the outer border of the shaded area compared to the inner border. 
Whereas for spin noise the value of ρ≈ -1 is strongly preferred. 





 §.§  Distribution of noise parameters as a mixture model


One may also be interested whether the distribution of pulsar noise parameters is more complex. 
For example, in Figure 1 in <cit.>, one may notice the clustering of noise parameters in two areas.
First, around γ≈ 1.
Second, a cluster of pulsars around γ=13/3 of the gravitational wave background. 
A similar clustering is observed in the PTA noise analysis by the <cit.>, where it is explained that pulsars with γ≈ 13/3 contribute to the common-spectrum process in full-PTA analysis. 
An analogous Figure is shown in <cit.>, although the presence of two over-densities is not obvious there.
Although we have taken steps to marginalize the contribution of the common-spectrum process and uncover properties of pulsar-intrinsic noise, it is nevertheless useful to consider the possibility of a mixture model of a Gaussian distribution 𝒩 and a Uniform distribution 𝒰, where a broader Uniform distribution may fit potential outliers. 
The density of such a distribution is then ν𝒩 + (1-ν)𝒰, where ν∈ [0,1] is the contribution of the normal distribution.



 §.§  Summary of the results


As part of our measurements of hyperparameters, we have also obtained Bayesian evidence values 𝒵.
With this, we perform model selection to determine which of the models fit the observed distribution of pulsar noise parameters best. 
First, we find that the standard `static' uniform priors are disfavored compared to hierarchical uniform priors with natural log Bayes factors of ≳ 100. 
Therefore, the use of standard uniform priors is not recommended for PTA data analysis. 
The rest of the results of our hyperparameter estimation and hierarchical model selection are shown in Table <ref>. 

Some ensemble properties of spin noise (SN) are similar to DM noise, whereas some properties are different. 
For spin noise, the normal distribution fits the data as well as the uniform distribution. 
However, the covariance between A_SN and γ_SN plays an important role. 
In Figures <ref> and <ref>, one may notice a diagonal trend which results in finding a covariance parameter ρ_SN≈-1. 
The model with a lack of covariance, ρ_SN=0 is disfavored with lnℬ=6. 
The covariance between A_SN and γ_SN can be explained by noticing that it follows a line of equal noise power. 
In contrast, for DM noise, we do not find strong evidence for covariance between A_DM and γ_DM. 
With a lack of covariance, the uniform model is strongly preferred over the Gaussian model with a log Bayes factor of 6. 
The mixture model is disfavored by our observations of both spin noise and DM noise, indicating the lack of outlying noise terms. 





























































































§ CAVEATS OF THE CIRCULAR ANALYSIS
 

Circular analysis or “double dipping” refers to the usage of data to inform on the model which is then applied to analyse the same data. 
Such a practice may lead to the underestimation of measurement uncertainty and systematic errors. 
This point is stressed by van Haasteren (2024, in prep.). 
The author argues against the approach of determining which pulsar spin noise terms to include in the full-PTA analysis based on (a lack of) evidence for these terms in single-pulsar noise model selection. 
Indeed, double dipping is not a proper statistical treatment, it should be abandoned in favour of the proper model selection and model averaging performed in one global fit, as recommended by the author.
However, to give a fair overview of the problem, we also show that circular analysis does not always yield incorrect results in analyses of PTA data. 

Before we discuss the application of hierarchical inference in the subsections below, we would like to point out that it is not clear from van Haasteren (2024, in prep.) that the aforementioned approach of selecting spin noise terms for subsequent full-PTA analyses has led to underestimation of measurement uncertainties reported by PTAs. 
For model averaging, the prior odds between the existence or absence of common red noise are proposed to be 50%-50% in van Haasteren (2024, in prep.).
For the case of including a red noise term to every pulsar, the prior odds will be arbitrary and vary from pulsar to pulsar which is clear from the equations in the Appendix in van Haasteren (2024, in prep.).
However, data informs on these odds, so they should become a model (hyper-)parameter to avoid prior misspecification, the same problem discussed in this work.
These parametrised odds will act as model selection. 
In the limit where there is no measurable red noise in pulsars and the data informs on it sufficiently well, this will lead to the elimination of the contribution of one of the models. 
Because different noise models lead to different measurement uncertainties for ( A_CP,γ_CP), it is possible that the uncertainty introduced by enforcing red noise to every pulsar is suboptimal.













 §.§  A comparison between circular analysis, incorrect priors, and a proper analysis


One example of circular analysis in the context of our hierarchical inference is finding best-fit Λ and using it for the gravitational wave search. 
It also leads to resolving prior misspecification but potentially at the cost of a reduced measurement uncertainty for ( A, γ) of the gravitational wave background. 
The reduction may come thanks to simply ignoring a part of the intrinsic measurement uncertainty for Λ, so a regular noise fluctuation is more likely to render our estimate of ( A, γ) to be inconsistent with the true value. 
In Figure <ref>, we show three measurements of ( A, γ) of the common-spectrum process in the EPTA data. 
Blue dashed contours correspond to a posterior obtained with the standard uniform spin noise priors, π(θ|Λ)=𝒰(θ|∅). 
Red dotted contours correspond to π(θ|Λ)=𝒩(θ |μ_0,σ_0,ρ=0), where (μ_0,σ_0) are the values obtained for this model from Table <ref>. 
Therefore, the red dotted line corresponds to a circular analysis. 
Filled green contours correspond to a proper (global-fit) measurement obtained based on Equation <ref>, it is discussed in more detail in the companion paper by Goncharov et al. (2024, in prep.). 
The measurement uncertainty of ( A, γ) obtained with circular analysis matches that obtained with a full hierarchical analysis. 
It is even 5-8% larger, suggesting that in the current setting circular analysis overestimates the intrinsic measurement uncertainty. 
Note, however, that the circular analysis shifts the maximum-aposteriori value of ( A, γ), so the use of it is still not recommended.























































 §.§  A toy example: circular analysis of common red noise


Keeping in mind the caveats above, for the test purposes we also perform the estimation of hyperparameters Λ for the mixture model of two Gaussian distributions where hyper-priors on the second component of the mixture correspond to EPTA measurement of the common-spectrum process parameters <cit.>. 
This exercise remains useful to explore the interplay between pulsar-intrinsic noise and the common-spectrum process and to understand the performance of the model. 
For the two-Gaussian mixture model, π( A_SN,γ_SN | Λ_2𝒩), and Λ_2𝒩=(μ_ A_SN^1,2, σ_ A_SN^1,2, μ_γ_SN^1,2, σ_γ_SN^1,2, ρ_SN^1,2, and ν_𝒩).
Indices (1,2) refer to each of the two components in a mixture, while ν corresponds to a fraction of the first mixture component in the total prior probability. 
For the second component to correspond to the common-spectrum process in the EPTA data, we impose a Gaussian hyperprior determined by the measurement uncertainty on ( A,γ) reported by the EPTA for our data.

We find that the fraction of the second normal component is mostly consistent with one and inconsistent with zero. 
However, hyperparameters of the second mixture component – mean and (co-)variance – are well-constrained. 


Interestingly, we also find a marginal excess posterior density at (σ_ A_SN^2,σ_γ_SN^2)=(0,0). 
This result is related to <cit.>'s measurement of σ_ A to be consistent with zero.
When removing the common-spectrum term with γ=13/3 from our pulsar-intrinsic noise model, the posterior density increases, as expected.
This is shown in Figure <ref>.
















§ CONCLUSIONS
 

We performed inference of ensemble properties of spin noise and DM noise in the 25-year version of the second data release of the European Pulsar Timing Array (EPTA). 
Overall, we find that the standard uniform priors used in the previous analysis are too wide and not representative of the observations. 
However, a parametrised uniform prior distribution is a good fit. 
We, therefore, recommend correctly accounting for ensemble pulsar noise properties using Equation <ref> in future full-PTA analyses to avoid systematic errors (parameter estimation and model selection biases). 



§ ACKNOWLEDGEMENTS


We thank Rutger van Haasteren for insightful discussions about hierarchical inference and circular analysis. 
Some of our calculations were carried out using the OzSTAR Australian national facility (high-performance computing) at Swinburne University of Technology. 





§ DATA AVAILABILITY


 
The data used for our study, EPTA DR2 by the <cit.>, is available at zenodozenodo.
The code to reproduce the results of this analysis and the analysis by <cit.> is available in the  branch at https://github.com/bvgoncharov/pta_gwb_priorsgithub.com/bvgoncharov/pta_gwb_priors. 
Other data may be provided by the corresponding author upon request. 
The PTA likelihood is incorporated in enterprise <cit.> and posterior sampling is performed using ptmcmcsampler <cit.> and dynesty <cit.>. 








mnras
 



















§ NUMERICAL MARGINALISATION OVER HYPERPARAMETERS


In this Section, we derive Equation <ref> which describes a hierarchical posterior marginalised over hyperparameters Λ. 
After taking an integral over Λ, Equation <ref> becomes

    𝒫(θ|δ t) =  𝒵^-1ℒ(δ t|θ) ∫π(θ|Λ) π(Λ) dΛ.

Following the formalism of importance sampling which was used to obtain Equation <ref>, we refer to π(θ|Λ) as the target distribution. 
Similarly, we introduce the proposal distribution, the standard PTA noise prior that does not depend on hyperparameters, π(θ|∅). 
Next, we multiply Equation <ref> by unity and expand the unity on the right-hand side as a proposal distribution divided by itself.
Rearranging the multipliers,

    𝒫(θ|δ t) = ℒ(δ t|θ) π(θ|∅)/𝒵∫π(θ|Λ)/π(θ|∅)π(Λ) dΛ.

Next, we use the expression for the expectation value of a probability density f(x) given the known probability density p(x)[Equation <ref> is also used to derive Equation <ref>. There, p(x) is taken to be a posterior distribution instead.]:

    ⟨ f(x) ⟩_p(x) = ∫ f(x) p(x) dx ≈1/n_s∑_i^n_sf(x_i),

where x_i are n_s samples from p(x). 
Taking π(Λ) as p(x), we arrive at Equation <ref>. 












