text:
# An excess power statistic for detection of burst sources of gravitational radiation

## I Introduction and Summary

### A Background and Motivation

The inspiral, merger, and ringdown of binary black hole systems may be the most important source of gravitational radiation for detection by the kilometer-scale interferometric gravitational wave detectors such as LIGO , and VIRGO . The importance of these sources is twofold :

1. A large amount of gravitational radiation is expected to be emitted by the merger of two black holes. For intermediate mass ($`10M_{}\text{}1000M_{}`$) black hole binaries this radiation will be in the frequency band of highest sensitivity for LIGO and VIRGO.***While the relative abundance of such systems is still a very open question, we are encouraged by two recent developments in the astrophysics literature: (i) evidence suggesting that black holes in this mass range may exist  and (ii) a globular cluster model that suggests LIGO I may expect to see about one black hole coalescence event during the first two years of operation . These sources should therefore be amongst the brightest in the sky, and visible to much greater distances than other sources. The detection rate for coalescing binary black holes could therefore be higher than for any other source.
2. The radiation emitted from the merger of black holes probes the strong field regime of a purely gravitational system. This radiation should therefore provide a sensitive test of general relativity.

These benefits can only be realized, however, if the gravitational radiation from black hole mergers can be detected.

The best understood and most widely developed technique for detection of gravitational waves with interferometric detectors is matched filtering . Matched filtering is the optimal technique if the entire waveform to be detected is accurately known in advance (up to a few unknown parameters). Unfortunately, the gravitational radiation from black hole mergers results from highly non-linear self-interaction of the gravitational field. This makes it extremely difficult to obtain gravitational waveforms. Efforts to do so have met with only limited success thus far. Binary black hole mergers will therefore not amenable to detection by matched filtering, at least for the first gravitational wave searches in the 2002-2004 time frame.

Similarly, there are other classes of sources, like core-collapse of massive stars in supernovae, or the accretion induced collapse of white dwarfs, for which the physics is too complex to allow computation of detailed gravitational waveforms. For these sources, as for binary black hole mergers, we must seek alternative signal detection methods. These methods are often called “blind search” methods.

One class of search methods is based on time-frequency decompositions of detector data. (For an exploration of a variety of other methods, see Ref. .) Time-frequency strategies have become standard in many other areas of signal analysis . There is also a growing literature on the application of time-frequency methods to gravitational waves . For binary black hole mergers, it is possible to make crude estimates of signals’ durations and frequency bands , although these estimates need to be firmed up and refined by numerical relativity simulations. This suggests that one should look only in the relevant time-frequency window of the detector output.

Flanagan and Hughes  (FH) have suggested a particular time-frequency method for blind searches. The method uses only knowledge of the duration and frequency band of the signal: one simply computes the total power within this time-frequency window, and repeats for different start times. The method detects a signal if there is more power than one expects from detector noise alone. Thus, we call it the *excess power* search method. Similar methods have been discussed elsewhere in the gravitational wave literature. Schutz  investigated the method in the context of the cross-correlation of outputs from different detectors. An autocorrelation filter for unrestricted frequencies was published by Arnaud et al.  shortly after and independently of FH. A generalization of the excess power filter has also been discussed in the signal analysis literature , where it has recently been “attracting considerable interest” . Finally a method closely related to the excess power method has recently been explored by Sylvestre .

The excess power method distinguishes itself for the detection of signals of known duration and frequency band by a single compelling feature: in the absence of any other knowledge about the signal, *the method can be shown to be optimal*. Furthermore, it can be shown that for mergers of a sufficiently short duration and narrow frequency band, it performs nearly as well as matched filtering.

The essence of the power filter is that one compares the power of the data in the estimated frequency band and for the estimated duration to the known statistical distribution of noise power. It is straightforward to show that if the detector output consists solely of stationary Gaussian noise, the power in the band will follow a $`\chi ^2`$ distribution with the number of degrees of freedom being twice the estimated time-frequency volume (i.e., the product of the time duration and the frequency band of the signal). If a gravitational wave of sufficiently large amplitude is also present in the detector output, an excess of power will be observed; in this case, the power is distributed as a non-central $`\chi ^2`$ distribution  with non-centrality parameter given by the signal power. The signal is detectable if the excess power is much greater than the fluctuations in the noise power which scales as the square-root of the time-frequency volume. Thus, the viability of the excess power method depends on the expected duration and bandwidth of the gravitational wave as well as on its intrinsic strength. For instance, the method is not competitive with matched filtering in detecting binary neutron star inspirals, since the time-frequency volume for such signals is very large, $`10^4`$.

To implement this method, one needs to decide the range of frequency bands and durations to search over. For initial LIGO, the most sensitive frequency band is $`100\text{}300\text{Hz}`$, and it makes sense to search just in this band. For binary black hole mergers, signal durations might be of order tens or hundreds of milliseconds, depending on the black hole masses and spins . Thus, the time-frequency volume of a merger signal can be as large as $`100`$, and its power would need to be more than one tenth as large as the noise power for detectability with the excess power method.

One can also establish operational lower bounds on the time durations and frequency bands of interest. Because the largest operational frequency bandwidth is $`200\text{Hz}`$ for the initial LIGO interferometers, the shortest duration of signal that need be considered is $`5\text{ms}`$ (for a minimum time-frequency volume of unity). Similarly, for a maximum duration of $`0.5\text{s}`$, the smallest bandwidth that needs to be considered is $`2\text{Hz}`$. The excess power in any of the allowed bandwidths and durations can thus be obtained by judiciously summing up power that is output from a bank of one hundred $`2\text{Hz}`$ band-pass filters (spanning the $`200\text{Hz}`$ of peak interferometer sensitivity) for the required duration.

Having established the statistic and its operational range of parameters, the following simple algorithm for implementing the excess power method emerges naturally:

1. Pick a start time $`t_s`$, a time duration $`\delta t`$ (containing $`N`$ data samples), and a frequency band $`[f_s,f_s+\delta f`$\].
2. Fast Fourier transform (FFT) the block of (time domain) detector data for the chosen duration and start time.
3. Sum the power in each of the $``$ one hundred 2 Hz bands spanning the peak sensitivity region of the detector.
4. Further sum the power in the 2 Hz bands which correspond to the chosen frequency band.
5. Calculate the probability of having obtained the summed power from Gaussian noise alone using a $`\chi ^2`$ distribution with $`2\times \delta t\times \delta f`$ degrees of freedom.
6. If the probability is significant, record a detection.
7. Repeat the process for all allowable choices of start times $`t_s`$, durations $`\delta t`$, starting frequencies $`f_s`$ and bandwidths $`\delta f`$.

This procedure, which must be repeated for every possible start time, can lead to moderately-large computational requirements. We find that the computational efficiency of this implementation, which we call the short FFT algorithm, can be improved upon by considering data segments much longer than the longest signal time duration. In this case, after summing over the chosen band, we must FFT the data back into the time domain. This implementation, which we call the long FFT algorithm, is more efficient by at least a factor of $`4`$ over the parameter space of interest.

The most significant drawback of the filter outlined above is that the $`\chi ^2`$ statistic is appropriate only to Gaussian noise. Real detector noise will contain significant non-Gaussian components. There are likely to be transient bursts of broad band noise that have characteristics very similar to black hole merger signals.

The non-Gaussianity of real detector noise leads us to two considerations. First, like most blind search methods, the excess power method will likely be a useful tool for characterizing and investigating the non-Gaussian components of the noise. In particular, it can provide a simple and automated procedure for garnering statistical information about noise bursts. This is a useful and important feature of the excess power method, even though we focus almost exclusively on signal detection in this paper. Second, since the method cannot distinguish between noise bursts and signals in any one detector, it will be essential to use multiple-detector versions of the power statistic for actual signal detections. In Sec. V we derive the optimal multi-detector generalization of the excess power statistic under the assumption of Gaussian noise. It will be important in the future to generalize this analysis to allow for (uncorrelated) non-Gaussian noise components in individual detectors.

The layout of this paper is as follows: in Sec. I B we begin with an overview of the filter and some of its properties. This is done with an eye toward implementation, so that readers whose primary interest is in applying the filter need not concern themselves with mathematical aspects of the statistical theory of receivers. Subsequently we discuss properties of the excess power statistic in Sec. II, its derivation from a Bayesian framework in Sec. III, an efficient implementation of the statistic in Sec. IV, and the generalization of the power statistic to multiple detectors in Sec. V.

### B Overview

The output $`h(t)`$ of the gravitational wave detector is sampled at a finite rate $`1/\mathrm{\Delta }t`$ to produce a time series $`h_j=h(j\mathrm{\Delta }t)`$, where $`j=0,1,2\mathrm{}`$. This output can be written as

$$h_j=n_j+s_j$$
(1)
where $`n_j`$ is the detector noise and $`s_j`$ is a (possibly absent) signal. For most of this paper we assume that the noise is stationary and Gaussian. Under these assumptions, the components of the Fourier transform of a segment containing $`N`$ samples of noise,

$$\stackrel{~}{n}_k=\underset{j=0}{\overset{N1}{}}n_je^{2\pi ijk/N},$$
(2)
can be taken to be independent (we discuss the extent to which this is true in Sec. IV). The one-sided power spectrum $`S_k`$ of the noise is defined by:

$$\stackrel{~}{n}_k\stackrel{~}{n}_k^{}=\frac{1}{2}S_k.$$
(3)
Here, $``$ indicates an average over the noise distribution and $``$ denotes complex conjugation.

Consider a situation where all possible signals have a fixed time duration $`\delta t`$ and are band-limited to a frequency band $`[f_s,f_s+\delta f]`$, but no other information is known about them. Then, as we show in Sec. III the *optimal* statistic for detection of this class of signals is the *excess power statistic*

$$=4\underset{k_1k<k_2}{}|\stackrel{~}{h}_k|^2/S_k.$$
(4)
The sum in Eq. (4) is over the positive frequency components $`k_1k<k_2`$ that define the desired frequency band. The number of frequency components being summed is therefore equal to the time-frequency volume $`V=\delta t\delta f=k_2k_1`$.

In practice, we search over every time-frequency window that is consistent with a *range* of possible time durations and bandwidths. The number of such windows per start time is

$$N_{\text{windows}}=\frac{1}{2}N_{\text{channels}}(N_{\text{channels}}+1)(N_{\text{max}}N_{\text{min}}+1)$$
(5)
where $`N_{\text{max}}=\delta t_{\text{max}}/\mathrm{\Delta }t`$ is the number of samples in the longest expected signal duration $`\delta t_{\text{max}}`$, and $`N_{\text{min}}=\delta t_{\text{min}}/\mathrm{\Delta }t`$. Here $`N_{\text{channels}}=\delta f_{\text{max}}/\delta f_{\text{min}}`$ is the ratio of the largest bandwidth searched over $`\delta f_{\text{max}}`$ to the shortest bandwidth $`\delta f_{\text{min}}`$. One strategy for this search was outlined above. A flowchart for this algorithm is presented in Fig. 1. We call this algorithm the *short FFT method*. We have also considered a second algorithm which we call the *long FFT method*, and its flowchart is shown in Fig. 2.

Under the conditions that (i) the number $`M`$ of time domain data points being filtered is large and (ii) one searches over many different time-frequency volumes, the long FFT method is computationally more efficient than the short FFT method: the long FFT eliminates the redundancy of Fourier transforming data more than once when it falls into overlapping time-frequency volumes. In Sec. IV we estimate the computational costs of the two methods. We find that the short FFT method requires

$$C_{\text{short}}\frac{V^2}{2\alpha _{\text{max}}}\left(\frac{3\mathrm{log}_2V}{\alpha _{\text{max}}}+\frac{V}{3}\right)$$
(6)
floating-point operations per start time where $`V=\delta t_{\text{max}}\delta f_{\text{max}}`$ is the maximum time-frequency volume and $`\alpha _{\text{max}}=\delta f_{\text{max}}\mathrm{\Delta }t`$ is the maximum dimensionless bandwidth. (The frequency band from DC to Nyquist corresponds to $`\alpha _{\text{max}}=\frac{1}{2}`$.) The long FFT method requires only

$$C_{\text{long}}\alpha _{\text{max}}^1V^2\mathrm{ln}V$$
(7)
floating-point operations per start time; the computational improvement is a factor of

$$\frac{C_{\text{short}}}{C_{\text{long}}}\frac{3}{2\mathrm{ln}2}\frac{1}{\alpha _{\text{max}}}+\frac{1}{6}\frac{V}{\mathrm{ln}V}.$$
(8)
The first term shows that there is at least a factor of $`4`$ to be gained by the long FFT method; in addition to this, the computational gain increases with the total time-frequency volume to be searched. For $`V=100`$, the value of the second term is also $`4`$.

Having computed the total power for the various time-frequency windows of interest, one must decide which, if any, of those windows might contain a signal. A signal increases the expected power in a window, so we seek windows containing a statistically significant excess of power. If one knows the distribution of the noise data, one can derive the distribution for the noise power and thus set detection thresholds on the power.

The power $``$ is distributed as a $`\chi ^2`$ distribution with $`2V`$ degrees of freedom in the absence of a signal. When a signal is present, $``$ is distributed as a non-central $`\chi ^2`$ distribution with $`2V`$ degrees of freedom , where the non-central parameter is the signal power $`A^2`$:

$$A^2=4\underset{k_1k<k_2}{}|\stackrel{~}{s}_k|^2/S_k.$$
(9)
The quantity $`A`$ also represents the signal-to-noise ratio that one would expect to achieve if a matched filter were used to detect the signal.

In Fig. 3 we show the central and non-central $`\chi ^2`$ distributions for several choices of parameters. It is straightforward to use the $`\chi ^2`$ distributions plotted in Fig. 3 (a) to set a frequentist threshold for the excess power statistic so that a desired false alarm probability is achieved; then the false dismissal probability can be computed as a function of signal amplitude using the non-central $`\chi ^2`$ distributions plotted in Fig. 3 (b). Alternatively, on can fix both false alarm and false dismissal probabilities, and then use the $`\chi ^2`$ and non-central $`\chi ^2`$ distributions to determine the expected signal-to-noise ratio ($`A`$) of the signal which achieves these probabilities. A curve demonstrating this for a false alarm probability of $`10^9`$ and and false dismissal probability of $`0.01`$ is show in Fig. 3 (c).

A derivation of the statistical properties of the excess power statistic (both with and without a signal) can be found in Sec. II. In Sec. III we show that the method is optimal and unique under suitable assumptions using a Bayesian analysis.

The excess power statistic requires minimal information about the signals to be detected, making it a useful statistic for poorly modeled sources. Nevertheless, it is useful to compare the detection efficiency of the method to that of matched filtering. We characterize the excess power filter by the time-frequency volume $`V`$ and the matched filter bank by the number $`𝒩_{\text{eff}}`$ of effectively independent filters. For given false alarm and false dismissal probabilities, we denote the amplitude of the weakest signals detectable using the excess power filter or a bank of matched filters by $`A_{\text{min}}^{\text{EP}}`$ and $`A_{\text{min}}^{\text{MF}}`$ respectively. The relative effectiveness $`\eta `$ of the two search methods is given by the ratio of these amplitudes: $`\eta =A_{\text{min}}^{\text{EP}}/A_{\text{min}}^{\text{MF}}`$. In other words, the excess power statistic can detect a source at a fraction $`\eta `$ of the distance to which a bank of matched filters can. The relative effectiveness $`\eta `$ is plotted as a function of the time-frequency volume $`V`$ and the number of templates $`𝒩_{\text{eff}}`$ is shown in Fig. 4. This figure shows that for time frequency volumes $`V`$ less than $`100`$, and for all values of the effective number of templates $`𝒩_{\text{eff}}`$, the relative effectiveness $`\eta `$ is greater than $`1/2`$.

## II The search method in a Single Interferometer

In this section we define the search method in the context of a single interferometer, and derive its operating characteristics from the frequentist statistical framework.

### A Definition of method for a single time-frequency window

Consider stretches of discretely sampled detector data $`𝐡=\{h_0,h_1,\mathrm{},h_{N1}\}`$ consisting of $`N`$ data points. We will denote by $`𝒱`$ the $`N`$-dimensional vector space of all such data stretches. We assume that the detector output consists of a stationary, zero-mean, Gaussian noise component $`n_j`$, plus a possible signal $`s_j`$, so that $`h_j=n_j+s_j`$. Under these assumptions, the statistical properties of the noise are characterized by the $`N\times N`$ correlation matrix

$$R_{ij}n_in_j=C_n(|ij|\mathrm{\Delta }t).$$
(10)
Here $`C_n(t)`$ is the correlation function of the noise and $`\mathrm{\Delta }t`$ is the sampling time. This correlation matrix determines a natural inner product on $`𝒱`$ given by

$$(𝐚,𝐛)=\underset{i,j=0}{\overset{N1}{}}a_iQ_{ij}b_j.$$
(11)
where $`𝐐=𝐑^1`$.

We now discuss the notion of time-frequency projections. Consider the time-frequency window

$$𝒯=\{t_s,\delta t,f_s,\delta f\}$$
(12)
defined by the frequency interval $`[f_s,f_s+\delta f]`$, where $`f_s`$ is a starting frequency and $`\delta f`$ is a bandwidth, and the time interval $`[t_s,t_s+\delta t]`$, where $`t_s`$ is a starting time and $`\delta t`$ is a duration. Suppose that we want to focus attention on that portion of the data that lies inside the time-frequency window $`𝒯`$, to the extent that this is meaningful. One obvious thing to do is to truncate the data in the time domain, perform a discrete Fourier transform into the frequency domain, and then throw away the data points outside the frequency band of interest. One then obtains the quantities

$$\stackrel{~}{H}_K=\underset{J=0}{\overset{N_t1}{}}e^{2\pi iJK/N_t}h_{j+J},$$
(13)
where $`j=t_s/\mathrm{\Delta }t`$, $`N_t=\delta t/\mathrm{\Delta }t`$, and $`K`$ runs over the range $`f_s\delta tK(f_s+\delta f)\delta t`$. We denote by $`𝒲_𝒯`$ the vector space of the projected data $`\stackrel{~}{H}_K`$. The dimension of this vector space over the real numbers is

$$\mathrm{dim}𝒲_𝒯=2\delta t\delta f=2V_𝒯,$$
(14)
where $`V_𝒯\delta t\delta f`$ is the time-frequency volume of the time-frequency window $`𝒯`$.

Of course, there are many other methods of attempting to pick out the portion of the data in the time-frequency window $`𝒯`$.For example, one could FFT the entire data segment, truncate it in the frequency domain, FFT back to the time domain, and then truncate it again in the time domain. The lack of a preferred unique method is due to the uncertainty principle. In many circumstances the differences between different reasonable choices will be relatively unimportant. For the remainder of this section we will assume that we have picked some reasonable projection method.The choice of a projection method corresponds mathematically to the choice of a $`2V_𝒯`$-dimensional subspace of the dual space $`𝒱^{}`$ of $`𝒱`$. When one specifies in addition the detector noise spectrum, the projection method determines a $`2V_𝒯`$ dimensional subspace of $`𝒱`$. We can write the projected data in general as

$$\overline{h}_J=\underset{j=0}{\overset{N1}{}}A_J^jh_j,$$
(15)
where $`A_J^j`$ is a real $`2V_𝒯\times N`$ matrix, the quantities $`\overline{h}_J`$ are real, and $`J`$ runs over $`0J2V_𝒯1`$.

We define the power statistic associated with $`𝒯`$ and with a choice of projection method to be

$$_𝒯(𝐡)\underset{I,J=0}{\overset{2V_𝒯1}{}}Q_{IJ}\overline{h}_I\overline{h}_J,$$
(16)
where $`_JQ_{IJ}R_{JK}=\delta _{IK}`$ and

$$R_{JK}=\overline{h}_J\overline{h}_K=\underset{j,k}{}A_J^jA_K^kR_{jk}$$
(17)
is the correlation matrix of the projected data. The quantity $``$ is, roughly speaking, just the total power in the data stream within the given time-frequency window, where power is not the physical power but is measured relative to the detector noise (i.e., it is the conventional power of the pre-whitened data stream).

The statistic can also be described geometrically as follows. The linear mapping defined by Eq. (15) has a kernel $`\left\{h_i\right|\overline{h}_I=0\text{ for all }I\}`$. The set of all vectors $`𝐡`$ perpendicular to all elements of this kernel with respect to the inner product (11) form a subspace $`𝒱_𝒯`$ of $`𝒱`$ which can be naturally identified with $`𝒲_𝒯`$. Any element $`𝐡`$ in $`𝒱`$ can be decomposed as

$$𝐡=𝐡_{}+𝐡_{},$$
(18)
where $`𝐡_{}`$ lies in $`𝒱_𝒯`$ and $`𝐡_{}`$ is perpendicular to all elements of $`𝒱_𝒯`$. The statistic (16) is the squared norm of the parallel component:

$$_𝒯(𝐡)=(𝐡_{},𝐡_{}).$$
(19)
For the simple time-frequency truncation method (13) discussed above, one can obtain a simple approximate formula for the statistic. The correlation matrix of the quantities $`\stackrel{~}{H}_J`$ of Eq. (13) is given by, to a first, crude approximation,

$`\stackrel{~}{H}_J\stackrel{~}{H}_K`$ $`=`$ $`0,`$ (20)
$`\stackrel{~}{H}_J\stackrel{~}{H}_K^{}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\delta _{JK}S_K,`$ (21)
where $`0J,KN_t/2`$,

$$S_K=\delta tS_h(K/\delta t)/(\mathrm{\Delta }t)^2$$
(22)
and $`S_h(f)`$ is the conventional one-sided power spectral density of the detector noise. The expressions (20) and (21) are accurate only when $`f_s\delta t1`$ and when $`S_h(f)`$ does not vary substantially on scales $`1/\delta t`$; more accurate expressions can be computed if desired from Eqs. (10) and (13). The definition (16) now yields

$$_𝒯(𝐡)4\underset{K=f_s\delta t}{\overset{(f_s+\delta f)\delta t}{}}|\stackrel{~}{H}_K|^2/S_K,$$
(23)
cf., Eq. (4) of the Introduction. We show in Sec. IV below that the expression (23) is an adequate approximation to $`_𝒯(𝐡)`$ for most purposes.

The search method consists of searching over time-frequency windows $`𝒯`$, and selecting as possible events only those windows $`𝒯`$ for which $`_𝒯`$ exceeds a suitable threshold $`^{}`$. We discuss further how to search over time-frequency windows in Sec. IV below. For the remainder of this section we assume that the time-frequency window $`𝒯`$ is fixed and known, and discuss the performance of the statistic.

### B Operating characteristics of the statistic

When a signal is not present in the data stream, the statistic $`_𝒯(𝐡)`$ is the sum of the squares of $`2V`$ independent, zero-mean, unit-variance Gaussian random variables.§§§To see this, note that there is a basis of $`𝒲_𝒯`$ in which the correlation matrix $`R_{JK}`$ is equal to the $`2V\times 2V`$ identity matrix, and use Eqs. (16) and (17). Thus $``$ follows a $`\chi ^2`$ distribution with $`2V`$ degrees of freedom; the upper-tail cumulative probability is

$$Q_0(^{})=P(>^{})=\frac{\mathrm{\Gamma }(V,^{}/2)}{\mathrm{\Gamma }(V)},$$
(24)
where $`\mathrm{\Gamma }(a,x)=_x^{\mathrm{}}e^tt^{a1}𝑑t`$ is the incomplete Gamma function. The quantity $`Q_0(^{})`$ is the *false alarm* probability for the detection threshold $`^{}`$. The distribution (24) is plotted in Fig. 5 for several values of $`V`$. An approximate expression for $`Q_0`$ in the regime $`Q_01`$ is

$`Q_0(^{})`$ $`=`$ $`\sqrt{{\displaystyle \frac{2V}{\pi }}}{\displaystyle \frac{1}{A_{}^2}}\left(1+{\displaystyle \frac{A_{}^2}{2V}}\right)^Ve^{A_{}^2/2}`$ (26)
$`\times \left[1+O\left({\displaystyle \frac{1}{V}}\right)+O\left({\displaystyle \frac{1}{A_{}^2}}\right)+O\left({\displaystyle \frac{V}{A_{}^4}}\right)\right],`$
where $`A_{}^2^{}2V`$.

We next consider the case when a signal is present, so that $`𝐡=𝐧+𝐬`$. The formula for the statistic $``$ given by Eqs. (15) and (16) becomes

$$=\underset{I,J=0}{\overset{2V1}{}}Q_{IJ}(\overline{n}_I+\overline{s}_I)(\overline{n}_J+\overline{s}_J),$$
(27)
where $`\overline{n}_I=_jA_I^jn_j`$ is the projected noise and $`\overline{s}_I=_jA_I^js_j`$ is the projected signal. We define the amplitude $`A`$ of the signal byNote that $`A`$ is the signal-to-noise ratio that would be obtained by matched filtering if prior knowledge of the waveform shape allowed one to perform matched filtering, and if the signal $`𝐬`$ were confined to the time-frequency window $`𝒯`$.

$$A^2=(𝐬_{},𝐬_{})=\underset{I,J=0}{\overset{2V1}{}}Q_{IJ}\overline{s}_I\overline{s}_J,$$
(28)
where we use the notation of Eq. (18). The expression (27) can be simplified by (i) choosing the basis of $`𝒲_𝒯`$ so that $`Q_{IJ}=\delta _{IJ}`$ (which is roughly equivalent to whitening the detector output for a large class of time-frequency windows) and (ii) further specializing the choice of basis so that the signal vector is $`(\overline{s}_0,\overline{s}_1,\mathrm{},\overline{s}_{2V1})=(A,0,\mathrm{},0)`$. The result is

$$=(\overline{n}_0+A)^2+\underset{I=1}{\overset{2V1}{}}\overline{n}_I^2$$
(29)
where $`\overline{n}_0,\overline{n}_1,\mathrm{},\overline{n}_{2V1}`$ are independent Gaussian random variables with zero mean and unit variance.

From Eq. (29) one can compute the moment generating function for the random variable $``$. The result is

$$e^t=\frac{\mathrm{exp}[A^2t/(12t)]}{(12t)^V}.$$
(30)
The probability distribution for $``$ can be now obtained by taking the inverse Laplace transform, which is accomplished by expanding the argument of the exponential in Eq. (30) as a power series in $`A`$. The result is a weighted sum of $`\chi ^2`$ probability distribution functions:

$$p(|A,V)=\underset{n=0}{\overset{\mathrm{}}{}}\frac{e^{A^2/2}(A^2/2)^n}{n!}\frac{e^{/2}(/2)^{n+V1}}{\mathrm{\Gamma }(n+V)}.$$
(31)
This is the non-central $`\chi ^2`$ probability distribution with non-centrality parameter $`A^2`$ discussed in Sec. 26.4 of Ref. . A closed form expression for the probability distribution is

$$p(|A,V)=\frac{1}{2}e^{(+A^2)/2}(^{1/2}/A)^{V1}I_{V1}(A^{1/2}),$$
(32)
where $`I_n(x)`$ is the modified Bessel function of the first kind of order $`n`$.

The upper-tail cumulative probability distribution for $``$

$$Q_A(^{},A,V)=P(>^{}|A,V)=_{^{}}^{\mathrm{}}p(|A,V)𝑑$$
(33)
is the *true detection* probability for a given threshold $`^{}`$ and a given signal amplitude $`A`$. Figure 6 shows this true detection probability $`Q_A`$, expressed as a function of signal strength $`A`$ and false alarm probability $`Q_0`$ \[via Eq. (24)\], evaluated at $`Q_0=0.01`$, for several different values of the time-frequency volume $`V`$.

Some qualitative insight into the detectability of a signal can be obtained from the first two moments of the distribution for $``$. The expected (mean) value is $`=2V+A^2`$, while the variance is $`\mathrm{Var}=^2^2=4V+4A^2`$. For large values of $`V`$ the probability distributions are nearly Gaussian within a few sigma of the expected value, so we can imagine setting the threshold $`^{}`$ to be a few times $`\sqrt{4V}`$ above the mean noise level $`2V`$ in order to achieve the required false alarm probability. Thus, a signal will be detectable when $`2V>(\text{a few})\times \sqrt{4V}`$. In other words, the signal power $`A^2`$ can be small compared to the mean noise power $`2V`$ and still be detectable; it need only be comparable to the much smaller fluctuations $`\sqrt{4V}`$ in the noise power.

Once one specifies the time-frequency volume $`V`$ and desired values of the false alarm probability $`Q_0`$ and true detection probabilities $`Q_A`$, there is a minimum signal amplitude $`A_{\text{min}}`$ that can be detected with the excess power method. To compute this amplitude, one first inverts Eq. (24) to obtain the required threshold $`^{}`$ as a function of $`Q_0`$ and $`V`$: $`^{}=^{}(Q_0,V)`$. Second, one inverts Eq. (33) to obtain the amplitude $`A`$ as a function of $`^{}`$, $`Q_A`$ and $`V`$: $`A=A(^{},Q_A,V)`$. The minimum signal amplitude is then given by

$$A_{\text{min}}(Q_0,Q_A,V)=A[^{}(Q_0,V),Q_A,V].$$
(34)
This quantity is plotted as a function of the time frequency volume in Fig. 3 (c), for various values of $`Q_0`$ and for $`Q_A=0.99`$.

### C Comparison of performance to that of matched filtering.

The performance of the excess power statistic should be compared with that of a matched filtering search for the same class of signals. Of course, matched filtering searches will not be possible for the classes of signals we are interested in (for example supernovae) due to the lack of theoretical templates; nevertheless the comparison is useful as benchmark of the excess power method.

We start by discussing how the performance of matched filtering depends on the set of signals being searched for. Suppose that a given class of signals have a known duration and frequency band, so that they all lie inside a fixed time-frequency window $`𝒯`$ with time-frequency volume $`V`$. Let $`A_{}(Q_0)`$ be the signal-to-noise ratio threshold for the matched filtering search necessary to give a false alarm probability of $`Q_0`$ for each starting time. Now if the bank of matched filters consisted of $`𝒩`$ statistically independent filters, then the threshold $`A_{}`$ would be given by the formula

$$\mathrm{erfc}(A_{}/\sqrt{2})=1(1Q_0)^{1/𝒩}Q_0/𝒩.$$
(35)
In the more realistic case where the templates are not all statistically independent, the formula (35) motivates us to define an effective number of independent templates $`𝒩_{\text{eff}}=𝒩_{\text{eff}}(Q_0)`$ by the relationIn Ref.  the quantities $`A_{}`$, $`Q_0`$ and $`𝒩_{\text{eff}}`$ were denoted $`\rho _{}`$, $`ϵ/𝒩_{\text{start-times}}`$, and $`𝒩_{\text{shapes}}`$, respectively.

$$\mathrm{erfc}[A_{}(Q_0)/\sqrt{2}]=Q_0/𝒩_{\text{eff}}(Q_0).$$
(36)
This effective number of templates depends on the false alarm probability $`Q_0`$, or, equivalently, on the detection threshold $`A_{}`$ via the relation $`A_{}=A_{}(Q_0)`$. For a given class of signals (e.g., inspiralling binaries), it should be possible to estimate $`𝒩_{\text{eff}}`$ by a modification of the method of Ref.  wherein one does not eliminate the intrinsic parameters or the signal amplitude and one chooses a minimal match of order $`0.3`$ say to give approximately statistically independent templates. We suspect that $`𝒩_{\text{eff}}`$ will not differ too significantly from the actual number of templates used in a search.\**\**\**The grid of search templates used will be determined by having a minimal match  of $`0.97`$ say instead of $`0.3`$, which tends to make the actual number of templates larger than $`𝒩_{\text{eff}}`$. On the other hand, a template grid needs only 2 templates to cover all possible signal amplitudes and phases (when the other parameters are fixed), whereas the number of statistically independent templates that can be generated by varying the amplitude and phase can be much greater than 2 for small $`Q_0`$. In any case, for a given matched filtering search, the detection threshold and the resulting effective number of independent filters will be determined by Monte Carlo simulations.

An illustrative special case of matched filtering is when the signal manifold is a linear subspace $`𝒮`$ of the space of all possible signals. In this case the maximum over all templates of the signal-to-noise ratio squared is simply $`(𝐡_{},𝐡_{})`$, where $`𝐡_{}`$ is the perpendicular projection of the detector output $`𝐡`$ into $`𝒮`$ . Since this quantity has a $`\chi ^2`$ distribution, it follows from the approximate formula (26) and from the definition (36) that for this special case we have

$`𝒩_{\text{eff}}(A_{})`$ $`=`$ $`{\displaystyle \frac{\sqrt{2d}}{A_{}}}\left[1+{\displaystyle \frac{A_{}^2}{d}}\right]^{d/2}`$ (37)
$`\times `$ $`\left[1+O\left({\displaystyle \frac{1}{d}}\right)+O\left({\displaystyle \frac{1}{A_{}^2}}\right)+O\left({\displaystyle \frac{d}{A_{}^4}}\right)\right],`$ (38)
where $`d=\mathrm{dim}(𝒮)`$ is the number of signal parameters. In Ref.  this relation was used to define an effective dimension for any space of signals, and that effective dimension was used instead of $`𝒩_{\text{eff}}(Q_0)`$ to characterize the signal space. Here however we will instead parameterize our comparisons directly in terms of $`𝒩_{\text{eff}}`$.

We also note that for this special case of a linear signal subspace, we have

$$𝒩_{\text{eff}}(A_{})2^{I_s(A_{})},$$
(39)
where $`I_s(A_{})`$ is the number of bits of information about the source carried by a detected signal with signal-to-noise ratio $`A_{}`$, as defined in Ref. . We conjecture that this relation might be approximately valid for general signal manifolds.

We now turn to the relative performance of the excess power and matched filtering search methods. As explained in the previous subsection, once we specify the false alarm probability $`Q_0`$ and true detection probability $`Q_A`$ we can compute the minimum signal amplitude $`A_{\text{min}}`$ necessary for detection via the excess power method, as a function of $`Q_0`$, $`Q_A`$ and $`V`$. Let us denote that value here as $`A_{\text{min}}^{\text{EP}}`$:

$$A_{\text{min}}^{\text{EP}}=A_{\text{min}}(Q_0,Q_A,V).$$
(40)
Similarly, we can compute the minimum amplitude necessary for detection with matched filtering with $`𝒩_{\text{eff}}`$ independent templates; the result is

$$A_{\text{min}}^{\text{MF}}=A_{\text{min}}(Q_0/𝒩_{\text{eff}},Q_A,1/2).$$
(41)
In other words, one simply uses the same formula with $`Q_0`$ replaced by $`Q_0/𝒩_{\text{eff}}`$, and with the number of degrees of freedom $`2V`$ being unity. We define the relative effectiveness $`\eta `$ of the excess power method relative to the bank of filters by

$$\eta (Q_0,Q_A,𝒩_{\text{eff}},V)=A_{\text{min}}^{\text{EP}}/A_{\text{min}}^{\text{MF}}.$$
(42)
The factor by which the event rate for the excess power method is smaller than that for the bank of matched filters is $`\eta ^3`$. The relative effectiveness is plotted in Fig. 4 as a function of $`V`$ and $`𝒩_{\text{eff}}`$. For $`V100`$ we see that $`\eta 0.6`$ always, showing that the excess power method performs almost as well as matched filtering.

When the time-frequency window $`𝒯`$ is not known in advance, one must search over time frequency windows. This reduces the efficiency of the excess power method compared to matched filtering. An approximate parameterization of this reduction can be obtained by replacing in Eq. (40) the false alarm probability $`Q_0`$ by $`Q_0/𝒩_\text{w}`$, where $`𝒩_\text{w}`$ is the number of statistically independent time frequency windows searched over per start time. We have performed Monte Carlo simulations with white Gaussian noise which suggest that $`𝒩_\text{w}100V_{\text{max}}`$, where $`V_{\text{max}}`$ is the largest time-frequency volume searched over. The resulting change in the relative efficiency $`\eta `$ is not very large.

Further insight into the relation between the excess power and matched filtering methods can be obtained as follows. First, an approximate formula for the function (34) is obtained by approximating the distribution of $``$ to be a Gaussian:

$`A_{\text{min}}(Q_0,Q_A,V)^2=A_{}(Q_0,V)^2`$ (43)
$`+2\sqrt{2}\mathrm{erf}^1(2Q_A1)\sqrt{V+A_{}(Q_0,V)^2},`$ (44)
where $`A_{}(Q_0,V)`$ is obtained by inverting Eq. (26). This formula is typically accurate to a few percent for $`0.5Q_A0.99`$. When $`Q_01Q_A`$ (which will typically be the case), we can neglect the second term in Eq. (44) in comparison with the first, so that

$$A_{\text{min}}^{\text{EP}}(Q_0,Q_A,V)A_{}(Q_0,V).$$
(45)
Second, the quantity $`A_{\text{min}}^{\text{MF}}(Q_0,Q_A)`$ will similarly depend weakly on $`Q_A`$ and will be well approximated by the quantity $`A_{}(Q_0)`$ obtained from Eq. (36). Combining these approximations together with Eq. (26), we see that the excess power method is equivalent (in terms of detection thresholds) to matched filtering with an effective number of templates of

$$𝒩_{\text{eff}}=\left(1+\frac{A_{}^2}{2V}\right)^V\frac{\sqrt{V}}{A_{}}\left(1+\frac{A_{}^2}{2V}\right)^V.$$
(46)
The quantity (46) was shown in Ref.  to be the total number of distinguishable signals within the given time-frequency window with signal-to-noise $`A_{}`$. In other words, it is the maximum possible value for the effective number of independent templates $`𝒩_{\text{eff}}`$ for *any* manifold of signals inside the time-frequency window $`V_𝒯`$, as a function of the time-frequency volume $`V`$. Hence, we can understand the excess power method as a limiting case of matched filtering: the case where the manifold of signals becomes so large (perhaps curving back and intersecting itself) that (when smeared out by the noise) it effectively fills up the entire space $`V_𝒯`$ of signals within the given time-frequency window. The equivalence can also be seen from the fact, noted above, that the excess power statistic coincides with the matched filtering statistic when dependence of the signal on its parameters is linear and the number of parameters coincides with the dimension $`2V`$ of $`W_𝒯`$ .

## III BAYESIAN ANALYSIS OF SIGNAL DETECTION

In this section we show how our proposed search method arises naturally from an analysis of the detection of signals from a Bayesian point of view. Section III A defines the class of signals under consideration in terms of a prior probability density function (PDF). Section III B derives the excess power statistic, and Sec. III C compares detection criteria based on a false alarm rate to criteria based on the probability that a signal is present in the data.

### A The space of signals

The signals of interest (e.g., black hole mergers) are poorly understood. We characterize our knowledge in terms of a prior PDF $`p(𝐬)`$ for signals $`𝐬`$ in the vector space $`𝒱`$. In this subsection we explain how to encode knowledge of the expected bandwidth and duration of the signals in the PDF.

Suppose we know that the signal $`𝐬`$ lies approximately within some time-frequency window $`𝒯=[t_s,t_s+\delta t]\times [f_s,f_s+\delta f]`$, but that nothing else is known about the signal. Then we know that $`𝐬`$ belongs to a subspace $`𝒱_𝒯`$ of $`𝒱`$. Of course, there are several slightly inequivalent choices of such a subspace, as discussed in Sec. II A above, but we will assume that these differences are unimportant, and pick one choice of $`𝒱_𝒯`$.

For any vector $`𝐡`$ in $`𝒱`$, we can write $`𝐡=𝐡_{}+𝐡_{}`$, where $`𝐡_{}`$ is the projection of $`𝐡`$ into $`𝒱_𝒯`$ and $`𝐡_{}`$ is perpendicular to $`𝒱_𝒯`$. We assume the following form for the prior PDF $`p(𝐬|𝒯)`$ given the time-frequency window $`𝒯`$:

$`p(𝐬|𝒯)d^N𝐬`$ $`=`$ $`\delta ^{(N2V)}(𝐬_{})d^{(N2V)}𝐬_{}`$ (48)
$`\times p_1(A)d^{2V}𝐬_{},`$
where $`A^2=(𝐬_{},𝐬_{})`$ and $`N`$ is the dimensions of $`𝒱`$. Here the first factor consisting of the ($`N2V`$)-dimensional $`\delta `$ function restricts $`𝐬`$ to lie in $`𝒱_𝒯`$. The second factor depends only on the magnitude $`A`$ of $`𝐬_{}`$, which means that we assume all directions in the vector space $`𝒱_𝒯`$ are equally likely when one measures lengths and angles with the inner product (11).††††††It would be more realistic to make this assumption with respect to an inner product on $`𝒱`$ whose definition did not depend on the noise spectrum, but if the noise spectrum does not vary too rapidly within the bandwidth of interest, the distinction is not too important and our assumption will be fairly realistic. We can rewrite the prior PDF (48) as

$`p(𝐬|𝒯)d^N𝐬`$ $`=`$ $`\delta ^{(N2V)}(𝐬_{})d^{(N2V)}𝐬_{}`$ (50)
$`\times {\displaystyle \frac{\mathrm{\Gamma }(V)}{2\pi ^V}}d^{(2V1)}\mathrm{\Omega }_{}p(A)dA`$
where $`d^{(2V1)}\mathrm{\Omega }_{}`$ is the $`(2V1)`$-dimensional element of solid angle and where $`p(A)dA`$ is the probability that the signal amplitude lies between $`A`$ and $`A+dA`$. We discuss the choice of $`p(A)`$ in Sec. III C below.

So far we have assumed that the time interval $`[t_s,t_s+\delta t]`$ and frequency interval $`[f_s,f_s+\delta f]`$ are known. In a real search, however, one must account for ignorance of these parameters. An appropriate prior which does this is

$$p(𝐬)d^N𝐬=p(𝐬|𝒯)d^N𝐬p_𝒯(𝒯)d^4𝒯,$$
(51)
where $`p_𝒯(𝒯)=p_𝒯(t_s,\delta t,f_s,\delta f)`$ is a prior PDF on the time-frequency window parameters. The PDF $`p_𝒯(𝒯)`$ should be uniform in $`t_s`$, but its dependence on the parameters $`\delta t`$, $`f_s`$ and $`\delta f`$ will depend on the class of sources under consideration. We will see below that our analysis depends only weakly on the choice of PDF $`p_𝒯(𝒯)`$, as long as it is a slowly varying function of its parameters.

### B Derivation of the search method

In the Bayesian approach to signal detection there is a unique and optimal method to search the data stream for signals if the statistical properties of the detector noise and the prior probability distribution for signals are known. One computes the probability $`p_s(𝐡)`$ that some signal $`𝐬`$ is present in the measured data $`𝐡`$. The signal-detection criterion is that the probability $`p_s(𝐡)`$ exceeds some threshold value. This is the starting point for our analysis; more details can be found in Wainstein and Zubakov  or Finn and Chernoff .

The prior PDF given Eqs. (50) and (51) describes our state of knowledge about the signals to be searched for. Now let $`p_{s0}`$ denote the a priori probability that gravitational waves exist (or that our signal model is correct), for which an appropriate value for the first searches might be $`p_{s0}=1/2`$. The signal PDF (51) then gets modified to

$$(1p_{s0})\delta ^N(𝐬)d^N𝐬+p_{s0}p(𝐬)d^N𝐬.$$
(52)
It then follows that the posterior probability $`p_s(𝐡)`$ that a signal is present in the data $`𝐡`$ is given by

$$\frac{p_s(𝐡)}{1p_s(𝐡)}=\mathrm{\Lambda }(𝐡)\frac{p_{s0}}{1p_{s0}},$$
(53)
where the likelihood function $`\mathrm{\Lambda }(𝐡)`$ is

$`\mathrm{\Lambda }(𝐡)`$ $`=`$ $`{\displaystyle \mathrm{\Lambda }(𝐡;𝐬)p(𝐬)d^N𝐬}`$ (54)
$`=`$ $`{\displaystyle \mathrm{\Lambda }(𝐡;𝐬)p(𝐬|𝒯)p(𝒯)d^N𝐬d^4𝒯}`$ (55)
and

$$\mathrm{\Lambda }(𝐡;𝐬)=\frac{p(𝐡|𝐬)}{p(𝐡|𝐬=\mathrm{𝟎})}.$$
(56)
In Eq. (56) the quantity $`p(𝐡|𝐬)`$ is the probability of measuring the time series $`𝐡`$ when the signal $`𝐬`$ is present, and $`p(𝐡|𝐬=\mathrm{𝟎})`$ is the corresponding probability when no signal is present. For stationary Gaussian noise the likelihood ratio $`\mathrm{\Lambda }(𝐡;𝐬)`$ is .

$$\mathrm{\Lambda }(𝐡;𝐬)=\mathrm{exp}[(𝐡,𝐬)(𝐬,𝐬)/2].$$
(57)
Equation (53) shows that the probability $`p_s(𝐡)`$ increases monotonically with increasing $`\mathrm{\Lambda }(𝐡)`$. Consequently, thresholding on $`\mathrm{\Lambda }(𝐡)`$ to detect signals is equivalent to thresholding on the probability that a signal is present in the data stream. This is also the optimal signal detection strategy in the Neyman-Pearson sense of maximizing the detection probability for a fixed false alarm probability .

The integral (54) includes a integral over all possible time-frequency windows $`𝒯`$, which can be approximated as a sum:

$$\mathrm{\Lambda }(𝐡)\frac{1}{𝒩_{\text{windows}}}\underset{𝒯}{}\mathrm{\Lambda }(𝐡;𝐬)p(𝐬|𝒯)d^N𝐬.$$
(58)
Here $`𝒩_{\text{windows}}`$ is the number of grid points in a grid on the four dimensional space of time-frequency windows used to approximate the integral. Now if a signal is present, the summand in Eq. (59) will be a sharply peaked function of $`𝒯`$. If the grid spacing is chosen to approximately coincide with the width of the peak (so that $`𝒩_{\text{windows}}`$ is the number of statistically independent time-frequency windows) then the sum will be dominated by the largest term, and we obtain

$$\mathrm{\Lambda }(𝐡)\frac{1}{𝒩_{\text{windows}}}\underset{𝒯}{\mathrm{max}}\mathrm{\Lambda }(𝐡;𝐬)p(𝐬|𝒯)d^N𝐬.$$
(59)
Thus it is sufficient to consider only a single time-frequency region in the remainder of this section with the understanding that the signal detection will be based on the maximum of the likelihood function over all relevant time-frequency windows. Also we can factor $`𝒩_{\text{windows}}`$ as $`𝒩_{\text{windows}}=𝒩_{\text{st}}𝒩_\text{w}`$, where $`𝒩_{\text{st}}`$ is the number of statistically independent starting times $`t_s`$ in the search, and $`𝒩_\text{w}`$ is the number of statistically independent time-frequency windows per start time.

The evaluation of the integral in Eq. (59) with the prior PDF given by Eq. (50) can be done in stages. First we integrate over the delta-function to restrict the possible signals to the vector space $`𝒱_𝒯`$. This essentially replaces $`𝐬`$ by $`𝐬_{}`$ in Eq. (57). Next we use the definition (19) of $``$ to write the inner product appearing in Eq. (57) as

$$(𝐡,𝐬_{})=A^{1/2}\mathrm{cos}\theta ,$$
(60)
where $`\theta `$ is the angle between the vectors $`𝐡_{}`$ and $`𝐬_{}`$. We then obtain the formula

$$\mathrm{\Lambda }(𝐡)=\frac{1}{𝒩_{\text{st}}𝒩_\text{w}}\mathrm{\Lambda }(𝐡;A)p(A)𝑑A$$
(61)
with

$`\mathrm{\Lambda }(𝐡;A)`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Gamma }(V)e^{A^2/2}}{\pi ^{1/2}\mathrm{\Gamma }(V1/2)}}{\displaystyle _0^\pi }e^{A^{1/2}\mathrm{cos}\theta }\mathrm{sin}^{2V2}\theta d\theta `$ (62)
$`=`$ $`\mathrm{\Gamma }(V)e^{A^2/2}(A^{1/2}/2)^{1V}I_{V1}(A^{1/2})`$ (63)
$`=`$ $`p(|A,V)/p(|A=0,V)`$ (64)
\[cf. Eq. (32)\]. Here $`I_n(x)`$ is the modified Bessel function of the first kind of order $`n`$, and maximization over time frequency windows $`𝒯`$ is understood.

The quantity (62) is a monotonically increasing function of the power $``$. Hence thresholding on $`\mathrm{\Lambda }`$ is equivalent to thresholding on $``$, and so $``$ is the optimal (in the Bayesian sense) statistic for the detection of the class of signals we have considered.

### C Bayesian thresholds

Frequentist detection thresholds $`^{}`$ are set by specifying a false alarm rate, and can be computed using Eq. (24). As is well known, if such a threshold is exceeded it does not necessarily mean that a signal is present with high probability, even for low false alarm probabilities . To determine how likely it is that a signal is actually present in the data stream requires the use of Bayesian methods.

For a Bayesian detection strategy, one sets a threshold on the posterior probability $`p_s(𝐡)`$ that a signal is present given the data. This probability is related to the likelihood function $`\mathrm{\Lambda }(𝐡)`$ by Eq. (53). In general, the integral (61) required to compute $`\mathrm{\Lambda }(𝐡)`$ must be performed numerically. Since $`\mathrm{\Lambda }(𝐡)`$ depends on the data $`𝐡`$ only through $`(𝐡)`$, one can determine a Bayesian threshold for $``$ from the value of $`p_s`$.

Consider a search characterized by $`𝒩_{\text{st}}`$ statistically independent start times and $`𝒩_\text{w}`$ statistically independent time-frequency windows. The frequentist false alarm probability $`Q_0`$ of the previous section \[Eq. (24)\] is the false alarm probability for a given start time and a given time-frequency window. Hence the false alarm probability for the entire search is

$$p_{\text{fa}}(^{})=Q_0(^{})𝒩_{\text{st}}𝒩_\text{w}.$$
(65)
It is natural, in comparing frequentist and Bayesian thresholds, to set $`p_s=1p_{\text{fa}}`$. For example, for “$`99\%`$ confidence” one would choose $`p_s=0.99=1p_{\text{fa}}`$. We emphasize that this means “$`99\%`$ confidence that events will be due to signals” for the Bayesian, while it means “$`99\%`$ confidence that there will be no false events” for the frequentist; since these are different statistical statements, the frequentist and the Bayesian will obtain different thresholds.

We first discuss approximate evaluation of Bayesian thresholds. The integral (61) can be approximately evaluated in the the regime $`V1`$ by using the Laplace approximation, if the prior PDF $`p(A)`$ does not vary too rapidly. The result is

$`\mathrm{\Lambda }(𝐡)`$ $`=`$ $`{\displaystyle \frac{\sqrt{2\pi V}}{\widehat{A}}}{\displaystyle \frac{p(\widehat{A})}{𝒩_{\text{st}}𝒩_\text{w}}}\left(1+{\displaystyle \frac{\widehat{A}^2}{2V}}\right)^V\mathrm{exp}(\widehat{A}^2/2)`$ (67)
$`\times \left[1+O\left({\displaystyle \frac{1}{V}}\right)+O\left({\displaystyle \frac{1}{\widehat{A}^2}}\right)+O\left({\displaystyle \frac{V}{\widehat{A}^4}}\right)\right],`$
where $`\widehat{A}=\widehat{A}(𝐡)`$ is defined by

$$(𝐡)=2V+\widehat{A}(𝐡)^2.$$
(68)
If we now compare Eqs. (26), (53), (65) and (67) and use $`p_{s0}=1/2`$ and $`1p_s1`$, we see that the Bayesian threshold $`\widehat{A}`$ and frequentist threshold $`A_{}`$ are related by

$$\left(1+\frac{\widehat{A}^2}{2V}\right)^Ve^{\widehat{A}^2/2}=\left(1+\frac{A_{}^2}{2V}\right)^Ve^{A_{}^2/2},$$
(69)
where the factor $``$ is

$$=\frac{\widehat{A}A_{}^2}{2Vp(\widehat{A})}.$$
(70)
Clearly the two thresholds coincide when $`=1`$. However, typically the factor $``$ can be quite significant and can cause the Bayesian and frequentist thresholds to differ substantially.

It is useful to consider a specific example. Suppose that we are searching for black hole mergers which we expect to produce short (a few ms) broad band signatures with a time-frequency volume of $`V=100`$. We want to be $`99\%`$ sure of our detection, so we set $`p_s=0.99=1p_{\text{fa}}`$. Suppose that the search duration is 1/3 of a year, so that the number of independent start time is $`𝒩_{\text{st}}=10^{10}`$ say, and that the number of statistically independent windows is $`10^4`$. For $`99\%`$ confidence that there will be no false alarms, we should choose $`Q_0=10^{16}`$, from Eq. (65). This gives from Eq. (24) a frequentist threshold of $`^{}=411.3`$ corresponding to a signal-to-noise threshold of $`A_{}=14.5`$. The corresponding Bayesian threshold depends on the specification of prior PDF for signal amplitudes $`A`$. A reasonable choice of prior is

$$p(A)=\{\begin{array}{cc}3A_c^3/A^4\hfill & AA_c\text{,}\hfill \\ 0\hfill & A<A_c\text{.}\hfill \end{array}$$
(71)
This is just the distribution that would be expected for sources distributed uniformly in time and space, except that it is cutoff in an approximate way at small $`A`$ in order to ensure correct normalization. The parameter $`A_c^3`$ is prior probability per start-time of an event being present with signal-to-noise ratio exceeding unity. Based on population estimates such Ref. , we optimistically assume a prior probability of order unity for approximately one merger event per year with $`A>1`$, which translates into $`A_c^310^{10}`$. We can now compute a Bayesian threshold by combining Eqs. (53), (61), and (62). The result is $`=539`$ corresponding a signal-to-noise ratio threshold of $`\widehat{A}=18.4`$, which is substantially higher than the frequentist value of $`14.5`$.

## IV Implementation

As discussed in Sec. II A, one will not know in advance the start time $`t_s`$, duration $`\delta t`$ and frequency band $`[f_s,f_s+\delta f]`$ of signals in a real search, and thus one must perform a search over these four parameters. One needs to compute the excess power $`_𝒯(𝐡)`$ for each possible time-frequency window, and record as possible events all of those windows for which $`_𝒯`$ is above threshold.‡‡‡‡‡‡Note that different thresholds will be required for each window $`𝒯`$, but the false alarm probability $`Q_0`$ will be the same for each window. We assume that we wish to search over all values of $`\delta t`$ in the range

$$\delta t_{\text{min}}\delta t\delta t_{\text{max}},$$
(72)
and over all $`f_s`$ and $`\delta f`$ with

$$\begin{array}{ccc}\delta f_{\text{min}}\hfill & \hfill & \delta f\delta f_{\text{max}}\hfill \\ f_{\text{min}}\hfill & \hfill & f_s\hfill \\ f_s+\delta f\hfill & \hfill & f_{\text{min}}+\delta f_{\text{max}}\hfill \end{array}\}.$$
(73)
It is clear that the computational cost can quickly grow to unreasonable proportions, so it is important to achieve an efficient implementation of the search technique.

There are (at least) two different ways to implement a search over a pre-specified set of time-frequency windows. The first uses many FFTs of data segments with durations in the range (72) as suggested by the derivation in Sec. II, and for each FFT computes $``$ for all frequency bands in the range (73). This process is then repeated for every possible start time. We call this procedure the *short FFT method*. The second method partitions the time series into long data segments each containing $`M`$ samples, and for each of these segments computes its FFT. That FFT is then partitioned into $`\delta f_{\text{max}}/\delta f_{\text{min}}`$ non-overlapping frequency bands each of width $`\delta f_{\text{min}}`$, and for each one the FFT is bandpass filtered to that frequency band and then inverse Fourier transformed. The result is $`\delta f_{\text{max}}/\delta f_{\text{min}}`$ different timeseries, which we call channels, each containing particular frequency information. The elements of these time series are then squared. One obtains in this way a time-frequency plane in which each pixel represents the total power in a time-frequency volume of order $`1`$. Finally, one computes the total power in the various rectangles in this time-frequency plane. We call this procedure the *long FFT method*. We now consider the computational cost of each method in turn and argue that the second method is more computationally efficient.

### A Sufficiency of approximate version of statistic

In Sec. II A above we discussed how to compute an approximate version of the excess power statistic \[cf. Eq. (23)\]. Namely, for a given start time $`t_s`$, perform a Fourier transform of the $`K`$ point time series $`\{h_j\}`$ corresponding to the time window $`\delta t`$, where $`K=\delta t/\mathrm{\Delta }t`$. Denote the discrete Fourier Transform (DFT) by $`\stackrel{~}{h}_k`$ where $`0kK/2`$. Identify the frequency components $`k_1k<k_2`$ of $`\stackrel{~}{h}_k`$ which belong to the frequency band $`\delta f`$, and construct the statistic

$$=4\underset{k_1k<k_2}{}|\stackrel{~}{h}_k|^2/S_k,$$
(74)
where $`S_k`$ is the noise power spectrum defined in Eq. (22).

The quantity (74) differs from the exact statistic $``$ due to the fact that the expectation value $`\stackrel{~}{h}_k\stackrel{~}{h}_k^{}^{}`$ is not diagonal. (It becomes effectively diagonal only in the limit $`\delta t\mathrm{}`$.) Consequently, the expression (74) is not a sum of squares of independent unit-variance Gaussian random variables, and so its distribution could in principle differ from the non-central $`\chi ^2`$ distribution. However, in practice, if the power spectrum of the noise is a slowly varying function of frequency, then the correlations introduced by using the expression (74) are small. To confirm this, we have examined the behavior of the statistic (74) computed from the DFT of colored, Gaussian noise. We generated colored noise according to the correlation generating scheme

$$n_j=\left(m_j0.8m_{j1}+1.2944n_{j1}0.64n_{j2}\right)/1.3145$$
(75)
where $`m_j`$ are uncorrelated Gaussian deviates and $`n_j=m_j=0`$ for $`j<0`$. To determine detection statistics, we used the signal model

$$s_j=S\mathrm{exp}[16(j/2N1)^2]\mathrm{cos}(2\pi jf_0)$$
(76)
with $`N=4096`$ samples in the signal. The central frequency of the signal was $`f_0=600/4096`$ and the constant $`S`$ was chosen to give the required value of signal amplitude $`A`$. We found the operating characteristics of $``$ were not significantly affected by using the approximate formula (74) rather than the exact formula. This is demonstrated in Figs. 5 and 6 where we have overlaid simulated false alarm and true detection probabilities on top of the distributions computed in Sec. II B. We calculated the goodness of fit using a $`\chi ^2`$ test for a few of the curves in these figures. In each case, the reduced $`\chi ^2`$ value was $`1.03`$, indicating that it is unlikely that the simulated data is drawn from distributions other than those presented in Sec II B. We therefore conclude that we can use the approximate formula (74) without significantly modifying the behavior of the statistic $``$.

### B The short FFT method

The algorithm is: (1) pick a start time, (2) pick a time duration $`\delta t`$, (3) FFT the selected data and compute the power in each frequency bin, (4) sum the power in the bands of interest, (5) loop over steps (2)–(4) until all time durations are used, and (6) repeat steps (1)–(5) for all start times.

The computational cost for steps (3)–(4) can be estimated as follows. The number of data points in segment of data of duration $`\delta t`$ is $`N=\delta t/\mathrm{\Delta }t`$, so each FFT requires $`3N\mathrm{log}_2N`$ floating point operations. Since the relevant frequency band has a dimensionless bandwidth\**\**\**By dimensionless bandwidth we mean number of frequency bins, i.e., bandwidth multiplied by $`\mathrm{\Delta }t`$. Note that the dimensionless bandwidth of the entire frequency band up to the Nyquist frequency is $`1/2`$. $`\alpha _{\text{max}}=\delta f_{\text{max}}\mathrm{\Delta }t`$, the total cost to compute the power in each frequency bin in the band is $`3N\alpha _{\text{max}}`$ operations. Now, for the $`k`$th frequency bin, it costs $`(N\alpha _{\text{max}}k)`$ floating point operations to compute the power in all frequency intervals whose have lowest frequency component is in the $`k`$th bin; the number of operations required to do this for all frequencies $`k`$ is

$$\underset{k=1}{\overset{N\alpha _{\text{max}}1}{}}(N\alpha _{\text{max}}k)=\frac{N\alpha _{\text{max}}}{2}\left(N\alpha _{\text{max}}1\right).$$
(77)
These steps must be repeated for each $`N`$ with $`N_{\text{min}}NN_{\text{max}}`$, where $`N_{\text{min}}=\delta t_{\text{min}}/\mathrm{\Delta }t`$ and $`N_{\text{max}}=\delta t_{\text{max}}/\mathrm{\Delta }t`$. Thus, the total computational cost per start time $`C_{\text{short}}`$ is

$$C_{\text{short}}=\underset{N=N_{\text{min}}}{\overset{N_{\text{max}}}{}}N[3\mathrm{log}_2N+3\alpha _{\text{max}}+\frac{1}{2}\alpha _{\text{max}}(N\alpha _{\text{max}}1)].$$
(78)
One will typically have $`N_{\text{min}}1`$ and $`N_{\text{max}}\alpha _{\text{max}}=V1`$, where $`V`$ is the total time-frequency volume to be searched. In this case, a useful approximation to the computational cost per start time is

$$C_{\text{short}}\frac{V^2}{2\alpha _{\text{max}}}\left(\frac{3\mathrm{log}_2V}{\alpha _{\text{max}}}+\frac{V}{3}\right).$$
(79)
The total computational cost in flops (floating-point operations per second) can be obtained by multiplying $`C_{\text{short}}`$ by the sampling rate.

### C The long FFT method

As discussed above, in the long FFT method one constructs a time-frequency plane consisting of $`N_{\text{channels}}=\delta f_{\text{max}}/\delta f_{\text{min}}`$ different channels. The power in any time-frequency window can then be computed by summing the power in that region of the time-frequency plane. The data stream is broken into chunks of length $`M`$ points, each chunk is FFTed, and the requisite number of channels are produced by bandpass filtering, Fourier transforming back into the time domain, and squaring the time samples. For each chunk, the computational cost of this step is

$$C_1=M[3(1+N_{\text{channels}})\mathrm{log}_2M+N_{\text{channels}}].$$
(80)
To search over the time-frequency plane, we first pick the frequency interval $`\delta f`$ and construct the $`\delta f_{\text{max}}/\delta f`$ channels of this bandwidth; this requires $`M(\delta f/f_{\text{min}}1)`$ additions. For each of these new channels, we sum up the power in various time intervals. This step requires $`N_{\text{max}}`$ operations per start time, of which there are $`MN_{\text{min}}`$. Thus the total cost at this stage of the search is

$$C_2=\underset{j=1}{\overset{N_{\text{channels}}}{}}\frac{N_{\text{channels}}}{j}[M(j1)+N_{\text{max}}(MN_{\text{min}})].$$
(81)
Since there are approximately $`M`$ different start times, the cost per start time is given by the approximate formula

$$C_{\text{long}}\alpha _{\text{max}}^1V^2\mathrm{ln}V.$$
(82)
### D Comparison of the two methods

The space of time-frequency windows to search over was delineated in the Introduction for the initial interferometers in LIGO. We adopt the corresponding parameter values $`\delta f_{\text{min}}=2\text{Hz}`$, $`\delta f_{\text{max}}=200\text{Hz}`$, $`\delta t_{\text{min}}=0.005\text{s}`$, and $`\delta t_{\text{max}}=0.5\text{s}`$. The computational power required using the long FFT method is $`0.3\text{GFlops}`$, which saves a factor of $`14`$ over the short FFT method if $`\mathrm{\Delta }t=0.001`$ seconds.

In general, the computational gain afforded by the long FFT method over the short FFT method is given approximately by

$$\frac{C_{\text{short}}}{C_{\text{long}}}\frac{3}{2\mathrm{ln}2}\frac{1}{\alpha _{\text{max}}}+\frac{1}{6}\frac{V}{\mathrm{ln}V}.$$
(83)
The first term shows that there is at least a factor of $`4`$ to be gained by the long FFT method; in addition, the computational gain increases with the total time-frequency volume $`V`$. For $`V=100`$, the second term is also $`4`$.

There is a further benefit to the long FFT technique. It allows finer frequency resolution in the choice of starts $`f_s`$ and ends $`f_s+\delta f_s`$ of the frequency bands to be explored (although the above estimates of computational cost were for a search equivalent to the short FFT search). Moreover, as part of a hierarchical search, the long FFT method has a further advantage in that it allows follow ups to be made without significant further computations. The next stage of a hierarchical search might involve techniques other than the excess-power method, e.g., Hough transforms or other line tracking algorithms.

## V Multiple detectors

The network of gravitational wave detectors under construction around the world brings benefits that a single instrument cannot. This is especially true for “blind” search techniques, such as the power statistic. Since these techniques do not require the signal to have a specific form, random noise glitches are much more likely to meet the detection criteria than is the case for signal-specific searches such as matched filtering. Multiple-detector statistics will be much more efficient at rejecting such false alarms than single-detector statistics . In this section we consider the construction of the optimal detection strategy for a network of detectors. The derivation requires further formal development. For maximum clarity, we introduce most of our notation in Sec. V A. We derive the multi-instrument detection statistic for a network of aligned detectors in Sec. V B. The two LIGO interferometers at the Hanford site form such a network. In addition, if we ignore the slight misalignment that arises from curvature of the earth, we can also include the interferometer in Livingston to form a three interferometer network. The general case when not all instruments are aligned is treated in Sec. V C.

Our analysis is based on the formalism of Ref.  which followed earlier work of Ref. . We assume that the noise of the detector network is Gaussian. Even though we allow correlations between noise in different instruments, the assumption of Gaussian noise is a serious limitation since the main benefit of having several detectors is to combat non-Gaussian noise. It should be possible to adapt the theoretical models of non-Gaussian noise given in Ref.  in order to derive robust multi-detector statistics. However, it is necessary to understand first the Gaussian case.

### A Notation and terminology

Suppose the detector network consists of $`n_d`$ detectors. Denote the output of the entire network by the vector of time series

$`\stackrel{}{𝐡}`$ $`=\{`$ $`h_j^A\}`$ (84)
$`=\{`$ $`h^1(0),h^1(\mathrm{\Delta }t),\mathrm{},h^1[(N1)\mathrm{\Delta }t],`$ (86)
$`h^2(0),\mathrm{},h^2[(N1)\mathrm{\Delta }t],\mathrm{},h^{n_d}[(N1)\mathrm{\Delta }t]\}`$
where $`A\{1,2,\mathrm{},n_d\}`$ and $`j\{0,1,\mathrm{},N1\}`$. We assume that the noise $`\stackrel{}{𝐧}`$ of the network follows a multivariate Gaussian distribution which is determined by the $`(Nn_d)^2`$ correlations

$$n^A(i\mathrm{\Delta }t)n^B(j\mathrm{\Delta }t)=C^{AB}(|ij|\mathrm{\Delta }t)=R_{ij}^{AB}$$
(87)
where $``$ denotes an ensemble average. In general, the elements of the network correlation matrix are

$$(𝐑)_{AN+i,BN+j}=R_{ij}^{AB}.$$
(88)
The convention (88) for combining the capital and lower case indices to form in an $`Nn_d\times Nn_d`$ matrix is used from here on. The probability density function of the noise is given by

$$p(\stackrel{}{𝐧})=[(2\pi )^{Nn_d}det𝐑]^{1/2}\mathrm{exp}[\frac{1}{2}(\stackrel{}{𝐧},\stackrel{}{𝐧})]$$
(89)
where the inner product is given by

$$(\stackrel{}{𝐩},\stackrel{}{𝐪})=\underset{A,B=1}{\overset{n_d}{}}\underset{i,j=0}{\overset{N1}{}}p_i^AQ_{ij}^{AB}q_j^B$$
(90)
and

$$Q_{ij}^{AB}=(𝐑^1)_{AN+i,BN+j}.$$
(91)
For later convenience, we note that the inner product (90) can be written in terms of the discrete Fourier transforms of the detector time series, that is

$$(\stackrel{}{𝐩},\stackrel{}{𝐪})\underset{A,B=1}{\overset{n_d}{}}4\mathrm{R}\mathrm{e}\underset{k=0}{\overset{N/2}{}}\stackrel{~}{p}_k^A(𝐒_k^1)^{AB}\stackrel{~}{q}_k^B$$
(92)
where

$$𝐒_k^{AB}\delta _{kk^{}}2\underset{j,j^{}=0}{\overset{N1}{}}e^{2\pi ijk/N}R_{jj^{}}^{AB}e^{2\pi ij^{}k^{}/N}.$$
(93)
The notation $`N/2`$ denotes the greatest integer less than or equal to $`N/2`$. These relations are strictly speaking valid only in the continuum and infinite time limits $`\mathrm{\Delta }t0`$ and $`N\mathrm{\Delta }t\mathrm{}`$. Nevertheless, they are sufficiently accurate for most practical applications. Finally, we note that the likelihood ratio $`\mathrm{\Lambda }(\stackrel{}{𝐡};\stackrel{}{𝐬})`$ is given by

$$\mathrm{\Lambda }(\stackrel{}{𝐡};\stackrel{}{𝐬})=\frac{p(\stackrel{}{𝐡}|\stackrel{}{𝐬})}{p(\stackrel{}{𝐡}|\stackrel{}{\mathrm{𝟎}})}=\mathrm{exp}[(\stackrel{}{𝐡},\stackrel{}{𝐬})\frac{1}{2}(\stackrel{}{𝐬},\stackrel{}{𝐬})].$$
(94)
### B Aligned detectors

The simplest type of multi-instrument network to analyze is a network consisting of instruments which all respond to the same polarization component of the gravitational wave field. The two LIGO interferometers in Hanford form such a detector, and if we ignore the slight misalignment arising from the curvature of the earth (a $`10\%`$ correction effect; see Table I) the third LIGO interferometer in Livingston can also be included.

The signal at any detector is simply a time-delayed version of the signal that would be detected at the coordinate origin, which for simplicity we take to be at the center of the earth. Thus, the signal at detector $`A`$ is

$$s^A(t)=s(t+\tau _A),$$
(95)
where $`\tau _A`$ is the time of flight for a gravitational wave between detector $`A`$ and the coordinate origin, and $`s(t)`$ is the signal at the coordinate origin. The time delays $`\tau _A`$ depend on the direction to the source: if $`𝐦`$ is a unit three vector in the direction of propagation of the gravitational wave (i.e., opposite to the direction to the source) and $`𝐱^A`$ is the location of detector $`A`$, then

$$\tau _A=𝐦𝐱^A/c$$
(96)
where $`c`$ is the speed of light. Finally, the DFT of the signal at detector $`A`$ can be written as

$$\stackrel{~}{s}_k^A=e^{2\pi ik\mathrm{\Delta }_A(𝐦)/N}\stackrel{~}{s}_k$$
(97)
where $`\mathrm{\Delta }_A(𝐦)=\tau _A/\mathrm{\Delta }t`$ and $`\stackrel{~}{s}_k`$ is the DFT of the signal at the origin of coordinates.

A convenient description of the multi-instrument network response is the effective strain $`\stackrel{~}{h}_k^{\text{(eff)}}`$ defined by

$$\stackrel{~}{h}_k^{\text{(eff)}}=S_k^{\text{(eff)}}\underset{A,B=1}{\overset{n_d}{}}e^{2\pi ik\mathrm{\Delta }_A(𝐦)/N}(𝐒_k^1)^{AB}\stackrel{~}{h}_k^B$$
(98)
where

$$1/S_k^{\text{(eff)}}=\underset{A,B=1}{\overset{n_d}{}}e^{2\pi ik\mathrm{\Delta }_A(𝐦)/N}(𝐒_k^1)^{AB}e^{2\pi ik\mathrm{\Delta }_B(𝐦)/N}.$$
(99)
We note that the effective strain depends on the direction $`𝐦`$ to the putative source through the time delays $`\mathrm{\Delta }_A(𝐦)`$. Generally the same is true for $`S_k^{\text{(eff)}}`$, although not if there are no correlations between the instrumental noise at separated sites .\*†\*†\*†It might be reasonable to make this assumption for the three detectors at LIGO for example. We can now write the likelihood ratio (94) as

$$\mathrm{\Lambda }(\stackrel{}{𝐡};𝐬)=\mathrm{exp}[((\stackrel{~}{h}^{\text{(eff)}},\stackrel{~}{s}))\frac{1}{2}((\stackrel{~}{s},\stackrel{~}{s}))].$$
(100)
where

$$((\stackrel{~}{p},\stackrel{~}{q}))=4\mathrm{R}\mathrm{e}\underset{k=0}{\overset{N/2}{}}\stackrel{~}{p}_k\stackrel{~}{q}_k^{}/S_k^{\text{(eff)}}.$$
(101)
The posterior probability $`p_s(\stackrel{}{𝐡})`$ that a signal is present given the data $`\stackrel{}{𝐡}`$ is determined by integrating the likelihood against prior probabilities densities for the signal and for the source direction $`𝐦`$. Thus

$$\frac{p_s(\stackrel{}{𝐡})}{1p_s(\stackrel{}{𝐡})}=\mathrm{\Lambda }(\stackrel{}{𝐡})\frac{p_{s0}}{1p_{s0}}$$
(102)
where

$$\mathrm{\Lambda }(\stackrel{}{𝐡})=p(\theta ,\varphi )d^2\mathrm{\Omega }p(𝒯)𝑑𝒯p(𝐬|𝒯)d^N𝐬\mathrm{\Lambda }(\stackrel{}{𝐡};\stackrel{}{𝐬})$$
(103)
and $`p_{s0}`$ is the a priori probability that any gravitational wave sources exist. The mechanism of how information about the signals is encoded in the prior probabilities $`p(𝐬|𝒯)`$ and $`p(𝒯)`$ is treated in Sec. III A, and applies directly to the current context with only one modification: the inner product $`(,)`$ should be replaced by $`((,))`$. In particular, the probability distribution $`p(𝐬|T)`$ is given by Eq. (50). The function $`p(\theta ,\varphi )`$ is the expected distribution of source directions. For sources that are mostly further than $`30\text{Mpc}`$, the distribution should be uniform on the sphere.

The integral in the expression (103) for the likelihood function includes a sum over all time-frequency windows $`𝒯`$ and all source directions $`𝐦`$. However, it is nearly equivalent, and much easier, to adopt the *maximum* term in the sum as an approximation to the likelihood function, since the largest term will dominate the sum when a signal is present. It is therefore sufficient to consider only a single time-frequency region $`𝒯`$ and fixed direction $`𝐦`$ in the remainder of this section, with the understanding that the detection statistic will include a maximization over these variables .

Using arguments similar to those in Sec. III B, we can perform the integral over signals $`𝐬`$. In particular, the prior probability in Eq. (50) restricts $`𝐬`$ to lie in a vector space $`𝒱_𝒯`$ which contains only signals in the time-frequency window $`𝒯`$. This has the effect of replacing $`𝐬`$ by $`𝐬_{}`$ in the inner products. Moreover, the inner product $`((,))`$ induces a natural inner product on the subspace $`𝒱_𝒯`$. The integral over angles can be performed as in Eq. (62) to show that

$$\mathrm{\Lambda }(\stackrel{}{𝐡})=\mathrm{\Gamma }(V)e^{A^2/2}(A^{1/2}/2)^{1V}I_{V1}(A^{1/2})p(A)𝑑A$$
(104)
where

$$=((\stackrel{~}{h}_{}^{\text{(eff)}},\stackrel{~}{h}_{}^{\text{(eff)}}))=4\mathrm{R}\mathrm{e}\underset{k_1k<k_2}{}\stackrel{~}{h}_k^{\text{(eff)}}\stackrel{~}{h}_k^{\text{(eff)}}/S_k^{\text{(eff)}}.$$
(105)
Here $`\stackrel{~}{h}_{}^{\text{(eff)}}`$ is the projection of the effective strain into the time-frequency subspace $`𝒱_𝒯`$; the second equality holds for a particular time-frequency window in which the signal has duration $`N\mathrm{\Delta }t`$ and is localized to the frequency band $`k_1k<k_2`$. The amplitude $`A`$ is defined by the $`A^2=((\stackrel{~}{s}_{},\stackrel{~}{s}_{}))`$. Since the right hand side of Eq. (104) is a monotonically increasing function of $``$, the Neyman-Pearson theorem tells us that $``$ provides the optimal multi-instrument detection statistic. Note that, as mentioned above, this detection statistic includes an implicit maximization over all source directions $`𝐦`$ and time-frequency windows $`𝒯`$.

### C General networks of detectors

When the network contains at least one instrument with a different orientation to the others, it is necessary to discuss the two degrees of freedom, or polarizations, of the gravitational wave signal. We denote these two independent signals as $`s^+(t)`$ and $`s^\times (t)`$, where the definition is with respect to a radiation coordinate system associated with the gravitational waves. (See Appendix B for a detailed discussion of these and the other coordinate systems relevant to this section.) For the $`A`$th detector in the network, the gravitational wave strain is

$$s^A(t)=F_+^As^+(t+\tau _A)+F_\times ^As^\times (t+\tau _A)$$
(106)
where $`F_+^A`$, $`F_\times ^A`$ are the detector beam pattern functions and $`\tau _A=𝐦𝐱^A/c`$ is the time delay between the origin of earth-fixed coordinates and the detector $`A`$ (located at $`𝐱^A`$) for a wave propagating in the direction $`𝐦`$.

The concept of effective strain is particularly useful in the formal development of the multi-instrument detection statistic for gravitational waves. To introduce this concept, consider the inner product

$$(\stackrel{}{𝐬},\stackrel{}{𝐬})=\underset{A,B=1}{\overset{n_d}{}}4\mathrm{R}\mathrm{e}\underset{k=0}{\overset{N/2}{}}\stackrel{~}{s}_k^A(𝐒_k^1)^{AB}\stackrel{~}{s}_k^B.$$
(107)
The DFT of the signal at the $`A`$th detector is

$$\stackrel{~}{s}_k^A=e^{2\pi i\mathrm{\Delta }_A(𝐦)k/N}(F_+^A\stackrel{~}{s}_k^++F_\times ^A\stackrel{~}{s}_k^\times )$$
(108)
where $`\mathrm{\Delta }_A(𝐦)=\tau _A/\mathrm{\Delta }t`$ is the discrete time delay, $`\mathrm{\Delta }t`$ is the sampling rate, and $`\stackrel{~}{s}^{+,\times }`$ are the DFTs of the plus and cross polarizations of the signal defined in Eq. (106). The inner product can be rewritten as

$$(\stackrel{}{𝐬},\stackrel{}{𝐬})=\underset{\alpha ,\beta =+,\times }{}4\mathrm{R}\mathrm{e}\underset{k=0}{\overset{N/2}{}}\stackrel{~}{s}_k^\alpha \mathrm{\Theta }_{\alpha \beta }^k\stackrel{~}{s}_k^\beta $$
(109)
where

$$\mathrm{\Theta }_{\alpha \beta }^k=\underset{A,B=1}{\overset{n_d}{}}e^{2\pi i(\mathrm{\Delta }_A\mathrm{\Delta }_B)k/N}F_\alpha ^A(𝐒_k^1)^{AB}F_\beta ^B.$$
(110)
We can now introduce a pair of effective strains, corresponding to the $`+`$ and $`\times `$ gravitational wave signals for the network, by

$$\stackrel{~}{h}_k^\alpha =\underset{\beta =+,\times }{}\mathrm{\Theta }_k^{\alpha \beta }\underset{A,B=1}{\overset{n_d}{}}F_\alpha ^Ae^{2\pi i\mathrm{\Delta }_Ak/N}(𝐒_k^1)^{AB}\stackrel{~}{h}_k^B$$
(111)
where $`_{\beta =+,\times }\mathrm{\Theta }_k^{\alpha \beta }\mathrm{\Theta }_{\beta \gamma }^k=\delta _\gamma ^\alpha `$. In terms of the effective strains $`\stackrel{~}{h}_k^\alpha `$ and the signals $`s_{+,\times }`$, the likelihood ratio is

$$\mathrm{\Lambda }(\stackrel{}{𝐡};\stackrel{}{𝐬})=\mathrm{exp}[((\stackrel{~}{h},\stackrel{~}{s}))\frac{1}{2}((\stackrel{~}{s},\stackrel{~}{s}))]$$
(112)
where the inner product is now defined as

$$((\stackrel{~}{p},\stackrel{~}{q}))=\underset{\alpha ,\beta =+,\times }{}4\mathrm{R}\mathrm{e}\underset{k=0}{\overset{N/2}{}}\stackrel{~}{p}_k^\alpha \mathrm{\Theta }_{\alpha \beta }^k\stackrel{~}{q}_k^\beta .$$
(113)
The effective strains and the inner product (113) depend on the direction $`𝐦`$ to the putative source. Consequently the probability that a signal with plus polarization $`s^+`$ and cross polarization $`s^\times `$ is present in the data stream is given by Eqs. (102) and (103) with $`\mathrm{\Lambda }(\stackrel{}{𝐡};\stackrel{}{𝐬})`$ defined by Eq. (112) where the measure $`p(𝐬|𝒯)`$ on the space of signals is defined as follows.

The signal $`\{s^+,s^\times \}`$ now belongs to a $`4V`$-dimensional vector space which is the tensor product of two copies of $`𝒱_𝒯`$. Within this vector space, all directions can also be considered to be equally likely. These assumptions about the signal $`\{s_+,s_\times \}`$ reduce to the assumption in Sec. V B when the detectors are aligned. Moreover, the reasoning from that section can be readily applied here provided one understands that the vector space of signals is now $`4V`$-dimensional and that all angles and lengths are measured using the inner product (109). Thus, the integrals can be carried out in much the same way to arrive at the excess power statistic for a multiple-instrument network:

$$=((\stackrel{~}{h}_{},\stackrel{~}{h}_{}))=\underset{\alpha ,\beta =+,\times }{}4\mathrm{R}\mathrm{e}\underset{k_1k<k_2}{}\stackrel{~}{h}_k^\alpha \mathrm{\Theta }_{\alpha \beta }^k\stackrel{~}{h}_k^\beta .$$
(114)
As before there is an implicit maximization over time-frequency windows and source directions.

Since the effective strains are linear combinations of the outputs of each of the detectors in the network, the statistic (114) is a bilinear function of the outputs of all the detectors, containing both auto-correlation terms from each detector individually and cross-correlation terms between each pair of detectors. It is the optimal statistic in Gaussian noise. When the noise in the instruments is non-Gaussian, it remains to be seen what is the best strategy. One obvious strategy is to simply omit the auto-correlation terms in Eq. (114) and retain only the cross-correlation terms; the resulting statistic will share many of the nice features of $``$ and be more robust against non-Gaussian noise bursts. The real challenge is to derive the optimal statistic in the presence of uncorrelated noise bursts which are Poisson distributed in time. It is likely that the model introduced in Ref.  can be used to address this issue.

###### Acknowledgements.

We are grateful to Bruce Allen, Sam Finn, Soumya Mohanty, Julien Sylvestre and Kip Thorne for helpful discussions. This work was supported by National Science Foundation awards PHY 9970821, PHY 9722189, PHY 9728704, PHY-9407194, and Phy-9900776. EF acknowledges the support of the Alfred P. Sloan foundation.

## A Related detection statistics

In this appendix we discuss some detection statistics that can be obtained from the Bayesian formalism discussed in Sec. III starting from different prior PDFs for signals $`𝐬`$.

### 1 Known signal spectrum

Suppose that one knows, in addition to the duration and frequency band, the spectrum of the expected signal, but that one does not know the phase evolution. Let us adopt the Fourier basis (assuming that the autocorrelation matrix is reasonably close to diagonal in this basis) and assume that the noise is stationary and Gaussian. Then the likelihood ratio is

$$\mathrm{\Lambda }(𝐡;𝐬)=\mathrm{exp}\left[4\underset{k=0}{\overset{N/2}{}}(|\stackrel{~}{h}_k||\stackrel{~}{s}_k|\mathrm{cos}\varphi _k\frac{1}{2}|\stackrel{~}{s}_k|^2)/S_k\right]$$
(A1)
where $`\varphi _k`$ represents the relative phase difference between the data and the expected signal in the $`k`$th frequency bin. Since the signal phases are considered unknown, we should integrate out these angles to obtain the integrated likelihood ratio

$$\mathrm{\Lambda }(𝐡;\{P_k\},A)=\underset{k=0}{\overset{N/2}{}}2\pi e^{A^2P_k/S_k}I_0(2^{3/2}AP_k^{1/2}|\stackrel{~}{h}_k|/S_k)$$
(A2)
where the one-sided signal spectrum is given by $`\frac{1}{2}A^2P_k=|\stackrel{~}{s}_k|^2`$ with $`2_{k=0}^{N/2}P_k/S_k=1`$.

In the limit of weak signal amplitude $`A`$, we can approximate the likelihood ratio of Eq. (A2) by its expansion in powers of $`A`$. The first non-trivial term is the *locally optimal*  detection statistic

$$\frac{d^2\mathrm{ln}\mathrm{\Lambda }(𝐡;\{P_k\},A)}{dA^2}|_{A=0}=\text{(const)}+4\underset{k=0}{\overset{N/2}{}}P_k|\stackrel{~}{h}_k|^2/S_k^2.$$
(A3)
This statistic is the weighted average of the detector output power in each frequency bin. Unfortunately, it is not possible to get simple expressions for the false alarm and false dismissal probabilities for this statistic; one needs to use numerical methods to obtain these given a known signal power spectrum $`\{P_k\}`$.

### 2 Non-Gaussian noise

It is unlikely that a gravitational wave detector will produce purely stationary and Gaussian noise. In the case that the detector noise distribution is known, we can obtain a detection statistic for unknown signals using the Bayesian methodology. Unfortunately the most general noise distribution contains many free functions and will not be known in practice. However, constructing simple analytic non-Gaussian noise models and the associated detection statistics can give us insight into what kind of statistics to try out with real detectors.

One such simple model is as follows. We assume that the detector noise is stationary, and, as before, that each frequency bin in the Fourier basis is uncorrelated. Let make the additional assumption that the power in each frequency bin is independentally distributed, while the phases of each frequency bin are uniform and independent. Then

$$p(𝐧)d^N𝐧=\underset{k=1}{\overset{(N1)/2}{}}f_k(|\stackrel{~}{n}_k|^2/S_k)\frac{d|\stackrel{~}{n}_k|^2}{S_k}\frac{d\mathrm{arg}\stackrel{~}{n}_k}{2\pi }$$
(A4)
where $`f_k(x)`$ are known non-exponential probability distribution functions (exponential functions would correspond to Gaussian noise). Here we have omitted the DC (and a possible Nyquist) component. The likelihood ratio is

$`\mathrm{\Lambda }(𝐡;𝐬)`$ $`=`$ $`{\displaystyle \underset{k=1}{\overset{(N1)/2}{}}}{\displaystyle \frac{f_k(|\stackrel{~}{h}_k\stackrel{~}{s}_k|^2/S_k)}{f_k(|\stackrel{~}{h}_k|^2/S_k)}}`$ (A5)
$`=`$ $`{\displaystyle \underset{k=1}{\overset{(N1)/2}{}}}\{12{\displaystyle \frac{\mathrm{Re}(\stackrel{~}{h}_k\stackrel{~}{s}_k^{})}{S_k}}{\displaystyle \frac{f_k^{}(|\stackrel{~}{h}_k|^2/S_k)}{f_k(|\stackrel{~}{h}_k|^2/S_k)}}`$ (A8)
$`+[2\left({\displaystyle \frac{\mathrm{Re}(\stackrel{~}{h}_k\stackrel{~}{s}_k^{})}{S_k}}\right)^2{\displaystyle \frac{f_k^{\prime \prime }(|\stackrel{~}{h}_k^2/S_k)}{f_k(|\stackrel{~}{h}_k^2/S_k)}}`$
$`+{\displaystyle \frac{|\stackrel{~}{s}_k|^2}{S_k}}{\displaystyle \frac{f_k^{}(|\stackrel{~}{h}_k|^2/S_k)}{f_k(|\stackrel{~}{h}_k^2/S_k)}}]+O(|\stackrel{~}{s}_k|^3)\}.`$
We have expanded the likelihood ratio in powers of the (presumed small) signal in order to construct the locally optimal detection statistic .

To compute the integrated likelihood function we need to integrate over our prior knowledge of signals. Let us suppose that we do not know the signal phase evolution; then we can integrate over the unknown phases $`\mathrm{arg}\stackrel{~}{s}_k`$ in each frequency bin. We find

$$\mathrm{\Lambda }(𝐡;\{P_k\},A)=1+\frac{A^2}{2}\underset{k=1}{\overset{(N1)/2}{}}\frac{|\stackrel{~}{s}_k|^2}{S_k}g_k(|\stackrel{~}{h}_k|^2/S_k)+O(A^4)$$
(A9)
where $`\frac{1}{2}A^2P_k=|\stackrel{~}{s}_k|^2`$ and

$$g_k(x)=[xf_k^{\prime \prime }(x)+f_k^{}(x)]/f_k(x).$$
(A10)
For Gaussian noise, $`f_k(x)=e^x`$ and $`g_k(x)=x1`$ for all $`k`$, which gives essentially the same detection statistic as in Eq. (A3). For a probability distribution with tails that decrease more slowly in the $`k`$th bin, e.g., $`f_k(x)(1x/2)^2`$, then we have $`g_k(x)=(x1)/(1x/2)^2`$, which increases with $`x`$ up to $`x=2`$, and then decreases for larger values of $`x`$. Thus, large amounts of excess power in the $`k`$th bin are *suppressed*.

When the signal is known to be band-limited to frequency bins $`k_1k<k_2`$, but we have no reason to believe that any particular bin in the band will contribute more to the overall signal-to-noise ratio than any other bin, then we obtain the locally optimal statistic by assuming a uniform weighting of the terms in Eq. (A9). Thus the locally optimal statistic is

$$\underset{k_1k<k_2}{}g_k(|\stackrel{~}{h}_k|^2/S_k).$$
(A11)
In the case of Gaussian noise this is the excess power statistic; for noise models with larger tails, the components of the sum are attenuated if they have large power.

## B Multidetector amplitudes

In Sec. V C, we discussed the detection of burst signals using multiple detectors. When the detectors are not aligned, one needs the response functions of each detector in the network to a gravitational wave signal from a given sky position. Here we define reference coordinates to which we refer each detectors response. Consider a coordinate system fixed at the center of the earth. In terms of latitude and longitude $`\{\phi ,\lambda \}`$, the coordinate axes are oriented so that the $`x`$-axis pierces the earth at $`\{000,000\}`$, the $`y`$-axis pierces the earth at $`\{000,090^{}\text{E}\}`$, and the $`z`$-axis pierces the earth at $`\{090^{}\text{N},000\}`$. We denote the location of a source on the celestial sphere by standard spherical polar coordinates $`\{\theta ,\varphi \}`$ measured with respect to this earth fixed frame. A fiducial signal $`𝐬`$ comes from a sky position with right ascension $`\alpha =\varphi +\text{GMST}`$ (GMST is the Greenwich mean sidereal time of arrival of the signal), declination $`\delta =\pi /2\theta `$, and has polarization angle $`\psi `$—the angle (counter-clockwise about the direction of propagation) from the line of nodes to the $`𝐗`$-axis of the signal coordinates. In particular, this gravitational wave signal can be represented by a tensor

$$s_{ij}=s^+(𝐞_+)_{ij}+s^\times (𝐞_\times )_{ij}$$
(B1)
where the polarization tensors are given by

$`(𝐞_+)_{ij}`$ $`=`$ $`(𝐗𝐗𝐘𝐘)_{ij}`$ (B2)
$`(𝐞_\times )_{ij}`$ $`=`$ $`(𝐗𝐘+𝐘𝐗)_{ij}.`$ (B3)
The vectors $`𝐗`$ and $`𝐘`$ are the axes of the wave-frame, given explicitly by

$`𝐗`$ $`=`$ $`(\mathrm{sin}\varphi \mathrm{cos}\psi \mathrm{sin}\psi \mathrm{cos}\varphi \mathrm{cos}\theta )𝐢`$ (B5)
$`(\mathrm{cos}\varphi \mathrm{cos}\psi +\mathrm{sin}\psi \mathrm{sin}\varphi \mathrm{cos}\theta )𝐣+\mathrm{sin}\psi \mathrm{sin}\theta 𝐤`$
$`𝐘`$ $`=`$ $`(\mathrm{sin}\varphi \mathrm{sin}\psi \mathrm{cos}\psi \mathrm{cos}\varphi \mathrm{cos}\theta )𝐢`$ (B7)
$`+(\mathrm{cos}\varphi \mathrm{sin}\psi \mathrm{cos}\psi \mathrm{sin}\varphi \mathrm{cos}\theta )𝐣+\mathrm{sin}\theta \mathrm{cos}\psi 𝐤`$
where the polarization angle $`\psi `$ is defined above, and $`𝐢`$, $`𝐣`$ and $`𝐤`$ are unit vectors along the $`x`$, $`y`$ and $`z`$-axes respectively. Note, we use a right handed coordinate system in which the vector $`𝐙=𝐗𝐘`$ points in the direction from the source towards the detector. The waveforms in Refs.  are referred to these coordinates; Thorne uses a different definition in Ref. .

One can characterize the response of an interferometer on the surface of the earth to the impinging gravitational wave using another tensor $`𝐃`$ given by

$$𝐃_{ij}=\frac{1}{2}(𝐧^x𝐧^x𝐧^y𝐧^y)_{ij}$$
(B8)
where $`𝐧^x`$ and $`𝐧^y`$ are unit vectors along the $`x`$ and $`y`$ arms of the interferometer respectively. For a given interferometer $`A`$, it is now straightforward to compute the response

$$s^A=\underset{i,j=1}{\overset{3}{}}D_{ij}^As_{ij}$$
(B9)
and to extract the response functions $`F_{+,\times }^A`$ by comparing the result with the formula:

$$s^A=F_+^As^++F_\times ^As^\times .$$
(B10)
For a detector having its arms aligned with the coordinate axes at the center of the earth, we find

$`F_+`$ $`=`$ $`{\displaystyle \frac{1}{2}}(1+\mathrm{cos}^2\theta )\mathrm{cos}2\varphi \mathrm{cos}2\psi \mathrm{cos}\theta \mathrm{sin}2\varphi \mathrm{sin}2\psi `$ (B11)
$`F_\times `$ $`=`$ $`+{\displaystyle \frac{1}{2}}(1+\mathrm{cos}^2\theta )\mathrm{cos}2\varphi \mathrm{sin}2\psi \mathrm{cos}\theta \mathrm{sin}2\varphi \mathrm{cos}2\psi `$ (B12)
Finally, the response $`𝐃`$ of the various detectors around the world can be determined using the latitude North $`\phi `$, longitude East $`\lambda `$, and arm orientations $`(\psi _{x,y},\omega _{x,y})`$ where $`\psi _{x,y}`$ are the azimuths (North of East) of the $`x`$ and $`y`$ arms and $`\omega _{x,y}`$ are the tilts of the $`x`$ and $`y`$ arms above the horizontal defined by the WGS-84 earth model . This model is an oblate ellipsoid with semi-major axis $`a=\mathrm{6\hspace{0.17em}378\hspace{0.17em}137}\text{m}`$ and semi-minor axis $`b=\mathrm{6\hspace{0.17em}356\hspace{0.17em}752.314}\text{m}`$. The position $`𝐱=x𝐢+y𝐣+z𝐤`$ of a detector at a given latitude $`\phi `$, longitude $`\lambda `$, and elevation $`h`$ above (normal to) the surface is given by

$`x`$ $`=`$ $`[R(\phi )+h]\mathrm{cos}\phi \mathrm{cos}\lambda `$ (B13)
$`y`$ $`=`$ $`[R(\phi )+h]\mathrm{cos}\phi \mathrm{sin}\lambda `$ (B14)
$`z`$ $`=`$ $`[(b^2/a^2)R(\phi )+h]\mathrm{sin}\phi `$ (B15)
where $`R(\phi )=a^2(a^2\mathrm{cos}\phi +b^2\mathrm{sin}\phi )^{1/2}`$ is the local radius of the earth. At this position, the unit vectors pointing East, North, and Up are

$`𝐞_\lambda `$ $`=`$ $`\mathrm{sin}\lambda 𝐢+\mathrm{cos}\lambda 𝐣`$ (B16)
$`𝐞_\phi `$ $`=`$ $`\mathrm{sin}\phi \mathrm{cos}\lambda 𝐢\mathrm{sin}\phi \mathrm{sin}\lambda 𝐣+\mathrm{cos}\phi 𝐤`$ (B17)
$`𝐞_h`$ $`=`$ $`\mathrm{cos}\phi \mathrm{cos}\lambda 𝐢+\mathrm{cos}\phi \mathrm{sin}\lambda 𝐣+\mathrm{sin}\phi 𝐤`$ (B18)
respectively. The unit vector along the $`x`$ arm is then given by

$$𝐧^x=\mathrm{cos}\omega _x\mathrm{cos}\psi _x𝐞_\lambda +\mathrm{cos}\omega _x\mathrm{sin}\psi _x𝐞_\phi +\mathrm{sin}\omega _x𝐞_h$$
(B19)
and similarly for the $`y`$ arm. For completeness, we list these vectors $`𝐧^x`$ and $`𝐧^y`$ for each of the interferometers in Table I. For the two LIGO interferometers, these vectors are provided in . For the other interferometers we used the values in note  (with tilt angles $`\omega =0`$), or the values given in Ref.  (with elevations $`h=0`$ and tilt angles $`\omega =0`$).


format: text

text:
# Twin Peaksaafootnote aPresented by J. Smit.

The on-shell imaginary part of the retarded selfenergy of massive $`\phi ^4`$ theory in 1+1 dimensions is logarithmically infrared divergent. This leads to a zero in the spectral function, separating its usual bump into two. The twin peaks interfere in time-dependent correlation functions, which causes oscillating modulations on top of exponential-like decay, while the usual formulas for the decay rate fail. We see similar modulations in our numerical results for a mean field correlator, using a Hartree ensemble approximation.

In our numerical simulations of 1+1 dimensional $`\phi ^4`$ theory using the Hartree ensemble approximation $`^\mathrm{?}`$ we found funny modulations in a time-dependent correlation function. Fig. 1 shows such modulations on top of a roughly exponential decay. The correlation function is the time average of the zero momentum mode of the mean field, $`F_{\mathrm{mf}}(t)=\overline{\phi (t)\phi (0)}\overline{\phi (t)}\overline{\phi (0)}`$, where the over-bar denotes a time average, $`\overline{X(t)}=_{t_1}^{t_2}𝑑t^{}X(t+t^{})/(t_2t_1)`$, taken after waiting a long time $`t_1`$ for the system to be in approximate equilibrium. This equilibrium is approximately thermal and $`F_{\mathrm{mf}}(t)`$ is analogous to the symmetric correlation function of the quantum field theory at finite temperature, $`F(t)=\frac{1}{2}\{\widehat{\phi }(t),\widehat{\phi }(0)\}_{\mathrm{conn}}.`$ A natural question is now, does $`F(t)`$ also have such modulations?

The function $`F(t)`$ can be expressed in terms of the zero momentum spectral function $`\rho (p^0)`$,

$$F(t)=_{\mathrm{}}^{\mathrm{}}\frac{dp^0}{2\pi }e^{ip^0t}\left(\frac{1}{e^{p^0/T}1}+\frac{1}{2}\right)\rho (p^0),$$
(1)
and the latter in turn in terms of the retarded selfenergy $`\mathrm{\Sigma }(p^0)`$,

$$\rho (p^0)=\frac{2\mathrm{I}\mathrm{m}\mathrm{\Sigma }(p^0)}{[m^2(p^0+iϵ)^2+\mathrm{Re}\mathrm{\Sigma }(p^0)]^2+[\mathrm{Im}\mathrm{\Sigma }(p^0)]^2}.$$
(2)
The selfenergy can be calculated in perturbation theory. The one and two loop diagrams in the imaginary time formalism which have nontrivial energy-momentum dependence are shown in Fig. 2. Diagrams not shown give only rise to an effective temperature dependent mass, which we assume to be the mass in the propagators of the diagrams in Fig. 2, after adding a counterterm that sets the real part of $`\mathrm{\Sigma }`$ to zero at $`p^0=m`$.

The one loop diagram is present only in the ‘broken phase’ (for which $`\widehat{\phi }0`$; there is really only a symmetric phase in 1+1 dimensions, but this is due to symmetry restoration by nonperturbative effects which will not obliterate the one-loop damping.) The corresponding selfenergy has been calculated in $`^\mathrm{?}`$, for example. It only leads to damping for frequencies $`p_0^2>4m^2`$, which are irrelevant for the quasiparticle damping at $`p_0^2=m^2`$. So from now on we concentrate on the two-loop diagram. After analytic continuation to real time one finds that it is given by the sum of two terms, $`\mathrm{\Sigma }_1+\mathrm{\Sigma }_2`$ (see e.g. $`^\mathrm{?}`$). The first has an imaginary part corresponding to $`13`$ processes requiring $`p_0^2>9m^2`$, so it does not contribute to plasmon damping. The second is given by

$`\mathrm{\Sigma }_2`$ $`=`$ $`{\displaystyle \frac{9\lambda ^2}{16\pi ^2}}{\displaystyle \frac{dp_2dp_3}{E_1E_2E_3}\frac{(1+n_1)n_2n_3n_1(1+n_2)(1+n_3)}{p^0+iϵ+E_1E_2E_3}}`$ (3)
$`+\left[(p^0+iϵ)(p^0+iϵ)\right],`$
where $`\lambda `$ is the coupling constant (introduced as $`_1=\lambda \phi ^4/4`$), and $`E_1=\sqrt{m^2+(p_2+p_3)^2}`$, $`E_i=\sqrt{m^2+p_i^2}`$, $`i=2,3`$, $`n_i=[\mathrm{exp}(E_i/T)1]^1`$, $`i=1,2,3`$. Its imaginary part corresponds to $`22`$ processes, which contribute in the regions near $`p_0=\pm m`$.

Now the usual definition of the thermal plasmon damping rate (at zero momentum) in terms of the retarded selfenergy,

$$\gamma =\mathrm{Im}\mathrm{\Sigma }(m)/2m,$$
(4)
leads to a divergent answer (a collinear divergence). A natural way out of this difficulty may be to continue the selfenergy analytically into the lower half of its second Riemann sheet, $`p^0mi\gamma `$, and replace (4) by the improved definition

$$m^2(mi\gamma )^2+\mathrm{\Sigma }(mi\gamma )=0.$$
(5)
The analytic continuation of the selfenergy into the region Im $`p^0<0`$ poses the puzzle how to deal with the logarithmic branch point coming from the collinear singularity at $`p^0=m`$. However, the ambiguity is present only in the real part of $`\mathrm{\Sigma }`$. For weak coupling $`\lambda /m^21`$ we get from (5) the equation

$$\frac{\gamma }{m}=\frac{9\lambda ^2}{16\pi m^4}\frac{e^{m/T}}{\left(e^{m/T}1\right)^2}\left[\mathrm{ln}\frac{m}{\gamma }+c(T)\right].$$
(6)
The constant $`c`$ has to be determined by matching a numerical evaluation of $`\mathrm{\Sigma }`$ to the logarithmic singularity at $`p^0=m`$.

We evaluated $`\mathrm{\Sigma }_2`$ in (3) for $`T=m`$ by numerical integration with $`ϵ/m=0.02`$, $`0.01`$ and linear extrapolation $`ϵ0`$, giving $`c0.51`$. For example, Eq. (6) now gives $`\gamma /m=0.061`$, for $`\lambda /m^2=0.4`$.

To see how well this $`\gamma `$ describes the decay of the correlator $`F(t)`$ we evaluated this function directly from (1) and (2). The divergence in $`\mathrm{Im}\mathrm{\Sigma }(p^0)`$ at $`p^0=m`$ leads to a zero in the spectral function $`\rho (p^0)`$. So is there a peak at all in $`\rho (p^0)`$? Fig. 3 shows what happens: the ‘usual’ peak has separated into two twins! Fig. 4 shows the resulting $`F(t)`$. The effect of the double peak is indeed an oscillating modulation on top of the roughly exponential decay. The decay corresponding to $`\mathrm{exp}(\gamma t)`$, with $`\gamma `$ given by (6), is also indicated in the plot: it does not do a good job in describing the average decay beyond the first interference minimum. The ‘Twin Peaks’ phenomenon implies that the usual definition of damping rate (5) is unreliable in 1+1 dimensions.

Acknowledgements. We thank Gert Aarts for useful conversations. This work is supported by FOM/NWO.