Description and Application of the Surfing Effect

Abstract

The standard technique for very low-frequency gravitational wave detection is mainly based on searching for a specific spatial correlation in the variation of the times of arrival of the radio pulses emitted by millisecond pulsars with respect to a timing model. This spatial correlation, which in the case of the gravitational wave background must have the form described by the Hellings and Downs function, has not yet been observed. Therefore, despite the numerous hints of a common red noise in the timing residuals of many millisecond pulsars compatible with that expected for the gravitational wave background, its detection has not yet been achieved. By now, the reason is not completely clear and, from some recent works, the urgency to adopt new detection techniques, possibly complementary to the standard one, is emerging clearly. Of course, this demand also applies to the detection of continuous gravitational waves emitted by supermassive black hole binaries populating the Universe. In the latter case, important information could, in principle, emerge from the millisecond pulsars considered individually in a single-pulsar search of continuous GWs. In this context, the surfing effect can then be exploited, helping to select the best pulsars to carry out such analysis. This paper aims to clarify when the surfing effect occurs and describe it exhaustively. A possible application to the case of the supermassive black hole binary candidate PKS 2131–021 and millisecond pulsar J2145–0750 is also analyzed.

1. Introduction

Since the first direct observation of gravitational waves (GWs), carried out by the interferometers of the LIGO/VIRGO collaboration [1], gravitational radiation has proved to be a precious messenger of information for astrophysics, allowing not only to confirm the predictions of general relativity [2,3] but also to study extraordinary phenomena, such as the merger of neutron stars [4], and to discover several black holes [5]. However, gravitational astronomy is still in its infancy and to date, only the tip of the iceberg of the GW spectrum has been scratched. Additionally, low-frequency GWs are not accessible by ground-based interferometers, which are currently only sensitive to high-frequency GWs (in the range 10–${10}^4$ Hz) emitted in the final stages of the coalescence of compact objects [6]. To explore the low-frequency side of the GW spectrum (in the range ${10}^{-5}$–1 Hz), space-based gravitational interferometers, such as Laser Interferometer Space Antenna (LISA) [7], are required and, in addition, for very low-frequency GWs (in the range ${10}^{-10}$–${10}^{-6}$ Hz), pulsar timing arrays (PTAs) must be used, which are currently the only instruments sensitive to them [8].

PTAs are arrays of extremely regular millisecond pulsars (MSPs), constantly monitored by radio telescopes to measure the variations in the times of arrival (ToAs) of the emitted radio pulses, commonly referred to as timing residuals (TRs). Since GWs can influence the radio pulse ToAs, their signature can be found by studying the TRs [9,10,11], once cleaned from other effects that can affect them. At present, the main purpose of PTAs is to reveal the gravitational wave background (GWB) due to the superposition of GWs emitted by the supermassive black hole binary systems (SMBHBs) that, according to the current cosmological model, populate the Universe [12]. The GWB should induce spatially correlated TRs in all the MSPs of PTAs [13]. Since this correlation should take the form described by the Hellings and Downs function [14], the latter is considered the smoking gun for the GWB [15]. Although from the data collected by the main PTA collaborations (the European Pulsar Timing Array (EPTA) [16], the Indian Pulsar Timing Array (InPTA) [17], the North American Nanohertz Observatory for Gravitational Waves (NANOGrav) [18], and the Parkes Pulsar Timing Array (PPTA) [19], which join their efforts as the International Pulsar Timing Array (IPTA) [20]) it has been possible to obtain strong constraints on both the frequency and the amplitude of the GWB, and despite the emergence from them of a common red noise compatible with that expected for the GWB, the smoking gun is still missing, and the first GW detection by PTAs has not yet been claimed. Furthermore, the first detection of continuous GW emitted by a SMBHB seems even more distant, as it may be subject to knowledge of the GWB.

This article is motivated by the exigency to look for additional instruments and tools for detecting very low-frequency GWs [21]. This necessity is made more urgent by recent research works, which show the limits and criticalities of the investigation principles adopted up to now [22]. For this purpose, the potential of the surfing effect for the detection of the continuous GWs emitted by a SMBHB has been analyzed. Although the surfing effect is present in some works on PTAs, the literature to date is lacking in articles dedicated to its discussion in the context of general relativity, and in which a possible concrete application is presented. For this reason, this paper is structured as follows: Section 2 is dedicated to a detailed description of the surfing effect, Section 3 describes an application of the surfing effect to the SMBHB candidate PKS 2131–021, in Section 4 the signatures of the GWs possibly emitted by the SMBHB candidate PKS 2131–021, have been searched in the averaged narrowband and wideband post-fit whitened TRs obtained by NANOGrav collaboration on the MSP J2145–0750, and finally, the main conclusions are presented Section 5.

2. An Overview on the Surfing Effect

The surfing effect occurs when the angle $\theta$ seen by an observer on Earth between the travel direction of the radio pulses emitted by a MSP and the travel direction of the GWs emitted by a single GW source is sufficiently small. In this case, the radio pulses behave as if they were surfing the GWs, producing an increased value of the pulsar TRs. That can be shown by considering the function $R\left(t,n_i\right)$ which describes the pulsar TRs induced by the GWs emitted by a single source [9,10,11]:

$$R\left(t,n_i\right)=\sum _{A=+,\times }{F}^A\left[r_e^A\left(t,n_i\right)-r_p^A\left(t,n_i\right)\right]$$

(1)

where t is the time, $n_i$ is the versor oriented toward the source of the GWs1, A is the GW polarization state index, which is + in the case of the plus-polarization state, and × in the case of the cross-polarization state, and $F^A$ is the antenna pattern function, defined as2:

$$F^A=\frac{1}{2}\frac{p^i{p}^j{e}_{ij}^A}{(1-n_i{p}^i)}$$

(2)

Here, $p^i$ is the versor oriented toward the MSP, and $e_{ij}^A$ is the GW polarization tensor:

$$e_{ij}^{+}=\left(\begin{array}{ccc}1& 0& 0\\ 0& -1& 0\\ 0& 0& 0\end{array}\right)$$

(3)

$$e_{ij}^{\times }=\left(\begin{array}{ccc}0& 1& 0\\ 1& 0& 0\\ 0& 0& 0\end{array}\right)$$

(4)

Moreover, $r_e^A\left(t,n_i\right)$ and $r_p^A\left(t,n_i\right)$ are the Earth and the pulsar terms defined respectively by3:

$$r_e^A\left(t,n_i\right)={\int }_0^{t}t{h}_e^A\left({\omega }_{e}t\right)$$

(5)

$$r_p^A\left(t,n_i\right)={\int }_0^{t}t{h}_p^A\left({\omega }_{p}t-{\omega }_p{t}_p\left(1-n_i{p}^i\right)\right)$$

(6)

where ${\omega }_e$ and ${\omega }_p$ are the GW angular frequency, evaluated at the Earth and MSP positions, respectively, $t_p$ is the time distance from the MSP to the Earth, and $h_e^A\left({\omega }_{e}t\right)$ and $h_p^A\left({\omega }_{p}t-{\omega }_p{t}_p\left(1-n_i{p}^i\right)\right)$ are the perturbations on the flat space-time metric induced by the GWs, evaluated at the Earth and MSP positions, respectively4.

Equation (1) is particularly useful to describe the pulsar TRs induced by the GWs emitted by a circularized SMBHBs in the continuous emission regime. In the literature, this assumption is often made in order to simplify the description of the principles on which PTAs are based. That is reasonable because, according to the most widely accepted models, SMBHBs form during the collision of galaxies by dynamic friction, which, as it is widely known, tends to circularize their orbits [23].

In this case, it is convenient to adopt a coordinate system $Oxyz$ with its origin on the Earth5, so that:

$$n_i\equiv \left(0,0,1\right)$$

(7)

$$p^i\equiv \left(sin{\theta }_{p}cos{\varphi }_p,sin{\theta }_{p}sin{\varphi }_p,cos{\theta }_p\right)$$

(8)

where ${\theta }_p$ is the angle between the versor $p^i$ and the positive z-axis, and ${\varphi }_p$ is the angle between the x-axis and the projection on the $xy$-plane of the versor $p^i$ (see Figure 1). With that choice, the product $n_i{p}^i$, which appears in Equations (2), (5), and (6), is simply $cos{\theta }_p$.

The antenna pattern function in Equation (2) can be rewritten in this coordinate system, taking into account the Equations (3) and (4), for both the polarization states:

$$F^{+}=\frac{1}{2}\frac{{sin}^2{\theta }_p}{1-cos{\theta }_p}cos\left(2{\varphi }_p\right)$$

(9)

$$F^{\times }=\frac{1}{2}\frac{{sin}^2{\theta }_p}{1-cos{\theta }_p}sin\left(2{\varphi }_p\right)$$

(10)

Moreover, the perturbations on the flat space-time metric induced by the GWs evaluated at the Earth and MSP positions are [24]:

$$h_e^A\left({\omega }_{e}t\right)={\mathcal{F}}^A\mathcal{A}sin\left({\omega }_{e}t+{\alpha }^A\right)$$

(11)

$$h_p^A\left({\omega }_{p}t-{\omega }_p{t}_p\left(1-cos{\theta }_p\right)\right)={\mathcal{F}}^A\mathcal{A}sin\left({\omega }_{p}t-{\omega }_p{t}_p\left(1-cos{\theta }_p\right)+{\alpha }^A\right)$$

(12)

where ${\mathcal{F}}^A$ is a factor that accounts for the orbital inclination angle $\iota$ of the SMBHB [25]:

$${\mathcal{F}}^{+}=1+{cos}^2\iota$$

(13)

$${\mathcal{F}}^{\times }=-2cos\iota$$

(14)

and $\mathcal{A}$ is the perturbation amplitude, given by:

$$\mathcal{A}=2\frac{{\mathcal{M}}^{5/3}}{D}{\left(\frac{\omega }{2}\right)}^{2/3}$$

(15)

where $\mathcal{M}$ is the chirp mass of the SMBHB, D is the luminosity distance, and $\omega$ is the angular frequency of GWs, and ${\alpha }^A$ is the initial phase, with ${\alpha }^{+}=0$ and ${\alpha }^{\times }=\pi /2$. Note that all the quantities in Equation (15) are red-shifted [24].

The difference between the Earth term, in Equation (5), and the pulsar term, in Equation (6), can be determined for both the polarization states by taking into account Equations (11)–(15), obtaining:

$$\begin{array}{cc}\hfill & r_e^{+}\left(t,n_i\right)-r_p^{+}\left(t,n_i\right)=-\frac{\mathcal{A}}{{\omega }_e}\left(1+{cos}^2\iota \right)\left[cos\left({\omega }_{e}t\right)-1\right]+\hfill \\ \hfill & +\frac{\mathcal{A}}{{\omega }_p}\left(1+{cos}^2\iota \right)\left[cos\left({\omega }_{p}t-{\omega }_p{t}_p\left(1-cos{\theta }_p\right)\right)-cos\left(-{\omega }_p{t}_p\left(1-cos{\theta }_p\right)\right)\right]\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill & r_e^{\times }\left(t,n_i\right)-r_p^{\times }\left(t,n_i\right)=-\frac{\mathcal{A}}{{\omega }_e}\left(-2cos\iota \right)\left[-sin\left({\omega }_{e}t\right)\right]+\hfill \\ \hfill & +\frac{\mathcal{A}}{{\omega }_p}\left(-2cos\iota \right)\left[-sin\left({\omega }_{p}t-{\omega }_p{t}_p\left(1-cos{\theta }_p\right)\right)+sin\left(-{\omega }_p{t}_p\left(1-cos{\theta }_p\right)\right)\right]\hfill \end{array}$$

(17)

Then, the full expression of the function $R\left(t,n_i\right)$ is found by substituting the Equations (9), (10), (16) and (17) in Equation (1). The result obtained above can be simplified by introducing some additional working hypotheses. First of all, the SMBHB can be assumed to be face-on, which means that $\iota =\pi /2$. In this case, as can be seen from Equations (13) and (14), ${\mathcal{F}}^{+}=1$ and ${\mathcal{F}}^{\times }=0$. Secondly, since only one MSP is under consideration, the angle ${\theta }_p$ coincides with the angle $\theta$, so ${\theta }_p=\theta$, the coordinate system is rotated to have ${\varphi }_p=0$, and $t_p$ is replaced by T to simplify the notation. In this case, as can be seen by Equations (9) and (10), $F^{+}={sin}^2\theta {\left[2\left(1-cos\theta \right)\right]}^{-1}$ and $F^{\times }=0$. Thirdly, the SMBHB can be assumed to be very far from its coalescence, as most of the SMBHBs observable with PTAs should be. Under this assumption, the variation of the angular frequency of the GWs [26] during the interval of time $\Delta t=T\left(1-cos\theta \right)$, which denotes the time delay between the Earth epoch and the pulsar epoch [27], can be neglected, so ${\omega }_e\approx {\omega }_p$ and the angular frequency of GWs can be denoted by just $\omega$. Using all these assumptions, one gets:

$$\begin{array}{cc}\hfill & R\left(t,\theta \right)=-\frac{1}{2}\frac{{sin}^2\theta }{1-cos\theta }\frac{\mathcal{A}}{\omega }\left[cos\left(\omega t\right)-1\right]+\hfill \\ & -\frac{1}{2}\frac{{sin}^2\theta }{1-cos\theta }\frac{\mathcal{A}}{\omega }\left[-cos\left(\omega t-\omega T\left(1-cos\theta \right)\right)+cos\left(-\omega T\left(1-cos\theta \right)\right)\right]\hfill \end{array}$$

(18)

where the dependence on $n_i$ has been replaced by the dependence on $\theta$. As can be seen from Equation (18), the way the function $R\left(t,\theta \right)$ depends on the angle $\theta$ makes its use for the detection of GWs rather problematic. For example, it is interesting to note that $R\left(t,\theta \right)$ has a discontinuity for $\theta =0$, which can be removed by placing $R\left(t,0\right)\equiv {lim}_{\theta \to 0}R\left(t,\theta \right)=0$. Moreover, for values of $\omega$ and T of the same order of magnitude of the average angular frequency of the GWs emitted by SMBHB observable by PTAs and of the average time distance of the MSPs included in PTAs, respectively, the function $R\left(t,\theta \right)$ oscillates extremely rapidly with $\theta$.

Since the function $R(t,\theta )$ carries all the information about the GWs emitted by the SMBHB, it is crucial to have a clear understanding of its behavior. The easiest way to do this is by considering, instead of $R\left(t,\theta \right)$, the function $S(\zeta ,\eta )$, which, using the substitutions:

$$\zeta \equiv \chi \frac{t}{T}$$

(19)

$$\eta \equiv \chi \left(1-cos\theta \right)$$

(20)

where $\chi \equiv \omega T$, $\zeta \in \left[0,\infty \right[$ for $t\in \left[0,\infty \right[$ and $\eta \in \left[0,2\chi \right]$ for $\theta \in \left[0,\pi \right]$, can be defined by:

$$S\left(\zeta ,\eta \right)\equiv 2\frac{\chi }{\mathcal{A}T}\left|R\left(\zeta ,\eta \right)\right|$$

(21)

The function defined in Equation (21) can be written more explicitly as:

$$S\left(\zeta ,\eta \right)=\left(2-\frac{\eta }{\chi }\right)\left|cos\zeta -1-cos\left(\zeta -\eta \right)+cos(-\eta )\right|$$

(22)

where the first factor in Equation (22) has taken out the absolute value since it is always positive in the defined domain. The TRs induced by GWs are enhanced for the values which maximize the function $S\left(\zeta ,\eta \right)$ in Equation (22). The function $S\left(\zeta ,\eta \right)$ is a periodic function characterized by an amplitude that depends only on $\eta$ through the first linear factor:

$$a\left(\eta \right)\equiv 2-\frac{\eta }{\chi }$$

(23)

and its shape depends on both the variables $\zeta$ and $\eta$ through the second oscillating factor:

$$f\left(\zeta ,\eta \right)\equiv \left|cos\zeta -1-cos\left(\zeta -\eta \right)+cos(-\eta )\right|$$

(24)

so that:

$$S\left(\zeta ,\eta \right)=a\left(\eta \right)f\left(\zeta ,\eta \right)$$

(25)

Equation (25) allows to find the global maximum of the function $S\left(\zeta ,\eta \right)$ by studying the maxima of the functions $a\left(\eta \right)$ and $f\left(\zeta ,\eta \right)$. The function $a\left(\eta \right)$ is maximized when:

$$a\left(\eta \right)=2\mathrm{for}\eta =0$$

(26)

that occurs when:

$$\theta =0$$

(27)

while the function $f\left(\zeta ,\eta \right)$ is maximized when:

$$f\left(\zeta ,\eta \right)=4\mathrm{for}\left(\zeta ,\eta \right)=\left(\pi +2\pi n_{\zeta },\pi +2\pi n_{\eta }\right)$$

(28)

that occurs when:

$$t=\frac{\pi T}{\chi }\left(1+2n_{\zeta }\right)\equiv t_{n_{\zeta }}$$

(29)

$$\theta =arccos\left(1-\frac{\pi }{\chi }\left(1+2n_{\eta }\right)\right)\equiv {\theta }_{n_{\eta }}$$

(30)

Therefore, the results in Equations (27) and (30) imply that, for each value of t, the global maximum of the function $R\left(\zeta ,\eta \right)$ can be found placing $n_{\eta }\equiv 0$. Then, since for $n_{\eta }=0$ results ${\theta }_0\ll 1$, this property is referred to as surfing effect, and ${\theta }_0=arccos\left(1-\pi /\chi \right)$ is the surfing effect angle.

3. An Application of the Surfing Effect to the Case of the Supermassive Black Hole Binary Candidate PKS 2131–021

The existence of SMBHBs is not only predicted by most of the hierarchical structure formation models [28] but also supported by observational evidence. In some rare cases, the inferred orbital parameters suggest that, if the data analysis is correct, the SMBHB should emit GWs observable with PTAs. In these cases, as has already been conducted with the SMBHB candidate in the Seyfert galaxy 3C 66B, the data collected by PTAs can be used to rule out the GW emission [27,29,30].

An interesting recent study on the blazar PKS 2131–021 opens the possibility for an application of the surfing effect [31]. In fact, according to the data analysis of the blazar PKS 2131–021 radio emission, it can be identified as an SMBHB candidate located at a redshift $z=1.285$ [32], with an observed orbital period of $1760.4_{-5.3}^{+5.3}$ d and a non-red-shifted chirp mass $\lesssim 5.4\times {10}^9{M}_{\odot }$ (see ref. [31]). Adopting the most recent values for cosmological constants determined by the Planck collaboration [33], the luminosity distance of the SMBHB candidate is $D=9.2$ Gpc, while the angular frequency and the amplitude of the expected GW emission are $\omega =8.{262}_{-0.025}^{+0.025}\times {10}^{-8}$ Hz and $\mathcal{A}\lesssim 2.4\times {10}^{-15}$, respectively.

An optimal strategy to confirm these results, by taking advantage of the surfing effect, is to perform a single-pulsar search of continuous GWs focusing on the MSP with the smallest angular distance from the SMBHB candidate PKS 2131–021 among the ones currently included in PTAs. For this reason, the best MSP is the MSP J2145–0750, which lies at the angle $\theta =0.1156$, corresponding to about $6.{6234}^{\circ }$, from the SMBHB candidate PKS 2131–021 [34,35]. The MSP J2145–0750 is a well-studied MSP, lying in a binary system with a white dwarf, at a distance of $0.{62}_{-0.02}^{+0.00}$ kpc [35] from the Earth and timed for about $12.5$ years by the NANOGrav collaboration [36,37]. The choice of the MSP J2145–0750 is also convenient because it allows neglecting the GW angular frequency variation between the Earth and the pulsar terms due to the SMBHB PKS 2131–021 orbital evolution. In fact, by considering the equation describing the GW angular frequency variation [26]:

$$\dot{\omega }=\frac{12}{5}{2}^{1/3}{\mathcal{M}}^{5/3}{\omega }^{11/3}\left(1+z\right)$$

(31)

and substituting in the integral of Equation (31) the time delay between the Earth epoch and the pulsar epoch, which is $\Delta t=13.5\pm 0.4$ yr, the GW angular frequency variation between the Earth and the pulsar terms turns out to be $\Delta \omega \lesssim 0.042\times {10}^{-8}$ Hz.

With all these data in hand, Equation (30) can be used to determine the function $S\left(\zeta ,\eta \right)$ and then to find out the surfing effect angle (see Figure 2).

Once t has been chosen to satisfy Equation (29), by setting, for example, $n_{\zeta }\equiv 0$, it is possible to define the function $R\left(\theta \right)\equiv R\left(t_0,\theta \right)$, useful for verifying whether the angular separation $\theta$ between the MSP and the SMBHB candidate is such as to maximize the TRs. In this case, $\omega =8.{262}_{-0.025}^{+0.025}\times {10}^{-8}$ Hz and $T={2024}_{-65}^{+16}$ yr, therefore $\chi ={5274}_{-186}^{+58}$, and at the angle ${\theta }_0=0.0345$, corresponding to about $1.{9778}^{\circ }$, one has $R\left({\theta }_0\right)\lesssim 0.116$ μs. Even if the actual angle between the SMBHB candidate PKS 2131–021 and the MSP J2145–0750, which is $\theta =0.1156$, is larger than ${\theta }_0$, it can be found that this is very close to ${\theta }_5=0.1145$, which identifies the sixth peak of the function $R\left(\theta \right)$. In fact, for $\theta =0.1156$ results $R\left(\theta \right)\lesssim 0.104$ μs (see Figure 3 and Figure 4).

We emphasize that this is a lucky coincidence and allows one to say that if the chirp mass of the SMBHB candidate PKS 2131–021 coincides with the estimated upper limit, the TRs induced by the GWs emitted by the SMBHB candidate PKS 2131–021 on the MSP J2145–0750 might be observable even at the current PTA sensitivity level due to the surfing effect. In addition, as can be easily deduced from the Figure 3 and Figure 4, if $\theta$, whether small or large, is in a node of the function $R\left(\theta \right)$, where it vanishes, in the TRs of the considered MSP there will be no trace of the influence of GWs. Eventually, it is necessary to stress the fact that the larger is $\theta$, the higher the precision required for $\chi$. In fact, as it is shown in Figure 5, for the same value of $\theta$ the function $R\left(\theta \right)$ can assume sensibly different values within the uncertainty on $\chi$. However, if $\theta$ is in the first peak, the role played by the error on $\chi$ is marginal.

4. A Look to NANOGrav Data

The NANOGrav collaboration has recently published the $12.5$ yr Data Set, which is publicly available (see ref. [38]), making it possible to examine the TRs of the MSP J2145–0750. In particular, the datasets useful for searching GWs are the averaged narrowband [36] and wideband [37] post-fit whitened TRs (see Figure 6 and Figure 7) since they are characterized by a low weighted root mean squared value ${\sigma }_w$, which has been evaluated by using the following definition:

$${\sigma }_w=\sqrt{\frac{{\sum }_i^n\frac{{\left(x_i-\overline{x}\right)}^2}{{\sigma }_i^2}}{{\sum }_i^n\left(\frac{1}{{\sigma }_i^2}\right)}}$$

(32)

which is valid for any set of quantities $x_1,x_2,\dots ,x_n$ characterized by the errors ${\sigma }_1,{\sigma }_2,\dots ,{\sigma }_n$ and by an average $\overline{x}$. In this case, using Equation (32) it has been found ${\sigma }_w=0.333$ μs for narrowband, and ${\sigma }_w=0.276$ μs for wideband.

In Figure 6 and Figure 7, the origin of time is at Modified Julian Date 53267, corresponding to 19 September 2004, from which the observation started. Equation (18) implies that the TRs induced by the GWs emitted by a SMBHB are periodic, with an angular frequency that is the same as the angular frequency $\omega$ of the GWs. This periodic signature can be searched in TRs using periodic analysis algorithms, such as the Lomb-Scargle (LS) algorithm [39,40].

The LS periodograms of the TRs plotted in Figure 6 and Figure 7 have been obtained and plotted in Figure 8 and Figure 9, respectively.

In both cases (see Figure 8 and Figure 9), among the periodicities with the highest normalized amplitude, the LS algorithm finds one associated with an angular frequency close to $\omega =8.{262}_{-0.025}^{+0.025}\times {10}^{-8}$ Hz, expected for GWs possibly emitted by the SMBHB candidate PKS 2131–021. Specifically, it finds $\omega =8.006\times {10}^{-8}$ Hz in the case of the narrowband dataset and $\omega =8.0139\times {10}^{-8}$ Hz in the case of the wideband dataset. At this point, a shuffling test was conducted to verify that these values are not artifacts and emerge from true periodicities in the datasets. To this aim, several datasets were created by shuffling the starting datasets, and, for each of the shuffled datasets, the LS periodogram was obtained. After doing so, the periodogram associated with the starting dataset was compared with each periodogram associated with the shuffled datasets to establish which is characterized by the largest normalized amplitude at the frequency value under consideration. This procedure was iterated 10,000 times, with the result that the former normalized amplitude was larger with respect to the latter the ≃87% of cases for narrowband and the ≃93% of cases for wideband.

Although it can be read as an encouraging result, it is important to stress that the two periodograms have been obtained using the complete datasets, provided by the NANOGrav collaboration. However, from Figure 6 and Figure 7, it can be noted that, in the first observation years, the cadence of TRs was sporadic and irregular, and the error bands were significantly large than in the last observation years. Therefore, it may be convenient to reanalyze the datasets after applying some filters. Based on the previous considerations, one possibility is to remove both the TRs falling in the first six observation years and those associated with an uncertainty larger than $3{\sigma }_w$ (see Figure 10 and Figure 11). Applying this filter also assures that, as shown in ref. [31], the dataset is limited to a time interval from 2010 onwards when the SMBHB candidate should have been in the continuous emission regime.

The LS periodograms of the TRs plotted in Figure 10 and Figure 11 have been obtained and plotted in Figure 12 and Figure 13, respectively.

In both cases (see Figure 12 and Figure 13), the LS algorithm allows finding the periodicity with the highest normalized amplitude associated with an angular frequency even closer to $\omega =8.{262}_{-0.025}^{+0.025}\times {10}^{-8}$ Hz. Specifically, it is found $\omega =8.147\times {10}^{-8}$ Hz in the case of the narrowband dataset, and $\omega =8.205\times {10}^{-8}$ Hz in the case of the wideband dataset. Even if these values are so close to the angular frequency of the GW emission expected for the SMBHB candidate PKS 2131–021, this is not enough for claiming its detection. Further, more sophisticated analyses must be conducted in order of ensuring the solidity of this result and the nature of this periodicity.

5. Conclusions

In this article, a possible investigation method to search for continuous GWs, based on the surfing effect, has been proposed. The surfing effect, described in Section 2, is usually ignored in most of the studies on the GWB since it is not relevant when the entire array of MSPs is used for its detection [14,41,42,43]. However, when performing a single-pulsar search of GWs emitted by a circularized SMBHBs in the continuous emission regime, it has to be taken into consideration. In fact, in this case, as shown in Section 2 and Section 3, the global maximum of the function $R\left(t,\theta \right)$, which describes the TRs induced by GWs, corresponds, at any time, to an angular separation between the SMBHB and the MSPs almost null. However, it is important to keep in mind that, at any time, the function $R\left(t,\theta \right)$ is characterized by very rapid angular oscillations, the more frequent the greater the $\chi$ factor, as shown in Figure 3 and Figure 4. Therefore, the function $R\left(\theta \right)$ represents a sort of map, which can be used to verify if a MSP is more or less suitable for performing a single-pulsar search of continuous GWs, so MSPs that fall on one of the peaks will be preferred, and those that fall into one of the nodes will be discarded.

As an example, in Section 3 and Section 4, the case of the SMBHB candidate PKS 2131–021 was considered. This system, based on the orbital parameters determined by the radio analysis of its emission [31], should be an emitter of GWs satisfying the hypotheses made in Section 2 for surfing effect. By a lucky coincidence, as shown in Figure 3 and Figure 4, among the MSPs currently included in PTAs, the angularly closest MSP to the SMBHB candidate PKS 2131–021, which is the MSP J2145–0750, is such that its angular separation falls into a peak of the function $R\left(\theta \right)$. This makes it worthy of attention since in its TRs a periodic signature induced by the GWs emitted by the SMBHB candidate PKS 2131–021 may be present. For this reason, in Section 5, the LS algorithm was used to search for such periodic signature. The periodograms plotted in Figure 8 and Figure 9 have been obtained from the complete datasets, provided by the NANOGrav collaboration (see Figure 6 and Figure 7), while the periodograms plotted in Figure 12 and Figure 13 have been obtained from the filtered datasets (see Figure 10 and Figure 11), discussed in Section 4. Interestingly, the LS algorithm finds in each dataset a periodicity associated with an angular frequency very close to the angular frequency of the GW emission expected for the SMBHB candidate PKS 2131–021. Although this result does not constitute proof of the SMBHB candidate PKS 2131–021 GW emission, it is certainly worthy of further investigation. A more detailed and rigorous analysis of the MSP J2145–0750 TRs, using the NANOGrav $12.5$ yr Data Set in conjunction with the IPTA Data Release 2 [44], will be proposed in a future paper entirely dedicated to searching for GWs potentially emitted by the SMBHB candidate PKS 2131–021. There, it will also be discussed the use of other MSPs not too angularly far from the SMBHB candidate PKS 2131–021, as well as the possible improvements that will be brought on this kind of search by next-gen PTAs, such as Square-Kilometer Array (SKA) [45,46,47,48].