The NANOGrav 12.5 yr Data Set: A Computationally Efficient Eccentric Binary Search Pipeline and Constraints on an Eccentric Supermassive Binary Candidate in 3C 66B

The radio galaxy 3C 66B has been hypothesized to host a supermassive black hole binary (SMBHB) at its center based on electromagnetic observations. Its apparent 1.05 yr period and low redshift (∼0.02) make it an interesting testbed to search for low-frequency gravitational waves (GWs) using pulsar timing array (PTA) experiments. This source has been subjected to multiple searches for continuous GWs from a circular SMBHB, resulting in progressively more stringent constraints on its GW amplitude and chirp mass. In this paper, we develop a pipeline for performing Bayesian targeted searches for eccentric SMBHBs in PTA data sets, and test its efficacy by applying it to simulated data sets with varying injected signal strengths. We also search for a realistic eccentric SMBHB source in 3C 66B using the NANOGrav 12.5 yr data set employing PTA signal models containing Earth term-only as well as Earth+pulsar term contributions using this pipeline. Due to limitations in our PTA signal model, we get meaningful results only when the initial eccentricity e 0 < 0.5 and the symmetric mass ratio η > 0.1. We find no evidence for an eccentric SMBHB signal in our data, and therefore place 95% upper limits on the PTA signal amplitude of 88.1 ± 3.7 ns for the Earth term-only and 81.74 ± 0.86 ns for the Earth+pulsar term searches for e 0 < 0.5 and η > 0.1. Similar 95% upper limits on the chirp mass are (1.98 ± 0.05) × 109 and (1.81 ± 0.01) × 109 M ☉. These upper limits, while less stringent than those calculated from a circular binary search in the NANOGrav 12.5 yr data set, are consistent with the SMBHB model of 3C 66B developed from electromagnetic observations.


Introduction
Routine detections of gravitational waves (GWs) from coalescing stellar-mass compact object binaries by the ground-based LIGO-Virgo-KAGRA observatories have opened a new window to our Universe (e.g., Abbott et al. 2023).While the ground-based observatories are sensitive to GW signals in the tens of hertz to kilohertz frequency range, pulsar timing array (PTA : Sazhin 1978;Detweiler 1979;Foster & Backer 1990) experiments are sensitive to GWs in the nanohertz frequency range.PTAs achieve this by routinely monitoring ensembles of millisecond pulsars using some of the world's most sensitive radio telescopes to use them as accurate celestial clocks (Hobbs et al. 2019).There are six PTA collaborations active at present: the North American Nanohertz Observatory for Gravitational Waves (NANOGrav: Demorest et al. 2012), the European Pulsar Timing Array (EPTA: Desvignes et al. 2016), the Parkes Pulsar Timing Array (PPTA: Manchester et al. 2013), the Indian Pulsar Timing Array (InPTA: Tarafdar et al. 2022), the MeerKAT Pulsar Timing Array (MPTA: Miles et al. 2023), and the Chinese Pulsar Timing Array (CPTA: Xu et al. 2023).The International Pulsar Timing Array (IPTA: Verbiest et al. 2016) consortium combines data and resources from a subset of these collaborations and fosters scientific discourse to advance the prospects of nanohertz GW detection and post-discovery science.
Recently, evidence for the presence of a low-frequency stochastic GW background (GWB) in their respective PTA data sets was reported independently by NANOGrav (Agazie et al. 2023d), EPTA+InPTA (Antoniadis et al. 2023b), PPTA (Reardon et al. 2023), and CPTA (Xu et al. 2023).This signal manifests as a common-spectrum red noise process with a cross-pulsar spatial correlation that is consistent with the expected Hellings-Downs overlap reduction function (Hellings & Downs 1983).These exciting results have ushered in the era of nanohertz GW astronomy and the upcoming IPTA Data Release 3 is expected to strengthen the prospects of PTA 63 Sloan Fellow. 64NASA Hubble Fellowship: Einstein Postdoctoral Fellow. 65Infinia ML, 202 Rigsbee Avenue, Durham NC, 27701, USA. 66NANOGrav Physics Frontiers Center Postdoctoral Fellow. 67Deceased. 68NSF Astronomy and Astrophysics Postdoctoral Fellow.
Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
science further.The source(s) of this low-frequency GWB remain(s) inconclusive, and the observed GWB is consistent with many different physical processes in our universe, e.g., a population of supermassive black hole binaries (SMBHBs) emitting GWs, cosmic inflation, first-order phase transitions, and cosmic strings (e.g., Afzal et al. 2023;Agazie et al. 2023b;Antoniadis et al. 2023a).Nevertheless, an ensemble of inspiralling SMBHBs present in active galactic nuclei (AGNs) is believed to be the most prominent source of the nanohertz GWB (Burke-Spolaor et al. 2019).
SMBHBs are expected to form via galaxy mergers.Observations suggest that nearly every galaxy contains a supermassive black hole (SMBH) at its center (Richstone et al. 1998).When two galaxies merge, the pair of SMBHs interact with the broader merger environment during its early evolution; dynamical friction with the surrounding stars and gas helps the SMBH pair to efficiently sink toward the center of the merger remnant, eventually forming a bound SMBHB system (Begelman et al. 1980).Many different mechanisms have been proposed over the years to evolve the binary further down to separations 0.1 pc, such as three-body interactions with a dense galactic stellar core and torques from a circumbinary disk, and the exact nature of such dynamic mechanisms is commonly known as the final parsec problem (De Rosa et al. 2019).Once the SMBHBs reach a sub-parsec separation, they inspiral by emitting GWs in the PTA frequency range and eventually merge.Such sources have previously been predicted to be detected first as a stochastic GWB formed via the incoherent superposition of GWs coming from a population of SMBHBs, and then as individual sources that stand out above the background (Sesana et al. 2009;Rosado et al. 2015;Kelley et al. 2017Kelley et al. , 2018;;Mingarelli et al. 2017;Pol et al. 2021;Bécsy et al. 2022).While the detection of the GWB with accurate spectral characterization could inform us about the average properties of the cosmic population of SMBHBs in our universe (Taylor 2021), the detection of continuous GWs from an individual SMBHB will be a strong indicator of SMBHBs' strength of contribution to the observed GWB (Kelley et al. 2019;Charisi et al. 2022).Further, a continuous GW detection can lead to persistent multi-messenger nanohertz GW astronomy if we can identify its electromagnetic counterpart, which is likely to be an AGN.Therefore, it is highly desirable to detect GWs from individual SMBHBs in the PTA data sets in the near future.
Over the years, PTA experiments have searched for and put progressively more stringent constraints on the presence of continuous GWs from individual SMBHBs in their data sets (Jenet et al. 2004;Zhu et al. 2014;Babak et al. 2016;Perera et al. 2018;Aggarwal et al. 2019;Feng et al. 2019;Arzoumanian et al. 2020bArzoumanian et al. , 2021Arzoumanian et al. , 2023;;Agazie et al. 2023e;Antoniadis et al. 2023c;Falxa et al. 2023).These searches can be divided into two different categories: all-sky searches where no prior information about the source of the GWs is taken into consideration (e.g., Babak et al. 2016;Aggarwal et al. 2019), and targeted searches where GWs from a particular SMBHB candidate, identified through electromagnetic observations, is searched for (e.g., Jenet et al. 2004;Arzoumanian et al. 2020b).For a targeted search, the sky location and the luminosity distance of the continuous GW source are fixed to their known values for the SMBHB candidate.Furthermore, the binary period of the candidate, if obtained from electromagnetic observations, can be used to narrow down the prior for the GW frequency.It has been shown that a targeted search can improve both the GW upper limits (Arzoumanian et al. 2020b) and the detection probability (Liu & Vigeland 2021) by a factor of a few as compared to an all-sky search.
Historically, most PTA searches for continuous GWs have been limited to individual SMBHBs in inspiralling quasi-circular orbits (e.g., Aggarwal et al. 2019;Arzoumanian et al. 2020b).However, many proposed solutions for the final parsec problem rely on the binaries entering the nanohertz frequency regime while retaining non-negligible eccentricities (e.g., Ryu et al. 2018).Furthermore, a few SMBHB candidates such as Spikey (Hu et al. 2020) and OJ 287 (Dey et al. 2019) have electromagnetic signatures whose interpretation requires eccentric orbits.Attempts at searching for eccentric SMBHBs in PTA data sets have been stymied by the prohibitive computational cost of evaluating the PTA signals induced by inspiral GWs from such sources.This resulted in most of the past PTA searches for eccentric SMBHB sources being restricted to using a single pulsar observation and suboptimal periodogram-based methods (Jenet et al. 2004;Feng et al. 2019).Falxa et al. (2023) and Antoniadis et al. (2023c) searched for GWs from eccentric SMBHBs to follow up on the candidate sources identified in their circular SMBHB searches, although they did not report any detection or constraint pertaining to such sources.The computational challenges involved in PTA searches for eccentric SMBHBs were investigated and an efficient computational framework for evaluating the PTA signals was presented in Susobhanan (2023), influenced by Susobhanan et al. (2020).A pilot single-pulsar Bayesian all-sky search was also presented in Susobhanan (2023) as a proof of concept.The development of the above framework into a PTA-based eccentricity search was initially explored in Cheeseboro (2021).In this paper, we use the framework of Susobhanan (2023) to perform, for the first time, a full-PTA-scale multi-messenger targeted search for continuous GWs from an eccentric SMBHB in 3C 66B.We carry out this search using the NANOGrav 12.5 yr (NG12.5)narrowband data set (Alam et al. 2020).
3C 66B (R.A.: 02 h 23 m 11 4112, decl.: 42 59 31. 385+  ¢  ) is an elliptical radio galaxy with an estimated redshift of 0.02092.This redshift corresponds to a luminosity distance (D L ) of 94.17 Mpc assuming a flat ΛCDM model with H 0 = 67.66km/(Mpc s), Ω m0 = 0.30966, and Ω Λ0 = 0.68885 (Aghanim et al. 2020).The presence of an SMBHB in 3C 66B was first hypothesized by Sudou et al. (2003) using Very Long Baseline Interferometry (VLBI) measurements.They attributed the apparent elliptical motion of 3C 66B's radio core to the orbital and precessional motion of the smaller black hole's jet.Based on this model, they estimated an orbital period of 1.05 ± 0.03 yr and a chirp mass of 1.3 × 10 10 M ☉ .Given the proximity of this source to our Galaxy, the estimated binary parameters implied a GW amplitude that should have been well within the detection capabilities of pulsar timing experiments.Jenet et al. (2004) used 7 yr of timing data of a single pulsar, PSR B1855+09, obtained using the Arecibo telescope to search for a signal consistent with the proposed orbital period.In the absence of a detection, they put an upper limit of 7 × 10 9 M ☉ on the chirp mass of the binary at zero eccentricity, ruling out the initially proposed model for the VLBI observations.Iguchi et al. (2010) reported flux variations of 3C 66B at 3 mm with a periodicity of 93 days, which they attributed to the Doppler-shifted flux modulation due to the orbital motion of the SMBHB.Assuming the SMBHB to be in a circular orbit with an orbital period of 1.05 yr as in Sudou et al. (2003), they estimated a chirp mass of 7.9 × 10 8 M ☉ using this revised model, an order of magnitude less than the upper bound set by Jenet et al. (2004).The PTA upper limits on the GWs emitted by this source were refined further by Arzoumanian et al. (2020b) and Arzoumanian et al. (2023), and the latest upper limit on the chirp mass of the 3C 66B system is 1.41 × 10 9 M ☉ .This upper limit is calculated using the NG12.5 data set (Alam et al. 2020) while assuming a circular binary and also accounting for the presence of a common-spectrum red noise process in the data.
While Iguchi et al. (2010) assumed a circular orbit to explain the periodic flux variations in 3C 66B, an eccentric SMBHB can also explain all the observed variations that suggest the presence of a binary in 3C 66B.In fact, limits on the maximum possible eccentricity of the postulated binary in 3C 66B for a few fixed chirp masses of the system were also calculated in Jenet et al. (2004) using single-pulsar data.However, a PTAbased search for continuous GWs from an eccentric binary in 3C 66B using multiple pulsars has never been reported in the literature.In this paper, we present a targeted search for PTA signals, induced by an eccentric SMBHB system in the radio galaxy 3C 66B, in NG12.5 data set using a Bayesian framework.For this targeted search, we have used the orbital frequency of the proposed binary measured by Sudou et al. (2003), the known sky location of the radio galaxy core, and the luminosity distance to the source (derived from its measured redshift) as informative priors.
This paper is arranged as follows.In Section 2, we provide a brief overview of how we model the PTA signal induced by GWs from an eccentric SMBHB, including the expressions for GW-induced pulsar timing residuals R(t) (Section 2.1), a description of the general relativistic orbital dynamics of eccentric binaries required to calculate R(t) (Section 2.2), and the PTA signal parameterization for the targeted search (Section 2.3).In Section 3, we describe our data analysis methods, including the data set we used for the GW search (Section 3.1), the pulsar timing and noise models (Section 3.2), the prior distributions for the parameters (Section 3.3), the model comparison method (Section 3.4), and software pipeline used for performing the search (Section 3.5).The application of our method on a simulated data set and the corresponding results are discussed in Section 4. Finally, we present the results of our targeted search for GWs from eccentric binaries in NG12.5 data set in Section 5 along with a summary and discussion in Section 6.

PTA Signal Model for Eccentric Binaries
We begin by briefly describing the PTA signal induced by an SMBHB inspiralling along a relativistic eccentric orbit.The computation of the PTA signal usually involves modeling the relativistic motion of the binary using a post-Newtonian (PN) formalism, where general relativistic corrections to the Newtonian binary dynamics are expressed in powers of (v/c) 2 ∼GM/(c 2 r), 69 where M, r, and v are the total mass, relative separation, and the orbital velocity of the binary, respectively, and c is the speed of light in vacuum (Blanchet 2014).The two GW polarization amplitudes h +,× (t) can be expressed as functions of orbital variables, and the GW-induced pulsar timing residual (the PTA signal) R(t) involves the time integrals of h +,× (t).
The PTA signal originating from such a source is described in Section 2.1 and the corresponding orbital motion is described in Section 2.2, following Susobhanan (2023).Different approaches to modeling such PTA signals may be found in Jenet et al. (2004), Taylor et al. (2016), andSusobhanan et al. (2020).

PTA Signal
When a GW travels across our line of sight to a pulsar, it modulates the observed times of arrival (TOAs) of the pulses from that pulsar.Such modulations are known as the PTA signal and can be expressed for a given pulsar as (Detweiler 1979) where h is the GW strain (see below for precise definition), t and t¢ are coordinate times measured at the solar system barycenter (SSB), t 0 is an arbitrary fiducial time.The geometric time delay Δ p is given by where D p is the distance to the pulsar and μ is the angle between the lines of sight to the pulsar and the GW source.
The GW strain h(t) can be written as a linear combination of the two polarization amplitudes h +,× (t): where ψ is the GW polarization angle and F +,× are the antenna pattern functions that depend on the sky locations of the pulsar and the GW source (see, e.g., Taylor et al. (2016) and the contributions to R(t) arising from s +,× (t) and s +,× (t − Δ p ) are known as the Earth term and the pulsar term, respectively.The quadrupolar-order expressions for s +,× (t), derived assuming slow advance of periapsis and orbital decay (i.e., the timescales of the advance of periapsis and the orbital decay are much longer than the orbital period), are given by Jenet et al. (2004), and Susobhanan (2023) A term appearing at (( ) )  v c n 2 is usually referred to as an nPN correction.
where , where ι is the orbital inclination; ω(t) is the argument of periastron, and e t (t) and u(t) are the time eccentricity and the eccentric anomaly, respectively.The PTA signal amplitude is given by S(t) = H(t)/n(t), where n(t) = 2π/P(t) is the mean motion of the binary and H(t) = GMη x(t)/(D L c 2 ) is the dimensionless GW strain amplitude; P(t) is the binary period, G is the universal gravitational constant, M = m 1 + m 2 is the total mass, η = m 1 m 2 /M 2 is the symmetric mass ratio, m 1 and m 2 are the black hole masses, and D L is the luminosity distance to the binary.The dimensionless PN parameter x is given by ( where k(t) is the advance of pericenter per orbit.M and η relate to the chirp mass M ch as M ch = η 3/5 M. Here, we have divided the time-varying amplitude of the signal into a constant amplitude S 0 = S(t 0 ) and a time-varying part ς(t) = S(t)/S 0 for the convenience of parameterization.The orbital parameters appearing in Equations (5b) and (6c) will be explained in greater detail in the following section.

PN-accurate Orbital Motion
It is evident from Equations (4)-(6c) that we need to calculate the temporal evolution of various binary orbital parameters in order to determine R(t).The orbital evolution of a relativistic binary can be split into two parts: the conservative part that takes into account the advance of periapsis of the orbit, and the reactive evolution due to emission of GWs from the system.We use the PN-accurate quasi-Keplerian parameterization to describe the conservative dynamics of our relativistic eccentric binary systems (Damour & Deruelle 1985;Memmesheimer et al. 2004).We begin our description by defining the mean anomaly l as and expressing the eccentric anomaly u implicitly as a function of l via the PN-accurate Kepler equation (Memmesheimer et al. 2004;Boetzel et al. 2017) , 8 t t

= -+
where ( ) F u t is a periodic function of u that first appears at the 2PN order.The angular coordinate in the orbital plane (orbital phase) f can be written as where f is the true anomaly given by k is the advance of periastron per orbit, ϖ is the periastron angle defined by e f is known as the angular eccentricity, and ( ) F u f is a periodic function of u first appearing at the 2PN order.The argument of periastron ω = f − f, and is not to be confused with ϖ. 70 We use the 3PN accurate expressions for k, e f , F t , and F f involving x, e t , and η found in (Boetzel et al. 2017, supplementary material).
In the absence of galactic environmental influences, the orbital period and eccentricity of the binary decrease with time due to the emission of GWs from the system, and this reactive evolution of the orbit can be incorporated into our framework by allowing n and e t to vary slowly as functions of time.Accurate up to the leading quadrupolar (2.5PN) order, the orbital evolution can be expressed as a system of four coupled differential equations (Damour et al. 2004 In this work, we solve the above system of differential equations using the analytic solution provided in Susobhanan et al. (2020) to get the temporal evolution of n, e t , ϖ, and l, starting with certain initial conditions (n 0 , e 0 , ϖ 0 , l 0 ) defined at t 0 .Thereafter, we use Equations ( 8)-( 10) to calculate the evolution of u and ω as a function of time.Now we have all the inputs to calculate the PTA signal R(t) using Equations (4), (5b), and (6c).

Parameterizing the PTA Signal Model for a Targeted Search
To calculate the PTA signal due to GWs from an individual SMBHB for a particular pulsar, we need information about various quantities related to the GW source and the pulsar.These include the sky position of both the GW source (R.A. gw , decl.gw ) and the pulsar (R.A. psr , decl.psr ), required to calculate the antenna pattern functions F +,× .The other parameters for the eccentric SMBHB GW model, gleaned from Sections 2.1 and 2.2, are the PTA signal amplitude (S 0 ) at a reference time t 0 , the projection angles (ψ, ι), the binary mass parameters (M, η), and the initial condition for the orbital evolution (n 0 , e 0 , ϖ 0 , l 0 ) given at t 0 .Further, the pulsar distance D p also enters the model if the pulsar term contributions are included, where p is an index denoting the pulsar.
In the case of a targeted search, the GW source coordinates and the redshift (z) are usually well known from electromagnetic observations, and the luminosity distance D L to the source can be computed from z by assuming a cosmological model.For some SMBHB candidates (e.g., 3C 66B, OJ 287; Dey et al. 2018), the orbital period P 0 (and therefore n 0 ) at some t 0 is also accurately measured from the electromagnetic data.Additionally, given S 0 , n 0 , e 0 , η, and D L , one can compute M by numerically solving the equation 71 We can also rewrite this equation as where M ch = Mη 3/5 is the chirp mass of the binary.Therefore, the unknown or free independent model parameters related to the GW source for a targeted search include (S 0 , ψ, ι, η, e 0 , ϖ 0 , l 0 ).
Regarding the pulsar-related parameters, the pulsar sky positions are very accurately known from pulsar timing.The pulsar distances are usually not well constrained by electromagnetic observations and hence we treat the pulsar distance for each pulsar as an unknown parameter.In addition, the initial orbital phase parameters l p and ϖ p for the pulsar term are treated as independent from l 0 and ϖ 0 due to the poor constraint on the pulsar distances and the errors arising from restricting the reactive orbital evolution to the leading PN order.To summarize, we use, along with the seven GW source parameters mentioned above, (D p , ϖ p , l p ) for each pulsar as free model parameters for our targeted search.Obviously, the pulsar-specific parameters are omitted in the case of an Earth term-only search, where we consider only the Earth term contributions to the PTA signal and ignore the pulsar term contributions.

Data Analysis Methods
In this section, we begin by describing briefly the pulsar timing data set used in this paper, along with the full model we used to describe the pulsar timing residuals.The priors we have used for different free model parameters are also discussed in detail.The software pipeline used in our analysis is also discussed.

NANOGrav 12.5 yr Narrowband Data Set
As mentioned earlier, we have used the NANOGrav 12.5 yr (NG12.5)narrowband data set to perform our targeted search for continuous GWs from an eccentric SMBHB in 3C 66B.This data set consists of TOAs and timing models for 47 pulsars obtained using the Green Bank Telescope and the Arecibo Telescope from 2004-2017 (Alam et al. 2020).We excluded from our analysis PSRs J1946+3417 and J2322 +2057 due to their short observation span (less than 3 yr), and PSR J1713+0747 due to the presence of two chromatic timing events within the data span that were identified to be of non-GW (possibly pulsar magnetospheric or interstellar medium) origin (Lam et al. 2018).The earliest and latest TOAs in this data set are MJDs 53216 (2004 July 30) and 57933 (2017 June 29) respectively, making the data span ∼12.92 yr.Each TOA is transformed to the SSB frame using the DE438 solar system ephemeris. 72The dispersion measure (DM) variations present in the TOAs are corrected by applying DMX parameters (see Alam et al. 2020), which provide a piecewise constant model thereof.

Model for the Timing Residuals
Pulsar timing residuals represent the deviations in the observed TOAs from the ones predicted by a timing model.The timing residuals r in the vicinity of the maximum likelihood values of the timing model parameters can be expressed as a sum of different components as follows: where M is the pulsar timing design matrix and ò is a vector that represents the deviations in timing model parameters from their maximum likelihood values arising due to the presence of signal components that were unaccounted for in the initial timing model.The WN term n WN is modeled using three phenomenological parameters EFAC, EQUAD, and ECORR, each for every telescope receiver/backend combination.The EFAC parameters scale the TOA measurement uncertainties by a multiplicative factor, the EQUAD parameters add to the measurement uncertainties in quadrature, and the ECORR parameters model the correlated WN between narrowband TOAs derived from the same observation (Agazie et al. 2023c).
The IRN for each pulsar is modeled as a Gaussian red noise process with a power-law spectral density with amplitude A p and spectral index γ p , where p is an index denoting the pulsar and f yr = 1 yr −1 .We model the IRN of a pulsar as a Fourier sum consisting of 30 linearly spaced frequency bins ranging from 1/T span to 30/T span , where T span represents the total observation time span of the data set.
In addition to the IRN, Arzoumanian et al. (2020a) showed the presence of a common-spectrum red noise process n CURN in the NG12.5 data set, albeit without any conclusive detection of spatial correlations.Recently, this common process was shown to be an early hint of the GWB, which manifests as an HDcorrelated common-spectrum process in the NANOGrav 15 yr data set (Agazie et al. 2023d).In this work, we treat n CURN as a spatially uncorrelated common-spectrum process since spatial 71 We apply the Newton-Raphson method.This requires the derivative , where k 0 denotes the value of k at t 0 . 72https://naif.jpl.nasa.gov/pub/naif/JUNO/kernels/spk/de438s.bsp.lblcorrelations are not detected in the NG12.5 data set.We model n CURN as a Gaussian process with a power-law spectrum as where A CURN is the amplitude and γ CURN is the spectral index.We model the CURN as a Fourier sum consisting of five linearly spaced frequency bins ranging from 1/T span to 5/T span to be consistent with the NG12.5 GWB analysis (Arzoumanian et al. 2020a).Overall, the models we use for the WN, IRN, and the CURN are same as those from the NG12.5 GWB analysis.

Priors
For IRN and the CURN parameters, we use priors that are identical to those used in the NG12.5 GWB analysis (Arzoumanian et al. 2020a).This corresponds to log uniform prior for the amplitudes and uniform prior for the spectral indices, namely where A and γ represent the amplitudes and spectral indices appearing in Equations ( 16) and (17).Further, we fix the WN parameters at their maximum likelihood values obtained from single-pulsar analysis as is standard practice in PTA analyses (e.g., Arzoumanian et al. 2023).This amounts to 90 free noise parameters in total: two IRN parameters per pulsar and two CURN parameters.
The prior distributions for the eccentric SMBHB parameters used in our analysis are listed in Table 1 along with the fixed parameters known from the binary model of 3C 66B (Sudou et al. 2003;Iguchi et al. 2010).We fix the sky position of the GW source to the known optical location of 3C 66B.The uncertainty in the sky position of 3C 66B should not significantly affect our results as it is very small compared to the spatial resolution sensitivity of our data.This can be inferred from previous all-sky searches for continuous GWs (e.g., Agazie et al. 2023e;Arzoumanian et al. 2023) where it is seen that the upper limits do not change rapidly with sky position.The orbital period P 0 of the binary in 3C 66B estimated from electromagnetic observations is 1.05 ± 0.03 yr.We treat P 0 as a constant in our analysis because the abovementioned uncertainty in the period corresponds to a GW frequency uncertainty of ∼0.9 nHz, which is smaller than the frequency resolution of ∼2.5 nHz for the NG12.5 data set.The luminosity distance D L of the GW source is estimated using the redshift (z = 0.02092) assuming a flat Universe with cosmological parameters given in Aghanim et al. (2020), available as Planck18 in astropy (Price-Whelan et al.

2018).
We use two types of priors for the PTA signal amplitude S 0 following the standard practice adopted in PTA searches: the detection analysis is performed using a uniform prior on S log 10 0 and the upper limit analysis is performed using a uniform prior on S 0 (denoted as "LinearExp" for S log 10 0 in Table 1).Although η ä (0, 0.25] and e 0 ä [0, 1) by definition, we set their lower bounds to small non-zero numbers to avoid division by zero errors.We also restrict e t < 0.8 because, for large values of e t , the angular eccentricity e f may go above 1 for highly relativistic systems, and this renders the quasi-Keplerian parameterization invalid.This is discussed further in Section 5.The upper bound of ψ is set to π as it enters the PTA signal as sin 2y and cos 2y, and similar priors apply for ϖ 0 and ϖ p also.
To set the priors for D p for each pulsar, we use the pulsar distance obtained from the parallax measurement if available.Otherwise, we use the NE2001 galactic DM model (Cordes & Lazio 2002) to estimate the distance from the DM of the pulsar (see Section 2.3.4 of Arzoumanian et al. (2023) for a detailed discussion).The pulsar distances and uncertainties used in our analysis can be found in Table 2 of Arzoumanian et al. (2023).We use a truncated Normal distribution to ensure that D p > 0, and scale the uncertainties of DM-based D p estimates by a  (2020).The fiducial time t 0 is an approximate reference time for the end of Sudou et al. (2003) observations."PsrDistPrior" represents normal distributions truncated at 0 obtained using parallax and DM-based distance estimates (see the text for details).In addition to the above single-parameter priors, the joint prior distribution of S 0 , η, and e 0 is modified by certain validation criteria ( ) V S e log , , 10 0 0 h as described in the text.
factor of 2 to account for systematic errors inherent in such measurements.This is denoted as "PsrDistPrior" in Table 1.The pulsar sky locations are accurately known from pulsar timing.We treat them as fixed parameters while computing the PTA signals since their uncertainties are too small to affect our results.
In addition to the above single-parameter priors, we set the prior to zero for the combinations of parameters that fail either of the following criteria.
1.The binary should not merge within the data span (see Susobhanan et al. (2020) for an expression for a quadrupolar-order approximation for the merger time).2. The post-Keplerian parameters k < 0.5 and e f < 1 throughout the evolution of the binary within the data span.
As the binary approaches the merger, the quadrupolar approximation becomes inadequate and a more sophisticated inspiral-merger-ringdown description of the PTA signal becomes necessary.Since we lack such a description currently, we must exclude this regime from our analysis.Furthermore, the quasi-Keplerian description breaks down as the binary becomes more and more relativistic, and we exclude this regime from our analysis by employing the somewhat arbitrary condition k < 0.5.For a given P 0 and D L of the SMBHB, the evolution of the binary is uniquely determined by S 0 , e 0 , and η as the total mass M can be calculated using Equation (13).
Therefore, we define a validation function ( ) V S e log , , 10 0 0 h such that it is set to 1 if the combination of S 0 , e 0 , and η values satisfies the criteria mentioned above, otherwise 0. We multiply this validation function with the otherwise uniform prior distributions for S log 10 0 , e 0 , and η, as shown in Table 1, to make sure the prior is zero in case our description of the PTA signal is inadequate.The effect of these extra criteria is that they alter the joint prior distribution of S 0 , η, and e 0 , and this is shown as a corner plot in Figure 1.
It is evident from Figure 1 that the upper limit of the valid prior space for S log 10 0 depends on e 0 and η, and it is lower for large e 0 and small η values.This is expected because for the same PTA signal amplitude S 0 (which is proportional to M ch 5 3 at the leading order; see Equation ( 14)), a system with higher eccentricity is more relativistic due to the lower relative separation at the periapsis.Further, for the same PTA signal amplitude S 0 , a binary with a very small value of η demands a very large value of the total mass (as evident from Equation (13)), which in turn results in a large value of the rate of advance of periapsis k.These properties of the valid joint prior distributions can affect our analysis, especially the calculation of the upper limit of S 0 in case of non-detection, in certain regions of the parameter space.This will be further discussed in Section 5. We stress that these criteria represent the limitations of our PTA signal model rather than being physical.We also note in passing that a similar constraint on the prior distribution based on the innermost stable circular orbit frequency was employed in the context of circular binary searches in Aggarwal et al. (2019) and Arzoumanian et al. (2023).

Model Comparison
In order to determine the significance of the detection of a PTA signal due to an eccentric binary, we compare the model  1 = timing model + WN + IRN + CURN + eccentric PTA signal against the model  0 = timing model + WN + IRN + CURN by computing the Bayes factor using the Savage-Dickey formula (Dickey 1971).We are able to apply the Savage-Dickey formula here because  1 and  0 are nested models, with  0 being the special case of  1 when S 0 = 0.The Bayes factor is given in this case by where  represents the data, and = represent the values of the marginalized prior distribution and the posterior distribution at S 0 = 0, respectively.The uncertainty in this estimate is given by  N 10 0 where N 0 is the number of samples in the lowest amplitude bin.

Search Pipeline Implementation
We use the GWecc.jl73package (Susobhanan 2023) to evaluate the PTA signal and the ENTERPRISE74 package (Ellis et al. 2019) to evaluate the single-parameter prior distributions and the likelihood function.In GWecc.jl, the eccentric_pta_signal_target function computes the PTA signal parametrized for a targeted search, the gwecc_-target_block function provides a corresponding deterministic signal object that can be used with ENTERPRISE, the PsrDistPrior class implements the informative pulsar distance priors, and validate_params_target function provides the criteria for defining the modified joint prior distribution as described in Section 3.3.The Hermite splinebased method described in Susobhanan (2023) is used to accelerate the PTA signal computation.
We use PTMCMCSampler75 (Ellis & van Haasteren 2017) to draw samples from the posterior distribution.PTMCMCSampler implements Adaptive Metropolis (AM), Single Component Adaptive Metropolis (SCAM), and differential evolution (DE) proposal distributions in combination with parallel tempering to effectively draw samples from a given posterior distribution.PTMCMCSampler also supports user-defined proposal distributions along with tunable weights for each proposal distribution.We use the following proposals available in the enterprise_extensions76 package (Taylor et al. 2021): draws from the single-parameter priors (JumpProposal.draw_from_prior and Jump-Proposal.draw_from_par_prior), and draws from empirical distributions (JumpProposal.draw_fro-m_empirical_distr) for IRN and CURN parameters.In addition, we use parallel tempering with four geometrically spaced temperatures.ENTERPRISE, enterprise_extensions, and PTMCMCSampler are described in Johnson et al. (2023).The search and visualization scripts used in this work can be found at https://github.com/lanky441/NG12p5_3C66B_GWecc.

Simulations
Before applying our methods to the NG12.5 data set, we perform targeted searches for eccentric PTA signals in simulated data sets to test our methods' performance and reliability.Inspired by Charisi et al. (2023), we carry out two separate analyses: an Earth term-only (E-only) analysis (pulsar term is excluded from the searched for signal model) and an Earth+pulsar term (E+P) analysis (pulsar term is included in the searched for signal model).
The simulated data sets used in our study correspond to the ten pulsars with the highest dropout factors (Aggarwal et al. 2019) as estimated in the NG12.5 GWB search (Arzoumanian et al. 2020a), and they are generated using the libstempo package77 (Vallisneri 2020), a python wrapper for tempo278 (Hobbs et al. 2006).We begin by creating an ideal noise-free realization of the pulsar TOAs by subtracting the timing residuals from the measured TOAs taken from the NG12.5 narrowband data set.We inject a WN realization with the same measurement uncertainties as in the NG12.5 narrowband data set, with EFAC = 1 and EQUAD = 0. We also inject a CURN realization with A CURN = 2.4 × 10 −15 (the GWB amplitude estimated from the NANOGrav 15 yr data set Agazie et al. 2023d), and γ CURN = 4.33 as expected of a GWB from a population of circular SMBHBs undergoing GW emissioninduced orbital evolution (Phinney 2001).We do not include IRN realizations in the simulated data for the sake of simplicity and to ensure that the analysis finishes within a reasonable time.Thereafter, we inject an eccentric PTA signal from a 3C 66B-like SMBHB, i,e., an SMBHB at the same sky location of 3C 66B with the same proposed binary parameters but a different distance.Both the Earth and the pulsar terms contributions are included while injecting the PTA signal.Finally, we fit the original timing model to the signal-injected TOAs in a maximum likelihood way and save the resulting post-fit timing model and TOAs as par and tim files.Note that we use the same simulated data set for both the E-only and E+P analyses since a physical PTA signal will include both the Earth term and the pulsar term contributions.
The aforementioned injected eccentric SMBHB signal is computed using the binary period and sky location as listed in Table 1, and the binary mass estimates given in Iguchi et al. (2010).In order to explore the performance and reliability of our pipeline for different signal strengths, we created simula- , where D L true is the luminosity distance value listed in Table 1.
We performed the E-only and E+P searches in the simulated data sets using our pipeline described in Section 3.5 and the priors listed in Table 1 (we use the uniform prior on S log 10 0 here).We fix the WN parameters to their injected values and search for a CURN and an eccentric PTA signal while analytically marginalizing the linearized timing model.
The posterior samples obtained from these searches, for the simulated data set with injected luminosity distance = D 5 L true , are shown as corner plots in Figures 2 (E-only) and 3 (E+P).In both the E-only and E+P cases, we find that the recovered CURN parameter values are consistent with the injected values within the 2σ level.We find that the estimated A CURN and γ CURN values in the E-only case show an offset as compared to the injected value whereas they are in closer agreement in the case of the E+P case.Whether this is a general feature of E-only searches caused by the unmodeled pulsar terms will be investigated in a future work.
The injected signal is detected in both the E-only and E+P searches; the S log 10 0 posterior distribution is consistent with the injected value although it has a long tail to the left.The l 0 , cos i, and η posterior distributions are broad but informative and consistent with the injected values.The e 0 posterior distributions for both the E-only and E+P searches are also informative and consistent with the injected values but are also consistent with e = 0.This is especially true in the case of the E +P search, where the posterior becomes flat as e → 0. Further, in both cases, ϖ 0 has a nearly flat posterior and ψ has a bimodal posterior where one of the peaks is consistent with the injected value.The bimodal structure in the posterior distribution of ψ is a consequence of the π-periodicity of the PTA signal expression in this parameter (see Equation (4)).Overall, we see that both the E-only and E+P searches are able to detect the injected signal but are not able to estimate many of the model parameters precisely.
In the case of the E+P search, we found that the posterior distributions for D p , ϖ p , and l p for all the pulsars are very similar to the prior distributions of the respective parameters, suggesting that the data does not provide any new insights or information about these parameters.Further, we find that the marginalized joint posterior distribution of e 0 and η has a geometry reminiscent of the prior distribution plotted in   4 that the upper limit on e 0 seen in Figure 3 is informed by the data, whereas the lower limit on η is an artifact of the prior.We have found that similar conclusions can be drawn for the posterior distributions for e 0 and η for the E-only search.
We repeated the above-described simulation studies with different injected values of D L , namely D L true and D 10 L true .Injecting the original value of the D L (D L true ) results in a nondetection where we are able to put an upper limit on S log 10 0 in both the E-only and E+P cases.It turned out that the posterior distribution in this case does not gain any significant information from the data for certain combinations of e 0 and η, where the S log 10 0 upper limit is not physically meaningful (this is also seen in the D 5 L true case, as can be seen in Figure 4).We defer a detailed discussion on this caveat to Section 5. Injecting a lower value of the D L (namely D 10 L true ) resulted in a detection in the E-only search, and the geometry of the posterior distribution turned out to be very similar to Figure 2. Unfortunately, the corresponding E+P search failed to adequately explore the parameter space.We suspect that this is due to a combination of its high dimensionality arising from the pulsar term parameters and complex geometry arising from the strong signal present in the data, i.e., our search method is inadequate in the strong signal regime in the E+P case.
The limited number of simulations we performed indicates that an E-only search can detect a PTA signal from an eccentric SMBHB, which is consistent with the results found in  Charisi et al. (2023) for circular binaries.However, a more rigorous simulation study is required to determine whether an E-only search will introduce a bias in the parameter estimations (or the upper limits in case of a non-detection) compared to an E+P search for eccentric binaries in PTA data sets.

Results
We now present the results of our search for continuous GWs from an eccentric SMBHB in 3C 66B in the NG12.5 data set.We run both an E-only search and an E+P search using the uniform priors for S log 10 0 for the detection analysis.The full signal model including a linearized timing model, WN, IRN, CURN, and eccentric PTA signal, as described in Equation (15), is used.The Bayes factors estimated using the Savage-Dickey formula (see Section 3.4) turned out to be 0.99 ± 0.01 for the E-only model and 1.01 ± 0.01 for the E+P model.It is evident that we did not detect any GWs from an eccentric SMBHB in 3C 66B as these numbers do not favor the presence of an eccentric PTA signal in our data.
In the absence of a detection, we now turn our attention to constraining the PTA signal amplitude S 0 and the chirp mass M ch for a possible eccentric SMBHB present in 3C 66B.The posterior samples for the upper limit analysis are obtained by reweighting the posterior samples from the detection analysis as described in Hourihane et al. (2023).Before calculating the upper limits, let us recall that, for certain parts of the parameter space, the posterior distribution of S 0 may only be informed by the prior distribution, determined in part by the validation criteria described in Section 3.3 and visualized in Figure 1 (this is seen to some extent in the injection study results shown in Figures 3 and 4).Therefore, the upper limit on S 0 (and M ch ) calculated for those parts of the parameter space can be misleading.Since e 0 , η, and S log 10 0 are not independent parameters in the prior distribution, it is more meaningful to explore the posterior distributions of S log 10 0 as functions of (e 0 , η) rather than marginalize over these parameters.
We divide the prior ranges for both e 0 and η into eight equally spaced bins and calculate the upper limit on S 0 for each e 0 and η bin from the posterior distribution while marginalized over all other parameters.The 95% upper limits on S log 10 0 for all pairs of bins from the E-only search are plotted using color maps in the upper panel of Figure 5. Recall that the upper limits are obtained after reweighting the samples generated by the detection analysis run.Further, we calculate the 95% upper limit on S log 10 0 from the modified joint prior distributions (Figure 1) for each (e 0 , η) pixel following the same way.The values of the 95% upper limits on S log 10 0 derived from the posterior distribution are quoted in each (e 0 , η) pixel along with the corresponding 95% upper limits obtained from the prior distribution within parentheses.We also plot the 95% upper limit on M log 10 ch for the E-only search in the lower panel of Figure 5 and the corresponding values for each pixel are also quoted.We see from panel (a) of Figure 5 that for pixels with high e 0 and/or low η values, the 95% upper limits obtained The underlying simulated data is the same as in Figure 2. We see that the upper limit of the posterior distribution for e 0 is not limited by the limitations in the prior.However, the lower limit of the posterior distribution of η may be limited by the prior.bottom) obtained from the Earth term-only search, binned in e 0 and η, plotted using color maps.S 0 and M ch are expressed in units of seconds and M ☉ , respectively.In each (e 0 , η) pixel of each panel, the 95% upper limit value is quoted.In the top panel, the 95% upper limit obtained from the prior distribution is quoted in parentheses.The pixels where the S log 10 0 upper limit may be restricted by the prior distribution and not determined by the data are highlighted with red font color and do not represent physically relevant results (they arise due to the limitations of our PTA signal model, see the text for discussion).We are able to obtain astrophysically meaningful constraints for e 0 < 0.5 and η > 0.1.from the posterior distribution of S log 10 0 are very close to those obtained from the prior.This is an indication that the upper limits of S log 10 0 for those pixels is determined by the prior, which we adopted to avoid regions of parameter space where the signal model is not valid, rather than the data.Hence, our upper limits for such pixels will not be astrophysically meaningful due to being affected by the limitations of our PTA signal model.We discern the pixels where the upper limit may be arising due to the restrictions imposed by our joint prior distribution and not by the data, using the following criterion: is less than a conservative tolerance of 5%.These upper limits do not represent physically meaningful results, arising due to the limitations of our PTA signal model, and we highlight such pixels using red font color in Figure 5.Note that, in principle, the upper limits can be limited by the prior distribution even in cases where the posterior distribution gains significant information over the prior from the data.This is why we base the above criterion on the upper limits rather than on a statistic such as KullbackLeibler divergence (Kullback & Leibler 1951) that treats the prior and posterior distributions more holistically.
Upper limits on S log 10 0 and M log 10 ch for the Earth+pulsar term search are calculated in the same way and shown in Figure 6.It can be seen from Figures 5 and 6 that, for the subset of the parameter space where e < 0.5 and η > 0.1, the upper limits are informed by the data and not just the restrictions of the prior.This is consistent with the shape of the joint prior distribution seen in Figure 1.We have also computed the Savage-Dickey Bayes factor in each pixel, and they remain close to 1 indicating non-detections in both the E-only and E+P searches.
We plot the marginalized posterior distributions for the CURN parameters and S log 10 0 , e 0 , and η obtained from the E-only and E+P detection analyses in Figures 7 and 8, respectively.We omit the samples with e 0 > 0.5 or η < 0.1 in these plots since they are significantly affected by the prior distribution as seen above.The 95% upper limits (calculated after reweighting) on S 0 are 88.1 ± 3.7 ns and 81.74 ± 0.86 ns, respectively, for the E-only and E+P analyses while restricting e 0 < 0.5 and η > 0.1.Similar 95% upper limits on M ch are (1.98 ± 0.05) × 10 9 M ☉ and (1.89 ± 0.01) × 10 9 M ☉ .The uncertainties in the upper limits are calculated using the bootstrap method.We also overplot the posterior distributions for the CURN parameters obtained from the NG12.5 GWB search (Arzoumanian et al. 2020a) in Figures 7 and 8.We observe that the addition of the eccentric SMBHB signal does not alter the posterior distribution of the CURN parameters appreciably in both the E-only and E+P searches.This indicates the robustness of our search against the leakage of power between the CURN and the eccentric SMBHB signal.Moreover, for the E+P search, we found that the posterior distributions for D p , ϖ p , and l p for all the pulsars closely align with their respective prior distributions.
We now turn our attention to the robustness of our search against the leakage of power between IRN and the eccentric SMBHB signal.We do this by comparing the posterior distributions of the IRN parameters obtained from our E-only and E+P detection analyses against those obtained from the NG12.5 GWB search using the Raveri-Doux tension statistic, available in the tensiometer package (Raveri & Doux 2021).Figure 9 shows a histogram of the tensiometer statistics computed for different pulsars, and we see that the IRN parameter posteriors from our analyses agree with those from the NG12.5 GWB search (Arzoumanian et al. 2020a) within the 0.04σ level for all pulsars.This gives us confidence that the search for the eccentric SMBHB signal is not affected significantly by the IRN.

Summary and Discussion
We have developed a pipeline for performing a targeted search for continuous GWs from individual eccentric SMBHB in a PTA data set.We calculate the pulsar timing residuals induced by GWs from an eccentric SMBHB using the semianalytic approach presented in Susobhanan et al. (2020) and Susobhanan (2023).We tested our pipeline by applying it to simulated data sets, and we have performed an E-only search where contributions from the pulsar term are neglected and an E+P search where contributions from both terms are included.
In both cases, we have found that the marginalized posterior distributions for our free model parameters are consistent with the injected values for the simulated data.Thereafter, we performed both E-only and E+P targeted searches for GWs from an eccentric SMBHB in the radio galaxy 3C 66B in the NANOGrav 12.5 yr data set.In addition to the eccentric SMBHB signal, our model incorporates a CURN process, detected as a precursor to the GWB in the NG12.5 data set (Arzoumanian et al. 2020a), IRN processes for each pulsar, and WN processes for each pulsar.We fix the WN parameters for each pulsar to the values obtained from its single-pulsar noise analysis.This is the first time a multi-messenger targeted search for GWs from an eccentric SMBHB is performed on a fullscale PTA data set.We did not find any evidence for continuous GWs from an eccentric SMBHB in 3C 66B in our search.Hence, we calculated the upper limits on the GW signal amplitude (S 0 ) and the chirp mass (M ch ) of the SMBHB candidate in 3C 66B as a function of the initial eccentricity (e 0 ) and symmetric mass ratio (η) of the binary.We found that for certain combinations for e 0 and η of the binary, the S log 10 0 upper limit is limited by the prior and is not informed by the data.This arises due to the limitations of our PTA signal model and such upper limits should not be considered astrophysically meaningful.For the parameter range of e 0 < 0.5 and η > 0.1, the upper limits are not influenced by the limitations of our PTA signal model and are informed by the data.In that regime, the 95% upper limits obtained from the E-only search are 88.1 ± 3.7 ns for S 0 and (1.98 ± 0.05) × 10 9 M ☉ for M ch .Similar 95% upper limits obtained from the E+P search are 81.74 ± 0.86 ns for S 0 and (1.89 ± 0.01) × 10 9 M ☉ for M ch .We note that the peculiar motion of 3C 66B can give rise to an uncertainty in the luminosity distance estimate from the redshift of that galaxy.This may lead to a systematic bias in the estimated upper limits on M ch whereas the upper limits on S 0 should be unaffected.We also found that the posterior distributions of the CURN and IRN parameters in both our searches are consistent with those obtained from the NG12.5 GWB search (Arzoumanian et al. 2020a).
The is consistent with the upper limits we obtain from our analysis.We also see that the upper limits on the chirp mass from our targeted search for GWs from eccentric binary in 3C 66B are higher than that obtained from the targeted search from circular binaries in the NG12.5 data set (Arzoumanian et al. 2023).This may be because of the higher number of free parameters present in our eccentric search as compared to the circular search.Another reason for this may be that, as the eccentricity pushes the signal power to higher frequencies, the dominant frequency might lie at high frequencies where the data is less sensitive due to WN.Further, the power of an eccentric signal is spread across multiple frequencies and the signal at some of these frequencies may lie below the noise threshold.
We do not find any clear trend in the upper limits on S log 10 0 as a function of e 0 or η for both the E-only and E+P searches, as evident from the 2D posterior distributions shown in Figures 7 and 8.This is not surprising because of the parameterization used for our targeted search (discussed in Section 2.3), where the overall PTA signal amplitude at the Newtonian order is just given by S 0 .Although different values of e 0 and η can lead to different temporal evolution of the PTA signal amplitude S(t) through ς(t) (due to different evolution of n(t) and k(t) appearing at higher PN order), we suspect that for a non-detection these do not affect the upper limits significantly.This is also evident from Figure 10, where we show the crest-to-trough amplitudes of the gravitational waveform (top panel) and the PTA signal (not to be confused with S 0 ) induced by GWs from an SMBHB (bottom panel) as a function of the eccentricity when we keep all the other binary parameters fixed.We see that although the gravitational waveform amplitude increases with the eccentricity, the amplitude of the PTA signal does not vary significantly with the eccentricity.Therefore, we expect that the upper limit on S 0 should also not depend significantly on the eccentricity of the binary.Similarly, upper limits on M ch also should not depend significantly on e 0 and η as we can see from Equation (14) that at the Newtonian order, the PTA signal amplitude depends only on M ch and not on e 0 or η.Please note that Figure 10 does not imply anything about the measurability of eccentricity in the case of a detection, which would only depend on the quality and sensitivity of the data set as well as the strength of the signal.
In this work, we have used a Gaussian red noise process with a 30-component power-law spectrum for modeling IRN and a DMX model to account for DM variations, following standard NANOGrav practices (Arzoumanian et al. 2020a(Arzoumanian et al. , 2023;;Agazie et al. 2023d).However, different PTAs make different choices for modeling IRN and DM variations in their data.Recently, Agazie et al. (2023a) compared the GWB search results from different PTAs and showed that, despite making different modeling choices, there is no significant difference in the GWB parameters measured by different PTAs.Further, for the majority of the pulsars, the IRN parameters are consistent among different PTAs.This suggests that the choice of individual pulsar noise models should not significantly affect the upper limits obtained in our analysis.
We find that certain aspects of our methods and the search pipeline can be improved further, especially in order to perform targeted PTA searches for GWs from eccentric binaries based on a catalog of SMBHB candidates.For many such candidates, the orbital period of the binary is not very accurately known, and therefore using a Gaussian prior for the initial binary orbital period (P 0 ) around the proposed values instead of fixing it should lead to more realistic results.Further, due to the limitation of our PTA signal model, we were able to calculate meaningful upper limits only for a certain regime of the parameter space (of e 0 and η) and that regime depends on the reference time t 0 and the total span of our data.A later reference time t 0 for the binary parameters will extend that regime whereas a longer data span will shrink it.Therefore, mitigating the limitations of our PTA signal model will be crucial for exploring the whole (e 0 , η) parameter space and Figure 10.Maximum−minimum (crest-to-trough) amplitude of the gravitational waveform (top panel) and PTA signal (bottom panel) for PSR J1909-3744 due to an SMBHB as a function of orbital eccentricity.The crest-totrough amplitude is the difference between the maximum and minimum values of the waveform for a given time span and is not to be confused with S 0 .We have only considered the Earth term.The masses of the binary are fixed to the 3C 66B binary model values given in Iguchi et al. (2010): m 1 = 1.2 × 10 9 M ☉ , m 2 = 7.0 × 10 8 M ☉ .We have used cos 1 i = , ψ = 0, l 0 = 0, and γ 0 = 0, and a time span of 30 yr.The sky position and luminosity distance of the binary are taken from Table 1.We do not find any strong variation in the PTA signal amplitude with eccentricity, whereas the waveform amplitude is an increasing function of eccentricity.Note that the overall weak decreasing trend seen in the bottom plot is not universal and changes depending on binary parameter values.
obtaining meaningful upper limits on the PTA signal amplitude.This will be especially relevant for CHIME, FAST, and SKA-era PTA data sets, where a high observation cadence will make the data sensitive to higher GW frequencies, which correspond to highly relativistic SMBHBs.Ways to extend the validity of our PTA signal model closer to the merger event include the effective one-body formalism (Hinderer & Babak 2017) and replacing the currently unmodeled merger and ringdown parts of the waveform with a generic burst signal model.
A new method to search for individual sources in PTA data sets, named QuickCW (Bécsy et al. 2022), was introduced recently and was applied to great effect to search for circular SMBHBs in the NANOGrav 15 yr data set (Agazie et al. 2023e).QuickCW exploits the mathematical structure of the PTA signal expression to accelerate the likelihood computation in the case of projection parameter updates (e.g., in our model, the projection parameters are S 0 , ψ, cos i, ϖ 0 , and ϖ p for each pulsar).Extending this method for eccentric SMBHB searches will be a promising avenue to explore in order to keep future searches computationally feasible in the face of growing data volumes.Further, Charisi et al. (2023) showed that ignoring the pulsar term contributions does not significantly affect the parameter estimation in the case of targeted PTA searches for circular SMBHBs.Performing detailed simulation studies to confirm that the same is true for targeted searches of eccentric SMBHBs for both detection and upper limit analyses will help us efficiently perform such searches in future PTA data sets by ignoring the pulsar term.It will also be interesting to explore the possibility of using Hamiltonian Monte Carlo (HMC) to search for SMBHB signals in PTA data sets, given its performance advantages over other types of Markov Chain Monte Carlo (MCMC) methods in high-dimensional parameter spaces (see Freedman et al. 2023 for a PTA application of HMC).These considerations will be especially relevant for conducting an eccentric SMBHB search in the upcoming IPTA Data Release 3, which is expected to be much larger and more sensitive than the NG12.5 data set.
support of the Department of Atomic Energy, Government of India, under project identification No. RTI 4002.We thank the anonymous referee for their comments and suggestions.
tions with different GW source luminosity distance values, namely D L true ,

Figure 1 .
Figure1.Corner plot of the joint prior distribution of S 0 , η, and e 0 after the validation criteria ( ) V S e log , , 10 0 0 h have been applied.The samples are drawn from the prior using rejection sampling.

Figure 2 .
Figure 2. Posterior distributions for the CURN parameters (γ CURN and A log 10 CURN ) and the eccentric SMBHB parameters, marginalized over other parameters, for our Earth term-only search for GWs from eccentric binary in a simulated data set.The sky location and period of the simulated GW source are taken from Table1, the binary masses are taken fromIguchi et al. (2010), and the luminosity distance is taken to be one-fifth the value given in Table1.The green lines indicate the true values of the parameters for the injected signals used to create the simulated data set.

Figure 1 .
Figure 1.To see if the (e 0 , η) joint posterior provides more information than the corresponding prior distribution, we overplot them in Figure4.To provide a fair comparison, we restrict the prior samples in Figure4to those where the S log 10 0 values fall within the 95% credible interval of S log 10 0 derived from the posterior distribution.It is clear from Figure4that the upper limit on e 0 seen in Figure3is informed by the data, whereas the lower limit on η is an artifact of the prior.We have found that similar conclusions can be drawn for the posterior distributions for e 0 and η for the E-only search.We repeated the above-described simulation studies with different injected values of D L , namely D L true and D 10

Figure 3 .
Figure 3. Posterior distributions of the parameters for our Earth+pulsar term search in simulated data.See the caption of Figure 2 for details on the underlying simulated data set and plotting conventions.

Figure 4 .
Figure 4. Comparison between the valid joint prior distribution and the 2D marginalized posterior distribution of e 0 and η for the E+P search of GWs in simulated data.The underlying simulated data is the same as in Figure2.We see that the upper limit of the posterior distribution for e 0 is not limited by the limitations in the prior.However, the lower limit of the posterior distribution of η may be limited by the prior.

Figure 5 .
Figure 5. 95% upper limits on S log 10 0 (panel (a), top) and M log 10 ch (panel (b), bottom) obtained from the Earth term-only search, binned in e 0 and η, plotted using color maps.S 0 and M ch are expressed in units of seconds and M ☉ , respectively.In each (e 0 , η) pixel of each panel, the 95% upper limit value is quoted.In the top panel, the 95% upper limit obtained from the prior distribution is quoted in parentheses.The pixels where the S log 10 0 upper limit may be restricted by the prior distribution and not determined by the data are highlighted with red font color and do not represent physically relevant results (they arise due to the limitations of our PTA signal model, see the text for discussion).We are able to obtain astrophysically meaningful constraints for e 0 < 0.5 and η > 0.1.
upper limits obtained from the prior and posterior distributions, respectively.This implies that we treat the upper limit for a pixel as invalid if the fractional difference between [

Figure 6 .
Figure 6.95% upper limits on S log 10 0 (panel (a), top) and M log 10 ch (panel (b), bottom) obtained from the Earth+pulsar term search, binned in e 0 and η, plotted using color maps.The plotting conventions are identical to those in Figure5.We are able to obtain physically relevant constraints for e 0 < 0.5 and η > 0.1.For comparison, the 95% upper limit from the Earth+pulsar term targeted search of a circular binary in 3C 66B using NG12.5 data set is M log 9.15 10 ch <(Arzoumanian et al. 2023).

Figure 7 .
Figure 7. Posterior distribution for the CURN parameters (γ CURN and A log 10 CURN ) and the eccentric SMBHB signal parameters (e 0 , η, and S log 10 0 ) for the Earth term- only search marginalized over all other parameters (plotted in blue).We only include samples with e < 0.5 and η > 0.1 in this plot as the posterior distribution is significantly affected by the prior outside this region of the parameter space.The CURN parameters obtained from the NANOGrav 12.5 yr GWB search are overplotted in green for comparison.
SMBHB model of 3C 66B proposed by Iguchi et al. (2010) with a chirp mass of

Figure 8 .
Figure 8. Posterior distribution for the CURN parameters (γ CURN andA log 10 CURN ) and the eccentric SMBHB signal parameters (e 0 , η, and S log 10 0 ) for the Earth +pulsar term search marginalized over all other parameters (plotted in orange).We only include samples with e < 0.5 and η > 0.1 in this plot as the posterior distribution is significantly affected by the prior outside this region of the parameter space.The CURN parameters obtained from the NANOGrav 12.5 yr GWB search are overplotted in green for comparison.

Figure 9 .
Figure 9.The Raveri-Doux tension between marginalized IRN posterior distributions obtained from the NANOGrav 12.5 yr GWB search and this work for different pulsars, plotted as a histogram.The blue histogram corresponds to the Earth term-only search and the orange histogram corresponds to the Earth +pulsar term search.The tension remains less than 0.04σ for all pulsars in both searches.This implies that our searches are not significantly affected by the IRN. )

Table 1
Sudou et al. (2003) for the PTA Signal Parameters Note.The source coordinates of our GW source 3C 66B are taken from NASA/IPAC Extragalactic Database (http://ned.ipac.caltech.edu/).P 0 is taken fromSudou et al. (2003)and D L is estimated from the redshift assuming a flat Universe with cosmological parameters given inAghanim et al.