= =22.0truecm =15.5truecm =0.8truecm =0.2truecm ø[$\bullet~$]{} =cmbx12 =cmr9 =cmr8

[LISA data analysis: The monochromatic binary]{}

[detection and initial guess problems]{}

Ronald W. Hellings

Department of Physics, Montana State University, Bozeman MT 59715

ABSTRACT: We consider the detection and initial guess problems for the LISA gravitational wave detector. The detection problem is the problem of how to determine if there is a signal present in instrumental data and how to identify it. Because of the Doppler and plane-precession spreading of the spectral power of the LISA signal, the usual power spectrum approach to detection will have difficulty identifying sources. A better method must be found. The initial guess problem involves how to generate [*a priori*]{} values for the parameters of a parameter-estimation problem that are close enough to the final values for a linear least-squares estimator to converge to the correct result. A useful approach to simultaneously solving the detection and initial guess problems for LISA is to divide the sky into many pixels and to demodulate the Doppler spreading for each set of pixel coordinates. The demodulated power spectra may then be searched for spectral features. We demonstrate that the procedure works well as a first step in the search for gravitational waves from monochromatic binaries.

The signal coming down from a spaceborne gravitational wave detector will contain signals from many close compact binaries, but, to look at the signal, one would see nothing but noise. This is because the instrument noise, integrated over the bandwidth of the detector, far exceeds the strength of any individual binary signal. This situation is common in data analysis, and the usual procedure for detecting weak coherent signals in the presence of large incoherent noise is to generate the power spectrum of the signal. Since the binary star signal is monochromatic, its power will be concentrated in a single frequency bin of the power spectrum, while the noise is spread over all frequencies. The signal will then be seen as a single spike rising above the neighboring noise spectrum. However, the usefulness of this power spectrum method for signal detection is limited for the spaceborne gravitational wave detectors, due to the motion of the detector around the sun and to the change in the orientation of the detector. As a result of these motions, the gravitational wave signal from a particular source will be strongly phase modulated, spreading the total power from the source over many frequency bins. It is therefore quite possible that none of the frequency components of the signal will be seen above the noise.

To illustrate, we show in Fig. 1 the power spectrum of the gravitational wave signal from a binary star, where the signal frequency is 0.01 Hz and the amplitude is $3.5\times10^{-22}$. The signal has been sampled at 10 s intervals for one year. This is the power spectrum that would be seen if the detector were at rest in the solar system. The power spectrum as seen in a moving detector is shown in Fig. 2. To generate this figure, the detector has been assumed to be the LISA detector \[1\], in which three free-flying spacecraft move on heliocentric trajectories in such a way as to remain at roughly constant distances from each other and to remain in a plane inclined by $60\deg$ to the plane of the ecliptic, the plane precessing about the ecliptic pole once per year. The spreading of the power over about 100 frequency bins can be clearly seen in the figure. If a noise spectrum, characteristic of the proposed LISA instrument noise, is added to the signal, the power spectrum in Fig. 3 results. As may be seen, the spreading of the binary signal power over many frequency bins has produced a power spectrum that cannot now be seen above the noise.

One of the goals of this paper is to show how this “detection problem” may be solved.

The ultimate goal of gravitational wave detectors is to use the gravitational wave signals to determine the astrophysical properties of the gravitational wave source. Papers by Cutler and Vecchio \[2,3\] and by Moore and Hellings \[4\] have addressed the problem of what information is available in the signal and how well the parameters of the source may be determined from a signal with a particular signal-to-noise ratio. The data analysis technique to be used to extract the information from the signal is called linear least-squares parameter estimation. This technique accounts for the possible correlation of the parameters with each other and gives the best-fit values of the parameters and their formal standard deviation, with the assumption of a background of Gaussian noise. However, the “linear” in “linear least-squares” reminds us that this procedure will only work in a totally linear problem or in the limit that the problem may be linearized by writing it in terms of small quantities, the differences between the final best-fit values for the parameters and some [*a priori*]{} values.

The signal from a monochromatic binary is characterized by seven independent parameters. These are:

øthe frequency $f$ of the gravitational wave signal (a parameter that depends on the masses of the two members of the binary and on the radius of their mutual orbit),

øthe amplitude $h$ (which depends of the binary masses, the orbital radius, and on the distance from the binary system to the Earth),

øthe two-parameter location of the binary in the sky (ecliptic coordinates $\theta$ and $\phi$),

øthe two-parameter orientation of the orbit plane of the binary (taken to be the inclination $i$ of the orbit plane normal to the line-of-sight to the binary and the angle $\psi$ by which the major axis of the binary’s apparent ellipse in the sky is rotated relative to the ecliptic plane, around the line-of-sight), and

øthe phase $\phi_0$ of the binary orbit at some initial time $t_0=0$.

The received gravitational wave signal is linear in only one of these seven parameters, the amplitude $h$. For all the rest, initial guesses must be found that are sufficiently close to the correct values of the parameters to allow the linear least-squares analysis to converge to the proper place. Often, in parameter estimation processes like this, the initial guesses are known from some other aspect of the observations. For example, in a similar problem to the one we are facing, the times-of-arrival of pulsar pulses are used to determine the position of the pulsar in the sky and the parameters characterizing pulsar’s rotational dynamics. The initial guess for the position, however, is not found from the time series itself, but from the direction in which the antenna is pointing when the pulsar data is collected. The spaceborne gravitational wave detector is an omnidirectional instrument, so there will be no such ancillary information on source direction; all information must be extracted from the time series itself.

A critical question for LISA data analysis, therefore, is how does one take the time series and generate initial guesses for the six non-linear parameters. One method would be to cross-correlate the time series with a series of templates that represent all possible values of all possible parameters and look for the templates that have a strong correlation with the signal, the parameters of the template becoming the initial guess. The spacing between templates is set by the need to get an initial guess that is within the linear regime of the final answer. Ultimately, this is the only method that can be employed in the absence of data from outside the time series. However, let us consider the nature of the six non-linear parameters to see if a more step-by-step method may be found.

The frequency $f$ of a signal buried in noise is ordinarily found by spectral analysis. The reason why that avenue is not available here is evident by inspection of Figs. 2 and 3; the spectral power is spread out over many frequency bins. The contribution to this spreading by the motion of the detector around the sun is given by the usual Doppler formula $$f_r = f\left[1-{\Omega R\over c}\sin\theta \sin(\Omega t-\phi)\right]
\eqno(1)$$ where $f_r$ is the frequency observed in the detector, $f$ is the frequency that would be observed at the solar system barycenter, $\Omega$ is the mean motion of the detector around the sun, and $R$ is the radius of the detector orbit. The magnitude of the Doppler spreading depends on the gravitational wave frequency $f$ and on $\theta$. The phase of the Doppler depends on $\phi$. For gravitational waves at frequencies higher than $f=c/(2R)(\approx10^{-3}\,$Hz for LISA) this Doppler shift dominates over the spreading produced by the precession of the detector plane and is the major impediment to determination of the frequency by spectral analysis. Now, if the frequency and position of a source were to be known, the template matching could be a much simpler procedure, with only three parameters left to determine – $i$, $\psi$, and $\phi_0$. What is more, these three parameters have little effect on the signal in the detector and are much less well-determined in the final least-squares analysis than are $\theta$ and $\phi$. Therefore the resolution needed to find initial values for these parameters that lie within the linear regime, where next-higher-order terms are not greater than the final uncertainty, is very coarse, and the number of templates needed to span the parameter space will be small. A method to independently find the position and frequency would thus greatly simplify the template matching procedure.

The second goal of this paper is to show how to generate such a step-by-step approach to the ”initial guess” problem by approximating the frequency and location of the source independently of the other parameters of the system.

Both problems we have discussed, the detection and initial guess problems, may be solved by first undoing the position-dependent Doppler modulation the signal undergoes. The Doppler modulation depends only on $\theta$ and $\phi$ and is given simply by Eq. 1. The modulation produced by the detector-plane precession is much more complicated (see Moore and Hellings \[4\]) and depends on all four angular parameters, $\theta$, $\phi$, $i$, and $\psi$, so demodulating the Doppler only is both a simple task to undertake and one that will undo the greatest part of the spectral spreading at frequencies where the spreading is most pronounced.

Given a position on the sky, and given a time series from the detector, a new time series is generated which is the time series that would be received from that direction if the detector were at rest at the solar system barycenter. The plane of this fictitious detector will continue to precess, however, at the yearly period, just like the true LISA detector. If the time series at the detector is $y(t_n)$, we seek a time series $z(t_n)$ that represents the signal that the barycentric detector would have received. The algorithm for generating $z(t_n)$ is $$z(t_n)=y(t_n-\tau), \eqno(2)$$ where $\tau$ is the time of flight from when the wavefront passed the detector to when it arrives at the barycenter, as seen in Fig. 4. The value of $\tau$, given $t_n$, comes from the solution of the transcendental equation $$\tau={\hat n\cdot\vec r(t_n-\tau)\over c},  \eqno(3)$$ where $\hat n$ is the unit vector toward the source and $\vec r(t)$ is the vector from the barycenter to the detector at time $t$. At time $t_n$, the value of $\tau$ is found by solving Eq. 3, and the barycentric time series is generated using Eq. 2.

This procedure, of course, assumes that one knows the direction $\hat n$ to the source. When the direction is known, and the procedure applied, a time series is generated whose Fourier spectrum is shown in Fig. 5. As may be seen in comparison with Fig. 3, the spread spectrum that lay beneath the noise in the uncorrected time series is now seen strongly above the noise. The breadth of the spectral “line” in Fig. 5 is a result of the uncorrected modulation due to the precession of the detector plane. Indeed, when we produce a demodulated spectrum for the OMEGA detector \[4\], where the detector plane does not precess, the entire signal power is found in a single line, as seen in Fig. 6.

Clearly, the remaining problem is that of determining the direction to the source. With no [*a priori*]{} knowledge, the only course of action is to consider all directions one at a time, to demodulate for that direction, and to examine the power spectrum of the resulting demodulated time series. This has the effect of solving the detection problem and the initial guess problem simultaneously. The demodulated signals will collect the power into a few frequency bins where they will protrude above the noise spectral power, the frequency at which this occurs will be the frequency of the source, and the particular direction that produces this visible spectral feature will be the direction to the source.

To illustrate this procedure we have simulated signals from three binary systems, as they would be observed in the LISA detector. We have combined these signals together, as they would appear in the detector, and then added random Gaussian noise at a level consistent with the LISA detector in this frequency band. The time series is sampled every 10 s and accumulated for a year ($3.16\times10^6$ data points). The sky is then divided into pixels of one square degree, for a total of 20,582 pixels. The pixel coordinates are used one at a time to demodulate the data according to Eqs. 2-3, and the power spectrum is formed for each case by fast Fourier transform of the demodulated time series. The power spectrum is then searched for features that are a chosen factor above the average power in the spectrum. Since the spectral lines are still spread over several adjoining bins by the detector plane precession, the detection procedure takes three-point running averages of the power spectrum to enhance such features as it searches. The result of this procedure is seen in Fig. 7.

The limits of the sky plot in Fig. 7 are from the ecliptic plane ($\theta=90\deg$) to the ecliptic pole ($\theta=0\deg$) and over 360$\deg$ in ecliptic longitude. The negative ecliptic lattitudes ($\theta>90\deg$) are the mirror image of the positive ones due to the symmetry of $\sin\theta$ in Eq. 1. The resolution of the plot is one square degree. The three sources, at $\{\theta,\phi\}$ coordinates $\{60,120\}$, $\{78,20\}$ and $\{12,280\}$ are seen clearly above the noise in the plots. However, the size and complicated structure of each region make it difficult to select a value for the coordinates of the source to the same accuracy as the resolution of the plot. In each of the three regions, a search for the brightest pixel for each source will give values that are within $\pm10\deg$ of the correct values, but this is not sufficient to be an initial guess for $\theta$ and $\phi$. This may be seen by considering the Doppler formula in Eq. 1. The change in the Doppler shift of the signal due to a change in $\theta$, for example, is proportional to $$\delta(\sin\theta)=
\cos\theta\,\delta\theta-{1\over2}\sin\theta\,\delta\theta^2 + {\rm
Higher}\, {\rm Order}\, {\rm Terms} \eqno (4)$$ The second term on the right represents the non-linear error introduced into the linear least-squares fit. For the source at $\theta=60\deg$, this non-linear term would be $1.3\times10^{-2}$ radians for $\delta\theta=10\deg$. This should be compared with the error in the first term due to the final uncertainty in the linear least-squares parameter fit, $\cos\theta\,\sigma_\theta=(0.5)(5\times10^{-3})=2.5\times10^{-3}$ (from the least-squares fits of Moore and Hellings \[4\]). Thus, when a linear least-squares fit tries to make a $10\deg$ change from initial guess to final value for $\theta$, it will use only the first term in Eq. 4 and will make an error five times its calculated formal uncertainty. The situation for $\phi$ is identical. From this, we conclude that one cannot go directly from the position estimates taken from the demodulated signal to a final least-squares fit.

So what is to be done? At this point we need to remember that our goal in the demodulation procedure was not to generate initial guesses directly, but to reduce the range of templates that would be needed for a template-matching procedure. Indeed, the complicated structure that keeps us from determining a single accurate position via the demodulation is a result of the precession of the detector plane that depends on $\theta$, $\phi$, $i$, and $\psi$. For the non-precessing OMEGA detector, the demodulation alone is sufficient to get the initial position guess to the required accuracy, as may be seen by examination of Fig. 8. For LISA, the result of the demodulation has only been to find an exact value for the frequency and approximate values for $\theta$ and $\phi$. However, this will enable one to reduce the required template space from ($10^6$ frequencies)$\times$($2\times10^4$ positions) to less than 100 position templates. A five-parameter template matching procedure ($f$ is now known and $h$ is always linear) over a restricted set of values for $\theta$ and $\phi$ will be needed to produce a sufficient best guess to be supplied to the final least-squares estimator.

Several other things may also be noted in Fig. 7. First, a single bright source will dominate the observed demodulated power over a large range of values of $\theta$ and $\phi$, possibly masking weaker nearby sources. Part of this is an artifact of the method we have used to generate Fig.7, the picking out of the strongest spectral feature to represent each pixel. Perhaps the wisest use of the figure is to eliminate regions where there is no significant demodulated power as being a region free of sources that reach above the noise at all. One may then concentrate on the regions with significant power and examine spectra from many points in the region to look for secondary sources with different frequencies. A second thing to be noticed Fig. 7 is the significant leaking of power from sources near the ecliptic plane ($\theta=90\deg$) down towards the ecliptic plane. This is a result of the diminished position accuracy for the least squares fits for monochromatic binaries as one approaches the ecliptic (see Moore and Hellings \[4\], Fig. 3). Looking at it another way, this is a result of the fact that the phase of the Doppler modulation is the same for sources at the same ecliptic longitudes; only the amplitude of the modulation is different.

Finally, let us consider some of the questions that remain to be addressed. First, the details of the final template matching need to be studied. Cornish and Larson \[5\] have generated an efficient code for template matching that works in the frequency domain and that is capable of fast correlation of signals over the entire parameter space. They have also generated a Doppler demodulation code\[6\], that similarly works in the frequency domain, whose efficiency should be compared with the time-domain code that we have used in this paper. Lastly, we should point out that among the test cases we have run are some in which there are two or more sources nearby each other in frequency and in position. When this occurs, there seems to be a pulling of the best demodulated pixels away from both source positions to some intermediate position. This has implications not only for the initial guess problem, but also for the problem of parameter estimation in the presence of non-Gaussian noise. These are issues we intend to address in future papers.

It is a pleasure to thank Neil Cornish and Thomas Moore for helpful discussions. This study was accomplished with support from NASA grants NAG5-11469 and NCC5-579.

=cmbx10

=0.2in =1 \[1\] K. Danzmann [*et al.,*]{} [*LISA Pre-Phase A Report (second edition)*]{}, MPQ 233 (1998).

=0.2in =1 \[2\] C. Cutler and A. Vecchio, [*Proceedings of the 2$^{nd}$ LISA Symposium (AIP Conference Proceedings 456)*]{}, 95 (1998).

=0.2in =1 \[3\] A. Vecchio and C. Cutler, [*Proceedings of the 2$^{\rm nd}$ LISA Symposium (AIP Conference Proceedings 456)*]{}, 101 (1998).

=0.2in =1 \[4\] T.A. Moore and R.W. Hellings, Phys. Rev. [D65]{} 062001 (2002).

=0.2in =1 N.J. Cornish and S.L. Larson, manuscript in preparation (2002).

=0.2in =1 N.J. Cornish and S.L. Larson, manuscript in preparation (2002).

CAPTIONS:

Figure 1. A portion of the power spectrum of the gravitational wave from a binary star, as seen by a stationary, non-precessing detector at the barycenter. The parameters characterizing the gravitational wave source are $h=3.5\times10^{-22}$, $f=1.0\times10^{-2}\,$Hz, $\theta=60\deg$, $\phi=120\deg$, $i=10\deg$, $\psi=84\deg$, and $\phi_0=0$.

Figure 2. The same gravitational wave as in Fig. 1, as seen by the LISA detector. The power is now spread over about one hundred frequency bins.

Figure 3. The same gravitational wave in the LISA detector as in Fig.2, but with noise added at a level characteristic of the LISA detector in this frequency band. The spread spectrum cannot now be detected above the noise.

Figure 4. The geometry of the revision of the wavefront timing to the solar system barycenter. In the first picture, the LISA detector receives the signal at time $t=t_n-\tau$. In the second picture, a time $\tau$ later, the same signal is received at the barycenter, while LISA has moved to a new position.

Figure 5. The signal from Fig. 3, after the revision to the solar system barycenter has been performed for the particular direction of this source. The signal at 0.01 Hz is clearly visible above the noise. The remaining breadth of the spectral peak is due to the modulation produced by the precessing plane of the LISA detector.

Figure 6. The same as Fig. 5, but for the non-precessing OMEGA detector. The entire gravitational wave power is now concentrated in a single peak.

Figure 7. The all-sky map of peak spectral power in the demodulated signal for one-degree pixels. The signal being analyzed represents the received LISA time series for three binary gravitational wave sources plus random noise at the predicted LISA level. The parameters of the three sources are: SOURCE 1 — $h=3.3\times10^{-22}$, $f=3.0\times10^{-3}\,$Hz, $\theta=60\deg$, $\phi=120\deg$, $i=10\deg$, $\psi=84\deg$, $\phi_0=0$; SOURCE 2 — $h=4.0\times10^{-22}$, $f=4.0\times10^{-3}\,$Hz, $\theta=78\deg$, $\phi=20\deg$, $i=40\deg$, $\psi=0\deg$, $\phi_0=0$; and SOURCE 3 — $h=3.7\times10^{-22}$, $f=3.5\times10^{-3}\,$Hz, $\theta=12\deg$, $\phi=280\deg$, $i=0\deg$, $\psi=84\deg$, $\phi_0=0$.

Figure 8. The same type of all-sky map for the same signals as Fig. 7, but for the non-precessing OMEGA detector. The relative weakness of the power near the ecliptic (the $\theta=78\deg$, $\phi=20\deg$ location) is due to the non-uniform sensitivity of the OMEGA detector – more sensitive than LISA near the ecliptic pole and less sensitive near the ecliptic plane.
