© The Authors 2025

1. Introduction

Pulsar timing array (PTA) collaborations search for gravitational waves (GWs) with frequencies in the nanoHertz (nHz) band (10⁻⁹ − 10⁻⁷ Hz). The primary target of PTA experiments is a stochastic gravitational wave background (GWB) generated by the incoherent superposition of the GWs emitted by all the super-massive black hole binaries (SMBHBs) that populate our universe (Rajagopal & Romani 1995; Jaffe & Backer 2003; Wyithe & Loeb 2003; Sesana et al. 2008). The signal manifests itself in the data as a red noise process common to all the pulsars in the array, with a specific angular correlation, first derived by Hellings & Downs (1983). Recently, the PTA community found evidence for a signal with both these characteristics: the collaborations reporting this finding were the European PTA collaboration, along with the Indian PTA collaboration (EPTA, InPTA; EPTA Collaboration et al. 2023a), North American Nanohertz Observatory for Gravitational Waves collaboration (NANOGrav, Agazie et al. 2023a), Parkes PTA collaboration (PPTA, Reardon et al. 2023) and Chinese PTA collaboration (CPTA, Xu et al. 2023).

While the Hellings and Downs (HD hereinafter) correlation is distinctive of a GW origin of the signal, separating it from possible correlated noise of different nature (Tiburzi et al. 2016), its source has not beeen confidently established yet. In fact, besides SMBHBs, a nHz cosmic GWB can also be produced by a number of physical processes occurring in the early universe, such as cosmic strings networks (Damour & Vilenkin 2000), primordial curvature perturbations (Tomita 1967), or QCD phase transitions in the early universe (Kosowsky et al. 1992), all of which can provide a viable explanation of the observed signal (Afzal et al. 2023; EPTA Collaboration et al. 2024). Moreover, the HD correlation is not unique of stochastic GWBs; as shown in Cornish & Sesana (2013) and in Allen (2023), the signal generated by a single SMBHB and observed by an array of pulsars results in an angular correlation with the HD expectation value, albeit with a somewhat different variance. We note that the quality of the current data does not allow for a description of the signal variance; therefore, the possibility that a particularly loud source (e.g. a massive and nearby SMBHB) could be mistaken by a GWB is real. The term ‘loud source’ is used here to describe one that can be singled out from the noise and (when present) from the GWB by means of a deterministic, template-based search. Since the SMBHBs producing such signals are far from mergers and persistent in the data, their signal is generally referred to as continuous gravitational wave (CGW, Babak & Sesana 2012).

Searches for CGWs have been performed on the latest data release of EPTA (DR2new, Antoniadis et al. 2024) and on the NANOGrav 15 yr dataset (Agazie et al. 2023b). These studies both found evidence supporting the presence of a CGW with a 2.5–3σ significance, compared to a model where only pulsar noise is present in the data. This value is comparable to the evidence in favour of the presence of a GWB, which has a 3-to-4σ significance. Some simulations carried out in Antoniadis et al. (2024) already showed how a CGW can be misinterpreted as a GWB and (vice versa) a GWB can be misinterpreted as a CGW. In both the EPTA and NANOGrav analyses, when a GWB is also searched for, evidence for a CGW significantly decreases, indicating that its presence is (at the very least) not necessary to explain the observations. However, the degeneracy between the pulsar response to the two signals does not allow us to conclusively rule out the presence of a CGW.

The presence of these resolvable sources in the PTA data on the one hand poses a number of interesting challenges, on the other hand opens the door to invaluable opportunities. In fact, while the analysis pipelines currently used to search for a GWB rely on the assumption that the background is a Gaussian, isotropic, and stationary process characterised by a power-law Fourier spectrum (Lentati et al. 2015; Arzoumanian et al. 2016), the presence of a CGW can result in an excess of power at specific frequencies; these produces some level of anisotropy and non Gaussianity (e.g. Gardiner et al. 2024; Sah et al. 2024). Thus, the standard analysis may depend on inapplicable assumptions and yield biases on the estimated GWB parameters (Bécsy et al. 2023; Valtolina et al. 2024). This is a recognised issue within the PTA community and analysis techniques are being developed to go beyond the power-law, isotropic Gaussian assumption (e.g. Taylor & Gair 2013; Gair et al. 2014; Taylor et al. 2020; Ali-Haïmoud et al. 2020; Xin et al. 2021; Bécsy et al. 2022) and to draw GWBs and CGWs out of those data simultaneously (Bécsy & Cornish 2020; Freedman & Vigeland 2024).

Within this context, it is of paramount importance to study the interplay of CGWs and GWB when both are present in the data and what challenges arise in the interpretation of the analysis results. From an exquisitely scientific perspective, conversely, the prospects of detecting CGWs is particularly appealing. In fact due to their deterministic nature, PTAs are able to estimate the source parameters, in particular the sky location via triangulation (Babak & Sesana 2012; Ellis et al. 2012) and (possibly) the distance to the source (especially if the pulsar term of the induced timing delay can be measured, Ellis 2013). The detection of one or more resolvable sources (which is expected in the SKA PTA era, Truant et al. 2025) would be very valuable, as it would allow us to track the orbital evolution of the binary and thus obtain direct information on the massive black holes (MBHs) residing in the environment (e.g. Kocsis & Sesana 2011; Sesana 2013a; Ravi et al. 2014). This information cannot be well constrained by a measure of the GWB power spectral density, because it is strongly degenerate with the MBH mass and eccentricity distribution. In addition, CGWs offer the exciting opportunity to search for electromagnetic counterparts of the GW signal, opening up the prospects for low-frequency GW multimessenger astronomy (Sesana et al. 2012; Burke-Spolaor 2013; Kelley et al. 2019; Liu & Vigeland 2021; Charisi et al. 2022). Although PTAs can realistically localise the CGW source only within several deg² (Sesana & Vecchio 2010; Goldstein et al. 2018) they must be originated by particularly massive and/or nearby SMBHBs, making the search of their galaxy host viable (Goldstein et al. 2019; Petrov et al. 2024). Identifying the host and detecting accretion activity onto the binary will shed light about the physics governing the formation and (thermo)dynamical evolution of SMBHBs and the accretion discs surrounding them.

In this paper, we carry out the first detailed investigation of the interplay between stochastic GWBs and deterministic CGWs in realistic PTA datasets, assessing the performance of currently used analysis pipelines when both signals are present in the data. We started by building state-of-the-art mock datasets capturing the complexity of real data and injecting them with as realistic as possible GW signals, mimicking the ones found in the recent EPTA analyses (EPTA Collaboration et al. 2023a,b; Antoniadis et al. 2024). Realistic GWBs were injected as the incoherent sum of individual waves consistent with astrophysical observations (Rosado et al. 2015), and CGW are modelled according to the candidate identified in Antoniadis et al. (2024). The datasets were then analysed with the standard models implemented in the Enhanced Numerical Toolbox Enabling a Robust PulsaR Inference SuitE (ENTERPRISE, Ellis et al. 2019; Taylor et al. 2021), which models the GWB as a common red process characterized by a single power law spectrum. We also tested CGW models that include only the Earth term or the full Earth plus pulsar components.

The paper is organized as follows. In Sect. 2, we describe in detail the simulations: the properties of the adopted datasets, the injected signals, and the model used for the recovery. In Sect. 3, we present the main results of our analysis. Using the methods described in the previous section, we carried out several investigations. We first assess the importance of including the pulsar term in the CGW model (Sect. 3.1). We then study the degeneracy of signal recovery on datasets containing only a GWB or a CGW (Sect. 3.2), showing that datasets injected only with a CGW can produce evidence of a GWB (Sect. 3.2.1) and vice versa (Sect. 3.2.2). We also show that a joint search is able to pin down the correct nature of the signal (Sect. 3.2.3). The performance of a combined analysis of a GWB and a CGW on datasets containing both signals is then evaluated in Sect. 3.3. We discuss our main findings in Sect. 4. Throughout the paper, we assume a Λ cold dark matter cosmology with parameters (H₀, Ω_m, Ω_Λ) = (70 kms⁻¹ Mpc⁻¹, 0.3, 0.7).

2. Methods

In this section we describe both the simulation structure, along with the models and methods used in the analysis.

2.1. Simulations

Our simulations are based on the second data release of EPTA EPTA Collaboration et al. (2023b). In particular, we have reproduced the DR2full dataset: 25 pulsars with a maximum observation time of 24.8 yr. The main properties of the real times of arrival (ToAs) for each pulsar are reproduced: the observation time, the number of ToAs, the timing model (TM), the multiband measurements (note: the frequency coverage is identical to that of the EPTA DR2full pulsars), and the noise levels. The uneven sampling of the ToAs is also reproduced, but large gaps (which are common in the real data streams) are not present in the simulations (see Fig. 1 for a comparison between real and simulated residuals). All the properties of the arrival times were injected in the dataset using the libstempo package (Vallisneri 2020).

Fig. 1. Pulsar PSR J1022+1001. Top: Simulated data. Bottom: Real data from EPTA DR2full. The properties that are exactly reproduced are the time span, number of observations, frequency coverage, and the noise levels. The main difference between the two datasets is in the distribution of the ToAs; while in our simulations, ToAs are unevenly sampled, they are still uniformly distributed over the observational timespan, meaning that sparse observations and long gaps are absent (as can be clearly seen in a visual comparison of the two panels).

2.1.1. Noise injection

The libstempo package is also used to inject the three components of the noise budget: white noise (WN), achromatic (independent of the observation radio frequency) red noise (RN) and chromatic red noise (due to dispersion measure variations and customarily referred to as DM).

To model the white noise, the ToAs are sampled from a Gaussian distribution centred at the value predicted by the TM and with a width given by:

$$\begin{array}{c}\hfill {\sigma }_i=\sqrt{{\mathrm{EFAC}}_i^2{\sigma }_{\mathrm{ToA}}^2+{\mathrm{EQUAD}}_i^2},\end{array}$$(1)

where σ_ToA is the error bar associated with each ToA according to the template-fitting errors obtained in EPTA DR2, the EFAC takes into account for the ToA measurement errors and the EQUAD accounts for any putative extra source of white noise. The two parameters EFAC and EQUAD are specific for each observing backend. In the simulations, each pulsar’s backend, as well as the ToAs’ uncertainties, are identical those of DR2full; therefore, each pulsar has ToAs taken from ∼12 different backends, which means 12 different couples of values for EFAC and EQUAD. The values used for EFAC and EQUAD are the maximum likelihood values obtained by the single pulsar noise analysis of EPTA DR2full (EPTA Collaboration et al. 2023b).

In addition to the white noise, achromatic and chromatic red noise contribute to the single-pulsar noise budget. These are time-correlated noise components modelled as a stationary Gaussian process with a power-law power spectral density (PSD) of the form

$$\begin{array}{c}\hfill P(f;A,\gamma )=\frac{A^2}{12{\pi }^2}(\frac{f}{{\mathrm{yr}}^{-1}}{)}^{-\gamma }{\mathrm{yr}}^3.\end{array}$$(2)

For achromatic red noise, in each ToA, the delay induced in the residuals is proportional to the red noise PSD,

$$\begin{array}{c}\hfill {\mathrm{\Delta }}_{\mathrm{RN}}\propto \sqrt{P(f;A_{\mathrm{RN}},{\gamma }_{\mathrm{RN}})}.\end{array}$$(3)

This noise component accounts for the long-term variability of the pulsar spin along with other sources, such as unmodelled objects orbiting the pulsar. For chromatic red noise, on the other hand, the induced delay also depends on the observing radio frequency because it is due to the interaction of the radio signal with the ionised interstellar medium (IISM), the interplanetary medium in the Solar System, and the Earth’s ionosphere. This effect is taken into account during the observations and inside the timing model which considers its value at a reference epoch together with its first and second derivatives. However, the inhomogeneous and turbulent nature of the IISM also induces stochastic variations in the DM value, which are modelled as chromatic red noise. The residual DM delay is modelled as

$$\begin{array}{c}\hfill {\mathrm{\Delta }}_{\mathrm{DM}}\propto \sqrt{P(f;A_{\mathrm{DM}},{\gamma }_{\mathrm{DM}})}\times {\nu }^{-2},\end{array}$$(4)

where ν is the observing radio frequency. The values of A_RN, γ_RN, A_DM, and γ_DM values used in the simulations are reported in Table 1, together with the number of Fourier components used in the modelling of these processes. These values, namely, the amplitude, slope, and number of Fourier components, correspond to maximum likelihood values obtained from the single pulsar noise analysis of EPTA DR2full. An example of simulated data compared to the real EPTADR2 dataset is shown in Fig. 1.

2.1.2. Gravitational wave signal injection

Having specified the properties of the adopted PTA, we now turn to the description of the injected GW signals. The signal produced by a SMBHB is composed by two terms, customarily referred to as Earth term (ET) and pulsar term (PT), which are respectively associated with the GW strain hitting the Earth at the ToA of the considered pulse and the strain hitting the pulsar at the time of emission of the same pulse.

Throughout the paper we will always assume circular SMBHBs, evolving under the effect of GW backreaction only (i.e. we neglect any possible coupling with the local stellar or gaseous environment). Under these assumptions, the delays induced in the ToAs at an epoch t are fully defined by 8 parameters: the chirp mass ℳ¹, the observed frequency f_gw, the redshift z, the inclination angle ı, the sky location (θ, ϕ) – respectively polar and azimuthal angle –, the initial phase ϕ₀ and the polarization angle ψ. Following Ellis (2013), we model them as:

$$\begin{array}{c}\hfill s(t,\widehat{\mathrm{\Omega }})=F^{+}(\widehat{\mathrm{\Omega }})\mathrm{\Delta }{s}_{+}(t)+F^{\times }(\widehat{\mathrm{\Omega }})\mathrm{\Delta }{s}_{\times }(t),\end{array}$$(5)

where

$$\begin{array}{c}\hfill \mathrm{\Delta }{s}_A(t)=s_A(t_p)-s_A(t_e),\end{array}$$(6)

the functions $F^A(\widehat{\mathrm{\Omega }})$ – with A indexing the polarization and $\widehat{\mathrm{\Omega }}$ being the unit vector pointing from the GW source to the Solar System Barycentre (SSB) – are the antenna pattern functions that depend on the relative position of the pulsar and the source (see Ellis 2013, for full details). s_A(t_p) and s_A(t_e) are the PT and ET, respectively. They both take the form:

$$\begin{array}{cc}\hfill s_{+}(t)& =\frac{{\mathcal{M}}^{5/3}}{d_L\omega {(t)}^{1/3}}(-\text{sin}(2(\mathrm{\Phi }(t)-{\mathrm{\Phi }}_0))(1+{\text{cos}}^2ı)\text{cos}(2\psi )\hfill \\ \hfill & -2\text{cos}(2(\mathrm{\Phi }(t)-{\mathrm{\Phi }}_0))\text{cos}ı\text{sin}(2\psi ))\hfill \\ \hfill s_{\times }(t)& =\frac{{\mathcal{M}}^{5/3}}{d_L\omega {(t)}^{1/3}}(-\text{sin}(2(\mathrm{\Phi }(t)-{\mathrm{\Phi }}_0))(1+{\text{cos}}^2ı)\text{sin}(2\psi )\hfill \\ \hfill & +2\text{cos}(2(\mathrm{\Phi }(t)-{\mathrm{\Phi }}_0))\text{cos}ı\text{cos}(2\psi )),\hfill \end{array}$$(7)

where

$$\begin{array}{c}\hfill \mathrm{\Phi }(t)={\mathrm{\Phi }}_0+\frac{1}{32{\mathcal{M}}^{5/3}}({\omega }_0^{-5/3}-\omega {(t)}^{-5/3}),\end{array}$$(8)

and

$$\begin{array}{c}\hfill \omega (t)=({\omega }_0^{-8/3}-\frac{256}{5}{\mathcal{M}}^{5/3}t{)}^{-3/8}.\end{array}$$(9)

The only difference between PT and ET is that while the latter is the delay induced at the time of observation on Earth t_e, the former has to be evaluated at the time of emission of the radio pulse at the pulsar $t_p=t_e-L(1+\widehat{\mathrm{\Omega }}·\widehat{p})$, where L is the pulsar distance and $\widehat{p}$ is the pulsar position in the sky. We inject in our data two different GW components: a stochastic GWB and a resolvable CGW. Since the evolution of the signal over the observation time is small, the phase approximation is used both in the injection and in the recovery: ω_Earth = ω₀ at each epoch t.

Stochastic GWB. The GWB is injected as a realistic superposition of GWs coming from a cosmic population of SMBHBs. The injection pipeline starts from a catalogue of SMBHBs, computes the delays induced by each binary in each pulsar ToA according to Eq. (5) and then sums all the delays in the time domain. The distribution of the sources’ ℳ, f_gw, z and ı determines the observed power spectrum generated by the SMBHBs population: the binned total characteristic strain, where the frequency bin width is Δf = 1/T with T the time span of the pulsar array, can be evaluated as (Rosado et al. 2015)

$$\begin{array}{c}\hfill h_c^2(f_i)={\displaystyle \sum _{j\in \mathrm{\Delta }{f}_i}}\frac{h_j^2(f_r)f_r}{\mathrm{\Delta }{f}_i},\end{array}$$(10)

where f_i is the central frequency of the bin Δf_i, and the sum runs over all the binaries for which the observed frequency f = f_r/(1 + z)∈Δf_i, where f_r is the GW frequency in the binary rest frame. The strain of each GW signal is given by

$$\begin{array}{c}\hfill h(f_r)=2\sqrt{\frac{1}{2}(a{(ı)}^2+b{(ı)}^2)}\frac{{(G\mathcal{M})}^{5/3}{(\pi f_r)}^{2/3}}{c^{4}d(z)},\end{array}$$(11)

where a(ı) = 1 + cos²ı, b(ı) = − 2cosı and d(z) is the comoving distance of the binary, which is computed as

$$\begin{array}{c}\hfill d(z)=\frac{c}{H_0}{\displaystyle {\int }_0^z}\frac{1}{\sqrt{{\mathrm{\Omega }}_m{(1+z)}^3+{\mathrm{\Omega }}_{\mathrm{\Lambda }}}}dz.\end{array}$$(12)

For the population of circular GW driven SMBHBs, assumed here, the characteristic strain given by Eq. (10) can be approximated as a power law. It is therefore customary to describe the signal as

$$\begin{array}{c}\hfill h_c(f)=A_{\mathrm{GW}}{\left(\frac{f}{{\mathrm{yr}}^{-1}}\right)}^{\alpha },\end{array}$$(13)

where α = (3 − γ)/2 = −2/3 being γ the PSD induced in the timing residual (cf. Eq. (2)).

The power spectrum $h_c^2(f_i)$ thus depends on the number of emitting sources per unit redshift, mass and frequency d³N/(dzdℳdf_r). A list of ∼120 k binaries is constructed by extracting ℳ, f_gw, z as Monte Carlo samples of a numerical distribution d³N/(dzdℳdf_r), obtained from the observation-based models described in Sesana (2013b) and Rosado et al. (2015). In those papers, the authors construct hundreds of thousands of SMBHB populations in agreement with current constraints on the galaxy merger rate, the relation between SMBHs and their hosts, the efficiency of SMBH coalescence and the accretion processes following galaxy mergers. From this set of models, we selected 100 populations producing a GWB amplitude A ∼ 2.4 × 10⁻¹⁵, consistent with the signal observed by EPTA (EPTA Collaboration et al. 2023a). The other parameters of the binaries are then randomly sampled such that sources are isotropically distributed in the sky – cos θ ∈ Uniform(−1, 1) and ϕ ∈ Uniform(0, 2π) –, their orbital angular momenta are uniformly distributed on the unit sphere – cos ι ∈ Uniform(−1, 1) –, the phase and the polarization angle are random – ϕ₀ ∈ Uniform(0, 2π), ψ ∈ Uniform(0, π).

Resolvable CGW. In some of the dataset, an individually resolvable CGW is injected. The waveform is again injected in each pulsar in the time domain residual according to Eq. (5). The only difference, compared to the SMBHBs contributing to the unresolved GWB is that the amplitude of the signal is large enough that the source can be identified above the pulsar noise and the unresolved GWB by a signal search employing a deterministic filter to match the signal (see the next section). We stress that we always inject both the ET and the PT for every SMBHB, regardless of whether it provides an unresolvable contribution to the GWB or it is a resolvable CGW.

Datasets. For the purpose of this investigation, seven different datasets were produced, with the exact same noise realisation, but different gravitational signals:

• Dataset injGWB01, injGWB02, and injGWB03: These three datasets contain a GWB-only. They are produced injecting three background selected from the 100 populations mentioned above. The three populations are chosen among those displaying a smooth behaviour in the lowest frequency bins, without particularly loud sources.

• Dataset injCGW5 and injCGW20: These two datasets only contain a single loud CGW-only. The injection procedure is similar to the GWB datasets, but the population file used contains only one binary. The two injected sources have Earth term GW frequencies of 5 nHz (injCGW5) and 20 nHz (injCGW20). The parameters of the 5 nHz source are chosen to approximately reproduce frequency and amplitude of the candidate CGW found in the EPTA DR2new, Antoniadis et al. (2024), which in fact has a frequency of ∼5 nHz and an amplitude of h ∼ 1.1 × 10⁻¹⁴. The source is then located in the Virgo cluster; this choice determines the sky location and distance, and consequently the chirp mass. Inclination, initial phase and polarisation angles are randomly extracted as in the GWB case. The source injected in injCGW20 has identical parameters except for the higher frequency (20 nHz instead of 5 nHz), thus resulting in a higher strain (h ∼ 2.7 × 10⁻¹⁴). The parameters of the two CGWs used are listed in Table A.1. The delays induced in pulsar PSR J1022+1001 by these two CGWs are shown in Fig. 2.

  Fig. 2. CGWs injected in Pulsar PSR J1022+1001. Top: Injected delay induced by the 5 nHz CGW. Bottom: Injected delay induced by the 20 nHz CGW. Earth term, pulsar term and total delay are shown.

• Dataset injGWB03+CGW5 and injGWB03+CGW20: These two datasets are injected with a GWB and a CGW. The GWB is in both cases the one resulting from population 03, while the CGWs are either at 5 nHz or 20 nHz.

2.2. Signal modelling and analysis methods

All the analyses are carried using the ENTERPRISE (Ellis et al. 2019; Taylor et al. 2021) package. The starting point of the analysis is the likelihood function of the timing residuals δt

$$\begin{array}{c}\hfill P(\delta \mathbf{t}|\xi ,\zeta )=\frac{\text{exp}[-\frac{1}{2}{(\delta \mathbf{t}-\mathbf{M}ϵ-\mathbf{s})}^T{\mathbf{C}}^{-1}(\delta \mathbf{t}-\mathbf{M}ϵ-\mathbf{s})]}{\sqrt{{(2\pi )}^l\text{det}\mathbf{C}}}·\end{array}$$(14)

The term Mϵ is the timing model error expressed as a linear function of the offset from the nominal timing model parameters, ϵ. ξ are the timing model parameters and ζ are the stochastic components parameters. The likelihood function sampled to search for the gravitational signal is marginalized over the timing model parameters as described in Chalumeau et al. (2022). C = C(ζ) is the covariance matrix that contains all the contributions from stochastic processes, which are:

• WN: modelled according to Eq. (1). All the white noise parameters are fixed at the injected values in the GW searches;

• RN: modelled as a Gaussian process with amplitude of the Fourier components given by Eq. (2) (Lentati et al. 2014);

• DM: modelled as a Gaussian process with amplitude of the Fourier components given by Eq. (2) and amplitude depending on the observed radio frequency (see Eq. (4)). As for RN, the DM is not present in all the pulsars and, if present, the number of Fourier bins used in the Gaussian process depends on the pulsar and it is equal to the number of components used in the injection. For a summary of the noise model characteristics, see Table 1.

• Common red noise (CRN), if searched for: A red noise process which is common to all the pulsars can be accounted for in the covariance matrix. If included in the analysis, it is also modelled as a Gaussian process. The amplitude of the components are modelled in two ways: as a power law function of the frequency and, therefore, the parameters are the slope γ_CRN and the amplitude A_CRN; or as a binned spectrum, and thus the parameters are the amplitudes of the common process in each frequency bin f_i. In both cases, the number of components used in the Gaussian process is 30.

Although we injected the GWB as a superposition of unresolved SMBHBs, we followed the standard PTA analysis practice and modelled it in our likelihood function as a CRN characterised by an amplitude and a spectral index, as in Eq. (2), which can be directly translated in a characteristic strain of the form given by Eq. (13). We note that depending on the dataset, we use two different flavours of CRN. In most searches, the GWB is modelled as a common uncorrelated red noise (CURN). This means that inter-pulsar correlations induced by the GWB are neglected and the GWB posterior distribution is then reconstructed a posteriori with the re-weighting technique Hourihane et al. (2023), when possible. The advantage of this model is that it is much faster since it keeps the correlation matrix in the likelihood function diagonal and it naturally provides the Bayes factor of the reconstructed model versus the sampled model, besides the posterior distribution of both models. However, the re-weighting technique is not always efficient and (as we explain later in this work), using a CURN to model a GWB is not always appropriate; therefore, when needed, we use the HD-correlated CRN model which we referred to as GWB model. In this case, GWB induced inter-pulsar spatial correlation are also modelled, according to the HD curve (Hellings & Downs 1983)

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{\mathrm{ab}}=\frac{3}{2}\frac{1-cos{\theta }_{\mathrm{ab}}}{2}ln\left(\frac{1-cos{\theta }_{\mathrm{ab}}}{2}\right)-\frac{1}{4}\frac{1-cos{\theta }_{\mathrm{ab}}}{2}+\frac{1}{2},\end{array}$$(15)

where a and b are indexes identifying two pulsars in the array, and θ_ab is their angular separation in the sky. We note that searching for an HD-correlated CRN is much more computationally expensive, since the correlation matrix is not diagonal anymore.

The last term in the likelihood function is s, which is the delay induced by an individually resolvable CGW, and is analytically described by Eqs. (5)–(9) which are implemented in the CGW model of ENTERPRISE. An approximation that can be adopted is to consider in the CGW model the ET only, neglecting the PT. This is computationally very convenient, since the former only depends on the 8 source parameters. Conversely, the latter also depends on the distance of each pulsar in the array. Since those distances are generally known with large uncertainties (of the order of 20%), they must be considered as additional free parameters in the search. Therefore, the complete model has at least one more parameter per each pulsar with respect to the ET only model. Moreover, as prescribed in Ellis (2013), in order to ease the search of the best parameters, one extra parameter is added per each pulsar, which is an initial phase at the pulsar, Φ_p. We assess the performance of an ET only search on realistic signals including ET and PT in Sect. 3.1.

We have mainly adopted Bayesian techniques in our analysis, although we also exploit the frequentist optimal statistics, when useful. The priors used in the search are listed in Table B.1 and are all uninformative, except for the pulsar distances’ priors that exploit information obtained from independent electromagnetic observations. The prior distributions have a width which is taken as the measurement error listed in Verbiest et al. (2012); for the pulsars that are not included in this paper, the relative error is fixed at 20%. The mean of the prior is chosen such that the pulsar distance used for the injection is within 1σ of the prior distribution for each pulsar. Parameter estimation is then performed by sampling the posterior distribution with the Markov chain Monte Carlo (MCMC) sampler PTMCMCSampler, Ellis & van Haasteren (2017). For the 25 pulsar PTA we employ, the dimensionality of the posterior distribution is 62 (2 parameters (slope and amplitude) for each red noise and dispersion measurement process included in the model) plus the GW model dimension (2 for the GWB power-law model, 8 or 58 for an ET only or a full ET + PT CGW model search respectively). Model selection is done by computing the Bayes factor, either with the re-weighting technique, when it is efficient, or with the product space method implemented in ENTERPRISE_EXTENSIONS (for details on the method see Hee et al. 2015; Taylor et al. 2020, for how it is implemented in PTA).

3. Results

Using the methods described in the previous section, we carried out several investigations. In Sect. 3.1 we assess the importance of modelling a CGW including both the Earth and pulsar terms. In Sect. 3.2, we investigate the degeneracy between the analysis results obtained with a GWB-only dataset and with a CGW-only dataset. In Sect. 3.2.1 it is shown how CGW-only datasets can produce strong evidence in favour of an HD-correlated GWB. Vice versa, in Sect. 3.2.2 it is shown how a GWB-only dataset can be misinterpreted as a continuous wave. This degeneracy is solved in Sect. 3.2.3, by means of a combined analysis that searches simultaneously for a GWB and a CGW on both GWB-only and CGW-only datasets. Finally, in Sect. 3.3, we evaluate the performance of a combined analysis of a GWB and a CGW on datasets containing both signals.

3.1. Importance of the pulsar term in CGW modelling

As described in Sect. 2.2, the signal imprinted by a CGW in the ToAs is the sum of two terms: the Earth term and the pulsar term. The ET is common to all pulsars and only depends on the 8 parameters describing the source. Conversely, the PT depends on the distance of each pulsar and the phase of the GW signal at the pulsar’s location. This means that modelling of the full CGW ET + PT in a search performed on an array of 25 pulsars (like the one we are considering here) requires 50 extra free parameters, that must be added to the eight source parameters and to the 62 parameters of the individual pulsars noise components. It is therefore interesting to compare searches of CGW modelled as ET only and ET + PT on realistic datasets, to see whether the computational efficiency of the ET only search bears any significant drawback.

The two templates have been tested on both CGW-only datasets (injCGW5 and injCGW20), and the main results are shown in Fig. 3. From the figure, it is clear that the performance of the ET only template depends on the CGW frequency. In fact, when applied to the injCGW20 dataset (right panel of Fig. 3), the ET only model performs well on all parameters that it is sensitive to, i.e. all source parameters except the chirp mass and the luminosity distance. Conversely, the same model applied to the injCGW5 dataset results in some non-negligible biases in the recovery. As it can be seen from the left panel of Fig. 3, the injected values for the frequency and sky location are not compatible with the 90% C.I. of the recovered posterior distributions.

Fig. 3. Marginalised posterior distributions of the 8 source parameters (6 in the corner plot plus the 2 sky coordinates in the sky map) and of the CGW amplitude (obtained as combination of distance, frequency and chirp mass posterior distributions, stand-alone panel) for the CGW injected in the dastaset injCGW5 (left panel) and injCGW20 (right panel). The signal is searched with a CGW-only model with ET-only (orange distributions) and with the complete ET + PT template (purple distributions). Injected values are shown as black solid lines.

The different performance of the ET only template in the two cases is due to the different degree of separability of the ET and PT terms for the two injected sources. For the injected system chirp mass, ℳ = 1.16 ⋅ 10⁹ M_⊙, the average frequency difference between ET and PT is $\overline{\mathrm{\Delta }f}\sim 3.8\times {10}^{-6}$ nHz yr${}^{-1}\overline{\mathrm{\Delta }t}\sim 1.0\times {10}^{-2}\mathrm{nHz}\sim 0.02\mathrm{\Delta }{f}_{\mathrm{bin}}$ for the 5 nHz CGW and $\overline{\mathrm{\Delta }f}\sim 1.4\times {10}^{-3}$ nHz yr${}^{-1}\overline{\mathrm{\Delta }t}\sim 4.0\mathrm{nHz}\sim 3.2\mathrm{\Delta }{f}_{\mathrm{bin}}$ for the 20 nHz CGW. Therefore, while in the 20 nHz CGW, ET and PT falls at very distinct frequencies, they blend together in the same frequency bin for the 5 nHz CGW. In the former case, the ET only template only picks the common power produced by the ET at its frequency, disregarding the PT term; while in the latter, the ET only template is trying to fit a superposition of ET and PT essentially overlapping at the same frequency. This concept is demonstrated in Fig. 4, where the reconstructed waveforms are compared to the injected ones. In the 20 nHz case (injCGW20 dataset), the ET only template almost perfectly fits the ET of the full injected waveform, while in the 5 nHz case (injCGW5), the ET only template absorbs power from both ET and PT and the reconstructed waveforms sits, for many pulsars, somewhat in between the ET and the full injected signal and is also way more uncertain as shown by the much larger credible region.

Fig. 4. Signal recovered with the ET-only template applied to the 5 nHz CGW (dataset injCGW5, top panel) and to the 20 nHz CGW (data injCGW20, bottom panel). Median and 90% credible interval of the signal recovery are shown for pulsar PSR J1022+1001 (top) and PSR J1713+4707 (bottom). Solid and dashed lines are respectively the injected total and the ET-only signals.

Conversely, the full ET + PT template is able to correctly reconstruct the signal and, because it is sensitive to the frequency evolution, it also allows us to constrain both the source chirp mass and the luminosity distance. The performance is very different for the 5 nHz CGW and the 20 nHz CGW: for the former the chirp mass recovery is $\mathcal{M}=0.{52}_{-0.48}^{+9.98}\times {10}^8{M}_{\odot }$, while for the latter it is $\mathcal{M}=1.{01}_{-0.06}^{+0.07}\times {10}^9{M}_{\odot }$, where we give the median and the 68% credible interval. The superior recovery of the chirp mass in the 20 nHz case is due to the much larger frequency evolution of the source between ET and PT, allowing us to place a much tighter constrain on ḟ and thus on ℳ. Another consequence of the large frequency evolution of the 20 nHz CGW is that the model is sensitive to the pulsar distances. This is shown in Fig. 5 where the reconstructed waveform (left) and the posterior distribution of the pulsar distance (right) are shown for PSR J1713+0747 for the injCGW20 (top row) and injCGW5 (bottom row) datasets. In the first case, the posterior distribution on the distance is narrower than the prior and centred at the injected value: this is allowed thanks to an appropriate modelling of the frequency evolution, which is reflected in the excellent reconstruction of the waveform. Conversely, for the 5 nHz case, the posterior distribution of the pulsar distance is almost indistinguishable from the prior, but signal reconstruction still matches the injection, as the frequency evolution is too small for the search to be sensitive to it. We note that the GW signal reconstruction and distance determination is not always as good for the 20 nHz CGW. This is shown in the middle row of Fig. 5 for pulsar PSR J1918-0642. In this case the posterior distribution resembles the prior, and while the reconstructed signal does not strongly disagree with the injection, the median fails to capture the signal shape accurately. In fact, only about 13 pulsars, whose evolution is accurately reconstructed, significantly contribute to the chirp mass recovery.

Fig. 5. Injected CGW residuals (black solid line) and the reconstructed residuals with median and 90% C.I. (solid line and shaded area) are shown in the first column. The second column shows the pulsar distance prior distribution (dashed line) and recovered posterior distribution (histogram) and the injected value of the distance (solid black line). The first two rows correspond to the results obtained by analysing injCGW20 dataset and are shown for pulsar PSR J1713+0747 (first row) and pulsar PSR J1918-0642 (second row); the last row shows the results obtained with injCGW5 dataset displayed for pulsar PSR J1713+0747.

In addition to obtaining a better and unbiased recovery of the CGW parameters, there is another rationale for preferring the complete CGW template over the ET-only template: in a joint CGW + GWB search, the contribution to the CGW signal due to the pulsar term, if not taken into account by the CGW model, is absorbed by the GWB model, since it is a common red noise process. It will therefore add extra spurious power to the GWB, biasing the recovery of its parameters. An example of this behaviour is shown in Fig. 6, where the GWB + CGW search is applied to the dataset injCGW5: when the CGW is searched with the full template, the posterior distribution of the GWB parameters is unconstrained, when the search is done with the ET-only template, the posterior distribution of the GWB parameters remains well constrained because it absorbs all the power coming from the unmodelled pulsar term.

Fig. 6. Posterior distribution of the GWB parameters obtained from dataset injCGW5 with the GWB + CGW ET-only (purple) and GWB + CGW ET + PT (orange) template.

We therefore conclude that while the ET-only template can be used to search for possible CGW candidates and provides a reasonable estimate of the source parameters (see also Charisi et al. 2024), the full ET + PT model must be preferred, especially in a joint search with a stochastic GWB. Despite the much lower computational efficiency, we therefore chose to always model our CGWs as ET + PT according to the procedures described below.

3.2. Degeneracy between GWB and CGW

In this section, searches for GWB, CGW and combined searches of the two signals together are applied to the three GWB realisations (injGWB01, injGWB02, injGWB03) and to the two CGW realisations (injCGW5, injCGW20). The aim is to study the degeneracy between the two signals and how it can be solved.

3.2.1. CGW misinterpreted as a GWB

We started by performing a GWB analysis on the three injGWB and the two injCGW datasets. To this end, we made use of both Bayesian and frequentist techniques, following the procedure below.

• A power-law-shaped and HD-correlated stochastic process is assumed for the GWB, then the slope γ_GWB and the amplitude log₁₀A_GWB of the powerlaw process are sampled, together with the noise parameters to obtain the posterior distribution;

• The Bayes factor, ${\mathcal{B}}_{\mathrm{CURN}}^{\mathrm{GWB}}$, is evaluated to compare the evidence for an HD-correlated and an uncorrelated red noise process. The Bayes factor is computed either with the re-weighting method or with the product space when re-weighting is not efficient;

• The averaged angular correlation between the pulsars is estimated with the optimal statistics method averaging over ten angular bins, each one containing 30 pulsar pairs. A power-law-shaped background with variable γ is assumed;

• The S/N of the HD cross-correlation marginalized over the noise parameters is obtained with the optimal statistics method, assuming a powerlaw-shaped background with variable γ.

The results of the power-law spectrum parameter estimation are listed in the last two columns of Table 2 for all five datasets and two examples of the recovered GWB parameters’ posterior distributions are shown in Fig. 7 for datasets injGWB02 and injCGW5. In general, the GWB search returns well constrained posteriors regardless of the actual content of the dataset, which in practice means that a CGW can easily be detected as (and mistaken for) a GWB. However, the signal recovered from the two CGW-only datasets (injCGW5 and injCGW20) is higher and flatter than the one obtained from the background datasets. This can be clearly seen both from the A and γ values of Table 2 and from the 2D contours in Fig. 7. This is expected since the extra power due to the two single sources is at frequencies which are higher than the loudest GWB frequency bins that dominate the powerlaw spectrum recovery in the GWB-only datasets.

Fig. 7. Posterior distribution of the powerlaw GWB parameters. Purple empty posterior is obtained from dataset injGWB02, orange full posterior is obtained from dataset injCGW05. The solid line is the ‘nominal’ GWB injection (γ ∼ 13/3, log10AGWB ∼ = − 14.66).

The Bayes factors reported in Table 2 show that an HD-correlated common process (i.e. a GWB) is generally strongly favoured regardless of the content of the data. In fact, dataset injCGW5 produces a Bayes factor ${\mathcal{B}}_{\mathrm{CURN}}^{\mathrm{GWB}}$ of ∼10 and HD signal-to-noise ratio (S/N) of 3.0${}_{-0.5}^{+0.5}$ (fixing γ to 3), while injCGW20 gives ${\mathcal{B}}_{\mathrm{CURN}}^{\mathrm{GWB}}$ of ∼60 and S/N = $3.2_{-0.6}^{+0.6}$ (fixing γ to 2). Interestingly, injGWB03 does not favour an HD correlated process, displaying ${\mathcal{B}}_{\mathrm{CURN}}^{\mathrm{GWB}}$ of ∼0.15. This is consistent with Valtolina et al. (2024), who found that, due to the stochastic nature of the signal, individual realisations of the same SMBHB population can result in Bayes factor that differ by several orders of magnitude.

Although the Bayes factors obtained on the injCGW5 and injCGW20 are compatible with a GWB detection, the resulting HD S/N is significant only in a small range of values of γ that are consistent with the posterior distributions obtained by the Bayesian search, as shown in Fig. 8. By fixing γ to the nominal value of 13/3, the HD S/N drops to 1.3${}_{-0.4}^{0.4}$ and to 1.1${}_{-0.9}^{1.0}$ for the injCGW5 and injCGW20 datasets respectively. Conversely, the injGWB01 and injGWB02 datasets (also shown in Fig. 8) result in distributions of the HD S/N that are consistently high for all values of γ ≥ 4. Moreover, at higher values of γ, the CGW-only dataset displays increasing evidence for a dipolar-correlated signal whose S/N surpasses that of the HD. This behaviour is not observed in the GWB-only datasets, for which the dipole S/N always stays lower than the HD S/N. Finally, two examples of HD reconstruction are shown in Fig. 9 for datasets injGWB02 and injCGW5. While the quadrupolar pattern can be seen in both datasets, the scatter around the expected relation appears to be much larger when a CGW only is present in the data, consistent with Cornish & Sesana (2013) and Allen (2023). Moreover, the HD correlation for the injCGW5 dataset is reconstructed assuming γ = 3, while it tends to disappear for higher values of γ.

Fig. 8. Optimal statistics S/N for an HD-correlated (empty violins) and dipole-correlated (filled violins) process as a function of the assumed γ power law index. The two upper panels show in purple the results for datasets injGWB01 (top) and injGWB02 (bottom), while the two lower panels show in orange the results for datasets injCGW5 (top) and injCGW20 (bottom).

Fig. 9. Overlap reduction function obtained with the optimal statistic. Purple circular and orange squared points indicate the results for the injGWB02 and injCGW5 datasets respectively. The correlation coefficients for each pair of pulsars are grouped in 10 bins of 30 pulsar pairs each, the points represent the weighted average in each bin. The solid line shows the expectation value of the HD correlation. Note that the purple rightmost point lie below the orange one.

3.2.2. GWB misinterpreted as a CGW

We went on to search for a CGW in the same five datasets described in the previous section. In the search, we estimated the posterior distributions of all the noise parameters (except for the WN, which is fixed) and the CGW parameters (eight source parameters, plus 50 additional parameters for the pulsars’ distances and phases). We then used hypermodelling to evaluate the Bayes factor between a model containing a CGW versus a model featuring intrinsic pulsar noise only.

We can see from Table 3 that in this case as well, regardless of whether the dataset contains a CGW or a GWB, a CGW is recovered with high confidence, with a Bayes factor of ${\mathcal{B}}_{\mathrm{PSRN}}^{\mathrm{CGW}}>1000$ in all cases. Table 3 also shows that the CGW frequency and amplitude are fairly well constrained. The performance of the CGW search on the injCGW datasets has already been assessed in Sect. 3.1. As shown in Fig. 3 when a full ET + PT template is used for the CGW, the model can correctly recover the injected signals and estimate their parameters, although with differences depending on the frequency of the source. Here it is particularly interesting to check the CGW recovery on the three injGWB datasets, which is shown in Fig. 10: all source parameters are fairly well constrained, including the sky location, which is different in each dataset. These parameters, however, do not correspond to any specific source present in the datasets, as can be seen in the three panels of Fig. 11. The recovered CGW frequency always corresponds to the frequency bin that is the most constrained by the HD-correlated free spectrum search; moreover, the amplitude recovered is always compatible with the GWB amplitude in that frequency bin. Thus, it can be concluded that the CGW parameters are constrained because the CGW model absorbs part of the HD-correlated common red noise produced by the incoherent superposition of all the SMBHBs composing the GWB.

Fig. 10. Posterior distributions of the CGW parameters, inclination angle, chirp mass, luminosity distance, frequency, initial phase, and polarization angle, recovered with a CGW-only model (considering both Earth and pulsar term) applied to the three GWB only datasets. The inset sky map shows the reconstruction of the (fake) CGW sky location. The black stars in the sky map represent the positions of the 25 EPTA pulsars, while the two red dots are the Virgo and the Fornax clusters, there as a reference.

Fig. 11. Injected and recovered signals in the three injGWB datasets. In each panel, the green stars represent the contributions of each binary of the injected population to the GWB, the black line is the total power spectral density per frequency bin. The orange violin plot is the HD-correlated free spectrum recovery and the purple violin is the posterior distribution of the recovered CGW amplitude at the location of the recovered CGW frequency.

Finally, we note that the CGW amplitudes recovered from the three injGWB datasets are respectively $1.2_{-0.7}^{+1.5}\times {10}^{-14}$, $1.4_{-0.5}^{+0.6}\times {10}^{-14}$ and $1.8_{-0.4}^{+0.5}\times {10}^{-14}$ (median and 68% credible interval). These values are all well compatible with the amplitude of the DR2new CGW candidate, which is $1.0_{-2.6}^{+0.5}\times {10}^{-14}$ (median and 90% credible interval) (Antoniadis et al. 2024).

3.2.3. Breaking the degeneracy: combined search for a GWB and a CGW

Fig. 12. Posterior distributions obtained from GWB-only, CGW-only, CURN + CGW and GWB + CGW searches performed on the GWB-only dataset injGWB02 (left) and on the CGW-only dataset injCGW5 (right).

We now perform a joint CGW + GWB search on the five datasets considered so far. We perform two flavours of combined searches: (i) common uncorrelated red noise plus CGW (CURN + CGW), and (ii) HD-correlated common red noise plus CGW (GWB + CGW). The CGW is always modelled as the sum of ET and PT, and together with the signal parameters, all the RN and DM parameters are also searched over.

The results of these two searches are reported in Table 4 for all datasets, while the posterior distributions of selected parameters of the model components are shown in Fig. 12 for datasets injGWB02 and injCGW5. In this figure, the orange and light-purple distributions are obtained respectively from the GWB + CGW and CURN + CGW searches, and results of GWB-only and CGW-only searches are also shown for comparison. To see how searching for both type of signals simultaneously allows us to properly characterise the actual content of the data, let us discuss the injGWB and injCGW separately.

In the injGWB cases, when a common red noise component is added to the model alongside the CGW, the evidence for the latter strongly decreases. This can be seen both from the drop of the Bayes factors in the first three lines of Table 4, and from the posterior distributions of the CGW parameters in the left panel of Fig. 12. In Sect. 3.2.2, when searching for a CGW on a dataset containing only a GWB, we found ${\mathcal{B}}_{\mathrm{PSRN}}^{\mathrm{CGW}}>1000$. This is because the GWB produces a common noise that the model prefers to assign to the CGW signal (common to all pulsars) rather than to the individual pulsars’ RN components. When a CURN is included in the model, ${\mathcal{B}}_{\mathrm{CURN}}^{\mathrm{CURN}+\mathrm{CGW}}$ drops significantly but is still slightly larger than unity, suggesting that the presence of CGW component is still preferred. This can also be seen in the posterior distribution of the CGW amplitude shown in purple in the left panel of Fig. 12, which (despite featuring a long tail that extends to the lower end of the prior range) still shows a significant peak around log A ≈ −14.3. This behaviour is shared by all three injGWB datasets and is due to the fact that the GWB manifests itself primarily as a common red noise process, with a weaker, emerging ‘off-diagonal’ HD-correlated component. The CURN models it as a common red noise process, but is unable to account for the emerging HD correlation, which the model then assigns to the CGW. Conversely, when a GWB is added to the model, we get ${\mathcal{B}}_{\mathrm{GWB}}^{\mathrm{GWB}+\mathrm{CGW}}<1$ (first three lines of Table 4), strongly suggesting that the GWB component is now able to describe the full content of the data. This is also clear by looking at the CGW posterior parameters in the left panel of Fig. 12. The orange distributions show that both the frequency and amplitude posteriors are almost perfectly flat, with the latter displaying a cut-off, which we find to be at a similar value for all three realisations. From this, we can extract a 95% upper limit of the CGW amplitude which is ≈4.0 × 10⁻¹⁵, ≈3.8 × 10⁻¹⁵ and ≈4.6 × 10⁻¹⁵ for datasets injGWB01, injGWB02 and injGWB03 respectively. We note that these upper limits are lower than the GWB level, since the total GW content of the dataset is correctly accounted for by the GWB component. Finally, the recovered source is not localised in the sky, with the posterior distributions of the sky location almost uniform across the sphere, consistent with the prior.

In contrast to the CGW parameters, the GWB parameters’ posterior distributions are not significantly affected by the addition of the CGW component to the model: the GWB posterior distributions obtained from the GWB + CGW search and with the GWB-only search are well consistent with each other (see left panel of Fig. 12 for the injGWB02 dataset and Table 4 for all the datasets) confirming again that the GWB component of the model is correctly accounting for all the HD-correlated power. Only a small variation can be noticed in the CURN + CGW recovery of the background parameters: the posterior distributions are slightly wider with a longer low amplitude tail. This is because the CURN cannot absorb any HD correlated power, which is (spuriously) assigned to the CGW component.

Conversely, when the same analysis is applied to the injCGW datasets, the posterior distributions of the CGW parameters are always well constrained, both with the CURN + CGW and with the GWB + CGW search, as shown in the right panel of Fig. 12. Searching for CGW-only, CURN + CGW or GWB + CGW, in this case, gives consistent results for the CGW parameter, with no significant changes in the posterior distributions’ shape and width (see right panel of Fig. 12). Moreover, GWB + CGW is able to recover the sky location with a precision which is comparable to the one obtained with the CGW-only search and the model comparison strongly favours the GWB + CGW model over the GWB-only model (see Table 4). As expected, the GWB posterior distribution is instead completely unconstrained when the model applied in the analysis contains both a common red noise and a CGW: the power that generates the GWB signal recovered with the GWB-only search is now correctly completely absorbed by the CGW component (see Fig. 12).

3.2.4. Noise-only datasets

As a further check, we tested the analysis pipelines on a dataset containing only the noise and no GW signal, neither deterministic nor stochastic. The noise realisation is the one contained in all the datasets analysed. In this case, neither a GWB nor a CGW is detected. As for the GWB, we find a value of the Bayes factor ${\mathcal{B}}_{\mathrm{CURN}}^{\mathrm{HD}}\sim 1.5$ and the posterior distributions of A and γ are poorly constrained, with a 95% upper limit on the GWB amplitude of ∼4.47 × 10⁻¹⁵. As for the CGW, the Bayes factor ${\mathcal{B}}_{\mathrm{PSRN}}^{\mathrm{CGW}}$ is ∼0.5 and the CGW parameters are totally unconstrained, with the amplitude showing a 95% upper limit of ∼4.5 × 10⁻¹⁵. Although many more noise realisations should be tested to draw robust statistical conclusions, these tests show on fully realistic datasets that current PTA detection pipelines are unlikely to return a spurious detection in absence of an actual GW signal in the data.

3.3. Joint GWB + CGW analysis on combined CGW + CGW datasets

In this section we present the results of simultaneous searches for a GWB and a CGW, applied to datasets injGWB03+CGW5 and injGWB03+CGW20. We remind the reader that both datasets contain the same noise and GWB realisation from dataset injGWB03. In each dataset a loud CGW is superimposed to the GWB: in one case the CGW is at 5 nHz (injGWB03CGW5), in the second one it is at 20 nHz (injGWB03CGW20). This is done in order to see how the joint analysis performs in presence of a slowly evolving source (the 5 nHz CGW) and on a fast evolving one (the 20 nHz CGW), and how the two of them interact with the GWB.

3.3.1. Bias induced in the GWB recovery by the presence of a resolvable source in the SMBHBs population

As already discussed by Bécsy et al. (2023) and Valtolina et al. (2024), the presence of a particularly loud source among the systems contributing to the GWB, if not accounted for in the analysis, can significantly bias the recovery of the overall GWB signal. As a preliminary step in our investigation, we therefore perform a GWB only search on the injGWB03CGW5 and injGWB03CGW20 datasets. As expected, the analysis produces a biased recovery of the GWB parameters: the results are shown in the right panels of Figs. 13 and 14: in both cases, the recovered background is flatter and higher than the injected signal due to the extra power introduced by the loud source. In the two following sections, we demonstrate how this bias can be solved by performing a joint analysis of the GWB and the CGW components.

Fig. 13. CGW (left) and GWB (right) recovered by the joint search applied on dataset injGWB03CGW5. Left panel: Posterior distributions of the CGW source parameter using the full ET + PT template. The panel is organized as the ones in Fig. 3, with 6 source parameters shown in the triangle plot, the sky localization error shown in the sky map and the recovered signal amplitude h shown in the stand-alone inset. Injected values are marked by the black solid lines. Right panel: Posterior distributions of the GWB parameters. The GWB recovery from the joint search on the injGWB03+CGW5 dataset is shown in orange. For comparison, the recovery of a GWB only search on the injGWB03CGW5 is shown in red while the recovery of a GWB only search on the GWB-only injGWB03 dataset is shown in purple. Solid black lines represent the injected values. In both panels we mark for each parameter the median and 68% credible interval of the recovery.

3.3.2. Results for the injGWB03+CGW5 dataset

The results presented here are obtained by sampling the posterior distributions of the CURN + CGW model and then re-weighting them to obtain the posteriors of the GWB + CGW model. The results of the combined search applied to the dataset injGWB03+CGW5 provides good overall results, with the majority of the CGW parameters recovered within 2σ of the injected values.

The posterior distributions of the CGW parameters are shown in the left panel of Fig. 13, while the recovery of the background parameters is shown in the right panel. The recovery of the CGW parameters is successful: the frequency and the sky localization of the source are well constrained – the width of the posterior distributions is ≲10% of the prior width for all three parameters – and compatible within 68% with the injected values: the injected frequency is f = 5.02 nHz, recovered frequency is $f=5.{04}_{-0.04}^{+0.04}$ nHz; injected sky location is (12^h19^min50^s, 12.38°), recovered sky location is (${12}^{\mathrm{h}}{33}_{-19}^{+20\min }$, ${10}_{-5}^{+5}$ °). However, the 3D source location cannot be properly constrained, since the luminosity distance posterior distribution is extremely wide due to its strong covariance with the chirp mass. Both posterior distributions are roughly bimodal, with a narrow peak close to the injected value – in particular for the chirp mass – and a long, bumpy tail extending to low values of luminosity distance and chirp mass. We note that the posterior of the luminosity distance spans the interval (0.02, 24) Mpc within the 90% credible region; for such a nearby putative source, ancillary information coming from galaxy catalogues can effectively enable the identification of the source host (Petrov et al. 2024). The inclination angle, initial phase and polarisation angle are well constrained and compatible with the injected values.

Focussing now on the GWB parameters, the posterior distribution in the right panel of Fig. 14 shows how accounting for the resolvable source in the model solves the bias in the GWB parameters’ recovery. In fact, the posterior distribution obtained with the combined search is well consistent with the one obtained with the GWB-only search applied to the GWB-only dataset injGWB03: they peak at similar values and have comparable widths, showing that the parameters of the stochastic and of the deterministic components are not significantly correlated and the interplay between the two signals does not affect the search.

3.3.3. Results for the injGWB03+CGW20 dataset

The joint search applied to the injGWB03CGW20 dataset unveils some new interesting behaviour of the recovery pipeline. In fact, in this case, when a joint CURN + CGW (ET + PT) search is performed and the GWB is obtained with the re-weighting technique, the estimate of the GWB parameters is completely off, as demonstrated by the filled purple distributions shown in the right panel of Fig. 14. The result is indeed quite similar to a GWB search only (shown in red), meaning that a significant fraction of the CGW power is absorbed in the common background component and the CURN + CGW search followed by re-weighting fails in this case, producing strong biases in the results.

Fig. 14. Same as Fig. 13 but for dataset injGWB03+CGW5. The only difference is that in the right panel, we show the joing GWB + CGW analysis either modelling the GWB as a CURN and then using the re-weighting technique (filled purple) or by directly modelling the HD correlations in the GWB search (filled orange). All other distributions and markers have the same colour style and meaning as in Fig. 13.

The reason behind this failure is explained in the following: as shown in Sect. 3.1, due to the significant difference in the 20 nHz CGW between the Earth and the pulsar term GW frequency, in order to correctly account for both terms with the CGW model, the majority of the 50 PT parameters (i.e. 25 pulsar distances and 25 individual PT GW phases) have to be constrained. This is a difficult task for the sampler, which in the CURN + CGW case converges to an easier, alternative solution: the pulsar terms are all absorbed by the common uncorrelated red noise component, which is possible because they are indeed a common sinusoidal process that is uncorrelated, since the pulsars are all at different sky locations and distances from Earth. This justifies the much higher and flatter recovered background shown in the right panel of Fig. 14: the recovered amplitude is $1.{83}_{-0.26}^{+0.28}\times {10}^{-14}$, while the injected value is 2.4 × 10⁻¹⁵.

This interpretation is confirmed by the recovered CGW parameters (not shown). Most of them are biased and in particular the chirp mass, which determines the frequency evolution of the signal, is much lower than the injected value: the obtained value is $6.9_{-5.0}^{+18.0}\times {10}^6{\mathrm{M}}_{\odot }$, which corresponds to a 95% upper limit on the frequency derivative of 4.2 × 10⁻⁶ nHz⋅yr⁻¹, ∼300 times smaller than the ḟ of the injected SMBHB. This is consistent with the fact that the values of the CGW parameters are such that the recovered time domain signal matches the injected Earth term and completely neglects the pulsar term contribution. This is confirmed by the shape of the reconstructed CGW shown in Fig. 15. In the top panel we can see that in the CURN + CGW search the recovered CGW waveform for pulsar PSR J1713+0747 reproduces exactly the injected ET contribution. In fact, being HD-correlated, in presence of a CGW component in the model, the Earth term is not absorbed by the CURN component. Indeed, the recovered background is slightly lower than the one obtained by a GWB-only search on the same dataset (shown by the red distributions in the right panel of Fig. 14), because in that case both the Earth and pulsar terms are absorbed by the background component.

Fig. 15. Median and 90% credible region of the recovered CGW signal for pulsar PSR J1713+4707 from the joint analysis of dataset injGWB03+CGW5. The top panel shows the result obtained with a CURN + CGW (ET + PT) search and bottom panel shows the outcome of a GWB + CGW (ET + PT) search. Solid ad dashed line are respectively the total and the ET-only signal injected. In the first case, the CGW model neglects the pulsar term, resulting in a biased CGW parameters recovery (detail explanation in the main text).

We stress that this bias indeed corresponds to the maximum of the posterior distribution. This can be explained by the fact that the CURN model does not exactly match the injected common red noise, which is HD correlated, and therefore it prefers to absorb the pulsar term, since it is much cheaper to describe it as a two-parameter CURN signal, rather than to determine the 2 × 25 parameters describing each individual sinusoidal pulsar term.

In this dataset, the injected signal components are properly recovered only by a direct GWB + CGW (ET + PT) search, which is more computationally expensive due to the off-diagonal HD terms that must be computed when evaluating the correlation matrix in Eq. (14). In fact, in a GWB + CGW search the pulsar term is not absorbed by the HD-correlated common red noise component and both signals are properly reconstructed. In the right panel of Fig. 14 the GWB recovery obtained by this search is shown in filled purple: it agrees very well with that obtained with the GWB-only template applied to the GWB-only dataset (empty purple), although it is slightly larger, which can be explained by the fact that the GWB + CGW model contains 58 more parameters describing the single source. A tiny correlation with this high number of parameters can easily explain the less precise recovery of the background parameters.

The correct recovery of the GWB parameters is a direct consequence of the correct recovery of the CGW parameters. In the left panel of Fig. 14, the posterior distributions of the 8 source parameters are shown and in Fig. 15 the median and 90% credible region of the recovered CGW signal is shown for pulsar PSR J1713+4707. From the recovery, it can be seen how the model is able to account for the frequency evolution of the signal: the frequency recovered is $19.{99}_{0.18}^{+0.21}$ nHz (injected is 20 nHz) and the chirp mass $1.{14}_{0.10}^{+0.11}\times {10}^9{M}_{\odot }$ (injected is 1.16 × 10⁹ M_⊙). This allows for the luminosity distance to be much better constrained than in the low frequency case, with a median and 68% credible interval of $32.7_{-9.3}^{+15.1}$ Mpc. Since the injected value is 21.6 Mpc, the value obtained is a mild overestimate, which is probably due to the anti-correlation with the inclination angle and the initial phase, which are both slightly underestimated. The 2D sky location recovery is also improved with respect to the 5 nHz case, resulting in a declination angle of ${12}^{\mathrm{h}}{38}_{-16}^{+17\min }$ and a right ascension of $11.5_{-3.4}^{+3.0}$ °.

4. Discussion and conclusions

In this paper, we carry out the first comprehensive study of the interplay between stochastic GWBs and deterministic CGWs in realistic PTA datasets, assessing the performance of current analysis pipelines. To this end, we constructed mock datasets featuring all the complexity of EPTA DR2full (EPTA Collaboration et al. 2023c); we simulated up to 24.8 years of unevenly sampled, multi-frequency ToAs for the 25 pulsars of DR2full and injected WN, RN, and DM values that are consistent with the maximum likelihood values found in the DR2full single pulsar analysis (EPTA Collaboration et al. 2023b). The injected GW signals were also as realistic as possible, tailored to mimic the findings of the recent EPTA analyses (EPTA Collaboration et al. 2023a). A GWB with nominal amplitude of A_GW = 2.4 × 10⁻¹⁵ was injected as the incoherent sum of individual waves coming from observationally based cosmic populations of SMBHBs (EPTA Collaboration et al. 2024), while the deterministic CGW was modelled according to the candidate signal identified in the analysis of EPTA DR2new (Antoniadis et al. 2024). In all cases, SMBHBs were assumed to be circular and both the ET and PT of the GW-induced time delays were injected. These GW signals were searched with the standard models implemented in the ENTERPRISE analysis suite. The GWB was therefore modelled as a common red process characterized by a single power law spectrum, while the model CGW matched the one used for injecting the signal (i.e. circular, GW driven SMBHB), either including only ET or ET + PT. With this set-up we carried a number of investigations, whose results are summarised in the following.

We first compared the performance of ET-only and ET + PT templates on two datasets containing only a CGW either at 5 nHz (injCGW5) or at 20 nHz (injCGW20). In both cases, the ET + PT template was able to fully recover the signal, providing an unbiased estimate of all its parameters, including chirp mass, distance and sky location. However, things were not as straightforward when the ET only template was employed. In this case, the 5 nHz CGW dataset resulted in a biased recovery of the frequency and the sky location; whereas in the 20 nHz case, the recovery was unbiased. In fact, in the latter case Earth and pulsar term can be easily disentangled and the ET-only template correctly matches the ET component of the signal. Conversely, in the 5 nHz case, ET and PT fall in the same frequency bin and the ET template tries to fit a mixture of the two, leading to a bias in the signal reconstruction. In both cases, the ET-only model is not sensitive to the chirp mass and therefore to the luminosity distance, preventing 3D source localisation. We also stress that in both cases, modelling ET only leaves behind some unmodelled residual (e.g. the full PT in the 20 nHz case), which will likely be absorbed by the GWB template in a joint search, leading to a biased recovery of the background parameters. These are the reasons why the complete CGW template (ET + PT) must be preferred to the ET-only, even if this means having two times the number of pulsars extra parameters (50 in our case) and, thus, much greater computational costs for the sampling.

Therefore, we adopted the full ET + PT template for the CGW and turned to the investigation of datasets where only one of the two types of signal was injected (either a GWB or a CGW); these analyses were performed separately with one model at the time (eiter a GWB or a CGW). We analysed five datasets: three realisations of a realistic GWB (dataset injGWB01, injGWB02, and injGWB03) and two realizations of single resolvable sources (datasets injCGW5 and injCGW20). All datasets contained the same noise realisation. These searches showed that a GWB can easily be misinterpreted as a CGW and vice versa. In fact, all sources returned evidence in favour of either GW signal and well constrained recovery of the model parameters. However, when a combined search is performed, meaning that the model used for the recovery contains both an HD-correlated GWB and a CGW, this degeneracy is broken: when only a GWB is present in the data, the posterior distributions of the CGW parameters obtained with the combined search are unconstrained and the CGW amplitude shows an upper limit; on the contrary, when a CGW is the only signal injected, the combined search gives unconstrained posterior distributions for the GWB parameters. We find that no CGW evidence is absorbed by the HD-correlated background if the single source is present in the data. Therefore, the combined search is able to correctly recover the signal that is actually present in the data. Focussing on GWB-only datasets, we highlight that the CGW evidence is totally absorbed when an HD-correlated GWB is present in the model, while it does not completely disappear if only a CURN component is modelled. In this case, even if poorly constrained, the CGW parameters posterior distributions still peak at specific values. Therefore, modelling the common red noise as HD-correlated is necessary to totally absorb the spurious CGW. These results are similar to the ones obtained with the real EPTA DR2new data (Antoniadis et al. 2024): CURN + CGW searches peaking at a CGW amplitude of ∼10⁻¹⁴, with the peak disappearing once the background is modelled as HD-correlated and the Bayes factor ${\mathcal{B}}_{\mathrm{GWB}}^{\mathrm{GWB}+\mathrm{CGW}}$ is of the order of unity. Our simulations show that this behaviour arises when a CGW is indeed not present in the data, strengthening the conclusion reached in the EPTA DR2 analysis.

Next, we analysed datasets containing both a stochastic GWB and a CGW, either at 5 nHz (injGWB03+CGW5) or at 20 nHz (injGWB03+CGW20). The low-frequency CGW is an almost monochromatic source, while the high-frequency one is drifting and the signal frequency significantly changes between the PT and the ET. Performing a combined search on these types of datasets is particularly important because a single resolvable source, if not properly taken into account in the search, will induce a bias in the GWB parameters recovery (Bécsy et al. 2023; Valtolina et al. 2024). The results of the combined searches indicate that the bias in the GWB recovery induced by the presence of the single source is resolved when the CGW is correctly accounted for by the deterministic CGW model, so long as it contains both ET and PT. This result is obtained with both datasets, injGWB03CGW5 and injGWB03CGW20. However, in the injGWB03CGW5 case, the same unbiased result can be obtained with two methods: by performing a CURN + CGW search and then re-weighting to obtain the HD + CGW posterior distribution; or by directly sampling the HD + CGW posterior distribution. The first method is computationally advantageous, as the CURN model is considerably faster to compute than the HD one. In the injGWB03CGW20 case, however, the first method (re-weighting) is unable to recover the unbiased posterior distribution, as the CGW template solely accounts for the ET contribution, while the PT is more easily absorbed by the CURN. The posterior distribution on the CURN parameters is therefore strongly biased with respect to the GWB injection. Performing a direct sampling of the HD + CGW model, on the contrary, gives an unbiased recovery of the GWB parameters: the posterior distributions of the CURN and HD parameters are consequently disjoint, which makes the re-weighting technique extremely inefficient. Therefore, while computationally much more expensive, the only method that is able to converge to unbiased results is the direct HD + CGW search. As for the GWB, the CGW recovery is successful, with most of the parameters constrained and compatible with the injection.

Given the recent PTA results indicating the presence of a gravitational signal in the data, it is crucial to be able to distinguish between the two types of signal that can be produced by SMBHBs. This work represents a first thorough study of this problem. Despite the limited parameter exploration, which is constrained by the high computational cost of the analyses, the investigation decisively points towards the importance of performing combined searches. Indeed, in the cases analysed here, combined searches were able to both identify which signal was present in the data, without misinterpretation, and to guarantee an unbiased recovery of the GWB parameters in the presence of a loud individual source in the SMBHB population. Looking ahead, it will be important to extend this analysis to richer and more sensitive datasets, investigating the possibility of disentangling multiple sources, also in the presence of eccentricity or environmental coupling. In fact, the high amplitude of the currently observed signal implies that the future SKA might be able to resolve several tens of individual SMBHBs at S/N > 5 (Truant et al. 2025). Given the complexity and high dimensionality of these searches, it is of paramount importance to design more efficient search techniques (Bécsy 2024; Freedman & Vigeland 2024; Vallisneri et al. 2024) capable of handling multiple CGW signals. In fact, a full-sky search over the entire frequency spectrum, such as the one performed in this paper, might struggle to handle multiple sources due to the increasing number of parameters. A better strategy might be to gather some prior information about putative CGW candidates from anisotropy searches or from quick CGW identification algorithms (Bécsy & Cornish 2020). This would reduce the computational cost of the analysis and also prevent CGW blending. Having the ability to successfully identify and disentangle those signals will have invaluable benefits with respect to our understanding of SMBHB evolution and dynamics in the context of multi-messenger astrophysics.

Acknowledgments

I.F. thanks Stanislav Babak for his hospitality and mentoring at the APC, where this investigation was initiated. We thank the B-Massive group at Milano-Bicocca University for useful discussions and comments. A.C., A.S., G.M.S. acknowledge the financial support provided under the European Union’s H2020 ERC Consolidator Grant “Binary Massive Black Hole Astrophysics” (B Massive, Grant Agreement: 818691). A.C. acknowledges financial support provided under the European Union’s Horizon Europe ERC Starting Grant “A Gamma-ray Infrastructure to Advance Gravitational Wave Astrophysics” (GIGA).

References

Appendix A: CGW injection details

In this appendix, the values of the two injected continuous gravitational waves are listed.

Appendix B: MCMC details

In this appendix, the functional shape and the range of the prior distributions used in the MCMC sampling are listed.

All Tables

All Figures

Fig. 1. Pulsar PSR J1022+1001. Top: Simulated data. Bottom: Real data from EPTA DR2full. The properties that are exactly reproduced are the time span, number of observations, frequency coverage, and the noise levels. The main difference between the two datasets is in the distribution of the ToAs; while in our simulations, ToAs are unevenly sampled, they are still uniformly distributed over the observational timespan, meaning that sparse observations and long gaps are absent (as can be clearly seen in a visual comparison of the two panels).

In the text

	Fig. 2. CGWs injected in Pulsar PSR J1022+1001. Top: Injected delay induced by the 5 nHz CGW. Bottom: Injected delay induced by the 20 nHz CGW. Earth term, pulsar term and total delay are shown.
In the text

Fig. 3. Marginalised posterior distributions of the 8 source parameters (6 in the corner plot plus the 2 sky coordinates in the sky map) and of the CGW amplitude (obtained as combination of distance, frequency and chirp mass posterior distributions, stand-alone panel) for the CGW injected in the dastaset injCGW5 (left panel) and injCGW20 (right panel). The signal is searched with a CGW-only model with ET-only (orange distributions) and with the complete ET + PT template (purple distributions). Injected values are shown as black solid lines.

In the text

	Fig. 4. Signal recovered with the ET-only template applied to the 5 nHz CGW (dataset injCGW5, top panel) and to the 20 nHz CGW (data injCGW20, bottom panel). Median and 90% credible interval of the signal recovery are shown for pulsar PSR J1022+1001 (top) and PSR J1713+4707 (bottom). Solid and dashed lines are respectively the injected total and the ET-only signals.
In the text

Fig. 5. Injected CGW residuals (black solid line) and the reconstructed residuals with median and 90% C.I. (solid line and shaded area) are shown in the first column. The second column shows the pulsar distance prior distribution (dashed line) and recovered posterior distribution (histogram) and the injected value of the distance (solid black line). The first two rows correspond to the results obtained by analysing injCGW20 dataset and are shown for pulsar PSR J1713+0747 (first row) and pulsar PSR J1918-0642 (second row); the last row shows the results obtained with injCGW5 dataset displayed for pulsar PSR J1713+0747.

In the text

	Fig. 6. Posterior distribution of the GWB parameters obtained from dataset injCGW5 with the GWB + CGW ET-only (purple) and GWB + CGW ET + PT (orange) template.
In the text

	Fig. 7. Posterior distribution of the powerlaw GWB parameters. Purple empty posterior is obtained from dataset injGWB02, orange full posterior is obtained from dataset injCGW05. The solid line is the ‘nominal’ GWB injection (γ ∼ 13/3, log10AGWB ∼ = − 14.66).
In the text

	Fig. 8. Optimal statistics S/N for an HD-correlated (empty violins) and dipole-correlated (filled violins) process as a function of the assumed γ power law index. The two upper panels show in purple the results for datasets injGWB01 (top) and injGWB02 (bottom), while the two lower panels show in orange the results for datasets injCGW5 (top) and injCGW20 (bottom).
In the text

Fig. 9. Overlap reduction function obtained with the optimal statistic. Purple circular and orange squared points indicate the results for the injGWB02 and injCGW5 datasets respectively. The correlation coefficients for each pair of pulsars are grouped in 10 bins of 30 pulsar pairs each, the points represent the weighted average in each bin. The solid line shows the expectation value of the HD correlation. Note that the purple rightmost point lie below the orange one.

In the text

Fig. 10. Posterior distributions of the CGW parameters, inclination angle, chirp mass, luminosity distance, frequency, initial phase, and polarization angle, recovered with a CGW-only model (considering both Earth and pulsar term) applied to the three GWB only datasets. The inset sky map shows the reconstruction of the (fake) CGW sky location. The black stars in the sky map represent the positions of the 25 EPTA pulsars, while the two red dots are the Virgo and the Fornax clusters, there as a reference.

In the text

Fig. 11. Injected and recovered signals in the three injGWB datasets. In each panel, the green stars represent the contributions of each binary of the injected population to the GWB, the black line is the total power spectral density per frequency bin. The orange violin plot is the HD-correlated free spectrum recovery and the purple violin is the posterior distribution of the recovered CGW amplitude at the location of the recovered CGW frequency.

In the text

	Fig. 12. Posterior distributions obtained from GWB-only, CGW-only, CURN + CGW and GWB + CGW searches performed on the GWB-only dataset injGWB02 (left) and on the CGW-only dataset injCGW5 (right).
In the text

Fig. 13. CGW (left) and GWB (right) recovered by the joint search applied on dataset injGWB03CGW5. Left panel: Posterior distributions of the CGW source parameter using the full ET + PT template. The panel is organized as the ones in Fig. 3, with 6 source parameters shown in the triangle plot, the sky localization error shown in the sky map and the recovered signal amplitude h shown in the stand-alone inset. Injected values are marked by the black solid lines. Right panel: Posterior distributions of the GWB parameters. The GWB recovery from the joint search on the injGWB03+CGW5 dataset is shown in orange. For comparison, the recovery of a GWB only search on the injGWB03CGW5 is shown in red while the recovery of a GWB only search on the GWB-only injGWB03 dataset is shown in purple. Solid black lines represent the injected values. In both panels we mark for each parameter the median and 68% credible interval of the recovery.

In the text

Fig. 14. Same as Fig. 13 but for dataset injGWB03+CGW5. The only difference is that in the right panel, we show the joing GWB + CGW analysis either modelling the GWB as a CURN and then using the re-weighting technique (filled purple) or by directly modelling the HD correlations in the GWB search (filled orange). All other distributions and markers have the same colour style and meaning as in Fig. 13.

In the text

Fig. 15. Median and 90% credible region of the recovered CGW signal for pulsar PSR J1713+4707 from the joint analysis of dataset injGWB03+CGW5. The top panel shows the result obtained with a CURN + CGW (ET + PT) search and bottom panel shows the outcome of a GWB + CGW (ET + PT) search. Solid ad dashed line are respectively the total and the ET-only signal injected. In the first case, the CGW model neglects the pulsar term, resulting in a biased CGW parameters recovery (detail explanation in the main text).

In the text