Directed searches for continuous gravitational waves from binary systems: parameter-space metrics and optimal Scorpius X-1 sensitivity

Directed searches for continuous gravitational waves from binary systems: parameter-space metrics and optimal Scorpius X-1 sensitivity

Paola Leaci Max-Planck-Institut für Gravitationsphysik, Albert-Einstein-Institut, D-14476 Golm, Germany Dip. di Fisica, Università di Roma “Sapienza”, P.le A. Moro, 2, I-00185 Rome, Italy    Reinhard Prix Max-Planck-Institut für Gravitationsphysik, Albert-Einstein-Institut, D-30167 Hannover, Germany
July 18, 2019

We derive simple analytic expressions for the (coherent and semi-coherent) phase metrics of continuous-wave sources in low-eccentricity binary systems, for the two regimes of long and short segments compared to the orbital period. The resulting expressions correct and extend previous results found in the literature. We present results of extensive Monte-Carlo studies comparing metric mismatch predictions against the measured loss of detection statistic for binary parameter offsets. The agreement is generally found to be within . As an application of the metric template expressions, we estimate the optimal achievable sensitivity of an Einstein@Home directed search for Scorpius X-1, under the assumption of sufficiently small spin wandering. We find that such a search, using data from the upcoming advanced detectors, would be able to beat the torque-balance level Wagoner (1984); Bildsten (1998) up to a frequency of , if orbital eccentricity is well-constrained, and up to a frequency of for more conservative assumptions about the uncertainty on orbital eccentricity.

04.80.Nn, 95.55.Ym, 95.75.-z, 97.60.Gb, 07.05.Kf

I Introduction

Continuous gravitational waves (CWs) are a promising class of signals for the second-generation detectors currently under construction: advanced LIGO (aLIGO) Harry (2010), advanced Virgo The Virgo Collaboration (2009) and KAGRASomiya (2012). These signals would be emitted by spinning neutron stars (NSs) subject to non-axisymmetric deformations, such as quadrupolar deformations (“mountains”), unstable oscillation modes (e.g. r-modes) or free precession. For a general review of CW sources and search methods, see for example Prix (2009).

A particularly interesting type of potential CW sources are NSs in low-mass X-ray binaries (LMXBs), with Scorpius X-1 being its most prominent representative Watts et al. (2008). The accretion in these systems would be expected to have spun up the NSs to the maximal rotation rate of  G.B. Cook, S.L. Shapiro, and S.A. Teukolsky (1994a). All the observations performed to date, however, show they only spin at several hundred , so there seems to be something limiting the accretion-induced spin-up. A limiting mechanism that has been suggested is the emission of gravitational waves Wagoner (1984); Bildsten (1998); Ushomirsky et al. (2000), which would result in steady-state CW emission where the accretion torque is balanced by the radiated angular momentum. For a discussion of alternative explanations, see  Bildsten (1998); Ushomirsky et al. (2000). The resulting torque-balance CW amplitude increases with the observed X-ray flux (independent of the distance of the system, e.g. see Eq.(4) in Bildsten (1998)), therefore the X-ray brightest LMXB Scorpius X-1 is often considered the most promising CW source within this category Watts et al. (2008); Bulten et al. (2014).

Several searches for CW signals from Scorpius X-1 have been performed (without any detections) on data from initial LIGO Abbott et al. (2007); Aasi et al. (2014a, b), and several new pipelines have been developed and are currently being tested in a Scorpius X-1 Mock Data Challenge (MDC) Bulten et al. (2014). So far, all current methods fall into either one of two extreme cases: highly coherent with a short total time baseline (6 hours in Abbott et al. (2007), 10 days in Aasi et al. (2014b)), or highly incoherent with a long total baseline but very short (hours) coherent segments Aasi et al. (2014a); Chung et al. (2011); Aasi et al. (2014c).

In order to increase sensitivity beyond these methods, and to be able to effectively absorb large amounts of computing power (such as those provided by Einstein@Home Ein (); Aasi et al. (2013)), it is necessary to extend the search approach into the realm of general long-segment semi-coherent methods, by stacking coherent -statistic segments Brady and Creighton (2000); Krishnan et al. (2004); Prix and Shaltev (2012). Such methods have already been employed over the past years for both directed and all-sky searches for CWs from isolated NSs Pletsch and Allen (2009); Aasi et al. (2013, 2013), but they have not yet been extended to the search for CWs from binary systems.

Key ingredients required for building a semi-coherent “StackSlide” search are the coherent and semi-coherent parameter-space metrics Owen (1996); Balasubramanian et al. (1996); Brady et al. (1998); Brady and Creighton (2000). The coherent binary CW metric was first analyzed in Dhurandhar and Vecchio (2001), and this was further developed and extended to the semi-coherent case in Messenger (2011). Here we will largely follow the approach of Messenger (2011), considering only low-eccentricity orbits, and focusing on two different regimes of either long coherent segments, or short segments compared to the orbital period, both of which admit simple analytic results.

Plan of this paper

Sec. II provides a general introduction of the concepts and notation of semi-coherent StackSlide methods, parameter-space metric and template banks. Secs. IIIV build on the work of Messenger (2011), rectifying some of the results and extending them to the case of a general orbital reference time. In Sec. VI we present the results of extensive Monte-Carlo software-injection tests comparing the metric predictions against the measured loss in statistics due to parameter-space offsets. Sec. VII applies the theoretical template-bank counts to compute optimal StackSlide setups for a directed search for Scorpius X-1 and estimates the resulting achievable sensitivity. Finally, Sec. VIII presents concluding remarks.

Notation. Throughout the paper we denote a quantity as when referring to the coherent case, and as when referring to the semi-coherent case.

Ii Background

The strain in a given detector due to a continuous gravitational wave is a scalar function , where is the time at the detector. The set of signal parameters denotes the four amplitude parameters, namely the overall amplitude , polarization angles and , and the initial phase . The set of phase-evolution parameters consist of all the remaining signal parameters affecting the time-evolution of the CW phase at the detector, notably the signal frequency , its frequency derivatives (also known as “spindown” terms), sky-position of the source, and binary orbital parameters. We will look in more detail at the phase model in Sec. III.

The observed strain in the detector affected by additive Gaussian noise can be written as , and the detection problem consists in distinguishing the (pure) “noise hypothesis” of from the “signal hypothesis” .

ii.1 Coherent detection statistic

As shown in Jaranowski et al. (1998), matched-filtering the data against a template , and analytically maximizing over the unknown amplitude parameters , results in the coherent statistic , which only depends on the data and on the template phase-evolution parameters . This statistic follows a (non-central) distribution with degrees of freedom, and a non-centrality parameter , which depends on the signal parameters and the template parameters . The quantity (which depends linearly on ) is generally referred to as the signal-to-noise ratio (SNR) of the coherent -statistic. The expectation value of this statistic over noise realizations is given by Jaranowski et al. (1998)


In the case of unknown signal parameters within some parameter space , the number of required template that need to be searched is typically impractically large. The maximal achievable sensitivity by this coherent statistic is therefore severely limited Brady et al. (1998) by the required computing cost. It was shown that semi-coherent statistics generally result in better sensitivity at equal computing cost (e.g. see Brady and Creighton (2000); Prix and Shaltev (2012)).

ii.2 Semi-coherent detection statistic

Here we focus on one particular semi-coherent approach, sometimes referred to as StackSlide, which consists of dividing the total amount of data into segments of duration , such that (in the ideal case of gapless data). The coherent statistic is then computed over all segments , and combined incoherently by summing111This form of the “ideal” StackSlide statistic is not directly used for actual searches. For reasons of computational cost, the coherent statistics across segments are computed on coarser template banks, and are combined on a fine template bank by summing across segments using interpolation on the coarse grids. This is discussed in detail in Prix and Shaltev (2012), but is not relevant for the present investigation.


This statistic follows a (non-central) -distribution with degrees of freedom and a non-centrality parameter , which is found to be given by


in terms of the per-segment coherent SNRs . Note that (contrary to the coherent case) can not be regarded as a semi-coherent SNR. The expectation for is


ii.3 Template banks and metric mismatch

In order to systematically search a parameter space using a statistic (which here can refer either to the coherent or the semi-coherent ), we need to select a finite sampling of the parameter space, commonly referred to as a template bank, and compute the statistic over this set of templates, i.e. .

A signal with parameters will generally not fall on an exact template, and we therefore need to characterize the loss of detection statistic as a function of the offset from a signal. This is generally quantified using the expected statistic of Eqs. (1) and (4), namely (by removing the bias ) as the relative loss in the non-centrality :


which defines the measured -statistic mismatch function. This mismatch has a global minimum of at and is bounded within .

We now use two standard approximations: (i) Taylor-expand this up to second order in small offsets , and (ii) neglect the dependence on the amplitude parameters (the effect of which was analyzed in detail in (Prix, 2007a)), which leads to the well-known phase metric , namely


with a positive-definite symmetric matrix, and implicit summation over repeated indices is used. The quality of this approximation will be quantified in a Monte-Carlo study in Sec. VI. The metric mismatch has a global minimum of at and (contrary to ) is semi-unbounded, i.e. . Previous studies testing against Prix (2007a); Wette and Prix (2013); Wette (2014) for all-sky searches have shown that the phase-metric approximation Eq. (6) works reasonably well for observation times and for mismatch values up to . Above this mismatch regime, the metric mismatch increasingly over-estimates the actual loss . For example, in Prix (2007a) an empirical fit for this behavior was given as .

The metric formalism was first introduced in Owen (1996); Balasubramanian et al. (1996) (in the context of compact binary coalescence signals), and was later applied to the CW problem in Brady et al. (1998), where it was shown that the coherent phase metric can be expressed explicitly as


where , with the CW signal phase, and denotes the time average over the coherence time , i.e. .

The corresponding semi-coherent phase metric was first studied in Brady and Creighton (2000) and was found to be expressible as the average over all per-segment coherent phase metrics222This assumes constant signal SNR over all segments , which should be a good approximation for stationary noise and sufficiently long segments, such that diurnal antenna-pattern variations have averaged out.) , i.e.


where we defined the average operator over segments as


The metric is useful for constructing template banks, by providing a simple criterion for how “close” we need to place templates in order to limit the maximal (relative) loss in detection statistic, which can be written as


This states that the worst-case mismatch to the closest template over the whole template bank should be bounded by a maximal value . Note that each template “covers” a parameter-space volume given by


which describes an -dimensional ellipsoid, where is the number of parameter-space dimensions. The template bank can therefore be thought of as a covering of the whole parameter space with such template ellipsoids, such that no region of remains uncovered. One can show Brady et al. (1998); Prix (2007b); Messenger et al. (2009) that the resulting number of templates in such a template bank is expressible as


where is the number of template-bank dimensions, and is the center density (also known as normalized thickness) of the covering lattice, which quantifies the number of lattice points per (metric) volume for a unit mismatch. The center density is a geometric property of the lattice, and for the typical covering lattices used here is given by Prix (2007b)


An important caveat for using Eq. (13) is that the metric determinant and integration must only include fully resolved template-bank dimensions, which require more than one template to cover the extent of the parameter space along that dimension. We can estimate the per-dimension template extents at given maximal mismatch as333Note that the extra factor of here compared to Watts et al. (2008) is required to account for the total extent of the ellipse, not just the extent from the center. This can also be seen by considering that the center density in Eq. (13) for the 1-dimensional case.


where are the elements of the inverse matrix to the metric . We can therefore define the metric template-bank “thickness” along dimension in terms of the corresponding effective number of templates along that direction as


For dimensions with we must exclude this coordinate from the bulk template-counting formula Eq. (13), as it would effectively contribute “fractional templates” and thereby incorrectly reduce the total template count.

Iii The binary CW phase

iii.1 The general phase model

The general CW phase model assumes a slowly spinning-down NS with rotation rate and a quadrupolar deformation resulting in the emission of CWs. The phase evolution can therefore be expressed as a Taylor series in the NS source frame as


where denotes the reference time and are the CW frequency and spindown parameters. These are defined as


where is a model-dependent constant relating the instantaneous CW frequency to the NS spin frequency . For example, for a quadrupolar deformation rotating rigidly with the NS (“mountain”) we have  Bildsten (1998); Ushomirsky et al. (2000); Cutler (2002); Melatos and Payne (2005); Owen (2005), while for r-modes  Bildsten (1998); Owen et al. (1998); Andersson et al. (1999) and for precession  Jones and Andersson (2002); Van Den Broeck (2005).

In order to relate the CW phase in the source frame to the phase in the detector frame needed for Eq. (8), we need to relate the wavefront detector arrival time to its source emission time , i.e. , such that . Neglecting relativistic wave-propagation effects444These effects are taken into account in the actual matched-filtering search codes, but are negligible for the calculation of the metric. (such as Einstein and Shapiro delays, e.g. see Hobbs et al. (2006); Edwards et al. (2006) for more details), we can write this as


where is the vector from solar-system barycenter (SSB) to the detector, is the speed of light, is the (generally unknown) distance between the SSB and the binary barycenter (BB), is the radial distance of the CW-emitting NS from the BB along the line of sight, where means the NS is further away from us than the BB, and is the unit vector pointing from the SSB to the source. In standard equatorial coordinates with right ascension and declination , the components of the unit vector are given by .

The projected radial distance along the line of sight can be expressed (e.g. see Roy (2005)) as


where is the distance of the NS from the BB, is the inclination angle between the orbital plane and the sky, is the argument of periapse, and is the true anomaly (i.e. the angle from the periapse to the current NS location around the BB). We further approximate the orbital motion as a pure Keplerian ellipse, which can be described as


in terms of the semi-major axis and the eccentricity . The ellipse can be written equivalently in terms of the eccentric anomaly , namely


and the dynamics is described by Kepler’s equation, namely


which provides a (transcendental) relation for . Combining Eqs. (20), (21) and (22), we can rewrite the projected radial distance in terms of , namely


where we defined . Combining this with Eq. (23) fully determines (albeit only implicitly) the functional relation required for the timing model of Eq. (19).

Dropping the unknown distance to the BB (which is equivalent to re-defining the intrinsic spindown parameters), and defining the SSB wavefront arrival time as


we can rewrite the timing relation Eq. (19) as


Here we are interested only in binary systems with known sky-position , and we can therefore change integration variables from to in the phase for the metric integration of Eq. (8). Furthermore, given that the Earth’s Rømer delay is bounded by s, and , this difference will be negligible for the metric and can be dropped. This is equivalent to effectively placing us into the SSB, which is always possible for known . In order to simplify the notation, we now simply write .

Plugging this timing model into the phase of Eq. (17), we obtain555The sign on (and equivalently on ) here agrees with Eq.(2.26) in Blandford and Teukolsky (1976), but differs from Eq.(2) in Messenger (2011), which is incorrect.


with . The binary systems we are interested in have semi-major axis of order s, and binary periods of order of several hours. Hence, the change in , and therefore , during the time will be negligible, and so we can approximate , namely


For the purpose of calculating the metric using Eq. (8), we can further approximate the phase model in the following standard way (e.g. see Jaranowski et al. (1998); Wette and Prix (2013)) as


where we expanded the factors and kept only the leading-order terms in , keeping in mind that and , where denotes the coherence time. In the last term we replaced the instantaneous intrinsic CW frequency as a function of by a constant parameter , namely we approximated . This “frequency scale” of the signal could be chosen as the average (or the largest) intrinsic CW frequency of this phase model over the coherence time . Given that this only enters as a scale parameter in the metric, and the changes of intrinsic frequency over the observation time will generally be small, this introduces only a negligible difference.

iii.2 The small-eccentricity approximation ()

We follow the approach of Messenger (2011) and consider only the small-eccentricity limit of the metric. Hence, by Taylor expanding Eqs. (24) and (28) up to leading order in , i.e. inserting into Kepler’s equation Eq. (28), we obtain


where is the mean orbital angular velocity. Plugging this into Eq. (24), we obtain the Rømer delay of the binary to leading order in as


where a constant term was omitted, which is irrelevant for the metric. We use the standard Laplace-Lagrange parameters defined as


and the mean orbital phase


measured from the time of ascending node , which (for small ) is related to by Hobbs et al. (2006)


and which (contrary to the time of periapse ), remains well-defined even in the limit of circular orbits.

The small-eccentricity phase model is therefore parametrized by the 5 binary parameters (referred to as the “ELL1” model in Hobbs et al. (2006); Edwards et al. (2006)), and in the circular case () this reduces to the 3 binary parameters with .

Iv The binary CW metric

iv.1 Phase derivatives

Following Dhurandhar and Vecchio (2001); Messenger (2011) we restrict our investigation to (approximately) constant-frequency CW signals, i.e. . This is motivated by the assumed steady-state torque-balance situation in LMXBs, which are our main target of interest. However, the corresponding fluctuations in the accretion rate are expected to cause some stochastic frequency drift, and one will therefore need to be careful to restrict the maximal coherence time in order to limit the frequency resolution. This will be discussed in more detail in the application to Scorpius X-1 in Sec. VII. The total phase-evolution parameter space considered here is therefore spanned by the following 6 coordinates


The small-eccentricity phase model used here can now be explicitly written as


with given by Eq. (35). The frequency parameter is treated as a constant in the phase derivative (which corresponds to neglecting the small correction ), but numerically we have in the present constant-frequency case without spindowns. Hence, we obtain the phase derivatives


which are inserted into Eq. (8) in order to obtain the coherent metric. The semi-coherent metric is obtained by averaging the coherent metrics over segments according to Eq. (9).

The resulting analytic expressions for these metrics in the general case are quite uninstructive and unwieldy, while it is straightforward to compute them numerically for any case of interest. However, as noticed in previous investigations Dhurandhar and Vecchio (2001); Messenger (2011), it is instructive to focus on two limiting regimes that yield particularly simple analytical results, namely the long-segment limit () where , and the short-segment limit () where .

When taking these limits on the metric , in order to decide whether a particular off-diagonal term is negligible or not, it is useful to consider the diagonal-rescaled metric


which is dimensionless and has unit diagonal, . In this rescaled metric we can then naturally neglect off-diagonal terms if they satisfy , as their corresponding contribution to the metric mismatch of Eq. (12) will then be negligible.

iv.2 The long-segment (LS) regime ()

iv.2.1 Coherent metric

As noted in Messenger (2011), it is convenient to use the discrete gauge freedom in the choice of the orbital reference epoch , as we can choose the time of ascending node during any orbit. Therefore we can add any integer multiple of a period and redefine


without changing the system. An (infinitesimal) offset (or “uncertainty”) changes under such a transformation in the presence of an (infinitesimal) offset on the period, i.e. if , namely


All other coordinate offsets are unaffected by this change of orbital reference epoch. This needs to be carefully taken into account in the metric when redefining . Namely, given that this is a pure “relabeling” of the same physical situation, the corresponding metric mismatch must be invariant, i.e.


where the only nonzero components of the Jacobian of this coordinate transformation Eq. (41) are


In the long-segment limit, these discrete steps are assumed to be small compared to the segment length . We can therefore choose the time of ascending node to be approximately at the segment midtime , i.e. we consider the special gauge choice , in which the phase derivatives Eq. (39) are time-symmetric around , which simplifies the metric calculation. Introducing the time offset


we therefor consider first the special case .

Keeping only leading-order terms in , and up to first order in (i.e. first order in and ), we obtain the diagonal metric


in agreement with the result Eq. (30) of Messenger (2011).

The general case of can be obtained by applying the Jacobian666A direct but somewhat more complicated calculation of the metric at yields the same result of Eq. (44) with , and so and , resulting in the general coherent long-segment metric with components


while all other off-diagonal terms are (approximately) zero in this limit. The nonzero off-diagonal term shows that in general (i.e. for ) there are correlations between offsets in and , introduced by the relation Eq. (42). This generalizes the result in Messenger (2011) to arbitrary choices of orbital reference epoch .

iv.2.2 Semi-coherent metric

In order to compute the semi-coherent metric, we only need to average the per-segment coherent metrics of the previous sections over segments according to Eq. (9). Here it is crucial to realize that all segment metrics must use the same coordinates in order for this averaging expression to apply, in particular we can only fix the gauge on once for all segments, and so generally . This means we cannot use the special form of of Eq. (46) for the average over all segments, as done in Messenger (2011), which would result in the erroneous conclusion that the semi-coherent metric would be identical to the coherent one. Instead, using Eq. (9), we find by averaging the per-segment metrics :


It is interesting to consider the ideal case of regularly-spaced segments without gaps, i.e. , and


resulting in the first two moments


where is the average mid-time over all segments. In this case we can write the variance of as


If we use the gauge freedom and choose , then , and , so that in this case the metric is again diagonal, and the only component different from the coherent metric is


corresponding to a refinement in the coordinate. Therefore the semi-coherent template-bank spacing in needs to be finer by a factor of compared to coherent spacing in order to achieve the same mismatch. The existence of this metric refinement in in the long-segment limit had already been noticed earlier Pletsch (2014), and was also used in the discovery of a binary pulsar in Fermi LAT data Pletsch et al. (2012).

iv.3 The short-segment (SS) regime ()

iv.3.1 Coherent metric

In the case of very short coherent segments compared to the period, i.e. , as pointed out in Messenger (2011), strong parameter-space degeneracies render the metric in these coordinates quite impractical for template-bank generation, and we therefore follow their approach of Taylor-expanding the phase around the mid-point of the coherent observation time, namely


omitting a constant term.

The -coordinates are defined as the k-th time-derivatives of the phase at , namely


Note that the phase model of Eq. (54) is formally identical to the phase in terms of the frequency and spindowns of an isolated NS (with playing the role of the reference time ), as can be seen in Eq. (29). The corresponding “spindown metric” has a well-known analytical form, which can be expressed most conveniently in terms of the rescaled dimensionless -coordinates,


resulting in the -coordinate metric


correcting the incorrect expression Eq. (21) of Messenger (2011).

Applying Eq. (55) to the small-eccentricity phase of Eq. (38), we obtain the following expressions


with , which differ in the sign of compared to Eqs. (A2-A5) of Messenger (2011), due to the sign error on in the phase model discussed in Sec. III.1.

Note that for the general elliptical case we expect independent -coordinates and in the circular case. The nonlinear coordinate transformation Eqs. (58) can be analytically inverted to obtain , as shown explicitly in Appendix A.

We follow Messenger (2011) in estimating the range of validity of this Taylor approximation up to order by considering the order of magnitude of the first neglected term of order , namely


where we defined the fraction of an orbit sweeped during the segment duration, i.e. . In order to link this phase error to a mismatch, we observe that here of Eq. (56), and the corresponding metric element of Eq. (IV.3.1) gives for and therefore


Plugging in typically values used in the numerical tests later, e.g. Hz, s, , we see that the mismatch grows very rapidly from at to at . We therefore expect this approximation with to typically hold for segment durations up to a to of an orbit.

iv.3.2 Semi-coherent metric

While the segments are assumed to be short compared to the period , we only consider the case where the total observation time is long compared to i.e. in addition to we also assume , which also implies .

We cannot directly use the coherent per-segment metrics in -coordinates to compute the semi-coherent metric as an average using Eq. (9), because a fixed physical parameter-space point would have different -coordinates in each segment.

Instead, as first shown in Messenger (2011), one can go back to physical coordinates and use the fact that to replace the discrete sum in Eq. (9) over segments by an integral over time, namely


in terms of coherent per-segment metrics expressed as a function of the segment mid-times , and where is the mid-time of the whole observation. In order to compute this, we first calculate the coherent metric of Eq. (8) for each segment by using the phase derivatives of Eq. (39). As before, we keep only first-order terms in (i.e. first-order in ), then Taylor-expand the results in up to second order, obtaining . This is integrated over the total observation time according to Eq. (IV.3.2), and only leading-order terms in are kept. We only keep off-diagonal terms where the corresponding diagonal-rescaled elements Eq. (40) are not . The resulting semi-coherent short-segment long-observation limit metric elements are found as


which generalizes the result in Messenger (2011) (for gauge choice , i.e. ) to the general case of .

Note that a more general form encompassing the semi-coherent metric in both the short-segment and long-segment limits has recently found in Whelan et al. (2015).

V Number of templates

In order to express the explicit template counts based on Eq. (13), we will assume a simple parameter space bounded by , , , , , and