Open Access
Issue
A&A
Volume 690, October 2024
Article Number A118
Number of page(s) 14
Section Astrophysical processes
DOI https://doi.org/10.1051/0004-6361/202348568
Published online 03 October 2024

© The Authors 2024

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

1. Introduction

The population of supermassive black hole binaries (SMBHBs) in the relatively local Universe is the most promising astrophysical source of gravitational waves (GWs) at nanohertz frequencies, probed by pulsar timing array (PTA) observations. The signal is generated by binaries in wide orbits with periods of months to years. Each binary is far from a merger and evolving slowly, so the emitted GWs are almost monochromatic. However, the incoherent superposition of GWs from many binaries creates a stochastic GW background (SGWB) signal with a characteristic broad red-noise-type spectrum. A search of the second data release (DR2) of the European Pulsar Timing Array (EPTA) for an SGWB was reported in EPTA Collaboration and InPTA Collaboration (2023a). This analysis reported increasing evidence of an SGWB, based on seeing a red noise process with a common spectral shape in all pulsars and seeing evidence that the correlation of the signal between pairs of pulsars was consistent with the forecasted Hellings-Downs (HD) correlation curve that is expected from an SGWB. The statistical significance reported in EPTA Collaboration and InPTA Collaboration (2023a) is not yet high enough to claim a detection, but the data are starting to show some evidence of an SGWB.

The EPTA DR2 includes 25 pulsars selected to optimise for detection of the HD correlations, based on the methods described in Speri et al. (2023). The analysed data were collected with six EPTA telescopes: the Effelsberg Radio Telescope in Germany, the Lovell Telescope in the UK, the Nançay Radio Telescope in France, the Westerbork Synthesis Radio Telescope in the Netherlands, the Sardinia Telescope in Italy and the Large European Array for Pulsars. For this paper, we used the DR2new subset of all of the data from EPTA Collaboration and InPTA Collaboration (2023a). It includes only the last 10.3 years of data, which were collected with new-generation wide-band backends.

The analysis of the EPTA data DR2new has reported evidence of an HD-correlated background as opposed to a common red noise process (EPTA Collaboration and InPTA Collaboration 2023a). However, this evidence may have resulted from a stochastic GW signal or from one or a few bright continuous gravitational wave (CGW) sources in the data (e.g. see Allen 2023 for details). Therefore, it is important to determine whether such CGW sources are present. The aim for this paper was to search for individual CGW sources in the EPTA data DR2new.

For this work, we used frequentist and Bayesian approaches to search for a CGW source. We adopted the model of a single binary system in a circular orbit. We analysed the data using both the Earth-term-only and a full signal (Earth + pulsar terms) model. For each pulsar, we assumed the custom-made noise model reported in EPTA Collaboration and InPTA Collaboration (2023b). We also allowed for the presence of a common red noise (CRN) component in the data. Evidence of a CRN was reported in the analysis of the reduced EPTA DR2 dataset comprising the six pulsars with the best timing accuracy (Chen et al. 2021). The six-pulsar dataset was not informative on the nature of this CRN signal, but the more recent 25-pulsar analysis reported in EPTA Collaboration and InPTA Collaboration (2023a) favours an SGWB origin for this noise. In this analysis, we consider models that include a deterministic CGW signal and one of three different noise models: individual pulsar noises only (PSRN), PSRN plus a common uncorrelated red noise (CURN) process, or PSRN plus an SGWB with Hellings-Downs (GWB) correlations between the pulsars. Simple power-law power spectral densities represent all common noises.

We conducted a Bayesian search for CGWs across a wide frequency band by splitting the dataset into sub-bands of width Δlog10fgw = 0.05. We followed up on the most significant candidate from this search with a detailed analysis. In addition, we performed an analysis on simulated datasets generated with noise properties consistent with the posterior distribution inferred from the actual data. The evaluated Bayes factors for an SGWB model versus SGWB with the addition of a CGW are close to unity. Therefore, the evidence provided by the data is equally consistent with both the inclusion or the absence of a CGW source. In other words, there is no preference for one hypothesis over the other based on the data, even though the CGW model has 58 additional parameters. This indicates that more data are needed to make a conclusive decision.

The paper is organised as follows. In Sect. 2 we describe the model used to describe the data and the frequentist and Bayesian methods we employed for our analysis. In Sect. 3 we present the results of the analysis of the EPTA DR2 dataset. In Sect. 4 we discuss results from a simulation study that we undertook to understand the results of the EPTA DR2 analysis. Finally, in Sect. 5 we summarise our results and current findings.

2. Methods

2.1. Noise model

We adopt the model for the noise in a single pulsar described in EPTA Collaboration and InPTA Collaboration (2023b), in which timing residuals are written as

δ t = M ϵ + n WN + n RN + n DM + n Sv PSRN + n CRN Common Red Noise + s CGW . $$ \begin{aligned} \delta t = \underbrace{\mathrm{M} \epsilon + n_{\mathrm{WN}} + n_{\mathrm{RN}} + n_{\mathrm{DM}} + n_{\mathrm{Sv}}}_{\rm PSRN} + \,\underbrace{n_{\mathrm{CRN}}}_{\rm Common \, Red \, Noise} + \underbrace{s}_{\rm CGW}. \end{aligned} $$(1)

The timing model error, Mϵ, is represented by a linear model based on the design matrix M and an offset from the nominal timing model parameters, ϵ. The white noise component nWN is described by two parameters for each backend, which apply multiplicative (EFAC) and additive (EQUAD) corrections to the estimated timing uncertainty. The pulsar red noise, nRN, dispersion measure variations, nDM, and scattering variations, nSv, are each represented by an incomplete Fourier basis defined at i/Tobs frequency bins (i is integer). The amplitudes are assumed to be generated by a stationary Gaussian process (van Haasteren & Vallisneri 2014), with PSD described by a power-law, characterised by spectral index, γ, and the amplitude, A, at reference frequency fref = 1/year. The noise models, including the number of frequencies in the Fourier basis, are customised for each pulsar, as described in EPTA Collaboration and InPTA Collaboration (2023b). We call the model that includes all of the aforementioned noise components the PulSaR Noise (PSRN) model.

We also allow for the presence of a common red noise (CRN), nCRN, affecting all the pulsars, that can take the form of an uncorrelated noise among pulsars (CURN) or a gravitational wave background (GWB) with a correlation described by the HD curve. We model the properties of the CRN similarly to the individual pulsar red noises, using an incomplete Fourier basis, with amplitudes described by a Gaussian process with a power-law PSD. In the Bayesian analysis below, we have used either three or nine Fourier bins for describing the CRN, and we have adopted the same priors for the pulsar noise components as presented in EPTA Collaboration and InPTA Collaboration (2023a). We refer the reader to EPTA Collaboration and InPTA Collaboration (2023b), Chalumeau et al. (2022) for a more complete description of the noise models.

The final component of the model for the timing residuals is a continuous gravitational wave (CGW) signal, s. This is described in the next section.

2.2. Continuous gravitational wave model

A supermassive black hole binary (SMBHB) system in a circular orbit produces monochromatic and quasi-non-evolving GWs (Arzoumanian et al. 2023; Falxa et al. 2023). Such signals induce pulsar timing residuals s a ( t , Ω ̂ ) $ s_a(t, \hat{\Omega}) $ of the form:

s a ( t , Ω ̂ ) = A = + , × F A ( Ω ̂ ) [ s A ( t ) s A ( t τ a ) ] , $$ \begin{aligned} s_a(t, \hat{\Omega }) = \sum _{A=+,\times } F^A(\hat{\Omega })[s_A (t) - s_A (t - \tau _a)], \end{aligned} $$(2)

where sA(t) and sA(t − τa) are referred to as the Earth term (ET) and the Pulsar term (PT), F A ( Ω ̂ ) $ F^A(\hat{\Omega}) $ are the antenna pattern functions that characterise how each GW polarisation, +, ×, affects the residuals as a function of the sky location of the source, Ω ̂ $ \hat{\Omega} $ is the direction of propagation of the GW and τa is a delay time between the source and pulsar a. The full expressions for sA(t) are:

s + ( t ) = M 5 / 3 d L ω ( t ) 1 / 3 { sin [ 2 Φ ( t ) ] ( 1 + cos 2 ι ) cos 2 ψ 2 cos [ 2 Φ ( t ) ] cos ι sin 2 ψ } , $$ \begin{aligned} s_+(t) &= \frac{\mathcal{M} ^{5/3}}{d_L \omega (t)^{1/3}} \left\{ -\sin \left[2\Phi (t)\right](1+\cos ^2 \iota ) \cos 2\psi \right. \nonumber \\& \qquad \qquad \qquad \left.- 2\cos \left[2\Phi (t)\right]\cos \iota \sin 2\psi \right\} ,\end{aligned} $$(3)

s × ( t ) = M 5 / 3 d L ω ( t ) 1 / 3 { sin [ 2 Φ ( t ) ] ( 1 + cos 2 ι ) cos 2 ψ + 2 cos [ 2 Φ ( t ) ] cos ι sin 2 ψ } , $$ \begin{aligned} s_\times (t) &= \frac{\mathcal{M} ^{5/3}}{d_L \omega (t)^{1/3}} \left\{ -\sin \left[2\Phi (t)\right](1+\cos ^2 \iota ) \cos 2\psi \right. \nonumber \\& \qquad \qquad \qquad \left. + 2\cos \left[2\Phi (t)\right]\cos \iota \sin 2\psi \right\} , \end{aligned} $$(4)

with ℳ the chirp mass, dL the luminosity distance to the source, ω(t) = πfgw(t) the time dependent frequency of the GW, ι the inclination, ψ the polarisation angle, and Φ(t) the time dependent phase of the GW. The amplitude of the GW is given by:

h = 2 M 5 / 3 d L ( π f gw ) 2 / 3 . $$ \begin{aligned} h = 2\frac{\mathcal{M} ^{5/3}}{d_L} (\pi f_{gw})^{2/3}. \end{aligned} $$(5)

For a slowly evolving binary, ω(t) is considered constant (ω(t) = ω0) over the duration of PTA observations of ∼10 years, giving for the Earth and Pulsar phases:

Φ ( t ) = Φ 0 + ω 0 t , $$ \begin{aligned} \Phi (t)&= \Phi _0 + \omega _0 t,\end{aligned} $$(6)

Φ ( t τ a ) = Φ 0 + Φ a + ω ( t τ a ) t . $$ \begin{aligned} \Phi (t-\tau _a)&= \Phi _0 + \Phi _a + \omega (t-\tau _a) t. \end{aligned} $$(7)

Nonetheless, the difference in frequency between Earth and Pulsar terms can be significant. The frequency of the pulsar term can be computed using the leading order radiation reaction evolution:

ω ( t ) = ω 0 [ 1 256 5 M 5 / 3 ω 0 8 / 3 ( t t 0 ) ] 3 / 8 . $$ \begin{aligned} \omega (t) = \omega _0 \bigg [ 1 - \frac{256}{5} \mathcal{M} ^{5/3} \omega _0 ^{8/3} (t - t_0) \bigg ]^{-3/8}. \end{aligned} $$(8)

This difference in ω(t) is determined by the time delay τa given by:

τ a = L a ( 1 + Ω ̂ · p ̂ a ) , $$ \begin{aligned} \tau _a = L_a (1 + \hat{\Omega } \cdot \hat{p}_a), \end{aligned} $$(9)

where La the distance between the Earth and pulsar a and p ̂ a $ \hat{p}_a $ is a unit vector pointing to pulsar a. If the SMBHB has significantly evolved during the time τa, the Earth term will have a higher frequency than the Pulsar term. This will usually be the case for frequencies above ∼10 nHz. For binaries at lower frequencies, binary evolution is typically negligible, and both terms will have the same frequency (within the resolution of the PTA) but different phases. The characterisation of the pulsar term can be difficult because the distance La is known with poor accuracy. Consequently, the pulsar distance La and phase Φa must be treated as free parameters that are fitted while searching for the signal (Corbin & Cornish 2010). In our analysis, we use a Gaussian prior on the distances La with the measured mean, μa, and uncertainty, σa, from Verbiest et al. (2012)1. For the pulsars not included in that paper, we use a mean of 1 kpc and an error of 20%.

2.3. Frequentist analysis

We analyse the data using a frequentist approach based on the Earth term-only ℱe-statistic (Babak & Sesana 2012; Ellis et al. 2012). The ℱe detection statistic is the log-likelihood maximised over the ‘extrinsic’ CGW parameters (h, ι, ψ, ϕ0) for a fixed set of intrinsic parameters (θ, ϕ, fgw). If the residuals are Gaussian, the null distribution is expected to be a χ2 distribution with 4 degrees of freedom. In the presence of a signal, ℱe is distributed as a non-central χ2-distribution with non-centrality parameter related to the square of the signal-to-noise-ratio (s(t)|s(t))2 (see Ellis et al. 2012 for further details). However, to calculate ℱe, we need to make assumptions about the noise properties. We take two different approaches: (i) we use the posterior distributions obtained from fitting the noise parameters to obtain a distribution of ℱe for each set of intrinsic parameters; (ii) we fix the noise parameters to their maximum likelihood estimates, as is often done for the EFAC and EQUAD parameters. The second approach is standard within frequentist analysis, but we also use (i) for the red noise components to emphasise that the inferred parameters have rather significant uncertainties. Varying the noise parameters generates a distribution of the optimal statistic for each choice of intrinsic parameters and thus brings an element of Bayesian approaches into this frequentist analysis.

We want to evaluate the significance of the highest ℱe measured on the observed data by computing the p-value, which is a statement about how improbable it would be to draw the observed data if no signal was present. Calculating the true p-value requires the true distribution of ℱe in the absence of signal (the null distribution), which we do not have. Building a proper null distribution for the ℱe statistics preserving the properties of the noise (of observed data) requires an extensive research programme and is outside the scope of this paper. Instead, we propose two simplified ways of estimating the significance: (i) using the theoretical null distribution of 2ℱe, assuming that our noise model is correct and complete. The theoretical distribution behaves as a χ2 with 4 degrees of freedom, (ii) shuffling the pulsars on the sky but keeping the “source” position fixed to its maximum likelihood value. The second approach requires further explanations. The signal-to-noise ratio for the candidate event is dominated by several pulsars, which implies that a simple shuffle of the pulsars will not destroy a monochromatic GW but will change its sky position and extrinsic parameters (initial phase, inclination, etc). By shuffling the pulsars on the sky, we ask, “how likely to obtain the observed value of ℱe by random position of pulsars”. This approach has the advantage of making no distributional assumptions about the pulsars’ noise properties. Further, it requires no minimal match or any orthogonality condition between shuffled pulsar positions (see discussion on this topic in Taylor et al. 2017; Marco et al. 2023) The sky shuffled distribution was obtained as follows:

  • (i) We produce 3000 shuffles by drawing random pulsar positions uniformly in the sky.

  • (ii) For each of the 3000 draws, we evaluate ℱe for 1000 noise parameters drawn from the posterior distributions obtained in EPTA Collaboration and InPTA Collaboration (2023b) and take the mean value.

  • (iii) We produce a histogram of the 3000 ℱe statistic values, representing the background distribution.

  • (iv) We repeat (i–iii) 20 times, obtaining a slightly different histogram each time due to differences in the 3000 shuffles and the noise realisations. This allows us to estimate the uncertainty in the background distribution and, therefore, in the computed p-value.

We will apply this procedure to the CGW candidate in Sect. 3.

2.4. Bayesian analysis

We also conducted a Bayesian analysis to obtain posterior probability distributions for the noise and signal parameters in the model described in Sect. 2. We make use of MCMC samplers PTMCMC (Ellis & van Haasteren 2017), QuickCW (Bécsy et al. 2022a) and Eryn (Karnesis et al. 2023) to explore the parameter space. The Bayesian analysis allows us to perform parameter inference and model selection. The latter is quantified by evaluating the Bayes factor: the ratio of the marginal posterior distributions (or evidence) for two different models. The marginal posterior is a quantity challenging to compute and can be estimated numerically using parallel tempering or nested sampling (Skilling 2004).

In this paper, we use ENTERPRISE (Ellis et al. 2020; Taylor et al. 2021) to evaluate the posterior probability for a given model. We compute the Bayes factors using the product-space method (Hee et al. 2016) implemented in ENTERPRISE and through Reversible Jump Markov chain Monte Carlo (RJMCMC), as implemented in Eryn. In both approaches, at each step of the Markov chain, either the parameters within the current model can be updated, or a switch to a different model can be proposed. The acceptance rule for the model switch is defined to ensure that detailed balance is maintained, thus ensuring that the stationary distribution of the Markov chain is the desired posterior distribution over models. The sampler will spend more time exploring the model with the highest marginal posterior probability. The Bayes factor ℬA/B between models A and B can then be calculated as the ratio between the final number of chain samples corresponding to each model.

In the product-space method, the chain samples are in a hypermodel space, which is a union of all the parameters of all the models being considered. An additional parameter determines which model is active within each sample, while inactive parameters undergo a random walk during the within-model steps. The effect is that the product space method retains some memory of where it had been exploring the other models, which can increase the probability that a proposed switch back to the other model is accepted. In RJMCMC, the chain typically only samples in the parameters of the currently active model and does not retain any memory. This can lead to lower model-switch acceptance rates but guarantees a more complete exploration of the parameter spaces of the different models. It is important to note that RJMCMC methods applied to PTA data analysis have been previously explored in Bécsy & Cornish (2020).

3. Results of data analysis

3.1. Frequentist analysis

Within the frequentist approach, we want to maximise the detection statistic (ℱe in our case) over all intrinsic parameters of the model. We perform the search using the noise models described in Sect. 2 for 100 logarithmically spaced frequencies from 1 to 100 nHz dividing the sky into 3072 different pixels using healpix (Zonca et al. 2019; Górski et al. 2005).

To account for the fact that the noise model has broad posteriors, we use the posterior samples of the noise parameters obtained in EPTA Collaboration and InPTA Collaboration (2023b) to calculate ℱe for fixed CGW parameters and average ℱe over 1000 randomly drawn samples of the PSRN model. The maximum of ℱe is found at 4.64 nHz, consistent with the results of the Bayesian analyses described in Sect. 3.2.

The sky distribution of ℱe at this frequency is given in Figure 1. The region of high statistic value (bright yellow) is relatively sparse and inconclusive with regard to the localisation of the CGW candidate. The maximum ℱe is depicted by a black star and corresponds to a region of the sky where we lack pulsars and, hence where the array is expected to be less sensitive.

thumbnail Fig. 1.

e-statistic of the candidate source at fgw = 4.64 nHz averaged over the noise uncertainties for the custom PSRN model. The black star shows the position of highest ℱe, whereas the red stars show the positions of the pulsars. The Fornax and Virgo clusters are shown as black dots.

The analysis was repeated by including the CURN component in the noise model, and results are presented in Figure 2. We show two distributions of ℱe. These are both evaluated at the optimal sky position and at the GW frequency 4.64 nHz and are obtained by varying the noise parameters (random draw) with (orange histogram) and without (blue histogram) the CURN component. Inclusion of the CURN slightly reduces the significance of the CGW candidate.

thumbnail Fig. 2.

Distribution of the ℱe-statistic over the noise uncertainties without CURN (blue) and with CURN (orange) at 4.64 nHz. The null distributions of the ℱe are obtained from the analysis of the EPTA DR2new data with shuffled sky positions (grey shaded region) and from the theoretical formula of a χ42-distribution (black solid line).

To evaluate the p-value, we compute the distribution of sky-shuffled pulsars according to the steps outlined in Sect. 2.3 at the CGW candidate parameter values (maximising the noise averaged ℱe). These results are indicated by the grey shaded distribution in Figure 2. The theoretical χ2 null distribution is shown as a black curve. The sky-shuffled distribution of 2ℱe is close to the theoretical χ42 distribution but not identical. This could be due to (i) non-gaussian noise present in the array; (ii) we are asking different questions in two approaches.

We compute p-values using the obtained distributions and the measured mean values of the orange and blue ℱe distributions. The results are summarised in Table 1, the top row is for the theoretical distribution and the second row is for the sky-shuffled distribution (with uncertainty). The obtained p-value for ℱe corresponds to about 3.5σ while ℱe, CURN corresponds to about 2.5σ.

Table 1.

Statistical significance of the candidate source at 4.64 nHz.

In order to apply the full detection pipeline, we would need to set the detection threshold based on the chosen false alarm probability. This, in turn, would require estimating the number of independent cells on the sky-frequency parameter space where the ℱe was evaluated (Babak et al. 2016; Ellis et al. 2012). Here, we decided to use the ℱe only to identify potential candidates and produce sky maps of significance without a precise statement of detection within the frequentist approach.

Estimating the p-value using ℱe. In the procedure that we describe, we select the pixel with highest ℱe, which is the same as maximising over sky location. Hence, the correct way to produce the background distribution would be to also select the highest ℱe pixel for each scramble. But when doing so, we always find a high ℱe value. This is because the S/N of the observed signal is dominated by a few pulsars (∼3 pulsars). One needs at least 3 pulsars for computing ℱe, based on a simple counting of the measurements and unknowns. Any additional pulsar in the array will improve S/N and the source’s sky localisation (ruling out some parts of the sky). For this reason, shuffling the pulsar position around does not help; we always have a solution based on the 3 best pulsars, and the overall value will depend on the relative contribution of other pulsars in the array. This is not the case when all pulsars are equally good timers as demonstrated in Babak & Sesana (2012). In Figure 3 we demonstrated, that sky scramble of the pulsar position still delivers a significant maximum ℱe. Note that it does not make ℱe useless, it still indicates a presence of a signal, it rather indicates that the sky scramble is not an appropriate mechanism to obtain the “null” dataset.

thumbnail Fig. 3.

Pulsar-shuffled version of the ℱe-statistic of the candidate source at fgw = 4.64 nHz averaged over the noise uncertainties for the custom PSRN model. The black star shows the position of highest ℱe, whereas the red stars show the positions of the pulsars. We see that the highest ℱe pixel changed position but is still significant. The Fornax and Virgo clusters are shown as black dots.

The reported p-value answers the question: how likely to get an observed ℱe at given point in the parameter space assuming the Gaussian noise and a random sky position (shuffling) of pulsars in the array. That said, the produced distribution is indeed not the real null distribution. In the text, we used the word sky-shuffled instead of scramble to avoid confusion with previous works. The p-values appearing in Table 1 report the highest possible significance that is achievable using the real positions of pulsars. What we are essentially doing is searching for mis-modelled CGW-like features, but this procedure does not settle whether this feature is only due to noise.

3.2. Bayesian analysis

We perform a Bayesian CGW search by splitting the frequency range 10−9 − 10−6.5 Hz into 50 logarithmically spaced subsegments and assuming Earth term only. We have computed the Bayes factor in each sub-band using PTMCMC and the product-space method. The results are presented in Figure 4 for the noise model described in Sect. 2. We find a Bayes factor above 100 around 4–5 nHz, and perform a model comparison analysis in a restricted frequency range fgw ∈ [10−9, 10−7.85] with the priors defined in Table 2. We use Bayes factors as the decision-maker between models. We use Eryn (Karnesis et al. 2023) as our fiducial sampler and we crosscheck the results using PTMCMC sampler (Ellis & van Haasteren 2017) and QuickCW (Bécsy et al. 2022a). The computation of the Bayes factors is performed using RJMCMC and confirmed with the product-space method. Our findings are summarised in Table 3, and in the following, we comment on them. For all models with CGW described below and quoted in the table, we have assumed a circular binary described in Sect. 2.2 with Earth and pulsar term, using the priors given in Table 2 unless otherwise stated. The validity of these assumptions is discussed at the end of the section.

thumbnail Fig. 4.

Bayes factor for the model comparison PSRN+CURN+CGW (Earth term) versus PSRN+CURN for 50 logarithmically spaced frequency sub-bands in the region fgw ∈ [1.5, 320] nHz.

Table 2.

Parameters of the continuous gravitational wave model with their respective priors and ranges.

Table 3.

Bayes factors obtained for different models.

The simplest considered data model includes only the custom pulsar noise (PSRN), therefore, no CRN is included. The PSRN model is used as a null hypothesis, and the alternative is given by the pulsar noise plus CGW (PSRN+CGW) considering Earth and pulsar terms. The Bayes factor for the model comparison PSRN+CGW vs PSRN is 4000.

Next, we include a CRN to the custom noise PSRN: a CURN or a GWB correlated according to the HD pattern. These become the new null hypotheses (CURN+PSRN) and (GWB+PSRN). We also consider two descriptions for the CRN: one using the three lowest Fourier harmonics (3 bins) and one using the nine lowest Fourier harmonics (9 bins) as in EPTA Collaboration and InPTA Collaboration (2023a). For reference, the maximum resolvable frequency for three and nine bins are log10(3/(10.3 yr))  ≈   − 8.03 and log10(9/(10.3 yr))  ≈   − 7.56, respectively. In practice, we should let the data choose the number of frequency bins similarly to what was done in EPTA Collaboration and InPTA Collaboration (2023b). However, for practical reasons, we present the results only for 3 and 9 bins. The maximum resolvable frequencies smaller than the one over-a-year frequency dictate the choice of nine frequency bins. Since the power-law model adopted to describe the background has two parameters, the minimum number of frequencies we can fit is three. Therefore, we also present the model selection results for a choice of three bins. Note that the free spectrum analysis presented in EPTA Collaboration and InPTA Collaboration (2023a) shows that the first, fourth, and ninth bins are constrained with a tail extending to low amplitudes. Only the second bin is constrained with zero support at low powers. Since the CGW candidate is located close to the second Fourier bin, showing the results for three bins can help to single out the red noise components of the spectrum that might be potentially affected by the other high-frequency noises. The Bayes factors of PSRN+CURN+CGW vs PSRN+CURN are 4 and 12, for 9 and 3 bins, respectively. The choice of the number of bins affects the spectral properties of the CRN and, consequently, also the Bayes factors. In fact, the slope of the CURN model becomes steeper when using 3 bins, allowing a possible CGW to emerge from the noise. When including the HD correlations, the Bayes factors of PSRN+GWB+CGW vs PSRN+GWB drop to 0.7 and 1, for 9 and 3 bins, respectively.

As already was pointed out in EPTA Collaboration and InPTA Collaboration (2023a), the HD component of the noise absorbs most of CGW signal and this can be clearly seen in the drop of the Bayes factor and in the posterior distributions shown in Figure 5. When the CRN is a CURN, the CGW model absorbs the power of the background around log10fgw ∈ [ − 8.5, −8.2] and yields an amplitude log10A posterior distribution with tails extending up to the lowest end of the prior range (see correlation in Figure 5 for parameters log10fgw, log10A).

thumbnail Fig. 5.

Posterior distributions of the CGW search in the second data release of the EPTA DR2new. The posteriors are obtained using a CGW model with Earth and pulsar term, the custom PSRN model and a CRN with either CURN or HD (GWB) correlations represented by 9 frequency bins. We show the posterior distribution for the gravitational wave frequency and amplitude fgw, h of the CGW, and the common noise spectral index and amplitude γ and A. The contours indicate the 1,2,3-σ Gaussian contours.

For the model CURN+PSRN+CGW with 9 bins, the log-frequency fgw is measured to be 4.61−2.98+1.11 nHz and the log-amplitude log10h is measured to be −14.0−2.6+0.5 (median and symmetric 90% credible interval). The chirp-mass posterior is uninformative and the sky localisation posterior is shown in Figure 6 where we also show the Virgo and Fornax clusters which are 16.5 Mpc and 19.3 Mpc from the Earth Jordan et al. (2007), respectively. If we use the median values of the amplitude and frequency to estimate the luminosity distance, we obtain dL ≈ 16.6 (ℳ/109 M)5/3 Mpc, using Eq. (5).

thumbnail Fig. 6.

Posterior distribution of the sky localisation obtained by searching for a CGW in the second data release of the EPTA DR2new. The posteriors are obtained using a CGW model with Earth and pulsar term, the inclusion of the custom pulsar noise (PSRN) and a common uncorrelated red noise (CURN) represented by 9 frequency bins. For reference, we show the position of the analysed pulsars and the Virgo and Fornax clusters.

We compute the sky marginalised 95% upper limit on strain amplitude, h95(fgw), across the studied frequency range for the model PSRN+CURN+CGW with 9 bins. For this analysis, we used a uniform prior on h in the range [10−18, 10−11] instead of the uniform prior on log10h used for the search (see Arzoumanian et al. 2023; Falxa et al. 2023). The strain upper limit was converted into a horizon distance, DH, (i.e. the distance up to which SMBHB systems should produce detectable CGW signals) using Eq. (5):

D H = 2 M 5 / 3 h 95 ( π f gw ) 2 / 3 . $$ \begin{aligned} D_H = 2 \frac{\mathcal{M} ^{5/3}}{h_{95}}(\pi f_{gw})^{2/3}. \end{aligned} $$(10)

We plot DH as a function of fgw in Figure 7 for three values of chirp mass ℳ = [108M, 109M, 1010M]. The highest DH is recovered around 20 nHz. The closest galaxy-cluster candidates (Fornax and Virgo) that could host a SMBHB lie at distances larger than 10 Mpc, meaning that we need binary systems with chirp masses larger than 109M in order for them to be detectable. One of the known galaxies hosting a supermassive black hole in the Fornax cluster is the radio galaxy Fornax A (NGC 1316). The radio observations suggest the presence of a powerful AGN, and thus a central supermassive black hole, whose estimated mass is M = 1.5 × 108M (Nowak et al. 2008). The mass and distance of such a black hole put it outside of the horizon distance of our PTA experiment.

thumbnail Fig. 7.

Horizon luminosity distance, DH, obtained from the sky averaged 95% upper limit on strain amplitude h using the PSRN+CURN+CGW (Earth term + pulsar term) model. The horizon distance is calculated with Eq. (10) for three chirp masses: 108M, 109M and 1010M. For reference, the Virgo and Fornax clusters are at 16.5 Mpc and 19.3 Mpc from the Earth, respectively.

For the model GWB+PSRN+CGW with 9 bins, we cannot constrain the CGW parameters, so we set a 95% upper limit on log10h95% = −13.75. We constrain the spectral properties of the GWB background to be 2.66−1.02+1.43 and −13.95−0.62+0.25, respectively for (γ and log10A) (median and symmetric 90% credible interval). In Figure 8 we show the posteriors for the same model GWB+PSRN+CGW, but this time with 9 and 3 bins. The background’s spectral properties (γ and A) are affected by the number of bins. The median log-amplitude decreases from –13.95 to –14.34 and the median slope from 2.66 to 3.832, following the typical γ − log10A correlation. The steeper slope allows the CGW to emerge from the noise and its posteriors show two clear peaks, one at 4.64 nHz and one at 12.6 nHz. Since the GWB model with 3 bins extends in frequency up to 9.33 nHz, the appearance of the second peak at 12.6 nHz might be due to unmodelled noise.

thumbnail Fig. 8.

Posterior distributions of the CGW search in the second data release of the EPTA DR2new. The posteriors are obtained using a CGW model with Earth and pulsar term, the custom PSRN model and an HD-correlated background (GWB) represented by 9 and 3 frequency bins. We show the posterior distribution for the gravitational wave frequency and amplitude fgw, h of the CGW, and the common noise spectral index and amplitude γ and A. The contours indicate the 1,2,3-σ Gaussian contours.

The inclusion of the HD correlation absorbs the presence of the CGW signal. We still need to check whether the peak frequency at 4.6 nHz for a CURN background is unaffected by the sampler choice, the pulsar term or the inclusion of eccentricity. To do this, we consider the model CURN+PSRN+CGW. The narrow green posterior contours in Figure 9 correspond to using Eryn to sample the model with the Earth-term only and the broad blue posterior is inferred with the model including the pulsar term. We have overplotted similar results (including the pulsar term) obtained with the QuickCW sampler in orange. The results are in broad agreement independently of the sampler choice and the exclusion of the pulsar term. To check if the circular CGW model is appropriate, we carried out a separate analysis including orbital eccentricity in the model following Taylor et al. (2016) but using only the Earth term in the analysis. This analysis inferred a low eccentricity (e < 0.2) for the CGW candidate, indicating that the analysis performed with the circular binary model is adequate.

thumbnail Fig. 9.

Inference of the amplitude h and frequency fgw of CGW using the PSRN+CURN+CGW model. The results obtained with Eryn are shown as green and blue histograms, for a Earth term only, and full CGW model, respectively. The blue contours are to be compared with the orange posterior obtained with QuickCW. The shown contours are the 90% normal credible regions.

3.3. Optimal statistics

The Bayesian analysis presented in the previous subsection indicates that the results about the nature of the observed signal are inconclusive. This subsection starts a long investigation process attempting to answer the question: “what is it that we see?” and provides a transition to the next section which describes analysis of the simulated data.

We compute the signal-to-noise ratio (S/N) of the CRN using the optimal statistic approach (Chamberlin et al. 2015; Vigeland et al. 2018) implemented within enterprise_extensions software package (Taylor et al. 2021). We estimate the S/N assuming quadrupolar correlation (HD). Following the procedure outlined in (EPTA Collaboration and InPTA Collaboration 2023a; Vigeland et al. 2018), we vary the pulsar noise parameters (using the custom noise models from EPTA Collaboration and InPTA Collaboration 2023b) and get as a result a distribution of S/N. The pulsar noise parameters are drawn from the posterior samples of a Bayesian analysis, which includes, in addition to single pulsar noises, a CURN, using a 9 Fourier bin basis. The solid orange line in Figure 10 reproduces the findings reported in EPTA Collaboration and InPTA Collaboration (2023a).

thumbnail Fig. 10.

Distributions of the signal-to-noise-ratio (S/N) for the HD correlations for a common red process with γ = 13/3, with (dashed) or without (solid) adding CGW to the data model. The orange lines correspond to the results obtained on DR2new and blue are on the simulated data with only CGW (no CRN).

Next we include a CGW in the model: we use all the posterior samples obtained in the Bayesian analysis which preserve the correlation between the noise parameters and the CGW (instead of only the noise parameters), and re-evaluate the optimal statistic. As a result, for each posterior sample, the S/N is computed, given the sample pulsar noise parameters, after subtracting the deterministic CGW residuals from the pulsar residuals. The CGW residuals are computed using the sample CGW parameter values. The resulting S/N distribution of HD correlations is given in Figure 10 as dashed orange line. This result implies that the data (minus CGW) does not show any sign of quadrupolar (GWB) correlation, in other words, a CGW alone can explain the HD feature observed in the data.

We corroborate our results obtained on the EPTA DR2new by repeating the same analysis on a simulated dataset. We produce a fake PTA based on the real (EPTA DR2new) pulsars and the noise estimation in which we inject only one CGW and no GWB (see Sect. 4 for a detailed description of the simulation). The blue solid line in Figure 10 indeed resembles the result obtained on the DR2new (orange line). As we will discuss in detail in the next section, a single CGW signal could be interpreted as a GWB (see Allen 2023). As expected, the subtraction of the CGW from the timing residuals removes the quadrupolar correlation (blue dashed line) and reproduces the previously obtained results on the EPTA data.

4. Simulation

We perform a simulation campaign to try to reproduce the features observed in the analysis of DR2new. We generate a fake array with the same time of arrivals (TOA)s and pulsar positions as in the real dataset. We inject noises using the maximum a posteriori of the noise parameter posterior obtained in EPTA Collaboration and InPTA Collaboration (2023b). We use a Gaussian process to simulate the noise components and consider different realisations in order to reproduce the observed results (see Appendix A). Using the simulated array as a basis, we propose two cases to study:

  • PSRN+CGW: A simulated analogue of DR2new with only one circular CGW injected at 4.8 nHz with sky location at (3h38, –35°27) as if it was in the Fornax cluster with a chirp mass of 109.2M and amplitude h = 10−13.6, without any CRN.

  • PSRN+GWB: A simulated analogue of DR2new with a gaussian and isotropic GWB as CRN with a powerlaw spectrum corresponding to A = 10−14.5 and spectral index γ = 13/3 (without any CGW).

Each simulation is analysed with the custom PSRN model and either with CGW (using the Earth term only) or GWB with a powerlaw spectrum.

We have considered the PSRN+GWB simulated data and analysed it with a single CGW source (no GWB). In Figure 11 we show that we can recover the CGW even if we have injected an isotropic GWB. We have repeated the analysis on 10 simulated datasets, in all cases the recovered “CGW” was centred at the lowest Fourier bin (1/Tobs ∼ 3 nHz) and located in the close vicinity of pulsars, in many cases around J1713+0747. Recovering the lowest frequency bin is compatible with the injected powerlaw spectrum for which it is the loudest bin. We present the posterior frequency of one case as a orange histogram in Figure 11. We have also analysed this data using GWB model and the inferred posterior is given in orange in Figure 12.

thumbnail Fig. 11.

Posterior distribution of the gravitational wave frequency fgw of a CGW fitted to two simulated PTAs: one with a single injected CGW (PSRN+CGW blue histogram), and one with an injected GWB (PSRN+GWB orange histogram). The posterior distribution is obtained with a MCMC analysis for a PSRN+CGW (Earth term) model. The dashed lines are the medians of the distributions.

thumbnail Fig. 12.

Posterior distribution of amplitude log10A and spectral index γ of the GWB for two simulated PTAs: one realistic PTA with only one CGW injected (PSRN+CGW) and one realistic PTA with a GWB injected (PSRN+GWB). The posterior distribution is obtained with a MCMC analysis for a PSRN+GWB model. The dashed lines are the medians of the distributions.

Next, we consider the PSRN+CGW simulated data and analyse it with the model of isotropic GWB with a power-law spectrum. As expected, this model gives a constrained posterior. Indeed, the excess of power at low frequency gives support to a power-law spectrum, though not the best description because we try to fit a single frequency process with a continuous spectrum. The HD correlations could be reproduced as shown in Figure 10. A single source can generate the HD pattern by averaging over the pulsar pairs (see Allen 2023; Cornish & Sesana 2013 for details). In the present case, the S/N is dominated by a few pulsars so the averaging is strongly biased and we would not necessarily expect a strong HD component. The analysis of this PSRN+CGW dataset with a CGW model is shown in Figure 11 as a blue histogram. Analysis of the same data with GWB is given in Figure 12 as a blue posterior. One can see that it has a higher amplitude and is shallower, similar to what was presented in EPTA Collaboration and InPTA Collaboration (2023a).

The analysis of the simulated data demonstrates how a CGW model can mimic a GWB and vice-versa. The point source also produces HD correlations, as previously shown in Cornish & Sesana (2013), Bécsy et al. (2022b), Allen (2023). Moreover, the anisotropic configuration of the current PTA (pulsars not uniformly distributed in the sky and having very different noise properties) produces an uneven response across the sky and the studied frequency range (see Sect. 3.1). The interactions between our configuration of pulsars and the signal models is a challenging aspect of PTAs that has been investigated in previous studies (Chalumeau et al. 2022; Speri et al. 2023). The inclusion of pulsars from the southern hemisphere (PPTA) would provide a better coverage of the sky allowing a better localisation of single sources, which is crucial to differentiate between a GWB or a CGW.

5. Summary

This paper presents an analysis of the EPTA DR2new dataset searching for continuous GW signals from super-massive black hole binaries in quasi-circular orbits. We performed a frequentist (based on the ℱe-statistic) and Bayesian (using the Bayes factor) analysis of the data, and, in both cases, found a significant CGW candidate at 4–5.6 nHz. The frequentist analysis gives a p-value of (5 × 10−4–6 × 10−3), equivalent to a 2.5–3σ significance level, depending on the evaluation procedure and whether or not a CURN is included in the noise model. Within the Bayesian analysis of the CGW candidate, we computed the Bayes factor to assess whether the data prefer the inclusion of a CGW source. We find weak evidence (Bayes factors ∼4 − 12) for the inclusion of a CGW on top of a CURN process. If the common red noise is assumed to have the HD correlation typical of a GWB, the Bayes factors for the inclusion of a CGW drop to ∼0.7 − 1. The data are equally well described by a model including both a GWB and CGW and a model including GWB only. Even though the CGW model depends on 58 parameters and therefore comes with a large dimensionality penalty, the Bayes factor is close to unity. More data will allow tests of these hypotheses to be more conclusive.

In an attempt to understand if the observed signal is due to a GWB or a CGW, we performed a simulation campaign. We simulated data based on the noise parameters inferred in EPTA Collaboration and InPTA Collaboration (2023b) and injected a GW signal. Our main finding is that simulated data with only an isotropic GWB injected can be fitted with a CGW model, and vice versa; a GWB model can explain simulated data containing only a single injected CGW. Therefore, we cannot conclusively distinguish between the presence of a single continuous gravitational wave or a gravitational wave background. In EPTA Collaboration and InPTA Collaboration (2024), considering models that produce a GW signal consistent with the one present in DR2new, the probability of detecting a single source with an S/N larger than 3 is estimated to be 50%. We hope that an analysis of the combined IPTA data (Data Release 3) will help to confirm the presence or not of a CGW signal and shed light on its nature.


1

https://www.atnf.csiro.au/research/pulsar/psrcat/

2

(s|s) denotes the noise weighted inner product, (s|s) = sC−1s, with C the covariance matrix of the noise model.

Acknowledgments

The European Pulsar Timing Array (EPTA) is a collaboration between European and partner institutes, namely ASTRON (NL), INAF/Osservatorio di Cagliari (IT), Max-Planck-Institut für Radioastronomie (GER), Nançay/Paris Observatory (FRA), the University of Manchester (UK), the University of Birmingham (UK), the University of East Anglia (UK), the University of Bielefeld (GER), the University of Paris (FRA), the University of Milan-Bicocca (IT), the Foundation for Research and Technology, Hellas (GR), and Peking University (CHN), with the aim to provide high-precision pulsar timing to work towards the direct detection of low-frequency gravitational waves. An Advanced Grant of the European Research Council allowed to implement the Large European Array for Pulsars (LEAP) under Grant Agreement Number 227947 (PI M. Kramer). The EPTA is part of the International Pulsar Timing Array (IPTA); we thank our IPTA colleagues for their support and help with this paper and the external Detection Committee members for their work on the Detection Checklist. Part of this work is based on observations with the 100-m telescope of the Max-Planck-Institut für Radioastronomie (MPIfR) at Effelsberg in Germany. Pulsar research at the Jodrell Bank Centre for Astrophysics and the observations using the Lovell Telescope are supported by a Consolidated Grant (ST/T000414/1) from the UK’s Science and Technology Facilities Council (STFC). ICN is also supported by the STFC doctoral training grant ST/T506291/1. The Nançay radio Observatory is operated by the Paris Observatory, associated with the French Centre National de la Recherche Scientifique (CNRS), and partially supported by the Region Centre in France. We acknowledge financial support from “Programme National de Cosmologie and Galaxies” (PNCG), and “Programme National Hautes Energies” (PNHE) funded by CNRS/INSU-IN2P3-INP, CEA and CNES, France. We acknowledge financial support from Agence Nationale de la Recherche (ANR-18-CE31-0015), France. The Westerbork Synthesis Radio Telescope is operated by the Netherlands Institute for Radio Astronomy (ASTRON) with support from the Netherlands Foundation for Scientific Research (NWO). The Sardinia Radio Telescope (SRT) is funded by the Department of University and Research (MIUR), the Italian Space Agency (ASI), and the Autonomous Region of Sardinia (RAS) and is operated as a National Facility by the National Institute for Astrophysics (INAF). The work is supported by the National SKA programme of China (2020SKA0120100), Max-Planck Partner Group, NSFC 11690024, CAS Cultivation Project for FAST Scientific. This work is also supported as part of the “LEGACY” MPG-CAS collaboration on low-frequency gravitational wave astronomy. JA acknowledges support from the European Commission (Grant Agreement number: 101094354). JA and SCha were partially supported by the Stavros Niarchos Foundation (SNF) and the Hellenic Foundation for Research and Innovation (H.F.R.I.) under the 2nd Call of the “Science and Society – Action Always strive for excellence – Theodoros Papazoglou” (Project Number: 01431). AC acknowledges support from the Paris Île-de-France Region. AC, AF, ASe, ASa, EB, DI, GMS, MBo acknowledge financial support provided under the European Union’s H2020 ERC Consolidator Grant “Binary Massive Black Hole Astrophysics” (B Massive, Grant Agreement: 818691). GD, KLi, RK and MK acknowledge support from European Research Council (ERC) Synergy Grant “BlackHoleCam”, Grant Agreement Number 610058. SB, HQ acknowledge funding from the French National Research Agency (grant ANR-21-CE31-0026, project MBH_waves). This work is supported by the ERC Advanced Grant “LEAP”, Grant Agreement Number 227947 (PI M. Kramer). KLi acknowledges support from the Shanghai Astronomical Observatory, Chinese Academy of Sciences, 80 Nandan Road, Shanghai 200030, China. AV and PRB are supported by the UK’s Science and Technology Facilities Council (STFC; grant ST/W000946/1). AV also acknowledges the support of the Royal Society and Wolfson Foundation. JPWV acknowledges support by the Deutsche Forschungsgemeinschaft (DFG) through thew Heisenberg programme (Project No. 433075039) and by the NSF through AccelNet award #2114721. NKP is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – Projektnummer PO 2758/1–1, through the Walter–Benjamin programme. ASa thanks the Alexander von Humboldt foundation in Germany for a Humboldt fellowship for postdoctoral researchers. APo, DP and MBu acknowledge support from the INAF Large Grant 2022 “GCjewels” (P.I. Andrea Possenti) approved with the Presidential Decree 30/2022 (Italy). RNC acknowledges financial support from the Special Account for Research Funds of the Hellenic Open University (ELKE-HOU) under the research programme “GRAVPUL” (K.E.-80383/grant agreement 319/10-10-2022: PI N. A. B. Gizani). EvdW, CGB and GHJ acknowledge support from the Dutch National Science Agenda, NWA Startimpuls - 400.17.608. BG is supported by the Italian Ministry of Education, University and Research within the PRIN 2017 Research Program Framework, n. 2017SYRTCN. LS acknowledges the use of the HPC system Cobra at the Max Planck Computing and Data Facility. The Indian Pulsar Timing Array (InPTA) is an Indo-Japanese collaboration that routinely employs TIFR’s upgraded Giant Metrewave Radio Telescope for monitoring a set of IPTA pulsars. BCJ, YG, YM, SD, AG and PR acknowledge the support of the Department of Atomic Energy, Government of India, under Project Identification # RTI 4002. BCJ, YG and YM acknowledge support of the Department of Atomic Energy, Government of India, under project No. 12-R&D-TFR-5.02-0700 while SD, AG and PR acknowledge support of the Department of Atomic Energy, Government of India, under project no. 12-R&D-TFR-5.02-0200. KT is partially supported by JSPS KAKENHI Grant Numbers 20H00180, 21H01130, and 21H04467, Bilateral Joint Research Projects of JSPS, and the ISM Cooperative Research Program (2021-ISMCRP-2017). AS is supported by the NANOGrav NSF Physics Frontiers Center (awards #1430284 and 2020265). AKP is supported by CSIR fellowship Grant number 09/0079(15784)/2022-EMR-I. SH is supported by JSPS KAKENHI Grant Number 20J20509. KN is supported by the Birla Institute of Technology & Science Institute fellowship. AmS is supported by CSIR fellowship Grant number 09/1001(12656)/2021-EMR-I and T-641 (DST-ICPS). TK is partially supported by the JSPS Overseas Challenge Program for Young Researchers. We acknowledge the National Supercomputing Mission (NSM) for providing computing resources of ‘PARAM Ganga’ at the Indian Institute of Technology Roorkee as well as ‘PARAM Seva’ at IIT Hyderabad, which is implemented by C-DAC and supported by the Ministry of Electronics and Information Technology (MeitY) and Department of Science and Technology (DST), Government of India. DD acknowledges the support from the Department of Atomic Energy, Government of India through Apex Project - Advance Research and Education in Mathematical Sciences at IMSc. The work presented here is a culmination of many years of data analysis as well as software and instrument development. In particular, we thank Drs. N. D’Amico, P. C. C. Freire, R. van Haasteren, C. Jordan, K. Lazaridis, P. Lazarus, L. Lentati, O. Löhmer and R. Smits for their past contributions. We also thank Dr. N. Wex for supporting the calculations of the galactic acceleration as well as the related discussions. The EPTA is also grateful to staff at its observatories and telescopes who have made the continued observations possible. Author contributions. The EPTA is a multi-decade effort and all authors have contributed through conceptualisation, funding acquisition, data-curation, methodology, software and hardware developments as well as (aspects of) the continued running of the observational campaigns, which includes writing and proofreading observing proposals, evaluating observations and observing systems, mentoring students, developing science cases. All authors also helped in (aspects of) verification of the data, analysis and results as well as in finalising the paper draft. Specific contributions from individual EPTA members are listed in the CRediT (https://credit.niso.org/) format below. We especially thank K. Grunthal for having found a bug in the implementation of the F-statistics. InPTA members contributed in uGMRT observations and data reduction to create the InPTA data set which is employed while assembling the DR2full+ and DR2new+ data sets. Additionally, InPTA members contributed to GWB search efforts with DR2full+ and DR2new+ data sets and their interpretations. Further, they provided quantitative comparisons of GWB posteriors that arise from these data sets and multiple pipelines. For this work specifically, SB, MF, LS equally share the correspondence of the paper. CRediT statement: Conceptualization Conceptualization: AShe, av, GMS, HQL, LS, LSe, MF, MK. Methodology Methodology: HQL, IF, JG, JWMo, LS, LSe, MF, MK. Software Software: AC, GMS, HQL, JWMo, LS, LSe, MF, MJK, SC. Validation Validation: DI, HQL, IF, JWMo, LS, LSe, MF. Formal Analysis Formal Analysis: GMS, HQL, IF, JG, LS, LSe, MF. Investigation Investigation: DP, HQL, I, LS, LSe, MBM, MF, SAS. Resources Resources: AC, AP, GJu, GT, I, JPWVi, KL, LS, LSe, MF, MJK, MK. Data Curation Data Curation: AC, DP, GMS, HH, I, KL, LS, LSe, MBa, MBM, MF, MJK, SAS, SC. Writing - Original Draft Writing - Original Draft: GMS, HQL, LS, LSe, MF. Writing - Review & Editing Writing - Review & Editing: AC, AFo, av, GMS, HQL, JG, KL, LS, LSe, MF, SC. Visualization Visualization: HQL, KL, LS, LSe, MF. Supervision Supervision: AP, AShe, av, GT, JG, JPWVi, LS, LSe, MF. Project Administration Project Administration: av, GT, JPWVi, LS, LSe, MF, MK. Funding Acquisition Funding Acquisition: AP, AShe, av, GT, I, MJK, MK.

References

  1. Allen, B. 2023, Phys. Rev. D, 107, 043018 [NASA ADS] [CrossRef] [Google Scholar]
  2. Arzoumanian, Z., Baker, P. T., Blecha, L., et al. 2023, ApJ, 951, L28 [NASA ADS] [CrossRef] [Google Scholar]
  3. Babak, S., & Sesana, A. 2012, Phys. Rev. D, 85, 044034 [NASA ADS] [CrossRef] [Google Scholar]
  4. Babak, S., Petiteau, A., Sesana, A., et al. 2016, MNRAS, 455, 1665 [NASA ADS] [CrossRef] [Google Scholar]
  5. Bécsy, B., & Cornish, N. J. 2020, Class. Quant. Grav., 37, 135011 [CrossRef] [Google Scholar]
  6. Bécsy, B., Cornish, N. J., & Digman, M. C. 2022a, Phys. Rev. D, 105, 122003 [CrossRef] [Google Scholar]
  7. Bécsy, B., Cornish, N. J., & Kelley, L. Z. 2022b, ApJ, 941, 119 [CrossRef] [Google Scholar]
  8. Chalumeau, A., Babak, S., Petiteau, A., et al. 2022, MNRAS, 509, 5538 [Google Scholar]
  9. Chamberlin, S. J., Creighton, J. D. E., Siemens, X., et al. 2015, Phys. Rev. D, 91, 044048 [NASA ADS] [CrossRef] [Google Scholar]
  10. Chen, S., Caballero, R. N., Guo, Y. J., et al. 2021, MNRAS, 508, 4970 [NASA ADS] [CrossRef] [Google Scholar]
  11. Corbin, V., & Cornish, N. J. 2010, arXiv e-prints [arXiv:1008.1782] [Google Scholar]
  12. Cornish, N. J., & Sesana, A. 2013, Class. Quant. Grav., 30, 224005 [NASA ADS] [CrossRef] [Google Scholar]
  13. Ellis, J., & van Haasteren, R. 2017, https://doi.org/10.5281/zenodo.1037579 [Google Scholar]
  14. Ellis, J. A., Siemens, X., & Creighton, J. D. E. 2012, ApJ, 756, 175 [NASA ADS] [CrossRef] [Google Scholar]
  15. Ellis, J. A., Vallisneri, M., Taylor, S. R., & Baker, P. T. 2020, https://doi.org/10.5281/zenodo.4059815 [Google Scholar]
  16. EPTA Collaboration and InPTA Collaboration (Antoniadis, J., et al.) 2023a, A&A, 678, A50 [CrossRef] [EDP Sciences] [Google Scholar]
  17. EPTA Collaboration and InPTA Collaboration (Antoniadis, J., et al.) 2023b, A&A, 678, A49 [Google Scholar]
  18. EPTA Collaboration and InPTA Collaboration (Antoniadis, J., et al.) 2024, A&A, 685, A94 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  19. Falxa, M., Babak, S., Baker, P. T., et al. 2023, MNRAS, 521, 5077 [NASA ADS] [CrossRef] [Google Scholar]
  20. Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759 [Google Scholar]
  21. Hee, S., Handley, W. J., Hobson, M. P., & Lasenby, A. N. 2016, MNRAS, 455, 2461 [NASA ADS] [CrossRef] [Google Scholar]
  22. Jordan, A., Blakeslee, J. P., Cote, P., et al. 2007, ApJS, 169, 213 [NASA ADS] [CrossRef] [Google Scholar]
  23. Karnesis, N., Katz, M. L., Korsakova, N., Gair, J. R., & Stergioulas, N. 2023, arXiv e-prints [arXiv:2303.02164] [Google Scholar]
  24. Marco, V. D., Zic, A., Miles, M. T., et al. 2023, ApJ, accepted [arXiv:2305.04464] [Google Scholar]
  25. Nowak, N., Saglia, R. P., Thomas, J., et al. 2008, MNRAS, 391, 1629 [NASA ADS] [CrossRef] [Google Scholar]
  26. Skilling, J. 2004, AIP Conf. Ser., 735, 395 [Google Scholar]
  27. Speri, L., Porayko, N. K., Falxa, M., et al. 2023, MNRAS, 518, 1802 [Google Scholar]
  28. Taylor, S. R., Huerta, E. A., Gair, J. R., & McWilliams, S. T. 2016, ApJ, 817, 70 [NASA ADS] [CrossRef] [Google Scholar]
  29. Taylor, S. R., Lentati, L., Babak, S., et al. 2017, Phys. Rev. D, 95, 042002 [NASA ADS] [CrossRef] [Google Scholar]
  30. Taylor, S. R., Baker, P. T., Hazboun, J. S., Simon, J., & Vigeland, S. J. 2021, https://github.com/nanograv/enterprise_extensions [Google Scholar]
  31. van Haasteren, R., & Vallisneri, M. 2014, Phys. Rev. D, 90, 104012 [NASA ADS] [CrossRef] [Google Scholar]
  32. Verbiest, J. P. W., Weisberg, J. M., Chael, A. A., Lee, K. J., & Lorimer, D. R. 2012, ApJ, 755, 39 [NASA ADS] [CrossRef] [Google Scholar]
  33. Vigeland, S. J., Islo, K., Taylor, S. R., & Ellis, J. A. 2018, Phys. Rev. D, 98, 044003 [NASA ADS] [CrossRef] [Google Scholar]
  34. Zonca, A., Singer, L., Lenz, D., et al. 2019, J. Open Source Softw., 4, 1298 [Google Scholar]

Appendix A: Simulation

The simulations were performed using the fakepta package (https://github.com/mfalxa/fakepta). The injected noises are assumed to be stationary gaussian processes with a power spectral density Sn(f, α) function of frequency f and hyper-parameters α according to EPTA Collaboration and InPTA Collaboration (2023b). In time domain, such noises n(t) can be decomposed on a Fourier basis of size N:

n ( t ) = k = 1 N X k sin ( 2 π k t / T obs ) + Y k cos ( 2 π k t / T obs ) = F a , $$ \begin{aligned} n({t}) = \sum _{k=1} ^N X_k \sin (2\pi k {t} / T_{obs}) + Y_k \cos (2\pi k {t} / T_{obs}) = \mathbf F {a}, \end{aligned} $$(A.1)

where we rewrote the sum as a matrix-vector multiplication with a a vector of (Xk, Yk) and F the design matrix of size NTOAs × N containing the sine and cosine terms sin(2πkt/Tobs),cos(2πkt/Tobs) per each TOA.

The covariance matrix Cn(t, t′) of the gaussian process n(t) is given by the expectation value ⟨⟩ of the noise at two times:

C n ( t , t ) = n ( t ) , n ( t ) = F a 2 F = F Φ F $$ \begin{aligned} \mathbf C _n ({t}, {t^{\prime }}) = \langle n({t}), n({t^{\prime }}) \rangle = \mathbf F ^\intercal \langle {a}^2 \rangle \mathbf F = \mathbf F ^\intercal {\Phi } \mathbf F \end{aligned} $$(A.2)

with Φ = diag{Sn(k/Tobs)/Tobs}.

Considering that the stochastic process n(t) follows a multivariate gaussian distributions 𝒩(0, Cn), a random realisation of n(t) corresponds to a random draw from the distribution 𝒩(0, Cn). The noises are injected by summing the obtained n(t) to the simulated pulsar residuals for the maximum a posteriori set of noise hyper-parameters α (entering in the Sn) inferred from the real data.

All Tables

Table 1.

Statistical significance of the candidate source at 4.64 nHz.

Table 2.

Parameters of the continuous gravitational wave model with their respective priors and ranges.

Table 3.

Bayes factors obtained for different models.

All Figures

thumbnail Fig. 1.

e-statistic of the candidate source at fgw = 4.64 nHz averaged over the noise uncertainties for the custom PSRN model. The black star shows the position of highest ℱe, whereas the red stars show the positions of the pulsars. The Fornax and Virgo clusters are shown as black dots.

In the text
thumbnail Fig. 2.

Distribution of the ℱe-statistic over the noise uncertainties without CURN (blue) and with CURN (orange) at 4.64 nHz. The null distributions of the ℱe are obtained from the analysis of the EPTA DR2new data with shuffled sky positions (grey shaded region) and from the theoretical formula of a χ42-distribution (black solid line).

In the text
thumbnail Fig. 3.

Pulsar-shuffled version of the ℱe-statistic of the candidate source at fgw = 4.64 nHz averaged over the noise uncertainties for the custom PSRN model. The black star shows the position of highest ℱe, whereas the red stars show the positions of the pulsars. We see that the highest ℱe pixel changed position but is still significant. The Fornax and Virgo clusters are shown as black dots.

In the text
thumbnail Fig. 4.

Bayes factor for the model comparison PSRN+CURN+CGW (Earth term) versus PSRN+CURN for 50 logarithmically spaced frequency sub-bands in the region fgw ∈ [1.5, 320] nHz.

In the text
thumbnail Fig. 5.

Posterior distributions of the CGW search in the second data release of the EPTA DR2new. The posteriors are obtained using a CGW model with Earth and pulsar term, the custom PSRN model and a CRN with either CURN or HD (GWB) correlations represented by 9 frequency bins. We show the posterior distribution for the gravitational wave frequency and amplitude fgw, h of the CGW, and the common noise spectral index and amplitude γ and A. The contours indicate the 1,2,3-σ Gaussian contours.

In the text
thumbnail Fig. 6.

Posterior distribution of the sky localisation obtained by searching for a CGW in the second data release of the EPTA DR2new. The posteriors are obtained using a CGW model with Earth and pulsar term, the inclusion of the custom pulsar noise (PSRN) and a common uncorrelated red noise (CURN) represented by 9 frequency bins. For reference, we show the position of the analysed pulsars and the Virgo and Fornax clusters.

In the text
thumbnail Fig. 7.

Horizon luminosity distance, DH, obtained from the sky averaged 95% upper limit on strain amplitude h using the PSRN+CURN+CGW (Earth term + pulsar term) model. The horizon distance is calculated with Eq. (10) for three chirp masses: 108M, 109M and 1010M. For reference, the Virgo and Fornax clusters are at 16.5 Mpc and 19.3 Mpc from the Earth, respectively.

In the text
thumbnail Fig. 8.

Posterior distributions of the CGW search in the second data release of the EPTA DR2new. The posteriors are obtained using a CGW model with Earth and pulsar term, the custom PSRN model and an HD-correlated background (GWB) represented by 9 and 3 frequency bins. We show the posterior distribution for the gravitational wave frequency and amplitude fgw, h of the CGW, and the common noise spectral index and amplitude γ and A. The contours indicate the 1,2,3-σ Gaussian contours.

In the text
thumbnail Fig. 9.

Inference of the amplitude h and frequency fgw of CGW using the PSRN+CURN+CGW model. The results obtained with Eryn are shown as green and blue histograms, for a Earth term only, and full CGW model, respectively. The blue contours are to be compared with the orange posterior obtained with QuickCW. The shown contours are the 90% normal credible regions.

In the text
thumbnail Fig. 10.

Distributions of the signal-to-noise-ratio (S/N) for the HD correlations for a common red process with γ = 13/3, with (dashed) or without (solid) adding CGW to the data model. The orange lines correspond to the results obtained on DR2new and blue are on the simulated data with only CGW (no CRN).

In the text
thumbnail Fig. 11.

Posterior distribution of the gravitational wave frequency fgw of a CGW fitted to two simulated PTAs: one with a single injected CGW (PSRN+CGW blue histogram), and one with an injected GWB (PSRN+GWB orange histogram). The posterior distribution is obtained with a MCMC analysis for a PSRN+CGW (Earth term) model. The dashed lines are the medians of the distributions.

In the text
thumbnail Fig. 12.

Posterior distribution of amplitude log10A and spectral index γ of the GWB for two simulated PTAs: one realistic PTA with only one CGW injected (PSRN+CGW) and one realistic PTA with a GWB injected (PSRN+GWB). The posterior distribution is obtained with a MCMC analysis for a PSRN+GWB model. The dashed lines are the medians of the distributions.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.