Exploiting orbital constraints from optical data to detect binary gamma-ray pulsars

It is difficult to discover pulsars via their gamma-ray emission because current instruments typically detect fewer than one photon per million rotations. This creates a significant computing challenge for isolated pulsars, where the typical parameter search space spans wide ranges in four dimensions. It is even more demanding when the pulsar is in a binary system, where the orbital motion introduces several additional unknown parameters. Building on earlier work by Pletsch&Clark (arXiv:1408.6962), we present optimal methods for such searches. These can also incorporate external constraints on the parameter space to be searched, for example from optical observations of a presumed binary companion. The solution has two parts. The first is the construction of optimal search grids in parameter space via a parameter-space metric, for initial semi-coherent searches and subsequent fully-coherent follow-ups. The second is a method to demodulate and detect the periodic pulsations. These methods have different sensitivity properties than traditional radio searches for binary pulsars, and might unveil new populations of pulsars.


INTRODUCTION
The Large Area Telescope (LAT; Atwood et al. 2009) on the Fermi satellite has helped to increase the known Galactic population of gamma-ray pulsars to more than 250 pulsars 1 (for a review see, e.g., Caraveo 2014). However, in the recent Fermi LAT Fourth Source Catalog (4FGL; Abdollahi et al. 2020) 1,525 out of 5,098 gamma-ray sources remain unassociated. Many of those are thought to be pulsars, perhaps in binary systems.
Gamma-ray pulsars may be detected in three ways. (a) A known (radio or X-ray) pulsar position and ephemeris guides a follow-up gamma-ray pulsation search within a nearby LAT source (e.g., Abdo et al. 2009a,b;Guillemot et al. 2012). (b) A similar gamma-ray pulsation search is done for a known pulsar, but without an obvious gamma-ray source being present (Smith et al. 2017). (c) A "blind" search hunts for gamma-ray pulsations around a LAT source where no pulsar has yet been identified.
Blind gamma-ray searches are the focus of this paper. Such searches have discovered more than 50 young pulsars (YPs) (e.g., Abdo et al. 2009c;Saz Parkinson et al. 2010;Clark et al. 2017), and three millisecond pulsars (MSPs) (Pletsch et al. 2012a;Clark et al. 2018). Many of these pulsars could not have been found via radio or X-ray emissions, which were not detected in extensive follow-up searches. Such systems are of particular interest because they constrain models of pulsar emission and beaming. Blind searches also have the potential to discover new populations of pulsar/neutron star objects.
So far, most blind gamma-ray searches have targeted isolated pulsars. The searches are a substantial computing effort, and have been carried out in campaigns or surveys that last several years. More recent surveys find new systems because the ongoing LAT operations provide additional data, which enables the detection of weaker pulsations (e.g., Clark et al. 2017). However, there is also a downside: the computing power required also increases quickly with longer observation time spans.
Until now, blind gamma-ray searches have only found one binary MSP, PSR J1311−3430. This is tantalizing because three quarters of the known MSPs in the Australia Telescope National Facility (ATNF) Pulsar Catalogue 2 (Manchester et al. 2005) are in binaries. So if search sensitivity were not limited by computing power, it might be possible to find many more. But even for isolated pulsars it is expensive to search for high (> 100 Hz) spin frequencies, and adding (at least three) additional orbital parameters makes it even more costly. By improving the techniques, the methods presented here are a first step towards finding more of these systems.
Much of our focus is on binary pulsars in so-called "spider" systems, in which the pulsar companion is being evaporated by an energetic pulsar wind. A typical example is the first "black widow" pulsar to be discovered, PSR B1957+20 (Fruchter et al. 1988). This was found in radio, where pulsations are eclipsed for a large fraction of the orbit, presumably by material ablated from the companion. Spider pulsars are categorized as black widows if the companion mass M c is very low (M c 0.1 M ) or as "redbacks" (another spider species) for larger companion masses (M c ∼ 0.15 − 0.4 M ) (e.g., Roberts 2013).
For many of the known MSPs in spider systems, the companions are visible in the optical. The light originates from nuclear burning, and/or from pulsar wind heating up the companion. The orbital motion of the companion then leads to a detectable modulation of the orbital brightness. The source of this modulation is not well understood. It might be that the side of the companion facing the pulsar is hotter than the other side, and is more visible at the companion's superior conjunction. The companion might also be tidally elongated into an ellipsoid, whose projected cross section onto the line of sight varies over the orbit.
The new blind search methods presented here are well suited to gamma-ray pulsars in spider systems, with nearly circular orbits (eccentricity e < 0.05) and for which optical observations of the pulsar's companion provide information about the orbital motion, and thus constrain the gamma-ray pulsation search space.
For concreteness, we present the search designs for two promising gamma-ray sources: (a) 3FGL J1653.6−0158, a likely MSP in a circular binary (Romani et al. 2014;Kong et al. 2014), and (b) 3FGL J0523.3−2528, a probable MSP in a slightly eccentric binary (Strader et al. 2014). These are ranked among the most likely pulsar candidates (Saz Parkinson et al. 2016). We demonstrate the feasibility of a search using the computing resources of the distributed volunteer computing project Einstein@Home (Allen et al. 2013). 2 http://www.atnf.csiro.au/research/pulsar/psrcat The paper is organized as follows. Section 2 reviews blind search methods for isolated gamma-ray pulsars and introduces the concepts required for such searches. Section 3 extends the methods to gamma-ray pulsars in circular orbit binaries and Section 4 further extends these to eccentric orbit binaries. In Section 5 our methods are compared with alternatives used in radio and gravitational-wave astronomy. Finally, in Section 6 we discuss the feasibility of future blind searches for binary gamma-ray pulsars and also consider some specific sources. This is followed by Appendices A, B, and C containing some technical details.
In this paper, c denotes the speed of light and G denotes Newton's gravitational constant.

BLIND SEARCHES FOR GAMMA-RAY PULSARS
Blind-search methods for isolated gamma-ray pulsars have been studied in detail by Pletsch & Clark (2014). Here, we summarize and extend their framework. The following Sections generalize the search methods to binary pulsars.
The search for gamma-ray pulsations begins with a list of N photons from a posited source, which we label with the index j = 1, . . . , N. The data available for these photons are their detector arrival time t j , their direction of origin, and their energy, spanning an observation interval T obs .
We are dealing with many sums and products in this paper. Sums and products over j, k, run from 1, . . . , N unless otherwise specified. Furthermore, we adopt the notation Not all photons are equally significant. Photons at low energies are less well-localized than those at higher energies, and cannot be so readily attributed to a target source. Photons whose energy is more consistent with a distributed background are less likely to come from the pulsar. Photons originating from a nearby point source might contaminate the data set. For such reasons, searches may be improved by modeling the spatial and energy distribution of the sources. This assigns a weight w j ∈ [0, 1] to each photon, which is the only place where the energy and arrival direction of the photons enter our analysis. The weight w j represents the probability that the j'th photon originated at the nominal pulsar (Bickel et al. 2008;Kerr 2011). These weights are used for noise suppression and to reduce computing cost by removing the lowest-weighted photons. In this paper, we assume that these weights have been determined in advance for each photon, so the only information available for the j'th photon is its arrival time t j in the detector and the weight w j .
The question that we need to answer is, are the arrival times of these photons random, or is there an underlying periodic-ity? To answer this question (in the statistical sense), we first need a model for the periodicity, which we assume is tied to the physical rotation of the pulsar.

Pulse profile and photon arrival probability
For now, assume that "in isolation" the pulsar would have a linearly-changing angular velocity. Using Φ to denote the rotational phase in radians where t psr is the time that would be measured by a fictitious observer freely falling with the center of mass of the pulsar, and t ref is a reference time. Note that detector time ticks at a different rate than t psr , because the detector is moving around the Earth and Sun, and because the pulsar might be orbiting a binary companion, or accelerating towards the Galaxy. Also note that without loss of generality we have set the phase at the reference time to zero. The parameters λ describe the pulsar. Here they are the spin frequency f and its first time derivativeḟ at reference time t ref . This second-order Taylor approximation holds for many pulsars and most MSPs, but for very young and "glitching" pulsars, additional higher order terms may be needed.
The flux of photons can be broken into three parts. The first does not come from the pulsar: it is a background which is uncorrelated with pulsar rotation. We call these unpulsed photons "background". The second part originates from the pulsar itself but is also uncorrelated with pulsar rotation. We call these "unpulsed source" photons. The last part is a periodically time-varying flux from the source, which we call "pulsed". We use p to denote the ratio of the number of pulsed photons to the total number of source photons (pulsed and unpulsed source).
The pulsed photon flux may be described with a periodic function F S (Φ) of the pulsar's phase around its rotational axis, Φ ∈ [0, 2π], and is time-stable for most pulsars. The normalized probability that a pulsed photon arrives in the phase interval [Φ, Φ + dΦ] is F S (Φ) dΦ. The function F S (Φ) has minimum value zero and encloses unit area in the interval [0, 2π].
We can now give the probability density function for the rotation phase associated with a given photon. This differs from one photon to the next because photons with small weight w j are more likely to have a phase-independent probability distribution. The probability that the j'th photon originates from a rotation phase interval [Φ j The first term (with probability 1 − w j ) describes the background photons, and the second and third terms (with prob-ability w j ) the unpulsed and pulsed source photons respectively. The probability distribution of pulsed photons may be expressed as the Fourier series The complex Fourier coefficients are Note that the Fourier coefficients γ n are constrained because F S has minimum value zero. Note also that for known gamma-ray pulsars |γ n | 2 decreases quickly with increasing index n (Pletsch & Clark 2014). In many cases the first 5 harmonics are sufficient to describe the pulse profile.
In principle, to detect gamma-ray pulsations, we assume a rotational model f ,ḟ and then compute the rotational phase associated with each photon. "Binning" these phases (mod 2π) with weights w j provides an estimate of F(Φ) = j w j F j (Φ)/ j w j , from which we can estimate F S (Φ) by shifting the minimum value to zero and re-scaling to unit area. If that function is compatible with zero (meaning: coefficients γ n are small), then no pulsations were detected. Conversely, if the γ n are large for some values of f andḟ , we have found pulsations.

Relationship of detector time t to t psr
The situation is slightly more complicated than described in the previous paragraph because computing t psr for each photon from its time of arrival at the Fermi satellite also requires the pulsar's sky position (right ascension α and declination δ). The sky position allows for " barycentric corrections", e.g., to account for Doppler shifts due to the LAT's movement around the Solar System Barycenter (SSB). Thus the photon's emission time t psr (t, α, δ) is a function of its arrival time t at the LAT and the putative pulsar's sky position. The pulsar's putative phase is a function of t and the four parameters λ = { f ,ḟ , α, δ}.
In blind searches the spin parameters are unknown. Although each photon is tagged with an arrival direction α, δ, these are not sufficiently precise to detect pulsations, so those location parameters must also be searched. Hence the parameter-space search volume Λ for isolated pulsars (λ ∈ Λ) is 4-dimensional. In Sections 3 and 4, the higherdimensional search spaces for binary pulsars in circular and elliptical orbits are discussed.

Searching for pulsations
For realistic searches the parameter space Λ is too large to search by the straightforward computational process described above. Instead Λ is explored with a multistage search based on several different test statistics (e.g., Meinshausen et al. 2009). This gives the greatest sensitivity at fixed computational cost (Pletsch & Clark 2014). The approach is hierarchical. In the first stage, a coarse grid covering the parameter space Λ is searched at low sensitivity using inexpensive test statistics. These are relatively insensitive to mismatch between tested parameters and pulsar parameters. In the following stages, smaller regions of Λ around the most promising candidates are searched at higher sensitivity. These use more expensive test statistics on finer, more closely spaced grids. Thus, a search is defined by a test statistic/grid hierarchy.
The spacing of the grids in parameter space is governed by the mismatch described above. For a given test statistic, we calculate a "metric", which is the fractional loss in the expected signal-to-noise ratio (S/N). The details of this are found later in this Section.
The search described in this paper has four stages, which employ detection statistics P 1 , S 1 , and H. Here we briefly describe the overall structure. The test statistics are defined and characterized later in this Section.
The first three stages search for significant power in the first harmonic |γ 1 | 2 . Each discards regions of parameter space which contain no signals; what remains is passed to the following stage. The first stage uses the "semi-coherent" test statistic S 1 with a low threshold. The second stage tests S 1 on a finer grid, with a higher threshold. The third stage uses the fully coherent test statistic P 1 . This searches coherently for power |γ 1 | 2 over the full observation span T obs with much greater sensitivity and a finer grid than before.
The fourth stage employs the expensive H statistic, which combines P 1 , . . . , P 5 . This coherently integrates over T obs to identify power in the first five harmonics |γ 1 | 2 , . . . ,|γ 5 | 2 . By searching around the surviving candidate points in parameter space with a still finer grid, this completes the hierarchy.

Coherent power test statistic P
The basis for all of our test statistics is the coherent Fourier power, evaluated over different periods of time. For the n'th harmonic, and including all of the photons, this is To simplify notation, from here on we use Φ(t j , λ) to denote Φ(t psr (t j , α, δ), f ,ḟ ), where t j is the photon arrival time measured at the LAT. The normalization constant is How does P n behave in the absence of pulsations and in the presence of pulsations?
To answer this question, we compute expectation values as shown in Appendix A. The power P n has an expected value (Eq. A5) and variance (in the absence of a pulsed signal, p = 0) The power P n is a detection statistic because it is sensitive to a non-vanishing pulse profile. If γ n is non-zero, then P n should be larger than two. It becomes larger as the fraction p of pulsed to source photons increases (which we cannot control). It also becomes larger as the number of photons (or equivalently, the observation time) grows. But to understand what values of P n correspond to statistically significant detections, we need to know about its statistical fluctuations, meaning the variance in P n .
Note that the diagonal-free double sum in these expressions can be re-expressed as ( j w 2 j ) 2 − j w 4 j . Thus the variance can be written If there are many photons from the source, and the weights are relatively uniformly distributed, then it follows that the numerator in Eq. 10 is O(N) and the denominator is O(N 2 ). Hence, the variance Var 0 [P n ] → 4 − O(1/N) approaches 4. In this limit, and with the statistical assumptions of Appendix A, P n has a non-central χ 2 -distribution with two degrees of freedom (Pletsch & Clark 2014). The non-centrality parameter is the second term appearing in Eq. (8). The expected S/N associated with P n is In the many-photon limit the quantity µ → j w 2 j /T obs is proportional to the mean weighted photon arrival rate.

Loss of P from parameter mismatch
In a real search, we compute detection statistics at a grid of discrete values of the signal parameters λ. If there is a signal present, its actual (true) parameters might be close to one of these discrete values, but will not match it exactly. There will always be some offset between the tested parameters and the true parameters. Here we quantify how much signal-to-noise is expected to be lost because of this mismatch.
Assume that the tested parameters λ are close to the true pulsar parameters λ psr , and introduce the notation dλ a = λ a − λ a psr (12) for the small parameter offsets. Here and elsewhere in the paper we index the parameter space dimension with lower-case Latin letters "a" and "b". These offsets change the pulsar rotation phase by where the notation is introduced and we neglect higher powers in dλ. We also adopt the Einstein summation convention that repeated parameter space indices are summed over all the dimensions of the parameter space. We now compute the fractional loss in expected S/N associated with this parameter mismatch. For the offset parameters the coherent power is where Φ j = Φ(t j , λ psr ) and ∆Φ j = ∆Φ(t j ). Following Appendix A, the expectation value of this is It follows that for the mismatched signal the expected S/N is The fractional loss in S/N (often called the "mismatch") is We need the mismatch to help set the spacings of the parameter space search grids, but for that purpose, approximations suffice. Assume there are many photons, and the weights are uniformly distributed in time (or at least slowly varying in a way that is not correlated with the pulsar rotation phase). The sums over the weights may then be replaced with simple integrals over time, giving Here we introduce "angle bracket" notation for an average over a time interval of length T centered around an arbitrary time t 0 . This takes an input function Q(t ) and outputs a new function of time t defined by which is the average of Q around the time t.

Parameter space metric g ab
Since the sensitivity of these searches is limited by available computing power, we need to construct a grid that covers the relevant parameter space with the smallest number of grid points. This means that the parameters λ psr of any possible pulsar should be close enough to a grid point that we do not lose too much S/N from the mismatch, but the grid should have as few points as possible.
The distance metric on the search space is a useful tool for such constructions (Balasubramanian et al. 1996;Owen 1996). It provides an analytical approximation to the mismatch. For example, the coherent mismatch in Equation (19) can be approximated by the "coherent metric" g ab m(λ, λ psr ) = n 2 g ab (λ) dλ a dλ b + O(dλ 3 ) for small coordinate offsets dλ a from the true pulsar parameters.
Expanding the exponential that appears in Eq. (19) to first order, one finds To evaluate the metrics, we need to account for the way in which the detected pulsar rotation phase depends upon the different pulsar parameters.

Evaluation of g ab for isolated pulsars
As seen by an observer freely falling at the center of mass of the pulsar, the rotation phase just depends upon the intrinsic frequency f and its derivativeḟ as given in Eq. 2. But as explained in Sec. 2.2, these must be converted to detector time.
For computing the metric, we do not need a conversion that is accurate to microseconds, but only one that takes into account the largest shifts between detector and pulsar time, of order ≈ 500 s, arising from the motion of the Earth around the Sun (Pletsch & Clark 2014). We denote the orbital angular frequency Ω E = 2π/yr, the orbital light-crossing time by r E = 1 AU/c, and the obliquity of the ecliptic by = 23.4 • .
If we choose a coordinate axis z along the line of sight to the pulsar, the the projected motion is where n x = cos α cos δ , n y = cos sin α cos δ + sin sin δ , and the sky location is given by the right ascension α, and the declination δ. The (arbitrary) choice for the origin of the time coordinate determines the constant ϕ ref , which is the Earth's orbital phase at that moment. Note that this simplified version of the Rømer delay does not account for the motion of the Fermi satellite around the Earth. It is not accurate enough to use in a search for pulsations, and is only used in the metric calculation.
For the purpose of computing the metric we can model the detected pulsar rotation phase as the sum of (2) and the additional phase cycles introduced by the Rømer delay (23): Here, the search parameters are λ = { f ,ḟ , n x , n y }, and the terms correcting the arrival times t have been neglected for theḟ summand. The metric for the coherent power P 1 follows from Eq. (22). The formulae are complicated but if we keep only the most significant terms then they simplify. To determine these, consider the relative size of the different quantities: T obs ≈ 10 yr ≈ 3 × 10 8 s , Most MSPs have parameters f andḟ in the given range. With these in mind, one finds diagonal metric components Most of the off-diagonal metric components are negligible. Determining if off-diagonal metric components are significant requires some care because they need to be compared to the corresponding diagonal components. This arises here, and in several other places in the paper. Here, we show in detail how this significance is determined. The same reasoning is used for the other cases that arise later, but is not elaborated.
Since the fundamental quantity of interest is the mismatch m, for fixed a and b (no Einstein summation convention), consider m = g aa (dλ a ) 2 + g bb (dλ b ) 2 + 2g ab dλ a dλ b . Re-scale the coordinates {λ a , λ b } to new coordinates {λ a = uλ a , λ b = wλ b } such that the two diagonal components of the metric in the new coordinates are both unity. (Here, u and w denote the re-scaling factors.) This implies that g aa (∂λ a /∂λ a ) 2 = g aa u −2 = 1 and g bb (∂λ b /∂λb ) 2 = g bb w −2 = 1. Then all offdiagonal metric components are of O(1/Ω E T obs ), apart from Note that all the off-diagonal terms may be neglected in the case that the integration time T obs 1 yr and the reference time t ref = t 0 .
For this case the diagonal "coherent metric" terms reduce to

Semicoherent power test statistic S
The coherent power P n in (6) provides a good statistical basis to find pulsations (meaning γ n nonzero) but is inefficient to compute. So the first two stages of our searches use the "semi-coherent" Fourier power S n . Its definition is similar to P n except that photons are only combined if their arrival-time difference is smaller than a coherence time, T coh T obs . This makes it less expensive to compute (but also less sensitive). The coherence time in a typical search in the first stage is T coh = 2 21 s ≈ 24 d, in the second stage is T coh = 2 22 s ≈ 48 d, and the observation span T obs (i.e. the operation time of the LAT) is more than 10 years.
For convenience the statistic S n differs from P n in one other way: we omit the diagonal j = k terms in the sum. This ensures that in the no-signal (p = 0) case the expected value of S n vanishes, with The rectangular window function restricts the sum to photons in which the arrival time difference τ jk = t j − t k (or "lag") is not larger than T coh : The semi-coherent normalization constant is chosen to bē which ensures that in the no-signal (p = 0) case S n has unit variance (Clark et al. 2017).
To characterize this detection statistic, we calculate the expectation value and variance with the calculational framework of Appendix A, obtaining: The expectation value is the same as the second term of P n in Eq. (8), except that the sum is restricted to the lag window. In fact, the formulae above hold for any choice of window function. The S/N for the semi-coherent Fourier power S n is simplified by assuming a rectangular window function (which equals its square). This gives The second line adopts the definition of µ given after Eq. (11) and makes the same assumptions of steady photon flux and large photon number. In practice, how large are these detection statistics? A typical gamma-ray pulsar might have a pulsed flux for which |γ 1 | 2 ≈ 0.2 and a 70%-fraction of pulsed photons for which p 2 ≈ 0.5. The weighted flux of source photons detected might be j w 2 j ≈ 500 over T obs = 10 yr, implying a rate µ ≈ 50/yr. With T coh = 24 d, this leads to coherent and incoherent S/Ns of order θ 2 P1 ≈ 50 and θ 2 S1 ≈ 4, significant at the 50σ and 4σ levels respectively.

Loss of S from parameter mismatch
We now turn to the metric for the semi-coherent statistic. To compute the mismatch for the semi-coherent detection statistic S n , with the same assumptions as above, we can replace the sums with integrals, obtaininḡ Note that the inner integral in the second line can include times outside the observation span t ∈ [t 0 − T obs /2,t 0 + T obs /2], going down to t = t 0 − T obs /2 − T coh /2 or up to t = t 0 + T obs /2 + T coh /2. In such cases the integrand should be set to zero, and normalized so that 1 = 1.

Parameter space metricḡ ab
We now evaluate these mismatches to lowest order, obtaining a distance metric on the parameter space. We evaluate the integrals in Eq. (37) naively, without setting the integrands to zero outside of the "valid data range". This gives rise to terms (complex or linear in dλ a ) which are not present in the exact expression. We assume that T coh T obs (typically T coh = 24 d and T obs > 10 yr). In that case, these terms are small, and we discard them.
The partial derivatives with respect to λ a ∈ { f ,ḟ , n x , n y }, under the assumption that T coh 1 yr T obs , can be approximated as as Pletsch & Clark (2014) did. (Here and in what follows, for readability, the time dependence of phase derivatives such as ∂ a Φ is not shown explicitly.) With these assumptions the semi-coherent mismatch Eq. (37) can be approximated by the semi-coherent metric where dλ a = λ a − λ a psr as earlier. Note that Eq. (39) has the same form as the coherent mismatch in Eq. (21).
The metric components arē where we have introduced which is exactly the coherent metric given in Eq. (22), but with T obs replaced by T coh 1 yr and t 0 replaced by t . Thus, terms of O(1/Ω E T coh ), similar to those appearing in Eq. (28), cannot be neglected.

Evaluation ofḡ ab for isolated pulsars
The non-vanishing semi-coherent metric components are: The semi-coherent metricḡ is diagonal for t 0 = t ref , as was the case for the coherent metric g in Eq. (30). Note that g nxnx and g nyny are not equal because the neglected terms of O(1/Ω E T obs ) have opposite signs. The metric componentḡ˙f˙f differs from that given by Pletsch & Clark (2014), but our results are identical in the limit of a large number of photons N homogeneously distributed over the observation span. This is the case, since we assumed it in deriving Eqs. (19) and (37).
Comparison of Eqs. (28) and (42) illustrates the benefits of the multistage search process described in Section 2.3. For grids with the same mismatch m =m, the ratio between the density of the coherent grid and the semi-coherent grid would be coherent grid density semi-coherent grid density For the timescales T coh and T obs given above, the ratio ∼ 10 6 . This is why the semi-coherent search stage is beneficial.

Multiple harmonic test statistic H
In the last and most sensitive stage of the multistage search, we adopt the widely used statistic which incoherently sums the coherent power from up to the first M max harmonics in the pulse profile. The H statistic provides a sensitive test for unknown (generic) pulse profiles. The original simulations by de Jager et al. (1989) recommended M max = 20, and to assess the false-alarm probability, carried out a numerical study of the distribution of H in pure noise.
Later results by Kerr (2011) show that the single-trial probability ρ of exceeding a value H threshold in pure noise is well modeled by ρ ≈ exp(−0.398 H threshold ) if the number of harmonics M max is very large. Obviously, if M max is reduced, then the single-trial probabilities are smaller than this, so exp(−0.4 H threshold ) is a reliable upper bound.
To avoid over-fitting, we generally use smaller limits M max = 3, 4, or 5 on the number of harmonics. Typical blindsearch gamma-ray pulsar detections have H values in the hundreds, corresponding to single-trial ρ values which must lie below 10 −30 .
Normally, the last search stage is not computationally limited. So we use a grid fine enough to secure power in the higher harmonics, while over-covering the search space for power in the lower harmonics. In practice, the grid is built using the coherent metric presented in Section 2.4.2 with n = M max .

Searches for isolated pulsars
Blind searches for isolated pulsars within gamma-ray data recorded by the LAT have been very successful (see, e.g., Clark et al. 2017). The key ingredients are the utilization of the powerful volunteer distributed computing system Ein-stein@Home (Allen et al. 2013) and searches which use these computing resources as efficiently as possible.
Most of the tools for constructing efficient searches have been presented in the earlier Sections. To discard unpromising regions in parameter space, the multistage approach is used as described in Section 2.3. For the first and computationally most crucial search stage, efficient grids covering the parametersḟ , n x , n y are built based on the distance metric, and f is searched using fast Fourier transform (FFT) algorithms (Frigo & Johnson 2005). In later search stages, f is also gridded with the metric but it is not efficient to use FFTs on the small ranges in f around the few most significant candidates from the semi-coherent search stage.

SEARCH METHOD: CIRCULAR BINARY ORBITS
The main problem in blind searches for pulsars is that the phase model from equation (2) depends on the (photon emission) time at the pulsar, while a gamma-ray detector records the time of arrival at the telescope. For binary pulsars the largest corrections to shift between these two times arise from the line-of-sight motion of the Fermi satellite around Earth and Sun r z,sky (t) and of the pulsar around its companion r z,cir (t psr ).
The line-of-sight motion of a binary pulsar in a circular orbit can be described via 3 parameters, which are usually taken to be: the orbital frequency Ω orb ; the projected semimajor axis x in seconds; and the epoch of ascending node T asc . With these, the two times are related by t psr + r z,cir (t psr ) = t + r z,sky (t) , where the corrections, also called Rømer delays, are expressed in seconds. The simplest expression of the pulsar's orbital line-of-sight motion r z,cir depends on the time measured at the pulsar t psr . In many cases, this time may be replaced with the detector time because and the quantity xΩ orb 1. In such cases This holds for most black widow and some redback systems with projected semimajor axes in the order of a few light seconds (see, e.g., ATNF Pulsar Catalogue 3 by Manchester et al. 2005). In all cases, it is accurate enough to compute the metric, and in many cases accurate enough for maintaining phase coherence in a search. The Rømer delay can be expressed in terms of the three orbital parameters as Here the orbital frequency Ω orb is connected to the orbital period P orb via P orb = 2π/Ω orb . In gamma-ray searches, in addition to the Rømer delay, we also have to correct for other effects like the Shapiro and Einstein delays. In contrast to radio observations, we do not have to account for the frequency-dependent dispersion caused by the Interstellar medium (ISM) because gamma rays are well above the plasma frequency of the ISM.
All of these effects are described by Lorimer & Kramer (2004) and Edwards et al. (2006). While these corrections must be included in gamma-ray searches, only the largest effects need to be included in the phase model for the derivation of a distance metric approximation.

Parameter-space metrics
In order to compute the metric, a simplified phase model can be used which accounts for the corrections (23) and (48): Here the search parameters are λ = { f ,ḟ , n x , n y , Ω orb , x, T asc } and the terms correcting the arrival times t have been neglected for theḟ summand. This phase model is not sufficient for searches because it would not maintain phase coherence 3 http://www.atnf.csiro.au/research/pulsar/psrcat with a true pulsar signal. However, it is sufficient to describe how varying the signal parameters leads to loss of S/N. The dominant components of the coherent metric for the orbital parameters are where we have assumed that the integration time span T obs is much larger than the orbital period P orb . Compared to the diagonal terms, as done in the text below Eq. (28), all other components are of O(1/Ω orb T obs ). The off-diagonal component g ΩorbTasc is vanishingly small if the epoch of the ascending node is close to the middle of the gamma-ray data set, T asc ≈ t 0 . In principle, T asc can be shifted forwards or backwards by an integer number N of orbital periods P orb to achieve this. However, when T asc is constrained, for example by optical observations, this is undesirable because it introduces uncertainties in the shifted value of T asc that grow linearly with N .
Even if T asc ≈ t 0 , our current searches ignore the offdiagonal term in the metric. The only negative consequence is that the grids are more closely spaced than needed, which reduces the efficiency of the search.
If we include the additional orbital parameters, the semicoherent mismatch (37) can still be written in metric form However, the assumptions made previously in (38) to calculate this only hold for the "isolated pulsar" parameter space coordinates λ iso = { f ,ḟ , n x , n y }. They do not hold for the additional orbital parameters λ orb = {Ω orb , x, T asc }.
If λ a ∈ λ orb is an orbital parameter and P orb T coh (typical coherence time T coh ≈ 24 d) the approximations are valid. By this, we mean that the ratio of the resulting metric to the correct metric is 1 + O(P orb /T coh ). With these assumptions the semi-coherent metricḡ ab is composed of three types of components. For the first type, the parameters λ a , λ b ∈ λ orb are orbital. For these components,ḡ giving the coherent result from Eq. (50).
For the second type, the parameters λ a , λ b ∈ λ iso are isolated. For these components which is the semi-coherent result found in Eq. (42). For the third type, one of a or b is in λ orb and the other is in λ iso . One obtains the same equation as for the second type. This vanishes by virtue of Eq. (52) and because In short, the non-vanishing semi-coherent metric components reduce to earlier results. For the orbital parameters, they are the same as the coherent metric components. For the isolated (spin and celestial) parameters, they are the same as the semi-coherent metric components for an isolated pulsar. To reiterate, the non-vanishing semi-coherent orbital metric components are: As before, for epoch of ascending node close to the middle of the dataset, i.e. T asc ≈ t 0 , the semi-coherent metric is diagonal.
At the end of Section 2.5.3, we discussed the relative densities of the coherent and semi-coherent grids for isolated sources. Now we have added three additional (orbital) dimensions to the parameter space. Because the metric factors into a product of a metric on the orbital parameters and a metric on the isolated parameters, the grid may also be constructed as a product of the grids on the corresponding subspaces. For the isolated parameters, the ratio between the density of the coherent grid and the semi-coherent grid is the same as for the search for isolated pulsars. For the orbital parameters, the number of grid points needed is the same as in the coherent case. Hence, the ratio of grid densities is the same as in Eq. (43).
In Figure 1 the mismatch and its coherent-metric approximation are compared for small parameter offsets, for a realistic simulated pulsar. The corresponding plot for the semicoherent mismatch looks very similar but has different fandḟ -scales. The mismatch and its metric approximation agree well for mismatch m ≤ 0.4. This is a typical value for a search: in Appendix B, we show that maximum sensitivity for a given computing resource is obtained for an average mismatchm = 0.383 (see Table 3).
The celestial parameters are not shown in Figure 1; for spider pulsars they are usually known to high precision from optical observations (e.g., from the Gaia DR2 Catalog; Gaia Collaboration et al. 2018), so no grid is required. For other pulsars where the sky position is less constrained a grid may be needed.
The search ranges for the orbital parameters are very large, and without further knowledge a full blind search is not possible. On the other hand, some searches are possible if the pulsar's companion is visible in optical/X-ray observations, which constrains the search parameters. In the next Section, we discuss a gamma-ray pulsar search design for 3FGL J1653.6−0158, which is thought to be an MSP in a circular-orbit binary (Romani et al. 2014;Saz Parkinson et al. 2016).

Search design for circular binary
This Section shows how to reduce the binary-pulsar search parameter space by exploiting orbital constraints from the companions.
We use the gamma-ray source 3FGL J1653.6−0158, which is predicted to be a spider pulsar (Romani et al. 2014;Kong et al. 2014), as an example. It was ranked second in Saz Parkinson et al.'s list (published 2016) of the most significant Fermi LAT Third Source Catalog (3FGL) unassociated sources predicted to be pulsars. The list also classifies it as a likely MSP. The gamma-ray source 3FGL J1653.6−0158 shows typical pulsar properties: a time-stable photon flux and a spectrum described by an exponential-cutoff power law.
The search ranges in spin frequency f and spin-down parameterḟ are guided by the known pulsar population and computational constraints. The search range is divided into YPs, with lower frequencies ( f < 44 Hz), and MSPs, with higher frequencies (44 Hz < f < 1500 Hz) 4 . Correspondingly, the spin-down lies between 0 and −10 −10 Hz s −1 , for YPs, and between 0 and −10 −13 Hz s −1 for MSPs.
The constraints for f andḟ define a region in parameter space that has to be searched. The frequency dimension can be efficiently scanned using the FFT algorithm (Frigo & Johnson 2005) as described by Pletsch & Clark (2014); Clark et al. (2016) and Clark et al. (2017) for isolated pulsars. Thė f -dimension can be covered by a uniformly spaced lattice. Special treatment for these parameters is possible: since their metric components are independent of the other parameters, so is the spacing.
In practice, the FFTs are computed in frequency intervals of bandwidth f BW = 8 Hz. These have f BW T coh frequency grid points, with frequency spacing 1/T coh . In the semi-coherent stage, for two points separated by half the grid spacing, this gives a worst-case metric mismatch m = π 2 /24 ≈ 0.411. (As . A comparison of the coherent metric approximation to the actual mismatch, for parameters of a simulated circular-orbit binary pulsar in 3FGL J1653.6−0158. Blue contours show the actual mismatch and red contours the metric approximation, at m = 0.2, and 0.4. As is generally the case (Allen 2019), the metric contours are conservative and lie inside the actual mismatch contours. discussed in Appendix B following Eq. (B7), this can be reduced by interpolation to a worst case value of m = 0.14, at no significant cost.) Thus, for one f BW -interval, the computing cost is the product of the cost of a single FFT multiplied with the number of parameter-space grid points in the other dimensions.
The sky position is tightly constrained because a likely optical and X-ray counterpart with significant light-curve modulation was found (Romani et al. 2014;Kong et al. 2014;Hui et al. 2015), and proposed to be an irradiated pulsar companion. At the time, the best estimate for the position of the likely optical counterpart was from the USNO B1.0 Catalog (Monet et al. 2003). Using this instead of the 3FGL position makes it possible to search 3σ-ranges of the sky parameters with only one semi-coherent sky grid point. At high frequencies extra sky grid points are needed only in the follow-up stages. The computing costs of these are negligible com-pared to the semi-coherent stage. The same optical source can now be identified in the Gaia DR2 Catalog (Gaia Collaboration et al. 2018), see Table 1. For this, the uncertainty in sky position is small enough that even at f = 1.5 kHz no extra sky points are needed.
The orbital parameters Ω orb and T asc are directly constrained by Romani et al. (2014) using optical observations of the companion. As shown in Table 1, they found a significant modulation at a period of P orb = 0.05194469 ± 1.0 × 10 −7 d, with epoch of ascending node T asc = 56513.48078 ± 5.2 × 10 −4 MJD.
Additional observations allow the third orbital parameter, the projected semimajor axis of the pulsar x = a 1 sin i/c (in units of light travel time), to be constrained. Here we denote the neutron star with subscript "1" and the companion with subscript "2". Measurements of the companion's velocity amplitude K 2 = 666.9 ± 7.5 km s −1 , together with the orbital period, imply that the pulsar mass function has the value where the mass ratio is q = M 2 /M 1 . This implies that the neutron star has mass M 1 > 1.60 ± 0.05 M . Since redback companions have masses M 2 0.4 M (Roberts 2013), this in turn implies q < 0.25. From Eq. (56), a mass ratio of q = 0.25 allows neutron-star masses up to 2.5 M for i = 90 • . (This is reassuringly conservative, since the most massive known neutron star (Cromartie et al. 2020) has mass 2.14 M .) Combining the mass function with Kepler's third law (a 1 + a 2 ) 3 = G(M 1 + M 2 )(P orb /2π) 2 and the center of mass definition a 1 M 1 = a 2 M 2 gives The upper limit for q then implies an upper limit x 0.2 s. It is challenging to build a search grid that covers the 3dimensional orbital parameter space with as few points as possible. This is because (as can be seen from the metric) the orbital parameter space is not flat, so a constant-spacing lattice is not optimal. A solution to this is presented by Fehrmann & Pletsch (2014), starting with "stochastic search grids" (Babak 2008;Harry et al. 2009). A stochastic grid is built by placing grid points with a random distribution that follows the expected distribution of metric distances, while ensuring a preset minimum distance between them. The resulting grid is then optimized by nudging grid points towards regions where neighboring grid points have higher than average separation. The resulting search grid is efficient, and has a well behaved mismatch distribution, which simplifies the S/N distribution in the absence of signals.
The minimum number of grid points needed to cover the orbital parameter search space at mismatch m can be estimated from the proper 3-volume Here, the integral is over the relevant range of orbital parameter space, g denotes the orbital metric from Eq. (50), and numerical factors of order unity related to the efficiency (technically "thickness"; see Appendix B) of the grid lattice have been dropped. To understand how this depends on parameters, note that the integral is proportional to where the search range for x is [x min , x max ]. ∆Ω orb and ∆T asc are the search ranges around the values of Ω orb and T asc estimated from the optical modeling. Furthermore, we make the assumption that ∆Ω orb Ω orb . The strong dependency of Derived search range Projected semimajor axis a , x (s) . . 0 -0.2 NOTE-The JPL DE405 solar system ephemeris has been used and times refer to TDB. a Assuming a mass ratio of q < 0.25, see text following Eq. (56).
N orb on x max and f means that searches for YPs (smaller f ) in tight binary orbits (smaller x max ) are computationally much cheaper than searches for MSPs in wide orbits. The latter are only possible if the orbital constraints are very narrow.
If the parameter space is small in a particular direction, this reduces the effective dimension of the parameter space and changes the formulae above. For example, denote the range of x by [x min , x max ]. Now consider the case where ∆x = x max − x min is small enough that g xx ∆x 2 m. Then only a single grid point is needed in the x-direction, and Eq. (58) must be replaced with a two-dimensional integral, and the exponent on m replaced with −1. Since the orbital metric components in Eq. (55) depend on the parameters, for example g xx = 2π 2 f 2 , this reduction in dimension can take place for certain ranges of parameters (here small frequency f ) and not for others.
We can estimate the computing cost of a search for 3FGL J1653.6−0158 by computing the number of grid points  Table 1 (no grid is needed over sky location). The frequency range is gridded in intervals of bandwidth f BW = 8 Hz as discussed earlier in this Section. The total computing cost is obtained by multiplying the cost of one FFT, the number ofḟ grid points, and the number of orbital grid points (which depends on the f interval), and then summing over the f intervals. Since the orbital grid depends on frequency, a new search grid is constructed for each frequency interval, using the metric at the maximum frequency of that interval.
A convenient way to express the computing cost is in terms of search duration on Einstein@Home, where we assume that the project provides 25,000 GPU-hrs/week. This is shown in Figure 2 as a function of the maximum frequency searched. Searching up to f = 1,500 Hz requires less than 80 d. Note that the search cost in one frequency step is proportional to the number of orbital grid points. To search 3σ ranges in T asc and Ω orb within a reasonable amount of time, either the maximum f or x need to be reduced.
We can also give a general estimate for the MSP-search duration. Since the semimajor axis is typically not well constrained, we assume x min = 0. We evaluate Eq. (59), using Kepler's third law to replace x max with the corresponding maximum searched mass ratio q max , obtaining where M 1 is the neutron star mass. As before, we assumeḟ ∈ [−10 −13 , 0] Hz s −1 for an MSP search. The search duration up to a maximum frequency f max is then where the dimensionless parenthetical factors are of order unity for typical systems of interest, and For redbacks (q < 0.3) one has B(q) < 0.02, whereas for black widows (q < 0.08) one has B(q) < 4×10 −4 . The time A depends on the details of the search and the available computing resources. A typical Einstein@Home search as described in this Section has A ∼ 10 d.
In summary, this Section has shown how the circular-orbit binary pulsar search for 3FGL J1653.6−0158 can be carried out. It is computationally expensive, but by exploiting the orbital constraints it is feasible, even for high MSP frequencies. In practice, a search would start at low frequencies, gradually working up to 1.5 kHz. To further reduce cost, the search should be stopped if a pulsar is found.
While here we have considered one specific example, these methods are more broadly applicable. With them, circular orbit binary pulsar searches are practical if there are good orbital constraints from optically visible companion stars and if the pulsar's projected semimajor axis is not too large.

SEARCH METHOD: ECCENTRIC BINARY ORBITS
For pulsars in eccentric binary orbits, the photon arrival times have to be corrected for the line-of-sight motion r z,ell (t), which is the projection of the eccentric orbit in the line-ofsight direction. In analogy with equation (47), we can approximate the photon emission time at the pulsar as up to O(xΩ orb ). Compared with the circular case, two extra parameters are needed to describe the projected line-of-sight motion, r z,ell (t). For now, we take these to be the orbital eccentricity e and the angle ω between the ascending node and the pericenter. We note that the approximation to O(xΩ orb ) is sufficient for the elliptical example source considered in this paper. If the value of x were larger, a higher order approximation in x would also be required .
YPs with main-sequence stars as companions can have very eccentric orbits. For small orbits the pulsars tidally deform the companion, which dissipates energy. This tidally locks the companion, so that the same side of the companion faces the pulsar, and over time circularizes the orbit (Phinney 1992). This explains why old, spun-up MSPs are usually found in binaries with small or unobservable eccentricity. Only a few exceptions are known (Knispel et al. 2015).
If the energy loss in a spider system is small for each orbit, the pulsar moves around a smaller ellipse and the companion around a larger ellipse. The fixed center of mass is a focus of both ellipses, and the separation vector between pulsar and companion also traces an ellipse.
The line-of-sight variation due to the elliptical motion, r z,ell (t), was derived by Blandford & Teukolsky (1976) and can be written as In this formula the label "ell" is replaced by "BT" to denote that this is the Blandford & Teukolsky model. The eccentric anomaly E is a parameter along the pulsar path which increases with time. If ψ is the angular position of the pulsar measured from the center of the ellipse, then tan ψ = (1 − e 2 ) 1/2 tan E. Equivalently, project the pulsar's position parallel to the semiminor axis, onto a circle whose radius is the semimajor axis, and whose center is the center of the ellipse. Then E is the angular position of that projected point on the circle. E obeys Kepler's equation where M is the mean anomaly. This is a linear function where T 0 = T asc + ω/Ω orb is the epoch of pericenter passage. Unfortunately, there are some problems with the BT model and this parameterization. Kepler's equation (65) cannot be solved in closed form to find E as a function of t. Furthermore, in small-eccentricity orbits, the pericenter is not welldefined and the mismatch arising from offsets in T 0 and ω does not take the simplest possible form. For these reasons, we shift to an uncorrelated set of parameters and Taylorexpand r z,BT as function of e.
A new set of parameters was suggested by Lange et al. (2001). These are the time of ascending node T asc and two Laplace-Lagrangian parameters 1 and 2 defined via The parameters {T 0 , e, ω} are given by With the old parameters, the region of constant mismatch around a grid point is an ellipsoid whose principal directions are not parallel to the {T 0 , e, ω} axes. In the next Section, we show that with the new parameters, the region of constant mismatch is a sphere. This simplifies the code used to optimize grid point locations.
The Rømer delay r z,BT for the pulsar's motion can be expanded to first order in e. Following convention, we use the label "ELL1" for this linear-in-e model: r z,BT = r z,ELL1 + O(e 2 ). This can be described using the parameters {T 0 , e, ω} or the parameters {T asc , 1 , 2 } as We have introduced which is similar to M in Eq. (66) but shifted from pericenter to ascending node. (Note that the term −3e sin ω/2 = −3 1 /2 is typically dropped, as it is time independent.) The ELL1 approximation to the BT model can accurately track the pulsar's rotational phase for eccentricities e below some threshold value. In Appendix C, we show how this threshold depends upon the spin frequency f and semimajor axis x.
Later in the paper, in Section 4.2, we design a search for 3FGL J0523.3−2528, which is a gamma-ray source predicted to harbor a redback pulsar in an eccentric orbit. For that case, the ELL1 model is insufficient and a third-order-in-e model is needed. In Appendix C, we derive higher-order-in-e approximations to r z,BT , and demonstrate how they improve the match (decrease the mismatch) to the BT model.

Parameter-space metrics
In this Section, we calculate the coherent and semicoherent parameter-space metric for the ELL1 model. Compared to the circular case, the parameter space has two extra dimensions.
Since the ELL1 model differs at first order in e from the circular model, the coherent metric also differs at first order. However, for the { f ,ḟ , n x , n y , Ω orb , x} metric components, the first order terms are of O(1/Ω orb T obs ) and can be neglected; the dominant difference is second order in e. Thus, the coherent metric components given in Eqs. (28) and (50) Note that the off-diagonal component g T0ω does not vanish. As described in the previous Section, this complicates the form of the mismatch. We now change to the parameters {T asc , 1 , 2 }, for which it is convenient to use Eq. (74). For these, the diagonal components are For small eccentricities e, the dominant metric components are given above. For completeness, we list the O(e 1 )corrections, which are all off-diagonal: The remaining off-diagonal components of the orbital metric are of O(1/Ω orb T obs ). These metric components have been found to be a good approximation even for higher eccentricities where the ELL1 model is not sufficient to track the rotational phase in a search and higher order models need to be used. This might be because many of the linear-in-e terms vanish from the metric.
The semi-coherent metric components are very similar to the coherent ones. The components associated with the non-eccentric parameters { f ,ḟ , n x , n y , Ω orb , x}, calculated in the circular case in Eqs. (42) and (55) remain valid; they have only second-order corrections in e. For the remaining orbital parameters {T asc , 1 , 2 } the semi-coherent metric components are the same as in the coherent case (this follows from Eq. (53)). Thus, the diagonal components for {T asc , 1 , 2 } areḡ TascTasc =2π 2 f 2 x 2 Ω 2 orb , where we omit terms of O(1/Ω orb T obs ). Thus, the semicoherent metric for the ELL1 model simply adds the components above to the semi-coherent metric for the circular model. In Figure 3, the mismatch and its coherent-metric approximation are compared for small parameter offsets, for a realistic simulated pulsar. Apart from different f -andḟ -scales, the corresponding plot for the semi-coherent mismatch looks very similar. The mismatch agrees well with its metric approximation for mismatch m ≤ 0.5, which is typical: in Appendix B, we show that the highest sensitivity at given computing cost for an elliptical search is obtained with an average mismatchm = 0.471 (see Table 3).
The sky position parameters {n x , n y } are not shown in Figure 3 because we assume that for spider pulsars they are known to high precision from optical observations. A full blind search for binary pulsars in elliptic orbits is computationally impossible. There are too many parameterspace dimensions -even for circular orbits with reasonable parameter ranges the grid has too many points. To make a search possible, one needs tight constraints derived from optical/X-ray observations of the pulsar's companion star. In the next Section, we will discuss constraints and the search design for the probable eccentric orbit binary gammaray pulsar in 3FGL J0523.3−2528 (Strader et al. 2014;Saz Parkinson et al. 2016).

Search design for low eccentricity binary
In this Section, we discuss how to reduce the search parameter space using orbital constraints for the gamma-ray source 3FGL J0523.3−2528, presumed to be a pulsar in an eccentric binary orbit. This is similar to the circular example of Section 3.2.
The gamma-ray source itself was investigated by Saz Parkinson et al. (2016) and ranked 9th highest in a list of most significant 3FGL unassociated sources predicted to be pulsars. It shows typical pulsar-like properties: the photon flux is stable over time and the spectrum is fit by an exponential cutoff power law. The source is not in the Galactic disc, which increases the odds that it hosts an MSP. Earlier optical observations identified a likely companion and indicate an orbit with small, but not negligible, eccentricity of e = 0.04 (Strader et al. 2014). In contrast to the previous paragraph, this suggests that the pulsar is a YP, because binary MSPs tend to have rather circular orbits (Phinney 1992).
The sky-position search range of the probable pulsar within 3FGL J0523.3−2528 is tightly constrained from the X-ray and optical observations of the likely companion discussed above (Strader et al. 2014). At the time, the best estimate for the optical position was from the USNO-B1.0 Catalog (Monet et al. 2003). It is now also identified in the Gaia DR2 Catalog (Gaia Collaboration et al. 2018), whose pointing is so precise (see Table 2) that even at f = 1.5 kHz no search over sky position is required.
The orbital-parameter search ranges shown in Table 2 come from the Strader et al. (2014) analysis of the photometric and spectroscopic optical data. The orbital period and eccentricity parameters are constrained by the periodic opticalflux modulation. They assume that this arises from viewing a tidally-locked and deformed (ellipsoidal) companion at different aspect angles. Hence, the orbital period is twice the observed modulation period. (Another possible explanation for the modulation would be irradiation, but spectroscopic data does not show the orbital-phase-dependent temperature change that would be expected.) The orbital period is constrained to P orb = 0.688134 ± 0.000028 d at epoch of superior conjunction T 0.5 = 56577.14636 ± 0.0037 MJD. The eccentric parameters {e, ω} fall in the ranges e = 0.040 ± 0.006 and ω = 214 ± 10 deg.
The semimajor axis x is constrained using Eq. (57). This is similar to our previous example in Section 3.2, but requires fewer assumptions because the mass ratio q = M 2 /M 1 is directly bounded from the observations. To do this, Strader et al. (2014) estimate the rotational velocity of the companion's Roche lobe from high-quality optical spectra. Combined with the companion's radial velocity K 2 = 190.3 ± 1.1 km/s, this constrains the mass ratio q = 0.61 ± 0.06. Returning to Eq. (57), this gives x = 3.66 ± 0.38.
For our search, we also need the epoch of ascending node T asc . However, the results of Strader et al. (2014) are given in terms of the epoch of superior conjunction T 0.5 . For circular orbits, T 0.5 and T asc differ by P orb /4, but for eccentric orbits the relation is more complicated. To second order in e, it is For 3FGL J0523.3−2528 with e = 0.04, this O(e 2 ) approximation is more accurate than the uncertainties in the measured quantities on the rhs. (Higher-order approximations in e would be required for pulsars in binary orbits with larger eccentricities or longer orbital periods.) The resulting T asc is given in Table 2.
A search for a pulsar in an eccentric orbit is very similar to one for a pulsar in a circular orbit. The only differences are that a more general model for the Rømer delay is required to track the pulsar phase, and the orbital grids need to cover five orbital dimensions. While the latter is much more complex, it can be done with the same optimized stochastic search grid construction methods that are used in the circular case.
To accurately track the rotational phase of the pulsar, requires a higher-order-in-e approximation to r z,BT than the ELL1 model, unless the eccentricity is very small. Such approximations are computed in Appendix C. There, we also determine which order in e is sufficient.
For the case of 3FGL J0523.3−2528, a model of O(e 3 ) is sufficient. In Figure 6, we show that the rotational-phase error is negligible for the constrained parameter ranges given above.
Analogously to Eq. (58), the minimum number of grid points for the orbital parameter space can be computed from the proper 5-volume N orb ≈ m −5/2 detḡ dλ orb .
Here, the metric has the 5 dimensions {x, Ω orb , T asc , 1 , 2 }. This integral is proportional to  where x ∈ [x min , x max ]. ∆Ω orb , ∆T asc , ∆ 1 , and ∆ 2 are the search ranges for the corresponding parameters, and we made the assumption that ∆Ω orb Ω orb . The number of orbital grid points and subsequently the computing cost depend even stronger on f and x in an eccentric search than in a circular one.
The computing cost of a search for 3FGL J0523.3−2528 is estimated based on the number of grid points. We assume search ranges in f andḟ as given earlier in this Section. The remaining parameter-space ranges are given in Table 2. The required total computing cost of the search is estimated by multiplying the cost of one FFT with the number ofḟ -grid points, and the f -dependent number of orbital grid points, and then summing over the f intervals.
To exemplify the computing cost of a search for 3FGL J0523.3−2528, we express it in terms of search duration on Einstein@Home, assuming that the project provides 25,000 GPU-hours per week. This is shown in Figure 4 as a function of the maximum searched frequency. For comparison, we also show the search duration for a circular binary search, i.e. setting e = 0 and not searching over { 1 , 2 }. An eccentric MSP search up 1.5 kHz would take more than 100 million years on Einstein@Home, and even a YP search would take more than 100 years. Circular searches for YPs or MSPs up to 400 Hz would take a few 100 days. Note that the search ranges are still the 1σ ranges, so searches within the 3σ range would be more computing intensive.
In summary, this Section has shown how computing intensive a search for 3FGL J0523.3−2528 would be. An eccentric MSP search even to low frequencies ∼ 100 Hz is not feasible with the current constraints, and a YP search would be very expensive. In the optical data, Strader et al. (2014) do not see evidence for a "false" eccentricity, but a circular search would be much less computing intensive than a eccentric one. With slightly tighter constraints searches up to 800 Hz could be feasible.

COMPARISON WITH OTHER METHODS
Similar and alternative methods are used to search for binary pulsars in data from radio telescopes and gravitationalwave detectors. In this Section, we will review these, compare them to the methods presented here, and discuss their applicability to searches for binary gamma-ray pulsars.
In addition to coming from diverse messengers and frequencies, the data have other key differences. The gammaray data is similar to the gravitational-wave data: the length of the data sets is months to years and the instruments simultaneously detect signals from a substantial fraction of the sky. In contrast, typical radio surveys collect data in stretches of minutes from tiny fractions of the sky. While gamma-ray data consist of discrete photon arrival times, radio and gravitational wave data are continuous. Therefore, it is not surprising that some pulsation search methods might work for one kind of data but not for the other.
For these other data sources, many methods have been employed by many individuals and groups. Here we are guided by reviews from Lorimer & Kramer (2004) for radio search methods and Messenger et al. (2015) for gravitational-wave methods. We exclude methods which require data from two detectors.

Acceleration searches
Time-domain "acceleration searches" have been very successful in finding new radio pulsars in binaries with orbital periods shorter than a day (see, e.g., Camilo et al. 2000). Fourier-domain acceleration searches have also been successfully used to discover binary radio pulsars (see, e.g., Ransom et al. 2001;Andersen & Ransom 2018). A similar approach to search for continuous gravitational waves is the "polynomial search" (van der Putten et al. 2010).
These searches do not use a model which describes periodic orbital motion. Instead, they assume constant acceleration along a straight line (see also Johnston & Kulkarni 1991). This accurately describes an orbiting system only if the data set is much shorter than one orbital period. Since the LAT data set is more than a decade long, acceleration searches would only find binary gamma-ray pulsars whose orbital periods were decades or longer.
It is straightforward to quantify the range of orbital periods an acceleration search is sensitive to. Assume that the data set is less than ∼ 10% of the orbital period and is near the superior or inferior conjunction, where the velocity is changing linearly with time (Johnston & Kulkarni 1991). An acceleration a along the line of sight ("los") towards Earth contributes an amountḟ to the observed spin-frequency derivative. The maximum acceleration at inferior or superior conjunction is for a circular orbit a = cx Ω 2 orb , and for an eccentric orbits a = cx Ω 2 orb (1 + e)/(1 − e). Therefore, searches would be sensitive if the sum of the intrinsic pulsar spin-down and this line-of-sight contribution to the spin-down is within the search range. Since the intrinsic spin-down is usually negative, this is most likely if the acceleration towards Earth is positive, i.e. if the pulsar is near the superior conjunction.
Current blind search surveys for isolated gamma-ray pulsars are a form of acceleration search because they scan over spin-down (Clark et al. 2017). For YPs they search down toḟ = −10 −9 Hz/s and for MSPs down toḟ = −10 −13 Hz/s. In principle, these searches are sensitive to pulsars like the young ( f ≈ 7 Hz) binary pulsar PSR J2032+4127, which is in a 45 − 50-yr orbit around its companion (Ho et al. 2017). It was found in an isolated gamma-ray search (Abdo et al. 2009c), and only afterwards it was discovered to be in a binary system . The orbit is highly eccentric (e ≈ 0.93 − 0.99) with x ≈ 7000 − 20000 s. The maximum spin-down contribution should therefore be of order | max{ḟ los }| = 10 −10 Hz/s. This is in the search range if the pulsar is near superior conjunction during the mission time.
Blind searches which assume linear acceleration, i.e. they search over constantḟ , are only sensitive to binary pulsars with P orb 10T obs . To become sensitive to shorter orbital periods, higher-order frequency derivatives must be searched. "Jerk" searches, which include the second-order frequency derivativef , improve the sensitivity for pulsars with orbital periods in the range P orb ∈ [7T obs , 20T obs ] and have been successfully used in a radio pulsar search (Andersen & Ransom 2018). Alternatively, the full orbital motion may be taken into account, as in Allen et al. (2013).

Stack/slide search
The "stack/slide" method has been used in radio pulsar searches like the Parkes Multibeam Pulsar survey to account for binary motion (Faulkner et al. 2004). This led to the discovery of the double neutron star system PSR J1756−2251 with an orbital period of 7.7 hr (Faulkner et al. 2005). (The words "stack/slide" are used in continuous gravitationalwave searches, not to account for binary pulsar motion but rather to remove the effects of the Earth rotation and motion around the SSB (Brady & Creighton 2000;Riles 2017). That is also the case for the semi-coherent searches we describe in this paper to account for the LAT's motion around the SSB.) In a stack/slide search the data set is broken into subsets of length T coh , corresponding to frequency bins of width ∆ f = 1/T coh . T coh is chosen to be small enough that the Doppler modulation induced by motion of the detector around the SSB, or of the pulsar around the binary center of mass, remains within a single bin. For circular binary motion, provided that T coh is a factor of a few smaller than P orb , this im-plies f xΩ 2 orb T coh < 1/T coh .
Each of these subsets is then Fourier transformed. The resulting power spectra are added (stacked) together after the Doppler modulation is compensated by shifting the frequency (slide) in each of the spectra; sources give rise to peaks in the stacked spectra. This technique is only sensitive if the subsets are much shorter than the Doppler modulation period. This technique is useless for spider gamma-ray pulsars because detection statistics are constructed from the differences of photon arrival times. Spider pulsars have typical orbital periods of P orb 1 d, so data subsets would have to be shorter than a few hours. Most data subsets would contain no photons. A few would contain one photon. Almost none would contain enough photons to compute the differences of photon arrival times.
Stack/slide could be used for gamma-ray pulsars in orbits where P orb is too small for an acceleration search, but is much larger than the T coh ≈ 24 d used in this paper. Using Kepler's third law, the condition of Eq. (84) can be written where M 1 is the pulsar mass and q = M 2 /M 1 is the mass ratio.
(In fact this applies provided T coh P orb .) This shows that with our choice of T coh , stack/slide methods might be able to find gamma-ray pulsars with planetary companions, with orbital periods longer than ∼ 1 yr and masses up to O(10) Earth masses.

Power spectrum search
The basic assumption of a "power spectrum search" is that the data set can be broken into subsets short enough that the observed spin frequency is constant in each one. This is the same assumption as in a stack/slide search. That technique is based on visual inspection and has been used to discover binary radio pulsars (see, e.g., Lyne et al. 2000).
To carry out the search, power spectra are computed for each subset. The spectra are binned in frequency and plotted with a frequency-versus-time color map. The colors show the power and make it easy to visually identify peaks in the power spectrum. A binary pulsar signal appears as a peak whose frequency varies sinusoidally with time.
The method "TwoSpect" uses a similar method to perform all-sky searches for continuous gravitational waves from sources in binary systems. The visual inspection is replaced by a second Fourier transform (hence the name TwoSpect; Goetz & Riles 2011). While no continuous gravitational waves have been detected, this technique has been used to put upper limits on continuous gravitational-wave emission from the low-mass X-ray binary Scorpius X-1 (Aasi et al. 2014).
The power spectrum search is not suitable for detecting gamma-ray spider pulsars for the same reasons as the stack/slide method.

Sideband search
"Sideband searches" have found many binary radio pulsars within globular clusters (Lorimer & Kramer 2004). The method has also been adapted to search for continuous gravitational waves from sources in binary systems (Messenger & Woan 2007;Sammut et al. 2014). One first carries out a search for isolated systems, as if there were no binary motion, and then looks for a characteristic structure in the results of that isolated search.
If a binary is present, orbital motion produces sidebands around a central peak at the spin frequency of the pulsar (Ransom et al. 2003). Since the isolated search does not remove the effects of the binary motion, a pulsar's power is spread over many Fourier bins (also called sidebands). This reduces the sensitivity compared to a matched-filter search.
The method is particularly useful for tight-orbit-binary pulsars where the orbital period is much smaller than the observation time span, which is the case of interest for spider pulsars. After detecting a signal, the binary parameters can be inferred from the locations and magnitudes of the sidebands and the central peak.
To see how this works, we compute the S/N of the coherent detection statistic P n for an isolated-pulsar template, with parameters {ν,ḟ , n x , n y , 0, 0, 0}, arising from a circular-binary pulsar with parameters { f ,ḟ , n x , n y , x, f orb , T asc }, where f orb = Ω orb /2π. This S/N is given by Eq. (17), which depends on the rotational phase difference due to the parameter mismatch: One can think of ν as denoting the pulsar frequency in the isolated search. Our derivation closely follows Ransom et al. (2003).
To compute the detection statistic P n , we evaluate Eq. (17) with the phase mismatch (86). We first re-express e in∆Φ using the Jacobi-Anger expansion with z = 2πn f x and ϑ = 2π f orb (t − T asc ), where J m is a Bessel function of the first kind. Multiplying this by e i2πn(ν− f )(t−tref) gives where, without loss of generality, we have set t ref = T asc . Since the S/N only depends on the modulus of e in∆Φ , we may also set T asc = 0. We assume that there is a large number of photons from the hypothetical pulsar, which have equal weights and arrive at uniformly spaced intervals in time. The double sum j =k in Eq. (17) may then be replaced by an integral over time, since On the rhs we have included the diagonal j = k term which is absent on the lhs, but is negligible in the limit where the number of photons N is large. The integral over time is with F = n(ν − f ) + m f orb . For observation times that include many orbits, the rhs of Eq. (90) is unity for F = 0 and is negligible otherwise. Thus, the only terms in Eq. (89) that survive are those for which ν = f − m f orb /n. When that is satisfied, we have where m is constrained by F = 0. Thus, the double sum in Eq. (89) vanishes at all frequencies ν except for the "sideband frequencies" ν = ν m = f − m f orb /n, where m takes on all integer values. We now evaluate the S/N θ 2 Pn (ν) from Eq. (17) by substituting in Eq. (91), assuming that the weights w j are constant. For the reasons just given, θ 2 Pn (ν) vanishes except at the discrete sideband frequencies ν m = f − m f orb /n. We obtain The quantity θ 2 Pn that appears on the rhs is given by Eq. (11). It is the S/N that the pulsar would have in an isolated search if the binary motion were absent. It is also the S/N that the pulsar would have in a binary-pulsar search at the true signal parameter values.
The structure in frequency space ν is evident from Eq.(92). As described by Ransom et al. (2003), the S/N is spread over equally-spaced sidebands around the pulsar frequency f , whose spacing is commensurate with the orbital frequency. The sideband width is ∼ 1/T obs , as can be seen from Eq. (90).
In comparison with a binary pulsar search, the isolatedpulsar search has lost some S/N, since J 2 m ≤ 1. To recover some of the lost S/N within the isolated-pulsar search, we introduce a new test statistic that sums over the first m orb sidebands around the central pulsar frequency. This cumulative sideband power may be written as with the detection statistic P n (ν) appropriate to an isolatedpulsar search with parameters {ν,ḟ , n x , n y , 0, 0, 0}. (A test statistic weighing the m'th sideband in Eq. (93) by J 2 m (2πn f x) would be more sensitive, but for simplicity it is not considered here.) The S/N for the cumulative sideband power B n is easily calculated. It is defined as where p is the pulsed fraction defined in Eq.
(3). The numerator of Eq. (94) can be found from Eq. (92), which implies that . Summing this over m gives the numerator: The denominator of Eq. (94) is defined in the absence of a signal, with p = 0. It is easily calculated if the noise at the different frequencies which contribute to the sum are independent. Since Poisson noise is stationary, these contributing terms will be independent if they are spaced more than one frequency bin apart, where the bins have width 1/nT obs . Since the sideband frequencies are separated by f orb /n, these different terms will be independent if there are many orbits in the observation time: f orb T obs 1. Each term in the denominator then has variance 4, so the sum yields E 0 [B 2 n ] − E 2 0 [B n ] = 4(2m orb + 1). Thus, the S/N for B n is To maximize this S/N, what is the optimal number of sidebands m orb to include? As shown by Ransom et al. (2003), the optimal number of sidebands to include depends upon where square brackets denotes "integer part". To see this, consider the sum which appears in Eq. (96): For m orb < M orb this sum grows (approximately linearly) with increasing m orb . But the addition theorem for Bessel functions ensures that Eq. (98) stops growing and approaches unity as soon as m orb exceeds M orb . Since the denominator of the S/N in Eq. (96) has a term that grows like √ 2m orb + 1, the S/N is maximized for m orb = M orb . For this number of sidebands, one thus obtains for the expected S/N of the cumulative sideband power. The behavior we have just described, considered alongside the definition (94) of the S/N, shows the main weakness of sideband searches. The numerator grows (approximately) linearly as we include more sidebands, meaning that we can recover all of the signal power. But, in the absence of a signal, B n undergoes a random walk as sidebands are included, and so the denominator of Eq. (94) (the root-mean-squared of B n in the absence of a signal) increases as √ 2M orb + 1. Thus, in comparison with an optimal matched filter, the incoherent summation over sidebands loses a factor of √ 2M orb + 1 in the S/N. This is explicit in Eq. (99), and makes sideband searches ineffective if there are many sidebands, as is often the case. For example, consider the potential circular-binary pulsar in 3FGL J1653.6−0158 and the potential eccentricbinary pulsar in 3FGL J0523.3−2528 discussed earlier in this paper. Their estimated parameter ranges in f and x give rise to large numbers of sidebands.
This means that sideband searches work best if only a few sidebands are expected, meaning that 2πx f , the total rotational phase arising from the orbital modulation, is small. This is the case for black-widow systems, which have very light companions. The small companion mass means that the pulsar orbits very close to the center of mass, so the projected semimajor axis x is extremely small. Note that the modulation can be small even for the high frequencies f typically found for black widows. Figure 5 illustrates this, for example, for the black-widow pulsar PSR J1311−3430, which would have been a candidate for a sideband search. The figure shows the expected optimal matched-filter S/N θ 2 P1 required to exceed a threshold in the expected cumulative sideband S/N θ 2 B1 > 100, which is a reasonable threshold for confident detection. From Eq. (99), this requires θ 2 P1 to exceed 100 √ 2M orb + 1. Hence, M orb is constant on the contour lines, which therefore denote boundaries of constant f x. Since the largest observed θ 2 P1 values for known pulsars are ∼ 1000, the region below and to the left to the contour line corresponding to θ 2 P1 = 1000 might be considered for sideband searches.
Sideband searches within gamma-ray binaries like LS 5039 and LS I +61 303 would also be justified. These systems contain a compact object: a black hole or neutron star. Comparison between the expected cumulative sideband S/N θ 2 B 1 and the expected optimal matched-filter S/N θ 2 P 1 . For given frequency f and semimajor axis x, the black contours show the θ 2 P 1 required to exceed a threshold θ 2 B 1 > 100. The crosses are at the locations of two known pulsars: PSR J1311−3430 and PSR J2339−0533. The red lines show four potential sideband search candidates. For the YP candidates LS 5039 and LS I+61 303, and the spider candidate 3FGL J0523.3−2528, the approximate values for the semimajor axes are known. The dashed line shows the maximum semimajor axis value for 3FGL J1653.6−0158.
Since both binaries are highly eccentric (0.3 < e < 0.6; Aragona et al. 2009), the compact objects could be YPs. These two candidate pulsars are both displayed in Figure 5. This is purely illustrative, since the sideband power B n defined here is only suitable for circular-binary pulsars. Eccentric pulsars will have additional sidebands (Ransom et al. 2003), and thus must have an even higher pulsed fraction to be detectable in a sideband search.
This Section has not discussed the implementation of a practical sideband search. We would need some constraints on the parameters f orb and x to hunt for the sidebands. If those are available from optical observations, then the sky position will be known precisely. This in turn would make a fully-coherent isolated pulsar search computationally feasible. The resulting test statistics could than be used to construct the sideband search statistic B n of Eq. (93).

Discussion
The methods discussed in this Section have little applicability to searches for gamma-ray pulsars in spider systems, which are the main focus of this paper. But they are of interest for other types of binary systems.
Acceleration searches could discover binary pulsars with orbital periods comparable to, or longer than our observation time T obs ∼ 10 yr. These binaries have pulsars whose companions are very-low-mass stars or planets, in wide or-bits. These pulsars might have been missed by isolated pulsar searches.
Stack/slide and power spectrum methods do not appear suitable for spider gamma-ray pulsar searches. They might potentially detect systems with orbital periods longer than our typical coherence time T coh ∼ 24 d and shorter than T obs ∼ 10 yr. However, these searches are very expensive computationally.
Sideband searches could be used to hunt for binary pulsars with low spin frequencies or in very close orbits. While these are computationally less expensive than the search methods discussed earlier in this paper, they are also considerably less sensitive.
All of these methods have a domain of applicability. Given prior knowledge and constraints on a specific target, one can investigate these different methods to determine which are feasible and to estimate which one is potentially the most sensitive.

CONCLUSIONS
This work presents computationally efficient methods to detect circular-and eccentric-orbit binary gamma-ray pulsars in blind searches. These generalize techniques that have been previously developed to search for isolated pulsars (Pletsch & Clark 2014).
We have presented all of the elements of this generalization. Physically, the central element is a model that accurately describes the rotational phase of a pulsar over time as would be observed at the Solar System Barycenter. In comparison with the isolated model, this must also account for the Rømer delay caused by the binary motion. A second key element are semi-coherent and coherent test statistics, along with their expected signal-to-noise ratios. The last key element are the metrics for these statistics, which measure the "distance" in parameter space between two different rotational phase models. This metric quantifies the expected fractional loss in signal-to-noise ratio, and enables the construction of efficient parameter-space grids for a search.
We have shown how these different elements can be used together to search for gamma-ray pulsars. This is analogous to the isolated pulsar case (Pletsch & Clark 2014): the most computationally efficient approach is a multistage search with several semi-coherent and coherent stages. The computing cost is proportional to the number of points in the parameter-space grid. We compute this from the metric and show how the computing cost depends on the search parameters. This in turn allows the grid spacing to be optimized, achieving the highest possible sensitivity at fixed computing cost. These methods have been very successful in discovering isolated gamma-ray pulsars (Pletsch et al. 2012b;Clark et al. 2015Clark et al. , 2016Clark et al. , 2017Clark et al. , 2018. A truly blind search for binary pulsars is computationally impossible. Because the parameter space has at least 7 of the 9 possible dimensions { f ,ḟ , α, δ, x, P orb , T asc , 1 , 2 }, too many grid points are needed to cover it. However, in some cases, the number of dimensions can be reduced and/or the corresponding search ranges can be tightly constrained by multiwavelength observations. Such searches may be characterized as "vision impaired" rather than "blind". This paper considers two illustrative examples of this type, drawn from potential spider pulsars. Here, analysis of optical observations constrains the orbital parameters, and we show that searches of reasonable sensitivity (in some cases limited to young pulsars) are feasible. This enables "blind" searches for binary gamma-ray pulsars that were previously not feasible. This is important because these pulsars might be impossible to detect in other wavebands.
The methods of this paper, particularly the metric in parameter space, have applications beyond blind searches. There are binary pulsars which are visible in radio, optical, or X-ray, for which gamma-ray pulsations have not yet been found. For recent discoveries, precise determination of their orbital and other parameters is often not possible, since it requires observations spanning several years. The methods here are useful in such cases, to carry out efficient followup searches to discover gamma-ray pulsations. This way, within days or weeks after radio pulsations are discovered, the pulsar's parameters can be precisely measured over the > 10 yr of elapsed LAT mission time. This approach led to the discovery of gamma-ray pulsations soon after the radio detection of the 707 Hz black-widow pulsar PSR J0952−0607 (Bassa et al. 2017;Nieder et al. 2019).
A significant shortcoming of this paper's methods is that the number of grid points, and hence the required computing resources, grow quickly with increasing frequency f and semimajor axis x. To make searches feasible, it might be necessary to balance a reduced search range (smaller maximum f and/or x) versus a reduced search sensitivity (wider grid spacing and/or shorter coherence time). Even with large computing resources like Einstein@Home, millisecond pulsar searches for binaries with x ∼ seconds are only feasible if the orbital parameters are precisely constrained.
The second significant shortcoming is that search sensitivity is lost if the pulsar's rotational phase does not match our model. This can happen for several types of pulsars and binary systems. This paper assumes that the intrinsic spin frequency f varies linearly with time. It does not include the time-dependent variations or the unpredictable frequency glitches often seen in young pulsars. This means that pulsars could be "detected" in the semi-coherent stages of a search, but are then discarded after the coherent stage, because they did not match the phase model well enough to produce a significant detection statistic (see, e.g., Clark et al. 2017). Phase model mismatch can also arise from time-dependent variations of the orbital period P orb , which seems to be common in redback systems (see, e.g., Pletsch & Clark 2015). For pulsars in short-orbital-period binaries with heavy companions, post-Keplerian gravitational corrections also have to be taken into account (see, e.g., Damour & Deruelle 1986;Edwards et al. 2006).
Because of these limitations, this paper also evaluates alternative search methods, which have previously been used in radio and gravitational-wave searches. While these may be applied to search for binary gamma-ray pulsars, only the sideband search methods appear to have some chance to detect tight-orbit spider pulsars, which are the main focus of this paper.
A more detailed study is necessary to make a fair sensitivity comparison between the sideband search and this paper's methods. Indeed, while the cumulative sideband power loses a lot of signal-to-noise ratio compared this paper's methods, it might be improved. Since the sideband structure follows a known form, one could obtain a larger S/N by assigning weights to the sidebands before summing them, rather than using equal weights as done here.
We have implemented the new methods developed in this paper in a mixture of C and Python codes. These have been tested using simulated pulsar signals, both with our own code and with the widely used TEMPO2 package . We are confident that these codes work correctly, in part because they have discovered new spider pulsars, soon to be published.
We are currently using these codes and methods to hunt for spider pulsars in the unassociated sources of the Fermi LAT Fourth Source Catalog Catalog. These "blind" searches are guided by orbital constraints from optical observations. The orbital grids are constructed on the computing cluster AT-LAS at the Albert Einstein Institute in Hannover. The first two (semi-coherent) stages and the third (coherent) stage are all done on Einstein@Home, whose volunteers provide a massive computing pool. The final, less compute-demanding (H statistic) follow-up stage is done on ATLAS. To increase the computing power available in the initial stages of the search, we ported the search codes to work on Ein-stein@Home-volunteer's GPUs. The ATLAS cluster is also used to carry out follow-up gamma-ray searches of newly discovered radio pulsars, to refine the parameters as discussed above and in Nieder et al. (2019). This paper has used the two gamma-ray sources 3FGL J1653.6−0158 and 3FGL J0523.3−2528 as examples, to show how a realistic search might be structured. Both of these searches are being carried out, and the results will be discussed in upcoming papers. A similar search for a pulsar within 3FGL J2039.6−5618 (Romani 2015;Salvetti et al. 2015) successfully detected pulsations using the methods presented here (C.J. Clark et al. 2020, in prep.).
The outlook for future searches is promising. The Gaia Catalog provides sky locations for the spider companions, which are precise enough so that no search in {α, δ} is required. In addition, since the Large Area Telescope mission is ongoing, data sets are getting longer. Current searches use T obs ∼ 11 yr of data, compared with initial searches with T obs ∼ 4 yr. Furthermore, our available computing power is also increasing with time. This means that current searches employ T coh ∼ 24 d in the first stage, compared with initial searches with T coh ∼ 12 d. Since search sensitivity scales with (T coh T obs ) 1/4 (Pletsch & Clark 2014), our current sensitivity has increased by more than 50%. We believe that O(10 − 30) of the unassociated sources in the 4FGL Catalog are undiscovered spider pulsars, and that we can find some of them.

A. EXPECTATION VALUES OF SIGNAL STATISTICS
Here we show how to calculate the expectation values of signal statistics. The statistics depend upon the j = 1, . . . , N modeled pulsar rotation phases at the photon arrival times t j . To simplify the language and notation, we suppose that the vector of parameters λ = { f ,ḟ , α, δ} is fixed, and denote the modeled rotation phases by Φ j = Φ(t j , λ) = Φ(t psr (t j , α, δ), f ,ḟ ). Sums and products over j, k, run from 1, . . . , N unless otherwise specified. Finally, we write "the phase of the j'th photon", rather than "the modeled pulsar rotational phase associated with the j'th photon".
Our key assumption is that the phase of each photon is an independent (hence uncorrelated) random variable. This is justified because the number of photons detected is much less than one per pulsar revolution. The phase Φ j of the j'th photon is drawn from the distribution F j (Φ j ) as given in Eq. (3). Thus, using Eq. (4), the probability distribution function of Φ j is where the Fourier coefficients γ n are defined by Eq. (5) for n > 0, by γ n = γ * −n for n < 0, and by γ 0 = 0 for n = 0.
The expectation value of any quantity Q(Φ 1 , . . . , Φ N ) is now given by where the statistical independence of the rotation phases allows the probability density to be written as a product. For example the expected value of exp(−inΦ j ) is where δ nm is the Kronecker delta, giving unity for n = 0. The expected value of the coherent power signal statistic Eq. (6) is In the product above, only two terms are nontrivial, for which either = k or = j. The integrand does not depend upon the other N − 2 integration variables, whose corresponding integrals give unity, since the probability density is normalized.
One obtains On the first line, the first sum comes from terms with j = k and the second sum from terms where j = k, and we have used Eq. A3 to simplify both terms.

B. MAXIMAL SENSITIVITY AT FIXED COMPUTING COST
The sensitivity of a search can be quantified via the pulse fraction p defined in Eq. (3). More sensitive searches can detect sources with smaller values of p.
If infinite computing power were available, we would employ the fully coherent detection statistics H or P 1 , and the sensitivity of a search would only be limited by the data. To determine that ultimate sensitivity, consider the expected S/N θ 2 P1 given in Eq. (11). A point in parameter space where θ 2 P1 exceeded some threshold θ 2 threshold (established by the desired false alarm and false dismissal probabilities) would be counted as a detection. A reasonable detection threshold might be θ 2 threshold = 50, corresponding to pulse fraction sensitivity p 2 > θ 2 threshold /|γ 1 | 2 µT obs . For typical values of µT obs = 500 effective photons and |γ 1 | 2 = 0.8, this gives an ultimate, data-limited sensitivity of p 2 > 0.13.
In practice, with limited computing power, we adopt the multi-stage hierarchical approach described in Section 2.3. A sensible choice is to use most of the computing power in the first, semi-coherent stage. Roughly speaking, this is because a signal will only be found if it rises above the detection threshold in the first stage of the search 5 . Hence we will assume that our sensitivity is limited by the first semi-coherent search stage.
The maximum possible sensitivity of the semi-coherent stage is determined by the threshold on the semi-coherent S/N, whose expected value is given in Eq. (36). The threshold is lower than before, typically θ 2 S1 > θ 2 threshold = 10. Using search parameters from Eq. (27) and later in that Section gives a minimum detectable pulse fraction of p 2 > θ 2 threshold /|γ 1 | 2 µ √ T obs T coh = 0.31. As before, this is the theoretical sensitivity that could be achieved with unlimited computing power, but employing the semi-coherent statistic.
In practice, we must take the computing cost into account. This cost is proportional to the number of grid points 5 Of course this depends upon the choice of threshold and the region of parameter space around a candidate which is searched in the subsequent stages. If the full parameter space is searched for each candidate, then the statement is false! in parameter space at which the detection statistic is calculated. Reducing the number of grid points (corresponding to a larger average mismatch) loses some S/N but the additional computing power may be used to increase the coherence time T coh , which increases the S/N. What compromise maximizes the search sensitivity for a given computing cost? To find the optimal balance between the worst-case grid mismatch m and the coherent integration time T coh , we maximize the sensitivity with the constraint that the compute power is fixed, as described in Prix & Shaltev (2012) and Pletsch & Clark (2014). What is important is the rate at which the number of grid points grows with increasing T coh , which in turn depends upon the dimension of the parameter space.
The number of dimensions d in the search parameter space is determined by our prior knowledge. To quantify that, we use n orb (possible values 3 or 5) for the number of orbital parameters searched, and n sky (possible values 0 or 2) for the number of sky dimensions searched, so d = 2 + n orb + n sky . In the case of an eccentric binary with poorly known position, we have the full parameter space discussed in the main text, { f ,ḟ , n x , n y , Ω orb , x, T asc , 1 , 2 }, so n orb = 5, n sky = 2, and d = 9. For an eccentric binary whose position is precisely known (for example from optical observations), {n x , n y } are omitted from the search, n orb = 5, n sky = 0, and d = 7. For a circular binary whose position is precisely known, { 1 , 2 } are also omitted, so n orb = 3, n sky = 0, and d = 5.
The smallest detectable pulse fraction (averaged over signal location in parameter space) may be written as Here,m represents the average (over parameter space) mismatch of the grid (Prix & Shaltev 2012). The construction of our parameter space grid is described following Eq. (55); its average mismatch may be estimated as follows. Within a given 8 Hz frequency interval, the grid is the direct product of an equally-spaced grid in the frequency direction, an equally-spaced grid in theḟ direction, a 2-dimensional hexagonal lattice in sky position {n x , n y }, and an optimized stochastic grid in the orbital parameters. Below, we call these "subgrids". To determine the computing cost, we need to count the number of grid points in these subgrids, and multiply them together.
Because the metric has no off-diagonal terms that couple the different subgrids, the average parameter space mismatcĥ m can be written aŝ m =m f +m˙f +m sky +m orb , wherem f is the average mismatch in the frequency dimension (if all other parameters are exactly matched to the signal), andm˙f ,m sky ,m orb are the corresponding average mis-matches in theḟ , sky, and orbital subgrids (if all other parameters are exactly matched to the signal). The frequency dimension is searched with an FFT whose frequency spacing d f = 1/T coh . For the worst case, which is two points separated by d f /2, the quadratic metric approximation predicts a mismatchḡ f f /(2T coh ) 2 = π 2 /24 = 0.411, and hence an average mismatchm f = 0.14. As is often the case, the quadratic approximation slightly overestimates the mismatch; the Spherical Ansatz of Allen (2019) predicts a worst-case m = sin 2 ( π 2 /24) ≈ 0.36 which agrees well with the numerically measured value given in Section 5.2 of Pletsch & Clark (2014). In fact, as described before Eq. (42) of that paper, we can reduce the average mismatch tom f = 0.075 at almost no extra computational cost, by interpolating the frequency spectrum.
Theḟ subgrid has uniform spacing dḟ , and is an example of a regular lattice. For regular lattices, the average mismatchm is related to the worst-case mismatch m viâ m = ξm, where ξ ∈ [0, 1] is a lattice-dependent dimensionless geometrical factor called "thickness' (Prix & Shaltev 2012). Here we have a (one-dimensional) hypercubic grid, for which ξ = 1/3, so the average mismatchm˙f = m˙f /3, where m˙f =ḡ˙f˙f (dḟ /2) 2 = π 2 T 2 coh T 2 obs dḟ 2 /288 is the maximal mismatch in theḟ dimension. (Since the differences are small, for simplicity we do not employ the Spherical Ansatz further.) The sky subgrid is a hexagonal lattice with thickness ξ = 5/12 ≈ 0.416. Hence,m sky = 0.416m sky , where m sky is the worst-case sky mismatch.
The orbital parameter grid has an average mismatch which is well-estimated during the process of its construction, and can be easily controlled via the parameter that determines when new points are added to the stochastic bank.
The computing cost is the product of the number of grid points in the non-frequency dimensions with the cost of a single FFT. The number of grid points can be estimated using arguments like those given in deriving Eq. (58). In each of the different subgrids, the number of grid points is proportional tom −D/2 wherem is the average mismatch in that subgrid, and D is the dimension of that subgrid. Hence the number of grid points in theḟ subgrid is proportional to T cohm −1/2 f and the number of grid points in the sky subgrid is proportional to T nsky cohm −nsky/2 sky . The number of grid points in the orbital subgrid is proportional tom −norb/2 orb and is independent of T coh . Since the cost of an FFT is proportional to T coh log T coh , this gives a total computing cost C Here C 0 is a constant, and following Pletsch & Clark (2014), we have omitted the slowly-varying logarithmic factor from the cost of the FFT.
The method of Lagrange multipliers can be used to maximize sensitivity p −2 S1 at fixed computing cost 6 . The quantity we extremize is where λ is the Lagrange multiplier, s = 2 + n sky , and c 1 and c 2 are constants (independent of the average mismatches and T coh ). Extremizing L with respect to the coherence time and the three different average mismatches gives where we have made use of Eq. (B7) to evaluate the derivatives ofm.
To find the average mismatches that maximize the sensitivity at fixed computing cost, combine the first equation in turn with the second or third or fourth: T coh drops out and one obtains a closed form for the corresponding average mismatch. The independence from coherence time T coh in the binary pulsar case was previously shown for the isolated pulsar case by Pletsch & Clark (2014). For example, to solve form˙f , multiply the first equation by T 1/2 coh , multiply the second equation by 2sm˙f T −1/2 coh , and add them. One obtains (1 −m)/2 − 2sm˙f = 0, whose solution ism˙f = (1 −m)/4s. Doing this for all three combinations yieldŝ m˙f = 1 −m 4(2 + n sky ) , m orb = 1 −m 4(2 + n sky ) n orb , and (B10) m sky = 1 −m 4(2 + n sky ) n sky .
(B11) Table 3. A comparison of computationally unlimited and optimal computationally-limited semi-coherent searches, showing mismatches and sensitivity. The columns show the average template bank mismatcĥ m, and the average mismatches in theḟ , sky and orbital subgrids. (Note that the average per-dimension mismatch is constant.) Then the corresponding maximumḟ and sky mismatch are listed with the (square of the) minimum detectable pulsed fraction p. The first row shows the ideal semi-coherent case where the grid points are infinitesimally spaced and the compute cost is infinite. The next three rows illustrate smaller and smaller binary-system parameter spaces. The final row is for an isolated pulsar with unknown sky position. Thus we havem˙f = 1 −m f 9 + n orb + 5n sky , m orb = 1 −m f 9 + n orb + 5n sky n orb , and (B12) m sky = 1 −m f 9 + n orb + 5n sky n sky , which in turn allows us to determine the average and maximum mismatch in each of the subgrids, and the corresponding search sensitivity compared with an extremely finely spaced (but computationally very expensive) semi-coherent search.
In practice, after setting the mismatch as given by this optimal point, one adjusts the coherence time T coh to be as long as allowed by the available computing resources. What does this imply about the sensitivity? Above, we showed that with reasonable assumptions, a semi-coherent search can detect a pulsed fraction p 2 > 0.31 if there are infinite computing resources. With finite computing resources, this is increased by a factor of 1/(1 −m) = (9 + n orb + 5n sky )/4(2 + n sky )(1 −m f ), as can be seen from Eq. (B6). The corresponding loss of sensitivity is shown in Table 3. The achievable pulsed-fraction sensitivity is not far from the ideal case.
This analysis extends previous work (Pletsch & Clark 2014), which assumed a grid with fixed thickness ξ = 1/3 in all dimensions. However, this is not the case for current searches. Here, we have considered a grid which is a product of subgrids, each of which can have different geometrical properties, as used in existing searches. If we assume fixed thickness, then our results and in particular the final line of Table 3 agree with Eq. (H2) from Pletsch & Clark (2014).

C. HIGH ORDER PHASE MODEL FOR ELLIPTICAL BINARIES
The main text uses a linear-in-e "ELL1" approximation to the correct "BT" line-of-sight motion in eccentric orbits.
Here, we consider higher orders in the eccentricity e. The BT model is given in Eq. (64): We express this as where α n (e) and β n (e) are power series in e. The goal here is to find these functions, and to determine the appropriate order needed for our searches. (Taff (1985) gives an expansion of sin E and cos E in powers of e, but does not give a similar expansion for the line-of-sight motion.) For the derivation of the power series we introduce the Bessel functions and some of their properties. For positive integers n the Bessel function can be expressed as the power series The relation is also needed. Following Taff (1985), we start with the Fourier expansion of cos E: where we have introduced the generalized binomial coefficient r k = r · (r − 1) · · · (r − (k − 1)) k! .
The Cauchy product of √ 1 − e 2 andβ gives β n (e) = The α n follow directly fromα n , and differ only for n = 0. We list the results to 11th order. (A similar calculation (Dhurandhar & Vecchio 2001) gives the coefficients to seventh order, but without a general formula.) The α's are given by: β 10 = 78125 145152 e 9 ,  Figure 6. Mismatch between the BT model and models truncated at orders e 0 , e 1 , e 2 , and e 3 , for the source 3FGL J0523.3−2528 with e = 0.04. This is computed on a grid of 100 × 100 simulated pulsar signals, with equally spaced log 10 f /Hz ∈ [0, 3], and log 10 x/s ∈ [−2, 1]. The gray, dashed line indicates the semimajor axis x = 3.66 of the likely pulsar in 3FGL J0523.3−2528. The slopes of the constant-mismatch contours are the same for different models because e is fixed.
A line-of-sight model accurate to O(e k ) requires retaining terms up to and including α k+1 and β k+1 . Depending on the search parameters { f , x, e} different orders of these Taylor series will be required.
Consider the source 3FGL J0523.3−2528. The expected eccentricity is e ∼ 0.04. To find the appropriate order in e, we simulated 10,000 realizations of a pulsar in 3FGL J0523.3−2528, with different spin frequencies f and semimajor axes x. Figure 6 shows the mismatches which arise from using approximations of different orders in e, compared to the full BT model. For high frequencies the mismatch m is significant m ∼ 0.3 (S/N loss of up to 30%) for the O(e 2 )-model. A sensible choice is the O(e 3 )-model, for which the mismatch is below 1% for frequencies f < 1 kHz.
For systems with different eccentricities, we can also provide guidance. Since most of the known spider pulsars are MSPs, we simulated 10,000 realizations of a 1 kHz-pulsar with different semimajor axes x and eccentricities e. Figure 7 shows the mismatches that arise up to sixth orders in e.   Figure 7. Same as Figure 6, but varying the eccentricity e with fixed frequency f = 1 kHz, and going up to e 5 . The mismatch is computed on a grid of 100 × 100 simulated pulsar signals, with equally spaced log 10 e ∈ [−3, 0], and log 10 x/s ∈ [−2, 1].