text:
# Time-frequency detection algorithm for gravitational wave bursts

## I Introduction

A number of large laser interferometers GWProjects  are approaching sensitivities to gravitational waves in the $`101000`$ Hz frequency band that could be sufficient for the detection of astrophysical events Thorne . The signals from these events will be buried deep in the instrumental noise, so that unambiguous detections will be possible only with highly efficient data processing algorithms.

The focus of this article will be on transient sources of gravitational radiation, which will be defined as sources that have a relatively short duration (milliseconds to tens of seconds) and a bandwidth which overlaps at least partially with that of the interferometric detectors. A significant number of such transient sources have been theorized, with various levels of sophistication. For instance, the inspiral portion of coalescing compact binaries is well-understood (by post-newtonian expansion techniques inspiral ), as is the ring-down portion if a black hole results from the coalescence Leaver , but the merger portion is understood at best only qualitatively merger ; BH . As the mass of the binary increases, the signal-to-noise ratio of the merger portion dominates that of the inspiral and ringdown portions; the coalescence of $`10M_{}1000M_{}`$ black hole binaries could be visible at large distances, provided that the merger waveform could be detected with sufficient efficiency BH ; BH2 . The collapse of the core of massive stars could also produce detectable signals fryer\_holz ; depending on the type of progenitor, bar modes, r-modes, fragmentation instabilities or black hole ringdowns could be important sources of gravitational radiation. The details of the waveforms for most of these different mechanisms are far from known; even in the best cases, only numerical simulations covering parts of the relevant physics are available in the literature ZM . Hence, as it can be seen from the preceding examples, the amount of information about the gravitational wave signal from various sources varies considerably, and this variability is obviously reflected in the efficiency of the algorithms that can be designed for each class of sources.

Only minimal assumptions about the signal will be made in this paper, and therefore the principal characteristic of the algorithm to be presented will be its robustness against poor modeling of the expected signal waveform. Stated differently, this algorithm will be moderately sensitive to a very large class of signals, by opposition to being very sensitive to only a few specific signals. It will correspondingly be useful to search for transient sources that do not have waveforms that are precisely predicted, and to characterize the non-Gaussian, transient component of the noise in the instruments.

It is important when designing a detection algorithm to compare its efficacy to that of other designs; the development process can stop when an algorithm that outperforms all others is discovered. Various measures of this efficacy can be adopted, and various techniques to obtain the optimal algorithm be used, as it is discussed in section II. It is argued that the algorithm that is presented in this article, the TFCLUSTERS algorithm, has a structure which is close to that of the optimal detector for a large variety of signal classes.

The TFCLUSTERS algorithm is explained in detail in section III. To summarize, it consists of the following four steps:

The data $`𝒚`$ from a gravitational wave detector are transformed into a time-frequency representation with fixed time and frequency resolutions $`T`$ and $`F`$, respectively; the instantaneous power at time $`t=iT`$ and frequency $`f=jF`$, estimated from this representation, is labeled $`P_{ij}(𝒚)`$.

A threshold $`\eta `$ is applied on the power, in order to retain only pixels with $`P_{ij}(𝒚)>\eta `$.

Clusters of pixels with power above threshold are formed by grouping pixels sharing a common side; clusters that do not conform to a fixed set of rules on their size or distance to other clusters are discarded.

A threshold on the sum of $`P_{ij}(𝒚)\eta `$ over the surviving clusters is applied. Data segments containing clusters satisfying this threshold lead to the acceptance of the hypothesis that they contain a signal.

In order to understand the operating characteristics of this algorithm, a simplified version \[without the clustering analysis, i.e. consisting only of steps (i), (ii) and (iv)\] is shown in section IV to maximize for all signals in $`𝖱^N`$ the signal-to-noise ratio among all detectors that are based on the estimation of a lower bound of the signal power. This simplified algorithm is especially efficient for signals with a sparse representation in the time-frequency domain. Since most signals are expected to have pixels that present a high degree of spatial correlation in the time-frequency domain, the clustering analysis of TFCLUSTERS is an interesting way to capitalize on that property to filter out a large portion of the noise, as it is shown in section V.

An analytical method for computing the false alarm rate associated with TFCLUSTERS is developed, with the details for the clustering analysis being presented in section V.1. Using the computer-generated enumeration of all the possible clusters of a certain size that can be formed, it is shown that large clusters are exponentially unlikely to occur when only noise is present. The rate of occurrence of pairs of clusters separated by a certain distance is computed in an analogous manner. One example of a complete analysis of the performances of TFCLUSTERS for a short, narrow-band signal, including the optimization of its efficiency at fixed false alarm rate, is presented in section V.2. The efficiency of TFCLUSTERS is compared to that of the (unrealistic) ideal power detector (which assumes a knowledge of the signal duration and central frequency) as a function of the signal-to-noise ratio of the signal; at fixed probability of detection, the reduction in signal-to-noise ratio for TFCLUSTERS is consistently less than about 30%.

Numerical simulations (section VI) confirm that the analytical model for the false alarm rate is accurate in the regime of operation relevant to TFCLUSTERS, with errors on the predicted false alarm rates around 1% in most situations, and with maximum errors around 10%, due to the exclusion of higher-order terms in the theoretical modeling, or to unmodeled finite size effects.

## II Transient Detection

The transient detection problem consists in choosing between two possibilities: the observed data $`y(t)`$, $`0t𝒯`$, consist of a signal term $`s(t)`$ and a noise term $`n(t)`$, where the noise is additive and is assumed to be white Gaussian with zero mean and unit variance whitening :

$$y(t)=s(t)+ϵn(t),$$
(1)
or the data are noise alone:

$$y(t)=ϵn(t).$$
(2)
In the simplest case, $`s(t)`$ is fixed and known. Perhaps the most natural optimality criterion is then to maximize the probability of detection when the signal is present for the probability of detection when no signal is present being smaller than some preassigned false alarm probability. This so-called Neyman-Pearson criterion leads to the use of the amplitude $`\widehat{Q}`$ of the signal, optimally estimated Fisher\_estimation  from the scalar product

$$\widehat{Q}=y(t),s(t),$$
(3)
as the statistic on which the decision between Eq. (1) and Eq. (2) is made.

In the case where $`s(t)W`$ for some function space $`W`$, one can integrate the signal distribution to reduce the problem to the simple case of discerning between two fixed probability distributions, as is done when $`s(t)`$ is fixed and known. This requires the knowledge of the prior distribution $`p[s(t)]`$, and results in using the likelihood ratio

$$\mathrm{\Lambda }[y(t)]=\frac{_{s(t)W}p[y(t)|s(t)]p[s(t)]𝑑s(t)}{p[y(t)|0]}$$
(4)
as the detection statistic. When the prior $`p[s(t)]`$ is unknown, the choice of an optimality criterion is not as simple as for the fixed signal case; the modified Neyman-Pearson criterion maximizes the minimum of the probability of detection over all $`s(t)W`$, for the false alarm probability being below some preassigned value, and is therefore an interesting “conservative” choice. It is well-known that the optimal algorithm in terms of the modified Neyman-Pearson criterion is obtained by choosing the prior $`p[s(t)]`$ that is least-favorable, i.e. the prior that minimizes the minimum over $`W`$ of the probability of detection when the signal is present, at fixed false alarm rate LeastFavorable . In particular, using the generalized likelihood ratio

$$\mathrm{\Lambda }^{}[y(t)]=\underset{s(t)W}{\mathrm{max}}\frac{p[y(t)|s(t)]}{p[y(t)|0]}$$
(5)
as the detection statistic does not guarantee optimality GLRT .

For signals with excellent theoretical models, such as binary neutron star or black hole inspirals, the function space $`W`$ is compact enough that the least-favorable prior is approximately uniform in $`s(t)`$, so that the generalized likelihood ratio \[Eq. (5)\] derives from Eq. (4) as the optimal detection statistic. This leads to a convenient implementation Owen  which simply thresholds on the maximum of the correlation \[Eq. (3)\] over a filter bank, which is a discretely sampled version of the function space $`W`$.

For the case where the signal space is not simple enough to allow matched filtering, that is when the signal is not well-modeled, a number of incoherent methods have been proposed in the literature. One of them is the so-called excess power detector epower , which basically threshold on the power integrated over a large number of different shapes at different positions in the time-frequency plane. In their discussion of the optimality of the excess power detector, the authors of epower  use a prior uniform in the “whitened” signal subspace of all waveforms of finite duration and bandwidth, which is imposed from physical arguments, and their result is therefore only a proof of optimality with respect to that particular prior. In particular, this prior is not necessarily least-favorable, and their result is therefore not a proof of optimality in the modified Neyman-Pearson sense. Similarly, the author of vicere  chooses a different prior, which is uniform in the “unwhitened” signal space, to derive a detector which is optimal with respect to this prior, and which is similar to the excess power detector, but is perhaps better adapted to colored instrumental noise.

In addition to these detectors, a number of ad hoc methods have been proposed to solve the transient detection problem for unmodeled transients. They are based on the analysis of patterns in the time-frequency plane tfplane ; tfplane1 , or on the time-domain analysis of the data with various filters tdomain . Only the authors of tdomain  discuss optimality, from a numerical point of view, i.e. by using a small number of simulated signals from ZM  that are injected into noise in numerical simulations in order to compare the performances of a few different detection algorithms.

When the noise is Gaussian, the likelihood ratio in Eq.(4) can be rewritten as kailath :

$$\mathrm{\Lambda }[y(t)]=\mathrm{exp}\left[_0^𝒯y(t)\widehat{s}(t)dt\frac{1}{2}_0^𝒯\widehat{s}^2(t)𝑑t\right],$$
(6)
where $`\widehat{s}(t)`$ is the causal minimum mean-square error estimator of $`s(t)`$, i.e.,

$$\widehat{s}(t)=E[s(t)|y(\tau )],$$
(7)
for $`E[|y(\tau )]`$ the expectation over the noise given the observation $`y(t)`$ for $`\tau <t`$, assuming the model of Eq.(1). The symbol $``$ represents the Itô stochastic integral Ito . Of course, the computation of $`\widehat{s}(t)`$ in Eq.(7) requires the knowledge of the prior $`p[s(t)]`$, but the structure of Eq.(6) suggests that an efficient approach to problems for which the integration in Eq.(4) cannot be carried out effectively (due to a complex or unknown prior, for instance) might be to develop efficient estimators of the signal, and to use them as filters to test for the presence of signals in the data (i.e., to use an estimator-correlator design).

It is shown in Donoho  that (i) transforming the (discretely sampled) observations $`\{y_i:i=1,\mathrm{},N\}`$ to a wavelet basis, (ii) applying on the transformed observations $`\{\stackrel{˘}{y}_i\}`$ the non-linearity

$$\stackrel{ˇ}{s}_i=\{\begin{array}{cc}0\hfill & \text{if }|\stackrel{˘}{y}_i|<\eta \hfill \\ \stackrel{˘}{y}_i\eta \hfill & \text{if }\stackrel{˘}{y}_i>\eta \hfill \\ \stackrel{˘}{y}_i+\eta \hfill & \text{if }\stackrel{˘}{y}_i<\eta ,\hfill \end{array}$$
(8)
for a threshold $`\eta ϵ\sqrt{2\mathrm{log}(N)/N}`$, and (iii) transforming the truncated observations $`\{\stackrel{ˇ}{s}_i\}`$ back to the time domain lead to an estimator $`\{\widehat{s}_i\}`$ that is optimal in the sense of giving the smallest value of the maximum of the expected least-square error over a wide class of signals. Using this estimator and preserving the structure of Eq.(6), a simple quantity to threshold on is

$$\underset{i=1}{\overset{N}{}}\left(y_i\widehat{s}_i\frac{\widehat{s}_i^2}{2}\right),$$
(9)
which can be rewritten as

$$\underset{i=1}{\overset{N}{}}\left(\frac{\stackrel{ˇ}{s}_i^2}{2}+\eta |\stackrel{ˇ}{s}_i|\right).$$
(10)
The algorithmic structure of TFCLUSTERS corresponds to this model; the two principal differences are that a short-time Fourier basis is used instead of a wavelet basis, and that an additional step involving the analysis of the correlations between the non-zero components of the estimator $`\stackrel{ˇ}{s}`$ is introduced.

## III Algorithm

It will be assumed that $`N`$ data corresponding to the gravitational wave strain are read and stored as real numbers into a vector $`𝒚`$. The data are time ordered and uniformly sampled with sampling frequency $`f_s`$, so the $`i^{\mathrm{th}}`$ component $`y_i`$ of the vector $`𝒚`$ is the value of the measured strain at time $`i/f_s`$: $`y_i=y(i/f_s)`$.

Step (i) of TFCLUSTERS is the construction of the time-frequency representation of the data from the spectrogram:

$$P_{ij}(𝒚)=|\stackrel{~}{y}_{ij}|^2,$$
(11)
where $`\stackrel{~}{y}_{ij}`$ is the $`j^{\mathrm{th}}`$ component of the discrete Fourier transform of the $`i^{\mathrm{th}}`$ segment of data of length $`M`$, $`𝒚_i=\{y_k:k=1+(i1)M,\mathrm{},iM\}`$. $`M`$ is assumed to be even and to be a factor of $`N`$, so that $`i=1,\mathrm{},N/M`$, and $`j=1,\mathrm{},M/2+1`$. $`j=1`$ and $`j=M/2+1`$ correspond to DC and to the Nyquist frequencies, and will not be considered anywhere below. Consequently, the maximum number of useful pixels in the time-frequency representation will be $`N_s=(M/21)(N/M)`$. The time resolution $`T`$ of the time-frequency representation is fixed and is given by $`T=M/f_s`$. The frequency resolution $`F`$ is simply $`F=1/T`$. In this time-frequency representation, the two hypotheses expressed in Eqs. (1) and (2) become $`P_{ij}(𝒚)=|\stackrel{~}{s}_{ij}+ϵ\stackrel{~}{n}_{ij}|^2`$ and $`P_{ij}(𝒚)=|ϵ\stackrel{~}{n}_{ij}|^2`$, respectively.

If the noise $`𝒏`$ is Gaussian and white, the power in different pixels is statistically independent. When no signal is present, $`P_{ij}`$ is the sum of the square of two independent Gaussian variables with zero mean and equal variance, and is therefore exponentially distributed:

$$p_{P_{ij}}(P)=\frac{e^{P/ϵ^2}}{ϵ^2}.$$
(12)
When a signal is present, the Gaussian variables have non-zero mean, and the probability density function (hereafter, pdf) of the power is spectra :

$$p_{P_{ij}}(P||\stackrel{~}{s}_{ij}|^2)=\frac{e^{(P+|\stackrel{~}{s}_{ij}|^2)/ϵ^2}I_0\left(2\sqrt{P|\stackrel{~}{s}_{ij}|^2}/ϵ^2\right)}{ϵ^2},$$
(13)
where $`I_0`$ is the modified Bessel function of order zero.

The spectrogram representation is of course not the only time-frequency representation available, and may not be optimal for most signal. It is however the simplest one to work with in this exploratory work. Wavelet bases share the independence property of the noise in the different pixels, as do any orthonormal basis, and they are known to offer better localization properties than Fourier transforms for many signals WavDonoho . However, their dyadic representation makes their analysis more complicated, due to the varying shape of the pixels with scale. Another classical way to improve the spectrogram representation is to use windows and to overlap the segments. This is a good way to reduce artifacts from Fourier transforms (e.g. edge effects) and to increase the time resolution of the time-frequency representation, but at the cost of destroying the statistical independence of the pixels. This independence is essential for the calculations involved in the clustering analysis presented in section V, and for the correct interpretation of some of the TFCLUSTERS thresholds. Any practical implementation of TFCLUSTERS can, however, deal equally well with spectrograms built with windows and overlapping segments as it does with the simpler type of spectrogram described above.

Step (ii) of TFCLUSTERS consists in applying a threshold on the power $`P_{ij}`$. Pixels with $`P_{ij}>\eta `$ are called black pixels, and other pixels are called white pixels. The probability for any given pixel to be black when no signal is present (the black pixel probability, $`p`$), is given by:

$$p=\mathrm{exp}(\eta /ϵ^2).$$
(14)
Figure 1 illustrates the result from Step (ii) and the effects of Step (iii) and (iv) below on simulated data.

Step (iii) of TFCLUSTERS considers the clustering of the black pixels. A cluster is defined as a set of black pixels containing all the black pixels that are the nearest neighbour of any pixel in the set. For a given pixel, its nearest neighbours are the pixels immediately to its left and to its right (time steps immediately before and after), and above and below it (frequency difference equal to the spectrogram resolution $`F`$). Two pixels touching only by a “corner” are called next nearest neighbours. The size $`S`$ of a cluster is simply defined as the number of black pixels it contains. The notion of distance $`d_c`$ between two clusters $`\mathrm{\Gamma }_1`$ and $`\mathrm{\Gamma }_2`$ is defined as the minimal distance $`d_p`$ between any two pixels $`p_1`$ and $`p_2`$ in the two clusters,

$$d_p(p_1,p_2)=|i_1i_2|+|j_1j_2|,$$
(15)
$$d_c(\mathrm{\Gamma }_1,\mathrm{\Gamma }_2)=\underset{p_1\mathrm{\Gamma }_1,p_2\mathrm{\Gamma }_2}{\mathrm{min}}d_p(p_1,p_2),$$
(16)
where $`(i_k,j_k)`$ are the coordinates of $`p_k`$, i.e., $`p_1`$ corresponds to $`P_{i_1j_1}(𝒚)`$, etc. Hence, nearest neighbour pixels have $`d_p=1`$, next-nearest neighbour pixels have $`d_p=2`$, and any two clusters must have $`d_c2`$. This choice of distance is made for convenience, but it has the implication of making the distance isotropic in the spectrogram representation of the time-frequency domain, irrespectively of the actual values of its time and frequency resolutions. Building a spectrogram with very long time bins and therefore very narrow frequency bins would have the effect of making the distances “longer” in the frequency direction than a spectrogram with short time bins and wide frequency bins.

Thresholds are applied both on the size of the cluster and on its distance to other clusters. The latter is easily motivated for physical signals; for example, although clusters of size two that are produced by fluctuations in the noise could be likely in a certain experiment, and therefore be below the size threshold, to have two or more such clusters close from each other, say being next nearest neighbours, could be rather unlikely. Hence, a signal with a well defined frequency that is slowly varying, so that its time-frequency track is a thin curve, could easily produce an archipelago of clusters of size two that is statistically unlikely to be produced by noise fluctuations alone.

If a cluster $`\mathrm{\Gamma }_1`$ has size $`S_1\sigma `$, it is immediately accepted as significant. Otherwise, its distance to all other clusters $`\mathrm{\Gamma }_2`$ with size $`S_2<\sigma `$ is compared to a threshold $`\delta _{S_1,S_2}`$ which depends explicitly on the size of the two clusters. All the clusters with $`d_c(\mathrm{\Gamma }_1,\mathrm{\Gamma }_2)\delta _{S_1,S_2}`$ are merged into a generalized cluster, which is declared significant. If $`\mathrm{\Gamma }_2`$ is already in a separate generalized cluster, $`\mathrm{\Gamma }_1`$ is added to that generalized cluster. If no cluster $`\mathrm{\Gamma }_2`$ with $`S_2<\sigma `$ satisfies the distance criterion, $`\mathrm{\Gamma }_1`$ is rejected. For a given choice of the minimum cluster size $`\sigma `$, there are $`\sigma (\sigma 1)/2`$ distance thresholds $`\delta _{S_1,S_2}`$, which will be organized below as a vector $`𝜹`$:

$$𝜹=[\delta _{S_1,S_2}]=[\delta _{1,1},\delta _{1,2},\mathrm{},\delta _{1,\sigma 1},\delta _{2,2},\mathrm{},\delta _{2,\sigma 1},\mathrm{},\delta _{\sigma 1,\sigma 1}].$$
(17)
Finally, step (iv) of TFCLUSTERS considers significant clusters from step (iii) and threshold on their excess power $`\widehat{Q}`$, which is defined as

$$\widehat{Q}=\underset{(i,j)\mathrm{\Gamma }}{}(P_{ij}(𝒚)\eta ),$$
(18)
for any given cluster $`\mathrm{\Gamma }`$ of size $`S`$. If no signal is present, i.e. if the cluster $`\mathrm{\Gamma }`$ results from a fluctuation of the noise, the pdf of $`P_{ij}(𝒚)`$ after the thresholding of step (ii) will be a truncated exponential for $`(i,j)\mathrm{\Gamma }`$:

$$p_{P_{ij}(𝒚)}(P)=\{\begin{array}{cc}\frac{e^{(P\eta )/ϵ^2}}{ϵ^2}\hfill & \text{if }P>\eta \hfill \\ 0\hfill & \text{otherwise.}\hfill \end{array}$$
(19)
The pdf of $`\widehat{Q}`$ will be the convolution of $`S`$ such distributions, i.e.

$$p_{\widehat{Q}}(P)=\{\begin{array}{cc}\frac{1}{(S1)!}\left(\frac{PS\eta }{ϵ^2}\right)^{S1}\frac{e^{(PS\eta )/ϵ^2}}{ϵ^2}\hfill & \text{if }P>S\eta \hfill \\ 0\hfill & \text{otherwise.}\hfill \end{array}$$
(20)
Hence, setting a size-dependent threshold $`𝒬_S`$ on $`\widehat{Q}`$, defined by the integral equation

$$\alpha =_{S\eta }^{𝒬_S}\frac{1}{(S1)!}\left(\frac{PS\eta }{ϵ^2}\right)^{S1}e^{(PS\eta )/ϵ^2}\frac{dP}{ϵ^2},$$
(21)
leads when only noise is present to the rejection of a fraction $`\alpha `$ of the clusters that survived step (iii), independently of the cluster size $`S`$.

Step (iv) is very similar to the scheme used in the excess power algorithm of epower . The essential difference is that TFCLUSTERS chooses the pixels included in the computation of the excess power by using a threshold on the individual pixel power and a clustering analysis, while epower  use a fixed set of “masks” of different shapes that are translated in the time-frequency plane, and over which the power for all pixels is integrated. The definition of these masks requires some prior expectations about the signal characteristics, and when this prior information is not sufficient to constrain the signal shapes to look for in the time-frequency plane, the selective approach of TFCLUSTERS is expected to become more efficient than the excess power approach of epower .

It should also be noted that other kinds of thresholds on the total cluster power could be used in step (iv). While the one presented here penalizes in a similar fashion clusters of all sizes (it reduces the false rate by $`\alpha `$ independently of $`S`$), a size independent threshold could also be used, for instance, in which case large clusters would be more likely to survive the last step of TFCLUSTERS.

## IV Operating Characteristics

Assessing the optimality of TFCLUSTERS in the modified Neyman-Pearson sense is a very involved mathematical task, as it can be inferred from similar problems treated in the literature Igster . Nevertheless, it is shown in this section that the algorithm presents interesting properties that suggest its near-optimality for a large class of signals. What will actually be discussed is a simplified version of the algorithm not involving the clustering analysis (step (iii) of TFCLUSTERS). It is assumed that the clustering analysis will improve the performances of this detector when signals that form clusters in the time-frequency plane are indeed present.

The binary test will be constructed by comparing the signal power estimate, $`\widehat{Q}(𝒚)=|\widehat{𝒔}|^2`$, to a certain threshold $`\zeta `$; $`\widehat{Q}>\zeta `$ will lead to the acceptance of Eq. (1). This estimator is constructed by summing the power in the spectrogram over the pixels that have power above a certain threshold:

$`|\widehat{\stackrel{~}{s}}_{ij}|^2=\{\begin{array}{cc}0\hfill & \text{if }P_{ij}(𝒚)<\eta \hfill \\ P_{ij}(𝒚)\eta \hfill & \text{otherwise}\hfill \end{array}`$ (24)
$`(P_{ij}(𝒚)\eta )_+,`$ (25)
$$\widehat{Q}(𝒚)=\underset{i,j}{}|\widehat{\stackrel{~}{s}}_{ij}|^2.$$
(26)
The sum in Eq. (26) is performed over the whole time-frequency plane, and therefore is over $`N/2`$ terms.

The following theorem, inspired from the work of Donoho  on signal estimation, is proved in appendix A:
Theorem 1:
Given the model Eqs. (24)-(26), for $`Q=|𝒔|^2`$, and provided that $`\eta =\beta ϵ^2\mathrm{log}N/2`$, $`\beta >1`$, there exists a series of numbers $`\pi _N`$ with $`\pi _N1`$ for $`N1`$ such that $`𝒔𝖱^N`$:
(i) $`Pr(\widehat{Q}Q)=\pi _N`$,
(ii) $`Pr(\widehat{Q}\widehat{q})=\pi _N`$, where $`\widehat{q}`$ is any power estimator with $`Pr(\widehat{q}Q^{})=\pi _N`$ $`𝒔^{}\mathrm{\Omega }(𝒔)`$, where $`\mathrm{\Omega }(𝒔)`$ is some neighbourhood of $`𝒔`$ defined in Eq. (48), and $`Q^{}=|𝒔^{}|^2`$.
(iii) $`Pr\left(\widehat{Q}Q_{i,j}\mathrm{min}(2\beta ϵ^2\mathrm{log}N/2,|\stackrel{~}{s}_{ij}|^2)\right)=\pi _N`$.

This theorem leads to a number of observations:

* Optimality: Theorem 1 shows the optimality of $`\widehat{Q}`$ in the subset of detectors for which the condition stated in part (i) is respected: $`\widehat{Q}`$ provides the tightest lower bound to $`Q`$ that can be constructed from the data, for all $`𝒔`$.
* False alarm probability: For a threshold $`\zeta `$ on $`\widehat{Q}`$, part (i) implies that the probability to detect a signal with power $`Q<\zeta `$ goes to zero asymptotically with $`N`$ getting large. In particular, part (i) implies that the false alarm probability goes to zero asymptotically, since when no signal is present $`Q=0`$.
* False dismissal probability: The direct interpretation of part (iii) is that any signal with $`Q>\zeta +_{i,j}\mathrm{min}(2\beta ϵ^2\mathrm{log}N/2,|\stackrel{~}{s}_{ij}|^2)`$, or equivalently with $`_{i,j}\left(|\stackrel{~}{s}_{ij}|^22\beta ϵ^2\mathrm{log}N/2\right)_+>\zeta `$, will be detected with a false dismissal probability approaching zero. Obviously, if $`\mathrm{max}|\stackrel{~}{s}_{ij}|^22\beta ϵ^2\mathrm{log}N/2`$, part (iii) implies that $`\widehat{Q}=0`$, giving a limit on how sensitive this detector can be.
* Sparse signals: Also as a consequence of part (iii), sparse signals that have only a small number of non-zero $`\stackrel{~}{s}_{ij}`$ are likely to be more easily detected, because they have a smaller value of $`_{i,j}\mathrm{min}(2\beta ϵ^2\mathrm{log}N/2,|\stackrel{~}{s}_{ij}|^2)`$ than signals with the same power spread over a larger number of pixels.

The scheme discussed above of thresholding on the power $`P_{ij}(𝒚)`$ integrated over pixels that have power above a certain threshold $`\eta `$ thus possesses certain optimal properties. It is very efficient for signals that have a sparse representation in the time-frequency domain. It is unlikely, however, that the optimal properties are preserved when the signal subspace is restricted to signals that form clusters in the time-frequency domain; it that case, intelligence about the spatial correlation of the signal pixels should certainly allow for more efficient algorithms. The idea behind TFCLUSTERS is that the ad hoc approach of merging the two-threshold scheme discussed above with a clustering analysis designed to effectively reject noise should lead to an efficient algorithm for such a signal subspace.

## V Clustering Analysis

On physical grounds, it can be expected for most transient sources that the pixels with excess power due to the signal will tend to cluster in the time-frequency domain. Short signals like black hole ringdowns or mergers, for instance, have durations roughly of the order of the inverse of the usable interferometer bandwidths ($`1`$ kHz), and therefore appear as connected clusters of duration $`T`$ and bandwidth equal to the search bandwidth. For longer signals that spend tens or hundreds of cycles in the interferometers band, the required stability of the source is likely to necessitate that the dominating mechanism for the emission of the waves is governed by rotation, with the wave instantaneous frequency being around $`J/\pi I`$, for $`J`$ the magnitude of the angular momentum along the principal axis of rotation, and $`I`$ the moment of inertia of the source. Complicated dynamics may lead to the spreading of the signal power over some finite frequency interval, or even the formation of sidebands well-separated from the principal frequency of the waves (in this case, each sideband will be considered an independent cluster). When the source is not too far away, it is argued below that the amount of angular momentum radiated by gravitational waves for sources that are near the sensitivity limit of the interferometric detectors is insufficient to cause a change in the rotation frequency that is rapid enough to produce disconnected pixels in the time-frequency representation of the signal.

Over a time interval $`T`$ corresponding to one time slice in the time-frequency representation, a source at a distance $`r`$ from the Earth that emits uniformly in all directions waves of frequency $`f`$ and characteristic strain amplitude $`h`$ will radiate a total amount of energy $`\mathrm{\Delta }E`$, where Thorne

$$\mathrm{\Delta }E\pi ^2\frac{c^3}{G}f^2h^2r^2T.$$
(27)
The change $`\mathrm{\Delta }J`$ in angular momentum magnitude corresponding to the emission of the waves is related to the amount of radiated energy by

$$|\mathrm{\Delta }J|\frac{\mathrm{\Delta }E}{\pi f}$$
(28)
if most of the radiation is quadrupolar angrel . On the other hand, for sources with dynamics dominated by rotation, the second time derivative of the mass quadrupole moment is bounded by $`\ddot{Q}\pi Jf`$, so for $`h=2G\ddot{Q}/rc^4`$,

$$J\frac{c^4}{2\pi G}\frac{hr}{f}.$$
(29)
The value of the characteristic amplitude $`h`$ is expressed in term of the noise power spectral density $`S_n(f)`$ so that the signal-to-noise ratio for a signal with bandwidth $`1/T`$ is unity, corresponding to a marginally detectable source:

$$h=S_n^{1/2}(f)/T^{1/2}.$$
(30)
Combining Eqs. (27)-(30):

$$\frac{|\mathrm{\Delta }J|}{J}310^4\left(\frac{r}{1\mathrm{M}\mathrm{p}\mathrm{c}}\right)\left(\frac{f}{100\mathrm{H}\mathrm{z}}\right)^2\left(\frac{S_n^{1/2}(f)}{510^{23}\mathrm{Hz}^{1/2}}\right)\left(\frac{T}{0.1\mathrm{s}}\right)^{1/2},$$
(31)
where the numbers in Eq. (31) for $`S_n(f)`$ and $`f`$ correspond approximately to the minimum of the noise spectral density of the interferometers being presently developed, and the value of $`T`$ is chosen to match the expected time resolution of the time-frequency representations to be used on actual data. For values of $`r`$ that are sufficiently small, a source does not need to radiate gravitational waves at a rate that produce variations in its angular momentum magnitude that are significant in order to be detectable. In order of magnitude, $`|\mathrm{\Delta }J|/J|\mathrm{\Delta }f|/f`$, for $`\mathrm{\Delta }f`$ the change in the wave frequency over a time $`T`$, provided that the source doesn’t change its moment of inertia by a large fraction, which is unlikely to happen over timescales of 0.1 s, except perhaps near the end of the gravitational wave signal (as e.g. in binary inspirals, at the innermost stable circular orbit), at which point the question of the amount of clustering becomes irrelevant. Hence, the wave frequency is not expected to vary enough over the time resolution $`T`$ to generate pixels in contiguous time slices that are disconnected. This can happen if

$$\frac{|\mathrm{\Delta }f|}{f}\frac{1}{Tf}.$$
(32)
Combining Eqs. (31) and (32) gives a necessary (but not sufficient) condition for pixels to be disconnected:

$$r300\text{ Mpc}\left(\frac{100\text{ Hz}}{f}\right)^3\left(\frac{510^{23}\mathrm{Hz}^{1/2}}{S_n^{1/2}(f)}\right)\left(\frac{0.1\text{ s}}{T}\right)^{3/2}.$$
(33)
In this context, any rotation dominated source at the detection limit of the detector that is closer than $`r`$ is expected to form a cluster in the time-frequency domain. Figure 2 shows the variation of $`r`$ with frequency, for the noise spectral density of the LIGO interferometers psdLIGO . Sources with $`f300`$ Hz will form detectable clusters even if they are as far as the Virgo cluster; those with $`f1250`$ Hz, even if they are in the galaxy, may not form clusters.

These considerations show that under certain restrictions, at least two broad classes of sources (short impulsive and rotation dominated at the detection limit) should lead to signals that form clusters in the time-frequency plane. As it will be shown below, Gaussian noise will have the opposite property, in the sense that black pixels will tend to fill the plane uniformly, without forming large clusters. This will be used as a powerful basis for denoising the thresholded spectrograms computed from the data.

### V.1 False Alarm Rate

The reader familiar with the mathematics of percolation theory in statistical physics has most likely recognized at this point the applicability of the results in that field to the problem under investigation here. In particular, for a black pixel probability $`p`$, the average number of clusters of size $`S`$ per pixel of an infinite image is:

$$n_S(p)=p^SD_S(1p),$$
(34)
where $`D_S`$ is the so-called perimeter polynomial, and is related to the number of shapes a cluster of size $`S`$ can have (counting shapes related only by a translation as identical), that is to the degeneracy $`g_{SR}`$ of a cluster of size $`S`$ and perimeter $`R`$:

$$D_S(q)=\underset{R}{}g_{SR}q^R.$$
(35)
As it can be seen from Eq. (35), the perimeter of a cluster is the number of white pixels having at least one black pixel in the cluster as a nearest neighbour. Coefficients of the perimeter polynomials for cluster sizes up to $`S=22`$ for nearest neighbours on the square lattice are tabulated in Mertens , and were mostly generated through the use of computer enumeration techniques. Note that $`n_S(p)`$ is simply the probability of any pixel to be in a cluster of size $`S`$, divided by $`S`$ Conway .

For low cluster densities, the expected number of clusters per unit time $`\lambda `$, i.e. the cluster false alarm rate, is related for a threshold $`\sigma `$ on the cluster size to the time and frequency resolutions, and to the bandwidth $`B`$ that is searched:

$$\lambda =\frac{B}{TF}\underset{S=\sigma }{\overset{\mathrm{}}{}}n_S(p).$$
(36)
In practice, because of the rapid decay of $`n_S(p)`$ with $`S`$, the sum can easily be truncated without major losses of precision. Figure 3 illustrates the nearly exponential decay of the cluster rate with the threshold on the cluster size, for different black pixel probabilities.

In analogy to Eq. (34), the average number per pixel of pairs of clusters of size $`S_1`$ and $`S_2`$ separated by a distance $`d`$ is defined as:

$$\nu _{S_1,S_2}(d)=p^{S_1+S_2}H_{S_1S_2}^d(1p),$$
(37)
where the polynomial $`H_{S_1S_2}^d`$ is related to the number of configurations $`k_{S_1S_2R}(d)`$ for a cluster of size $`S_1`$ to be within a distance $`d`$ from a cluster of size $`S_2`$, where the sum of the perimeters of the two clusters is $`R`$, and where configurations that are related only by a translation are again considered identical:

$$H_{S_1S_2}^d(q)=\underset{R}{}k_{S_1S_2R}(d)q^R.$$
(38)
As an example, for two clusters of size 2 separated by a distance of two, Eq. (37) becomes:

$$\nu _{2,2}(2)=\frac{56p^4q^{11}+40p^4q^{10}}{4}.$$
(39)
The factors in the numerator account for all the possible configurations for a cluster of size 2 to be a distance 2 from another cluster of size 2, and the factor of 4 in the denominator is the total number of pixels, necessary in order not to overcount clusters related by a simple translation. The general expression for $`d>2`$ is

$$\nu _{2,2}(d)=8(d+1)p^4q^{12},$$
(40)
and hence the associated false alarm rate is

$`\lambda _{2,2}={\displaystyle \frac{B}{TF}}\left[10p^4q^{10}+14p^4q^{11}+{\displaystyle \underset{d=3}{\overset{\delta _{2,2}}{}}}8(d+1)p^4q^{12}\right]`$ (41)
$`={\displaystyle \frac{B}{TF}}\left[10p^4q^{10}+14p^4q^{11}+4(\delta _{2,2}^2+3\delta _{2,2}10)p^4q^{12}\right],\text{ if }\delta _{2,2}3,`$ (42)
where $`\delta _{2,2}`$ is the threshold on the distance for two clusters of size 2 to be considered correlated.

Table 1 gives the coefficients $`k_{S_1S_2R}(d)`$ for small cluster sizes that were obtained by computer enumeration. General relations can easily be deduced from this table: on the square lattice, using the definition of distance from Eq. (15), the diameter of a circle of radius $`r`$ is $`4r`$, and therefore the number of points at distance $`d`$ from a certain cluster increases linearly with $`d`$. Moreover, when $`d>2`$, both clusters are guaranteed to not be sharing any perimeter white pixels, and therefore the number of configurations with fixed total perimeter consists of a “bulk” part (from the constant number of configurations occurring at pixels that are on the same row or column as a pixel in the cluster), and of a part growing linearly with $`d`$ (from the “diagonals” of the curve of constant distance from the first cluster). Table 2 gives the formula for $`\nu _{S_1,S_2}(d)`$ deduced from table 1, for $`d>2`$. Values for $`d=2`$ can be read directly from table 1. Figure 4 shows how the expected number of pairs of clusters of size $`S_1`$ and $`S_2`$ (i.e., $`_{d=2}^{\delta _{S_1,S_2}}\nu _{S_1,S_2}(d)`$) increases when the threshold $`\delta _{S_1,S_2}`$ is increased. As expected, it also shows that pairs of large clusters are more unlikely to form than pairs of smaller clusters, and that the expected number of pairs of clusters decreases when the black pixel probability is reduced.

It should be noted that TFCLUSTERS includes higher order terms ($`n`$-cluster configurations, $`n3`$) into its definition of generalized clusters. These complex configurations are built by merging all the clusters satisfying the distance criteria, and consequently the sum of Eq. (36) and of its equivalent for Eq. (37) overestimates slightly the true false alarm rate, especially when the distance thresholds $`\delta _{S_1,S_2}`$ are large for the density of small clusters.

### V.2 Efficiency

As shown above, the false alarm rate is computed analytically; when $`\alpha =0`$ (cf. step (iv) of TFCLUSTERS), the probability of detection of a certain signal can also be computed analytically. It is not possible at this time to model the detection efficiency when $`\alpha >0`$ because the equivalent to Eq. (20) when a non-null signal is present is not known analytically, and is required to compute the probability of detection of a given signal. The sensitivity and false alarm rate of TFCLUSTERS depend on a number of parameters: $`F,T,p,\sigma `$ and $`𝜹`$. Using some general information about the expected waveform, it is possible to run an optimization over all these parameters to maximize the probability of detection at fixed false alarm rate. More generally, it is possible to run multiple versions of the detector, each one for a different set of parameters, in order to cover as many different classes of signal as possible. The number of classes that can be considered is limited by the computational power available for the analysis, and is somewhat limited by the fact that results are not independent.

For comparison, the probability of detection using the ideal power detector can also be computed. This detector assumes that the bandwidth and duration of the signal are known and that their product is the time-frequency volume $`V`$, and it compares the excess power computed according to these parameters to a threshold $`\kappa `$. In order to get a false alarm rate rather than a false alarm probability, it is assumed that the detector is applied with frequency $`1/\tau `$ on segments of independent data of length $`\tau `$, and that all frequency bins in the search bandwidth $`B`$ are considered in disjoint groups having bandwidth $`\varphi `$. This detector is incoherent and is therefore not as efficient as matched filtering, but it still performs better than any detector that can be implemented when neither the duration nor the bandwidth is known. The false alarm rate is:

$$\lambda =\frac{B}{\varphi \tau }_\kappa ^{\mathrm{}}\frac{(P/ϵ^2)^{V1}e^{P/ϵ^2}}{ϵ^2(V1)!}𝑑P,$$
(43)
and the probability of detection is spectra

$$\beta =_\kappa ^{\mathrm{}}\left(\frac{P}{V\rho ^2ϵ^2}\right)^{(V1)/2}\mathrm{exp}\left(\frac{P+V\rho ^2ϵ^2}{ϵ^2}\right)I_{V1}\left(2\sqrt{\frac{V\rho ^2P}{ϵ^2}}\right)\frac{dP}{ϵ^2}.$$
(44)
The case where one or more of the bandwidth, duration, central frequency or arrival time parameters have to be searched over is more complicated to analyze: the probabilities are not statistically independent, because of the need to search over overlapping regions in the time-frequency representation of the signal, and therefore the computation of the false alarm rate is rather involved.

An example is now presented to make things more precise regarding how well TFCLUSTERS compares to the ideal power detector. The signal is taken to be of duration $`6T`$ and bandwidth smaller than $`F=1/T`$ (since the problem is symmetric under the interchange of the time and frequency axes, the results below are equivalent to the case of a short signal of duration $`T`$ and bandwidth $`6/T`$). It is assumed that the signal power is distributed uniformly over its full duration and bandwidth. As a definite example, let the bandwidth of the search be chosen to be $`B=512F`$, and the false alarm rate be $`\lambda =1/3600T`$. The central frequency is assumed to be such that the power is concentrated in a single constant frequency row of the spectrogram, while both the case of synchronized (the signal covers six pixels) and the case of random arrival time (the signal covers seven pixels) are considered. For the ideal power detector, $`V=6`$ or 7, and $`\tau `$ is set to $`6T`$. For the TFCLUSTERS algorithm, the various thresholds ($`\eta `$, $`\sigma `$ and $`\delta _{S_1,S_2}`$, $`\alpha `$ being zero) are optimized for every value of the signal-to-noise ratio in order to maximize the probability of detection (hereafter, POD) for the constraint $`\lambda 1/3600T`$. The details of the calculation are presented in appendix B. Figure 5 shows a comparison of the optimized POD as a function of signal to noise ratio $`\rho `$, for TFCLUSTERS and the ideal power detector. Only two sets of clustering analysis thresholds cover the whole range of signal-to-noise ratios of figure 5, illustrating the relative independence of the performances of TFCLUSTERS on its numerous parameters. Given the fact that the chosen signal was a line in the time-frequency plane, and was therefore the shape that is the easiest to “break” by changing a black pixel into a white pixel, the example presented here can be considered as a difficult situation; most other distributions of the power (i.e. other signal shapes) would make the performances of TFCLUSTERS and of the ideal power detectors only closer.

## VI Numerical Simulations

Extensive numerical simulations were carried out in order to confirm the validity of the analyses presented in the previous sections, and in order to explore properties of TFCLUSTERS that are hard to study analytically. All of the results below were produced using the same method: segments of Gaussian white noise of unit variance are produced using a random number generator RAN2 . The segments are of duration $`10^4`$ s, and are sampled at 2048 Hz. These numbers are chosen so that processing a single segment uses most of the RAM of the computers on which the code is running. On every segment, the implementation of TFCLUSTERS within the LIGO Scientific Collaboration Algorithm Library LAL  is used to generate a list of significant clusters, according to some pre-specified values of the various thresholds. The 90% central confidence interval for the rate of significant clusters in the data, assuming they form a Poisson process, is computed based on the number of detected clusters, using the standard Neyman construction Neyman . If the ratio of the width of the 90% confidence interval to its central value is smaller than 1%, the simulation terminates. Otherwise, a number of other $`10^4`$ s segments are processed, until the above termination criterion is met. For a true rate of $`\lambda `$, the termination criterion requires to process of order $`10^5/\lambda `$ worth of simulated data.

It should be noted that the timing distribution of clusters from white noise is indeed very well approximated by a Poisson distribution. This is confirmed in one specific case by figure 6 which shows that the distribution of the time delay $`\mathrm{\Delta }`$ between two successive clusters follows an exponential distribution, which is the defining property of a Poisson process.

The simulations were performed on a cluster of workstations with 1GHz Intel Pentium III “Coppermine” processors, 512 Mbytes of PC-133 RAM memory, with the Linux operating system. For a large range of parameter values, the processing time per CPU was on average 250 to 550 times shorter than the duration of the data segment. Most of the time was spent at grouping the black pixels into clusters, and quite logically the most important factor in determining the speed of TFCLUSTERS was the black pixel probability; independently of the other parameters, the ratio of processing to real time was around 300 for $`\eta /ϵ^2=2`$, around 500 for $`\eta /ϵ^2=4`$, and was increasing almost linearly with the power threshold $`\eta `$.

### VI.1 The Case $`𝜹=\mathrm{𝟎}`$, $`\alpha =0`$

When $`𝜹=\mathrm{𝟎}`$ and $`\alpha =0`$, the expected rate is given by Eq. (36). Figures 7 and 8 show the excellent agreement between rates from simulations and predictions from Eq. (36). The two agree to better than 0.5% most of the time, commensurably with the precision of the simulations. The sum in Eq. (36) is of course dominated by the first few terms, but values of $`n_S(p)`$ as predicted from Eq. (34) describe the simulated data also very well for large values of $`S`$, as shown in figure 9.

### VI.2 The Case $`𝜹\mathrm{𝟎}`$, $`\alpha =0`$

For $`𝜹\mathrm{𝟎}`$, there is a contribution to the cluster rate from both Eq. (36) and its equivalent for $`\nu _{S_1,S_2}(d)`$. Figures 10 and 11 show the good agreement between rates from simulations and predictions. When $`\eta /ϵ^22.8`$, the agreement is at the precision level of the simulations, although there is a systematic overestimation of the measured rate by the predictions for $`\eta /ϵ^23.7`$. This overestimation, however, reaches almost 20% for $`\eta /ϵ^2=2`$; this is expected, as explained in section V, because for $`\eta /ϵ^22.8`$ the cluster density is high enough that higher-order combinations above the 2-cluster one are likely to be produced.

To demonstrate that the error is indeed due to higher-order terms, a histogram of the contribution $`\lambda (S)`$ to the total rate of generalized clusters of size $`S`$ is presented in figure 12, together with the theoretical prediction based on 1- and 2-cluster configurations. The predictions systematically overestimate the measured rate for $`S<8`$, and underestimate it for $`S8`$. It should be noted that for $`𝜹=[0,0,0,0,0,0,2,3,4,4]`$, 3-cluster configurations can have sizes $`8S12`$, while 4-cluster configurations can have sizes $`10S16`$. Hence, one can expect a small error for $`S=8`$ corresponding to the 3-cluster $`(S_1,S_2,S_3)=(2,2,4)`$, and no errors from higher-order terms for $`S<8`$, i.e. for the terms dominating the sum leading to the prediction for $`\lambda `$. Overall, more small clusters than expected get merged into generalized clusters due to these high-order terms, and the total rate is consequently overestimated. Figure 13 quantifies the importance of high-order terms as a function of the cluster density, i.e. as a function of the threshold $`\eta `$; as expected, the relative importance of high-order terms becomes smaller than the 1% level around $`\eta /ϵ^2=3`$, in agreement with the results in figure 10.

### VI.3 The Case $`\alpha >0`$, and Finite-Size Effects

For $`\alpha >0`$, the rate is expected to be $`(1\alpha )\lambda _0`$, where $`\lambda _0`$ is the value of the rate for the same parameters, except that $`\alpha =0`$. Figure 14 illustrates the reduction in $`\lambda `$ as $`\alpha `$ is varied. The results are as expected from Eq. (21), within the errors from the simulations.

The simulations presented in figure 14 were carried out with a very asymmetrical spectrogram: $`T`$ was chosen to be 1/32 s, so $`F=32`$Hz, and the 992 Hz bandwidth was covered by only 31 pixels, while a segment had $`3210^4`$ bins in time. Because of that, significant finite size effects are expected; the rate predictions are based on the assumption that the time-frequency plane is infinite, but in the present case it is small enough in the frequency dimension that clusters are likely to be “clipped” and therefore to be smaller then expected. Consequently, it is expected that the predictions will again overestimate the measured rate, and this is what is observed in figures 14 and 15. As it can be seen from figure 15, the fractional error from the prediction is still smaller than about 10% when the cluster density is low enough.

## VII Conclusions

TFCLUSTERS, a new time-frequency detector for bursts of gravitational radiation in broad-band interferometric observatories, was described in some details in this paper. The behavior of the detector when applied to white Gaussian noise in the absence of signals was carefully analyzed, leading to a formalism for the computation of the false alarm rate of TFCLUSTERS for any values of its many parameters. The results from numerical simulations showed that this analysis was accurate in most situations: errors at the 1% level or better were obtained in “ideal” situations (low cluster density, large number of frequency bins), and errors of the order of 10% appeared in situations were the analysis was expected to be less accurate, due to high order terms not included in the sums over the cluster configurations, or due to finite size effects from the limited bandwidth of the search. In the case where the errors were large, the analysis presented in this paper was systematically overestimating the false alarm rates; should they prove from practical work to be necessary, more accurate estimates could therefore be obtained from more careful calculations, or from larger scale numerical simulations. While the false alarm properties of TFCLUSTERS are well-understood, the efficiency of the detector is subject to more uncertainties.

A calculation that was presented in Appendix B showed that the efficiency of TFCLUSTERS was comparable to that of the ideal power detector, unfortunately illustrating at the same time the mathematical complexity associated with producing such an estimate of the efficiency of the detector for a given signal. Nevertheless, the fundamental approach used by TFCLUSTERS, namely the use of an adaptive power integral over pixels with excessive power, was shown to be maximizing asymptotically the estimate of the power in the signal, provided it is not overestimating it, over all possible estimators of the same quantity. This naturally suggests optimal properties for TFCLUSTERS as a detector. However, this can be nothing more than a conjecture at this point, as the actual proof of the optimality of TFCLUSTERS in the (modified) Neyman-Pearson sense is of great difficulty.

Independently of the question of the optimality of TFCLUSTERS, the structure of the algorithm is particularly practical for its implementation for the analysis of actual data. The first power threshold on individual pixels (step (ii)) can be chosen to be frequency dependent in order to allow the analysis of colored Gaussian noise when whitening of the data is not convenient; the errors on the black pixel probability that are introduced by the non-zero correlations coloring the noise are generally negligible. Moreover, frequency bands containing spurious interferences are easily left out of the analysis, with minor modifications to the algorithm presented in this paper.

The search for gravitational waves will of course require the operation of TFCLUSTERS in coincidence at sites that are geographically separated and the use of information from auxiliary data that do not couple to the gravitational wave data, in order to avoid possible false detections due to non-Gaussian noise. A simple but complete analysis system using TFCLUSTERS and satisfying these requirements was developed in thesis  in order to compute upper limits on the rate of occurrence of events that are expected to radiate gravitational radiation. A comparison of event lists from the individual detectors was used to carry out this coincidence analysis. This approach may, however, not give the best network detection efficiency, in part because of the relatively coarse time resolution of TFCLUSTERS. Enhanced versions of the algorithm, designed to analyze coherently data from many detectors, might prove to be more efficient.

###### Acknowledgements.

The author wishes to thank Rainer Weiss for helpful discussions and comments on the work presented here, and the LSC group at the University of Wisconsin-Milwaukee for generous time allocation on their computer system. This work was supported by the National Science Foundation under cooperative agreement PHY-9210038 and the award PHY-9801158, and by the Fonds pour la Formation de Chercheurs et l’Aide à la Recherche of the Province of Québec, Canada.

## Appendix A Proof of theorem 1

Consider first the deterministic equivalent to Eq. (1):

$$𝒚=𝒔+\delta 𝒖,$$
(45)
where $`𝒖`$ is a nuisance parameter, so that the equivalent to Eq. (11) is

$$P_{ij}(𝒚)=|\stackrel{~}{s}_{ij}+\delta \stackrel{~}{u}_{ij}|^2,$$
(46)
where $`\stackrel{~}{u}_{ij}𝖢`$ and where by definition $`|\stackrel{~}{u}_{ij}|1`$.

Lemma 1:
For the model described by Eq. (45) with $`\delta ^2=\eta `$, $`𝒔𝖱^N`$ and $`𝒚`$ satisfying Eq. (45), $`\widehat{Q}(𝒚)Q`$.
Proof:
Consider Eq. (24). Trivially, if $`|\widehat{\stackrel{~}{s}}_{ij}|^2=0`$, then $`|\widehat{\stackrel{~}{s}}_{ij}|^2|\stackrel{~}{s}_{ij}|^2`$. For $`|\widehat{\stackrel{~}{s}}_{ij}|^2>0`$, $`|\widehat{\stackrel{~}{s}}_{ij}|^2=P_{ij}(𝒚)\eta `$. From Eq. (46), $`|\stackrel{~}{s}_{ij}|^2P_{ij}(𝒚)\delta ^2`$. Since $`\delta ^2=\eta `$, this gives $`|\stackrel{~}{s}_{ij}|^2|\widehat{\stackrel{~}{s}}_{ij}|^2`$ for all $`𝒔,𝒚`$. Summing over $`i`$ and $`j`$ gives $`\widehat{Q}Q`$.

Lemma 2:
Given Eq. (45) with $`\delta ^2=\eta `$, $`𝒔𝖱^N`$, and $`𝒚`$ respecting Eq. (45), $`\widehat{Q}(𝒚)\widehat{q}(𝒚)`$, where $`\widehat{q}(𝒚)`$ is any power estimator satisfying

$$\widehat{q}(𝒚^{})Q,𝒔^{}\mathrm{\Omega }(𝒔)\text{ and }𝒚^{}\text{ respecting Eq. (}\text{45}\text{)},$$
(47)
where

$$\mathrm{\Omega }(𝒔)=\{𝒔^{}𝖱^N:|𝒔𝒔^{}|^2\underset{ij}{}\mathrm{min}(4\delta ^2,2\delta |\stackrel{~}{s}_{ij}|\left[1+\frac{|\stackrel{~}{s}_{ij}|}{\delta }\left(1+\sqrt{1+\frac{2\delta }{|\stackrel{~}{s}_{ij}|}}\right)\right])\}.$$
(48)
Proof:
Suppose $`𝒚^0`$ such that $`\widehat{q}(𝒚^0)>\widehat{Q}(𝒚^0)`$, so $`\stackrel{~}{y}_{ij}^0=\stackrel{~}{s}_{ij}^0+\delta \stackrel{~}{u}_{ij}^0`$. One can construct a signal $`𝒔^{}`$ such that

$$\stackrel{~}{s}_{ij}^{}=\{\begin{array}{cc}\stackrel{~}{y}_{ij}^0\delta \stackrel{~}{u}_{ij}^{}\hfill & \text{if }|\stackrel{~}{y}_{ij}^0|>\delta \hfill \\ 0\hfill & \text{otherwise}.\hfill \end{array}$$
(49)
The freedom provided by $`\stackrel{~}{u}_{ij}^{}𝖢`$ is used to choose all values of $`\stackrel{~}{u}_{ij}^{}`$ so that $`\stackrel{~}{u}_{ij}^{}`$ and $`\stackrel{~}{s}_{ij}^{}`$ are orthogonal. If $`\stackrel{~}{s}_{ij}^{}=0`$, the direction of $`\stackrel{~}{u}_{ij}^{}`$ is unimportant. For all $`i,j`$, $`|\stackrel{~}{u}_{ij}^{}|=1`$. Note that this choice of $`\stackrel{~}{u}_{ij}^{}`$ and $`\stackrel{~}{s}_{ij}^{}`$ satisfies Eq. (45).

Also, since $`\stackrel{~}{s}_{ij}^{}\stackrel{~}{s}_{ij}^0=\delta (\stackrel{~}{u}_{ij}^0\stackrel{~}{u}_{ij}^{})`$, $`|\stackrel{~}{s}_{ij}^{}\stackrel{~}{s}_{ij}^0|^24\delta ^2`$. Now, $`|\stackrel{~}{y}_{ij}^0|^2=\delta ^2+|\stackrel{~}{s}_{ij}^{}|^2`$, so,

$`|\stackrel{~}{s}_{ij}^{}\stackrel{~}{s}_{ij}^0|^2=|\stackrel{~}{s}_{ij}^0|^2+|\stackrel{~}{s}_{ij}^{}|^22\stackrel{~}{s}_{ij}^0,\stackrel{~}{s}_{ij}^{}`$ (50)
$`2|\stackrel{~}{s}_{ij}^0|^2+2|\stackrel{~}{s}_{ij}^0|\delta +2|\stackrel{~}{s}_{ij}^0|^2\sqrt{1+{\displaystyle \frac{2\delta }{|\stackrel{~}{s}_{ij}^0|}}}.`$ (51)
Combining these inequalities and summing over $`i,j`$ shows that $`𝒔^{}\mathrm{\Omega }(𝒔)`$.

Then,

$`\widehat{Q}(𝒚^0)={\displaystyle \underset{i,j}{}}(P_{ij}(𝒚^0)\eta )_+`$ (52)
$`={\displaystyle \underset{i,j}{}}(|\stackrel{~}{s}_{ij}^{}+\delta \stackrel{~}{u}_{ij}^{}|^2\eta )_+`$ (53)
$`={\displaystyle \underset{i,j}{}}(|\stackrel{~}{s}_{ij}^{}|^2+\delta ^2\eta )`$ (54)
$`=|𝒔^{}|^2.`$ (55)
It follows that $`\widehat{q}(𝒚^0)>\widehat{Q}(𝒚^0)`$ implies $`\widehat{q}(𝒚^0)>|𝒔^{}|^2`$, contradicting Eq. (47).

Lemma 3:
For the model described by Eq. (45) with $`\delta ^2=\eta `$, $`𝒔𝖱^N`$ and $`𝒚`$ satisfying Eq. (45), $`\widehat{Q}(𝒚)Q_{i,j}\mathrm{min}(2\delta ^2,|\stackrel{~}{s}_{ij}|^2)`$.
Proof:
Obviously, $`|\stackrel{~}{s}_{ij}|^2|\widehat{\stackrel{~}{s}}_{ij}|^2|\stackrel{~}{s}_{ij}|^2`$. Also, from Eq. (46), $`P_{ij}(𝒚)|\stackrel{~}{s}_{ij}|^2\delta ^2`$, so from Eq. (24), $`|\widehat{\stackrel{~}{s}}_{ij}|^2=P_{ij}(𝒚)\delta ^2|\stackrel{~}{s}_{ij}|^22\delta ^2`$ when $`|\widehat{\stackrel{~}{s}}_{ij}|^2>0`$. When $`|\widehat{\stackrel{~}{s}}_{ij}|^2=0`$, $`\delta ^2>P_{ij}(𝒚)|\stackrel{~}{s}_{ij}|^2\delta ^2`$, so $`|\stackrel{~}{s}_{ij}|^2<2\delta ^2`$. Combining the three inequalities gives $`|\stackrel{~}{s}_{ij}|^2|\widehat{\stackrel{~}{s}}_{ij}|^2\mathrm{min}(2\delta ^2,|\stackrel{~}{s}_{ij}|^2)`$; summing over $`i,j`$ proves the lemma.

Proof of theorem 1:
In the noise model from Eq. (11), the power in the noise is distributed exponentially:

$$p_{|\stackrel{~}{n}_{ij}|^2}(P)=\frac{e^{P/\sigma ^2}}{ϵ^2}.$$
(56)
The corresponding cumulative probability distribution is

$$P_{|\stackrel{~}{n}_{ij}|^2}(P)=1e^{P/ϵ^2}.$$
(57)
Hence, the probability that the maximum of $`N`$ independent realizations of that random variable is less than a threshold $`C_N`$ is

$$\pi _NPr\left(\underset{i,j}{\mathrm{max}}|\stackrel{~}{n}_{ij}|^2<C_N\right)=\left(1e^{C_N/ϵ^2}\right)^N.$$
(58)
For $`N1`$, this can be rewritten as

$$\pi _N=\mathrm{exp}\left(e^{\mathrm{log}NC_N/ϵ^2}\right)$$
(59)
for $`C_N/ϵ^2>\mathrm{log}N`$. Hence, for $`C_Nϵ^2\mathrm{log}N`$, this probability goes to $`1/e`$ as $`N`$ becomes large. Moreover, any values of $`\pi _N>1/e`$ can be achieved with the right choice of $`C_N`$; for $`C_N=\beta ϵ^2\mathrm{log}N`$, $`\beta >1`$, $`\pi _N1`$ for $`N1`$.

It should be noted that the event

$$\underset{i,j}{\mathrm{max}}|\stackrel{~}{n}_{ij}|^2<\beta ϵ^2\mathrm{log}N$$
(60)
implies a realization of the noise which can be mimicked by the deterministic noise model \[Eq. (45)\] with $`\delta =ϵ\sqrt{\beta \mathrm{log}N}`$. Hence, the three statements of theorem 1 follow from the three lemmata proved in this appendix.

## Appendix B The Narrow-Band Signal Example

Given a $`u`$ by $`v`$ matrix $`𝑸`$ with elements representing the distribution of the power in the signal in a rectangular sub-region of the time-frequency plane, the signal-to-noise ratio is defined by:

$$\rho ^2=\frac{_{i,j}Q_{ij}}{uvϵ^2}.$$
(61)
The elements of the matrix $`𝑫`$ representing the black pixel probability corresponding to the matrix $`𝑸`$ are given by the integral of the density from Eq. (13):

$$D_{ij}=_\eta ^{\mathrm{}}p_{P_{ij}}(P|Q_{ij})𝑑P.$$
(62)
In general, not all pixels where some signal power is present will be black, and the signal will be detected only when a number of black pixels equal or larger than the threshold $`\sigma `$ form a connected cluster, or when a pair of smaller clusters are close enough. The contribution to the probability of detection of the signal of such a configuration will be the product of the black pixel probabilities ($`D_{ij}`$) and of the white pixel probabilities ($`1D_{ij}`$) for “holes”. Although noise fluctuations could in principle help the detection of a signal by forming a “bridge” over regions where no signal is present, summing over these possibilities involve $`2^{uv}`$ terms. A slight underestimate of the probability of detection is instead used by summing only over the $`n`$ pixels where signal is present ($`nuv`$), which reduces the number of terms to be considered to

$$\underset{n_H=0}{\overset{n\sigma }{}}\frac{n!}{(nn_H)!n_H!},$$
(63)
where $`n_H`$ is the number of holes. Of course, the enumeration process can be greatly simplified when the distribution of the power in the signal has some specific symmetries.

Consider now the example of section V.2. Under the assumption that the starting time of the signal matches exactly the binning of the spectrogram used to detect it, it will be represented by a row of 6 pixels, each with an equal probability $`p`$ to be black, neglecting effects such as power leakage. The columns labeled “POD ($`rp,s0`$)” in tables 3 and 4 give the contributions to the probability of detection (POD) of this signal from various thresholds. The problem is slightly more complex when the arrival time is taken to be random. In that case, the signal is spread in general over 7 pixels, with the central five having a black pixel probability $`p`$, and the leftmost and rightmost having smaller probabilities $`r`$ and $`s`$, respectively. $`r`$ and $`s`$ are given by:

$$r=_\eta ^{\mathrm{}}p_{P_{ij}}(P|P_s^{})𝑑P$$
(64)
and

$$s=_\eta ^{\mathrm{}}p_{P_{ij}}(P|P_sP_s^{})𝑑P,$$
(65)
where $`P_s`$ is the power in the central five pixels, $`P_s^{}`$ is the power in the leftmost pixel, and $`p_{P_{ij}}(P|Q)`$ is given by Eq. (13). The POD is then given by

$$\frac{1}{P_s}_0^{P_s}\mathrm{POD}(r(P_s^{}),s(P_sP_s^{}))𝑑P_s^{},$$
(66)
where $`\mathrm{POD}(r(P_s^{}),s(P_sP_s^{}))`$ is the probability of detection for the 7 pixels configuration; contributions to it from various thresholds can be found in the columns labeled “POD” in tables 3 and 4. For the ideal power detector with $`\lambda =1/3600T`$, $`\tau =6T`$ and $`B/\varphi =512`$, Eq. (43) gives $`\kappa 23.87ϵ^2`$ for $`V=uv=6`$ and $`\kappa 25.55ϵ^2`$ for $`V=uv=7`$. Eq. (44) is then used directly to generate the ideal power detector curves of figure 5.


format: text

text:
# Origin of non-keplerian motions of masers in NGC 1068

## 1 Introduction

Even though direct images are missing, evidences that Active Galactic Nuclei (AGN) host a super-massive black hole surrounded by a parsec-size accretion disc are compelling. In a few objects, the spatial resolution is just sufficient to unveil the outermost regions of the central engine at the sub-parsec scale thanks to the detection of water maser emission (Greenhill, 2002). The best example is NGC 4258 where masers have been associated with a thin disc in quasi-perfect keplerian rotation (Miyoshi al. 1995). In contrast, the azimuthal velocity of masers in NGC 1068 varies with the radius $`R`$ as $`R^{0.31\pm 0.02}`$ according to Greenhill et al. (1996) who mentioned the possible role of the disc/torus self-gravity. To our knowledge, this hypothesis has never been explored, even qualitatively. We show in section 2 that a sub-keplerian velocity exponent is theoretically expected close to the outer edge of a (low mass) disc when considering self-gravity. In section 3, we report possible parameters for the disc and for the black hole in NGC 1068 such that the computed velocity match the observed motion of maser spots. A few remarks and possible improvements are found in the last section.

## 2 Basic ideas

### 2.1 Gravity, orbiting velocity and velocity exponent

We consider the steady state system111In the following, quantities labeled with “bh” refer to the black hole, and those labeled with “disc” refer to the disc. schematically pictured in Fig. 1, made of a black hole with mass $`M_{\mathrm{bh}}`$ and Schwarzschild radius $`R_\mathrm{S}=\frac{2GM_{\mathrm{bh}}}{c^2}`$ surrounded by an axi-symmetrical, finite size gaseous disc with inner edge $`R_{\mathrm{in}}`$, outer edge $`R_{\mathrm{out}}`$, semi-thickness $`H`$, surface density $`\mathrm{\Sigma }`$ and mass $`M_{\mathrm{disc}}=qM_{\mathrm{bh}}`$. Assuming that gravity balances the centrifugal force in the disc, then the rotation velocity $`v_\varphi `$ of gas at cylindrical distance $`R`$ from the center obeys the Euler equation

$$v_\varphi ^2=Rg_R(>0),$$
(1)
where $`g_R=g_R^{\mathrm{bh}}+g_R^{\mathrm{disc}}`$ is the radial component of the total (black hole + disc) gravitational acceleration. Note that we neglect the motion of the system around its center of mass (although possibly significant) and consider the black hole as “fixed”. Further, the mass distribution in the disc is such that $`g_R<0`$ at the outer edge. We can define the total “velocity exponent” $`k`$ as

$$k\frac{\mathrm{ln}v_\varphi }{\mathrm{ln}R}=\frac{1}{2}\left[1+\frac{\mathrm{ln}(g_R)}{\mathrm{ln}R}\right],$$
(2)
so that the azimuthal velocity varies as $`R^k`$, locally. In particular, for a black hole only, one finds the keplerian value $`k^{\mathrm{bh}}=1\frac{3R^2}{2r^2}`$ where $`r=\sqrt{R^2+z^2}`$ is the spherical radius (and $`k^{\mathrm{bh}}=\frac{1}{2}`$ at the equatorial plane $`z=0`$).

### 2.2 Combining a black hole and a disc: edge effects

There is no reliable, regular expression222It is worth noting that any expression for $`g_R^{\mathrm{disc}}`$ established assuming a zero-thickness disc inevitably diverges at the edges (Mestel, 1963) and must not be used in the disc inside and in the neighborhood as well. This is why we shall employ here a numerical tri-dimensional approach which provides everywhere with a regular self-gravity field.  for the radial field $`g_R^{\mathrm{disc}}`$ due to a disc (as a tri-dimensional object) with given shape and density distribution, and thus no straightforward expression of the disc velocity exponent $`k^{\mathrm{disc}}`$ is available. For the present system however, the relation

$$k\frac{1}{2}=\left[k^{\mathrm{bh}}\frac{1}{2}\right]\frac{g_R^{\mathrm{bh}}}{g_R}+\left[k^{\mathrm{disc}}\frac{1}{2}\right]\frac{g_R^{\mathrm{disc}}}{g_R}$$
(3)
can easily be established from Eq.(2). We thus see that a departure to the keplerian motion (that is $`kk^{\mathrm{bh}}`$) can arise in two opposite situations: i) when the disc self-gravity exceeds the gravity due to the black hole, or ii) when the disc self-gravity is much smaller than the central gravity, but $`|k^{\mathrm{disc}}|1`$. In the extreme, the first case suggests that a non keplerian rotation can be produced by a disc only (without a black hole; i.e. $`q\mathrm{}`$). Here, we shall adopt the second point of view which rather corresponds to configurations with a mass ratio $`q1`$.

The point is that, depending on the mass distribution in a disc, large values of $`k^{\mathrm{disc}}`$ can occur at the edges, and especially at the outer edge of interest here, even when $`q1`$. The physical reason is that radial gravity is generally vanishingly small well in the middle of a disc because the gravitational contributions of matter located interior and exterior to any point essentially cancel out (unless very steep density gradients). Conversely, near an edge, there creates an imbalance between these two contributions, and the resulting field can have a quite large amplitude. Hence a gradient of radial self-gravity. An illustration of this ”edge effect” is given in Fig. 2 which displays the radial field and velocity exponent in a system with mass ratio $`q=0.1`$. In this example, we have adopted parameters typical for an AGN standard disc (e.g. Frank, King & Raine, 1992), namely black hole mass $`M_{\mathrm{bh}}=10^8`$ M, disc inner edge $`R_{\mathrm{in}}=3R_\mathrm{S}`$ and outer edge $`R_{\mathrm{out}}=10^4R_\mathrm{S}`$, aspect ratio $`\frac{H}{R}=0.01\times f(R)`$ and surface density $`\mathrm{\Sigma }R^{3/10}\times f(R)`$, where we have introduced the $`f`$-function in order to make the disc thickness and density gradually decrease to zero over a disc scale height typically (i.e. $`f1`$ over almost all the disc extent and $`f=0`$ outside $`[R_{\mathrm{in}},R_{\mathrm{out}}]`$). This is necessary to avoid a sharp discontinuity between the disc and the ambient medium which would not be realistic at all. Further, we impose that the vertical stratification of matter between the disc mid-plane and its surface is quadratic with the altitude. This conveniently mimics the Gaussian profile expected in a vertically isothermal disc as well as the solution of the plane Lane-Emden equation for a vertically self-gravitating disc (e.g. Ibañez & Sigalotti, 1984). Quite importantly, we point out that both the adopted prescription for $`f(R)`$ and the assumption regarding the vertical stratification have no noticeable influences on the global results within physically acceptable limits, so that issues raised here are not artefacts due to discontinuity effects. Finally, the radial gravitational acceleration due to the disc is computed numerically using the accurate Poisson 3D-solver described in Huré (2002).

We see that, although the disc is “light” with respect to the black hole and central gravity fully dominates, $`k`$ deviate markedly from the keplerian value near the outer edge. A general property (met here) is that the rotation is characterized by a “sub-keplerian velocity exponent” for $`RR_{\mathrm{out}}`$ in the sense that $`k\frac{1}{2}`$, whereas this is just the opposite for $`RR_{\mathrm{out}}`$.

## 3 The case of NGC 1068

### 3.1 Origin for the non-keplerian motions

NGC 1068 is probably one of the best studied333See for instance the proceedings of the workshop published in the ’97 special issue of Astrophysics and Space Science. among type-2 Seyfert galaxies. In this object with systemic velocity of $`1126`$ km/s and distance of $`15`$ Mpc, water maser emission has been detected at the parsec scale. Some emission features would originate from the surface of a disc, along a diameter perpendicular to the line of sight, in the range $`0.651.1`$ pc (Greenhill et al. 1996). This corresponds to a masing disc extending from $`9`$ to $`15`$ milli-arcseconds (mas) in the sky. From the rotation curve of water maser spots, the dynamical mass enclosed within $`1.1`$ pc is estimated to $`1.5\times 10^7`$ M and the best fit for the velocity exponent yields $`k_{\mathrm{obs}.}=0.31\pm 0.02`$.

The significant departure to Kepler law means that this dynamical mass may not be that of the black hole. Greenhill et al. (1996) suggested that disc self-gravity could be responsible for this non-keplerian motion. Two other, plausible explanations have be defended: radiation pressure effects (Pier & Krolik, 1992) and extra gravitational attraction by a stellar cluster (Kumar, 1999). However, the disc self-gravity hypothesis is quite natural, at least because NGC 1068, as an active galaxy, hosts a black hole and a disc with a certain mass. It is also corroborated by the arguments exposed in §2.2. From this point of view, the condition $`k_{\mathrm{obs}.}>k^{\mathrm{bh}}`$ (i.e. a sub-keplerian exponent) implies that the maser spots in this AGN rotate inside the disc (see Fig. 2), and not outside (which would mean a ”super-keplerian” exponent). We can thus conclude that the outer edge of the masing disc is larger than $`1.1`$ pc.

### 3.2 “Inverting” the rotation curve

Matching observational data means that one must find the right magnitude of the orbiting velocity, and not only the velocity exponent. As $`v_\varphi `$ depends on the distribution of matter in space (see Eq. (1)), we see that the problem is to determine, if possible, which mass distribution(s) can cause the observed motions — a common problem in Astrophysics. More formally, this means the inversion of some relation

$$v_\varphi v_\varphi (R,z;M_{\mathrm{bh}},R_{\mathrm{in}},R_{\mathrm{out}},\mathrm{\Sigma },H,\mathrm{}),$$
(4)
which obviously involves several parameters. In such a problem, not a single solution does exist but many, and finding all solutions would require the exploration of the whole parameter space. This task is complex and well out of the scope of this paper. However, preliminary investigations have revealed that i) the location of the inner edge has no influence provided $`R_{\mathrm{in}}R_{\mathrm{out}}`$, and it is set to $`3R_\mathrm{S}`$, ii) the disc thickness is not a critical parameter, provided the disc is geometrically thin (namely $`\frac{H^2}{R^2}1`$), and iii) the outer edge should stand beyond $`1.1`$ pc (see §3.1).

To simplify more the quest of solution, we further assume an uniform aspect ratio $`\frac{H}{R}`$, with $`0.3`$ as a firm upper limit (Gallimore et al., 1996), and we impose that the surface density is of the form $`\mathrm{\Sigma }=\mathrm{\Sigma }_0R^s`$ (we keep on using the smoothing $`f`$-function; see §2.2). Besides, we set $`z=\pm H`$, since we are interested in the velocity field at the disc surface where masers are assumed to be excited. Then, remains to be constrained mainly the central mass, the surface density and the location of the outer edge.

### 3.3 Properties of the outer disc and central mass

Two extreme solutions will be discussed here: a thick disc solution with the maximum value for aspect ratio, and a thin disc solution. Table 1 lists the parameters deduced for the disc and for the black hole. As Fig. 3 shows, both models fit very well observations, despite the total velocity index $`k`$ is far from being constant (but is of order of $`0.3`$ on average). It turns out that i) a single value for the black hole mass is required, ii) $`\mathrm{\Sigma }`$ decreases roughly as $`R^1`$, iii) the disc outer edge is larger than $`1.3`$ pc (or $`18`$ mas), and iv) the disc mass is close the black hole mass. Note that the data reported here concern the outer disc and must not extrapolated down to the central black hole, but only over one or two decade(s) in radius.

These results deserve two comments. First, the physical conditions required for the emission of water maser seem met (Neufeld, Maloney & Conger, 1994). Actually, at $`0.65`$ pc, adopting $`\mathrm{\Sigma }380`$ g/cm2 as a mean value (see Tab. 1), the typical number density is

$$n\frac{\mathrm{\Sigma }}{2H\mu m_\mathrm{H}}\frac{3\times 10^7}{\left(\frac{H}{R}\right)}\mathrm{cm}^3,$$
(5)
where $`\mu `$ is the mean mass per particle ($`2.3`$ in the molecular phase for cosmic abundances) in units of the mass of hydrogen $`m_\mathrm{H}`$. This is in the right density range (i.e. $`10^810^{10}`$ cm-3) provided the aspect ratio of the disc satisfies $`3\times 10^3\frac{H}{R}0.3`$. Further, for such densities and for temperatures appropriate for molecules to survive (i.e. $`3000`$ K), support by radiation pressure must be significant, even dominant (whatever the solution). This conclusion is akin to the analysis by Pier & Krolik (1992). Assuming a vertically self-gravitating disc in hydrostatic equilibrium, the mid-plane temperature is

$$T_{\mathrm{mid}}\left(\frac{3\pi Gc\mathrm{\Sigma }^2}{4\sigma }\right)^{1/4}2000\mathrm{K},$$
(6)
where $`\sigma `$ is the Stephan constant and $`c`$ is the speed of light. For a typical grain opacity of $`1`$ cm2/g, the vertical optical depth is $`\tau 200`$, implying a temperature at the disc surface of $`T_{\mathrm{mid}}/\tau ^{1/4}600`$ K, as required (Neufeld, Maloney & Conger, 1994).

The second remark concerns the gravitational stability of the disc. For both solutions, the Toomre $`Q`$-parameter is of the order of unity. Actually, again at $`0.65`$ pc, we have

$$Q\frac{\kappa c_\mathrm{s}}{\pi G\mathrm{\Sigma }}\frac{v_\varphi T_{\mathrm{mid}}^2}{\pi GR}\sqrt{\frac{8\sigma H}{3c\mathrm{\Sigma }^3}}2,$$
(7)
where $`\kappa \mathrm{\Omega }`$ is the epicyclic frequency and $`c_\mathrm{s}`$ is the sound speed. It means that the outer disc is in a marginally stable state, and does not resemble a standard disc (Kumar, 1999). This property is predicted by the standard model at the parsec scale, whatever the accretion rate in the inner disc. Note also that the fact that our solutions have $`\mathrm{\Sigma }R^1`$ and $`v_\varphi ^{\mathrm{disc}}cst`$ (see Fig. 3) is fully consistent with what we expect in an (radially) isothermal $`Qcst`$ disc (e.g. Mestel, 1963).

## 4 Concluding remarks

We have reminded that a non-keplerian motion can arise from the combination of the gravitational attraction of a black hole and a disc. We propose that the “sub-keplerian” profile traced by water masers in NGC 1068 results from such an effect, a possibility outlined already in Greenhill et al. (1996). We have reported two simple physical solutions producing approximately the observed mean velocity exponent, with, accordingly, the parameters for the black hole and for the outer disc in this active galaxy. As outlined, these solutions are not unique; for instance, it is possible to fit the feature near 12 mas (see Fig. 3), but $`\mathrm{\Sigma }`$ is no more a power law. From this point of view, it would be worthwhile to explore the whole parameter space.

The analysis discussed here can be applied to other AGN. Note however that the system considered here is very simple and would need refinements. In particular, there is no physical model for the disc and the present approach works with the most minimal conditions. For instance, Eq.(1) should be supplemented with non-gravitational terms (pressure gradient, advection, viscosity, etc.). Further, the gravitational attraction, possibly significant, of other components of the AGN (molecular torus, stellar cluster, galactic bulge, etc.) is not taken into account. Significant deviations to axi-symmetry (like a warp) could also change our results, quantitatively, and would make the problem much more tricky. With these restrictions, we have shown that it is feasible to “invert” the rotation curve in order to constrain the mass distribution in the central engine at the parsec scale. With more sophisticated or self-consistent investigations (e.g. Lodato & Bertin, 2002), one should be able to connect these outer solutions to the innermost ones where accretion theories are very uncertain.

###### Acknowledgements.

I am grateful to S. Collin, C. Boisson, F. Le Petit, and G. Lodato for stimulating discussions. I acknowledge James M. Moran, as the referee, for valuable comments and suggestions.


text:
# Evaluating matrix elements relevant to some Lorenz violating operators

## I introduction

In the recent literature, there have been considered a number of ways to modify the structure of space-time which can have experimental consequences. In one of the most popular scenarios, space-time is considered to become noncommutative at short distance scales, with space-time coordinates satisfying a commutation relation of the following form J.L. Hewett, F.J. Petriello, and T.G. Rizzo (2001); I. Hinchliffe and N. Kersting (2001); B. Jurco, L. Moller, S. Schraml, P. Schupp, and J. Wess (2001); Carlson et al. (2001); S. Godfrey and M.A. Doncheski (2002); A. Anisimov et al. (2002); Carlson et al. (2002a)

$$[\widehat{x}^\mu ,\widehat{x}^\nu ]=i\theta ^{\mu \nu },$$
(1)
where $`\widehat{x}^\mu `$ is a position four-vector promoted to an operator, and $`\theta ^{\mu \nu }`$ is a set of c-numbers antisymmetric in their indexes. The most striking effects of space-time non commutativity of the form (1) are the Lorentz violating effects appearing in field theories, which is a consequence of $`\theta ^{0i}`$ and $`\epsilon ^{ijk}\theta ^{ıȷ}`$ defining preferred directions in a given Lorentz frame.

Jurc̆o, Möller, Schraml, Schupp and Wess B. Jurco, L. Moller, S. Schraml, P. Schupp, and J. Wess (2001) have shown how to construct non-Abelian gauge theories in noncommutative spaces from a consistency relation. Using the same approach Carlson, Carone and Lebed Carlson et al. (2001) have derived the Feynman rules for consistent formulation of noncommutative QCD and they have computed the most dangerous, Lorentz-violating operator that is generated through radiative corrections. They have found that at the lowest order in perturbation theory, the formulation of noncommutative QCD that they have presented leads to Lorentz violating operators such as A. Anisimov et al. (2002)

$$\theta ^{\mu \nu }\overline{q}\sigma _{\mu \nu }q,\theta ^{\mu \nu }\overline{q}\sigma _{\mu \nu }\text{/}Dq\mathrm{and}\theta ^{\mu \nu }D_\mu \overline{q}\sigma _{\nu \rho }D^\rho q.$$
(2)
In Carlson et al. (2001) the phenomenological implications of the first of these operators were studied in detail. Noting that contributions from the space-space part of $`\theta ^{\mu \nu }`$ make $`\sigma _{\mu \nu }\theta ^{\mu \nu }`$ act like a $`\stackrel{}{\sigma }\stackrel{}{B}`$ interaction with $`\stackrel{}{B}`$ directly related to $`\theta ^{ij}`$, a limit was placed on the scale of non commutativity. One used the result of tests of Lorentz invariance in clock comparison experiments C.J. Berglund et al. (1995), which suggest that external $`\stackrel{}{\sigma }\stackrel{}{B}`$ like interactions are bounded at the $`10^7`$Hz level or few $`\times `$ $`10^{31}`$ GeV. Carlson et al. Carlson et al. (2001) concluded that

$$\theta \mathrm{\Lambda }^210^{29},$$
(3)
where $`\theta `$ is a typical scale for elements of the matrix $`\theta ^{\mu \nu }`$.

However, the effective Lorentz violating operator was obtained from a one loop correction to the quark propagator, and the operator proportional to $`\sigma _{\mu \nu }\theta ^{\mu \nu }`$ also contained a factor $`(\overline{)}pm)`$. With $`\stackrel{}{B}`$ constant, the evaluation of $`\stackrel{}{\sigma }\stackrel{}{B}`$ factors out from the evaluation of $`(\overline{)}pm)`$, and our discussion is focused on the later.

In Carlson et al. (2001) an ad hoc estimate was used for the matrix element of the operator $`(\overline{)}pm)`$, where $`m`$ is the current quark mass, in getting limit in equation (3). The matrix element $`\overline{)}pm`$ was estimated to be about $`M_N/3300`$ MeV, where $`M_N`$ is the nucleon mass. However, it has been argued that the expectation value of $`(\overline{)}pm)`$ could be much less than this naive estimate We thank M. Pospelov for discussion on this point. See also I. Mocioiu, M. Pospelov and R. Roiban (2002).

The aim of this paper is to calculate the matrix element of the operator $`(\overline{)}pm)`$, using variety of confinement potential models, so to evaluate the quality of the estimate made in Carlson et al. (2001).

The sample of potentials included four different confining potentials, two of them purely Lorentz scalar and two of them equal mixture of scalar and vector. The first scalar potential is a Bag-like potential

$$V(r)=\{\begin{array}{cc}V_0,\hfill & \text{if }rR\text{;}\hfill \\ 0,\hfill & \text{otherwise.}\hfill \end{array}$$
(4)
We also consider the one dimensional case for the nicety of the analytical result,

$$V(z)=\{\begin{array}{cc}V_0,\hfill & \text{if }z\frac{a}{2}\mathrm{or}z\frac{a}{2}\text{;}\hfill \\ 0,\hfill & \text{otherwise.}\hfill \end{array}$$
(5)
The $`V_0\mathrm{}`$ limit gives, of course, the MIT Bag model Chodos, A., R.L. Jaffe, K. Johnson, C.B Thorn, and V.F. Weisskopf (1974); Chodos, A., R.L. Jaffe, K. Johnson, and C.B Thorn, (1974) if one does not consider the Bag energy. We will consider models of vector+scalar confinement next, using in one case a linear spatial potential and on the other case a harmonic one,

$$V(r)=\frac{1}{2}(1+\gamma ^0)(V_0+\lambda r),$$
(6)
or

$$V(r)=\frac{1}{2}(1+\gamma ^0)Cr^2.$$
(7)
Finally we shall consider a purely scalar harmonic potential,

$$V(r)=Cr^2.$$
(8)
In the following sections it will be assumed that the current quark mass of $`510`$ MeV can be neglected compared to the quark eigenenergy of several hundred MeV.

## II Scalar Square-Well potential

For any given potential $`V`$, from the Dirac equation we have that

$$(\overline{)}pm)\psi =V\psi ,$$
(9)
therefore

$$\overline{)}pm=V.$$
(10)
In the three dimensional case, for the central potential $`V(r)`$ presented in (4), the solutions of the Dirac equation for the ground state, with $`m=0`$, in two regions $`𝐈.𝐫<𝐑,\text{and}\mathrm{𝐈𝐈}.𝐫>𝐑`$ (Fig. 1) have the following form

$$\psi _I(r)=N_I\left(\genfrac{}{}{0pt}{}{j_0(Er)}{i\sigma \widehat{r}j_1(Er)}\right)\chi ^{(s)},$$
(11)
$$\psi _{II}(r)=N_{II}\left(\genfrac{}{}{0pt}{}{h_0^{(1)}(ik_0r)}{\sigma \widehat{r}\sqrt{\frac{V_0E}{V_0+E}}h_1^{(1)}(ik_0r)}\right)\chi ^{(s)},$$
(12)
where $`k_0=\sqrt{V_0^2E^2}`$, and $`j_0,j_1`$ are the spherical Bessel functions, and $`h_0^{(1)},h_1^{(1)}`$ are the spherical Hankel functions of the first kind. The ground state energy can be found from the energy eigenvalue equation

$$j_1(ER)=j_0(ER)\left[\frac{1+k_0R}{(V_0+E)R}\right],$$
(13)
while for $`V_0\mathrm{}`$ the eigenvalue equation is $`j_1(ER)=j_0(ER)`$, as is familiar from the MIT Bag model Chodos, A., R.L. Jaffe, K. Johnson, C.B Thorn, and V.F. Weisskopf (1974); Chodos, A., R.L. Jaffe, K. Johnson, and C.B Thorn, (1974).

However, we know there are long range forces between baryons. If one wants to accommodate long range forces in this type of model, then one has to allow quarks to penetrate the walls of the potential well with some finite probability. Therefore the height of the potential, $`V_0`$, should be finite. A reasonable choice for $`V_0`$ and $`R`$ can be obtained by fitting the model parameters to obtain reasonable values, for example, for the mean square of charge radius of the nucleon $`r^2`$ and for the axial vector coupling constant $`g_A`$. We get a good fit by choosing $`R=1.12`$ fm and $`V_0=3`$ GeV for which we find $`r^2=0.64\mathrm{fm}^2`$ and $`g_A=1.15`$, as compared to experimental values of $`0.76\mathrm{fm}^2`$ and 1.27 respectively K. Hagiwara et al. (2002). Solving (13) for this choice of parameters for the ground state energy of a quark we find $`E=348`$ MeV.

Using solutions given in (11) and (12), we find

$$\overline{)}pm=21\text{MeV}.$$
(14)
Exploration of the integrals appearing in $`V(r)`$, shows that $`\overline{)}pm0`$ as $`1/V_0`$, when $`V_0\mathrm{}`$.

It may be of some pedagogic value to give the equivalent result for the 1D case (Fig. 2). The wave function for $`|z|<a/2`$ is just the free solution of the Dirac equation, and the solutions for $`|z|>a/2`$ are obtained from the free solution by the substitution $`EEV_0`$. We obtain

$$\overline{)}pm=2V_0_{\frac{a}{2}}^{\mathrm{}}\overline{\psi }\psi 𝑑z=\frac{E}{1+a\sqrt{V_0^2E^2}}.$$
(15)
One can note immediately that when the height of the potential $`V_0\mathrm{}`$ then $`\overline{)}pm0,`$ unless $`a0`$. For the choice of parameters made above, we obtain

$$\overline{)}pm=14\mathrm{MeV},$$
(16)
where for the ground state energy E we have used a value of $`260`$ MeV, from the energy eigenvalue equation.

## III Scalar + Vector Linear Confinement

Let us consider now the confinement problem of a spin 1/2 particle in a confining potential of the form

$$V(r)=\frac{1}{2}(1+\gamma ^0)(V_0+\lambda r).$$
(17)
This linear potential model for quark confinement was used in P.L. Ferreira, J.A. Helayel, and N. Zagury (1980) to calculate several properties of low-lying baryons. In P.L. Ferreira, J.A. Helayel, and N. Zagury (1980) the authors assumed nonzero quark masses. The straightforward modification of the wave functions for the case of vanishing current quark masses yields the following solution for the lowest energy eigenstate of the Dirac equation for the potential (17),

$`\mathrm{\Psi }(r)`$ $`=`$ $`N\left({\displaystyle \genfrac{}{}{0pt}{}{\mathrm{\Phi }(r)}{\sigma 𝐩/E\mathrm{\Phi }(r)}}\right)\chi ^{(s)},`$ (18)
$`\mathrm{\Phi }(r)`$ $`=`$ $`\sqrt{{\displaystyle \frac{K}{4\pi Ai^2(a1)}}}{\displaystyle \frac{1}{r}}Ai(Kr+a_1)),`$ (19)
where $`K=(\lambda E)^{1/3}`$. The energy eigenvalue $`E`$ and the normalization constant $`N`$ are given in (20)

$$E=V_0\frac{\lambda a_1}{K},N^2=\frac{3E}{4EV_0}.$$
(20)
In P.L. Ferreira, J.A. Helayel, and N. Zagury (1980) an analytic expression was obtained for the mean square charge radii of the barions and in P. L. Ferreira (1977) Ferreira obtained an analytic expression for the magnetic moment of the proton. We modified those expressions for the zero current quark mass case and used them together with the energy eigenvalue equation (20) to fit our model parameters $`V_0`$ and $`\lambda `$. We choose $`V_0=626`$ MeV and $`\lambda =0.98`$ GeV/fm to fit $`r^2`$ exactly and give the closest to the data value of $`\mu _p`$, obtaining

$$E=420\mathrm{M}\mathrm{e}\mathrm{V},r^2=0.76\mathrm{fm}^2\mathrm{and}\mu _p=2.44\mathrm{n}.\mathrm{m}.$$
(21)
For the above mentioned values of the model parameters we find that

$$\overline{)}pm=27\text{MeV}.$$
(22)
## IV Scalar + Vector Harmonic Confinement

Consider now potential of the form

$$V(r)=\frac{1}{2}(1+\gamma ^0)Cr^2.$$
(23)
The solution of Dirac equation with this potential is given in R. Tegen (1990). They write the lowest energy state Dirac spinor as

$$\psi (r)=\frac{1}{\sqrt{4\pi }}\left(\genfrac{}{}{0pt}{}{if(r)/r}{\sigma \widehat{r}g(r)/r}\right)\chi ^{(s)},$$
(24)
where $`\chi ^{(s)}`$ is a Pauli spinor, with the normalization $`\psi ^{}\psi d^3r=_0^{\mathrm{}}(f^2+g^2)𝑑r=1`$. Then the upper and lower components of the solution are

$`f(r)`$ $`=`$ $`N\left({\displaystyle \frac{r}{r_0}}\right)e^{r^2/2r_0^2},`$
$`g(r)`$ $`=`$ $`{\displaystyle \frac{N}{\sqrt{3}}}\left({\displaystyle \frac{r}{r_0}}\right)^2e^{r^2/2r_0^2},`$ (25)
$`N`$ $`=`$ $`\sqrt{8/(3r_0\sqrt{\pi })},r_0^2E_0^2=3,C={\displaystyle \frac{1}{9}}E_0^3,`$
where $`E_0`$ is the ground state eigenenergy and $`r_0`$ is a state dependent scale parameter.

Now we can calculate the matrix element of interest,

$`\overline{)}pm`$ $`=`$ $`{\displaystyle \overline{\psi }\frac{1}{2}(1+\gamma ^0)Cr^2\psi d^3r}`$
$`=`$ $`{\displaystyle _0^{\mathrm{}}}f(r)^2Cr^2𝑑r=\left({\displaystyle \frac{C}{3}}\right)^{1/3}={\displaystyle \frac{E_0}{3}}.`$
So, we can see that, in case of scalar+vector confinement of equal strengths, $`\overline{)}pm`$ is determined only by spin independent part of the Dirac spinor and is equal to one third of the ground state energy. In R. Tegen, R. Brockmann, and W. Weise (1982) it was also shown that for three massless quarks in their lowest 1$`s`$ orbit, with energy eigenvalues $`E_0`$ for each quark, the center-of-mass energy obtained with the potential (23) is just $`E_0`$, hence the nucleon mass in this model is: $`M_N=2E_0`$ (instead of $`M_N=3E_0`$ , as in non-relativistic and non-recoil models). Therefore, $`E_0=540`$ MeV and,

$$\overline{)}pm=180\mathrm{MeV}.$$
(27)
## V Pure Scalar Harmonic Potential

Tegen R. Tegen (1990) has considered scalar+vector harmonic confinement in calculating the weak neutron decay constant $`g_A/g_V`$ and found too small a value for $`g_A/g_V`$, compared to experiment. In R. Tegen, R. Brockmann, and W. Weise (1982) and R. Tegen, M. Schedl, and W. Weise (1983), a pure scalar harmonic potential $`V(r)=Cr^2`$ was studied numerically, and yielded more satisfactory results for $`g_A`$ and for the RMS charge radius. We find that

$$\overline{)}pm=C_0^{\mathrm{}}r^2\left(f(r)^2g(r)^2\right)𝑑r,$$
(28)
where $`f(r)`$ and $`g(r)`$ are defined as in eqn. (24).

We have fitted the numerical solution presented graphically in R. Tegen, R. Brockmann, and W. Weise (1982) with $`C=830\text{MeV/fm}^2`$ to calculate our integral of interest (28). The fitted wave functions are presented in Fig. 3, and as a benchmark for evaluation of the quality of the fit, we have calculated $`r^2`$ and $`g_A`$ and obtained values $`0.61`$ fm and $`1.26`$ respectively, as compared to $`r^2=0.64`$ fm and $`g_A=1.26`$ found in R. Tegen, R. Brockmann, and W. Weise (1982).

Thus we obtained, without any additional tuning, the following result:

$$\overline{)}pm=160\mathrm{M}\mathrm{e}\mathrm{V}.$$
(29)
## VI Summary

In this paper we have calculated, for the ground state of the quark in a nucleon, the matrix element of the operator $`(pm)`$, using variety of confinement potential models, under the assumption that the constituent quarks obey the Dirac equation. The motivation has been to solidify the estimates of the non-commutativity parameter of canonical (Lorentz violating) noncomutative QCD, where some of leading order Lorentz violating effects are proportional to factors of $`\overline{)}pm`$.

Interestingly, we found the following results,

$$\overline{)}pm=\{\begin{array}{cc}21\mathrm{MeV},\hfill & \text{for  3-D scalar central pot.}\hfill \\ 27\mathrm{MeV},\hfill & \text{for   scalar+vector linear pot.}\hfill \\ 180\mathrm{MeV},\hfill & \text{for   scalar+vector harm. pot.}\hfill \\ 160\mathrm{MeV},\hfill & \text{for   pure scalar harmonic pot.}\hfill \end{array}$$
(30)
We note that in the case of scalar central confinement as considered in section II, $`\overline{)}pm`$ vanishes as $`1/V_0`$ when $`V_0\mathrm{}`$, but it is different from zero in general. We note also that the value of $`\overline{)}pm`$ obtained for the scalar+vector linear confinement model is close to that obtained for scalar 3D potential well.

We have also shown that in case of scalar+vector harmonic confinement of equal strengths, $`\overline{)}pm`$ is determined only by spin independent part of the Dirac spinor and is equal to one third of the ground state energy.

For pure scalar harmonic confinement of the form $`V(r)=Cr^2`$, $`\overline{)}pm`$ was obtained using a fit to the numerical solution of the Dirac equation presented graphically in R. Tegen, R. Brockmann, and W. Weise (1982), and appears to have a value pretty close to that obtained for the scalar+vector harmonic confinement model.

Results obtained in this paper are within an order of magnitude agreement with the estimate made by Carlson et al. Carlson et al. (2001). The results obtained in Carlson et al. (2001) were used there to constrain the noncommutativity parameter in Lorentz violating noncommutative field theories. The constraints were very strong, and are still quite severe even if weakened by an order of magnitude. These results may be taken as a motivation to look for space-time noncomutativity in Lorentz-covariant ways Carlson et al. (2002b); Carlson and Carone (2002); Morita (2002); Berrino et al. (2002).

###### Acknowledgements.

The author expresses his gratitude to Prof. Carl Carlson for supervision of the current work. The author also expresses his appreciation to Prof. Christopher Carone for stimulating discussions and for careful reading of the manuscript.


text:
# Cosmological constant, dilaton field and Freund-Rubin compactification

## Acknowledgments

We would like to thank Kei-ichi Maeda and Daisuke Ida for fruitful discussions. We also give special thanks to James Overduin for reading this manuscript carefully. TS’s work is supported by Grant-in-Aid for Scientific Research from Ministry of Education, Science, Sports and Culture of Japan(No. 13135208, No.14740155 and No.14102004).