text:
# Detection methods for non-Gaussian gravitational wave stochastic backgrounds

## I Introduction and summary

Along with a new generation of gravitational wave detectors around the world ligo ; virgo ; geo ; tama , detection algorithms for a variety of sources are nearing completion. If the signals from these sources are detected, physicists stand to harvest unprecedented quantities of observational information concerning the nature of gravitation and the cosmos as a whole. The fruit of this harvest will be the outputs of detection algorithms. In this paper we introduce an algorithm designed for nearly optimal detection of a class of gravitational wave stochastic backgrounds. The non-Gaussian nature of this class of backgrounds causes the algorithm presented here to differ from the well studied cross-correlation based algorithms which are nearly optimal for Gaussian backgrounds.

### I.1 Gravitational wave stochastic backgrounds

Consider a large collection of similar gravitational wave sources. If we cannot resolve the individual signals produced by these sources and know only their statistical properties, the signals form a stochastic background. A wide variety of candidate sources of gravitational wave stochastic backgrounds have been studied (for an excellent general review see Ref. Allen Review ). These include high redshift supernovae gaussian supernovae ; non gaussian supernovae , the first stars or so-called population III objects first stars , rapidly rotating young neutron stars gaussian neutron stars 1 ; gaussian neutron stars 2 , early universe phase transitions and cosmic strings cosmic strings ; bubbles , inflation inflation , and high redshift compact binaries binaries .

Detecting a gravitational wave stochastic background produced by any one of these candidate sources could provide information on a variety of topics ranging from the evolution of the star formation rate Coward  to the numbers and sizes of posited extra dimensions Hogan . Because of this, stochastic backgrounds have long been thought to be among the most interesting possible types of gravitational radiation.

### I.2 Gaussian stochastic backgrounds

In order to develop detection methods, it is traditionally assumed that the individual events making up a background are uncorrelated and sufficiently frequent for the background to be Gaussian. That is, it is assumed that the conditions for applicability of the central limit theorem are satisfied.

Unlike electromagnetic waves, gravitational waves cannot be screened from a detector. Using a single gravitational wave detector, there is no practical way to distinguish between detector noise and a stochastic background of gravitational waves. As a consequence the sensitivity of a single detector to gravitational backgrounds is severely limited. By comparing the outputs of multiple detectors, sensitivity levels can be enhanced. Michelson Michelson  was the first to give a detailed description of such a detection method for a Gaussian stochastic background of gravitational waves in the presence of Gaussian detector noise. His detection strategy and its later refinements Christensen ; Flanagan ; Allen Romano  are often referred to as the cross-correlation method. Recently the cross-correlation method has been modified to treat more realistic detectors which themselves have sources of non-Gaussian noise robust gaussian ; robust gaussian II ; Klimenko and Mitselmakher .

We now briefly review the cross-correlation method. Consider two gravitational wave detectors. The output of each detector is a collection of dimensionless strain measurements. Suppose that $`N`$ such measurements are made with each detector at regular time intervals. Denote these measurements by a $`N\times 2`$ matrix $`h`$ with components $`h_i^k`$, where $`i=1,2`$ labels the detector, and $`k=1,2,\mathrm{},N`$ is a time index. To determine whether or not the data $`h`$ contains some desired signal, one usually compares the value of some detection statistic $`\mathrm{\Gamma }(h)`$ to a threshold value $`\mathrm{\Gamma }_{}`$. That is, if $`\mathrm{\Gamma }(h)>\mathrm{\Gamma }_{}`$ one concludes that a signal is present and otherwise one concludes that no signal is present. A detection statistic is said to be optimal if it yields the smallest probability of mistakenly concluding a signal is present (false alarm probability) after choosing a threshold which fixes the probability for mistakenly concluding a signal is absent (false dismissal probability).

Assume that the two detectors are collocated and aligned, and that each detector has white Gaussian noise with vanishing mean with no correlations between the two detectors. Then the standard cross-correlation detection statistic $`\mathrm{\Lambda }_{\text{CC}}`$ for a Gaussian signal is

$$\mathrm{\Lambda }_{\text{CC}}(h)=\frac{\widehat{\alpha }^2}{\overline{\sigma }_1\overline{\sigma }_2},$$
(1)
where

$`\widehat{\alpha }^2`$ $`=`$ $`\overline{\alpha }^2\theta (\overline{\alpha }^2),`$ (2)
$`\overline{\alpha }^2`$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{k=1}{\overset{N}{}}}h_1^kh_2^k,`$ (3)
$`\overline{\sigma }_i^2`$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{k=1}{\overset{N}{}}}\left(h_i^k\right)^2,`$ (4)
for $`i=1,2`$, and $`\theta (x)`$ is the Heaviside step function defined by

$$\theta (x)=\{\begin{array}{cc}1\hfill & \text{ if }x0\hfill \\ 0\hfill & \text{ if }x<0\hfill \end{array}.$$
(5)
This statistic is nearly optimal and can be derived from a maximum likelihood framework (see Sec. III.2). The subscript CC in $`\mathrm{\Lambda }_{\text{CC}}`$ denotes “cross correlation”. The generalization of this statistic to allow for colored noise and non-collocated, non-aligned detectors is discussed in Refs. Michelson ; Christensen ; Flanagan ; Allen Romano .

### I.3 Non-Gaussian stochastic backgrounds

A particular class of events will produce a Gaussian background if, on average, at any given moment, many individual events are arriving at the detector. However, if the ratio of average time between events to the average duration of events is large, then there are long stretches of “silence” or time during which no events arrive at the detector. The resulting stochastic background is non-Gaussian as the conditions for the applicability of the central limit theorem are not satisfied. Recent work has suggested that some candidate gravitational wave stochastic backgrounds, of both cosmological and astrophysical origin, may be non-Gaussian non gaussian supernovae ; cosmic strings ; first stars . However, predictions concerning the properties of most gravitational wave background sources rely heavily on theoretical arguments which extrapolate well beyond observational support. Such extrapolations are always in some sense speculative. It is conceivable that backgrounds predicted to be Gaussian may in fact turn out to be non-Gaussian, or vice versa.

In Sec. III.3 below, we apply a maximum likelihood framework to derive a detection statistic for a particular model of non-Gaussian stochastic background, which we now describe. Let $`h_i^k`$ be the outputs of two collocated aligned gravitational wave detectors with white, zero-mean, Gaussian noise with no correlations between the two detectors. The detector outputs $`h_i^k`$ consist of noise $`n_i^k`$ together with a common signal $`s^k`$:

$`h_1^k`$ $`=`$ $`n_1^k+s^k`$ (6)
$`h_2^k`$ $`=`$ $`n_2^k+s^k.`$
We wish to detect a non-Gaussian signal $`s^k`$ composed of long stretches of silence which separate short bursts whose amplitudes are Gaussianly distributed, and whose durations are smaller than the detector resolution time (see Fig. 4). We therefore assume that each signal sample $`s^k`$ is statistically independent with probability distribution \[cf. Eq. (49) below\]

$$p(s)=\xi \frac{1}{\sqrt{2\pi }\alpha }\mathrm{exp}\left[\frac{s^2}{2\alpha ^2}\right]+(1\xi )\delta (s).$$
(7)
The parameter $`\xi `$ is what we call the Gaussianity parameter of the stochastic background; it is the probability that, at any randomly chosen time, a burst is present in the detector. Thus $`\xi `$ takes values in the interval $`0\xi 1`$, and if $`\xi =1`$ then the background is Gaussian. The parameter $`\xi `$ can also be thought of as the duty cycle of the background. The parameter $`\alpha `$ in Eq. (7) is the rms amplitude of the bursts.

Our nearly-optimal detection statistic $`\mathrm{\Lambda }_{\text{ML}}^{\text{NG}}`$ for the signal model (7) is given by \[cf. Eq. (52) below\]

$`\mathrm{\Lambda }_{\text{ML}}^{\text{NG}}(h)`$ $`=`$ $`\underset{0<\xi 1}{\mathrm{max}}\underset{\alpha >0}{\mathrm{max}}\underset{\sigma _10}{\mathrm{max}}\underset{\sigma _20}{\mathrm{max}}{\displaystyle \underset{k=1}{\overset{N}{}}}\{{\displaystyle \frac{\overline{\sigma }_1\overline{\sigma }_2\xi }{\sqrt{\sigma _1^2\sigma _2^2+\sigma _1^2\alpha ^2+\sigma _2^2\alpha ^2}}}\mathrm{exp}[{\displaystyle \frac{\left(\frac{h_1^k}{\sigma _1^2}+\frac{h_2^k}{\sigma _2^2}\right)^2}{2\left(\frac{1}{\sigma _1^2}+\frac{1}{\sigma _2^2}+\frac{1}{\alpha ^2}\right)}}{\displaystyle \frac{\left(h_1^k\right)^2}{2\sigma _1^2}}{\displaystyle \frac{\left(h_2^k\right)^2}{2\sigma _2^2}}+1]`$ (8)
$`+`$ $`{\displaystyle \frac{\overline{\sigma }_1\overline{\sigma }_2}{\sigma _1\sigma _2}}(1\xi )\mathrm{exp}[{\displaystyle \frac{\left(h_1^k\right)^2}{2\sigma _1^2}}{\displaystyle \frac{\left(h_2^k\right)^2}{2\sigma _2^2}}+1]\}.`$
Here the quantities $`\overline{\sigma }_1`$ and $`\overline{\sigma }_2`$ are defined by Eq. (4). The values of $`\xi `$, $`\alpha ^2`$, $`\sigma _1^2`$ and $`\sigma _2^2`$ which achieve the maximum in Eq. (8) are, respectively, estimators of the signal’s Gaussianity parameter, the variance of the signal events, and the variances of the noise in the two detectors. If we calculate the quantity (8) at $`\xi =1`$, instead of maximizing over $`\xi `$, the result is a statistic which is equivalent to the standard cross-correlation statistic $`\mathrm{\Lambda }_{\text{CC}}`$.

The subscript ML on $`\mathrm{\Lambda }_{\text{ML}}^{\text{NG}}`$ stands for “maximum likelihood”, while the superscript NG stands for “non-Gaussian statistic”. The superscript NG does not necessarily mean that one is considering a non-Gaussian signal; both of the statistics $`\mathrm{\Lambda }_{\mathrm{CC}}`$ and $`\mathrm{\Lambda }_{\text{ML}}^{\text{NG}}`$ can be applied to data containing either a Gaussian signal or a non-Gaussian signal.

If the burst-amplitude parameter $`\alpha `$ is sufficiently large and the bursts are well separated in time, then the individual bursts can be seen in the detector output. In this case one could use, for example, the simple burst statistic 333In reality the statistic (9) would be especially susceptible to non-Gaussian noise bursts in the detector and so would not be used in practice; instead one would need search for events where $`|h_1^k|`$ and $`|h_2^k|`$ are simultaneously large. In this paper we restrict attention for simplicity to Gaussian detector noise; it will be important for future more general analyses to to allow for (uncorrelated) non-Gaussian noise components in the two detectors.

$$\mathrm{\Lambda }_\text{B}\underset{1kN}{\mathrm{max}}\left|h_1^k\right|.$$
(9)
on the data from detector 1 to detect the signal. The burst statistic (9) and the cross-correlation statistic $`\mathrm{\Lambda }_{\text{CC}}`$ are used as references for comparison for the maximum likelihood statistic below.

### I.4 Main results

There are two main results in this paper. The first result is the detection statistic $`\mathrm{\Lambda }_{\text{ML}}^{\text{NG}}`$ given by Eq. (8), which is derived in Sec. III.3. This statistic is nearly optimal for the detection of a class of non-Gaussian gravitational wave stochastic backgrounds incident on a pair of idealized detectors.

The second main result, summarized in Figs. 1 and 2, is a comparison of the performances of the maximum likelihood statistic $`\mathrm{\Lambda }_{\text{ML}}^{\text{NG}}`$, the cross-correlation statistic $`\mathrm{\Lambda }_{\text{CC}}`$, and the burst statistic $`\mathrm{\Lambda }_\text{B}`$.

That comparison is quantified in terms of the the minimum gravitational-wave energy density $`\mathrm{\Omega }_{\text{detectable}}`$ necessary for detection. The values of this quantity for the three different statistics $`\mathrm{\Lambda }_{\text{ML}}^{\text{NG}}`$, $`\mathrm{\Lambda }_{\text{CC}}`$ and $`\mathrm{\Lambda }_\text{B}`$ we will denote by $`\mathrm{\Omega }_{\text{detectable}}^{\text{ML}}`$, $`\mathrm{\Omega }_{\mathrm{detectable}}^{\text{CC}}`$, and $`\mathrm{\Omega }_{\text{detectable}}^\text{B}`$, respectively. Results for these three quantities obtained from Monte Carlo simulations are shown in Fig. 1, which gives $`\mathrm{\Omega }_{\text{detectable}}`$ as a function of $`\xi `$ for $`N=10^4`$ data points. The Monte Carlo simulations are described in Sec. IV.2 below. The figure shows that in the limit $`\xi 1`$ of Gaussian signals, the statistics $`\mathrm{\Lambda }_{\text{ML}}^{\text{NG}}`$ and $`\mathrm{\Lambda }_{\text{CC}}`$ perform approximately equivalently (the cross-correlation statistic is slightly better). As the Gaussianity parameter $`\xi `$ is decreased, the performance of $`\mathrm{\Lambda }_{\text{ML}}^{\text{NG}}`$ improves, until at $`\xi 10^{2.5}`$ it is better than that of $`\mathrm{\Lambda }_{\text{CC}}`$ by about a factor of $`3`$ in energy density. Finally, in the limit $`\xi 0`$, the individual bursts become visible and the burst statistic $`\mathrm{\Lambda }_\text{B}`$ becomes the best statistic.

Figure 1 also shows theoretical curves for the three quantities $`\mathrm{\Omega }_{\text{detectable}}^{\text{ML}}`$, $`\mathrm{\Omega }_{\mathrm{detectable}}^{\text{CC}}`$, and $`\mathrm{\Omega }_{\text{detectable}}^\text{B}`$. These curves are derived and discussed in Sec. IV below. For the burst and cross-correlation statistics, the theoretical curves should have a fractional accuracy $`1/\sqrt{N}`$. For the maximum likelihood statistic, the theoretical prediction is expected to be accurate to a few tens of percent. These expected accuracies are confirmed by the Monte Carlo simulations, as seen in Fig. 1.

The value $`N=10^4`$ of the number of data points is roughly appropriate for a space based detector like LISA, for which the duration of a measurement might be $`1`$ year and the effective bandwidth $`10^3`$ Hz. However, for year-long observations with ground based detectors, the effective bandwidth will be $`100`$ Hz and consequently the appropriate value of $`N`$ is $`10^9`$. We were unable to perform Monte Carlo simulations for this large value of $`N`$ due to limitations in available computing power. However, we show in Fig. 2 the theoretical curves for the three different statistics as functions of $`\xi `$ for $`N=10^9`$. In this case, the maximum likelihood statistic starts to outperform the cross-correlation statistic at $`\xi 10^3`$, and the maximum gain factor in energy density is of order $`2`$.

We next discuss the computational cost of the maximum likelihood statistic $`\mathrm{\Lambda }_{\mathrm{ML}}^{\mathrm{NG}}`$. As is well known, the computational cost of trying to detect a stochastic background using the cross-correlation statistic $`\mathrm{\Lambda }_{\text{CC}}`$ is negligible when compared to, say, matched-filter-based inspiral waveform searches. However, because of the non-trivial maximization in Eq. (8), the maximum likelihood statistic $`\mathrm{\Lambda }_{\text{ML}}^{\text{NG}}`$ is computationally intensive. In fact, every evaluation of the function to be maximized over the four parameters $`\xi `$, $`\alpha `$, $`\sigma _1`$, and $`\sigma _2`$ requires computing a length-$`N`$ sum or product, where $`N`$ is the number of data points, and takes longer than the entire cross-correlation detection method. Depending on the method of calculation, the computational cost of computing $`\mathrm{\Lambda }_{\text{ML}}^{\text{NG}}`$ is larger than that of computing $`\mathrm{\Lambda }_{\text{CC}}`$ by a factor anywhere from $`10^2`$ to $`10^4`$.

To summarize, under the idealized assumptions of this paper, if one searches for a stochastic background using the standard cross-correlation statistic, then one might not detect a signal that would have been detectable using our maximum likelihood statistic. This conclusion probably generalizes to realistic detector noise models and detector orientations.

### I.5 Outline of this paper

In Sec. II we introduce notation, review the general theory of signal detection and parameter measurement, and derive a general form of the maximum likelihood detection statistic. Then, in Sec. III, we derive the maximum likelihood statistics for both a Gaussian background (Sec. III.2) and for the model (7) of a non-Gaussian background (Sec. III.3), assuming two idealized detectors. In Sec. IV we discuss analytical calculations and Monte Carlo simulations comparing the performance of the maximum likelihood and cross-correlation detection statistics. Also in Sec. IV we show how the signal parameters $`\xi `$ and $`\alpha `$ can be estimated, with reasonable accuracy, for a strong non-Gaussian background. We conclude in Sec. V with a discussion of the results.

## II General theory of detection statistics and parameter estimation

In this section we review various formal aspects of the theory of signal detection and measurement. We derive a form of the maximum likelihood detection statistic that is more general than has been considered before in the context of gravitational wave data analysis Allen Romano ; general method ; excess power ; sam joe ; sam unpublished . The material in this section can be found in a variety of texts maximum likelihood ; we include this section for completeness and to introduce notation.

### II.1 Notational conventions

We use calligraphic letters $`𝒜,,𝒞,\mathrm{}`$ to denote random variables. As described in Sec. I.2, given $`D`$ detectors we can assemble an $`N\times D`$ detector output matrix $``$ with components $`_i^k`$ where $`k=1,2,\mathrm{},N`$ is a time index, and $`i=1,2,\mathrm{},D`$ labels the detector. We assume that the detector outputs are made up of noise $`𝒩`$ and signal $`𝒮`$ with components $`𝒩_i^k`$ and $`𝒮_i^k`$ respectively, such that

$$=𝒩+𝒮.$$
(10)
Specific realizations of random variables will be denoted by lower case Roman symbols. For example, $`h=n+s`$ is a specific realization of Eq. (10), where the components of $`h`$ are $`h_i^k`$.

Probability densities for random variables will always be denoted by a lowercase $`p`$ and will carry a subscript to indicate which random variable is being described. For example, $`p_𝒩(n)d^{ND}n`$ is the probability that $`n<𝒩<n+dn`$, where $`d^{ND}n`$ is the product

$$d^{ND}n=\underset{k=1}{\overset{N}{}}\underset{i=1}{\overset{D}{}}dn_i^k.$$
(11)
We write the normalization requirement for $`p_𝒩(n)`$ as

$$1=d^{ND}np_𝒩(n).$$
(12)
Unless otherwise specified, integrals are over $`^{ND}`$ where $``$ is the set of real numbers.

We assume a detector noise model with $`Q_n`$ parameters. Let $`𝒱_n`$ be a vector of length $`Q_n`$ whose components are the parameters characterizing the noise in the detectors. We denote by $`\mathrm{\Theta }_n`$ the space of all possible values of $`𝒱_n`$. Here the subscript $`n`$ is not an index; it is merely short for “noise”. We denote joint probabilities in the usual way. For example, $`p_{𝒩,𝒱_n}(n,𝐯_n)d^{ND}nd^{Q_n}v_n`$ is the probability that $`n<𝒩<n+dn`$ and $`𝐯_n<𝒱_n<𝐯_n+d𝐯_n`$, where $`d^{Q_n}v_n`$ is defined by

$$d^{Q_n}v_n=\underset{l=1}{\overset{Q_n}{}}dv_n^l,$$
(13)
and $`dv_n^l`$ is the $`l`$th component of $`d𝐯_n`$. We also use vertical bars to denote conditional probabilities. For example

$$p_{𝒩|𝒱_n}(n|𝐯_n)d^{ND}n=\frac{p_{𝒩,𝒱_n}(n,𝐯_n)d^{Q_n}v_n}{p_{𝒱_n}(𝐯_n)d^{Q_n}v_n}d^{ND}n$$
(14)
is the probability that $`n<𝒩<n+dn`$ given that $`𝐯_n<𝒱_n<𝐯_n+d𝐯_n`$.

We will often use the so-called total probability theorem Papoulis  to write probability densities for a specific random variable as an integral over the functional dependencies of that random variable. An example is

$$p_𝒩(n)=_{\mathrm{\Theta }_n}d^{Q_n}v_np_{𝒩|𝒱_n}(n|𝐯_n)p_{𝒱_n}(𝐯_n).$$
(15)
Expanding probability densities in this way allows us to treat parameters, such as the noise parameters $`𝒱_n`$ in Eq. (15), as unknowns. In fact, such a treatment of the noise parameters is the crucial difference between the derivations of this work and those in previous studies of gravitational wave data analysis techniques Allen Romano ; general method ; excess power ; sam joe ; sam unpublished .

We assume that the signal model contains $`Q_s`$ parameters, which we will treat as random variables like the noise parameters. We will denote by $`𝒱_s`$ the random vector of length $`Q_s`$ containing the signal parameters, and by $`\mathrm{\Theta }_s`$ the space of all possible values of $`𝒱_s`$.

We define the notions of “signal present” and “signal absent” in terms of a partition of the space $`\mathrm{\Theta }_s`$ of signal parameters into a disjoint union

$$\mathrm{\Theta }_s=\mathrm{\Theta }_{s0}\mathrm{\Theta }_{s1},$$
(16)
where $`\mathrm{\Theta }_{s0}`$ corresponds to the signal being absent, and $`\mathrm{\Theta }_{s1}`$ the signal being present. We define the random variable $`𝒯`$, taking values $`𝒯=0`$ or $`𝒯=1`$, according to

$$𝒯=\{\begin{array}{cc}1\hfill & \text{ if }𝒱_s\mathrm{\Theta }_{s1}\hfill \\ 0\hfill & \text{ if }𝒱_s\mathrm{\Theta }_{s0}\hfill \end{array}.$$
(17)
Thus $`𝒯=1`$ corresponds to a signal being present, and $`𝒯=0`$ to no signal being present. We define

$$p_{𝒮|𝒱_s,𝒯}(s|𝐯_s,0)=\{\begin{array}{cc}0\hfill & \text{ if }𝐯_s\mathrm{\Theta }_{s1}\hfill \\ \delta ^{ND}(s)\hfill & \text{ if }𝐯_s\mathrm{\Theta }_{s0}\hfill \end{array},$$
(18)
where $`\delta ^{ND}(s)`$ is the $`N\times D`$ dimensional Dirac delta function. We denote by $`p_{𝒯,}(t,h)d^{ND}h`$ the probability that $`𝒯=t`$ and that $`h<<h+dh`$, where $`t=0`$ or $`1`$. Similarly

$$p_{|𝒯}(h|t)d^{ND}h=\frac{p_{,𝒯}(h,t)}{P_𝒯(t)}d^{ND}h$$
(19)
is the probability that $`h<<h+dh`$ given that $`𝒯=t`$.

We denote probabilities (as opposed to probability densities) with an uppercase $`P`$. For example $`P_𝒯(1)`$ is the probability that a signal is present, and $`P_𝒯(0)`$ is the probability that a signal is absent.

Before examining the detector outputs, we may have some idea, say from previous experiments, of the probability that a signal will be present. We denote this prior probability by $`P^{(0)}`$. We denote by $`P^{(1)}`$ the posterior probability that the signal is present after examining $``$ in the context of all prior experiments etc. All posterior quantities have an implicit dependence on the detector outputs. To simplify the notation we will not explicitly show this dependence. For example, we write $`P^{(1)}`$ rather than the more cumbersome $`P^{(1)}()`$ for the posterior probability that a signal is present.

There are prior and posterior versions of all probability densities. When necessary we will append superscripts of $`(0)`$ and $`(1)`$ to distinguish priors and posteriors respectively. For example $`p_{𝒱_n}^{(1)}(𝐯_n)=p_{𝒱_n|}(𝐯_n|h)`$ is the posterior probability density for $`𝒱_n`$. The posterior distribution for the noise can be expanded in terms of $`p_{𝒱_n}^{(1)}(𝐯_n)`$ as

$$p_𝒩^{(1)}(n)=_{\mathrm{\Theta }_n}d^{Q_n}v_np_{𝒩|𝒱_n}(n|𝐯_n)p_{𝒱_n}^{(1)}(𝐯_n).$$
(20)
The conventions and symbols which have been introduced above are summarized in tables 1 and 2 respectively.

### II.2 Detection statistics

To detect a signal one uses a detection statistic, say $`\mathrm{\Gamma }=\mathrm{\Gamma }()`$, that is some function of the detector outputs $``$. A signal is said to have been detected when $`\mathrm{\Gamma }`$ exceeds some threshold value $`\mathrm{\Gamma }_{}`$.

Denote by $`P_{\text{FD}}(\mathrm{\Gamma }_{})`$ the probability of false dismissal, that is, the probability that we fail to detect a signal which is actually present. Similarly, let $`P_{\text{FA}}(\mathrm{\Gamma }_{})`$ be the probability that we claim to have detected a signal which in fact is absent—the probability of false alarm. For given signal and noise models and for a given statistic $`\mathrm{\Gamma }`$, the false alarm and false dismissal probabilities generate a curve in the $`P_{\text{FA}}`$-$`P_{\text{FD}}`$ plane parametrized by the threshold $`\mathrm{\Gamma }_{}`$. Such curves depend on the number of detectors $`D`$, the number of data points $`N`$, the signal parameters $`𝒱_s`$, and the noise parameters $`𝒱_n`$.

Suppose that the statistic $`\mathrm{\Gamma }`$ is bounded in the sense that there exist numbers $`\mathrm{\Gamma }_{\mathrm{min}}`$ and $`\mathrm{\Gamma }_{\mathrm{max}}`$ such that $`\mathrm{\Gamma }_{\mathrm{min}}<\mathrm{\Gamma }<\mathrm{\Gamma }_{\mathrm{max}}`$ for all $``$. Then it is clear that $`P_{\text{FD}}(\mathrm{\Gamma }_{\mathrm{min}})=0`$ and that $`P_{\text{FA}}(\mathrm{\Gamma }_{\mathrm{min}})=1`$. As the threshold $`\mathrm{\Gamma }_{}`$ increases toward $`\mathrm{\Gamma }_{\mathrm{max}}`$, $`P_{\text{FD}}(\mathrm{\Gamma }_{})`$ will increase while $`P_{\text{FA}}(\mathrm{\Gamma }_{})`$ decreases, until finally at $`\mathrm{\Gamma }_{}=\mathrm{\Gamma }_{\mathrm{max}}`$, $`P_{\text{FD}}=1`$, and $`P_{\text{FA}}=0`$. Thus, false dismissal-false alarm curves generally look something like those sketched in Fig. 3.

Note that if one uses a different statistic $`f(\mathrm{\Gamma })`$, where $`f`$ is any function, then the shape of the $`P_{\text{FA}}`$-$`P_{\text{FD}}`$ curve does not change as long as $`f`$ is monotonic in the sense that

$$\mathrm{\Gamma }>\mathrm{\Gamma }_{}f(\mathrm{\Gamma })>f(\mathrm{\Gamma }_{}).$$
(21)
Only the parametrization of the curve changes under such a transformation. Statistics related by transformations $`f`$ satisfying the monotonicity property (21) have identical false dismissal versus false alarm curves.

In 1933 Neyman and Pearson considered a simple signal detection scenario where the sets $`\mathrm{\Theta }_n`$, $`\mathrm{\Theta }_{s1}`$, and $`\mathrm{\Theta }_{s0}`$ each contain a single element Neyman and Pearson . They showed that for this scenario the detection statistic which minimizes $`P_{\text{FD}}`$ for any $`P_{\text{FA}}`$ is the so-called *likelihood ratio* $`\mathrm{\Lambda }`$, defined by

$$\mathrm{\Lambda }=\frac{p_{|𝒯}(h|1)}{p_{|𝒯}(h|0)}.$$
(22)
One notion of optimality for detection statistics is that the statistic should minimize the false dismissal probability at a fixed value of the false alarm probability. For the simple scenario above, this criteria, known as the Neyman-Pearson criteria, uniquely determines the likelihood ratio as the optimal statistic Ferguson . However in general, when any of $`\mathrm{\Theta }_n`$, $`\mathrm{\Theta }_{s1}`$, or $`\mathrm{\Theta }_{s0}`$ contains more than one element, the statistic selected by this criteria is a function of the unknown parameters $`𝒱_s`$ and $`𝒱_n`$. Thus, as is well known, the Neyman-Pearson criteria does not single out a unique statistic in such cases.

In this paper we will obtain our detection statistics from Bayesian considerations, but we will quantify their effectiveness using the Neyman and Pearson criteria of comparing false dismissal probabilities at fixed false alarm probabilities.

### II.3 Likelihood ratio and likelihood function

From a Bayesian point of view, a natural criterion for deciding that a signal is present is for the posterior probability $`P^{(1)}`$ to exceed some threshold Bayes . The posterior probability $`P^{(1)}`$ is related to the prior probability $`P^{(0)}`$ and to the likelihood ratio $`\mathrm{\Lambda }`$ defined by Eq. (22) by

$$\frac{P^{(1)}}{1P^{(1)}}=\mathrm{\Lambda }\frac{P^{(0)}}{1P^{(0)}}.$$
(23)
See appendix A for a derivation of Eq. (23) in the most general context where the sets $`\mathrm{\Theta }_n`$, $`\mathrm{\Theta }_{s1}`$, and $`\mathrm{\Theta }_{s0}`$ are all non-trivial. It follows from Eq. (23) that $`P^{(1)}`$ is a monotonic function of $`\mathrm{\Lambda }`$, so thresholding on $`P^{(1)}`$ is equivalent to thresholding on $`\mathrm{\Lambda }`$. This makes $`\mathrm{\Lambda }`$, or approximate versions of it, the natural choice for a detection statistic.

We derive in Appendix A the following general formula for the likelihood ratio as a function of the data $`=h`$:

$$\mathrm{\Lambda }=\frac{{\displaystyle _{\mathrm{\Theta }_{s1}}}d^{Q_s}v_s{\displaystyle d^{ND}s_{\mathrm{\Theta }_n}d^{Q_n}v_np_{𝒩|𝒱_n}(hs|𝐯_n)p_{𝒱_n}(𝐯_n)p_{𝒮|𝒱_s,𝒯}(s|𝐯_s,1)p_{𝒱_s|𝒯}(𝐯_s|1)}}{{\displaystyle _{\mathrm{\Theta }_n}}d^{Q_n}v_n^{}p_{𝒩|𝒱_n}(h|𝐯_n^{})p_{𝒱_n}(𝐯_n^{})}.$$
(24)
The various probability distributions that appear in Eq. (24) are (i) the prior distribution $`p_{𝒱_s|𝒯}(𝐯_s|1)`$ for the signal parameters $`𝐯_s`$; (ii) the distribution $`p_{𝒮|𝒱_s,𝒯}(s|𝐯_s,1)`$ for the signal $`s`$ given the signal parameters $`𝐯_s`$; (iii) the prior distribution $`p_{𝒱_n}(𝐯_n)`$ for the noise parameters $`𝐯_n`$; and (iv) the distribution $`p_{𝒩|𝒱_n}(h|𝐯_n)`$ for the noise $`n`$ given the noise parameters $`𝐯_n`$.

We can interpret Eq. (24) as follows. In the simple signal detection scenario, we choose between a pair of simple claims: (i) $`𝒱_s=𝐯_{s0}`$ or (ii) $`𝒱_s=𝐯_{s1}`$. In general we choose between a pair of complicated, or composite, claims: (i) $`𝒱_s\mathrm{\Theta }_{s0}`$ or (ii) $`𝒱_s\mathrm{\Theta }_{s1}`$, where both $`\mathrm{\Theta }_{s0}`$ and $`\mathrm{\Theta }_{s1}`$ contain many elements. Equation (24) says that the best way to chose between a pair of complicated claims is to first break the complicated pair of claims into pairs of simple claims, then compute the likelihood ratio for each pair of simple claims, and sum the results of each choice. That is, the likelihood ratio can be written as an integral over the parameters of the composite claims

$$\mathrm{\Lambda }=_{\mathrm{\Theta }_{s1}}d^{Q_s}v_s_{\mathrm{\Theta }_n}d^{Q_n}v_n\mathrm{\Lambda }(𝐯_s,𝐯_n),$$
(25)
where the integrand $`\mathrm{\Lambda }(𝐯_s,𝐯_n)`$, which we refer to as the *likelihood function*, can be read off from Eq. (24):

$$\mathrm{\Lambda }(𝐯_s,𝐯_n)=\frac{{\displaystyle d^{ND}sp_{𝒩|𝒱_n}(hs|𝐯_n)p_{𝒮|𝒱_s,𝒯}(s|𝐯_s,1)p_{𝒱_n}(𝐯_n)p_{𝒱_s|𝒯}(𝐯_s,1)}}{{\displaystyle _{\mathrm{\Theta }_n}}d^{Q_n}v_n^{}p_{𝒩|𝒱_n}(h|𝐯_n^{})p_{𝒱_n}(𝐯_n^{})}.$$
(26)
The likelihood function444There are two different conventions for the definition of the likelihood function. Some authors include the probability distributions for $`𝒱_s`$ and $`𝒱_n`$ in the definition of $`\mathrm{\Lambda }(𝐯_s,𝐯_n)`$ as we have in Eq. (26), while others leave these out of $`\mathrm{\Lambda }(𝐯_s,𝐯_n)`$ and would show these distributions explicitly in Eq. (25). $`\mathrm{\Lambda }(𝐯_s,𝐯_n)`$ can be used to compute the posterior probability density $`p_{𝒱_s,𝒱_n|𝒯}^{(1)}(𝐯_s,𝐯_n|1)`$ for the signal and noise parameters given that a signal is present, via the formula

$$\frac{P^{(1)}}{1P^{(1)}}p_{𝒱_s,𝒱_n|𝒯}^{(1)}(𝐯_s,𝐯_n|1)=\mathrm{\Lambda }(𝐯_s,𝐯_n)\frac{P^{(0)}}{1P^{(0)}}.$$
(27)
A derivation of Eq. (27) can be found in appendix A.

### II.4 Maximum likelihood detection statistics and parameter estimators

In many applications, it is impractical to compute the detection statistic (24) because of the multi-dimensional integrals involved Loredo . However, approximate versions of the statistic are often easier to compute and useful. If a signal is present with sufficiently large amplitude, then the integrand in the numerator of Eq. (24) will be sharply peaked. The integrand in the denominator of Eq. (24) will also be sharply peaked when there is sufficient data that the noise is well characterized. Under these circumstances, the integrals can be written as the values of the corresponding integrands at the peaks multiplied by “width factors”, where the width factors depend only weakly on the data $`h`$ and can be neglected without affecting much the performance of the statistic. \[The width factors from the integrals over the noise parameters will tend to cancel between the numerator and denominator\]. Also, frequently the prior distributions for $`𝒱_s`$ and $`𝒱_n`$ are slowly varying, and neglecting those distributions has a negligible effect on the performance of the statistic. Under these conditions the maximum likelihood detection statistic $`\mathrm{\Lambda }_{\text{ML}}`$ defined by

$$\mathrm{\Lambda }_{\text{ML}}=\frac{\underset{𝐯_s\mathrm{\Theta }_{s1}}{\mathrm{max}}\underset{𝐯_n\mathrm{\Theta }_n}{\mathrm{max}}{\displaystyle d^{ND}sp_{𝒩|𝒱_n}(hs|𝐯_n)p_{𝒮|𝒱_s,𝒯}(s|𝐯_s,1)}}{\underset{𝐯^{}\mathrm{\Theta }_n}{\mathrm{max}}p_{𝒩|𝒱_n}(h|𝐯_n^{})}$$
(28)
is a natural approximate version of $`\mathrm{\Lambda }`$ 555In the event that the priors for $`𝒱_s`$ and $`𝒱_n`$ restrict these parameters to regions $`\mathrm{\Theta }_{s1}^{}\mathrm{\Theta }_{s1}`$ and $`\mathrm{\Theta }_n^{}\mathrm{\Theta }_n`$, the bounds of the maximizations in Eq. (28) should be changed to $`\mathrm{\Theta }_{s1}\mathrm{\Theta }_{s1}^{}`$ and $`\mathrm{\Theta }_n\mathrm{\Theta }_n^{}`$.. The subscript ML denotes that (28) is the maximum likelihood approximate version of $`\mathrm{\Lambda }`$. See Ref. maximum likelihood  for further discussion of $`\mathrm{\Lambda }_{\text{ML}}`$ as an approximate version of $`\mathrm{\Lambda }`$ 666Note that $`\mathrm{\Lambda }_{\mathrm{ML}}`$ is an approximate version of $`\mathrm{\Lambda }`$ only in the sense that the false dismissal versus false alarm curves of the two statistics will be close to one another. The numerical values of $`\mathrm{\Lambda }_{\mathrm{ML}}`$ and $`\mathrm{\Lambda }`$ will in general differ significantly, due to the width factors and priors. Therefore the statistic $`\mathrm{\Lambda }_{\mathrm{ML}}`$ cannot be used in Eq. (23) to compute Bayesian thresholds for detection given a desired value of $`P^{(1)}`$..

A particular special case of the detection statistic (28), which is widely used, is the following. Assume that the noise parameters have some known values $`𝒱_n=𝐯_n`$. Then the noise priors and the $`\mathrm{\Theta }_n`$ integrals in Eq. (24) are trivial, and one obtains the detection statistic

$$\stackrel{~}{\mathrm{\Lambda }}_{\text{ML}}=\frac{\underset{𝐯_s\mathrm{\Theta }_{s1}}{\mathrm{max}}{\displaystyle d^{ND}sp_{𝒩|𝒱_n}(hs|𝐯_n)p_{𝒮|𝒱_s,𝒯}(s|𝐯_s,1)}}{p_{𝒩|𝒱_n}(h|𝐯_n)}.$$
(29)
See Ref. sam joe  for an exploration of the statistic (29) in the context of stochastic backgrounds. We will show below that for a Gaussian stochastic background, $`\mathrm{\Lambda }_{\text{ML}}`$ reduces to the standard cross-correlation statistic while the more specialized statistic $`\stackrel{~}{\mathrm{\Lambda }}_{\text{ML}}`$ does not. Thus for stochastic backgrounds, treating the noise parameters as unknowns is crucial robust gaussian II .

When the noise and signal parameters $`𝒱_n`$ and $`𝒱_s`$ can take on many values, one naturally would like to know which values are realized. Equation (27) suggests using the values $`\widehat{𝐯}_n`$ and $`\widehat{𝐯}_s`$ defined by

$$\mathrm{\Lambda }(\widehat{𝐯}_s,\widehat{𝐯}_n)=\underset{𝐯_s\mathrm{\Theta }_{s1}}{\mathrm{max}}\underset{𝐯_n\mathrm{\Theta }_n}{\mathrm{max}}\mathrm{\Lambda }(𝐯_s,𝐯_n).$$
(30)
The estimators $`\widehat{𝐯}_n`$ and $`\widehat{𝐯}_s`$ are known as maximum likelihood estimators. Note that $`𝐯_s=\widehat{𝐯}_s`$ and $`𝐯_n=\widehat{𝐯}_n`$ also maximize the numerator in Eq. (28). For the remainder of this paper we will use $`\mathrm{\Lambda }_{\text{ML}}`$, defined by Eq. (28), as our detection statistic, and $`\widehat{𝐯}_s`$ and $`\widehat{𝐯}_n`$, defined by Eq. (30), as parameter estimators.

## III Application to stochastic background searches

In this section we derive the maximum likelihood detection statistic (28) for a simplified model of the detection problem for stochastic gravitational waves, and for a specific simple model of a non-Gaussian stochastic background.

### III.1 Assumptions

We assume two detectors with outputs $`_i^k`$, where $`i=1,2`$ labels the detector and $`k=1,2,\mathrm{},N`$ is a time index. We assume that the noise in detector one is uncorrelated with the noise in detector two. We will require the noise in both detectors to have vanishing mean and to be both Gaussian and white, so that

$$p_{𝒩|𝒱_n}\left[n|(\sigma _1,\sigma _2)\right]=\underset{k=1}{\overset{N}{}}\frac{1}{2\pi \sigma _1\sigma _2}\mathrm{exp}\left[\frac{(n_1^k)^2}{2\sigma _1^2}\frac{(n_2^k)^2}{2\sigma _2^2}\right].$$
(31)
The parameters $`\sigma _1`$ and $`\sigma _2`$ in Eq. (31) are the square roots of the variances of the noise in the two detectors. For this model $`𝐯_n=(\sigma _1,\sigma _2)`$ and $`\mathrm{\Theta }_n=\left\{(\sigma _1,\sigma _2)\right|\sigma _10\text{ and }\sigma _20\}`$.

We assume that the detectors are collocated and aligned, so that the same signal is present in both detectors

$$𝒮_1^k=𝒮_2^k=𝒮^k.$$
(32)
Lastly we assume that the individual signal samples are uncorrelated and identically distributed, i.e., the signal is white, so that

$$p_𝒮(s)=\underset{k=1}{\overset{N}{}}p_{𝒮^k}(s^k).$$
(33)
Our assumptions (31)-(33) are unrealistic for both ground-based and space-based detectors: we expect the noise to be colored with significant non-Guasssian components, and in general detectors will not be co-located and aligned. Our analysis is therefore just a first step, and will need to be generalized. However, we expect that our central conclusion — the existence of statistics which outperform the standard cross-correlation statistic for nonGaussian signals — is robust, and will not be altered when these complications are taken into account.

We now derive a general formula for the maximum likelihood statistic (28), which we apply in both the Gaussian and non-Gaussian cases in the following two subsections. The denominator in Eq. (28) can be written, from Eq. (31), as

$$\underset{\sigma _10}{\mathrm{max}}\underset{\sigma _20}{\mathrm{max}}\left\{\left(2\pi \sigma _1\sigma _2\right)^N\mathrm{exp}\left[\frac{N}{2}\left(\frac{\overline{\sigma }_1^2}{\sigma _1^2}+\frac{\overline{\sigma }_2^2}{\sigma _2^2}\right)\right]\right\},$$
(34)
where $`\overline{\sigma }_1^2`$ and $`\overline{\sigma }_2^2`$ are defined by

$$\overline{\sigma }_i^2=\frac{1}{N}\underset{k=1}{\overset{N}{}}\left(h_i^k\right)^2$$
(35)
for $`i=1,2`$. It is easily shown that the maximum in Eq. (34) is achieved at $`\sigma _i=\overline{\sigma }_i`$. From Eq. (28) this yields

$$\mathrm{\Lambda }_{\text{ML}}=\frac{\underset{𝐯_s\mathrm{\Theta }_{s1}}{\mathrm{max}}\underset{𝐯_n\mathrm{\Theta }_n}{\mathrm{max}}{\displaystyle d^{ND}sp_{𝒩|𝒱_n}(hs|𝐯_n)p_{𝒮|𝒱_s,𝒯}(s|𝐯_s,1)}}{\left(2\pi \overline{\sigma }_1\overline{\sigma }_2\right)^N\mathrm{exp}\left(N\right)}.$$
(36)
Combining this with Eq. (33) yields the following final general expression for the maximum likelihood statistic:

$$\mathrm{\Lambda }_{\text{ML}}=\underset{𝐯_s\mathrm{\Theta }_{s1}}{\mathrm{max}}\underset{\sigma _10}{\mathrm{max}}\underset{\sigma _20}{\mathrm{max}}\underset{k=1}{\overset{N}{}}\frac{\overline{\sigma }_1\overline{\sigma }_2}{\sigma _1\sigma _2}_{\mathrm{}}^{\mathrm{}}𝑑s^kp_{𝒮^k|𝒱_s,𝒯}(s^k|𝐯_s,1)\mathrm{exp}\left[\frac{\left(h_1^ks^k\right)^2}{2\sigma _1^2}\frac{\left(h_2^ks^k\right)^2}{2\sigma _2^2}+1\right].$$
(37)
### III.2 Gaussian signal

We now consider the case where the signal is Gaussian and has a vanishing mean. We denote by $`\alpha ^2`$ the variance of the signal, so the prior for $`𝒮`$ is given by

$$p_{𝒮^k|𝒱_s,𝒯}(s^k|\alpha ,1)=\frac{1}{\sqrt{2\pi }\alpha }\mathrm{exp}\left[\frac{\left(s^k\right)^2}{2\alpha ^2}\right].$$
(38)
For this model $`𝐯_s=(\alpha )`$ has only one component, and $`\mathrm{\Theta }_{s1}=\{\alpha |\alpha >0\}`$.

Substituting the signal probability distribution (38) into the general expression (37) for $`\mathrm{\Lambda }_{\text{ML}}`$ yields a Gaussian integral which is straightforward to evaluate. The result is

$`\mathrm{\Lambda }_{\text{ML}}^\text{G}`$ $`=`$ $`\underset{\alpha >0}{\mathrm{max}}\underset{\sigma _10}{\mathrm{max}}\underset{\sigma _20}{\mathrm{max}}\{{\displaystyle \frac{\overline{\sigma }_1\overline{\sigma }_2}{\sqrt{\sigma _1^2\sigma _2^2+\sigma _1^2\alpha ^2+\sigma _2^2\alpha ^2}}}`$
$`\times `$ $`\mathrm{exp}[{\displaystyle \frac{\frac{\overline{\sigma }_1^2}{\sigma _1^4}+\frac{\overline{\sigma }_2^2}{\sigma _2^4}+\frac{2\overline{\alpha }^2}{\sigma _1^2\sigma _2^2}}{2\left(\frac{1}{\sigma _1^2}+\frac{1}{\sigma _2^2}+\frac{1}{\alpha ^2}\right)}}{\displaystyle \frac{\overline{\sigma }_1^2}{2\sigma _1^2}}{\displaystyle \frac{\overline{\sigma }_2^2}{2\sigma _2^2}}+1]\}^N,`$
where

$$\overline{\alpha }^2=\frac{1}{N}\underset{k=1}{\overset{N}{}}h_1^kh_2^k,$$
(40)
and we have appended a superscript G on $`\mathrm{\Lambda }_{\text{ML}}^\text{G}`$ to indicate the maximum likelihood detection statistic for a Gaussian signal.

One can show that the maximum in Eq. (III.2) is achieved at $`\alpha =\widehat{\alpha }`$, $`\sigma _1=\widehat{\sigma }_1`$, and $`\sigma _2=\widehat{\sigma }_2`$, where

$`\widehat{\alpha }^2`$ $`=`$ $`\overline{\alpha }^2\theta (\overline{\alpha }^2),`$ (41)
$`\widehat{\sigma }_i^2`$ $`=`$ $`(\overline{\sigma }_i^2\widehat{\alpha }^2)\theta \left(\overline{\sigma }_i^2\widehat{\alpha }^2\right),`$ (42)
for $`i=1,2`$, and $`\overline{\sigma }_1`$ and $`\overline{\sigma }_2`$ are given by Eq. (35). Here $`\theta (x)`$ is the step function (5). The quantities (41) and (42) are the maximum likelihood estimators for the variance $`\alpha ^2`$ of the signal and the variances $`\sigma _1^2`$ and $`\sigma _2^2`$ of the noise in the two detectors. The step functions in Eqs. (41) and (42) arise as a result of the bounds of the maximization in Eq. (37).

The corresponding detection statistic is, from Eq. (III.2) 777To simplify the formula for $`\mathrm{\Lambda }_{\text{ML}}^\text{G}`$ we assume that $`\overline{\sigma }_i^2\overline{\alpha }^2>0`$. This will be true for any realistic value of $`N`$ since $`\overline{\sigma }_i^2\overline{\alpha }^2=\sigma _{i,\mathrm{true}}^2+O(1/\sqrt{N})`$, where $`\sigma _{i,\mathrm{true}}`$ is the true value of $`\sigma _i`$ and the second term describes the statistical fluctuations.,

$$\mathrm{\Lambda }_{\text{ML}}^\text{G}=\left[1\frac{\widehat{\alpha }^4}{\overline{\sigma }_1^2\overline{\sigma }_2^2}\right]^{N/2}.$$
(43)
The cross-correlation statistic $`\mathrm{\Lambda }_{\text{CC}}`$ can be obtained from $`\mathrm{\Lambda }_{\text{ML}}^\text{G}`$ via a monotonic transformation which preserves false dismissal versus false alarm curves \[cf. Eq. (21) above\]:

$`\mathrm{\Lambda }_{\text{CC}}`$ $`=`$ $`\sqrt{1(\mathrm{\Lambda }_{\text{ML}}^\text{G})^{2/N}}`$ (44)
$`=`$ $`{\displaystyle \frac{\widehat{\alpha }^2}{\overline{\sigma }_1\overline{\sigma }_2}}.`$
Note that if we had assumed the noise parameters $`𝐯_n=(\sigma _1,\sigma _2)`$ were known, and derived a statistic from Eq. (29) rather than Eq. (28), we would have found instead the detection statistic $`\stackrel{~}{\mathrm{\Lambda }}_{\text{ML}}^\text{G}=\overline{\mathrm{\Lambda }}_{\text{ML}}^\text{G}\theta (\overline{\mathrm{\Lambda }}_{\text{ML}}^\text{G})`$, where

$$\overline{\mathrm{\Lambda }}_{\text{ML}}^\text{G}=\overline{\alpha }^2+\frac{1}{2}\left[\frac{\sigma _2^2}{\sigma _1^2}(\overline{\sigma }_1^2\sigma _1^2)+\frac{\sigma _1^2}{\sigma _2^2}(\overline{\sigma }_2^2\sigma _2^2)\right],$$
(45)
which is different from the standard cross-correlation statistic. This non-standard result is obtained because of the unrealistic assumption that the noise parameters $`𝐯_n=(\sigma _1,\sigma _2)`$ are known. Different derivations of the result (45) can be found in Refs. sam joe ; robust gaussian II .

It is often useful to characterize the “strength” of a stochastic background in terms of the signal-to-noise ratio of the cross-correlation statistic (44), which we now define. First note that for large $`N`$, the fractional fluctuations in $`\widehat{\alpha }^2`$ will be much larger than those in $`\overline{\sigma }_1\overline{\sigma }_2`$ 888This is true at fixed signal-to-noise ratio $`\rho `$.. For the purpose of defining the signal-to-noise ratio, we assume that $`N`$ is large enough that $`\overline{\sigma }_1`$ and $`\overline{\sigma }_2`$ in Eq. (44) can be taken to be independent of $`h`$, so that $`\mathrm{\Lambda }_{\text{CC}}`$ and $`\widehat{\alpha }^2`$ are equivalent detection statistics. We also use $`\overline{\alpha }^2`$ instead of $`\widehat{\alpha }^2`$ in the computations that follow, as is conventional when defining signal-to-noise ratios. If a signal is present, then the expected value of $`\overline{\alpha }^2`$ is, from Eqs. (10), (31)–(33), (38) and (40),

$$\overline{\alpha }^2=\alpha ^2.$$
(46)
If no signal is present, so that $`\alpha ^2=0`$, then the fluctuations in $`\overline{\alpha }^2`$ are given by

$$\mathrm{\Delta }\left(\overline{\alpha }^2\right)=\frac{\sigma _1\sigma _2}{\sqrt{N}}.$$
(47)
The signal-to-noise ratio $`\rho `$ is defined to be the ratio of these two quantities:

$$\rho =\frac{\alpha ^2\sqrt{N}}{\sigma _1\sigma _2}.$$
(48)
### III.3 Non-Gaussian signal

As mentioned in the introduction, the traditional assumption that a gravitational wave stochastic background will be Gaussian requires the individual events to be sufficiently frequent and uncorrelated. Our model for a non-Gaussian signal assumes instead that the events are infrequent.

Consider a collection of similar events generating a stochastic background $`𝒮`$. Let $`\xi `$ be the probability that, at any randomly chosen time, the waves from an event are arriving at the detectors. We assume that the time structure of individual events cannot be resolved by the detectors. That is, we assume that the events occur over timescales smaller than the detectors’ resolution time, as illustrated in Fig. 4.

We assume that the distribution of the amplitudes of the events is Gaussian with variance $`\alpha ^2`$. The probability distribution for the signal given the signal parameters $`(\xi ,\alpha )`$ is therefore given by

$`p_{𝒮^k|𝒱_s,𝒯}[s^k|(\xi ,\alpha ),1]`$ $`=`$ $`{\displaystyle \frac{\xi }{\sqrt{2\pi }\alpha }}\mathrm{exp}\left[{\displaystyle \frac{\left(s^k\right)^2}{2\alpha ^2}}\right]`$ (49)
$`+`$ $`(1\xi )\delta \left(s^k\right),`$
together with Eq. (33). Thus the signal model parameters are $`𝐯_s=(\xi ,\alpha )`$, which give respectively the “event probability” and “event variance” characterizing the stochastic background. The parameter space $`\mathrm{\Theta }_s`$ for this model is

$$\mathrm{\Theta }_s=\left\{(\xi ,\alpha )\right|0\xi 1\text{ and }\alpha 0\},$$
(50)
and the subset corresponding to a signal being present is

$$\mathrm{\Theta }_{s1}=\left\{(\xi ,\alpha )\right|0<\xi 1\text{ and }\alpha >0\}.$$
(51)
Note that our assumption that the time structure of events is not resolved by the detector is unrealistic. Detector resolution times can be as small as 0.1 ms in the case of ground-based detectors like LIGO 999For ground-based detectors, the effective resolution time in a cross-correlation between two detectors can be considerably longer than $`0.1`$ ms Allen Romano , which may help with this issue., and even supernova bursts are expected to have time scales $`10`$ ms waveform catalog ; new waveform catalog . It will be important for future studies to relax this assumption.

We now compute the maximum likelihood detection statistic $`\mathrm{\Lambda }_{\text{ML}}^{\text{NG}}`$ for our simple non-Gaussian signal model by substituting Eq. (49) into Eq. (37). This yields

$`\mathrm{\Lambda }_{\text{ML}}^{\text{NG}}`$ $`=`$ $`\underset{0<\xi 1}{\mathrm{max}}\underset{\alpha >0}{\mathrm{max}}\underset{\sigma _10}{\mathrm{max}}\underset{\sigma _20}{\mathrm{max}}{\displaystyle \underset{k=1}{\overset{N}{}}}\{{\displaystyle \frac{\overline{\sigma }_1\overline{\sigma }_2\xi }{\sqrt{\sigma _1^2\sigma _2^2+\sigma _1^2\alpha ^2+\sigma _2^2\alpha ^2}}}\mathrm{exp}[{\displaystyle \frac{\left(\frac{h_1^k}{\sigma _1^2}+\frac{h_2^k}{\sigma _2^2}\right)^2}{2\left(\frac{1}{\sigma _1^2}+\frac{1}{\sigma _2^2}+\frac{1}{\alpha ^2}\right)}}{\displaystyle \frac{\left(h_1^k\right)^2}{2\sigma _1^2}}{\displaystyle \frac{\left(h_2^k\right)^2}{2\sigma _2^2}}+1]`$ (52)
$`+`$ $`{\displaystyle \frac{\overline{\sigma }_1\overline{\sigma }_2}{\sigma _1\sigma _2}}(1\xi )\mathrm{exp}[{\displaystyle \frac{\left(h_1^k\right)^2}{2\sigma _1^2}}{\displaystyle \frac{\left(h_2^k\right)^2}{2\sigma _2^2}}+1]\}.`$
The values of $`\xi `$, $`\alpha ^2`$, $`\sigma _1^2`$, and $`\sigma _2^2`$ which achieve the maximum in Eq. (52) are, respectively, estimators of the signal’s Gaussianity parameter, the variance of the signal events, and the noise variances in the two detectors101010See Ref. MG9  for a derivation of a statistic similar to $`\mathrm{\Lambda }_{\text{ML}}^{\text{NG}}`$ and designed for the same non-Gaussian signals which is based on Eq. (29) rather than Eq. (28).. Note that if we evaluate Eq. (52) at $`\xi =1`$, rather than maximizing over $`\xi `$, we recover Eq. (III.2) and the statistic $`\mathrm{\Lambda }_{\text{ML}}^\text{G}`$.

We mention in passing an approximate version of the statistic (52) which is significantly easier to compute. Expanding the logarithm of the quantity to be maximized in Eq. (52) as a power series in $`\alpha ^2`$ to fourth order about $`\alpha ^2=0`$ yields the approximate statistic $`\widehat{\mathrm{\Lambda }}_{\text{ML}}^{\text{NG}}`$ given by

$`\mathrm{ln}\widehat{\mathrm{\Lambda }}_{\text{ML}}^{\text{NG}}`$ $`=`$ $`\underset{0<\xi 1}{\mathrm{max}}\underset{\alpha >0}{\mathrm{max}}\underset{\sigma _10}{\mathrm{max}}\underset{\sigma _20}{\mathrm{max}}{\displaystyle \underset{n=0}{\overset{4}{}}}{\displaystyle \underset{l=0}{\overset{8}{}}}{\displaystyle \underset{m=0}{\overset{8}{}}}\left({\displaystyle \frac{\alpha ^2}{\sigma _1\sigma _2}}\right)^n`$ (53)
$`\times `$ $`C_{nlm}(\xi ,\sigma _1,\sigma _2){\displaystyle \underset{k=1}{\overset{N}{}}}(h_1^k)^l(h_2^k)^m,`$
where the coefficients $`C_{nlm}(\xi ,\sigma _1^2,\sigma _2^2)`$ vanish unless $`l+m`$ is even and $`l+m8`$. In evaluating the statistic (53), one can first evaluate the 24 sums

$$\underset{k=1}{\overset{N}{}}(h_1^k)^l(h_2^k)^m$$
(54)
for the required values of $`l`$ and $`m`$, and subsequently numerically maximize over the parameters $`\xi `$, $`\alpha `$, $`\sigma _1`$, and $`\sigma _2`$. Thus the length-$`N`$ sums need only be performed once, rather than each time one tries a new set of values for $`\xi `$, $`\alpha `$, $`\sigma _1`$, and $`\sigma _2`$. Therefore the computational cost of $`\widehat{\mathrm{\Lambda }}_{\text{ML}}^{\text{NG}}`$ is only about an order of magnitude greater than that of the cross correlation statistic $`\mathrm{\Lambda }_{\text{CC}}`$, and this statistic may be useful to explore.

We now derive the signal-to-noise ratio $`\rho `$ for the cross-correlation statistic and for the non-Gaussian signal (49). If the signal is present, then from Eqs. (10), (33), (40), (41) and (49) the expected value of $`\overline{\alpha }^2`$ is

$$\overline{\alpha }^2=\xi \alpha ^2.$$
(55)
If no signal is present then the fluctuations in $`\overline{\alpha }^2`$ are given by

$$\mathrm{\Delta }\left(\overline{\alpha }^2\right)=\frac{\sigma _1\sigma _2}{\sqrt{N}}.$$
(56)
Therefore, taking the ratio of Eqs. (55) and (56), the signal-to-noise ratio $`\rho `$ is

$$\rho =\frac{\xi \alpha ^2\sqrt{N}}{\sigma _1\sigma _2}.$$
(57)
## IV Performance comparison

In this section we compare the performances of the cross-correlation statistic (44), the burst statistic (9), and the maximum likelihood statistic (52) for our model non-Gaussian signal described in Sec. III.3. The comparison is quantified in terms of the false alarm versus false dismissal curves, as discussed in Sec. II above. In Sec. IV.1 we discuss analytic predictions for these curves for the three different statistics. Section IV.2 describes our Monte Carlo simulation algorithm, and Secs. IV.3 and IV.4 describe the results.

### IV.1 Analytic computation of asymptotic behavior of statistics

We start by discussing the set of parameters on which the false dismissal versus false alarm curves can depend. As before, we assume two detectors with noise characterized by Eq. (31) with $`𝒱_n=(\sigma _1,\sigma _2)`$, and a non-Gaussian signal characterized by Eqs. (33) and (49) with $`𝒱_s=(\xi ,\alpha )`$. The curves for each statistic are given by some function

$$P_{\text{FD}}=P_{\text{FD}}(P_{\text{FA}},\xi ,\alpha ,\sigma _1,\sigma _2,N)$$
(58)
of the false alarm probability $`P_{\mathrm{FA}}`$, the Gaussianity parameter $`\xi `$, the rms amplitude $`\alpha `$ of events, the noise variances $`\sigma _1^2`$ and $`\sigma _2^2`$, and the number of data points $`N`$. We can simplify Eq. (58) by replacing $`\alpha `$ with the signal-to-noise ratio $`\rho `$ using the definition (57), and noting from dimensional analysis that $`P_{\text{FA}}`$ depends on $`\sigma _1`$ and $`\sigma _2`$ at fixed $`\rho `$ only through the ratio $`\sigma _1/\sigma _2`$. This gives

$$P_{\text{FD}}=P_{\text{FD}}(P_{\text{FA}},\xi ,\rho ,\sigma _1/\sigma _2,N).$$
(59)
For simplicity, we specialize to $`\sigma _1=\sigma _2`$ for the remainder of this paper. This implies that

$$P_{\text{FD}}=P_{\text{FD}}(P_{\text{FA}},\xi ,\rho ,N).$$
(60)
#### IV.1.1 Cross correlation statistic

The false dismissal versus false alarm curves for the cross-correlation statistic can be computed analytically in the large $`N`$ limit, as we now describe. Our derivation generalizes the analysis of Ref. Allen Romano  from Gaussian to non-Gaussian signals. For any detection statistic $`\mathrm{\Gamma }`$, we can express $`P_{\text{FA}}`$ and $`P_{\text{FD}}`$ in terms of the detection threshold $`\mathrm{\Gamma }_{}`$ as

$`P_{\text{FA}}(\mathrm{\Gamma }_{},\sigma _1,\sigma _2,N)`$ $`=`$ $`{\displaystyle _\mathrm{\Gamma }_{}^{\mathrm{}}}𝑑xp_{\mathrm{\Gamma }|𝒯}(x|0),`$ (61)
$`P_{\text{FD}}(\mathrm{\Gamma }_{},\xi ,\rho ,\sigma _1,\sigma _2,N)`$ $`=`$ $`1{\displaystyle _\mathrm{\Gamma }_{}^{\mathrm{}}}𝑑xp_{\mathrm{\Gamma }|𝒯}(x|1).`$
Here the definition of the random variable $`𝒯`$ is such that if $`𝒯=0`$ then no signal is present ($`\xi =\rho =0`$), and if $`𝒯=1`$ then a signal is present ($`\xi 0`$ and $`\rho 0`$); cf. Sec. II.1 above. Note that by eliminating $`\mathrm{\Gamma }_{}`$ between Eqs. (61) and (IV.1.1), we recover Eq. (58).

In the large $`N`$ limit, the distribution $`p_{\mathrm{\Lambda }_{\text{CC}}|𝒯}(x|t)`$ is a Gaussian by the central limit theorem, and the integrals (61) and (IV.1.1) can be evaluated analytically (see Appendix B) to give

$`P_{\text{FD}}(P_{\text{FA}},\xi ,\rho ,N)=1`$ (63)
$`{\displaystyle \frac{1}{2}}\mathrm{erfc}\left[{\displaystyle \frac{\mathrm{erfc}^1\left(2P_{\text{FA}}\right){\displaystyle \frac{\rho }{\sqrt{2}}}}{\sqrt{{\displaystyle \frac{\rho ^2}{N}}\left({\displaystyle \frac{3}{\xi }}1\right)+{\displaystyle \frac{2\rho }{\sqrt{N}}}+1}}}\right]+O\left({\displaystyle \frac{1}{\sqrt{N}}}\right).`$
Here the function $`\mathrm{erfc}(x)`$ (known as the compliment of the error function) is defined by

$$\mathrm{erfc}(x)=\frac{2}{\sqrt{\pi }}_x^{\mathrm{}}𝑑ye^{y^2},$$
(64)
and $`\mathrm{erfc}^1(x)`$ is the inverse of $`\mathrm{erfc}(x)`$. The formula (63) is valid only for $`P_{\mathrm{FA}}<1/2`$; $`P_{\text{FD}}`$ is undefined for $`1/2P_{\text{FA}}<1`$. In deriving Eq. (63), we assumed that the statistics $`\mathrm{\Lambda }_{\text{CC}}`$ and $`\widehat{\alpha }^2`$ are equivalent, and that the distribution for $`\overline{\alpha }^2`$ is Gaussian. Those assumptions are only valid up to fractional correction terms of order $`1/\sqrt{N}`$; hence the indicated correction term in Eq. (63).

In the regime where $`\rho ^2N\xi `$ in addition to $`N1`$, the result (63) simplifies to

$`P_{\text{FD}}(P_{\text{FA}},\xi ,\rho ,N)`$ $`=`$ $`1{\displaystyle \frac{1}{2}}\mathrm{erfc}\left[\mathrm{erfc}^1\left(2P_{\text{FA}}\right){\displaystyle \frac{\rho }{\sqrt{2}}}\right]`$
$`+`$ $`O\left({\displaystyle \frac{1}{\sqrt{N}}}\right)+O\left({\displaystyle \frac{\rho }{\sqrt{N}}}\right)+O\left({\displaystyle \frac{\rho ^2}{N\xi }}\right).`$
Note that the false dismissal versus false alarm relation (LABEL:specialized\_analytic) is independent of both $`N`$ and $`\xi `$. Sample curves from Eq. (LABEL:specialized\_analytic) are shown in Fig. 5.

The discontinuities at $`P_{\text{FA}}=1/2`$ are a result of the step functions in the definition (41) of $`\widehat{\alpha }^2`$.

#### IV.1.2 Burst statistic

By combining the definition (9) of the burst statistic together with the decomposition (6), the noise and signal distributions (31) and (49), and the change of variables (57) it is straightforward to derive the exact false alarm versus false dismissal relation. The result is given by

$`(1P_{\mathrm{FA}})^{1/N}=\mathrm{erf}\left({\displaystyle \frac{\mathrm{\Lambda }_{}}{\sqrt{2}}}\right)`$ (66)
and

$$P_{\mathrm{FD}}^{1/N}=\xi \mathrm{erf}\left[\frac{\mathrm{\Lambda }_{}}{\sqrt{2+\frac{2\rho }{\xi \sqrt{N}}}}\right]+(1\xi )\mathrm{erf}\left(\frac{\mathrm{\Lambda }_{}}{\sqrt{2}}\right),$$
(67)
where $`\mathrm{\Lambda }_{}`$ is the value of the threshold.

#### IV.1.3 Maximum likelihood statistic

We start by discussing the different regimes present in the space of signal parameters $`\xi `$, $`\rho `$ and $`N`$, treating the false alarm probability $`P_{\text{FA}}`$ as fixed. There are several different constraints on the three parameters $`\xi `$, $`\rho `$, and $`N`$ that define the regime in parameter space where we expect our maximum likelihood statistic to work well. First, it is clear that the total number of events $`\xi N`$ in the data set must be large compared to one:

$$\xi \frac{1}{N}.$$
(68)
Second, if the signal-to-noise ratio $`\alpha ^2/(\sigma _1\sigma _2)`$ of individual burst events is large compared to one, then one can detect the individual events using the burst statistic (9) and the method of this paper is not needed. From Eq. (48) we can write the constraint $`\alpha ^2/(\sigma _1\sigma _2)1`$ as

$$\xi \frac{\rho }{\sqrt{N}}.$$
(69)
A more precise version of this requirement can be obtained by noting that the detection threshold for the signal-to-noise ratio $`\alpha ^2/(\sigma _1\sigma _2)`$ is $`\sqrt{2\mathrm{ln}N}`$, since there are $`N`$ independent trials. This yields the constraint

$$\xi \frac{\rho }{\sqrt{2N\mathrm{ln}N}}.$$
(70)
The regime $`\xi \rho /\sqrt{2N\mathrm{ln}N}`$ is where the burst statistic $`\mathrm{\Lambda }_\text{B}`$ starts becoming as sensitive as the cross correlation statistic, as can be seen by combining Eqs. (LABEL:specialized\_analytic), (66) and (67) above. This behavior can also be seen in Figs. 1 and 2 above.

A third constraint on the space of signal parameters is derived as follows. Consider the statistic

$$\eta =\frac{1}{N}\underset{k=1}{\overset{N}{}}(h_1^k)^2(h_2^k)^2.$$
(71)
We can use this statistic to estimate the Gaussianity parameter $`\xi `$ in the following way. The mean value of $`\eta `$ when a signal is present is given by

$$\eta =3\xi \alpha ^4+\xi \alpha ^2(\sigma _1^2+\sigma _2^2)+\sigma _1^2\sigma _2^2,$$
(72)
and the variance when a signal is absent is

$$(\mathrm{\Delta }\eta )^2=\sigma _1^4\sigma _2^4\frac{8}{N}.$$
(73)
It follows from Eqs. (72), (73), and the relation $`\widehat{\alpha }^2=\xi \alpha ^2`$ that the estimator $`\widehat{\xi }`$ of $`\xi `$ defined by

$$\widehat{\xi }=\frac{3\widehat{\alpha }^4}{\eta \widehat{\alpha }^2(\widehat{\sigma }_1^2+\widehat{\sigma }_2^2)\widehat{\sigma }_1^2\widehat{\sigma }_2^2}$$
(74)
has a fractional accuracy of order

$$\frac{\mathrm{\Delta }\xi }{\xi }\frac{\xi \sqrt{N}}{\rho ^2}.$$
(75)
Now in the regime $`\mathrm{\Delta }\xi /\xi 1`$, we expect our maximum likelihood detection statistic to work well, since one’s first guess for a nonlinear statistic (71) can be used to detect the non-Gaussianity of the signal to high accuracy. In the regime $`\mathrm{\Delta }\xi /\xi 1`$, it is not obvious how the maximum likelihood detection statistic will perform, since it could have a performance much better than that of the statistic $`\eta `$. However, our Monte Carlo simulations \[Sec. IV.2 below\] and analytic computations \[Appendix C\] indicate that the maximum likelihood statistic does indeed perform poorly in the regime $`\mathrm{\Delta }\xi /\xi 1`$. Thus, our third constraint is $`\mathrm{\Delta }\xi /\xi 1`$, which from Eq. (75) can be written as

$$\xi \frac{\rho ^2}{\sqrt{N}}.$$
(76)
Our Monte Carlo simulations show that for $`\rho ^2/\sqrt{N}\xi 1`$, the maximum likelihood and cross-correlation statistics perform roughly equivalently, and that once $`\xi `$ becomes smaller than $`\rho ^2/\sqrt{N}`$, the maximum likelihood statistic starts to perform significantly better than the cross-correlation statistic; see Figs. 1 and 2 above.

In Appendix C we derive analytically the approximate expression (168) for the false dismissal probability for the maximum likelihood statistic, which we expect to be accurate up to corrections of order $`1/\rho ^4`$ or a few tens of percent. We also derive the expression (177) for the false alarm probability using a combination of analytical and numerical techniques. Combining these results gives the curves which are associated with the maximum likelihood statistic $`\mathrm{\Lambda }_{\text{ML}}^{\text{NG}}`$ and labeled “analytic” in Figs. 1, 2, 9, and 10.

### IV.2 Description of the Monte Carlo simulation algorithm

Next we describe our Monte Carlo simulations of the performances of the various statistics. We numerically estimate the false dismissal and false alarm probabilities $`P_{\text{FD}}`$ and $`P_{\text{FA}}`$ by conducting an ensemble of $`N_E`$ simulated experiments. For each experiment we simulate a detector output matrix, half of which have a signal present, and half of which do not. Since we know in advance whether or not a signal is present, we can easily estimate $`P_{\text{FA}}`$ and $`P_{\text{FD}}`$. More specifically, our algorithm for simulating false dismissal versus false alarm curves, for an arbitrary statistic $`\mathrm{\Gamma }`$, is as follows:

1. Choose values for $`\xi `$, $`\alpha `$, $`\sigma _1`$, $`\sigma _2`$, and $`N`$.
2. Choose the total number of trials $`N_E`$.
3. For $`r=1,2,\mathrm{},N_E/2`$:

1. Generate a data train $`h(\sigma _1,\sigma _2,N)`$ of noise only.
    2. Compute $`\mathrm{\Gamma }`$ and store result as $`\mathrm{\Gamma }_{r0}`$.
    3. Generate a data train $`h(\xi ,\alpha ,\sigma _1,\sigma _2,N)`$ which has a signal present.
    4. Compute $`\mathrm{\Gamma }`$ and store result as $`\mathrm{\Gamma }_{r1}`$.
4. Choose a discretization $`\mathrm{\Gamma }_j`$ of the set of thresholds, where $`j=1,2,\mathrm{},M`$.
5. Set $`P_{\text{FA}}(\mathrm{\Gamma }_j)=P_{\text{FD}}(\mathrm{\Gamma }_j)=0`$, for each $`j`$.
6. For $`r=1,2,\mathrm{},N_E/2`$:

1. for each $`j`$, if $`\mathrm{\Gamma }_{r0}>\mathrm{\Gamma }_j`$, increment $`P_{\text{FA}}(\mathrm{\Gamma }_j)`$ by $`2/N_E`$.
    2. for each $`j`$, if $`\mathrm{\Gamma }_{r1}\mathrm{\Gamma }_j`$, increment $`P_{\text{FD}}(\mathrm{\Gamma }_j)`$ by $`2/N_E`$.
7. Repeat steps 3-6 above several times to estimate the fluctuations in $`P_{\text{FA}}(\mathrm{\Gamma }_j)`$ and $`P_{\text{FD}}(\mathrm{\Gamma }_j)`$.

We use the above algorithm to simulate false dismissal versus false alarm curves for the three statistics $`\mathrm{\Lambda }_{\text{CC}}`$, $`\mathrm{\Lambda }_\mathrm{B}`$ and $`\mathrm{\Lambda }_{\text{ML}}^{\text{NG}}`$. The analytical expressions (63) and (66) – (67) for the cross-correlation and burst statistics are used as a check of the numerical method.

### IV.3 Simulation results

A family of simulated false dismissal versus false alarm curves for the cross correlation statistic $`\mathrm{\Lambda }_{\text{CC}}`$ and the maximum likelihood statistic $`\mathrm{\Lambda }_{\text{ML}}^{\text{NG}}`$ is shown in Fig. 6.

We see that at fixed $`\rho `$, as the Gaussianity $`\xi `$ of the signal decreases, $`\mathrm{\Lambda }_{\text{ML}}^{\text{NG}}`$ performs increasingly better than $`\mathrm{\Lambda }_{\text{CC}}`$. The curves for $`\mathrm{\Lambda }_{\text{CC}}`$ are almost indistinguishable from each other because $`\rho `$ is fixed, and the curves depend only on $`\rho `$ and not on $`\xi `$ for this detection statistic in the large $`N`$ limit \[cf. Eq. (LABEL:specialized\_analytic) above\].

If we maintain the same value for $`\rho `$ as in Fig. 6, but take $`\xi 0.03`$, the curves for $`\mathrm{\Lambda }_{\text{CC}}`$ and $`\mathrm{\Lambda }_{\text{ML}}^{\text{NG}}`$ cannot be distinguished from each other. We find in general that for *any* values of $`N`$, $`\sigma _1`$, $`\sigma _2`$, and $`\rho `$, as $`\xi 1`$, the false dismissal versus false alarm curves for $`\mathrm{\Lambda }_{\text{CC}}`$ and $`\mathrm{\Lambda }_{\text{ML}}^{\text{NG}}`$ cannot be distinguished from each other. Thus, the two statistics are nearly equivalent for Gaussian signals, as expected. However, for $`\xi 1`$, Fig. 6 demonstrates that $`\mathrm{\Lambda }_{\text{ML}}^{\text{NG}}`$ performs noticeably better than $`\mathrm{\Lambda }_{\text{CC}}`$.

We now discuss a comparison of the two statistics in terms of the minimum gravitational wave energy density necessary for detection, instead of in terms of the false dismissal versus false alarm curves. For a stochastic background with rms strain amplitude $`h_{\text{rms}}`$, we have $`\mathrm{\Omega }h_{\text{rms}}^2`$ Allen Review , where $`\mathrm{\Omega }`$ is the gravitational wave energy density. For our model signal (49) we have $`h_{\mathrm{rms}}^2\xi \alpha ^2`$, and comparing this with the formula $`\rho \xi \alpha ^2`$ from Eq. (57) shows that we can interpret the signal to noise ratio $`\rho `$ as the energy density in the stochastic background, even for non-Gaussian signals.

We compute the minimum detectable energy density or signal-to-noise ratio $`\rho _{\mathrm{detectable}}`$ as follows. First, we choose thresholds $`P_{\text{FA}}`$ and $`P_{\text{FD}}`$ for the false alarm and false dismissal probabilities. We refer to the pair $`(P_{\text{FA}},P_{\text{FD}})`$ as the *detection point*. For any statistic $`\mathrm{\Gamma }`$, the choice of detection point determines the detection threshold $`\mathrm{\Gamma }_{}`$, and inverting Eq. (60) gives the minimum detectable signal-to-noise ratio

$$\rho =\rho _{\text{detectable}}(P_{\text{FA}},P_{\text{FD}},\xi ,N),$$
(77)
as illustrated in Fig. 7.

For the cross-correlation statistic $`\mathrm{\Lambda }_{\text{CC}}`$ the result is, from Eq. (63),

$`\rho _{\text{detectable}}^{\text{CC}}`$ $`=`$ $`{\displaystyle \frac{2\sqrt{2}\gamma \left[1+\gamma \sqrt{2/N}\right]}{1+2\gamma ^2\left(1\frac{3}{\xi }\right)/N}}\left[1+O\left({\displaystyle \frac{1}{\sqrt{N}}}\right)\right]`$
$`=`$ $`2\sqrt{2}\gamma +O\left({\displaystyle \frac{1}{\sqrt{N}}}\right)+O\left({\displaystyle \frac{\gamma }{\sqrt{N}}}\right),`$
where $`\gamma =\mathrm{erfc}^1(2P_{\text{FA}})`$ and we have assumed that $`P_{\text{FA}}=P_{\text{FD}}`$. This relation is plotted in Fig. 8.

From the results of our simulations, we determine $`\rho _{\text{detectable}}(P_{\text{FA}},P_{\text{FD}},\xi ,N)`$ by numerically solving the equation

$$P_{\text{FD}}(P_{\text{FA}},\xi ,\rho ,N)P_{\text{FD}}=0$$
(80)
for $`\rho `$. Unfortunately, evaluating the function on the left hand side of Eq. (80) is computationally expensive. Each evaluation involves simulating the false dismissal versus false alarm curve which is itself a computationally intensive task. Moreover, it is only feasible for us to solve Eq. (80) for values of $`N`$ $`10^4`$ while a realistic detection scenario for ground based detectors would involve a year’s worth of data sampled at $`100`$ Hz for which $`N10^9`$. Therefore our conclusions about the applicability of the method to ground based detectors are based on our analytic results, as discussed in the Introduction.

Figure 9 shows the results obtained from numerically solving Eq. (80) for $`\rho _{\mathrm{detectable}}`$ for the parameter values $`\xi =0.02`$, $`P_{\text{FA}}=P_{\text{FD}}=0.1`$, and Fig. 10 shows the corresponding results for $`\xi =4.3\times 10^3`$. For the cross-correlation statistic, the results are in good agreement with the analytic prediction (LABEL:analytic\_detectable).

Figure 1 shows the minimum detectable energy density as a function of the Gaussianity parameter $`\xi `$ for $`N=10^4`$ (corresponding to space based detectors), for the cross-correlation and maximum likelihood statistics and also for the burst statistic (9). We again use the values $`P_{\text{FA}}=P_{\text{FD}}=0.1`$. The figure shows that the maximum likelihood statistic performs better than the other statistics by a factor which is roughly 3 for $`\xi `$ of order 1%. For smaller values of $`\xi `$, the maximum likelihood performs increasingly better than the cross-correlation statistic, but is eventually comparable to the burst statistic. Thus the maximum likelihood statistic gives an improvement in sensitivity to backgrounds composed of roughly $`10`$ to $`10^3`$ events per year.

Figure 2 is a similar plot, without the Monte Carlo simulation results, for $`N=10^9`$ (corresponding to ground based detectors). Here we use $`P_{\text{FA}}=P_{\text{FD}}=0.01`$. The results are similar to those in Fig. 1, except that here the gain in sensitivity occurs in the band $`10^5<\xi <10^3`$. This band corresponds to $`10^4`$-$`10^6`$ events per year.

### IV.4 Parameter estimation

The computation of the maximum likelihood statistic also serves to measure the parameters of the signal. The statistic $`\mathrm{\Lambda }_{\text{ML}}^{\text{NG}}`$, from Eq. (52), can be written as

$$\mathrm{\Lambda }_{\text{ML}}^{\text{NG}}=\underset{0<\xi 1}{\mathrm{max}}\underset{\alpha ^2>0}{\mathrm{max}}\underset{\sigma _1^20}{\mathrm{max}}\underset{\sigma _2^20}{\mathrm{max}}\lambda (\xi ,\alpha ^2,\sigma _1^2,\sigma _2^2).$$
(81)
The point $`(\widehat{\xi },\widehat{\alpha }^2,\widehat{\sigma }_1^2,\widehat{\sigma }_2^2)`$ where this maximum is achieved is the maximum likelihood estimator for $`(\xi ,\alpha ^2,\sigma _1^2,\sigma _2^2)`$. In Fig. 11 we show contours of the function $`\mathrm{ln}\lambda `$ for a strong ($`\rho =20`$) signal.

This figure shows that both $`\xi `$ and $`\alpha ^2`$ can be measured with good accuracy.

Note that the main benefit of using $`\mathrm{\Lambda }_{\mathrm{ML}}^{\mathrm{NG}}`$ is that it allows us to detect signals that are too weak to be seen using $`\mathrm{\Lambda }_{\text{CC}}`$. Using $`\mathrm{\Lambda }_{\text{ML}}^{\text{NG}}`$ also allows one to test if a detected signal is Gaussian, as obtained above, but this is not the main benefit of the method, as there are other, simpler, methods to test for non-Gaussianity.

## V Conclusions

The use of our maximum likelihood statistic in searches for a non-Gaussian background gives a gain in sensitivity over the standard cross-correlation statistic. Figures 1 and 2 show that the gain factor can be significant for sufficiently non-Gaussian signals. However, computing the maximum likelihood statistic requires significantly more computational power than the cross-correlation statistic.

The analysis presented here must be generalized in several ways before being usable in gravitational wave detectors. These generalizations, listed in order of importance, are:

* Our signal model (49) assumes a Gaussian distribution of amplitudes of the burst events. This assumption simplified our analysis and resulted in a statistic with the useful property of being nearly equivalent to the cross-correlation statistic in the Gaussian signal limit. In practice however, the distribution of the events should instead be based on the candidate sources. For example, a popcorn-like stochastic background produced by a spatially uniform distribution of standard-candle sources out to some maximum redshift would have a signal distribution of the form (49) with the Gaussian term replaced by a term proportional to $`s^4\theta (ss_{\mathrm{min}})`$, where $`\theta `$ is the step function and $`s_{\mathrm{min}}`$ is a cutoff signal strength.
* One should allow the burst durations to be longer than the detector resolution time. For this situation one possibility would be to preprocess the data with a lowpass filter, and then apply the techniques developed here. Another possibility would be to try to combine the analysis of this paper with the excess power detection method of Ref. excess power .
* Real detector noise always contains non-Gaussian components, so one needs to generalize the analysis to allow for this. Such a generalization for a Gaussian stochastic background can be found in Refs. robust gaussian ; robust gaussian II .
* It would be useful to consider a more general signal model which consists of a superposition of a Gaussian background and a non-Gaussian background, since the true gravitational wave background might consist of such a superposition.
* The analysis needs to be generalized to allow for colored detector noise, and separated, misaligned detectors. This generalization should be fairly straightforward.

###### Acknowledgements.

We thank Wolfgang Tichy, Tom Loredo, Teviet Creighton, and Bernard Whiting for helpful discussions, and the web site google.com for providing useful references on the generalized central limit theorem. The analytic computations in Appendix C were carried out using the software package Mathematica. This work was supported in part by National Science Foundation awards PHY-9722189 and PHY-0140209, the Alfred P. Sloan foundation, the Radcliffe Institute for Advanced Study, and the NASA/New York Space Grant Consortium.

## Appendix A General form of the likelihood ratio

In this appendix we give two derivations of the general formula (24) for the likelihood ratio. The first derivation is based on Eq. (22) while the second is based on Eq. (23). We also derive the formula (27) for the posterior probability density $`p_{𝒱_s,𝒱_n|𝒯}^{(1)}(𝐯_s,𝐯_n|1)`$.

### A.1 First derivation

We can derive Eq. (24) by using the total probability theorem to expand the distributions in the numerator and denominator of Eq. (22). Note that all distributions in this derivation are priors.

First expand $`p_{}(h)`$ just in terms of the random variable $`𝒯`$

$$p_{}(h)=P_𝒯(1)p_{|𝒯}(h|1)+P_𝒯(0)p_{|𝒯}(h|0).$$
(82)
Expanding $`p_{}(h)`$ in terms of all the degrees of freedom yields

$`p_{}(h)={\displaystyle \underset{t=0}{\overset{1}{}}}{\displaystyle _{\mathrm{\Theta }_s}}d^{Q_s}v_s{\displaystyle d^{ND}s_{\mathrm{\Theta }_n}d^{Q_n}v_n}`$ (83)
$`\times p_{|𝒯,𝒱_s,𝒮,𝒱_n}(h|t,𝐯_s,s,𝐯_n)p_{𝒯,𝒱_s,𝒮,𝒱_n}(t,𝐯_s,s,𝐯_n),.`$
The ratio of the coefficients of $`P_𝒯(1)`$ and $`P_𝒯(0)`$ in Eq. (83) will give the general expression for the likelihood ratio by Eq. (22).

The conditional distribution for $``$ in Eq. (83) can be translated into a conditional distribution for $`𝒩`$. From Eq. (10) it follows that

$$p_{|𝒮}(h|s)=p_{𝒩+𝒮|𝒮}(h|s)=p_{𝒩|𝒮}(hs|s),$$
(84)
and since $`𝒮`$ and $`𝒩`$ are statistically independent we obtain

$$P_{|𝒮}(h|s)=P_𝒩(hs).$$
(85)
Generalizing this argument gives

$$p_{|𝒯,𝒱_s,𝒮,𝒱_n}(h|t,𝐯_s,s,𝐯_n)=p_{𝒩|𝒱_n}(hs|𝐯_n),$$
(86)
since a priori $`𝒯`$, $`𝒱_s`$, and $`𝒮`$ are statistically independent of $`𝒩`$ and $`𝒱_n`$. For the same reason we can write the joint distribution that appears in Eq. (83) as

$$p_{𝒯,𝒱_s,𝒮,𝒱_n}(t,𝐯_s,s,𝐯_n)=p_{𝒯,𝒱_s,𝒮}(t,𝐯_s,s)p_{𝒱_n}(𝐯_n).$$
(87)
Substituting Eqs. (86) and (87) into Eq. (83) yields

$`p_{}(h)`$ $`=`$ $`{\displaystyle \underset{t=0}{\overset{1}{}}}{\displaystyle _{\mathrm{\Theta }_s}}d^{Q_s}v_s{\displaystyle d^{ND}s_{\mathrm{\Theta }_n}d^{Q_n}v_n}`$
$`\times `$ $`p_{𝒩|𝒱_n}(hs|𝐯_n)p_{𝒯,𝒱_s,𝒮}(t,𝐯_s,s)p_{𝒱_n}(𝐯_n).`$
We can also rewrite the distribution $`p_{𝒯,𝒱_s,𝒮}(t,𝐯_s,s)`$ as

$$p_{𝒯,𝒱_s,𝒮}(t,𝐯_s,s)=p_{𝒮|𝒱_s,𝒯}(s|𝐯_s,t)p_{𝒱_s|T}(𝐯_s,t)P_𝒯(t),$$
(89)
by Eq. (14). Substituting Eq. (89) into Eq. (A.1) and explicitly evaluating the sum over $`t`$ yields

$`p_{}(h)`$ $`=`$ $`P_𝒯(1){\displaystyle _{\mathrm{\Theta }_{s1}}}d^{Q_s}v_s{\displaystyle d^{ND}s_{\mathrm{\Theta }_n}d^{Q_n}v_np_{𝒩|𝒱_n}(hs|𝐯_n)p_{𝒱_n}(𝐯_n)p_{𝒮|𝒱_s,𝒯}(s|𝐯_s,1)p_{𝒱_s|𝒯}(𝐯_𝐬|1)}`$ (90)
$`+`$ $`P_𝒯(0){\displaystyle _{\mathrm{\Theta }_n}}d^{Q_n}v_np_{𝒩|𝒱_s}(h|𝐯_n).`$
Here we have used the following relations:

$`p_{𝒮|𝒱_s,𝒯}(s|𝐯_s\mathrm{\Theta }_{s1},0)`$ $`=`$ $`\delta ^{ND}(s)`$ (91)
$`p_{𝒱_s|𝒯}(𝐯_𝐬\mathrm{\Theta }_{s0}|1)`$ $`=`$ $`0`$ (92)
$`p_{𝒱_s|𝒯}(𝐯_𝐬\mathrm{\Theta }_{s1}|0)`$ $`=`$ $`0`$ (93)
$`{\displaystyle _{\mathrm{\Theta }_{s0}}}d^{Q_s}v_sp_{𝒱_s|𝒯}(𝐯_s|0)`$ $`=`$ $`1.`$ (94)
By comparing Eqs. (82) and (90) we can read off the distributions $`p_{|𝒯}(h|t)`$ and construct Eq. (24) from Eq. (22). Note that the expression (22) is independent of the space $`\mathrm{\Theta }_{s0}`$ of signal parameters corresponding to “no signal present”.

### A.2 Second derivation

Here we derive Eq. (24), and also Eq. (27), from Eq. (23). Consider the distribution

$$p_{𝒯,𝒱_s,𝒱_n|}(1,𝐯_s,𝐯_n|h)=\frac{p_{𝒯,𝒱_s,𝒱_n,}(1,𝐯_s,𝐯_n,h)}{p_{}(h)}.$$
(95)
We will justify Eq. (24) by the defining relation Eq. (23), which explicitly refers to priors and posteriors. Therefore we now append the appropriate superscripts as bookkeeping devices. Eq. (95) then reads

$$p_{𝒯,𝒱_s,𝒱_n}^{(1)}(1,𝐯_s,𝐯_n)=\frac{p_{𝒯,𝒱_s,𝒱_n,}^{(0)}(1,𝐯_s,𝐯_n,h)}{p_{}^{(0)}(h)}.$$
(96)
Using the expansion of $`p_{}(h)`$ given by Eq. (90), and what we will justify is the likelihood ratio $`\mathrm{\Lambda }`$ given by Eq. (24), we have

$$p_{𝒯,𝒱_s,𝒱_n}^{(1)}(1,𝐯_s,𝐯_n)=\frac{\left[\frac{p_{𝒯,𝒱_s,𝒱_n,}^{(0)}(1,𝐯_s,𝐯_n,h)}{_{\mathrm{\Theta }_n}d^{Q_n}v_n^{}p_{|𝒱_n}^{(0)}(h|𝐯_n)p_{𝒱_n}^{(0)}(𝐯_n)}\right]}{\mathrm{\Lambda }P^{(0)}+1P^{(0)}}.$$
(97)
Expanding the uppermost numerator in Eq. (97) over $`𝒮`$ by the total probability theorem gives

$`p_{𝒯,𝒱_s,𝒱_n,}^{(0)}(1,𝐯_s,𝐯_n,h)`$ $`=`$ $`{\displaystyle d^{ND}s}`$
$`\times `$ $`p_{𝒯,𝒱_s,𝒱_n,,𝒮}^{(0)}(1,𝐯_s,𝐯_𝐧,h,s),`$
and rewriting this gives

$`p_{𝒯,𝒱_s,𝒱_n,,𝒮}^{(0)}(1,𝐯_s,𝐯_𝐧,h,s)=p_{𝒩|𝒱_s}^{(0)}(hs|𝐯_s)`$
$`\times p_{𝒱_n}^{(0)}(𝐯_n)p_{𝒮|𝒱_s,𝒯}^{(0)}(s|𝐯_s,1)p_{𝒱_s|𝒯}^{(0)}(𝐯_s,1)P^{(0)}.`$ (99)
After putting Eq. (A.2) into Eq. (A.2), substitute the result into Eq. (97). Using $`\mathrm{\Lambda }(𝐯_s,𝐯_n)`$ given by Eq. (26) then yields

$$P^{(1)}p_{𝒱_s,𝒱_n|𝒯}^{(1)}(𝐯_s,𝐯_n|1)=\frac{\mathrm{\Lambda }(𝐯_s,𝐯_n)P^{(0)}}{\mathrm{\Lambda }P^{(0)}+1P^{(0)}}.$$
(100)
On the left hand side of Eq. (100) we have used

$$p_{𝒯,𝒱_s,𝒱_n}^{(1)}(1,𝐯_s,𝐯_n)=P^{(1)}p_{𝒱_s,𝒱_n|𝒯}^{(1)}(𝐯_s,𝐯_n|1).$$
(101)
Integrate Eq. (100) over $`\mathrm{\Theta }_n`$ and $`\mathrm{\Theta }_{s1}`$ using Eq. (25) and the normalization requirement

$$_{\mathrm{\Theta }_{s1}}d^{Q_s}v_s_{\mathrm{\Theta }_n}d^{Q_n}v_np_{𝒱_s,𝒱_n|𝒯}^{(1)}(𝐯_s,𝐯_n|1)=1$$
(102)
to get

$$P^{(1)}=\frac{\mathrm{\Lambda }P^{(0)}}{\mathrm{\Lambda }P^{(0)}+1P^{(0)}}.$$
(103)
Use Eq. (103) and Eq. (100) to form the ratio on the left hand side of Eq. (27) . This justifies Eq. (27).

Integrate Eq. (27) over $`\mathrm{\Theta }_n`$ and $`\mathrm{\Theta }_{s1}`$ using Eq. (25) and Eq.( 102) to see that the defining relation Eq. (23) is satisfied and thus Eq. (24) is justified.

## Appendix B Analytical expressions for false dismissal versus false alarm curves for cross-correlation statistic

This appendix derives the analytical form (63) of the false dismissal versus false alarm curves for the cross-correlation statistic $`\mathrm{\Lambda }_{\text{CC}}`$ in the large $`N`$ limit, for both Gaussian and non-Gaussian signals. A derivation for Gaussian signals can be found in Sec. IV of Ref. Allen Romano .

As noted in Sec. III.2, the statistics $`\mathrm{\Lambda }_{\text{CC}}`$ and $`\widehat{\alpha }^2`$ are equivalent in the large $`N`$ limit. Thus, in this limit, the false dismissal versus false alarm curves can be found by evaluating Eqs. (61) and (IV.1.1) with $`\mathrm{\Gamma }`$ replaced by $`\widehat{\alpha }^2`$. The relation (41) between the statistics $`\overline{\alpha }^2`$ and $`\widehat{\alpha }^2`$ implies the following relation between their probability distributions $`p_{\widehat{\alpha }^2|𝒯}(x|t)`$ and $`p_{\overline{\alpha }^2|𝒯}(x|t)`$:

$$p_{\widehat{\alpha }^2|𝒯}(x|t)=\theta (x)p_{\overline{\alpha }^2|𝒯}(x|t)+\delta (x)_{\mathrm{}}^0𝑑yp_{\overline{\alpha }^2|𝒯}(y|t).$$
(104)
Inserting this formula into Eqs. (61) and (IV.1.1) gives

$`P_{\text{FA}}(\widehat{\alpha }_{}^2)`$ $`=`$ $`\{\begin{array}{cc}{\displaystyle _{\widehat{\alpha }_{}^2}^{\mathrm{}}}𝑑xp_{\overline{\alpha }^2|𝒯}(x|0)\hfill & \text{ if }\widehat{\alpha }_{}^2>0\hfill \\ 1\hfill & \text{ if }\widehat{\alpha }_{}^20\hfill \end{array},`$ (107)
$`P_{\text{FD}}(\widehat{\alpha }_{}^2)`$ $`=`$ $`\{\begin{array}{cc}1{\displaystyle _{\widehat{\alpha }_{}^2}^{\mathrm{}}}𝑑xp_{\overline{\alpha }^2|𝒯}(x|1)\hfill & \text{ if }\widehat{\alpha }_{}^2>0\hfill \\ 0\hfill & \text{ if }\widehat{\alpha }_{}^20\hfill \end{array}.`$ (110)
In the large $`N`$ limit, the distribution $`p_{\overline{\alpha }_s^2|𝒯}(x|t)`$ must be Gaussian by the central limit theorem, and therefore this distribution is characterized entirely by its mean $`\overline{\alpha }_t^2`$ and variance $`[\mathrm{\Delta }(\overline{\alpha }_t^2)]^2`$. From Eqs. (10), (33), (40), (41) and (49), these are given by

$`\overline{\alpha }_0^2`$ $`=`$ $`0`$ (112)
$`\mathrm{\Delta }(\overline{\alpha }_0^2)`$ $`=`$ $`{\displaystyle \frac{\sigma _1\sigma _2}{\sqrt{N}}}`$ (113)
$`\overline{\alpha }_1^2`$ $`=`$ $`\xi \alpha ^2`$ (114)
$`\mathrm{\Delta }(\overline{\alpha }_1^2)`$ $`=`$ $`\sqrt{{\displaystyle \frac{\xi \alpha ^4(3\xi )+\xi \alpha ^2(\sigma _1^2+\sigma _2^2)+\sigma _1^2\sigma _2^2}{N}}}.`$
Substituting Gaussian distributions, with means and variances determined by Eqs. (112)-(LABEL:var1), into Eqs. (107) and (LABEL:Pfd) yields

$`P_{\text{FA}}(\widehat{\alpha }_{}^2,\sigma _1,\sigma _2,N)`$ $`=`$ $`\{\begin{array}{cc}{\displaystyle \frac{1}{2}}\mathrm{erfc}\left({\displaystyle \frac{\widehat{\alpha }_{}^2}{\sigma _1\sigma _2}}\sqrt{{\displaystyle \frac{N}{2}}}\right)\hfill & \text{ if }\widehat{\alpha }_{}^2>0\hfill \\ 1\hfill & \text{ if }\widehat{\alpha }_{}^20\hfill \end{array},`$ (118)
$`P_{\text{FD}}(\widehat{\alpha }_{}^2,\xi ,\alpha ,\sigma _1,\sigma _2,N)`$ $`=`$ $`\{\begin{array}{cc}1{\displaystyle \frac{1}{2}}\mathrm{erfc}\left[\left(\widehat{\alpha }_{}^2\xi \alpha ^2\right)\sqrt{{\displaystyle \frac{N}{2\left[\xi \alpha ^4(3\xi )+\xi \alpha ^2(\sigma _1^2+\sigma _2^2)+\sigma _1^2\sigma _2^2\right]}}}\right]\hfill & \text{ if }\widehat{\alpha }_{}^2>0\hfill \\ 0\hfill & \text{ if }\widehat{\alpha }_{}^20\hfill \end{array}.`$ (121)
If we now eliminate $`\widehat{\alpha }_{}^2`$ between Eqs. (118) and (LABEL:analytic2), change variables from $`\alpha `$ to $`\rho `$ using Eq. (57), and set $`\sigma _1=\sigma _2`$, we obtain Eq. (63).

## Appendix C Asymptotic behavior of maximum likelihood statistic

In this appendix we derive the large-$`N`$ behavior of the maximum likelihood statistic $`\mathrm{\Lambda }_{\mathrm{ML}}^{\mathrm{NG}}`$. From Eq. (52), we can write the statistic in the form

$$\mathrm{\Lambda }_{\mathrm{ML}}^{\mathrm{NG}}(h)=\mathrm{exp}\left[N(h)\right]$$
(123)
with

$$(h)=\underset{\sigma _1,\sigma _2,\xi ,\alpha }{\mathrm{max}}g(\sigma _1,\sigma _2,\xi ,\alpha ,h)$$
(124)
where

$$g=\frac{1}{N}\underset{k=1}{\overset{N}{}}g_k(\sigma _1,\sigma _2,\xi ,\alpha ),$$
(125)
and the function $`g_k=g_k(\sigma _1,\sigma _2,\xi ,\alpha )`$ is given by

$$e^{g_k}=\xi A_k(\alpha )+(1\xi )A_k(0)$$
(126)
with

$`A_k(\alpha )`$ $`=`$ $`\mathrm{exp}\left[{\displaystyle \frac{\left(\frac{h_1^k}{\sigma _1^2}+\frac{h_2^k}{\sigma _2^2}\right)^2}{2\left(\frac{1}{\sigma _1^2}+\frac{1}{\sigma _2^2}+\frac{1}{\alpha ^2}\right)}}{\displaystyle \frac{\left(h_1^k\right)^2}{2\sigma _1^2}}{\displaystyle \frac{\left(h_2^k\right)^2}{2\sigma _2^2}}+1\right]`$ (127)
$`\times {\displaystyle \frac{\overline{\sigma }_1\overline{\sigma }_2}{\sqrt{\sigma _1^2\sigma _2^2+\sigma _1^2\alpha ^2+\sigma _2^2\alpha ^2}}}.`$
We denote by $`\stackrel{~}{\sigma }_1`$, $`\stackrel{~}{\sigma }_2`$, $`\stackrel{~}{\xi }`$ and $`\stackrel{~}{\alpha }`$ the “true” parameters governing the distribution of the quantities $`h_1^k`$ and $`h_2^k`$ according to Eqs. (10), (31), (33), and (49), with untilded quantities replaced by the corresponding tilded quantities. \[These “true parameters” were denoted by $`\sigma _1`$, $`\sigma _2`$, $`\xi `$ and $`\alpha `$ in the body of the paper.\] We define $`\stackrel{~}{\rho }`$ to be the signal-to-noise ratio (57) with untilded quantities replaced by tilded quantities:

$$\stackrel{~}{\rho }\frac{\stackrel{~}{\xi }\stackrel{~}{\alpha }^2\sqrt{N}}{\stackrel{~}{\sigma }_1\stackrel{~}{\sigma }_2}.$$
(128)
For simplicity, in this appendix we restrict attention to the case $`\stackrel{~}{\sigma }_1=\stackrel{~}{\sigma }_2`$. Then, without loss of generality, we can take $`\stackrel{~}{\sigma }_1=\stackrel{~}{\sigma }_2=1`$ by rescaling our units of strain amplitude.

We discuss separately the computation of the false alarm and false dismissal probabilities, as different techniques are required to compute each.

### C.1 False dismissal probability

The false dismissal probability for the statistic (124) will be some function

$$P_{\mathrm{FD}}=P_{\mathrm{FD}}(_{},N,\stackrel{~}{\xi },\stackrel{~}{\rho })$$
(129)
of the threshold $`_{}`$ on $``$, the number of data points $`N`$, the Gaussianity parameter $`\stackrel{~}{\xi }`$ and signal-to-noise ratio $`\stackrel{~}{\rho }`$ of the signal. For applications to ground based detectors, we will have $`\stackrel{~}{\rho }`$ (a few), in order that the signal be detectable, $`N10^9`$, and $`10^3\stackrel{~}{\xi }1`$. Therefore it would be useful to find approximate analytic expressions for the false alarm probability in the limit of large $`N`$. There are actually several different, large $`N`$ regimes in the three dimensional parameter space $`(N,\stackrel{~}{\xi },\stackrel{~}{\rho })`$ that one might explore:

* The limit $`N\mathrm{}`$ with $`\stackrel{~}{\alpha }`$ and $`\stackrel{~}{\xi }`$ held fixed. This corresponds to fixing the stochastic background signal and going to a limit of long observation times. In this limit we have $`\stackrel{~}{\rho }\sqrt{N}`$ which diverges. This is not a very realistic limit to explore.
* The limit $`N\mathrm{}`$ with $`\stackrel{~}{\rho }`$ and $`\stackrel{~}{\xi }`$ held fixed. In this limit, the signal-to-noise ratio is held fixed, and correspondingly the amplitude $`\stackrel{~}{\alpha }`$ of the stochastic background signal goes to zero, from Eq. (128). This would be the most natural limit to explore. However, in this limit the statistical error $`\mathrm{\Delta }\stackrel{~}{\xi }`$ in our measurement of the Gaussianity parameter would diverge, from Eq. (75), and therefore in this limit we do not expect to be able to compute analytically the value of the parameter $`\xi `$ which achieves the maximum in Eq. (124). The analytic approximation methods which we discuss below do not work in this regime. \[In addition our Monte Carlo simulations show that the maximum likelihood statistic itself does not perform any better than the cross-correlation statistic in this regime, as discussed in the Introduction.\]
* The limit we actually explore is the limit $`N\mathrm{}`$ with $`\stackrel{~}{\xi }`$ fixed and $`\stackrel{~}{\rho }`$ scaling $`N^{1/4}`$, corresponding to $`\stackrel{~}{\alpha }N^{1/8}`$. The reason for our choosing to explore this particular limit is simply that it is amenable to analytic computations. Fractional corrections to our analytic results should scale like $`1/N`$ or as $`1/\stackrel{~}{\rho }^4`$. Since $`\stackrel{~}{\rho }`$ (a few) at the threshold for detection, the approximation should be good to $`10\%20\%`$ or so.

We now turn to a discussion of the computational technique. We write

$$\stackrel{~}{\alpha }=\stackrel{~}{\alpha }_0N^{1/8},$$
(130)
where $`\stackrel{~}{\alpha }_0`$ is independent of $`N`$. Correspondingly, from Eq. (49) we can write

$$s^k=N^{1/8}\widehat{s}^k,$$
(131)
where the distribution of $`\widehat{s}^k`$ is given by Eq. (49) with $`\xi `$ replaced by $`\stackrel{~}{\xi }`$ and $`\alpha `$ replaced by $`\stackrel{~}{\alpha }_0`$. In particular, the distribution of $`\widehat{s}^k`$ is independent of $`N`$. In computing the maximum over $`(\xi ,\alpha ,\sigma _1,\sigma _2)`$ in Eq. (124), it is useful change variables from $`\alpha `$ to $`\kappa `$ defined by

$$\kappa =\rho N^{1/4}=\frac{\xi \alpha ^2N^{1/4}}{\sigma _1\sigma _2},$$
(132)
which we expect to be independent of $`N`$ to leading order in the large $`N`$ limit. The value of the variable $`\kappa `$ that characterizes the signal is

$$\stackrel{~}{\kappa }=\stackrel{~}{\rho }N^{1/4}=\frac{\stackrel{~}{\xi }\stackrel{~}{\alpha }_0^2}{\stackrel{~}{\sigma }_1\stackrel{~}{\sigma }_2};$$
(133)
cf. Eqs. (130) and (132).

We now consider fixed realizations of the infinite sequences of random variables $`n_1^k`$, $`n_2^k`$ and $`\widehat{s}^k`$, and $`1k<\mathrm{}`$, and examine the limiting behavior of $`(h)`$ as $`N\mathrm{}`$. We compute this limiting behavior by substituting into the right hand side of Eq. (123) the relations

$$h_1^k=n_1^k+N^{1/8}\widehat{s}^kh_2^k=n_2^k+N^{1/8}\widehat{s}^k,$$
(134)
writing $`\alpha `$ in terms of $`\kappa `$ using Eq. (132), and expanding in powers of $`N^{1/8}`$. The result is an expression which can be written in terms of the sums $`Q_{abc}`$ defined by

$$Q_{abc}=\frac{1}{N}\underset{k=1}{\overset{N}{}}\left(\widehat{s}^k\right)^a\left(n_1^k\right)^b\left(n_2^k\right)^c,$$
(135)
where $`a`$, $`b`$, and $`c`$ are non-negative integers. From the central limit theorem we can write

$$Q_{abc}=\mu _{abc}+\frac{1}{\sqrt{N}}\mathrm{\Delta }_{abc},$$
(136)
where $`\mu _{abc}=Q_{abc}`$ are computable functions of $`\stackrel{~}{\xi }`$ and $`\stackrel{~}{\alpha }`$, and where the random variables $`(\mathrm{\Delta }_{100},\mathrm{\Delta }_{010},\mathrm{})`$ converge in distribution 111111See chapter 8 of Ref. Papoulis  for definitions of different notions of convergence for sequences of random variables.  as $`N\mathrm{}`$ to a multivariate Gaussian of zero mean whose variance-covariance matrix is independent of $`N`$. Thus, in particular the joint distribution of all $`\mathrm{\Delta }_{abc}`$’s is $`N`$-independent in limit that $`N\mathrm{}`$.

We define the vector

$$𝐯=(v^1,v^2,v^3,v^4)=(\xi ,\kappa ,\sigma _1^2,\sigma _2^2),$$
(137)
We denote the value of $`𝐯`$ that achieves the maximum in Eq. (124) by $`\widehat{𝐯}`$:

$$g(\widehat{𝐯})=\underset{𝐯}{\mathrm{max}}g(𝐯),$$
(138)
where $`\widehat{𝐯}=(\widehat{\xi },\widehat{\kappa },\widehat{\sigma _1^2},\widehat{\sigma _2^2})`$. These estimators satisfy a system of four equations 121212Here we are assuming that the maximum is achieved as a local maximum in the interior of the 4 dimensional parameter space. Cases when the maximum is achieved on the boundary are discussed below.

$$\frac{g}{v^l}|_{𝐯=\widehat{𝐯}}=0.$$
(139)
We solve Eq. (139) perturbatively. First assume that the estimators can be expanded in the form

$$\widehat{v^l}=\underset{j=0}{\overset{\mathrm{}}{}}\widehat{v^l}^{[j]}ϵ^j,$$
(140)
where for ease of notation we have defined $`ϵ=N^{1/8}`$. We define the expansion coefficients $`v^{l[j]}`$ analogously by an expansion of the form (140) but without the hats. Now using Eq. (136) the function $`g`$ can be expanded as a power series in $`ϵ`$ whose coefficients are functions of $`v^{l[k]}`$, $`\mu _{abc}`$, and $`\mathrm{\Delta }_{abc}`$:

$$g(𝐯)=\underset{j=0}{\overset{\mathrm{}}{}}g^{[j]}[v^{l[k]},\mu _{abc},\mathrm{\Delta }_{abc}]ϵ^j.$$
(141)
Substituting the expansions (140) and (141) into the condition (139) for a local extremum gives an infinite set of equations which must collectively be satisfied by the coefficients $`\widehat{v^l}^{[j]}`$

$$\frac{g^{[j]}}{v^{l[k]}}|_{v^{m[n]}=\widehat{v^m}^{[n]}}=0.$$
(142)
We solve these equations order by order to determine the coefficients $`\widehat{v^l}^{[j]}`$, and thereby justify a posteriori the ansatz (140).

We find that in order to compute the leading order expression for $``$, we must obtain the expansion for $`\widehat{\xi }`$ to zeroth order in $`ϵ`$, the expansion for $`\widehat{\kappa }`$ to fourth order in $`ϵ`$, and the expansions of $`\widehat{\sigma _1^2}`$ and $`\widehat{\sigma _2^2}`$ to sixth order in $`ϵ`$. The leading order results are

$`\widehat{\kappa }`$ $`=`$ $`\stackrel{~}{\kappa }+ϵ^2X+O(ϵ^3),`$ (143)
$`{\displaystyle \frac{1}{\widehat{\xi }}}`$ $`=`$ $`{\displaystyle \frac{1}{\stackrel{~}{\xi }}}+{\displaystyle \frac{Y}{\sqrt{6}\stackrel{~}{\kappa }^2}}+O(ϵ),`$ (144)
$`\widehat{\sigma _1^2}`$ $`=`$ $`1+O(ϵ^2),`$ (145)
$`\widehat{\sigma _2^2}`$ $`=`$ $`1+O(ϵ^2),`$ (146)
where

$`X=\mathrm{\Delta }_{011}`$ (147)
and

$`Y`$ $`=`$ $`{\displaystyle \frac{1}{8\sqrt{6}}}[4(\mathrm{\Delta }_{031}+\mathrm{\Delta }_{013})12(\mathrm{\Delta }_{002}+\mathrm{\Delta }_{020})`$ (148)
$`24\mathrm{\Delta }_{011}+\mathrm{\Delta }_{040}+\mathrm{\Delta }_{004}+6\mathrm{\Delta }_{022}].`$
Using Eqs. (31), (49), (135) and (136) one can show that the random variables $`X`$ and $`Y`$ are independent Gaussian random variables of zero mean and unit variance.

In deriving Eqs. (143) – (146) we assumed that the value of $`𝐯`$ which achieves the maximum in Eq. (124) corresponds a local maximum. However, if the right hand side of Eq. (143) is negative, the maximum will instead be achieved on the boundary of the parameter space at $`\widehat{\kappa }=0`$, since the variable $`\kappa `$ must be non-negative. Similarly, if the right hand side of Eq. (144) is less than 1, the maximum will be achieved at $`\widehat{\xi }=1`$, since $`1/\xi `$ must lie in the interval $`[1,\mathrm{})`$.

Substituting the results (143) – (146) \[together with the higher order corrections to those results which we have not shown\] into the expansion for the statistic $``$, and taking into account the various special cases discussed in the last paragraph, gives

$``$ $`=`$ $`[{\displaystyle \frac{1}{2}}(Y+\sqrt{6}q\stackrel{~}{\kappa }^2)^2ϵ^8\theta (Y+\sqrt{6}q\stackrel{~}{\kappa }^2)`$ (149)
$`+{\displaystyle \frac{1}{2}}(\stackrel{~}{\kappa }+ϵ^2X)^2ϵ^4\stackrel{~}{\kappa }^3ϵ^6+{\displaystyle \frac{7}{4}}\stackrel{~}{\kappa }^4ϵ^8`$
$`+\stackrel{~}{\kappa }Uϵ^7+\stackrel{~}{\kappa }Vϵ^8]\theta (\stackrel{~}{\kappa }+ϵ^2X)+O(ϵ^9).`$
Here $`\theta (x)`$ is the step function and

$`q`$ $`=`$ $`{\displaystyle \frac{1}{\stackrel{~}{\xi }}}1,`$ (150)
$`U`$ $`=`$ $`\mathrm{\Delta }_{101}+\mathrm{\Delta }_{110},`$ (151)
$`V`$ $`=`$ $`\mathrm{\Delta }_{200}{\displaystyle \frac{1}{2}}\stackrel{~}{\kappa }(\mathrm{\Delta }_{002}+\mathrm{\Delta }_{020})2\stackrel{~}{\kappa }\mathrm{\Delta }_{011}.`$ (152)
We note that the corresponding expression for the statistic $`(\mathrm{ln}\mathrm{\Lambda }_{\mathrm{ML}}^\mathrm{G})/N`$ \[which is equivalent to the cross-correlation statistic by Eq. (43)\] is given by Eq. (149) with the first term in the square brackets dropped.

Next we drop all the terms in the square bracket in Eq. (149) other than the first two terms. The reason is that these terms will give corrections that are smaller than the terms retained (both in expected value and in fluctuations) by a factor of

$$\stackrel{~}{\kappa }ϵ^2=\frac{\stackrel{~}{\rho }}{\sqrt{N}},$$

which will be small compared to unity for all cases we are interested in. This gives for the false dismissal probability the expression

$`P_{\text{FD}}`$ $`=`$ $`P(<_{})`$
$`=`$ $`{\displaystyle _{}}{\displaystyle \frac{dxdy}{2\pi }}\mathrm{exp}\left[{\displaystyle \frac{(xx_0)^2}{2}}{\displaystyle \frac{(yy_0)^2}{2}}\right],`$
where

$`x_0`$ $`=`$ $`\stackrel{~}{\kappa }/ϵ^2`$ (154)
$`y_0`$ $`=`$ $`\sqrt{6}q\stackrel{~}{\kappa }^2`$ (155)
$`r_0`$ $`=`$ $`\sqrt{2N_{}}.`$ (156)
Here the region $``$ in the $`x,y`$ plane is the union of the two regions

$`x`$ $``$ $`0`$
$`y`$ $``$ $`0`$
$`x^2+y^2`$ $``$ $`r_0^2`$ (157)
and

$`y`$ $``$ $`0`$
$`0x`$ $``$ $`r_0.`$ (158)
The integral over the region (158) is

$$P_{\mathrm{FD}}^{(1)}=𝒫(y_0)[𝒫(r_0x_0)𝒫(x_0)],$$
(159)
where

$$𝒫(x)1\frac{1}{2}\mathrm{erfc}(x/\sqrt{2})=_{\mathrm{}}^x𝑑t\frac{1}{\sqrt{2\pi }}\mathrm{exp}[t^2/2].$$
(160)
The integral over the region (157) can be written as

$`P_{\mathrm{FD}}^{(2)}`$ $`=`$ $`{\displaystyle \frac{1}{2\pi }}{\displaystyle _0^{\pi /2}}𝑑\theta {\displaystyle _0^{r_0}}𝑑rr`$
$`\times `$ $`\mathrm{exp}\left[{\displaystyle \frac{1}{2}}(r\mathrm{cos}\theta x_0)^2{\displaystyle \frac{1}{2}}(r\mathrm{sin}\theta y_0)^2\right].`$
The integrand in (LABEL:eq:pFD1) peaks at $`r\mathrm{cos}\theta =x_0`$, $`r\mathrm{sin}\theta =y_0`$. In order for $`P_{\mathrm{FD}}`$ to be small, its necessary that this peak occurs outside the domain of integration, at $`r>r_0`$. So we must have

$$x_0^2+y_0^2r_0^2.$$
(162)
The criterion $`x_0r_0`$ is, in order of magnitude, just the usual criterion for detectability with the cross-correlation statistic. The criterion $`y_0r_0`$ reduces to, in order of magnitude,

$$\xi \frac{\rho ^2}{\sqrt{N}}$$
(163)
which is what we claimed earlier to be the regime where the maximum likelihood statistic starts to work well, cf. Sec. IV.1.3 above.

Evaluating the integral (LABEL:eq:pFD1) using the Laplace approximation gives

$`P_{\mathrm{FD}}^{(2)}`$ $`=`$ $`{\displaystyle \frac{1}{r_0(\lambda 1)\sqrt{2\pi \lambda }}}\mathrm{exp}\left[{\displaystyle \frac{1}{2}}r_0^2(\lambda 1)^2\right]`$ (164)
$`\times \left[1+O\left({\displaystyle \frac{1}{r_0}}\right)\right],`$
where we define the variables $`\lambda `$ and $`\gamma `$ by

$$(x_0,y_0)=r_0\lambda (\mathrm{cos}\gamma ,\mathrm{sin}\gamma ).$$
(165)
However, the result (164) is not very accurate for small $`r_0`$. Alternatively we can integrate over $`r`$ in Eq. (LABEL:eq:pFD1) to obtain

$`P_{\mathrm{FD}}^{(2)}={\displaystyle _0^{\pi /2}}d\theta \{{\displaystyle \frac{1}{2\pi }}e^{\frac{r_{0}^{}{}_{}{}^{2}\left(1+\lambda ^2\right)}{2}}[e^{\frac{r_{0}^{}{}_{}{}^{2}}{2}}e^{r_{0}^{}{}_{}{}^{2}\lambda \mathrm{cos}(\gamma \theta )}]`$
$`+{\displaystyle \frac{r_0\lambda }{2\sqrt{2\pi }}}e^{\frac{r_{0}^{}{}_{}{}^{2}\lambda ^2}{4}\left[\mathrm{cos}(2\left\{\gamma \theta \right\})1\right]}\mathrm{cos}(\gamma \theta )[\mathrm{erf}\left({\displaystyle \frac{r_0\lambda \mathrm{cos}\{\gamma \theta \}}{\sqrt{2}}}\right)+\mathrm{erf}\left({\displaystyle \frac{r_0\left\{1\lambda \mathrm{cos}[\gamma \theta ]\right\}}{\sqrt{2}}}\right)]\},`$ (166)
where

$$\mathrm{erf}(x)=\frac{2}{\sqrt{\pi }}_0^x𝑑ye^{y^2}.$$
(167)
The integral (166) can be evaluated numerically. The false dismissal probability is then given by

$$P_{\mathrm{FD}}=P_{\mathrm{FD}}^{(1)}+P_{\mathrm{FD}}^{(2)},$$
(168)
with $`P_{\mathrm{FD}}^{(1)}`$ given by Eq. (159) and $`P_{\mathrm{FD}}^{(2)}`$ given by Eq. (166).

### C.2 False alarm probability

The false alarm probability is some function

$$P_{\mathrm{FA}}=P_{\mathrm{FA}}(_{},N)$$
(169)
of the threshold $`_{}`$ value of the detection statistic (124) and of the number of data points $`N`$. It does not depend on the signal parameters $`\stackrel{~}{\rho }`$ and $`\stackrel{~}{\xi }`$ because no signal is present. We would like to evaluate this quantity in the large $`N`$ limit.

We start by rewriting the statistic (124) in the form

$$=\underset{𝐯}{\mathrm{max}}\left\{\frac{1}{N}\underset{k=1}{\overset{N}{}}\mathrm{ln}A_k(0)+\frac{1}{N}\underset{k=1}{\overset{N}{}}\mathrm{ln}\left[1+\xi 𝒟_k(\alpha )\right]\right\},$$
(170)
where

$$𝒟_k(\alpha )=\frac{A_k(\alpha )}{A_k(0)}1.$$
(171)
Consider first the first term in Eq. (170). Using the definition (127) of $`A_k(\alpha )`$ and the definition (4) of $`\overline{\sigma }_1`$ and $`\overline{\sigma }_2`$ we can write this term as

$$\frac{1}{N}\underset{k=1}{\overset{N}{}}\mathrm{ln}A_k(0)=\frac{\mathrm{\Delta }\sigma _1^2}{\overline{\sigma }_1^2}\frac{\mathrm{\Delta }\sigma _2^2}{\overline{\sigma }_2^2}+O(\mathrm{\Delta }\sigma _1^3,\mathrm{\Delta }\sigma _2^3),$$
(172)
where $`\mathrm{\Delta }\sigma _1=\sigma _1\overline{\sigma }_1`$, $`\mathrm{\Delta }\sigma _2=\sigma _2\overline{\sigma }_2`$. Therefore the first term is maximized at $`\sigma _1=\overline{\sigma }_1`$, $`\sigma _2=\overline{\sigma }_2`$. Below we shall show that the second term in Eq. (170) is of order $`O(ϵ^2)`$, where in this subsection we define $`ϵ=1/\sqrt{N}`$. Therefore the values of $`\sigma _1`$ and $`\sigma _2`$ that achieve the maximum are

$`\widehat{\sigma }_1`$ $`=`$ $`\overline{\sigma }_1\left[1+O(ϵ^2)\right]`$
$`\widehat{\sigma }_2`$ $`=`$ $`\overline{\sigma }_2\left[1+O(ϵ^2)\right].`$ (173)
Moreover, in analyzing the second term it suffices to take $`\sigma _1=\overline{\sigma }_1`$, $`\sigma _2=\overline{\sigma }_2`$ in order to obtain the statistic to the leading $`O(ϵ^2)`$ order. Lastly, since we have assumed that $`\stackrel{~}{\sigma }_1=\stackrel{~}{\sigma }_2=1`$ and no signal is present, we have $`\overline{\sigma }_{1,2}=1+O(ϵ)`$. Hence, in analyzing the second term, it is sufficient to take $`\sigma _1=\sigma _2=1`$.

The statistic (170) therefore reduces to

$$=\underset{\alpha ,\xi }{\mathrm{max}}\frac{1}{N}\underset{k=1}{\overset{N}{}}\mathrm{ln}\left[1+\xi 𝒟_k(\alpha )\right]+O(ϵ),$$
(174)
where from Eqs. (127) and (171)

$$𝒟_k(\alpha )=\frac{1}{\sqrt{1+2\alpha }}\mathrm{exp}\left[\frac{w_k^2}{2+\frac{1}{\alpha }}\right]1.$$
(175)
Here $`w_k=(n_1^k+n_2^k)/\sqrt{2}`$, $`1kN`$, are independent Gaussian random variables of zero mean and unit variance.

It is straightforward to numerically compute the distribution of the statistic (174), by generating the Gaussian variables $`w_k`$ and numerically maximizing over $`\xi `$ and $`\alpha `$. The result is shown in Fig. 12. We find that at large $`N`$, the distribution of $`N`$ becomes independent of $`N`$, and is approximately given by

$$P(N>\xi )=\alpha _0e^{\beta _0\xi }$$
(176)
for $`\xi >0`$, where $`\alpha _00.42`$ and $`\beta _01.08`$. Therefore the false alarm probability is approximately given by

$$P_{\mathrm{FA}}=\alpha _0\mathrm{exp}\left[\beta _0N_{}\right].$$
(177)
Finally, we remark why it is plausible to expect the distribution of $`N`$ to be independent of $`N`$ in the large $`N`$ limit. The numerical maximizations over $`\xi `$ and $`\alpha `$ in Eq. (174) show that the maximum is nearly always achieved at $`\alpha 1`$ or $`\xi 1`$. In both these regimes, one can obtain some information about the $`N`$-dependence of the statistic.

Consider first the regime $`\xi 1`$. In this regime we can expand the expression (174) as a power series in $`\xi `$ to obtain

$$=\underset{\alpha ,\xi }{\mathrm{max}}\frac{1}{N}\underset{k=1}{\overset{N}{}}\left[\xi 𝒟_k(\alpha )\frac{1}{2}\xi ^2𝒟_k(\alpha )^2+O(\xi ^3)\right]+O(ϵ).$$
(178)
The generalized central limit theorem (reviewed in Appendix D) implies that

$$\frac{1}{N}\underset{k=1}{\overset{N}{}}𝒟_k(\alpha )=N^{\frac{1\gamma _1}{\gamma _1}}\left(\mathrm{ln}N\right)^{\delta _1}_N(\alpha ),$$
(179)
where for each fixed $`\alpha `$, the distribution of the random variable $`_N(\alpha )`$ becomes independent of $`N`$ in the large $`N`$ limit. Here

$$\gamma _1=\{\begin{array}{cc}2\hfill & 0<\alpha 1/2\hfill \\ 1+\frac{1}{2\alpha }\hfill & 1/2\alpha \hfill \end{array}$$
(180)
and

$$\delta _1=\{\begin{array}{cc}0\hfill & 0<\alpha 1/2\hfill \\ \frac{\alpha }{1+2\alpha }\hfill & 1/2\alpha .\hfill \end{array}$$
(181)
The limiting distribution is a Levy distribution with parameters $`p=1`$ and $`\gamma =\gamma _1`$. Similarly we have

$$\frac{1}{N}\underset{k=1}{\overset{N}{}}𝒟_k(\alpha )^2=N^{\frac{1\gamma _2}{\gamma _2}}\left(\mathrm{ln}N\right)^{\delta _2}𝒢_N(\alpha ),$$
(182)
where as $`N\mathrm{}`$ at each fixed $`\alpha `$ the distribution of the random variable $`𝒢_N(\alpha )`$ tends to a Levy distribution with parameters $`p=1`$ and $`\gamma =\gamma _2`$, with

$$\gamma _2=\{\begin{array}{cc}2\hfill & 0<\alpha 1/6\hfill \\ \frac{1+2\alpha }{4\alpha }\hfill & 1/6\alpha .\hfill \end{array}$$
(183)
and

$$\delta _2=\{\begin{array}{cc}0\hfill & 0<\alpha 1/6\hfill \\ \frac{2\alpha }{1+2\alpha }\hfill & 1/6\alpha .\hfill \end{array}$$
(184)
We now substitute the results (179) and (182) into the expression (178) for the statistic, and maximize analytically over the quadratic dependence on $`\xi `$. For $`\alpha 1/2`$, the value of $`\xi `$ which achieves the maximum goes to zero as $`N\mathrm{}`$, consistent with the assumption $`\xi 1`$, and the result is 131313For $`\alpha <1/2`$ this argument fails, which is why we must numerically verify that the distribution of $`N`$ is asymptotically independent of $`N`$.

$$N=\frac{1}{2}\underset{\alpha }{\mathrm{max}}\frac{_N(\alpha )^2}{𝒢_N(\alpha )}+O(ϵ).$$
(185)
In the regime $`\alpha 1`$, if we expand the expression (174) to quadratic order in $`\alpha `$, the result is an expression which is a linear function of $`1/\xi `$ at fixed $`\alpha \xi `$. Hence, when one maximizes over values of $`\xi `$ in the range $`0\xi 1`$, the maximum is always achieved either at $`\xi =0`$ or $`\xi =1`$. One can show that the maximum to this order is always achieved at $`\xi =1`$, and the resulting expression is

$$N=\frac{1}{4}𝒢^2+O(ϵ),$$
(186)
where

$$𝒢=\sqrt{N}\left[\frac{1}{N}\underset{k=1}{\overset{N}{}}w_k^21\right]$$
(187)
has a distribution that is independent of $`N`$ in the large $`N`$ limit.

## Appendix D Generalized central limit theorem

In this appendix we review the generalized central limit theorem that can be found on p. 574 of Ref. Feller . First we define a particular distribution function called the Levy distribution. It depends on 3 real parameters, a positive constant $`C`$, a parameter $`\gamma `$ in the range $`0<\gamma 2`$, and constant $`p`$ in the range $`0p1`$ 141414The parameter $`\gamma `$ is conventionally denoted by $`\alpha `$. We use $`\gamma `$ here to avoid confusion with the variable $`\alpha `$ defined in Eq. (7).. We say a random variable $`X`$ has a Levy distribution with parameters $`C`$, $`\gamma `$ and $`p`$ if the characteristic function of $`X`$ is given by

$`e^{i\zeta X}`$ $`=`$ $`\mathrm{exp}\{|\zeta |^\gamma {\displaystyle \frac{C\mathrm{\Gamma }(3\gamma )}{\gamma (\gamma 1)}}[\mathrm{cos}(\pi \gamma /2)`$ (188)
$`+i\mathrm{sgn}(\zeta )(pq)\mathrm{sin}(\pi \gamma /2)]\},`$
where $`q=1p`$. The corresponding probability distribution function is obtained by taking a Fourier transform and decays like $`x^{(1+\gamma )}`$ at large $`x`$ for $`\gamma <2`$ ($`\gamma =2`$ is the Gaussian case).

Consider now a random variable $`X`$ with probability distribution function $`f(x)`$ whose variance is infinite. Let

$$F(x)=_{\mathrm{}}^x𝑑yf(y)$$
(189)
be the cumulative distribution function and define

$$\mu (x)=_x^x𝑑yy^2f(y).$$
(190)
Suppose that the distribution satisfies the following conditions: (i) As $`x\mathrm{}`$ we have $`\mu (x)x^{2\gamma }L(x)`$, where $`0<\gamma 2`$, and $`L(x)`$ varies slowly in the sense that $`L(tx)/L(t)1`$ as $`t\mathrm{}`$ for all $`x>0`$. (ii) We have

$$\frac{1F(x)}{F(x)+1F(x)}p\frac{F(x)}{F(x)+1F(x)}q$$
(191)
as $`x\mathrm{}`$, where $`0p1`$, $`0q1`$ and $`p+q=1`$. (iii) For $`1<\gamma 2`$, we assume that the expected value $`𝑑xxf(x)`$ vanishes; this can be enforced by making a transformation of the form $`XX+\mathrm{constant}`$.

We define the sequence of random variables

$$S_N=\frac{1}{a_N}\underset{i=1}{\overset{N}{}}X_i,$$
(192)
where the $`X_i`$ are independent, identically distributed random variables with distribution function $`f`$, and the constants $`a_N`$ are chosen to satisfy

$$\frac{N\mu (a_N)}{a_N^2}C$$
(193)
as $`N\mathrm{}`$, where $`C`$ is a positive constant. Then, the distribution functions of the random variables $`S_N`$ converge to a Levy distribution with parameters $`C`$, $`\gamma `$ and $`p`$ as $`N\mathrm{}`$.