An Analytic Approximation to the Bayesian Detection Statistic for Continuous Gravitational Waves
AAn Analytic Approximationto the Bayesian Detection Statisticfor Continuous Gravitational Waves
John J. Bero
Center for Computational Relativity and Gravitation and School of Physics andAstronomy, Rochester Institute of Technology, 84 Lomb Memorial Drive, Rochester,NY 14623, USA
John T. Whelan
E-mail: [email protected]
Center for Computational Relativity and Gravitation and School of MathematicalSciences, Rochester Institute of Technology, 85 Lomb Memorial Drive, Rochester, NY14623, USA
Abstract.
We consider the Bayesian detection statistic for a targeted search forcontinuous gravitational waves, known as the B -statistic. This is a Bayes factorbetween signal and noise hypotheses, produced by marginalizing over the fouramplitude parameters of the signal. We show that by Taylor-expanding to first orderin certain averaged combinations of antenna patterns (elements of the parameter spacemetric), the marginalization integral can be performed analytically, producing a closed-form approximation in terms of confluent hypergeometric functions. We demonstrateusing Monte Carlo simulations that this approximation is as powerful as the full B -statistic, and outperforms the traditional maximum-likelihood F -statistic, for severalobserving scenarios which involve an average over sidereal times. We also show thatthe approximation does not perform well for a near-instantaneous observation, sothe approximation is suited to long-time continuous wave observations rather thantransient modelled signals such as compact binary inspiral. a r X i v : . [ g r- q c ] A ug nalytic B -stat Approximation
1. Introduction
The signal from a non-precessing source of gravitational waves (GWs) such as a rotatingneutron star or slowly-evolving binary system, can be described by phase-modulationparameters, which determine the shape of the signal, and amplitude parameters. In thecase where the phase-modulation parameters are assumed to be known, the likelihoodratio between models with and without signal is a function of the four amplitudeparameters. Jaranowski Kr´olak and Schutz [1] constructed a maximum-likelihoodstatistic (known as the F -statistic), which is the basis of many existing searches forcontinuous GWs. Prix and Krishnan [2] proposed a Bayesian alternative (the B -statistic)which instead marginalized the likelihood ratio over these parameteres, assuming ageometrically-inspired prior distribution. Exact evaluation of the B -statistic requiresintegration over the four-dimensional amplitude parameter space; Whelan et al [3]showed that two of the integrals can be done analytically, but a two-dimensionalnumerical integration remains. They also showed that the marginalization integrals canbe done exactly if the parameter-space metric (determined by averaged combinationsof antenna patterns) has a block-diagonal form. In this paper, we generalize this resultto produce an analytical approximation to the B -statistic by Taylor expanding to firstorder in the off-diagonal metric elements.This paper is laid out as follows: in section 2 we give a brief overview of thebackground information and formalism related to this topic, including a discussion ofGW signal analysis and a description of the two detection statistics which already exist.Section 3 contains the derivation of our approximation and in section 4 we test the powerof the approximation as a detection statistic. Section 5 concludes with a summary ofthe results and their practical implications.
2. Formalism
We follow the conventions and notation of [3], where more details can be found. Wesummarize the relevant expressions here. For a GW signal coming from a sky positionspecified by right ascension α and declination δ , we can define a propagation unit vector (cid:126)k pointing from the source to the solar-system barycenter (SSB). The tensor GW canthen be resolved in a basis of traceless tensors transverse to (cid:126)k as h ↔ ( τ ) = h + ( τ ) e ↔ + + h × ( τ ) e ↔× . (2.1)For a nearly periodic signal, as from a rotating neutron star (NS), the polarizationcomponents are h + ( τ ) ≡ h χ ) cos[ φ ( τ ) + φ ] and h × ( τ ) ≡ h χ sin[ φ ( τ ) + φ ] , (2.2) nalytic B -stat Approximation χ = cos ι is the cosine of angle between the line of sight and the neutron star’srotation axis, and h = 4 Gc | I xx − I yy | Ω d (2.3)is the amplitude in terms of the equatorial quadrupole moments { I xx , I yy } , the rotationfrequency Ω, and the distance d to the source. The preferred polarization basis tensorsare given by e ↔ + = ε ↔ + cos 2 ψ + ε ↔× sin 2 ψ and e ↔× = − ε ↔ + sin 2 ψ + ε ↔× cos 2 ψ . (2.4)where ε ↔ + = (cid:126)ı ⊗ (cid:126)ı − (cid:126) ⊗ (cid:126) and ε ↔× = (cid:126)ı ⊗ (cid:126) + (cid:126) ⊗ (cid:126)ı are the fiducial basis tensors definedusing unit vectors orthogonal to (cid:126)k , with (cid:126)ı pointing “West on the sky” in the directionof decreasing right ascension α , and (cid:126) pointing “North on the sky” in the direction ofincreasing declination δ . The polarization angle ψ measures the angle counter-clockwiseon the sky from (cid:126)ı to the NS’s equatorial plane.The phase evolution φ ( τ ) in terms of the arrival time τ at the SSB can be writtenin terms of NS rotation or spindown parameters, e.g., φ ( τ ) = 2 π (cid:18) f τ + f τ · · · (cid:19) , (2.5)although it may be more complicated, e.g., for NSs in binary systems.The strain, h , measured by an interferometric GW detector whose arms are parallelto the unit vectors (cid:126)p and (cid:126)p is given by h = h ↔ : d ↔ (2.6)where ‡ d ↔ = (cid:126)p ⊗ (cid:126)p − (cid:126)p ⊗ (cid:126)p (cid:126)a ⊗ (cid:126)b ) : ( (cid:126)c ⊗ (cid:126)d ) =( (cid:126)a · (cid:126)d )( (cid:126)b · (cid:126)c ). The GW strain can also be expressed as h = h + F + + h × F × , (2.8)where F + and F × are the detector antenna pattern functions which depend on the 3angles defining the source sky position and polarization basis relative to your detector,which in our case would be the right ascension α , the declination δ and the polarizationangle ψ . If we separate out their dependence on ψ , then the pattern functions have theform F + ( α, δ, ψ ) = a ( α, δ ) cos 2 ψ + b ( α, δ ) sin 2 ψ (2.9a) F × ( α, δ, ψ ) = − a ( α, δ ) sin 2 ψ + b ( α, δ ) cos 2 ψ , (2.9b) ‡ We limit attention in this section to the long-wavelength limit, where the detectors are assumed to besmall compared to the gravitational wavelength c/f , which is appropriate to most observations withground-based interferometric detectors. At higher frequencies, the detector tensor d ↔ ( f ) is frequency-dependent and complex. See e.g.,[3] for more details. nalytic B -stat Approximation a and b are amplitude modulation coefficients defined in terms of the detectortensor d ↔ as a ≡ ε ↔ + : d ↔ , (2.10a) b ≡ ε ↔× : d ↔ . (2.10b)These coefficients are defined with respect to the reference polarization basis anddepend both on the sky position of the GW source and the sidereal time at whichthe measurement is taking place.It is useful to divide the signal parameters into amplitude parameters { h , χ, ψ, φ } and phase-evolution parameters such as the sky position { α, δ } , and any parametersdescribing φ ( τ ). The dependence of the signal on the amplitude parameters can bewritten simply as[1, 3] h ↔ ( τ ; A , λ ) = A ˘ µ h ↔ ˘ µ ( τ ; λ ) , (2.11)where the Einstein summation convention implies the sum (cid:80) µ =1 over repeated indices.The amplitudes {A ˘ µ } are defined by §A ˘1 = A r cos φ r and A ˘2 = A r sin φ r (2.12a) A ˘3 = A l cos φ l and A ˘4 = A l sin φ l ; (2.12b)where A r = h (cid:18) χ (cid:19) and φ r = φ + 2 ψ ; (2.13a) A l = h (cid:18) − χ (cid:19) and φ l = φ − ψ (2.13b)are the amplitudes and phases of the right- and left-circularly-polarized components ofthe signal, respectively. If we denote the data recorded in the GW detector(s) as x , and assume those data toconsist of the signal A ˘ µ h ˘ µ plus Gaussian noise, the sampling distribution for the datawill be pdf( x |A ) ∝ exp (cid:18) − (cid:0) x − A ˘ µ h ˘ µ | x − A ˘ µ h ˘ µ (cid:1)(cid:19) (2.14)The log-likelihood ratio will thus beΛ( {A ˘ µ } ; x ) = ln pdf( x |A )pdf( x |
0) = A ˘ µ x ˘ µ − A ˘ µ M ˘ µ ˘ ν A ˘ ν (2.15) § Our coordinates {A ˘ µ } , introduced in [3], are related to the more familiar Jaranowski-Kr´olak-Schutz(JKS) coordinates {A µ } of [1] by A = A ˘1 + A ˘3 , A = A ˘2 − A ˘4 , A = −A ˘2 − A ˘4 , A = A ˘1 − A ˘3 . nalytic B -stat Approximation x ˘ µ ≡ ( x | h ˘ µ ) is the scalar product (see Appendix A) of the data with the templatewaveform, and {M ˘ µ ˘ ν } ≡ { ( h ˘ µ | h ˘ ν ) } = I L − K I K LL K J − K L J (2.16)forms a metric on parameter space.If we define {M ˘ µ ˘ ν } as the matrix inverse of {M ˘ µ ˘ ν } , we can write the maximum-likelihood values of the amplitude parameters {A ˘ µ } as (cid:98) A ˘ µ ( x ) = M ˘ µ ˘ ν x ˘ ν , (2.17)Since the maximum-likelihood parameters { (cid:98) A ˘ µ ( x ) } contain equivalent information tothe projections { x ˘ ν } (which form jointly sufficient statistics for the amplitude parameters A ), we can use { (cid:98) A ˘ µ } as a representation of the relevant part of the data. Their samplingdistribution can be written as the multivariate Gaussianpdf( (cid:98) A|A ) = (det 2 π M ) − / exp (cid:18) −
12 ( (cid:98) A ˘ µ − A ˘ µ ) M ˘ µ ˘ ν ( (cid:98) A ˘ ν − A ˘ ν ) (cid:19) (2.18)This is useful for conducting Monte Carlo simulations (as was done in [2]): one neednot simulate the full GW data, only generate draws of the four maximum-likelihoodparameters { (cid:98) A ˘ ν } representing the data.It is also convenient to write the log-likelihood ratio in terms of (cid:98) A as well:Λ( A ; (cid:98) A ) = A ˘ µ M ˘ µ ˘ ν (cid:98) A ˘ ν − A ˘ µ M ˘ µ ˘ ν A ˘ ν . (2.19)This is written explicitly in terms of the polar representation in Appendix B.The F -statistic[1] is defined as the maximized log-likelihood ratio, F = max A Λ( A ; x ) = Λ( (cid:98) A ; x ) = 12 (cid:98) A ˘ µ M ˘ µ ˘ ν (cid:98) A ˘ ν = 12 I (cid:98) A r + 12 J (cid:98) A l + (cid:98) A r (cid:98) A l (cid:104) K sin( (cid:98) φ r − (cid:98) φ l ) + L cos( (cid:98) φ r − (cid:98) φ l ) (cid:105) (2.20)The B -statistic[2] is defined as the Bayes factor between models with and without signal: B ( x ) = pdf( x |H s )pdf( x |H n ) = (cid:82) pdf( x |A ) pdf( A|H s ) d A pdf( x |
0) = (cid:90) e Λ( {A ˘ µ } ; x ) pdf( A|H s ) d A (2.21)The prior is taken to be uniform in χ ∈ ( − , ψ ∈ ( − π/ .π/
4) and φ ∈ (0 , π ), sothat pdf( h , χ, ψ, φ |H s ) = pdf( h |H s )2 π (2.22)The convention introduced in [2] is to use an improper prior pdf( h |H s ) = A , 0 < h < ∞ ,so that B ( x ) = A π (cid:90) π (cid:90) − (cid:90) π/ − π/ (cid:90) ∞ e Λ( {A ˘ µ } ; x ) dh dψ dχ dφ = A π (cid:90) π (cid:90) π (cid:90) ∞ (cid:90) ∞ e Λ( {A ˘ µ } ; x ) dA r dA l dφ r dφ l √ A r A l (2.23) nalytic B -stat Approximation
3. An approximate form for the B -statistic Previous work [3] showed that the B -statistic integral (2.23) can be exactly evaluated inthe case where K = 0 = L , so that the metric (2.16) becomes diagonal and the left- andright-circularly polarized subspaces decouple. We show in Appendix A that K and L can be small compared to I = J , especially in continuous-wave observations containingan average over sidereal times and/or detectors.When K and L are small compared to I and J , it is fruitful to consider a Taylorexpansion of the B-statistic integral (2.23), which we carry out in Appendix B, and find B (0) = A [Γ( )] / ( IJ ) / (3.1)andln B ( x ) B (0) ≈ ln F (cid:32) , , I (cid:98) A r (cid:33) +ln F (cid:32) , , J (cid:98) A l (cid:33) + (cid:104) K sin( (cid:98) φ r − (cid:98) φ l ) + L cos( (cid:98) φ r − (cid:98) φ l ) (cid:105) (cid:98) A r (cid:98) A l × F (cid:16) , , I (cid:98) A r (cid:17) F (cid:16) , , I (cid:98) A r (cid:17) + 14 F (cid:16) , , J (cid:98) A l (cid:17) F (cid:16) , , J (cid:98) A l (cid:17) − F (cid:16) , , I (cid:98) A r (cid:17) F (cid:16) , , I (cid:98) A r (cid:17) F (cid:16) , , J (cid:98) A l (cid:17) F (cid:16) , , J (cid:98) A l (cid:17) (3.2)where the terms omitted are second order and higher in K and/or L .We can compare this to several limiting cases and alternative forms. First, notethat if K = 0 = L , we recover the result of section 6.1 of [3]. [See equation (6.11) ofthat work.] Second, in the limit that (cid:98) A r and (cid:98) A l are both large, the asymptotic form ofthe confluent hypergeometric functions [see identity (13.5.1) of [4]] F (cid:32) , , I (cid:98) A (cid:33) (cid:98) A →∞ −→ ) (cid:32) I (cid:98) A (cid:33) − / e I (cid:98) A / (3.3a) F (cid:32) , , I (cid:98) A (cid:33) (cid:98) A →∞ −→ ) (cid:32) I (cid:98) A (cid:33) − / e I (cid:98) A / (3.3b)says thatln B ( x ) B (0) (cid:98) A r , (cid:98) A l →∞ −→ I (cid:98) A r + 12 J (cid:98) A l + (cid:104) K sin( (cid:98) φ r − (cid:98) φ l ) + L cos( (cid:98) φ r − (cid:98) φ l ) (cid:105) (cid:98) A r (cid:98) A l −
34 ln (cid:18) I (cid:98) A r (cid:19) −
34 ln (cid:18) J (cid:98) A l (cid:19) − (cid:18) (cid:19) = F −
32 ln( (cid:98) A r (cid:98) A l ) + const (3.4)which is the result in equation (5.37) of [3]. nalytic B -stat Approximation
4. Evaluation of Approximation
We evaluate the approximation for three cases of interest, which are further detailed inAppendix A:(i) The case originally considered in [2]: a T obs = 25 hr observation of a source atright ascension 2 radians, declination − . I = J = 0 . T obs S n ( f ) , K = − . T obs S n ( f ) , and L = − . T obs S n ( f ) .so K/I = − . L/I = − . (cid:107) (ii) An observation with perfect sidereal-time averaging of a source on the celestialequator (declination 0) using only H1. As shown in Appendix A, this is a worst-case long-observation scenario, for which I = J = 0 . T obs S n ( f ) , K = 0, and L = − . T obs S n ( f ) . so K/I = 0 and
L/I = − . − . I = J = 0 . T obs S n ( f ) , K = 0 . T obs S n ( f ) , and L = 0 . T obs S n ( f ) .so K/I = 0 .
236 and
L/I = 0 . B -statistic Integral To compare our approximate form of the B -statistic to its exact value, we have toevaluate the integral (2.23). It was shown in [3] that the log-likelihood ratio (2.15) canbe written in physical coordinates asΛ( {A ˘ µ } ; x ) = h ω ( x ; χ, ψ ) cos( φ − ϕ ( x ; χ, ψ )) − h [ γ ( χ, ψ )] h and φ integrals performed explicitly to reduce the B -statistic to a doubleintegral ¶ B ( x ) = A √ π (cid:90) − (cid:90) π/ − π/ I ( ζ ( x ; χ, ψ )) e ζ ( x ; χ,ψ ) γ ( χ, ψ ) dψ dχ , (4.2)where ζ ( x ; χ, ψ ) = [ ω ( x ; χ, ψ )] γ ( χ, ψ )] . (4.3) (cid:107) Note that this case is slightly less favorable than another realistic alternative with the same skyposition, which averages over the O1 segments from LIGO Hanford and LIGO Livingston, for which I = J = 0 . T obs S n ( f ) , K = − . T obs S n ( f ) , and L = − . T obs S n ( f ) . so K/I = − . L/I = − . ¶ A similar reduction to a two-dimensional integral appears in [5], with the integrand empiricallyestimated rather than evaluated analytically. nalytic B -stat Approximation γ ( χ, ψ ) and ω ( x ; χ, ψ ). (The form of ϕ ( x ; χ, ψ ) isirrelevant to the result of the integral.) From (B.2) we can see γ ( χ, ψ ) = A ˘ µ M ˘ µ ˘ ν A ˘ ν h = 1 h (cid:2) IA r + J A l + 2 A r A l [ K sin( φ r − φ l ) + L cos( φ r − φ l )] (cid:3) = I (cid:18) χ (cid:19) + J (cid:18) − χ (cid:19) + 2 (cid:18) χ (cid:19) (cid:18) − χ (cid:19) [ K sin(4 ψ ) + L cos(4 ψ )] . (4.4)while ω ( x ; χ, ψ ) cos( φ − ϕ ( x ; χ, ψ )) = A ˘ µ x ˘ µ h = (cid:18) χ (cid:19) ( x ˘1 cos φ r + x ˘2 sin φ r ) + (cid:18) − χ (cid:19) ( x ˘3 cos φ l + x ˘4 sin φ l )= U cos φ + V sin φ (4.5)so [ ω ( x ; χ, ψ )] = U + V , where U = cos 2 ψ (cid:34)(cid:18) χ (cid:19) x ˘1 + (cid:18) − χ (cid:19) x ˘3 (cid:35) + sin 2 ψ (cid:34)(cid:18) χ (cid:19) x ˘2 − (cid:18) − χ (cid:19) x ˘4 (cid:35) (4.6a) V = cos 2 ψ (cid:34)(cid:18) χ (cid:19) x ˘2 + (cid:18) − χ (cid:19) x ˘4 (cid:35) − sin 2 ψ (cid:34)(cid:18) χ (cid:19) x ˘1 − (cid:18) − χ (cid:19) x ˘3 (cid:35) (4.6b)The simulations that follow, we evaluate the integrals for the B -statistic using a 3000-point Monte Carlo integration on the space χ ∈ ( − , ψ ∈ ( − π/ , π/ ψ , we still estimatethe χ integral accurately. We compare our approximation to the numerically-evaluated exact B -statistic, and tothe F -statistic. Each statistic is a function of the four data values { x ˘ µ } . However, if weexpress it in terms of the maximum-likelihood parameters { (cid:98) A ˘ µ } , we see that all of thestatistics are independent of the combination (cid:98) φ r + (cid:98) φ l = 2 (cid:98) φ and depend on the angles (cid:98) φ r and (cid:98) φ l only in the combination (cid:98) φ r − (cid:98) φ l = 4 (cid:98) ψ . Thus we can consider the statistics on thethree-dimensional space parameterized by (cid:98) A r ≥ (cid:98) A l ≥
0, and (cid:98) φ r − (cid:98) φ l ∈ [0 , π ). Forvisualization purposes, we plot contours of constant statistic versus (cid:98) A r and (cid:98) A l on slices ofconstant (cid:98) φ r − (cid:98) φ l , in analogy to Figure 3 of [3], which considered a metric with K = 0 = L ,for which the statistics were independent of (cid:98) φ r and (cid:98) φ l . If we plot (cid:98) φ r − (cid:98) φ l = 0 in thefirst quadrant and (cid:98) φ r − (cid:98) φ l = π in the second, we are effectively plotting (cid:98) A r cos( (cid:98) φ r − (cid:98) φ l )versus (cid:98) A l on the slice sin( (cid:98) φ r − (cid:98) φ l ) = 0. Likewise, if we plot (cid:98) φ r − (cid:98) φ l = π in the first nalytic B -stat Approximation (cid:98) φ r − (cid:98) φ l = − π in the second, we are effectively plotting (cid:98) A r sin( (cid:98) φ r − (cid:98) φ l )versus (cid:98) A l on the slice cos( (cid:98) φ r − (cid:98) φ l ). Since the approximate B -statistic and the F -statisticboth depend on the combination K sin( φ r − φ l ) + L cos( φ r − φ l ), the former slice focuseson the impact of L and the second on the impact of K . Note that another choice ofslice would be to chose φ r − φ l = tan − (cid:0) − LK (cid:1) , so that the K -and- L -dependent part ofthe statistics vanished, or φ r − φ l = tan − (cid:0) KL (cid:1) , which would maximize the impact ofthis term. In practice, for the examples we chose, | L | is significantly larger than | K | , sothese slices would be similar to the ones we plot.We choose our contours for these plots to correspond to specific false-alarmprobabilities (estimated by drawing 10 random points { (cid:98) A ˘ µ } from a Gaussian withzero mean and variance-covariance matrix {M ˘ µ ˘ ν } ) rather than specific statistic values.In figure 1, we see that for the case (i), with K/I = − . L/I = − . B -statistic contours are nearlyindistinguishable. Figure 2 shows case (ii), for which K/I = 0 . L/I = − . (cid:98) ψ ≈
0, i.e., (cid:98) A r e i (cid:98) φ r ≈ (cid:98) A l e i (cid:98) φ l . Finally, in figure 3we show the case (iii), with K/I = 0 .
236 and
L/I = 0 . To evaluate the performance of our B -statistic approximation, we produced MonteCarlo simulations by drawing 10 sets of signal parameters, using a fixed value of h = 10 S n ( f ) T obs and drawing the parameters χ , ψ , and φ from uniform distributions.Each of these sets of parameters was converted into a point A ˘ µ , and then a signal (cid:98) A ˘ µ was generated by drawing from a Gaussian with mean A ˘ µ and variance-covariancematrix {M ˘ µ ˘ ν } . A receiver-operating-characteristic (ROC) curve was generated for eachstatistic by plotting the fraction of signal points above a signal threshold (detectionprobability) against the fraction of noise points (described in the previous section) abovethe same threshold. The latter fraction is known as false-alarm probability, Type I errorprobability, or, in the language of hypothesis testing, significance. A superior detectionstatistic will have a higher detection efficiency at a given false-alarm probability, andthus be found above and to the left of an inferior one. Note that while the Neyman-Pearson lemma states that the Bayes factor will be the optimal test statistic for aMonte Carlo using the same prior[6] this is not guaranteed to be the case here, sincethe delta-function prior on h is not the same as the uniform prior used in defining thestatistic.In figure 4 we show the ROC curve for case (i), in which our approximation wasshown to match the exact B -statistic well (see figure 1). As expected, the approximate B -statistic performs as well as the exact one, and both outperform the F -statistic, asshown in [2]. In figure 5 we show the ROC curve for case (ii), where our approximationwas shown in figure 2 to have some discrepancies with the exact B -statistic. Nonetheless, nalytic B -stat Approximation b A r √ I b A l √ J . . . . e - e - e - e - e - b φ r − b φ l = 0 , π for case (i) B -statapprox F -stat5 4 3 2 1 0 1 2 3 4 5 b A r √ I b A l √ J . . . . e - e - e - b φ r − b φ l = ± π for case (i) Figure 1.
Comparison of B -statistic (2.23) and approximation (3.2), along with F -statistic (2.20); using the assumptions of a 25-hour observation beginning 2004 Jan1 at 00:00 UTC (GPS time 756950413) [case (i)], for which K/I = − . L/I = − . (cid:98) A r , (cid:98) A l , and (cid:98) φ r − (cid:98) φ l . Top: the slice sin( (cid:98) φ r − (cid:98) φ l ) = 0, for which the L -dependent terms of the statistics are important; bottom: the slice cos( (cid:98) φ r − (cid:98) φ l ) = 0, forwhich the K -dependent terms of the statistics are important. The contours of constantexact and approximate statistic are nearly indistinguishable, indicating that this is agood approximation for these metric values. nalytic B -stat Approximation b A r √ I b A l √ J . . . . . e - e - b φ r − b φ l = 0 , π for case (ii) B -statapprox F -stat5 4 3 2 1 0 1 2 3 4 5 b A r √ I b A l √ J . . . . e - e - e - b φ r − b φ l = ± π for case (ii) Figure 2.
Comparison of B -statistic and approximation, along with F -statistic,assuming a source on the celestial equator and H1 observations which evenly samplesidereal time [case (ii)], for which K/I = 0 . L/I = − . B -statistic contours at low false alarm rate in the case of linear polarization (cid:98) A r ≈ (cid:98) A l . Note that the disagreement for this contour in other directions is becauseit is drawn at the same false alarm probability, so the approximate B -statistic contourmust be inside the exact B -statistic contour to compensate for the deformation in onedirection. nalytic B -stat Approximation b A r √ I b A l √ J . . . . . . . . e - e - b φ r − b φ l = 0 , π for case (iii) B -statapprox F -stat5 4 3 2 1 0 1 2 3 4 5 b A r √ I b A l √ J . . . . e - e - e - b φ r − b φ l = ± π for case (iii) Figure 3.
Comparison of B -statistic and approximation, along with F -statistic,assuming a single Greenwich sidereal time of 00:00 [case (iii)], for which K/I = 0 . L/I = 0 . B -statistic are quite far off of those of the exact B -statistic. In fact,the entirety of both plots lie above the median of the approximate B -statistic underthe no-signal hypothesis; the contour in the upper left of the top plot is a false alarmprobability of .
05, and the one in the center of the lower plot is . (cid:98) A r = 0 = (cid:98) A l is at the 98th percentile of the approximate B -statistic, but the minimumof the exact B -statistic. Thus the approximation is, as expected, inappropriate for avalue of √ K + L /I so close to unity. nalytic B -stat Approximation − − − − False alarm probability0 . . . . . . D e t ec t i o np r o b a b ili t y ROC curve for case (i) B -statapprox F -stat Figure 4.
ROC curves for B -statistic and approximation, along with F -statistic,using the metric from case (i) (see figure 1). In this case, the approximate B -statisticperforms identically to the exact one. Compare figure 3 of [2]. we see that it again performs as well as the exact B -statistic and better than the F -statistic. In figure 6 we show the ROC curve for case (iii), where our approximation wasshown in figure 3 to disagree considerably with the exact B -statistic. Unsurprisingly, wefind this approximation to be a poor detection statistic in this scenario, underperformingboth the exact B -statistic and the F -statistic.
5. Conclusions
We have produced an analytic approximation to the B -statistic, a Bayesian detectionstatistic for continuous gravitational waves based on a Bayes factor between signaland noise hypotheses. This approximation is based on a Taylor expansion in theparameters K/I and
L/I , which are related to observation-averaged combinations ofantenna patterns, and depend on the sky position of the source, detectors involvedin the observation, and distribution of the observations in sidereal time. For long-time observations which average over a range of sidereal times, these parameters tendto be small enough to produce a good first-order approximation, and we showedvia Monte Carlo simulations that the approximate statistic performed as well as theexact B -statistic, even for a case with an expansion parameter approaching 50%. The nalytic B -stat Approximation − − − − False alarm probability0 . . . . . . D e t ec t i o np r o b a b ili t y ROC curve for case (ii) B -statapprox F -stat Figure 5.
ROC curves for B -statistic and approximation, along with F -statistic, usingthe metric from case (ii) (see figure 2). Even though K/I = 0 . L/I = − . B -statistic, which is Taylor expanded in K/I and
L/I , still performsas well as the exact B -statistic (and better than the F -statistic) in this Monte Carlo. approximation is shown to break down for observations at a single sidereal time, whichindicates the approximation is not likely to be an appropriate statistic for transientmodelled signals such as compact binary inspiral.Unlike the exact B -statistic, which must be evaluated via a two-dimensionalnumerical integral, our approximation (like the maximum-likelihood F -statistic) canbe evaluated analytically, which should make it computationally more efficient. + This,combined with the better detection efficiency than the F -statistic at the same falsealarm rate, makes it a potentially useful replacement for, or alternative to, the F -statistic in a semicoherent search which combines F -statistic values at a range of signalparameters. One potential challenge is that the approximate B statistic is expressedin terms of confluent hypergeometric functions, which may be more time-consuming toevaluate than the algebraic functions involved in the F -statistic. Additionally, directevaluation of these confluent hypergeometric functions for large-amplitude signals can + For example, for the python code used to perform the Monte Carlo simulations for this paper,the signel code used to calculate both the F -statistic and the approximate B -statistic together took O (1 − µ s) per evaluation, while the numerical integration for the exact B -statistic took O (1 ms) perevaluation. nalytic B -stat Approximation − − − − False alarm probability0 . . . . . . D e t ec t i o np r o b a b ili t y ROC curve for case (iii) B -statapprox F -stat Figure 6.
Comparison of B -statistic and approximation, along with F -statistic, usingthe metric from case (iii) (see figure 3). Here K/I = 0 .
236 and
L/I = 0 . B -statistic performs poorly, considerably below boththe exact B -statistic and the F -statistic. produce overflow, even though the final approximation in terms of their logarithmsand ratios may be well-behaved. It may be necessary to supplement standard libraryfunctions with strategic use of asymptotic forms. Acknowledgments
We wish to thank the members of the LIGO Scientific Collaboration and VirgoCollaboration continuous waves group for useful feedback. This paper has been assignedLIGO Document Number LIGO-P1800060-v5.
Appendix A. Form and Behavior of the Metric Elements
Given some nearly monochromatic GW signal around frequency f , the multi-detectorscalar product of two time series x and y , used in the definition (2.16), can be expressedas ( x | y ) ≡ (cid:88) Xl S Xl ( f ) Re (cid:90) ∞ (cid:101) x X ∗ l ( f ) (cid:101) y Xl ( f ) df , (A.1) nalytic B -stat Approximation S Xl ( f ) is the one-sided noise power spectral density around the frequency f in detector X during time stretch l , and (cid:101) x Xl ( f ) , (cid:101) y Xl ( f ) are the corresponding Fourier-transforms of x X ( t ) , y X ( t ) restricted to the time stretch l . This assumes the data fromeach detector X has been divided into short stretches of data [ t l , t l + T sft ) of length T sft .The metric components can be written as I = A + B + 2 E and J = A + B − E and K = 2 C and L = A − B (A.2)where, in the long-wavelength limit, A = (cid:88) Xl T sft S Xl ( f ) (cid:0) a Xl (cid:1) and B = (cid:88) Xl T sft S Xl ( f ) (cid:0) b Xl (cid:1) and C = (cid:88) Xl T sft S Xl ( f ) a Xl b Xl (A.3)and E = 0 (so that I = J ). As shown in [3], the more general expression, with a complexfrequency-dependent detector tensor d ↔ ( f ) and amplitude-modulation coefficients a ( f )and b ( f ), the (real) metric components can be more generally written as ∗ I = (cid:88) Xl T sft S Xl ( f ) (cid:12)(cid:12) a Xl ( f ) − ib Xl ( f ) (cid:12)(cid:12) and J = (cid:88) Xl T sft S Xl ( f ) (cid:12)(cid:12) a Xl ( f ) + ib Xl ( f ) (cid:12)(cid:12) (A.4a) L + iK = (cid:88) Xl T sft S Xl ( f ) (cid:2) a Xl ( f ) − ib Xl ( f ) (cid:3) ∗ (cid:2) a Xl ( f ) + ib Xl ( f ) (cid:3) (A.4b)In this form, we can see that the Cauchy-Schwarz inequality implies that K + L = | L + iK | ≤ IJ ; (A.5)in the long-wavelength case, this becomes √ K + L ≤ I = J .Prix and Krishnan [2] give an example of a T obs = 25 hr observation of a sourceat right ascension 2 radians, declination − . A = 0 . T obs S n ( f ) , B = 0 . T obs S n ( f ) , and C = − . T obs S n ( f ) , which is equivalent to I = J = 0 . T obs S n ( f ) , K = − . T obs S n ( f ) ,and L = − . T obs S n ( f ) . or K/I = − . L/I = − . K/I is small ( < .
10) everywhere, while the ratio
L/I is smaller away from the celestial equator.As an alternative to the arbitrarily chosen 25-hour observing time of [2], we canconsider the idealization that a long observation will include roughly the same amountof data from each sidereal time, and construct the corresponding metric components forthis case. Under this idealization, the metric components will be independent of rightascension, allowing us to simply plot them versus declination. In figure A3 we plot themetric elements and their ratios versus declination. We find, as in figure A1, the ratio ∗ Note that equation (A.3b) of [3] has the formulas for K and L reversed. nalytic B -stat Approximation K/I for Prix & Krishnan (25 hrs, H1) -0.516217 0.516217
L/I for Prix & Krishnan (25 hrs, H1) -0.516217 0.516217 √ K + L /I for Prix & Krishnan (25 hrs, H1) − . − . . . . . . . . . . C u m u l a t i v e D i s t r i bu t i o n Prix & Krishnan (25 hrs, H1)
K/IL/I √ K + L /I Figure A1.
Plots of the metric element ratios
K/I , L/I , and √ K + L /I versus skyposition of targeted source, along with cumulative probability distributions of theseratios, assuming a randomly chosen sky location, using the assumption of a 25-hourobservation with LIGO Hanford Observatory (H1) beginning 2004 Jan 1 at 00:00 UTC(GPS time 756950413). L/I can approach 0 .
50 near the celestial equator. However, this is specific to the choiceof single-detector observations with H1 only. If we assume equal amounts of data fromLIGO Hanford (H1) and LIGO Livingston (L1), we find that
L/I (cid:46) .
15 over the entiresky. We also notice that K = 0 for this choice of observing time. This is a geometricalresult related to the symmetries of the quantity a X b X under rotations of the Earth.To give a more realistic example of a typical observing time, we consider the H1 andL1 segments associated with advanced LIGO’s first observing run (O1) (cid:93) , from the LIGOOpen Science Center[7]. We see that the ratios K/I and
L/I , plotted in figure A4, aresmall enough that a Taylor expansion should be promising.As a worst-case example (and an illustration of why this approximation is bettersuited to long continuous-wave observations than to transients), in figure A5, we showthe relevant metric component ratios for an observation at a single time, assumed tocorrespond to sidereal time 00:00 at the prime meridian. We see that in this case, thebound √ K + L ≤ I is nearly saturated for much of the sky. (cid:93) https://doi.org/10.7935/K57P8W9D nalytic B -stat Approximation − . − . − . . . . . S c a l e d M e t r i c C m p t Perfect averaging, H1
IKL -90 -60 -30 0 30 60 90Declination (degrees) − . − . − . . . . . M e t r i c C m p t R a t i o K/IL/I − . − . . . . . . . . . . C u m u l a t i v e D i s t r i bu t i o n Perfect averaging, H1
K/IL/I √ K + L /I Figure A2.
Left: Plots of metric elements I , K , and L , and the ratios K/I and
L/I versus declination of targeted source, assuming an observation using LIGO HanfordObservatory (H1) that results in a perfect average over sidereal time. The spacing indeclination is chosen to be proportional to sky area. Right: Cumulative probabilitydistributions of the metric element ratios
K/I , L/I , and √ K + L /I for this case,assuming a randomly chosen sky location. − . − . − . . . . . S c a l e d M e t r i c C m p t Perfect averaging, H1 and L1
IKL -90 -60 -30 0 30 60 90Declination (degrees) − . − . − . . . . . M e t r i c C m p t R a t i o K/IL/I − . − .
05 0 .
00 0 .
05 0 . . . . . . . C u m u l a t i v e D i s t r i bu t i o n Perfect averaging, H1 and L1
K/IL/I √ K + L /I Figure A3.
Left: Plots of metric elements I , K , and L , and the ratios K/I and
L/I versus declination of targeted source, assuming an observation using LIGO HanfordObservatory (H1) and LIGO Livingston Observatory (L1) that results in a perfectaverage over sidereal time. The spacing in declination is chosen to be proportionalto sky area. Right: Cumulative probability distributions of the metric element ratios
K/I , L/I , and √ K + L /I for this case, assuming a randomly chosen sky location. nalytic B -stat Approximation K/I for O1 segments -0.159683 0.159683
L/I for O1 segments -0.159683 0.159683 √ K + L /I for O1 segments − . − . − .
05 0 .
00 0 .
05 0 .
10 0 . . . . . . . C u m u l a t i v e D i s t r i bu t i o n O1 segments
Figure A4.
Plots of the metric element ratios
K/I , L/I , and √ K + L /I versus skyposition of targeted source, along with cumulative probability distributions of theseratios, assuming a randomly chosen sky location, for an observation corresponding tothe data segments (H1 and L1) from Advanced LIGO’s first observing run (O1). Appendix B. Derivation of Taylor Expansion
Here we collect the detailed derivation of the Taylor-expanding B -statistic.In terms of the polar representation, A ˘ µ M ˘ µ ˘ ν (cid:98) A ˘ ν = IA r (cid:98) A r cos( φ r − (cid:98) φ r ) + J A l (cid:98) A l cos( φ l − (cid:98) φ l )+ A r (cid:98) A l (cid:104) K sin( φ r − (cid:98) φ l ) + L cos( φ r − (cid:98) φ l ) (cid:105) + A l (cid:98) A r (cid:104) − K sin( φ l − (cid:98) φ r ) + L cos( φ l − (cid:98) φ r ) (cid:105) (B.1)and [see eqn (5.10) of [3]] A ˘ µ M ˘ µ ˘ ν A ˘ ν = IA r + J A l + 2 A r A l [ K sin( φ r − φ l ) + L cos( φ r − φ l )] (B.2) nalytic B -stat Approximation K/I for 00:00 GST (H1,L1) -1 1
L/I for 00:00 GST (H1,L1) -1 1 √ K + L /I for 00:00 GST (H1,L1) − . − . − . − .
25 0 .
00 0 .
25 0 .
50 0 .
75 1 . . . . . . . C u m u l a t i v e D i s t r i bu t i o n Single time, 00:00 GST (H1,L1)
K/IL/I √ K + L /I Figure A5.
Plots of the metric element ratios
K/I , L/I , and √ K + L /I versus skyposition of targeted source, along with cumulative probability distributions of theseratios, assuming a randomly chosen sky location, for a brief observation at Greenwichsidereal time 00:00. so thatΛ( A ; x ) = I (cid:18) − A r + A r (cid:98) A r cos( φ r − (cid:98) φ r ) (cid:19) + J (cid:18) − A l + A l (cid:98) A l cos( φ l − (cid:98) φ l ) (cid:19) + K (cid:16) − A r A l sin( φ r − φ l ) + A r (cid:98) A l sin( φ r − (cid:98) φ l ) − (cid:98) A r A l sin( φ l − (cid:98) φ r ) (cid:17) + L (cid:16) − A r A l cos( φ r − φ l ) + A r (cid:98) A l cos( φ r − (cid:98) φ l ) + (cid:98) A r A l cos( φ l − (cid:98) φ r ) (cid:17) = Λ r ( A r , φ r ; (cid:98) A r , (cid:98) φ r ) + Λ l ( A l , φ l ; (cid:98) A l , (cid:98) φ l )+ (cid:104) K sin( (cid:98) φ r − (cid:98) φ l ) + L cos( (cid:98) φ r − (cid:98) φ l ) (cid:105) × (cid:104) A r A l (cid:16) − cos( φ r − (cid:98) φ r ) cos( φ r − (cid:98) φ r ) + sin( φ r − (cid:98) φ r ) sin( φ l − (cid:98) φ l ) (cid:17) + A r (cid:98) A l cos( φ r − (cid:98) φ r ) + (cid:98) A r A l cos( φ l − (cid:98) φ l ) (cid:105) + (cid:104) K cos( (cid:98) φ r − (cid:98) φ l ) + L sin( (cid:98) φ r − (cid:98) φ l ) (cid:105) × (cid:20) A r A l (cid:16) cos( φ r − (cid:98) φ r ) sin( φ l − (cid:98) φ l ) + sin( φ r − (cid:98) φ r ) cos( φ l − (cid:98) φ l ) (cid:17) + A r (cid:98) A l sin( φ r − (cid:98) φ r ) + (cid:98) A r A l sin( φ l − (cid:98) φ l ) (cid:21) (B.3) nalytic B -stat Approximation e Λ( {A ˘ µ } ; x ) = e Λ r ( A r ,φ r ; (cid:98) A r , (cid:98) φ r )+Λ l ( A l ,φ l ; (cid:98) A l , (cid:98) φ l )+Λ ( A ; (cid:98) A ) ≈ e Λ r ( A r ,φ r ; (cid:98) A r , (cid:98) φ r ) e Λ l ( A l ,φ l ; (cid:98) A l , (cid:98) φ l ) (cid:16) ( A ; (cid:98) A ) (cid:17) (B.4)In this form, we can factor the integrals in each of the terms; they all reduce to one ofthree forms: (cid:90) π (cid:90) ∞ e − I A + IA (cid:98) A cos( φ − (cid:98) φ ) dA dφ √ A = 2 π (cid:90) ∞ e − I A I ( IA (cid:98) A ) dA √ A (B.5a) (cid:90) π (cid:90) ∞ e − I A + IA (cid:98) A cos( φ − (cid:98) φ ) cos( φ − (cid:98) φ ) √ A dA dφ = 2 π (cid:90) ∞ e − I A I ( IA (cid:98) A ) √ A dA (B.5b) (cid:90) π (cid:90) ∞ e − I A + IA (cid:98) A cos( φ − (cid:98) φ ) sin( φ − (cid:98) φ ) √ A dA dφ = 0 (B.5c)where I n ( x ) = i − n J n ( ix ) is the modified Bessel function, and we have used the Jacobi-Anger expansion[4], which tells us that e IA (cid:98) A cos( φ − (cid:98) φ ) = I ( IA (cid:98) A ) + 2 ∞ (cid:88) n =1 I n ( IA (cid:98) A ) cos( n [ φ − (cid:98) φ ]) . (B.6)Both of the remaining integrals can be done using equation (11.4.28) of [4], which says,in terms of the modified Bessel function, that, when Re( ν + µ ) > a ) > (cid:90) ∞ e − a t t µ − I ν ( bt ) dt = Γ (cid:0) ν + µ (cid:1) (cid:0) b a (cid:1) ν a µ Γ( ν + 1) F (cid:18) ν + µ , ν + 1 , b a (cid:19) (B.7)where F ( a, b, z ) = M ( a, b, z ) is the confluent hypergeometric function. We apply thiswith a = I/ b = I (cid:98) A , µ = 1 / /
2, and ν = 0 and 1, respectively, in (B.5a) and(B.5b), to get (cid:90) π (cid:90) ∞ e − I A + IA (cid:98) A cos( φ − (cid:98) φ ) dA dφ √ A = 2 π Γ (cid:0) (cid:1) / I / F (cid:32) , , I (cid:98) A (cid:33) (B.8a) (cid:90) π (cid:90) ∞ e − I A + IA (cid:98) A cos( φ − (cid:98) φ ) cos( φ − (cid:98) φ ) √ A dA dφ = 2 π Γ (cid:0) (cid:1) (cid:98) A / I / F (cid:32) , , I (cid:98) A (cid:33) (B.8b) nalytic B -stat Approximation B -statistic (2.23) as B ( x ) ≈ A π (cid:32) π Γ (cid:0) (cid:1) / I / (cid:33) (cid:32) π Γ (cid:0) (cid:1) / J / (cid:33) F (cid:32) , , I (cid:98) A r (cid:33) F (cid:32) , , J (cid:98) A l (cid:33) × (cid:104) K sin( (cid:98) φ r − (cid:98) φ l ) + L cos( (cid:98) φ r − (cid:98) φ l ) (cid:105) (cid:98) A r (cid:98) A l × F (cid:16) , , I (cid:98) A r (cid:17) F (cid:16) , , I (cid:98) A r (cid:17) + 14 F (cid:16) , , J (cid:98) A l (cid:17) F (cid:16) , , J (cid:98) A l (cid:17) − F (cid:16) , , I (cid:98) A r (cid:17) F (cid:16) , , I (cid:98) A r (cid:17) F (cid:16) , , J (cid:98) A l (cid:17) F (cid:16) , , J (cid:98) A l (cid:17) (B.9) Appendix C. Recovery of F -Statistic Our method expands the B -statistic to first order in the metric components K and L . Ithas been shown in [2] that the Bayes factor constructed with a prior uniform in the {A ˘ µ } is equivalent to the F -statistic, which we note in (2.20) has only zeroth- and first-orderterms in these quantities. This means that applying the Taylor-expansion method withthis prior should reproduce the exact F -statistic.If we replace the isotropic prior (2.22) with a uniform prior pdf( A ˘1 , A ˘2 , A ˘3 , A ˘4 |H f ) = C , the B -statistic integral (2.23) becomes B ( x ) = C (cid:90) ∞−∞ (cid:90) ∞−∞ (cid:90) ∞−∞ (cid:90) ∞−∞ e Λ( {A ˘ µ } ; x ) d A ˘1 d A ˘2 d A ˘3 d A ˘4 = C (cid:90) π (cid:90) π (cid:90) ∞ (cid:90) ∞ e Λ( {A ˘ µ } ; x ) A r A l dA r dA l dφ r dφ l (C.1)The Taylor expansion of the likelihood, and the angular integrals, proceed as inAppendix B, and the only difference is that the two principal integrals (B.5a) and(B.5b), become (cid:90) π (cid:90) ∞ e − I A + IA (cid:98) A cos( φ − (cid:98) φ ) A dA dφ = 2 π (cid:90) ∞ e − I A I ( IA (cid:98) A ) A dA (C.2a) (cid:90) π (cid:90) ∞ e − I A + IA (cid:98) A cos( φ − (cid:98) φ ) cos( φ − (cid:98) φ ) A dA dφ = 2 π (cid:90) ∞ e − I A I ( IA (cid:98) A ) A dA (C.2b)Using (B.7) with a = I/ b = I (cid:98) A , µ = 2 and 3, and ν = 0 and 1, respectively, we find,in place of (C.3), (cid:90) π (cid:90) ∞ e − I A + IA (cid:98) A cos( φ − (cid:98) φ ) A dA dφ = 2 πI F (cid:32) , , I (cid:98) A (cid:33) = 2 πI e I (cid:98) A (C.3a) (cid:90) π (cid:90) ∞ e − I A + IA (cid:98) A cos( φ − (cid:98) φ ) cos( φ − (cid:98) φ ) A dA dφ = 2 π (cid:98) AI F (cid:32) , , I (cid:98) A (cid:33) = 2 π (cid:98) AI e I (cid:98) A (C.3b) nalytic B -stat Approximation F ( a, a, z ) = e z . This then givesa statistic of B ( x ) ≈ C (cid:18) πI (cid:19) (cid:18) πJ (cid:19) e I (cid:98) A r + J (cid:98) A l (cid:110) (cid:104) K sin( (cid:98) φ r − (cid:98) φ l ) + L cos( (cid:98) φ r − (cid:98) φ l ) (cid:105) (cid:98) A r (cid:98) A l (1 + 1 − (cid:111) (C.4)So that, to first order in K and L ,ln B ( x ) B (0) ≈ I (cid:98) A r J (cid:98) A l (cid:104) K sin( (cid:98) φ r − (cid:98) φ l ) + L cos( (cid:98) φ r − (cid:98) φ l ) (cid:105) (cid:98) A r (cid:98) A l (C.5)Which is indeed the form given in (2.20) for the exact F -statistic. Appendix D. Relationship to High-SNR Approximation
Recent work[8] by Dhurandhar, Krishnan and Willis (hereafter DKW) contains adifferent approximate expression for the B -statistic, derived in the limit of high signal-to-noise ratio, but without assumptions on the form of the metric. In their notation,the approximate form is written [[8] equation (104)] B ( x ) ≈ (cid:18) π ζ − k ) (cid:19) (cid:34) e (cid:98) B † N (cid:98) B ( | (cid:98) B || (cid:98) B | ) (cid:35) (D.1)To make contact with our results, we collect here the conversion between DKW’snotation and ours. Their metric elements are ζ = I = J (they limit attention tothe long-wavelength limit) and κ = L + iK , with k = | κ | = √ K + L . They definecomplex amplitudes B = h e − iφ (1 + χ ) e − iψ = A r e − i (3 φ r + φ l ) = B ∗ (D.2a) B = h e − iφ (1 − χ ) e iψ = A l e − i ( φ r +3 φ l ) = B ∗ (D.2b)and a complex metric N ≡ { N µν } = 12 ζ κ ∗ κ ζ ζ κ ∗ κ ζ (D.3)from which we see12 (cid:98) B † N (cid:98) B = Re (cid:16) (cid:98) A r e i (3 (cid:98) φ r + (cid:98) φ l ) [ L − iK ] (cid:98) A l e − i ( (cid:98) φ r +3 (cid:98) φ l ) (cid:17) = 12 I (cid:98) A r + 12 J (cid:98) A l + (cid:98) A r (cid:98) A l (cid:104) K sin( (cid:98) φ r − (cid:98) φ l ) + L cos( (cid:98) φ r − (cid:98) φ l ) (cid:105) = F (D.4)and therefore their approximation can be written B ( x ) ≈ (cid:18) π IJ − K − L ) (cid:19) (cid:34) e F ( (cid:98) A r (cid:98) A l ) (cid:35) (D.5) nalytic B -stat Approximation A = 2 π .) References [1] Jaranowski P, Kr´olak A and Schutz B F 1998
Phys. Rev. D. Preprint gr-qc/9804014 )[2] Prix R and Krishnan B 2009
Class. Quant. Grav. Preprint arXiv:0907.2569 )[3] Whelan J T, Prix R, Cutler C J and Willis J L 2014
Class. Quant. Grav. Preprint )[4] Abramowitz M and Stegun I A 1964
Handbook of Mathematical Functions (National Bureau ofStandards)[5] Dergachev V 2012
Phys. Rev. D. Preprint )[6] Searle A C 2008
ArXiv E-Prints ( Preprint )[7] Vallisneri M, Kanner J, Williams R, Weinstein A and Stephens B 2015
Journal of Physics:Conference Series
Preprint arXiv:1410.4839 )[8] Dhurandhar S, Krishnan B and Willis J L 2017 (
Preprint1707.08163