The null hypothesis in Pulsar Timing Array PTA analyses includes assumptions about ensemble properties of pulsar timecorrelated noise. These properties are encoded in prior probabilities for the amplitude and the spectral index of the powerlaw 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.
['Astrophysics - High Energy Astrophysical Phenomena', 'Astrophysics - Instrumentation and Methods for Astrophysics']
Ensemble noise properties of the European Pulsar Timing Array
{'Ensemble noise properties of the European Pulsar Timing Array': 'Boris Goncharov 1 , 2 ★ , Shubhit Sardana 1 , 3  \n1 Max Planck Institute for Gravitational Physics (Albert Einstein Institute), D-30167 Hannover, Germany  \n- 2 Leibniz Universität Hannover, D-30167 Hannover, Germany  \n3 Department of Physics, IISER Bhopal, Bhauri Bypass Road, Bhopal, 462066, India  \n10 September 2024', 'ABSTRACT': 'The null hypothesis in Pulsar Timing Array (PTA) analyses includes assumptions about ensemble properties of pulsar timecorrelated noise. These properties are encoded in prior probabilities for the amplitude and the spectral index of the powerlaw 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.  \nKey words: pulsar timing arrays - pulsars - gravitational waves', '1 INTRODUCTION': "Pulsar Timing Arrays (PTAs, Foster & Backer 1990) are experiments that monitor pulse arrival times from galactic millisecond pulsars with a primary goal of detecting nanohertz-frequency gravitational waves (Sazhin 1978; Detweiler 1979). The most promising source of such gravitational waves is the stochastic superposition of inspiralling supermassive binary black holes in the nearby universe (Rosado et al. 2015). Thanks to spacetime metric perturbations from gravitational waves, pulse arrival times experience delays and advances (henceforth, delays). The power spectral density (PSD) of delays 𝑃 ( 𝑓 ) induced by stochastic gravitational waves manifests temporal correlations and corresponds to the background's characteristic strain spectrum ℎ c ( 𝑓 ) . The Fourier frequency of timing delays 𝑓 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, ℎ c ∝ 𝑓 -2 / 3 corresponding to 𝑃 ( 𝑓 ) ∝ 𝑓 -13 / 3 , highlighting signal prominence towards lower frequencies or, equivalently, longer observational time scales (Renzini et al. 2022). 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 (Keith et al. 2013) and the spin noise (Shannon & Cordes 2010). DM noise features an additional dependence 𝑃 ( 𝑓 ) ∝ 𝜈 -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.  \nA 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 Hellings & Downs (1983) function of pulsar angular separation.  \nContemporary 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, lg 𝐴 , 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 1 . So, for full-PTA analyses with multiple pulsars, a mismatch between our prior on ( lg 𝐴, 𝛾 ) and the observed distribution of these parameters may accumulate across pulsars and lead to biased inference. This case is known as prior misspecification .  \nImplicitly, noise prior misspecification is shown by Hazboun et al. (2020) 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 Goncharov et al. (2021b); Zic et al. (2022), pulsar spin noise with a wide range of ( lg 𝐴, 𝛾 ) 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  \nwave background (Goncharov et al. 2021b; NANOGrav Collaboration 2020; Chen et al. 2021).  \nGoncharov et al. (2022) 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 Goncharov et al. (2022) 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 van Haasteren (2024). Thedistribution 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.  \nIn this study, we introduce two new methods for modelling ensemble pulsar noise properties in PTA data. However, unlike van Haasteren (2024), (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 (DR2, EPTA Collaboration 2023) of the European Pulsar Timing Array (EPTA, Desvignes et al. 2016; Kramer & Champion 2013).  \nThe rest of the paper is organized as follows. In Section 2, 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 2.2.1 and 2.2.2. In Section 3, we report on our main results, where each subsection corresponds to a different model of a distribution of noise parameters ( lg 𝐴, 𝛾 ) . In Section 4, we test the robustness of our models to circular analysis. Finally, in Section 5, we draw conclusions.", '2 METHODOLOGY': "In this study, we analyze the second data release (DR2, EPTA Collaboration (2023)) 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 (Hobbs et al. 2006). 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 (Park et al. 2021).", '2.1 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 𝜹𝒕 (a vector of the measured pulse times of arrival,  \nToA) is a Gaussian distribution which is multivariate with respect to a number of observations,  \nL( 𝜹𝒕 | 𝜽 ) = exp GLYPH<18> -1 2 ( 𝜹𝒕 -𝝁 ) 𝑇 𝑪 -1 ( 𝜹𝒕 -𝝁 ) GLYPH<19> √︁ det ( 2 𝜋 𝑪 ) , (1)  \nwhere 𝜽 is a vector of parameters of models that describe the data. In other words, L( 𝜹𝒕 | 𝜽 ) is the time-domain likelihood. The model prediction for pulse arrival times as a function of time is 𝝁 ( 𝜽 ) , and 𝑪 ( 𝜽 ) is a covariance matrix that describes stochastic processes. Diagonal elements of 𝑪 , 𝝈 2 , correspond to temporally-uncorrelated 'white' noise. It is modelled as 𝝈 2 ( 𝒆 f , 𝒆 q ) = 𝝈 2 ToA 𝒆 2 f + 𝒆 2 q , where 𝝈 ToA is the measurement uncertainty on the data provided by the initial timing model fit and ( 𝒆 f , 𝒆 q ) are white noise model parameters 2 which are chosen to be separate for each telescope backend-receiver combination. Off-diagonal elements of 𝑪 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 (Lentati et al. 2013; van Haasteren & Vallisneri 2014):  \n𝝁 = 𝑭𝒂 + 𝑴𝝐 + 𝑼 𝒋 + 𝒅 ( 𝑡 ) . (2)  \nHere, 𝑭 are the Fourier sine and cosine basis function and 𝒂 are Fourier amplitudes of red processes. Timing model contributions are similarly modelled via the design matrix 𝑴 and the coefficients 𝝐 represent timing model parameters. Whereas a type of white noise corresponding to timing delays 𝒋 in units [s] that are the same for all ToAs obtained in the same observing epoch ('jitter' noise) are modelled using the basis 𝑼 . Vector 𝒅 ( 𝑡 ) corresponds to other deterministic signals as a function of pulsar observation time. Let us denote 𝒃 ≡ ( 𝒂 , 𝝐 , 𝒋 ) ∈ 𝜽 . Bayes theorem is used to obtain a posterior distribution of model parameters P( 𝜽 | 𝜹𝒕 ) from the likelihood and the prior 𝜋 ( 𝜽 ) :  \nP( 𝜽 | 𝜹𝒕 ) = L( 𝜹𝒕 | 𝜽 ) 𝜋 ( 𝜽 ) Z , (3)  \nwhere Z is the integral of the numerator over 𝜽 , it is termed Bayesian evidence.  \nHierarchical 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 𝒂 (two per frequency, one per sine term of 𝑭 and one per cosine term) it is more convenient to measure power-law parameters ( 𝐴, 𝛾 ) of the power spectral density of red processes:  \n𝑃 ( 𝑓 | 𝐴, 𝛾 ) = 𝐴 12 𝜋 2 GLYPH<18> 𝑓 𝑓 yr GLYPH<19> -𝛾 . (4)  \nThus, the prior on 𝒂 joins single pulsar likelihoods into a joint posterior:  \nP( 𝜽 ' , 𝜽 '' | 𝜹𝒕 ) = Z -1 ∫ L( 𝜹𝒕 | 𝒃 , 𝜽 ' ) 𝜋 ( 𝒃 | 𝜽 '' ) 𝜋 ( 𝜽 '' ) 𝜋 ( 𝜽 ' ) 𝑑 𝒃 , (5)  \nwhere 𝜽 ' are parameters which are not included in the reducedrank approximation in Equation 2, such that 𝜽 = ( 𝜽 ' , 𝒃 ) . Values 𝜽 '' are referred to as hyperparameters that include ( 𝐴, 𝛾 ) . The integral in Equation 5 is carried out analytically (Lentati et al. 2013; van Haasteren & Vallisneri 2014). The integral represents marginalisation over 𝒃 , and it is ultimately performed with ( 𝑭 , 𝑴 , 𝑼 ) and 𝒃  \nmodelled in 𝑪 . For more details, please refer to NANOGrav Collaboration (2016).  \nFull 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 ( 𝑨 , 𝜸 ) , with a pair of noise (hyper-)parameters ( 𝐴, 𝛾 ) per pulsar. We also have a pair of ( 𝐴 c , 𝛾 c ) to model PSD of each 'common' stochastic signal, i.e. , applicable to all pulsars. Using notations from Equation 2, an ' 𝒂 ' component of the term 𝜋 ( 𝒃 | 𝜽 '' ) for arbitrary pulsars ( 𝑎, 𝑏 ) and frequencies ( 𝑖, 𝑗 ) becomes  \n𝜋 ( 𝑎,𝑖 ) , ( 𝑏, 𝑗 ) ( 𝒂 | 𝐴 𝑎 , 𝛾 𝑎 , 𝐴 c , 𝛾 c ) = 𝑃 𝑎𝑖 ( 𝑓 𝑖 | 𝐴 𝑎 , 𝛾 𝑎 ) 𝛿 𝑎𝑏 𝛿 𝑖 𝑗 + + Γ 𝑎𝑏 𝑃 𝑖 ( 𝑓 𝑖 | 𝐴 c , 𝛾 c ) 𝛿 𝑖 𝑗 , (6)  \nwhere Γ 𝑎𝑏 is the overlap reduction function such as the HellingsDowns function. It encodes pulsar-to-pulsar correlations.", '2.2 Hierarchical inference of ensemble noise properties': 'As part of the standard PTA analysis routine described in Section 2.1, priors 𝜋 ( lg 𝐴, 𝛾 ) are assumed to be uniform distributions, U . Typically, 𝜋 ( lg 𝐴 ) = U(20 , -12 ) , 𝜋 ( 𝛾 ) = U( 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 ( lg 𝐴, 𝛾 ) ∈ 𝜽 , the posterior becomes  \nP( 𝜽 , 𝚲 | 𝜹𝒕 ) = Z -1 L( 𝜹𝒕 | 𝜽 ) 𝜋 ( 𝜽 | 𝚲 ) 𝜋 ( 𝚲 ) . (7)  \nThe approach of van Haasteren (2024) is to directly evaluate P( 𝜽 , 𝚲 | 𝜹𝒕 ) . 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 1.', '2.2.1 Marginalised posterior of hyperparameters': "To obtain P( 𝚲 ) 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 P( 𝜽 ' , 𝜽 '' | 𝜹𝒕 ) for every pulsar assuming a fixed prior on noise parameters 𝜋 ( 𝜽 𝑘 𝑖 |∅) . Here, ∅ means that hyperparameters are fixed as in the standard EPTA analysis (EPTA Collaboration & InPTA Collaboration 2023a). Second, one uses the resulting posterior samples and evidence values Z ∅ ,𝑖 to construct a likelihood marginalised over all parameters except 𝚲 , L( 𝜹𝒕 | 𝚲 ) . It is shown in Thrane & Talbot (2019) that the likelihood from Equation 7 can be written in the following form:  \nL( 𝜹𝒕 | 𝚲 ) = 𝑁 psrs GLYPH<214> 𝑖 Z ∅ ,𝑖 ( 𝜹𝒕 𝑖 ) 𝑛 𝑖 𝑛 𝑖 ∑︁ 𝑘 𝜋 ( 𝜽 𝑘 𝑖 | 𝚲 ) 𝜋 ( 𝜽 𝑘 𝑖 |∅) . (8)  \nIn the above equation, 𝑁 psrs is the number of pulsars in a PTA, 𝑛 𝑖 is the number of posterior samples obtained for a given pulsar, 𝜃 𝑘 𝑖 is the 𝑘 'th posterior sample for 𝑖 'th pulsar. Ensemble noise properties of the EPTA presented in this work are obtained based on Equation 8. 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 ( lg 𝐴, 𝛾 ) in all pulsars has to be imposed via a prior in individual pulsars. To model a commonspectrum process (NANOGrav Collaboration 2020; Goncharov et al.  \n2021b; Chen et al. 2021) associated with a gravitational wave background (NANOGrav Collaboration 2023a; EPTA Collaboration & InPTA Collaboration 2023b; Reardon et al. 2023) in our hierarchical likelihood defined by Equation 8, we impose an additional spin-noiselike 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.", '2.2.2 Marginalisation over hyperparameters': 'The posterior from Equation 7 marginalised over 𝚲 can be written in the following form:  \nP( 𝜽 | 𝜹𝒕 ) = L( 𝜹𝒕 | 𝜽 ) 𝜋 ( 𝜽 |∅) Z × 1 𝑛 p 𝑛 p ∑︁ 𝑘 𝜋 ( 𝜽 | 𝚲 𝑘 ) 𝜋 ( 𝜽 |∅) , (9)  \nwhere 𝚲 𝑘 are 𝑛 p samples from the prior. The reader may notice that the resulting posterior in Equation 9 is the product ( × ) of the standard PTA posterior and the weight factor. A derivation of Equation 9 is provided in the Appendix. Unlike Equation 8, Equation 9 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 9, and the rest of the analysis is performed in a standard way. Alternatively, the posterior P( 𝜽 | 𝜹𝒕 ) can be reweighted into a 𝚲 -marginalised posterior using rejection sampling.', '3 RESULTS': "In this work, we infer distributions for parameters governing pulsar spin noise, 𝜋 ( lg 𝐴 SN , 𝛾 SN ) , and DM noise, 𝜋 ( lg 𝐴 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 (as per the terminology of Goncharov et al. 2022). 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 2.2. It corresponds to modelling common signal in our data via a prior as per Table 1. Because we are agnostic about lg 𝐴 c , the approach imposes a suboptimal uncertainty for measuring ( lg 𝐴, 𝛾 ) 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 (EPTA Collaboration & InPTA Collaboration 2023a).", '3.1 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 ( lg 𝐴 DM ) , max ( lg 𝐴 DM ) , min ( 𝛾 DM ) , max ( 𝛾 DM ) . We apply the same principle to spin noise (SN) parameters. Pulsar-specific measurements of ( lg 𝐴, 𝛾 ) and the inferred boundaries of the uniform distribution for spin noise and DM noise  \nTable 1. Hierarchical Bayesian analysis methods for Pulsar Timing Arrays (PTAs). Parameters 𝜽 of a model describing PTA data include parameters of noise and signals such as the gravitational wave background (GWB). Hyperparameters 𝚲 determine the distribution of noise parameters in nature, 𝜋 ( 𝜽 | 𝚲 ) .  \nparameters - with measurement uncertainties - are shown in Figure 1. 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 P( lg 𝐴, 𝛾 ) from pulsars that do not have evidence for the respective noise term according to the EPTA Collaboration & InPTA Collaboration (2023a) 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 𝐴, 𝛾 (see Appendix E in Thrane & Talbot 2019). This would enable excess density of 𝜋 ( lg 𝐴 | 𝚲 ) towards low lg 𝐴 .', '3.2 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, 𝜋 ( lg 𝐴 SN , 𝛾 SN | 𝚲 N ) , and 𝚲 N = ( 𝜇 lg 𝐴 SN , 𝜎 lg 𝐴 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 lg 𝐴 or 𝛾 , as subscripts indicate. Parameter 𝜌 ∈ [-1 , 1 ] is the correlation coefficient between lg 𝐴 and 𝛾 . It is important to note that we obtained our posterior samples from individual pulsars on step one from Section 2.2 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 2. 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.', '3.3 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 Goncharov et al. (2021a), 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 NANOGrav Collaboration (2023b), 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 Reardon et al. (2023), although the presence of two over-densities is not obvious there. Although we have taken steps to marginalize the contribution of the commonspectrum process and uncover properties of pulsar-intrinsic noise, it is nevertheless useful to consider the possibility of a mixture model of a Gaussian distribution N and a Uniform distribution U , where a broader Uniform distribution may fit potential outliers. The density of such a distribution is then 𝜈 N + ( 1 -𝜈 )U , where 𝜈 ∈ [ 0 , 1 ] is the contribution of the normal distribution.', '3.4 Summary of the results': "As part of our measurements of hyperparameters, we have also obtained Bayesian evidence values Z . 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 2.  \nSome 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 lg 𝐴 SN and 𝛾 SN plays an important role. In Figures 1a and 2a, 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 B = 6. The covariance between lg 𝐴 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 lg 𝐴 DM and 𝛾 DM .  \nFigure 2. The inferred form of the normal distribution of pulsar red noise parameters. The colour and the black solid contours represent measurements of noise parameters in individual pulsars, as explained in Figure 1. Our hierarchical analysis yields hyperparameters which determine how pulsar noise parameters are distributed. Shaded ellipses correspond to 1 𝜎 credible levels for hyperparameters ( 𝝈 , 𝜌 ) which govern the width of the distribution and the tilt, respectively. Arrows correspond to 1 𝜎 credible levels for hyperparameters 𝝁 which govern the position of the distribution. Hyperparameter values corresponding to the shaded area and the dashed line represent a row N( 𝜽 | 𝝁 , 𝝈 , 𝜌 ) in Table 2.  \n<!-- image -->  \n(a) Spin noise (SN)  \n<!-- image -->  \n(b) Dispersion measure (DM) variation noise  \nFigure 1. The inferred form of the uniform distribution of pulsar red noise parameters. Black solid contours correspond to 1 𝜎 posterior levels of pulsar noise parameters when analysing single pulsar data, provided that the noise term is resolved in a pulsar based on model selection from EPTA Collaboration & InPTA Collaboration (2023a). Colour corresponds to the posterior density of pulsars where noise terms were not resolved. Our hierarchical analysis yields hyperparameters which determine how pulsar noise parameters are distributed. The shaded rectangular areas correspond to 1 𝜎 credible levels for hyperparameters ( min , max ) that govern limits of the uniform distribution based on our hierarchical analysis. Hyperparameter values corresponding to the shaded area and the dashed line represent a row U( 𝜽 | min , max ) in Table 2.  \n<!-- image -->  \n(b) Dispersion measure (DM) variation noise  \n<!-- image -->  \nWith 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.", '4 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 measure-  \nment 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.  \nBefore we discuss the application of hierarchical inference in the  \nTable 2. Results of hyperparameter estimation and hierarchical model selection for ensemble pulsar noise properties of the European Pulsar Timing Array. DMnoise and spin noise are parametrized by 𝜽 = ( lg 𝐴, 𝛾 ) . We present the results for the following models of pulsar noise parameter distributions: a uniform distribution U , a normal distribution N with and without covariance between lg 𝐴 and 𝛾 , a mixture model of a normal and a uniform distribution. Columns P( 𝚲 | 𝜹𝒕 ) contain maximum-aposteriori hyperparameter values with 1 𝜎 credible levels. Columns ln B contain log Bayes factors in favour of a model of interest against the uniform distribution model. The fact that they are negative means that the uniform model is the best fit.  \nsubsections 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  \nwill lead to the elimination of the contribution of one of the models. Because different noise models lead to different measurement uncertainties for ( lg 𝐴 CP , 𝛾 CP ) , it is possible that the uncertainty introduced by enforcing red noise to every pulsar is suboptimal.", '4.1 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 ( lg 𝐴, 𝛾 ) of the gravitational wave background. The reduction may come thanks to simply ignoring a part of the intrinsic measure-  \nFigure 3. A comparison between circular analysis (red dotted contours), incorrect broad uniform priors (blue dashed contours), and a proper hierarchical analysis (green contours) of the common-spectrum process associated with the gravitational wave background.  \n<!-- image -->  \nrtainty for 𝚲 , so a regular noise fluctuation is more likely to render our estimate of ( lg 𝐴, 𝛾 ) to be inconsistent with the true value. In Figure 3, we show three measurements of ( lg 𝐴, 𝛾 ) of the common-spectrum process in the EPTA data. Blue dashed contours correspond to a posterior obtained with the standard uniform spin noise priors, 𝜋 ( 𝜽 | 𝚲 ) = U( 𝜽 |∅) . Red dotted contours correspond to 𝜋 ( 𝜽 | 𝚲 ) = N( 𝜽 | 𝝁 0 , 𝝈 0 , 𝜌 = 0 ) , where ( 𝝁 0 , 𝝈 0 ) are the values obtained for this model from Table 2. Therefore, the red dotted line corresponds to a circular analysis. Filled green contours correspond to a proper (global-fit) measurement obtained based on Equation 9, it is discussed in more detail in the companion paper by Goncharov et al. (2024, in prep.). The measurement uncertainty of ( lg 𝐴, 𝛾 ) 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 ( lg 𝐴, 𝛾 ) , so the use of it is still not recommended.', '4.2 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 EPTA Collaboration & InPTA Collaboration (2023b). This exercise remains useful to explore the interplay between pulsar-intrinsic noise and the commonspectrum process and to understand the performance of the model. For the two-Gaussian mixture model, 𝜋 ( lg 𝐴 SN , 𝛾 SN | 𝚲 2 N ) , and 𝚲 2 N = ( 𝜇 1 , 2 lg 𝐴 SN , 𝜎 1 , 2 lg 𝐴 SN , 𝜇 1 , 2 𝛾 SN , 𝜎 1 , 2 𝛾 SN , 𝜌 1 , 2 SN , and 𝜈 N ) . 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 commonspectrum process in the EPTA data, we impose a Gaussian hyperprior determined by the measurement uncertainty on ( lg 𝐴, 𝛾 ) 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  \nand (co-)variance - are well-constrained. Interestingly, we also find a marginal excess posterior density at ( 𝜎 2 lg 𝐴 SN , 𝜎 2 𝛾 SN ) = ( 0 , 0 ) . This result is related to Goncharov et al. (2022)'s measurement of 𝜎 lg 𝐴 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 4.", '5 CONCLUSIONS': 'Weperformedinference 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 9 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 (highperformance computing) at Swinburne University of Technology.', 'DATA AVAILABILITY': 'The data used for our study, EPTA DR2 by the (EPTA Collaboration 2023), is available at zenodo. The code to reproduce the results of this analysis and the analysis by Goncharov et al. (2022) is available in the epta branch at github.com/bvgoncharov/pta\\_gwb\\_priors. Other data may be provided by the corresponding author upon request. The PTA likelihood is incorporated in /e.pc/n.pc/t.pc/e.pc/r.pc/p.pc/r.pc/i.pc/s.pc/e.pc (Ellis et al. 2020) and posterior sampling is performed using /p.pc/t.pc/m.pc/c.pc/m.pc/c.pc/s.pc/a.pc/m.pc/p.pc/l.pc/e.pc/r.pc (Ellis & van Haasteren 2017) and /d.pc/y.pc/n.pc/e.pc/s.pc/t.pc/y.pc (Speagle 2020).', 'REFERENCES': 'Chen S., et al., 2021, MNRAS, 508, 4970 Desvignes G., et al., 2016, MNRAS, 458, 3341 Detweiler S., 1979, ApJ, 234, 1100 EPTA Collaboration 2023, A&A, 678, A48 EPTA Collaboration InPTA Collaboration 2023a, A&A, 678, A49 EPTA Collaboration InPTA Collaboration 2023b, A&A, 678, A50 Ellis J., van Haasteren R., 2017, jellis18/PTMCMCSampler: Official Release, doi:10.5281/zenodo.1037579, https://doi.org/10.5281/zenodo. 1037579  \nEllis J. A., Vallisneri M., Taylor S. R., Baker P. T., 2020, ENTERPRISE: En- hanced Numerical Toolbox Enabling a Robust PulsaR Inference SuitE, Zenodo, doi:10.5281/zenodo.4059815, https://doi.org/10.5281/ zenodo.4059815  \nFoster R. S., Backer D. C., 1990, ApJ, 361, 300  \nGoncharov B., et al., 2021a, MNRAS, 502, 478  \nGoncharov B., et al., 2021b, ApJ, 917, L19  \nGoncharov B., et al., 2022, ApJ, 932, L22  \nHazboun J. S., Simon J., Siemens X., Romano J. D., 2020, ApJ, 905, L6  \nHellings R. W., Downs G. S., 1983, ApJ, 265, L39  \nHobbs G. B., Edwards R. T., Manchester R. N., 2006, MNRAS, 369, 655  \nKeith M. J., et al., 2013, MNRAS, 429, 2161  \nFigure 4. Modelling of the common-spectrum process as a narrow normally-distributed cluster of pulsar noise using data-informed hyperpriors 𝜋 ( 𝚲 ) . We assume that pulsar noise parameters lg 𝐴 and 𝛾 are distributed as a mixture model of two Gaussians, where hyperparameters of the second Gaussian are chosen to match the measurement uncertainty on the common spectrum process reported by the EPTA Collaboration & InPTA Collaboration (2023b). Parameters 𝜎 with respective indices correspond to standard deviations of the distribution that governs lg 𝐴 and 𝛾 . Consistency of 𝜎 with zero has first been shown in Goncharov et al. (2022) using a more rigorous although less flexible methodology. Parameter 𝜈 N corresponds to the fraction of the first, wider normal distribution in a model. Solid lines correspond to including the second 𝛾 = 13 / 3 red noise term to pulsars to account for the presence of the common-spectrum stochastic process, dashed lines correspond to not including it. As expected, the contribution of the common-spectrum process becomes more visible in the latter case at 𝜎 lg 𝐴,𝛾 = 0 and 𝜈 N ≈ 0 . 3.  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nN', 'Kramer M., Champion D. J., 2013, Classical and Quantum Gravity, 30, 224009': 'Lentati L., Alexander P., Hobson M. P., Taylor S., Gair J., Balan S. T., van  \nHaasteren R., 2013, Phys. Rev. D, 87, 104021  \nNANOGrav Collaboration 2016, ApJ, 821, 13  \nNANOGrav Collaboration 2020, ApJ, 905, L34  \nNANOGrav Collaboration 2023a, ApJ, 951, L8  \nNANOGrav Collaboration 2023b, ApJ, 951, L10  \nPark R. S., Folkner W. M., Williams J. G., Boggs D. H., 2021, AJ, 161, 105 Reardon D. J., et al., 2023, ApJ, 951, L6  \nRenzini A. I., Goncharov B., Jenkins A. C., Meyers P. M., 2022, Galaxies, 10, 34  \nRosado P. A., Sesana A., Gair J., 2015, MNRAS, 451, 2417  \nSazhin M. V., 1978, Soviet Ast., 22, 36  \nShannon R. M., Cordes J. M., 2010, ApJ, 725, 1607  \nSpeagle J. S., 2020, MNRAS, 493, 3132  \nThrane E., Talbot C., 2019, Publ. Astron. Soc. Australia, 36, e010  \nZic A., et al., 2022, MNRAS, 516, 410  \nvan Haasteren R., 2024, ApJS, 273, 23  \nvan Haasteren R., Vallisneri M., 2014, Phys. Rev. D, 90, 104012', 'APPENDIX A: NUMERICAL MARGINALISATION OVER HYPERPARAMETERS': 'In this Section, we derive Equation 9 which describes a hierarchical posterior marginalised over hyperparameters 𝚲 . After taking an integral over 𝚲 , Equation 7 becomes  \nP( 𝜽 | 𝜹𝒕 ) = Z -1 L( 𝜹𝒕 | 𝜽 ) ∫ 𝜋 ( 𝜽 | 𝚲 ) 𝜋 ( 𝚲 ) 𝑑 𝚲 . (A1)  \nFollowing the formalism of importance sampling which was used to obtain Equation 8, 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 A1 by unity and expand the unity on the righthand side as a proposal distribution divided by itself. Rearranging the multipliers,  \nP( 𝜽 | 𝜹𝒕 ) = L( 𝜹𝒕 | 𝜽 ) 𝜋 ( 𝜽 |∅) Z ∫ 𝜋 ( 𝜽 | 𝚲 ) 𝜋 ( 𝜽 |∅) 𝜋 ( 𝚲 ) 𝑑 𝚲 . (A2)  \nNext, we use the expression for the expectation value of a probability density 𝑓 ( 𝑥 ) given the known probability density 𝑝 ( 𝑥 ) 3 :  \n⟨ 𝑓 ( 𝑥 )⟩ 𝑝 ( 𝑥 ) = ∫ 𝑓 ( 𝑥 ) 𝑝 ( 𝑥 ) 𝑑𝑥 ≈ 1 𝑛 s 𝑛 s ∑︁ 𝑖 𝑓 ( 𝑥 𝑖 ) , (A3)  \nwhere 𝑥 𝑖 are 𝑛 s samples from 𝑝 ( 𝑥 ) . Taking 𝜋 ( 𝚲 ) as 𝑝 ( 𝑥 ) , we arrive at Equation 9.  \nThis paper has been typeset from a T E X/L A T E X file prepared by the author.'}