A Fourier Domain "Jerk" Search for Binary Pulsars

Binary pulsar systems produce extraordinary physics and astrophysics, yet detecting them can be a particularly difficult challenge. Isolated pulsars are detected by taking the fast Fourier transform (FFT) of an observational time series and searching through the resulting power spectrum. Because pulsars generally have narrow pulse profiles, resulting in many significant harmonics of the spin frequency in the power spectrum, those harmonics are summed to increase the significance of the final detection (e.g., Lorimer & Kramer 2012).

However, in the case of binary pulsars, detection is complicated by the orbital motion of the pulsar. If the observation duration T is longer than even a small fraction of the orbital period P_orb, Doppler shifting causes the apparent pulsar spin frequency to change with time. The Fourier power of the signal is "smeared" across neighboring frequency bins, leading to a decrease in the significance of the detection. Higher harmonics of the pulsar signal are impacted more than the fundamental, as Doppler smearing increases proportionally with harmonic number.

Acceleration searches account for this Doppler smearing by assuming that, over a small enough fraction of the orbit (T ≲ P_orb/10), the pulsar's acceleration is approximately constant, and therefore its measured spin frequency drifts linearly in time (e.g., Johnston & Kulkarni 1991). Blind acceleration searches iterate over many constant acceleration values, completing data manipulations in the time or Fourier domain to recover signal power that the orbital motion has smeared over several frequency bins.

The higher search dimensionality and computational expense of acceleration searches has precluded their use in wide pulsar surveys until this past decade. However, for high-priority targeted searches, such as those toward supernova remnants (e.g., Middleditch & Kristian 1984), globular clusters (e.g., Camilo et al. 2000), and Fermi unassociated sources (e.g., Ray et al. 2012), acceleration searches have uncovered more than 100 binary pulsars. The vast majority were found using a Fourier domain implementation of the acceleration search (Ransom et al. 2001, 2002), which has recently become known as the Fourier Domain Acceleration Search (FDAS; Dimoudi et al. 2018).³

While the FDAS is effective when T ≲ P_orb/10, search sensitivity to fainter pulsars improves with longer observing durations as $\sqrt{T}$. Acceleration searches are therefore limited to discovering only the brightest binary pulsars. To enable the detection of fainter and more accelerated systems, we can add another level of approximation to the acceleration search by assuming that the next derivative of orbital motion, the "jerk," is constant. Under this assumption, the measured spin frequency changes quadratically with time and the spin phase changes cubically.

Jerk searches allow for longer observation durations covering larger portions of the pulsar orbit than acceleration searches, while still retaining significant power in the Fourier domain (Bagchi et al. 2013). Therefore, although adding the jerk dimension substantially increases the runtime of the search, it also allows us to probe for systems previously missed with acceleration searches due to residual Doppler smearing. This capability will be especially useful in targeted searches for more exotic systems like Galactic Center pulsars, pulsar-black hole binaries, and the most compact and relativistic double neutron star systems.

We have developed and tested a Fourier domain jerk search as part of the PRESTO⁴ software package (Ransom 2001). This Letter briefly describes the details and performance of our implementation, and presents the detection of a new pulsar in the globular cluster Terzan 5 found using our jerk search.

Our jerk search functions as a fairly straightforward extension to the original FDAS implemented in PRESTO's accelsearch program. The mathematics and methodology of this FDAS are described in detail in Ransom et al. (2002) and Dimoudi et al. (2018). In the following section, we provide a brief summary.

2.1. Acceleration Search Review

In short, the FDAS is a matched filtering algorithm that corrects for Doppler smearing by assuming that the pulsar's acceleration, α, is roughly constant over the course of the observation. Under this assumption, each harmonic of the pulsar signal would experience a constant frequency derivative, $\dot{f}$, according to the relation

$$\begin{eqnarray}&&\alpha =\displaystyle \frac{\dot{{fc}}}{{hf}}=\displaystyle \frac{{zc}}{{{hfT}}^{2}},\end{eqnarray} \tag{ 1 }$$

where c is the speed of light, T is the observation duration, f is the frequency of the fundamental, and h is the harmonic number where h = 1 represents the fundamental. In this equation we also introduce the Fourier frequency derivative, z. Fourier frequencies (i.e., wavenumbers or FFT bin numbers) are defined as r = fT, and so $z=\dot{r}=\dot{f}{T}^{2}$ corresponds to the number of Fourier frequency bins that the signal drifts through over the course of the observation. Similarly, $w=\dot{z}=\ddot{f}{T}^{3}$ corresponds to the Fourier jerk, or the number of frequency derivative bins that the signal drifts through over the course of the observation. We use this r, z, w coordinate system throughout our acceleration and jerk search code, as it is computationally and intuitively advantageous when dealing with the properties of discrete Fourier transforms. When using these Fourier-based units defined over the whole observation, it is also convenient to switch to a normalized time coordinate, u, which represents the fraction of the observation complete at any given instant (such that 0 ≤ u ≤ 1).

Using the Convolution Theorem, the acceleration search correlates template Fourier domain amplitude and phase responses for a number of trial accelerations with the complex Fourier amplitudes from an FFT of the original input time series. The correlations are completed as a series of short FFTs using a computationally efficient "overlap-and-save" technique (see Section 3 of Dimoudi et al. 2018).

When stacked according to the trial $\dot{f}$ value, the resulting power spectra from these correlations form a 2D plane in Fourier frequency versus frequency derivative space. Figure 1 shows examples of such $f\mbox{--}\dot{f}$ planes for both simulated and actual pulsar signals. Once the $f\mbox{--}\dot{f}$ plane is constructed, the algorithm searches through it and identifies candidate pulsar signatures according to a pre-calculated power threshold. The search also incoherently sums the powers from a number of harmonics onto the fundamental to increase the probability of detecting narrow pulse profiles.

Zoom In Zoom Out Reset image size

For an acceleration search, the template responses can be calculated analytically (see Section 4.2.2 of Ransom et al. 2002, for a complete mathematical derivation). With a constant acceleration, we can approximate a harmonic of the pulsar signal as a sinusoid with a quadratically varying phase. Taking the Fourier transform of this signal yields an analytic expression for the template, composed of Fresnel integrals. The simulated signal in the bottom-left panels of Figure 1, which show the response of an accelerated (but jerk-corrected) pulsar across three harmonics, also effectively demonstrates the distinctive "X" shape that these templates form when displayed in the $f\mbox{--}\dot{f}$ plane. In the bottom-right panels, the z = 0 lines show the powers of the initial FFT of the time series, without any sort of acceleration or jerk correction.

2.2. Jerk Search

Just as the acceleration search is rooted in the approximation of constant acceleration via a linearly varying spin frequency, our jerk search approximates a constant jerk with linearly varying acceleration, or a quadratically varying spin frequency. A constant jerk corresponds to a constant second time derivative of the frequency,

$$\begin{eqnarray}&&\dot{\alpha }=\displaystyle \frac{\ddot{f}c}{{hf}}=\displaystyle \frac{{wc}}{{{hfT}}^{3}},\end{eqnarray} \tag{ 2 }$$

where $\dot{\alpha }$ is the jerk, $\ddot{f}$ is the second time derivative of the frequency, and w is the Fourier jerk.

The actual mechanics of the jerk search are very similar to what we have just described for the acceleration search. We use the same overall process of overlap-and-save matched filtering, threshold searching, and harmonic summing, except instead of generating an $f\mbox{--}\dot{f}$ plane of powers, we arrange the correlation results in an $f\mbox{--}\dot{f}\mbox{--}\ddot{f}$ (or r–z–w) volume, which we then search for candidates.

Generating this r–z–w volume requires us to calculate different template responses for correlation. Under the assumption of constant jerk, we approximate a harmonic of the pulsar signal n(u) as a sinusoid with a cubically varying phase

$$\begin{eqnarray}&&n(u)=a\cos \left[2\pi \left({r}_{0}u+\displaystyle \frac{1}{2}{z}_{0}{u}^{2}+\displaystyle \frac{1}{6}{{wu}}^{3}\right)+\phi \right],\end{eqnarray} \tag{ 3 }$$

where a is the amplitude of the signal, r₀ and z₀ are the Fourier frequency and frequency derivative at the start of the observation, and ϕ is a starting phase. To make our jerk search easier to implement, we want our collection of templates to be as compact in r, z, w space as possible around the template origin. To accomplish this, we can define our initial Fourier frequency and frequency derivative values, r₀ and z₀, in terms of the average frequency and average frequency derivative over the course of the observation, $\bar{r}$ and $\bar{z}$,

$$\begin{eqnarray*}&&{r}_{0}=\bar{r}-\displaystyle \frac{\bar{z}}{2}+\displaystyle \frac{w}{12}\ {\rm{and}}\ {z}_{0}=\bar{z}-\displaystyle \frac{w}{2}.\end{eqnarray*}$$

These relations come from averaging $r(u)=d{\rm{\Phi }}/{du}$ and $z(u)={d}^{2}{\rm{\Phi }}/{{du}}^{2}$ from u = 0 to u = 1, where Φ is the phase of the cosine in Equation (3).

Expanding Equation (3) into complex exponentials gives

$$\begin{eqnarray}n(u) & = & \displaystyle \frac{a}{2}\{{e}^{2\pi i\left[\left(\bar{r}-\tfrac{\bar{z}}{2}+\displaystyle \frac{w}{12}\right)u+\displaystyle \frac{1}{2}\left(\bar{z}-\displaystyle \frac{w}{2}\right){u}^{2}+\displaystyle \frac{1}{6}{{wu}}^{3}\right]}{e}^{i\phi }\\ & & +\ {e}^{-2\pi i\left[\left(\bar{r}-\tfrac{\bar{z}}{2}+\displaystyle \frac{w}{12}\right)u+\displaystyle \frac{1}{2}\left(\bar{z}-\displaystyle \frac{w}{2}\right){u}^{2}+\displaystyle \frac{1}{6}{{wu}}^{3}\right]}{e}^{-i\phi }\}.\end{eqnarray} \tag{ 4 }$$

The Fourier response of the jerked signal, at frequency r, is the Fourier transform of n(u),

$$\begin{eqnarray}{A}_{r} & = & N{\displaystyle \int }_{0}^{1}n(u){e}^{-2\pi {iru}}{du}\\ & = & \displaystyle \frac{{aN}}{2}{e}^{i\phi }{\displaystyle \int }_{0}^{1}{e}^{2\pi i\left[\left({\rm{\Delta }}r-\tfrac{\bar{z}}{2}+\displaystyle \frac{w}{12}\right)u+\displaystyle \frac{1}{2}\left(\bar{z}-\displaystyle \frac{w}{2}\right){u}^{2}+\displaystyle \frac{1}{6}{{wu}}^{3}\right]}\ {du},\end{eqnarray} \tag{ 5 }$$

where N is the number of samples in our original time series, ${\rm{\Delta }}r=\bar{r}-r$, and the second term in Equation (4) averages to zero. The correlation templates for the search are the complex conjugates of Equation (5) over a range of Δr.

Unlike the template integral for an acceleration search, Equation (5) has no analytic solution. However, as in the acceleration case, the template shape is independent of the absolute value of r₀ or $\bar{r}$. This independence allows us to simulate the portion of the binary orbit described by Equation (3) efficiently in a time series of only a few hundred or a few thousand points, depending on $\bar{z}$ and w. Taking an FFT of the time series then computes the Fourier integral numerically. The simulated signals in the bottom-right panel of Figure 1 show the uncorrected response of the first three harmonics of a jerked pulsar when displayed in the $f\mbox{--}\dot{f}$ plane (located at w = 0 in the $f\mbox{--}\dot{f}\mbox{--}\ddot{f}$ volume).

Once portions of the $f\mbox{--}\dot{f}\mbox{--}\ddot{f}$ volume are computed, they are typically converted to normalized powers, and then either incoherently added to other portions of the $f\mbox{--}\dot{f}\mbox{--}\ddot{f}$ volume at signal harmonics, or searched directly for significant signals (see e.g., Ransom et al. 2002). As per a normal Fourier search, the power-normalized incoherent harmonic sums are ${\chi }^{2}$-distributed with $2m$ degrees of freedom, where m is the number of summed harmonics. Determining the overall significance of a signal is beyond the scope of this Letter, however, for well-behaved data (i.e., a relatively uniformly sampled and white-noise-dominated input time series), the approximate number of independent Fourier bins searched must be accounted for (i.e., a search trials factor). For a jerk search, the approximate number of independent trials for a single choice of the number of harmonics to sum m is

$$\begin{eqnarray}&&{N}_{\mathrm{trials}}\sim \displaystyle \frac{{N}_{r}}{m}\left(\displaystyle \frac{{N}_{z}}{6.95}\right)\left(\displaystyle \frac{{N}_{w}}{44.2}\right),\end{eqnarray} \tag{ 6 }$$

where N_r, N_z, and N_w are the numbers of Fourier frequency bins, z bins, and w bins searched for the highest harmonic summed. The numerical values 6.95 and 44.2 are the Fourier widths at half-power of the Fourier response in the z and w directions, respectively.

When conducting searches, in order to help prevent "scalloping" of a signal's power due to the finite grid of computed r, z, and w points, accelsearch oversamples the volume in each of those directions. The grid spacing used is 0.5 (i.e., interbinning or Fourier interpolation), 2, and 20, for r, z, and w, respectively, given the measured half-widths of the peaks in each of those orthogonal directions (see above). To determine what range of z and w values to search, we simulated thousands of realistic pulsar binaries and "observed" them with integrations of tens of minutes to hours. Most systems are detected with $| z| \lt 200$ and $| w| \lt 600$. For larger values of z or w, given realistic orbits with stellar mass companions, regardless of T/P_orb, the constant acceleration or constant jerk assumption breaks down and we lose sensitivity.

The slices through the jerk volume displayed in Figure 1 illuminate other salient properties of the jerk search that are worth noting. First off, we can see that when the jerk volume is sliced through at the w of the pulsar signal (in this case 45.9), the response of the signal in the resulting $f\mbox{--}\dot{f}$ plane is the characteristic "X" shape that we expect from a standard acceleration search (see Ransom et al. 2002), indicating that the smearing due to jerk has been mitigated. For an animated projection through the r–z–w volume for a sinusoid, see Figure 2. Another important feature to note is that the higher harmonics of the signal are essentially "jerked" out of significance in the acceleration search. This is just as predicted by Equation (2), which shows that as the harmonic number (and therefore frequency, f) increases and the jerk remains constant, the Fourier w of the signal must increase proportionally. Thus, higher harmonics are more heavily affected by jerk. This effect is clearly illustrated in the left panels of Figure 1. As the frequency doubles from the fundamental to the second harmonic, the Fourier w also doubles from 45.9 to 91.7. As a result, the detection significance is much higher using the jerk search, as each of the harmonics has more power to contribute during harmonic summing.

Download figure:

Video Standard image High-resolution image

While developing and testing the algorithm and its implementation in accelsearch, we repeatedly searched two archival GBT observations of the globular cluster Terzan 5 taken 2005 May 15 and 2008 September 12. These data, described in Section 4 and in Ransom et al. (2005), were dedispersed at dispersion measures (DMs) of 238.00 and 238.72 pc cm⁻³, respectively. Each observation contains numerous binary millisecond pulsars (MSPs) with relatively short orbital periods (P_b ≲ 1 day) that are detectable with significant accelerations and jerks over integrations between 10 minutes to the ∼7 hr durations of the observations.

As an example of how detection significances can vary for real pulsars, we describe how five different Terzan 5 MSPs were detected in a single search of a 4096-s segment of the 2005 data, using accelsearch -numharm 4 -zmax 300 and with and without -wmax 900. Each of the reported detection significances are in σ (i.e., equivalent Gaussian significance) after correcting for the approximate number of independent trials searched (see Section 2.2), which is ∼41 times larger for the jerk searches than the acceleration searches, using these parameters.

For the isolated MSP Ter5L, the acceleration-only search detected the pulsar with a significance of 7.4σ, compared to 6.6σ for the jerk search. Similarly, the weakly accelerating long-period binary Ter5E was detected at 8.9 and 8.2σ, respectively, although the total summed power was larger in the jerk search. These results show, as expected, that you pay a penalty with a jerk search for weakly or unaccelerated pulsars due to the larger phase space searched. The situation was different for the compact binaries Ter5I, M, and N. Ter5I and Ter5M were detected with 1.2 and 3.0 extra sigma in the jerk searches, although Ter5N lost 0.9σ since w ∼ 0 during that portion of the pulsar's orbit.

These test searches, as well as the thorough analysis performed by Bagchi et al. (2013), show that jerk searches can improve the sensitivity to highly accelerating pulsars with T ∼ 0.05–0.15P_orb by a significant amount. The penalty is a slightly reduced sensitivity to weakly or unaccelerated pulsars due to the additional independent trials searched, as well as a substantially longer runtime that, in this case, increased by a factor of almost 80. In general, the runtime is roughly linear with the number of trials according to Equation (6).

While searching 11 overlapping segments of duration 4096 s from the 2008 observation, with -zmax 100 and -wmax 500, we detected a highly accelerated new pulsar with a spin period of 2.93 ms, z = −21, and w = −15 in one segment, using a four-harmonic sum.⁵ A similar search, using acceleration but no jerk, did not detect the pulsar, mostly because the higher harmonics had been "jerked" out of significance.

After determining a better DM for the pulsar, we searched other archival observations and detected it several additional times, allowing us to solve the compact and eccentric orbit, and eventually determine a fully coherent timing solution. Here we present the timing from the first ∼1000 days of Terzan 5 observations made with the GBT Pulsar Spigot (Kaplan et al. 2005). Most of the observations were taken at 2 GHz with a usable bandwidth of ∼600 MHz (out of 800 MHz total), a sample time of 81.92 μs, and 2048 frequency channels (see Ransom et al. 2005; Hessels et al. 2006, for details).

The timing solution for PSR J1748−2446am includes a strong detection of the advance of periastron (i.e., $\dot{\omega }$). If entirely relativistic, which is likely given a white dwarf companion, the total mass of the system is 1.85 ± 0.02 M_⊙, and the component masses can be constrained with the mass function by assuming random inclinations (see Table 1), in which case, the pulsar mass is ${1.649}_{-0.11}^{+0.037}$ M_⊙. A longer-duration timing solution, comprising ∼14 years of Spigot and coherently dedispersed Green Bank Ultimate Pulsar Processing Instrument (GUPPI) data will be presented elsewhere (S. M. Ransom et al. 2018, in preparation).

Table 1.  PSR J1748−2446am

Timing Parameters
Right Ascension (R.A., J2000)	17h48m048235(2)
Declination (decl., J2000)	−24°46'4721(9)
Pulsar Period (ms)	2.933819877244(2)
Pulsar Frequency (Hz)	340.8525546358(2)
Frequency Derivative (Hz s−1)	1.5893(4) × 10−14
Reference Epoch (MJD)	53700
Dispersion Measure (pc cm−3)	238.193(3)
Orbital Period (days)	0.80010926(2)
Projected Semimajor Axis (lt-s)	0.937815(5)
Orbital Eccentricity	0.204736(9)
Epoch of Periastron (MJD)	53700.440278(6)
Longitude of Periastron, ω (deg)	337.365(2)
Time Derivative of ω, $\dot{\omega }$ (deg yr−1)	0.454(4)
Span of Timing Data (MJD)	53193–54195
Number of TOAs	217
rms TOA Residual (μs)	38.5
Derived Parameters
Mass Function (M⊙)	0.00138336(2)
Total System Mass (M⊙)	1.85(2)
Min Companion Mass (M⊙)	≥0.15
Companion Mass (M⊙)	0.194 (+0.11, −0.023)
Pulsar Mass (M⊙)	1.649 (+0.037, −0.11)
Flux Density at 2 GHz (mJy)	∼0.015

Note. Numbers in parentheses represent 1-σ uncertainties in the last digit. The timing solution was determined using tempo with the DE436 Solar System Ephemeris and the DD binary model. The time system used is Barycentric Dynamical Time (TDB). The minimum companion mass was calculated assuming a pulsar mass of 1.4 M_⊙. The total system mass and 68% central confidence ranges on the masses of the pulsar and its companion were determined assuming that $\dot{\omega }$ is due completely to general relativity, and a random orbital inclination (i.e., probability density is uniform in $\cos \,i$).

Download table as:  ASCIITypeset image

In this Letter we have presented the implementation and performance of a new Fourier domain jerk search, which extends the PRESTO FDAS to account for linearly varying acceleration (i.e., constant jerk). Consistent with the analysis conducted by Bagchi et al. (2013), we find that our algorithm can significantly improve search sensitivity to pulsar binaries where T ∼ 0.05–0.15P_orb. Trade-offs of the algorithm include a significantly longer runtime and decreased sensitivity to unaccelerated systems due to the increased independent trials factor. While testing our algorithm on GBT observations of the globular cluster Terzan 5, we discovered Ter5am, an interesting compact eccentric binary that already shows relativistic periastron advance. Other relativistic effects (such as relativistic γ) will likely become detectable in coming years, allowing us to obtain precise mass measurements for both objects in the system.

Looking to the future, this jerk search technique will be especially useful in high-profile targeted searches for exotic systems such as Galactic Center and globular cluster pulsars, pulsar-black hole binaries, and the most compact and relativistic double neutron star systems. Because each of the template correlations in the r–z–w volume are independent of each other, this technique also has parallelization potential. While the current implementation makes use of OpenMP, a straightforward extension of the Dimoudi et al. (2018) methodology would allow us to implement the jerk search on graphics processor units. This could significantly reduce the runtime cost of the algorithm, opening it up to wider use and less focused searches.

The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. The Green Bank Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. S.M.R. is a CIFAR Senior Fellow and is also supported by the NSF Physics Frontiers Center award number 1430284. We thank Jason Hessels, Ingrid Stairs, and Paulo Freire for their help with GBT observations of Terzan 5 over the years, and Thankful Cromartie for a careful reading of the manuscript.

Facility: GBT. -

Software: PRESTO (Ransom 2001, 2011), Tempo (Nice et al. 2015).