© The Authors 2025

1. Introduction

Pulsar timing array (PTA, Foster & Backer 1990) experiments aim to detect ultra-low-frequency gravitational waves (GWs) in the nano- to milli-Hertz regime, where supermassive black hole binaries (SMBHBs) are a prime target, but cosmological sources are also possible. The Chinese Pulsar Timing Array (CPTA, Lee 2016) is a collaboration of various Chinese institutions using the Five-hundred-meter Aperture Spherical Telescope (FAST, Jiang et al. 2019) as the main instrument to observe millisecond pulsars with the highest single-dish sensitivity to date¹.

The detection of any GW signal with PTA data is challenging as it needs to be found amongst numerous other possible sources of noise. For a gravitational wave background (GWB), e.g. from a population of SMBHBs or from cosmological sources, the characteristic spatial correlation of the signal across a network of pulsars (Hellings & Downs 1983, henceforth referred to as HD) acts as a smoking gun. For individual GW sources, like a single binary, the resulting GW signal is deterministic. The signals from both, the GWB and individual sources, can be compared against a noise-only model. To be able to distinguish GW signals from pulsar related noise terms, robust noise models for each pulsar need to be employed in the analysis. Goncharov et al. (2021) and Chalumeau et al. (2022) have shown improvements in the detection significance of a GWB signal in PTA datasets using optimized noise models.

Globally, the efforts are coordinated in the International Pulsar Timing Array (IPTA, Verbiest et al. 2016; Perera et al. 2019) consortium, with regional member PTAs operating in Europe with the European Pulsar Timing Array (EPTA, Kramer & Champion 2013), North America with the North American Nanohertz Observatory for Gravitational waves (NANOGrav, McLaughlin 2013), Australia with the Parkes Pulsar Timing Array (PPTA, Manchester et al. 2013), and India with the Indian Pulsar Timing Array (InPTA, Tarafdar et al. 2022). Collaborations also exist in South Africa with the MeerKAT PTA (Miles et al. 2023) and in China with the CPTA. A series of coordinated papers searching for a GWB have been published by the member PTAs of the IPTA (Agazie et al. 2023a; EPTA Collaboration & InPTA Collaboration 2023a; Reardon et al. 2023a) and the CPTA (Xu et al. 2023) based on their most recent datasets. The MeerKAT PTA has also found evidence of a GW signal in its dataset (Miles et al. 2025a).

Different PTA collaborations have applied different methods to set and optimize the noise models used for their analysis (EPTA Collaboration & InPTA Collaboration 2023b; Reardon et al. 2023b; Agazie et al. 2023b; Miles et al. 2025b). Using the recent datasets, they find a consistent common signal with varying support for the HD correlation, which indicates some evidence of the GW origin of the measured signal. A detailed comparison between the NANOGrav, EPTA+InPTA, and PPTA results was performed by the IPTA (Agazie et al. 2024).

This paper is part of a series using the CPTA data release I (DR1) and expands upon different aspects of the overview work in Xu et al. (2023). Here, we analyse the single pulsar noise properties of 57 pulsars that form the CPTA DR1 with time spans up to 3.4 years and a weekly cadence (for details on the dataset see Xu et al., in prep.). Section 2 provides a summary of the analysis methods and different software packages used. The results of the analysis are presented in Section 3. We discuss comparisons with noise properties from other PTAs and conclude in Section 4.

2. Methods

2.1. Analysis

A more detailed description of the overall likelihood for the single pulsar noise model can be found in Chen et al. (2021) and Chalumeau et al. (2022), and references therein. The likelihood function with units of s⁻ⁿ can be written as

$$\begin{array}{c}\hfill L_{\mathrm{PSR}}=\frac{e^{-\frac{1}{2}{(\delta t_{\mathrm{post}})}^T{\mathbf{C}}^{-1}(\delta t_{\mathrm{post}})}}{\sqrt{{(2\pi )}^n|\mathbf{C}|}},\end{array}$$(1)

where the vector of the timing post-fit residuals, δt_post, and vector of the timing pre-fit residuals, δt_pre, have n elements with corresponding observation time, t_j, uncertainty, σ_j, and radio frequency, ν_j. They are related via the design matrix, M, of the timing model parameters, ξ (Hobbs et al. 2006): δt_post = δt_pre − Mξ. The timing model parameters can be analytically marginalized (van Haasteren et al. 2009; van Haasteren & Levin 2013; van Haasteren & Vallisneri 2014)². The noise covariance matrix, C, can be split into three components, white noise (WN), red noise (RN), and dispersion measure (DM) variation of the pulsar: C = C_WN + C_RN + C_DM.

Three parameters are applied to the initial estimate on the uncertainty, σ_j, of the times of arrival (TOAs) from the template matching method to test how well we can understand the white noise of the observing system and pulsar intrinsic properties. They are called EFAC, a multiplicative factor, EQUAD, an additional quadratic term, and ECORR, a white noise correlated between TOAs from a single observation. Together they result in the white noise covariance matrix

$$\begin{array}{c}\hfill {\mathrm{C}}_{\mathrm{WN},jk}={\mathrm{EFAC}}^2{\sigma }_j{\sigma }_k{\delta }_{\mathit{jk}}+{\mathrm{EQUAD}}^2{\delta }_{\mathit{jk}}+{\mathrm{ECORR}}^2{\mathrm{C}}_{\mathrm{ECORR},jk},\end{array}$$(2)

where δ_jk is the Kronecker delta symbol and the ECORR covariance matrix element C_ECORR, jk = 1, when the j-th and the k-th data points are of the same observing time but of different radio frequencies; and matrix element C_ECORR, jk = 0 otherwise. We note that ECORR is required for the CPTA data, as we produce 64 channelized TOAs. Thus, pulsar jitter can introduce white noise that is correlated between the TOAs of the same epoch (i.e. between different radio frequencies). Pulsar jitter noise should decrease as the observation time is increased averaging over more rotations. We tested this hypothesis by comparing the models of EQUAD and ECORR, whose values depend on the observation time length against the model where they do not. We find that most pulsars do not favour either model, thus we used the conventional time-independent model for the main analysis. More details can be found in Appendix A.

The average ⟨ ⟩ of many realizations of the timing residuals can be described by a stochastic process. Using the Wiener-Khinchin theorem, the covariance, C_X, jk, between two elements of the residual time series δt_X, j and δt_X, k at times t_j and t_k with radio frequencies ν_j and ν_k can be defined using a one-sided (f > 0) power spectral density (PSD) S_X(f, ν):

$$\begin{array}{c}\hfill {\mathrm{C}}_{\mathrm{X},jk}=⟨\delta t_{\mathrm{X},j}\delta t_{\mathrm{X},k}⟩={\displaystyle {\int }_0^{\infty }}{S}_{\mathrm{X}}(f,\nu )e^{2\pi if(t_j-t_k)}\mathrm{d}f,\end{array}$$(3)

where X can be either RN or DM and ν² = ν_jν_k. For RN, the PSD is independent of the radio frequency, while for DM it is not. If the PSD is constant across all frequencies, this equation can also describe WN.

For long-term variations in the residual time series, also known as red noise, the PSD S_RN(f) is usually modelled as a power law in the Fourier frequency domain:

$$\begin{array}{c}\hfill S_{\mathrm{RN}}(f)=\frac{A_{\mathrm{RN}}^2}{12{\pi }^2}{\left(\frac{f}{f_{\mathrm{yr}}}\right)}^{-{\gamma }_{\mathrm{RN}}}{\mathrm{yr}}^3,\end{array}$$(4)

with dimensionless amplitude A_RN and spectral index γ_RN sampled at frequencies f scaled by f_yr = 1/yr. The integral over an infinite frequency range in Eq. (3) can in practice be approximated by a sum over a limited number of frequency bins.

Two common methods are used to model the DM variations. DM variations can be measured in a piece-wise fashion between subsequent observations to produce a DM time series. These DM values can be applied to the residuals to subtract the effect of variations of the interstellar medium. This model, called DMX (NANOGrav Collaboration 2015), can be compared against a DM variation model that assumes a smooth variation over time with a Gaussian process (GP). Similarly to Larsen et al. (2024), we directly compared the results from a DM GP to those from a DMX-style data set. The DM GP models the DM variation over time after fitting for a constant, linear, and quadratic term.

Similar to the red noise, we used a power law to model the timing residuals induced by the remaining small DM variations. In the literature, there are two common definitions of the power-law model: TEMPONEST (TN, Lentati et al. 2014) and ENTERPRISE (EP, Ellis et al. 2020). In the TN convention, the PSD S_DM(f, ν) is written as

$$\begin{array}{c}\hfill S_{\mathrm{DM}}(f,\nu )=\frac{A_{\mathrm{DM},\mathrm{TN}}^2}{K^2{\nu }^4}{\left(\frac{f}{f_{\mathrm{yr}}}\right)}^{-{\gamma }_{\mathrm{DM}}}{\mathrm{yr}}^3,\end{array}$$(5)

where A_DM, TN has units of pc cm⁻³ s⁻¹ and K = 2.41 × 10⁻⁴ pc cm⁻³ MHz⁻² s⁻¹ is the DM constant. The power law can also be written in the EP convention as

$$\begin{array}{c}\hfill S_{\mathrm{DM}}(f,\nu )=\frac{A_{\mathrm{DM},\mathrm{EP}}^2}{12{\pi }^2}{\left(\frac{f}{f_{\mathrm{yr}}}\right)}^{-{\gamma }_{\mathrm{DM}}}{\left(\frac{\nu }{{\nu }_{\mathrm{ref}}}\right)}^{-4}{\mathrm{yr}}^3,\end{array}$$(6)

where A_DM, EP is dimensionless and ν_ref = 1400 MHz is a reference radio frequency. Comparing the two conventions gives a relation between the TN and EP amplitudes:

$$\begin{array}{c}\hfill A_{\mathrm{DM},\mathrm{EP}}=\sqrt{12{\pi }^2}{A}_{\mathrm{DM},\mathrm{TN}}/\left(K{\nu }_{\mathrm{ref}}^2\right).\end{array}$$(7)

Note that in this work, we employ the TN convention and denote A_DM, TN = A_DM in the following.

Additionally, we also model deterministic variations with a yearly period in the DM time series:

$$\begin{array}{c}\hfill \delta {\mathrm{DM}}_{\mathrm{yr},j}=A_{\mathrm{yr}}sin(2\pi f_{\mathrm{yr}}{t}_j+{\varphi }_{\mathrm{yr}}),\end{array}$$(8)

where A_yr in units of pc cm⁻³ is the amplitude of the yearly variations in the DM time series and ϕ_yr is the phase of the yearly sinusoid. The exact value of ϕ_yr depends on what time is chosen to correspond to ϕ_yr = 0.

The DM variations induced timing residuals δt_DMY in s and the DM time series δDM_Y in pc cm⁻³ are related via (Lentati et al. 2014)

$$\begin{array}{c}\hfill \delta {\mathrm{DM}}_{\mathrm{Y},j}=K{\nu }_j^2\delta t_{{\mathrm{DM}}_{\mathrm{Y}},j},\end{array}$$(9)

with Y denoting either the GP model or yearly variations.

To robustly model the red and DM noise as power laws in an efficient way, the Fourier frequency components need to be chosen carefully. While in theory a large enough number of frequency bins should provide an accurate fit, recent works have shown that in practice restricting the number of frequency bins can limit the covariance between white and red noise (Chalumeau et al. 2022; EPTA Collaboration & InPTA Collaboration 2023b). In addition, such a covariance could be falsely attributed to a GW signal, see e.g. Falxa et al. (2023). A small number of frequency bins will also reduce the computational cost. The conventional choice in PTA analysis has been 1/T, 2/T, up to n_high/T, linearly spaced frequency bins with width 1/T, up to a maximum number n_high. This approach has been shown to work well and efficiently, if there is a strong and dominant red signal in the lowest frequency bins. One can use logarithmically spaced frequency bins if the white noise is the dominant signal. A mixture between the two has been proposed by van Haasteren & Vallisneri (2015), which should result in an accurate model that is also efficient to evaluate computationally. However, that model requires knowledge of the transition frequency between red and white noise. For the CPTA analysis, we used n_bin logarithmically and evenly spaced frequency bins between 1/T and n_high/T. We applied the Romberg summation weights to the frequency bins to increase the accuracy of the model, details of which can be found in Appendix B.

Following the same principle as Chalumeau et al. (2022), we searched for the optimal frequency choice for the pulsar red and DM variation noise. We also used the Bayesian method to determine the optimal largest frequency bin, n_high. Additionally, we also optimized for the lowest total number of frequency bins, n_bin. The Bayes factor, ℬ, can be computed by comparing the evidence, 𝒵_i, between two models: ℬ₂₁ = 𝒵₂/𝒵₁. If ℬ₂₁ > 1, then model 2 is favoured over model 1, while model 1 is favoured over model 2, if ℬ₂₁ < 1. For a scale of significance see e.g. Kass & Raftery (1995).

2.2. Software

We used the TEMPO2 (Hobbs et al. 2006) package for the timing model design matrix and the residuals. To compute the likelihood function we employed four different software packages and samplers for the Bayesian analysis.

TEMPONEST (Lentati et al. 2014) is a Bayesian pulsar timing analysis software package written in the C programming language, which employs a nested-sampling approach (Feroz & Hobson 2008). As an integrated plugin of TEMPO2, it models the red noise and DM variation using either a power-law model or a model-independent method. Within our analyses, we employed MULTINEST (Feroz & Hobson 2008) as the sampler.

ENTERPRISE (Ellis et al. 2020) is a suite of python-based scripts that can be used with the timing model fit and residuals from TEMPO2 using libstempo (Vallisneri 2020). Signals and noises are modelled as Gaussian processes, typically with power-law spectra (van Haasteren & Vallisneri 2014). The likelihood computation from ENTERPRISE is passed to PTMCMCSAMPLER (Ellis & van Haasteren 2017) for the Bayesian parameter estimation via Monte Carlo Markov chain sampling.

42 (Caballero et al. 2016) is a software package written in python for Bayesian PTA data analysis. The red noise and DM terms can be modelled with power-law spectra or other models. The sampler MULTINEST (Feroz & Hobson 2008) is used for the single pulsar noise analysis.

42++ is a recently developed Bayesian PTA data analysis software written in the C++ language. It relies on the intel OneAPI package to implement high performance numerical linear algebra operations. The software package supports Metropolis–Hastings and parallel tempering sampling.

The prior choice for the parameters can be found in Table 1. For 42++, different priors for the power-law processes are adopted. Instead of the spectral amplitude A_X and index γ_X, 42++ uses the characteristic amplitude A_c, X and index α_X to specify the power-law model, such that the PSD becomes

$$\begin{array}{c}\hfill S_{\mathrm{X}}(f)=\frac{A_{\mathrm{c},\mathrm{X}}^2}{f}{\left(\frac{f}{f_{\mathrm{yr}}}\right)}^{-2{\alpha }_{\mathrm{X}}}.\end{array}$$(10)

The characteristic amplitude A_c, X has the intuitive physical meaning of the strength of the signal at a yearly timescale. For red noise, A_c, RN has the same unit as the timing residuals, i.e. second, while for DM, A_c, DM has the DM unit, i.e. pc cm⁻³.

The relations between the different amplitudes can be derived from comparing the PSDs from Eq. (4) with Eq. (10), and Eq. (5) with Eq. (10) by setting γ_X = 2α_X + 1 for both RN and DM³:

$$\begin{array}{c}\hfill A_{\mathrm{c},\mathrm{RN}}=A_{\mathrm{RN}}\mathrm{yr}/\sqrt{12{\pi }^2},A_{\mathrm{c},\mathrm{DM}}=A_{\mathrm{DM}}\mathrm{yr}.\end{array}$$(11)

3. Results

The noise optimization and selection was done in three steps:

1. The choice of the number of frequency bins and the high frequency cutoff was determined in a Bayesian model comparison.

2. The single pulsar noise analysis was run with these choices with our four software packages.

3. The final noise terms were selected based on the highest Bayes factor.

The optimal frequency bin choice is given in Table C.1, where n_bin is the number of logarithmically and evenly spaced bins between f₀ = 1/T and f_high = n_high/T, with T being the time span from the earliest to the latest TOA of a given pulsar’s dataset. In our optimization and selection processes, we modelled the high-frequency cutoff as one of the noise model parameters. We then ran the Bayesian inference in step 1 using 42++ only with different choices of number of frequency bins, where we tried 33, 65, 129, 257, and 513 (i.e. 2ⁿ + 1, where n is the order of the integration). These numbers satisfy the requirement of the Romberg integration. For all runs, we note that the values of the Bayesian evidence are nearly the same if the number of frequency bins is large enough. To reduce the computational cost, we chose the minimal number of frequency bins, at which the Bayesian evidence is ‘not significantly worse’ than the best run with the highest Bayesian evidence. Here, ‘not significantly worse’ means that the difference between the logarithmic Bayesian evidence is less than 1. The choice was applied to both the red and DM noise. In most cases, a small number is sufficient to describe pulsars without red or DM noise. Pulsars with noticeable red or DM noise need in general a slightly large number for n_high to reach into the high frequencies. In most pulsars, 33 bins are enough to describe the power law sufficiently well. A few pulsars favour the inclusion of up to 129 frequency bins.

Our standard noise model for step 2 includes three white noise, two red noise, two DM GP, and two parameters for the annual DM variation. This step was performed using all four software packages. However, the Romberg weights were only applied in the 42++ analysis. The tables in Appendix D show all these parameters for the four different pipelines. The Bayesian sampling is done with all nine (five) parameters for the DM GP (DMX) models simultaneously. Time series realizations of the red and DM GP noise can be found for two example pulsars in Figure 1. In general, all four pipelines give very consistent constraints or upper limits on the parameters, as demonstrated in Figure 2.

Fig. 1. Time series reconstructions for J0509+0856 in the top three and J1640+2224 in the bottom three panels. From top to bottom, each set of three panels shows the red noise for the DM GP and DMX datasets, as well as the DM GP reconstructions compared the DMX values (where linear and quadratic trends have been removed). The blue points indicate the median residual values and 1σ uncertainties.

Fig. 2. Recovered noise parameters from the four different software packages. The median values and central 68% credible regions are indicated by the points and error bars. The 95% upper limits are shown with arrows pointing down.

As several parameters are poorly constrained by the CPTA DR1 with its approximately three years of observations, the inclusion of the corresponding noise terms could increase the computational cost unnecessarily and possibly decrease the sensitivity of the GW signal search. Therefore, we performed another round of model comparisons investigating EQUAD, ECORR, DM and red noise to determine which of the noise models is the most favoured. The analysis was done simultaneously with fitting timing model parameters. Step 3 was run with TEMPONEST only. The most favoured combination of noise terms can be found in Table 2 of Xu et al. (in prep.).

3.1. DM GP versus DMX

We also performed steps 1 and 2 of the analysis for the dataset obtained with the DMX model. Figure 3 shows a comparison of the median and central 68% for the white and red noise parameter posteriors between the DMX and DM GP models.

Fig. 3. Comparison between the white and red noise parameter constraints from the DM GP and DMX model with ENTERPRISE. The median values and central 68% credible regions are indicated by the points and error bars. The 95% upper limits are shown with arrows pointing down.

The recovered EFAC and EQUAD white noise parameters are identical for all pulsars. There are, however, noticeable differences in the posteriors of the ECORR parameter. This could be due to power at high frequencies that the DMX model can fit, while the power-law model has a high frequency cutoff. On the other hand, the DMX model could also be overfitting and absorb high frequency noise, which could lead to lower Bayesian evidence.

Most pulsars show broadly consistent red noise constraints. However, a general trend of larger upper limits in the DMX model can be seen. It is possible that some of the pulsar intrinsic red noise has been absorbed by the DMX model or that high frequency white noise is wrongly attributed to red noise by the DM GP model.

We would also like to point out that in all cases the Bayesian evidence for the DMX dataset is significantly lower compared to the DM GP model. We note that for most of our pulsars the logarithmic Bayesian evidence difference Δlnℬ comparing the DM GP to the DMX model reaches the level of a few thousand. Thus, we focus on the analysis with the DM GP model⁴.

3.2. Parameter estimation

We summarize the noise properties of the 57 pulsars of the CPTA DR1. They are classified into four categories: 1. white noise only, 2. DM and white noise, 3. red and white noise, and 4. DM, red and white noise (see Table 2).

The EFAC white noise parameters are in general close to the expected value of 1, indicating that the initial estimate does not need to be modified much. Twenty pulsars have neither EQUAD nor ECORR constraints, thus we can only set upper limits for them (see Table 3). Out of the remaining 37 pulsars, 19 (28) have constraints on EQUAD (ECORR) respectively, with 10 having both. The appearance of the ECORR parameter is expected as the observations are channelized, i.e. several TOAs are created from the same observation by splitting the full radio frequency bandwidth into 64 channels.

About 14 pulsars have constrained parameter values for the red noise power law and 32 pulsars have constrained values for the DM GP. For the other pulsars we set upper limits, however, there are 15 (9) additional pulsars that have Bayesian evidence supporting the inclusion of the red noise (DM variations) in the model.

Five pulsars are found to favour only white noise models: J0509+0856, J0605+3757, J0732+2314, J1710+4923, and J2019+2425, where the first two only require EFAC and the latter three have EFAC and ECORR. For 23 pulsars we found significant evidence for a DM variation modelled with a power-law GP in addition to the white noise, but no evidence for pulsar red noise. Evidence to include red noise on top of white noise was found in 11 pulsars. Eighteen pulsars have the highest Bayes factor for the full noise model with white, red, and DM noise all being measured.

All four software packages produce consistent results for almost all pulsars, with a few notable exceptions. J1327+3423 has a bimodal posterior distribution for the ECORR parameter. While ENTERPRISE only finds the lower of the two peaks, the other three software packages produce a bimodal distribution, thus producing a larger upper limit. This could be either due to the differences in the implementation or the different sampling methods. Other examples include J1643−1224, J2010−1323, and J2033+1734, where small differences in the tails of the distributions lead to different upper limits. Additionally, some pulsars, like J1832−0836, J1903+0327, and J1910+1256, have constrained values in some software packages, while the other software packages produce upper limits.

3.3. Parameter covariances

Next, we investigate the covariances between different noise terms and parameters. The covariance between red noise and DM variations is expected as both cause long-term trends in the residuals. This covariance can be broken with multiple radio frequency coverage and longer observation time spans. For the amplitude parameter, posterior covariances usually result in distributions that have a peak, but also a long tail, resulting in the typical ‘L’-shape in the 2D posterior distributions. Table 4 summarizes the number of pulsars in the different categories detailed below. The full 2D posteriors can be found in the online supplementary material⁵.

3.3.1. DM GP – DM annual

For the CPTA DR1 with three years of observations, the annual DM signal is highly covariate with the DM GP power-law model. Although we only measure constrained values for the annual DM amplitude A_yr for about three pulsars (J0340+4130, J2033+1734, and J2322+2057), five more show an indication of a measurement (J0023+0923, J0824+0028, J1012+5307, J2234+0611, and J2234+0944), and 12 show a strong correlation between the DM GP power law and annual amplitude (J0154+1833, J0613−0200, J1643−1224, J1713+0747, J1738+0333, J1911+1347, J1918−0642, J1944+0907, J2017+0603, J2033+1734, J2145−0750, and J2317+1439). Longer observation time spans will help to break degeneracies.

3.3.2. Red noise – DM

The expected correlation between pulsar red noise and DM GP variations have been found in ten pulsars. However, five of these also show covariance with the ECORR white noise parameter (see next section). The five pulsars with purely red and DM noise covariance are J0030+0451, J0034−0534, J1853+1303, J1918−0642, and J1923+2515.

3.3.3. Red noise – DM – white noise

In addition to the covariance between red and DM noise, we also measured some degree of interplay between ‘long’-term variations and ‘short’-term white noise. Given that the time span is only three years, the distinction between red, DM, and white noise can be difficult. In total, 19 pulsars show difficulties in distinguishing the three noise terms. Five of these pulsars have a strong covariance between the ECORR parameter and both red and DM noise: J0154+1833, J0824+0028, J1903+0327, J2317+1439, and J2322+2057. The remaining 14 pulsars either show covariance between red and white noise (four pulsars) or DM and white noise (ten pulsars).

3.3.4. Red noise – white noise

As a sub-class of the red-DM-white noise covariance, ten pulsars only show interplay between the red and white noise. The DM parameters are well constrained for four pulsars (J0613−0200, J0621+1002, J1643−1224, and J2229+2643), while they are set as upper limits for the other six pulsars (J0645+5158, J1327+3423, J1453+1902, J1630+3734, J2010−1323, and J2022+2534).

3.3.5. DM – white noise

Another sub-class of the red-DM-white noise covariance are the DM-white noise pulsars, which are the opposite of red-white noise case. Only four pulsars show this type of interplay: J1910+1256, J1955+2908, J2033+1734, and J2214+3000. In all cases, the red noise term is not constrained.

3.3.6. No covariances

Thirty-three pulsars in the CPTA DR1 show no covariances between the noise terms. This is the case for the 26 pulsars that have well-constrained parameter distributions for red and/or DM noise, as well as seven pulsars that show no sign of red or DM noise terms. The seven pulsars with white noise only are J0509+0856, J0605+3757, J0732+2314, J1710+4923, J1944+0907, J2019+2425, and J2145−0750.

3.3.7. Covariances in the DMX dataset

As a comparison, we briefly summarize the constrained parameters and their covariances in the DMX model with only white and red noise. The figures showing the full 2D posterior distributions can be found in the online supplementary material⁶. Two pulsars show interplay between the EQUAD and ECORR parameters: J1955+2908 and J2150−0326. The similarity between the two parameters could explain the covariance. Eighteen pulsars have an indication for red noise, with ten of them with constrained parameters and eight of them sharing a covariance with the ECORR parameter.

4. Discussion and conclusion

We compared the parameter constraints for the red noise and DM GP power laws with the EPTA DR2 (EPTA Collaboration & InPTA Collaboration 2023b), NG15 (Agazie et al. 2023b), and PPTA DR3 (Reardon et al. 2023b) for a total of 22 overlapping pulsars, see Figure E.1. As each PTA uses different settings to produce the final noise models, the figures should be viewed with these caveats in mind. The EPTA DR2 has two versions of the dataset and only includes a noise term if it is significant. The PPTA DR3 noise analysis includes other noise terms that we have not considered. Finally, NANOGrav primarily uses the DMX model, thus there are no DM GP constraints to compare against. One can see very broad agreement between the CPTA results and those from the other PTAs. Agazie et al. (2024) provide a more detailed comparison between the noise models from the different PTAs.

The CPTA DR1 obtained using about three years of FAST observations is an unprecedentedly sensitive PTA dataset for GW searches. As such, we have investigated in detail the timing and noise analysis of each of the 57 pulsars in the dataset. The results for the timing model parameters and other details of the dataset can be found in Xu et al. (in prep.). Here, the results of the noise models have been presented. We followed a three-step process to determine the optimal frequency bin choice for the power-law model, verify the parameter posterior distributions with four different software packages, and finally use Bayesian analysis to select the optimal noise terms to be included in the final timing and noise model.

The analysis was performed for both DM GP and DMX style datasets. The DMX model allows for larger red noise amplitudes for the majority of the pulsars and has significantly lower Bayesian evidence compared to the DM GP model. The DM GP dataset has been chosen as our main dataset. In general, the four software packages, TEMPONEST, ENTERPRISE, 42, and 42++, produce very consistent constraints on all the parameters in the noise model.

For the majority of pulsars, we do not measure red or DM noise and therefore set upper limits on the amplitudes. Fourteen (32) pulsars show constraints for the red (DM) noise parameters, while an additional 15 (9) have Bayesian evidence for the respective noise term to be included (Xu et al., in prep.). We also tested the conventional description for the EQUAD and ECORR parameters against a time-dependent jitter noise description and found no conclusive evidence for either model. The hypothesis that ECORR should be linked to the pulsar jitter process is further tested in Wang et al. (in prep.).

We find no covariances between the parameters of the different noise terms for 33 pulsars, the majority of which have constrained posteriors for the parameters. The remaining 24 pulsars display covariances between red, DM, and white noise, notably ECORR. Given the short length of the dataset, it may be difficult to clearly distinguish between long-term red/DM variations and short-term white noise. While in this paper we focus on the noise terms with constrained parameters and their covariances, the final timing and noise models used for the GW search (Xu et al. 2023) were obtained from Bayesian model comparison in Xu et al. (in prep.).

The CPTA DR1 is opening a new frequency range that has not been accessible by the established PTA collaborations, with much increased sensitivity and reduced uncertainty using FAST observations. Being only about three years long, long-term variations cannot yet be well probed with the CPTA data; continued observation of the 57 pulsars will allow us to break the degeneracies found between different noise terms with future CPTA datasets. We plan to join the efforts in the IPTA to combine datasets from all PTA collaborations in the future. Such an IPTA dataset will significantly improve the GW science that can be done with a PTA.

Data availability

Supplementary material is available at https://doi.org/10.5281/zenodo.13828113

Acknowledgments

Observation of CPTA is supported by the FAST Key project. This work has used the data from the Five-hundred-meter Aperture Spherical radio Telescope (FAST, https://cstr.cn/31116.02.FAST). FAST is a Chinese national mega-science facility, operated by National Astronomical Observatories, Chinese Academy of Sciences (NAOC). This work is supported by the National SKA Program of China (2020SKA0120100), the National Key Research and Development Program of China No. 2022YFC2205203, the National Natural Science Foundation of China grant no. 12041303, 12250410246, 12173087 and 12063003, China Postdoctoral Science Foundation No. 2023M743518 and 2023M743516, Major Science and Technology Program of Xinjiang Uygur Autonomous Region No. 2022A03013-4, the CAS-MPG LEGACY project, and funding from the Max-Planck Partner Group. KJL acknowledges support from the XPLORER PRIZE and 20-year long-term support from Dr. Guojun Qiao. The data analysis are performed with computer clusters DIRAC and C*-SYSTEM of PSR@pku and computational resource provide by the PARATERA company. Software package 42++ is developed with INTEL oneAPI toolkits and the SCIENCE EXPLORER’S DEVELOPING GEARS from WEICHUAN TECHNOLOGY. We thank Dr. Jin Chang for providing the valuable long-term support to the CPTA collaboration. CPTA thanks Dr Duncan Lorimer, Dr Michael Kramer and Dr Richard Manchester for helping initialize the collaboration. SC would like to thank Dr Bruce Allen for his comments on the description of the analysis methods.

References

Appendix A: Comparison between jitter noise models

We compare two different models for the jitter noise, where, in model I, the white noise is modeled using EQUAD and ECORR as explained in the main text Section 2.1; and, in model II, the white noise is modeled as

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathrm{C}}_{\mathrm{WN},jk}=& {\mathrm{EFAC}}^2{\sigma }_j{\sigma }_k{\delta }_{\mathit{jk}}+{\mathrm{EQUAD}}^2{\left(\frac{1\mathrm{hr}}{t_{\mathrm{obs}}}\right)}^2{\delta }_{\mathit{jk}}\hfill \\ \hfill & +{\mathrm{ECORR}}^2{\left(\frac{1\mathrm{hr}}{t_{\mathrm{obs}}}\right)}^2{\mathrm{C}}_{\mathrm{ECORR},jk},\hfill \end{array}\end{array}$$(A.1)

i.e. both the uncorrelated jitter (EQUAD) and correlated jitter (ECORR) are modeled with their values being proportional to the observation time length $\sqrt{t_{\mathrm{obs}}}$. We list the Bayes factors between the two models in Table A.1. As one can see we can not tell the preference between the two models for most of the pulsars as the observation time length of each epoch is roughly a constant. For J1713+0747, J2033+1734 the time dependent jitter noise model is preferred with $ln\mathcal{B}\ge 3$, while the standard EQUAD-ECORR model is only preferred for J1643−1224. It is interesting to note that the jitter noise model is not systematically better than the EQUAD-ECORR model.

Appendix B: Numeric integration weights and covariance matrix computation

The covariance matrix C_j, k with j and k being the indices for time epochs and power spectral density function S(f) with f being in the frequency domain are related to each other via the Fourier transform

$$\begin{array}{c}\hfill {\mathrm{C}}_{\mathit{ij}}={\displaystyle \int S}(f)cos(2\pi f(t_i-t_j))\mathrm{d}f.\end{array}$$(B.1)

The above equation can also be computed by integrating over dlnf to get

$$\begin{array}{c}\hfill {\mathrm{C}}_{\mathit{ij}}={\displaystyle \int S}(f)fcos(2\pi f(t_i-t_j))\mathrm{d}lnf.\end{array}$$(B.2)

In order to speed up the numerical computation the PTA community usually use the frequency domain approximation of the covariance matrix, i.e.

$$\begin{array}{c}\hfill \mathbf{C}={\mathbf{F}\mathbf{S}\mathbf{F}}^{\mathrm{T}},\end{array}$$(B.3)

which is a consequence of approximating Eqn. B.2 using the following summation

$$\begin{array}{c}\hfill {\mathrm{C}}_{\mathit{jk}}={\displaystyle \sum _l}{S}_l(c_j{c}_k+s_j{s}_k),\end{array}$$(B.4)

and casting the equation into the matrix form of Eq B.3, where c_j ≡ cos2πf_lt_j, s_j ≡ sin2πf_lt_j, f_l is the discrete frequency, and S_l = S(f_l)Δf. Here, l increases up to 2, 3, 5, 9, ..., 2ⁿ + 1, where n is the order of integration. For the integration over dlnf, we have S_l = S(f_l)f_lΔlnf. However, the approximation in Eqn. B.4 is far from being accurate for most applications (Press et al. 2007). The difference between Eqn. B.2 and B.4 is on the order of O(Δf²). One widely used numerical method is to insert weights in the summation, such that the difference between the summation and integration can be of higher order of Δf. For sufficiently smooth functions the Romberg integration method (Press et al. 2007) can push the precision to the order of Δf²ⁿ using a slightly modified version of the summation from Eqn. B.4:

$$\begin{array}{c}\hfill {\mathrm{C}}_{\mathit{jk}}={\displaystyle \sum _l}{S}_l{W}_l(c_j{c}_k+s_j{s}_k),\end{array}$$(B.5)

where W_l are the Romberg weights (Press et al. 2007; Chapra & Canale 2014). They can be pre-computed using R(n, m) representing Romberg n-order integration with 2ⁿ + 1 bins and m-order extrapolation (m ≤ n). For the l-th bin, R_l(n, m) is defined recursively as

$$\begin{array}{c}\hfill R_l(n,0)=\{\begin{array}{cc}\frac{1}{2},\hfill & \mathrm{for}l=1,2^n+1;\hfill \\ 1,\hfill & \mathrm{for}l=2,...,2^n(\mathrm{if}n\ge 1);\hfill \end{array}\end{array}$$(B.6)

$$\begin{array}{c}\hfill R_l(n,m)=\frac{1}{4^m-1}[4^m{R}_l(n,m-1)-R_l(n-1,m-1)].\end{array}$$(B.7)

Extrapolating m to the n-order gives the Romberg weights for the l-th bin (up to 2ⁿ + 1), W_l = R_l(n, n). In 42++, all the covariance matrices are computed with Eqn. B.5. Our numerical experiments show that using Romberg weights can improve the precision of C_jk by two to three orders of magnitude for most of the number of frequency bins and spectral indices encountered in the application of PTAs, although the improvement of the likelihood evaluation precision can be limited by other factors. More details will be published elsewhere.

Appendix C: Frequency binning and properties of the data set

Appendix D: Posterior constraints on the noise parameters

Appendix E: Comparison with results from other PTAs

Fig. E.1. Comparison of the red noise and DM GP power law parameters between the CPTA DR1, EPTA DR2full and DR2new, NANOGrav 15yr and PPTA DR3 data sets.

All Tables

All Figures

	Fig. 1. Time series reconstructions for J0509+0856 in the top three and J1640+2224 in the bottom three panels. From top to bottom, each set of three panels shows the red noise for the DM GP and DMX datasets, as well as the DM GP reconstructions compared the DMX values (where linear and quadratic trends have been removed). The blue points indicate the median residual values and 1σ uncertainties.
In the text

	Fig. 2. Recovered noise parameters from the four different software packages. The median values and central 68% credible regions are indicated by the points and error bars. The 95% upper limits are shown with arrows pointing down.
In the text

	Fig. 3. Comparison between the white and red noise parameter constraints from the DM GP and DMX model with ENTERPRISE. The median values and central 68% credible regions are indicated by the points and error bars. The 95% upper limits are shown with arrows pointing down.
In the text

	Fig. E.1. Comparison of the red noise and DM GP power law parameters between the CPTA DR1, EPTA DR2full and DR2new, NANOGrav 15yr and PPTA DR3 data sets.
In the text