# Signal-to-noise-ratio and SNR-max detection statistics in template bank
searches for exotic physics transients with networks of quantum sensors

T. Daykin<EMAIL_ADDRESS>Department of Physics, University of Nevada, Reno,
89557, USA  C. Ellis  Department of Physics, University of Nevada, Reno,
89557, USA  A. Derevianko<EMAIL_ADDRESS>Department of Physics, University of
Nevada, Reno, 89557, USA

###### Abstract

Quantum sensor networks such as the existing networks of atomic clocks and
magnetometers offer intriguing capabilities in searches for transient signals
such as dark matter or fields sourced by powerful astrophysical events. The
common matched-filter technique relies on signal-to-noise (SNR) detection
statistic to probe such transients. For macroscopic dark matter objects, the
network would register a sweeping signal as the object propagates through at
galactic velocities. A potential event is registered when the SNR from a
specific template exceeds a threshold set by a desired false positive rate.
Generically, to span the continuous parameter space for the network-exotic-
physics encounter, one has to deal with multiple templates. In such template
bank searches, the natural generalization of the SNR statistic is the SNR-max
statistic, defined as the maximum of the absolute values of SNRs determined
from individual template matching. While individual SNR realizations are
Gaussian distributed, SNR-max probability distribution is non-Gaussian.
Moreover, as the individual template-bank SNRs are computed using the same
network data streams, SNRs become correlated between the templates. Cross-
template correlations have a sizable effect on the SNR-max probability and
cumulative distributions, and on the threshold SNR-max values. Computing
threshold SNR-max values for large template banks is computationally
prohibitive and we develop analytic approaches to computing properties of SNR-
max statistic. This is done for cases when the template bank is nearly
orthogonal (small cross-template correlations) and for banks with cross-
template correlation coefficients distribution “squeezed” about the most
probable cross-template correlation value. Since the cross-template
correlation coefficients quantify the similarity of templates, increasing
correlations tend to decrease SNR-max thresholds for specific values of false
positive rates. Increasing the number of templates in the bank increases the
SNR-max thresholds. Our derivations are carried out for networks that may
exhibit colored noise and cross-node correlations. Specific applications are
illustrated with a dark matter search with atomic clocks onboard satellites of
Global Positioning System (GPS) and with a “toy” planar network with cyclic
rotational symmetry.

##  I Introduction

As precision quantum sensors [1] become ubiquitous, it is natural to combine
them into networks. In this paper, we focus on time-domain searches for exotic
physics signals with networks of quantum sensors. Time-domain analysis is a
natural fit for transient non-oscillating signals. We refer the reader to
frequency-domain analysis [2] relevant to searches for oscillating dark matter
(DM) signals with networks or to combined time-frequency domain analyses [3,
4, 5] for detecting oscillating transients such as bursts of gravitational
waves [3, 4] or low-mass exotic field (ELFs) emitted by powerful astrophysical
events [5].

The spatial extent of networks can range from microscopic (atomic ensemble
scales) to continental or trans-continental scales. Specific examples include
a constellation of microwave atomic clocks on-board GPS satellites [6],
network of gravitational wave detectors [7, 8], trans-European network of
high-performance optical clocks [9], and the global network of atomic
magnetometers (GNOME) [10]. The nodes of the networks can be connected by
physical links (e.g., optical fiber or microwave/laser links). Alternatively,
synchronized cross-node comparison can be accomplished by a simple time-
stamping of sensor data relying on world-wide broadcast of GPS time. The
network can be composed of identical sensors (homogeneous network) or
different sensors (heterogeneous network). The position of the nodes can be
stationary or evolving in time (like in the case of orbiting GPS satellites).
Finally, individual sensors can use entangled or spin-squeezed ensembles that
operate below the quantum projection noise limit, see, e.g. [11, 12, 13]. The
formalism developed in this paper applies to all these cases and our results
hold for any network with stationary white or colored Gaussian noise and with
and without cross-node correlations.

This paper grew from our practical experience with the GPS.DM search [6] for
DM-induced transients using GPS network of atomic clocks. Multiple galactic
scale astrophysical observations indicate that DM comprises $85\%$ of all
matter in the universe, with only $15\%$ left to ordinary matter [14]. These
galactic scale observations only characterize the gravitational interaction of
ordinary matter and dark matter. Little is known of the microscopic
composition of DM or of its non-gravitational interactions with standard model
particles and fields. The primary target for current DM searches are weakly
interacting massive particles (WIMPs). With WIMP searches so far failing to
provide convincing evidence, alternative DM candidates are being considered
[14]. Self-interacting ultralight fields is one such alternative, where the DM
candidate can take a form of a macroscopic object that is coherent over large
scales [15, 16, 17, 18]. These macroscopic DM objects is the target of the
GPS.DM search.

In the context of DM models, ultralight fields are characterized by
macroscopic mode occupation numbers and can thus be described as classical
fields. The effect of exotic physics on quantum sensors depends on specific
interactions [18] (portals) between exotic fields and standard model particles
and fields. Particle detectors rely on energy-deposition. By contrast,
ultralight fields, through their portal interactions, gentle modulate atomic
(qubit) energy levels. Detecting such gentle modulations is the domain of
precision quantum sensors. For example, for atomic clocks [18, 19], atom
interferometers [20, 21, 22, 23, 24, 25] , and optical cavities [26, 22, 27] ,
DM objects can drive variation of fundamental constants of nature. The change
in fundamental constants in turn affects atomic or cavity frequencies causing
an atomic clock to speed up or slow down when the DM object overlaps with a
sensor. Alternatively, gradients of objects composed of axion-like fields act
as fictitious magnetic fields and can be detected with another class of
quantum sensors, magnetometers [28, 29] In both examples, as the DM object
sweeps through, the induced frequency shifts or fictitious magnetic fields
would appear as a transient perturbation sweeping through the network. A
geographically-distributed network can resolve the directionality and speed of
the transient signals. Further, this provides vetoing mechanism as the network
sweep speed and directionality must be consistent with standard halo model
(SHM) priors. Currently, there are several on-going DM searches that follow
these ideas [6, 30, 10, 23, 29, 25, 27, 24]

Our GPS.DM collaboration [6, 31], focuses on searching for DM transient
signatures by utilizing the network of atomic clocks aboard the GPS
satellites. The GPS data has an inherent advantage over building a new network
of quantum sensors as there is over two decades of publicly available archival
data [32], ready for data mining. DM signatures would consist of a correlated
propagation of atomic clock frequency perturbations through the GPS network at
galactic scale velocities ($\sim 300\,\mathrm{km/s}$). Previously, our GPS.DM
group has performed analysis of the archival GPS data in search for 2D DM
walls [6]. Although no DM signatures were found, prior astrophysical limits on
certain DM couplings to atoms were improved by several orders of magnitude.
Ref. [33] developed the application of the more sophisticated Bayesian search
techniques which may extend the DM discovery reach further by several orders
of magnitude in both the sensitivity and size/geometry of the DM objects. The
main focus of Ref. [31] was on the performance of an alternative frequentist
approach, the matched filter technique (MFT). That work obtained several
analytical results for an idealized network of white-noise sensors, including
cross-node correlation. Here we extend the results of Ref. [31] to template
banks.

The MFT is an ubiquitous technique, utilized for example by the Laser
Interferometer Gravitational Wave Observatory (LIGO) in their gravitational
wave detection [7, 8]. Also, MFT has applications in astrophysics [34, 35, 36,
37, 34], geophysics [38], and searches for exotic physics [39, 40, 41, 42,
43]. In applications, a large template bank determined by prior information on
the expected signal, is used to match data streams. Quantitatively, the degree
of the match is characterized by a signal-to-noise ratio (SNR). In Ref. [31],
the SNR for individual templates was shown to be a Gaussian normal random
variable. Further maximization of the template specific SNR variables
(ignoring dependence on sign), was performed over the template bank. In this
process, a Gaussian distribution was considered instead of the distribution
for the maximum of the absolute valued SNRs (SNR-max). As we discuss in this
paper, the SNR-max distribution is no longer a Gaussian distribution, but a
skewed (non-Gaussian) distribution. Moreover, there exist correlations in the
SNR statistics computed for different templates in the bank, which alter the
uncorrelated SNR-max distribution. Thresholds for specific false positive
rates determined in Ref. [31] overestimate the thresholds due to the effect of
correlated templates in a template bank. One of the goals of this paper is to
determine the thresholds correctly.

The structure of this paper is as follows. Sec. II reviews the SNR and SNR-max
statistics in the context of detecting transients and introduces notation used
in the rest of the paper. In Sec. III, we discuss various aspects of
generation of template banks in the network searches for transients. This is
illustrated with a pedagogical toy problem of a highly-symmetric planar
network and with a more practical GPS.DM search. Strategies for generating
optimal template banks are discussed in Sec. III.3. Sec. IV focuses on
correlations between SNRs for different templates and introduces template bank
covariance matrix. Computations of the template bank covariance matrix are
illustrated with the toy problem and with the GPS.DM search. Sec. V derives
the SNR-max probability distribution (PDF) and cumulative distribution (CDF)
distribution for several practical yet analytically treatable cases: for a
fully-orthogonal template bank, a nearly orthogonal template bank, and for
template bank ”squeezed” about the mean value of covariance matrix elements.
Sec. VI discusses the effect on thresholds when the template bank is fully-
orthogonal, a special case of $M=2$ template bank, a nearly-orthogonal bank,
and a “squeezed” bank. Lastly, Sec. VII draws the conclusions and addresses
the utility of networks with entanglement shared by geographically-separated
nodes. The paper contains several appendices where we present detailed
derivations for some of the paper results.

##  II SNR and SNR-max statistics

In the previous work [31] by our group, we explored the properties of the SNR
statistic in the context of networks of quantum sensors. In this section we
set up the problem, review relevant results from Ref. [31], and then introduce
the SNR-max detection statistic.

Throughout this paper, we assume that the network data streams are stationary
and Gaussian-distributed. The sensor noise can be colored and the noise can be
correlated across the network. For an atomic clock network, the cross-network
correlation comes from referencing every clock phase to that of “reference”
clock common to all clocks [6]. Notice our assumption of the noise
stationarity; for clocks this implies working with frequency data instead of
phase (time or time bias as recorded by the clock) data.

Network data streams $\mathbf{d}$ are comprised of sensor measurements
$d_{i}^{a}$, where the subscript $i$ refers to the time stamp (time $t_{i}$)
and superscript $a$ identifies a specific sensor with the total number of
sensors $N_{\rm S}$. The time grid is assumed to be uniform with a step size
$\Delta_{t}$. We will refer to a time interval $[t_{i},t_{i+1})$ as an epoch
$i$.

In the most basic search for transient signals, network data streams
$\mathbf{d}$ are partitioned into data windows of fixed length $J$ and the
data inside each window are matched to a bank of $M$ templates. We enumerate
windows with an index $w$ and refer to the portion of the data stream within a
specific data window as $\mathbf{d}^{w}$. The template bank contains multiple
signal templates, $\mathbf{s}_{k}$, $k=\overline{1,M}$. The templates must
reflect underlying physics of the sought signals. The length of each template
in time domain is equal to that of the data window $J$. For a network, each
template is represented by an $N_{\rm S}\times J$ matrix.

As a concrete example of physics-motivated generation of a signal template
bank, in Sec. III we discuss the bank used in the GPS.DM search. Briefly, in
the GPS.DM search we probe data streams for sweeps by walls of DM bubbles,
with the assumption that the radius of the DM bubble is much larger than the
network spatial extent. Then the network is swept by planar walls. These
hypothesized walls are to sweep through the network with velocities within a
range prescribed by the SHM [44]. Each template is distinguished by the
specific values of a DM wall incoming velocity and orientation. We typically
use $M=1024$ templates in the bank to ensure a sufficiently dense coverage of
the sweep parameters. There are nominally $N_{\rm S}=32$ atomic clocks on-
board GPS satellites populating orbits on a spherical shell of diameter $D\sim
5\times 10^{4}\,\mathrm{km}$. The atomic clock time readings (biases) are
sampled every $\Delta_{t}=30\,\mathrm{s}$ and the length of data window $J=61$
is determined by the ratio of the maximum network sweep time to the sampling
time, $(D/v_{\mathrm{min}})/\Delta_{t}$, where the minimum speed
$v_{\mathrm{min}}\approx 25\,\mathrm{km/s}$. The minimum value of the speed is
determined by the requirement that the spatial network configuration does not
change substantially over the duration of the sweep, since the orbital
velocity of the satellites is $\sim 5\,\mathrm{km/s}$.

The sought physics signal can be parameterized as $h\mathbf{s}_{k}$, with $h$
being the amplitude or the strength of the signal. $h$ is to be estimated from
the match to the template. Typically, $h$ encodes the coupling of the sensor
to the exotic physics. Formally, the signal strength estimator for a specific
window $w$ and a template $k$ reads [31]

$\hat{h}^{w}_{k}=\frac{\left(\mathbf{d}^{w}\right)^{T}\mathbf{E}^{-1}\mathbf{s}_{k}}{\mathbf{s}_{k}^{T}\mathbf{E}^{-1}\mathbf{s}_{k}}\,,$
(1)

where $\mathbf{E}^{-1}$ is the inverse of the network noise covariance matrix
$\mathbf{E}$, discussed in Sec. II.0.1.

Further, the standard deviation of the estimated signal amplitude, Eq. (1), is
[31]

$\sigma_{h,k}=\frac{1}{\sqrt{\mathbf{s}_{k}^{T}\mathbf{E}^{-1}\mathbf{s}_{k}}}\,.$
(2)

This uncertainty does not depend on the measurement data but it does depend on
a specific template. Finally, the SNR, $\rho^{w}_{k}$, in window $w$ for a
template $k$, is simply the signal strength estimator (1) divided by the
standard deviation of the signal strength (2),

$\rho^{w}_{k}\equiv\frac{\hat{h}^{w}_{k}}{\sigma_{h,k}}=\frac{\mathbf{s}_{k}^{T}\mathbf{E}^{-1}\mathbf{d}^{w}}{\sqrt{\mathbf{s}_{k}^{T}\mathbf{E}^{-1}\mathbf{s}_{k}}}.$
(3)

As shown in Ref. [31], in the absence of a signal (i.e. when
$\bm{d}^{w}\equiv\bm{n}^{w}$, noise of sensors), the SNR statistic is
Gaussian-distributed with zero mean and variance of 1,

$\rho^{w}_{k}\sim\mathcal{N}(0,1).$ (4)

Here $\mathcal{N}(\mu,\sigma)$ is the normal distribution with mean $\mu$ and
variance $\sigma^{2}$. This can be seen by examining the SNR definition (3).
The denominator, $\sqrt{\mathbf{s}_{k}^{T}\mathbf{E}^{-1}\mathbf{s}_{k}}$, is
a constant as the template and the network noise covariance matrix are fixed,
while the numerator $\mathbf{s}_{k}^{T}\mathbf{E}^{-1}\mathbf{d}^{w}$ is a
linear combination of Gaussian random variables. Thereby, the SNR is also a
Gaussian random variable with mean
$\left\langle\rho^{w}\right\rangle\propto\langle\mathbf{n}^{w}\rangle=0$. The
SNR variance is
$\left\langle(\rho^{w}_{k})^{2}\right\rangle=\mathrm{Var}((\mathbf{n}^{w})^{T}\mathbf{E}^{-1}\mathbf{s_{k}})/\mathbf{s}_{k}^{T}\mathbf{E}^{-1}\mathbf{s}_{k}$.
Ref. [31] (Appendix B) shows that
$\mathrm{Var}((\mathbf{n}^{w})^{T}\mathbf{E}^{-1}\mathbf{s}_{k})=(\mathbf{s}_{k})^{T}\mathbf{E}^{-1}\mathbf{s}_{k}$,
which results in an SNR variance
$\left\langle(\rho^{w}_{k})^{2}\right\rangle=1$.

Eq. (4) means that given multiple realizations of the signal-free data
streams, the values of SNR $\rho_{k}^{w}$ for a fixed template
$\mathbf{s}_{k}$ are distributed as

$p_{\mathrm{SNR}}(\rho_{k}^{w})=\frac{1}{\sqrt{2\pi}}e^{-(\rho_{k}^{w})^{2}/2}\,.$
(5)

The above definitions and properties of the SNR statistic hold for a general
case of colored noise and multi-sensor networks. If the subsequent windows
overlap, then the SNR values for the two windows are correlated as they are
based on partially overlapping data. Even then the SNR probability
distribution in individual windows still follows Eq. (5). This can be proven
by marginalizing joint probability distribution (multi-variate normal
distribution) over SNRs in all other overlapping windows.

So far we reviewed the properties of the SNR statistic for a single template.
Practical applications involve banks containing multiple templates. For a
given window $w$, all $M$ individual template SNRs are compared and the
maximum of the absolute values of the individual template SNRs is used as a
detection statistic,

$z^{w}\equiv\max(|\rho_{1}^{w}|,..,|\rho_{M}^{w}|)\,.$ (6)

Following the literature [45, 46, 47], we refer to this quantity as the SNR-
max statistic. To simplify notation, from this point forward we drop the
window superscript $w$, so that the above definition reads
$z\equiv\max(|\rho_{1}|,..,|\rho_{M}|)$.

Matching the same data stream to multiple templates induces cross-template
correlations in the SNR statistics, $\langle\rho_{i}\rho_{j}\rangle\neq 0$. In
addition, comparing SNRs and taking their maximum per Eq. (6) suppresses
smaller values of SNR. This suppression and cross-template correlations
substantially modify the single template SNR distribution (5). It is the focus
of our paper to quantitatively characterize the SNR-max probability and
cumulative distributions. Specifically, we are interested in a practically
important threshold value for false positives. We remind the reader that the
detection can be claimed when the value of chosen statistic exceeds the
detection threshold [48]. The threshold values for SNR-max are required to
ensure (at a given confidence level) that the observed SNR-max values
exceeding this threshold are due to probed physics and not due to randomness.
Our previous analysis [31] did not take into account that the detection
statistics may be correlated by templates or that our statistic is no longer
Gaussian-distributed. Here we present a more sophisticated analysis to improve
upon the treatment of Ref. [31]. Before proceeding to our analysis, below we
review relevant properties of the network noise covariance matrix
$\mathbf{E}$, entering the SNR definition (3).

####  II.0.1 Network noise covariance matrix $\mathbf{E}$

Network noise covariance matrix $\mathbf{E}$ is defined by matrix elements

$E_{(al)(bm)}=\left\langle n_{l}^{a}n_{m}^{b}\right\rangle\,.$ (7)

This definition involves averaging over realizations of signal-free data
$n_{l}^{a}$, i.e. intrinsic noise of the sensors. Without loss of generality,
we assume that the noise has zero mean, $\langle n_{l}^{a}\rangle=0$. In our
index convention, the letters at the beginning of alphabet $a,b,\ldots$
enumerate the sensors (range $\overline{1,N_{\rm S}}$) , while letters in the
middle of alphabet $l,m,\ldots$ — the epochs (index range $\overline{1,J}$ in
a data window). We also introduced a compound index notation $(al)$, which
organizes sensor-epoch data points into a vector, i.e. we unroll a $N_{\rm
S}\times J$ matrix into a 1D array. For example, for the first window,
$\mathbf{d}^{w=1}=(d_{1}^{1},d_{2}^{1},\ldots
d_{J}^{1},d_{1}^{2},d_{2}^{2},\ldots d_{J}^{2},\ldots d_{1}^{N_{\rm
S}},d_{2}^{N_{\rm S}},\ldots d_{J}^{N_{\rm S}})^{T}$. We will use Greek
letters for such compound indexes, e.g. Eq. (7) can be rewritten as
$E_{\alpha\beta}=\left\langle n_{\alpha}n_{\beta}\right\rangle$. We remind the
reader that the covariance matrix is symmetric and positive semi-definite.

As an illustration, for a network of white-noise sensors without a common
reference sensor, the noise covariance matrix is diagonal,

$E_{(al)(bm)}=\sigma^{2}_{a}\delta_{ab}\delta_{lm}\,,$ (8)

with $\sigma_{a}^{2}$ being the noise variance for sensor $a$. The matrix
elements of its inverse are

$\left(E^{-1}\right)_{(al)(bm)}=\frac{1}{\sigma_{a}^{2}}\delta_{lm}\delta_{ab}\,,$
(9)

i.e., $\mathbf{E}^{-1}$ is a diagonal matrix.

The same network but referenced to a common sensor (variance $\sigma_{R}^{2}$)
has the noise covariance matrix

$E_{(al)(bm)}=(\sigma^{2}_{a}\delta_{ab}+\sigma_{R}^{2})\delta_{lm}\,.$ (10)

Compared to the no-reference sensor case, Eq. (8), the reference sensor
introduces off-diagonal matrix elements. When all the sensors, except for the
reference sensor, are assumed to have the same variance $\sigma^{2}$, the
inverse of the covariance matrix $\mathbf{E}^{-1}$ can be found analytically
[31]

$\left(E^{-1}\right)_{(al)(bm)}=\frac{1}{\sigma^{2}}\delta_{lm}\left(\delta_{ab}-\frac{1}{N_{\rm
S}}\frac{\xi}{1+\xi}\right)\,,$ (11)

where the parameter

$\xi=N_{\rm S}\sigma_{R}^{2}/\sigma^{2}$ (12)

has the meaning of the square of the ratio of network-averaged sensitivity
($\propto\sigma/\sqrt{N_{\rm S}}$) to the reference sensor sensitivity
($\propto\sigma_{R}$). We will refer to $\xi$ as the relative reference-sensor
sensitivity parameter.

It is instructive to consider a limiting case of a single sensor ($N_{\rm
S}=1$) data stream with white noise of zero mean and variance $\sigma^{2}$. In
this case, the inverse of the covariance matrix is simply
$\mathbf{E}^{-1}=\sigma^{-2}\mathbf{I},$ where $\mathbf{I}$ is the $J\times J$
identity matrix. Further, Eqs. (1), (2), and (3) reduce to
$\hat{h}^{w}_{k}=\left(\mathbf{d}^{w}\cdot\mathbf{s}_{k}\right)/{\left(\mathbf{s}_{k}\cdot\mathbf{s}_{k}\right)}$,
$\sigma_{h}=\sigma/{\sqrt{\left(\mathbf{s}_{k}\cdot\mathbf{s}_{k}\right)}}$,
and
$\rho^{w}_{k}=\left(\mathbf{d}^{w}\cdot\mathbf{s}_{k}\right)/{(\sigma\sqrt{\left(\mathbf{s}_{k}\cdot\mathbf{s}_{k}\right)})},$
where we used the conventional definition of the dot product of two vectors.

##  III Generation of template banks in the searches for transients

Having defined the SNR and SNR-max statistics, now we formalize the procedure
for template bank generation. Template banks must reflect signals prescribed
by the underlying physics to be probed. While the signal templates are
application-specific, some generic requirements and challenges can be
illustrated with specific examples. First, in Sec. III.1 we consider a “toy”
2D problem of a network search for transient signals induced by sweeps of
randomly oriented lines through a circular network of high-degree symmetry.
This toy problem introduces a number of relevant issues specific for networks,
such as the angular and speed resolutions. After understanding the simplified
problem, in Sec. III.2, we discuss a more practical case of the GPS.DM search
[6, 31, 49, 50].

###  III.1 Toy problem: Circular topology network with $C_{N_{\rm S}}$
rotational symmetry

As an illustration of generating template banks and their properties, we
consider a planar network with a circular topology where identical point-like
sensors are distributed evenly on a circle of radius $R$, see Fig. 1.

Figure 1: (Color online) (Left panel) A planar network of $N_{\rm S}=10$
sensors with $C_{10}$ symmetry. (Right panel) Graphical representation of a
bank of $M=5$ templates. Arrows indicate direction of the line sweep
$\hat{u}_{i}$. The arrows $\hat{u}_{i}$ are normal to the sweeping line. For
example, template 1 corresponds to a vertically-oriented line sweeping from
the left.

Sensors’ polar angles are given by

$\phi_{a}=\frac{2\pi}{N_{\rm S}}(a-1),$ (13)

so that the sensor coordinates are

$\bm{r}_{a}=R(\cos\phi_{a},\sin\phi_{a})^{T}=R\hat{n}_{a},$ (14)

where we used the coordinate system with the origin at the circle center and
$\hat{x}$ pointing towards sensor 1. We also introduced a unit vector
$\hat{n}_{a}$ pointing towards sensor $a$. For example, for the network of
$N_{s}=10$ sensors, shown in Fig. 1, $\phi_{1}=0$ and $\phi_{2}=\pi/5$. This
arrangement of sensors has the $C_{N_{\rm S}}$ cyclic group symmetry, as
rotating the network by multiples of the angle $2\pi/N_{\rm S}$ about the
circle center maps the network onto itself. In this toy problem, we do not
include effects of the reference sensor.

We note that this 2D network topology is of relevance to the configuration of
satellites in the GPS constellation. The satellites in the GPS network are
arranged into six equally-spaced orbital planes with an inclination of
$55^{\circ}$, and each nearly-circular orbit is populated with 4-5 satellites.
Because the satellites are not evenly spaced apart within each orbit however,
the $C_{N_{\rm S}}$ symmetry is broken. Still, understanding our toy network
provides a number of useful insights applicable to general network topology.

With the toy network topology specified, now we generate a template bank for
sweeps of this 2D network by “lines”. These lines can be thought as 2D cross-
sections of thin DM walls. The lines sweep through the network, leaving a
geographically distributed imprint on the sensor signals. The sensor is
considered to be affected in a given epoch $l$, when the line crosses the
sensor in the time interval $[t_{l},t_{l}+\Delta_{t})$.

In our 2D example, the velocities $\bm{v}$ of the lines lie in the network
plane. However, the component of the velocity parallel to the line leads to an
unobservable translation of the line along itself. Different templates are
then uniquely specified by the direction and value of velocity component
$\bm{v}_{\perp}=v_{\perp}\hat{u}$, where $\hat{u}$ is the unit vector
perpendicular to the line and pointing in the direction of the sweep. We will
use polar angle $\theta$ in the network reference frame to characterize
incidence direction $\hat{u}$. For example, template bank in Fig. 1 samples
$M=5$ directions with $\theta_{k}=(2\pi/M)\,(k-1)$, $k=\overline{1,M}$. For
simplicity, we assume that all the directions $\hat{u}$ are equally probable.
The bank should provide an adequate sampling of both incidence directions
$\hat{u}$ and speeds $v_{\perp}$.

Templates can be generated based on the following argument. Suppose $t_{O}$
and $t_{a}$ are the moments of time when the line reaches the network center
and the sensor $a$, respectively. The difference between these two times is

$t_{a}-t_{O}=(\hat{u}\cdot\hat{n}_{a})\frac{R}{v_{\perp}}.$ (15)

Since our time sampling is finite, the sensor signal time series is simply
$s_{l}^{a}=\delta_{l,l_{a}}$, where the epoch $l_{a}$ is given by

$l_{a}=\left\lfloor(\hat{u}\cdot\hat{n}_{a})\frac{R}{v_{\perp}\Delta_{t}}\right\rfloor+l_{O}.$
(16)

Here $\lfloor x\rfloor$ is the floor function and $l_{O}=(J+1)/2$ is the
center of the data window.

Figure 2: (Color online) Construction of the template corresponding to a
vertically-oriented line sweeping from the left (polar angle $\theta=0$).
Sensor 6 is affected first and sensor 1 — last. Here we used
$R/(v_{\perp}\Delta_{t})=5.01$ and the number of epochs $J=15$. The center of
the network is crossed in the $l_{O}=7$ epoch.

Fig. 2 presents a construction of a template for a vertically-oriented line
sweeping our network from the left. This is template 1 from the template bank
shown in Fig. 1. The left panel shows slicing of the network by the sweeping
line with the spatial slice width $v_{\perp}\Delta_{t}$. Slice boundaries
correspond to the sweeping line snapshot at the beginning and the end of a
given epoch, e.g. the left-most slice is epoch 1.

Since $-1\leq\hat{u}\cdot\hat{n}_{a}\leq 1$ in Eq. (16), the constructed
template spreads over

$\Delta l=2R/(v_{\perp}\Delta_{t})$ (17)

epochs. Clearly, the number of epochs in the window $J$ should be greater than
$\Delta l$, otherwise a useful portion of the signal is cut out. For fixed
network geometry and sampling time $\Delta_{t}$, $\Delta l$ increases
indefinitely as $v_{\perp}$ approaches small values. Practical considerations,
however, introduce a cutoff at small $v_{\perp}$. For example, GPS satellites
positions change with time due to their orbital motion; this limits
$v_{\perp}$ from below, see Sec. II.

If $v_{\perp}>(v_{\perp})_{\mathrm{max}}=2R/\Delta_{t}$, all the sensors are
affected in the same epoch, so the template width $\Delta l=1$.
$(v_{\perp})_{\mathrm{max}}$ is the speed at which the line sweeps through the
entire spatial extent (diameter = $2R$) over the duration of the sampling
interval $\Delta_{t}$. In this case, the network losses directional and
velocity resolutions. That is for any incident line orientation and any
$v_{\perp}>(v_{\perp})_{\mathrm{max}}$, all the templates are identical,
collapsing into a vertical single-epoch column at $l_{a}=l_{O}$ in Fig. 2.
Underlying physics may constraint $(v_{\perp})_{\mathrm{max}}$ to lower
values, see, e.g., discussion of the GPS.DM search in Sec. III.2, where
$(v_{\perp})_{\mathrm{max}}$ is determined by the DM velocity distribution
cutoff.

As $v_{\perp}$ decreases, the template width $\Delta l$ increases and
different sensors begin to be perturbed at different epochs. On average, we
expect ${N_{\rm S}}/{\Delta l}$ sensors to be affected at the same epoch. We
refer to such groups of sensors as being degenerate and introduce the average
degeneracy factor

$\bar{g}=N_{\rm S}\frac{v_{\perp}\Delta_{t}}{2R}.$ (18)

This degeneracy factor is a rule of thumb. For example, for the bank and
network parameters in Fig. 2, from Eq. (18) we expect $\bar{g}=1$, while we
observe mostly two-fold degeneracies in the resulting template shown in Fig.
2. This is a consequence of the line propagation direction $\hat{u}$ being
aligned with one of the network axes of symmetry. Generically, in the limit of
small degeneracy factors, $\bar{g}\leq 1$, we expect all the sensors to be
perturbed at different epochs.

To generate the template bank we need to sample continuous parameter space of
directions $\hat{u}$ and normal speeds $v_{\perp}$. This leads to the question
of angular and velocity resolutions of our network.

We start with the angular resolution. Consider two sweeping lines which have
the same $v_{\perp}$ but differ by their angles of incidence: $\theta$ and
$\theta^{\prime}=\theta+\Delta_{\theta}$. Because both the number of sensors
and time sampling are finite, for sufficiently small relative tilt angles
$\Delta_{\theta}$, sweeps by both lines produce the same template. As we
increase the tilt, there is some critical value of $\Delta_{\theta}$ when the
two sweeps start mapping onto two distinct templates. Then the two incident
directions can be resolved in the network search. The critical value of the
tilt determines the network angular resolution.

As an illustration, consider the template construction in Fig. 2 for a line
sweeping from the left, $\theta=0$. Templates for tilted sweeps can be
generated by rotating the network sensors counter-clockwise by the angle
$\Delta_{\theta}$ about the circle center. Since the sensors 5 and 7 are near
the edge of epoch 3 slice (and sensors 2 and 10 are near the edge of epoch
10), even small tilts $\Delta_{\theta}>0$ move sensor 7 into epoch 4 and
sensor 2 into epoch 11. However, sensors being at the slice edge is an
exceptional situation. Generically, the sensors are more likely to be found in
the middle of epoch slices. Then the angular resolution is determined by the
maximally separated sensors within the same slice, i.e. sensors swept within
the central epoch. This leads to the worst-case-scenario angular resolution
estimate,

$\Delta_{\theta}\approx\sin^{-1}\left(\frac{v_{\perp}\Delta_{t}}{2R}\right)\,.$
(19)

Recalling the relation (17) for the template width $\Delta l$, we observe that
$\sin\Delta_{\theta}=1/\Delta l$. When the template collapses into a single
column due to high values of normal velocity, $\Delta l=1$, and the angular
resolution is lost ($\Delta_{\theta}\rightarrow\pi/2$) in agreement with our
previous qualitative statements. In the opposite limit of $\Delta l\gg 1$,
$\Delta_{\theta}\approx 1/\Delta l$. For the template of Fig. 2,
$\Delta_{\theta}\approx 1/12\,\mathrm{radians}$. Full, resolution-limited,
coverage of the incidence angles by the bank then requires
$M=2\pi/\Delta_{\theta}\approx 75$ templates. Having degenerate (affected in
the same epoch, but spatially separated) sensors improves the angular
resolution. It is worth emphasizing that the derived angular resolution is a
pure consequence of the line-network encounter geometry and it does not
account for the intrinsic sensor noise, see Refs. [33, 31], which will worsen
the angular resolution (19).

Angular resolution of a geographically-distributed network has been considered
in several papers. For example, in Ref. [5], there is a brief discussion on
the angular resolution for bursts of exotic fields propagating through a
network of quantum sensors. Their discussion suggests that the angular
resolution is roughly the ratio of the temporal resolution $\Delta_{t}$ to the
field burst propagation time through the sensor network ($2R/v_{\perp}$), in
agreement with our formula (19) in the limit of high sampling rates,
$\Delta_{t}\ll 2R/v_{\perp}$. We also refer the reader to papers [51, 3]
discussing “Rayleigh criterion” angular resolution for a network of
interferometers.

Now we find the resolution in normal velocities $v_{\perp}$. We apply the same
logic of finding the minimum increment $\Delta_{v}$ so that the sweeps with
two velocities $v_{\perp}$ and $v_{\perp}+\Delta_{v}$ map into two distinct
templates. Since the number of sensors and time sampling is finite,
vanishingly small variations in the normal velocity would produce identical
templates. A distinct new template is generated when the normal velocity is
increased by a critical value $\Delta v$,

$\Delta_{v}\approx\frac{v_{\perp}^{2}\Delta_{t}}{2R}.$ (20)

This velocity resolution improves for slower sweeps, as the template spreads
over more epochs, see Eq. (17).

Most of our observations are transferred to the GPS constellation in the
following section. We will return to this toy problem in the context of
computing template bank covariance matrices in Sec. IV.1.

###  III.2 GPS.DM search

Signal templates in the GPS.DM search are determined by the particular DM
model used (monopoles, strings, domain walls, etc…) and the parameters
associated with the event. For concreteness, we focus on DM wall events, which
are observed as a sweep through the geographically-distributed network, with a
given sensor perturbed when the wall overlaps with the sensor. The relevant
parameters for a particular event are the speed, incident direction, time of
the event, and the thickness of the wall. As in the toy problem, only the
normal component of the sweep velocity matters, however the velocity vector
now spans 3D. This parameter space is continuous and it is necessary to
discretize it to keep the number of templates finite. We discuss the
discretization strategies in Sec. III.3.

The standard halo model (SHM) dictates the necessary parameter prior
distributions. For DM objects the velocity distribution in the halo frame is
quasi-Maxwellian and isotropic with dispersion $v\approx 300\,\mathrm{km/s}$
[44] and with a cutoff set by the escape velocity of $v_{\mathrm{esc}}\approx
550\,\mathrm{km/s}$. Additionally, the solar system travels through the halo
at galactic velocities of $v_{g}\approx 200\,\mathrm{km/s}$ toward the Cygnus
constellation, implying that over $90\%$ of DM sweeps come from the forward
facing hemisphere [6]. In the GPS.DM search, we are interested in the normal
velocity of the DM wall in the Earth Centered Inertial (ECI) frame. In Fig. 3
we reproduce the normal velocity distribution for DM walls bound by the galaxy
in the ECI frame [6]. The maximum normal velocity, $v_{\perp}$, in this frame
is the addition of the DM escape velocity and the velocity of the solar
system, $v_{\perp}<v_{\mathrm{esc}}+v_{g}\approx 770\,\mathrm{km/s}$. As
previously stated, the minimum normal speed that can be resolved must be
larger than the velocities of the GPS satellites $v_{a}\approx
5\,\mathrm{km/s}$. This comes from the requirement that we can fix the
positions of orbiting satellites for a given time window. Similar to the toy
problem of Sec. III.1, the event parameters determine in which order the
sensors are swept and the overall sweep time. These characteristics
distinguish the templates within the signal template bank.

The GPS data-streams from satellite clocks are given as time measurement
biases with with respect to a fixed terrestrial reference sensor (clock) $R$.
DM-induced sensor bias (phase) of a given sensor $a$ is proportional to an
integral of the DM-induced frequency shift [18]. Since the clock biases are a
subject to a random-walk noise, we whiten the data by differencing the data
streams.

Both satellite and reference clocks are affected by the DM sweeps. In general,
sensitivity of the clock to the variation of fundamental constants depends on
the clock atom. The GPS constellation employs mostly Rb and a few Cs clocks. A
precision reference clock can be either an H-maser or a Rb clock (see Ref. [6]
for details) and one can choose among several terrestrial clocks to serve as
the reference clock. To simplify our consideration, we assume that we are only
dealing with a network of Rb clocks (including the reference clock). Then the
DM object affects each clock with the same signal amplitude $h$. This is the
limiting case of a homogeneous network. To bring out the essentials, we
further focus on “thin” walls characterized by an interaction time with each
sensor (clock) being much shorter than the sampling time interval
$\Delta_{t}$. This reduces the perturbation of atomic frequencies of each
clock to Kronecker deltas.

Figure 3:  (Color online) Normal speed $v_{\perp}$ probability distribution
[6] for galactically bound DM wall objects that cross the network. The
distribution has a maximum at $v_{\perp}=209\,\mathrm{km/s}$ and a cut-off at
$v_{\perp}=760\,\mathrm{km/s}$.

The terrestrial reference clock is centered in the network and we fix its
position in the center of our window of length $J$. Indeed, the reference
clock is located on the Earth, while the satellite clocks are orbiting
approximately 20 thousand kilometers above the Earth. $J$ is chosen to be an
odd number, so that the reference clock is perturbed in the epoch
$l_{R}=(J+1)/2$. Then the sought signal in the differenced data stream reads

$s_{l}^{a}=\delta_{l,l_{a}}+(-1)\delta_{l,l_{R}}\,,$ (21)

where $l$ enumerates epoch in the data window and $l_{a}$ and $l_{R}$ are the
epochs which the satellite clock and reference clock interact with the DM
object, respectively.

In general, the order in which the sensors are swept by the DM object depends
on the geometry of the network. In our setup, no matter the directionality of
the DM object some of the satellite sensors are affected prior to the
reference clock. However, for a terrestrial network, or when the reference
clock role is assigned to one of the satellite clocks, one can not ascertain
which sensors (reference or network) would be affected first in the window.
For this reason, we always set the reference sensor (on Earth or a satellite)
to be affected in the center on the window, thus mitigating any possibility of
missing information from a network sensor that is affected prior or after the
reference sensor.

The GPS.DM template construction can be visualized similar to the toy problem.
Fig. 4 shows a spherical shell of diameter $D=2R$ with sensors distributed on
the surface. A planar DM object passing through the network with normal
velocity $v_{\perp}$ and time resolution $\Delta_{t}$ segments the sphere into
$\Delta l=D/(v_{\perp}\Delta_{t})$ slices. As in our toy problem, $\Delta l$
has the meaning of template spread. Each segment has a spatial width
$L=v_{\perp}\,\Delta_{t}$. The normal velocity of thin walls determines the
number of equally spaced slices that cut the spherical shell.

Figure 4: (Color online) A spherical shell network sliced by a thin wall
propagating from the left. The wall is normal to the page. The shell is
segmented into $\Delta l$ slices and each slice segment of the sphere $S_{k}$
is defined by the polar angles $\theta_{k}$ and $\theta_{k+1}$ spanning its
width $L=v_{\perp}\Delta_{t}$. Angles are $\theta_{k}=\cos^{-1}(X_{k}/R)$ with
$X_{k}=\frac{2R}{\Delta l}(\frac{\Delta l}{2}-k)$. Reference clock is not
shown.

Fig. 5 shows a signal template for a wall sweep with normal speed
$v_{\perp}=185\,\mathrm{km/s}$. The network has $N_{\rm S}=25$ sensors and
data window contains $J=61$ epochs. The sensors are numbered in such a way
that the sensor $a=1$ is swept first and the sensors are swept in increasing
order. The slope of the blue line reflects the sweep speed. Notice, that once
the sensors are numbered, the sweep from a different direction would result in
a different pattern, with the straight red line converted into a dis-ordered
sequence of red tiles, see, e.g., Fig. 2. The blue vertical line would remain
the same with some of the tiles (corresponding to satellites degenerate with
the reference sensor) missing.

Figure 5: (Color online) Template for a thin wall signal sweeping through the
GPS constellation with the normal speed of $v_{\perp}=185\,\mathrm{km/s}$.
This network on January 3, 2012 has $N_{\rm S}=25$ satellite clocks and the
reference clock. Data window size is $J=61$. When the wall interacts with the
satellite clock, a positive spike of magnitude 1 (red tile) is observed. When
the wall interacts with the reference clock, there is a spike of opposite sign
(blue tile) at the epoch 31.

Our toy problem discussion (Sec. III.1) of template epoch spread, sensor
degeneracy, and angular and normal speed resolutions applies to the GPS
network search. The presence of the reference sensor requires refinement of
some of the considerations. For example, consider the limiting case when all
the satellite clocks are affected by the DM wall within the same $\Delta_{t}$
time interval (epoch) as the reference sensor. In this case, $l_{a}=l_{R}$ and
the signal template (21) collapses into a “null” template, where all elements
$s_{l}^{a}$ of the sought DM signal stream are zero. This effectively
eliminates detectable DM effects on the data stream. Even when only some of
the satellites are degenerate with the reference sensor, part of the network
sensitivity is effectively eliminated. This effect can be seen in Fig. 5,
where satellites 12 through 15 are degenerate with the reference clock and
lead to the lost (null) information from the satellites and reference sensor
(blank rows in Fig. 5).

Another peculiarity of the GPS network is that the network nodes move due to
the orbital motion of the satellites. The spatial positions of the orbiting
satellites is a part of the archival dataset [32] and this information is used
in the GPS.DM search. In analytical estimates, however, a more tractable
approach is to introduce a uniform occupancy approximation, where the sensors
populate a sphere of radius $R$ surrounding Earth with a random (but uniform
in probability) distribution of sensors on the sphere surface. In this
approximation, the probability for a sensor $a$ to reside on a spherical
surface segment $S_{k}$ segment is simply $P(a\in S_{k})=\Omega_{k}/4\pi$,
where solid angle $\Omega_{k}$ spans the segment $S_{k}$, see Fig. 4. Explicit
computations show that $\Omega_{k}=4\pi/\Delta l$, so that $P(a\in
S_{k})=1/\Delta l$, i.e. the probability is independent of the slice. The
average number of sensors in a given slice is the total number of network
sensor $N_{\rm S}$ multiplied by the probability $P(a\in S_{k})$ for a sensor
$a$ to reside on the segment $S_{k}$. Thereby, in this uniform occupancy
approximation, the degeneracy factor introduced in Sec. III.1 is the same as
for the 2D toy problem,

$\bar{g}=N_{\rm S}\frac{v_{\perp}\Delta_{t}}{2R}.$ (22)

As an illustration, template in Fig. 5 was generated for $N_{\rm S}=25$ sensor
network and sweep velocity $v_{\perp}=185\,\mathrm{km/s}$ with realistic
positions of GPS satellites. Based on Eqs. (17,22) we estimate a template
spread $\Delta l\approx 9$ with a degeneracy of $\bar{g}\approx 3$ sensors per
epoch. This is consistent with Fig. 5, which, depending on the epoch, exhibits
between 1 to 4 degenerate sensors. While not exact, the uniform occupancy is a
reasonable approximation for the true GPS network geometry.

The angular and speed resolutions remain the same as in the toy problem, Sec.
III.1. The angular resolution depends on just the normal speed of the DM wall,
while the normal speed resolution is inversely proportional to the square of
the normal speed,

$\displaystyle\Delta_{\theta}$
$\displaystyle\approx\sin^{-1}\left(\frac{v_{\perp}\Delta_{t}}{2R}\right)\,,$
(23) $\displaystyle\Delta_{v}$
$\displaystyle\approx\frac{v_{\perp}^{2}\Delta_{t}}{2R}.$ (24)

Because the velocity distribution is peaked at $v_{\perp}^{p}\approx
209\,\mathrm{km/s}$, we can introduce the most probable quantities, such as
the most probable template spread

$\overline{\Delta l}=2R/(v_{\perp}^{p}\Delta_{t})\,,$ (25)

which evaluates to $\overline{\Delta l}=8$ for GPS data.

###  III.3 Strategies in choosing template banks

Choosing templates in the template bank is a trade-off between computational
cost and the false-negative probability of missing the sought signal. Indeed,
one would like to provide dense coverage of the transient parameter space to
minimize false negatives, however increasing the number of templates $M$
raises computational cost of the search. We address this question in this
section and refine some of the statements further in the text.

There is a natural network resolution-limited choice of templates. For
example, in the GPS.DM search for thin walls, Sec. III.2, one may visualize
the template parameter space (normal speed, polar and azimuthal angles of
incidence for DM wall) as points in 3D space. Network resolutions
($\Delta_{v}$ and $\Delta_{\theta}$, Eqs. (24,23)) determine a small
“resolution volume” enclosing each point in this continuous space. Points
within the “resolution volume” lead to identical templates. Thereby, the
resolutions $\Delta_{\theta}$ and $\Delta_{v}$ effectively determine the grid
in this 3D parameter space. As an illustration, consider a thin wall sweeping
the GPS constellation at the most probable normal speed
$v_{\perp}=209\,\mathrm{km/s}$ (from Fig. 3). Then for
$\Delta_{t}=30\,\mathrm{s}$ and $2R\approx 5\times 10^{4}\,\mathrm{km}$, the
normal speed resolution is $\Delta_{v}\approx 26\,\mathrm{km/s}$ and an
angular resolution in both polar and azimuthal angles of incidence is
$\Delta_{\theta}\approx 7.2^{\circ}$. Complete resolution-limited coverage in
normal speed coordinate requires 68 points (with $25\,\mathrm{km/s}\leq
v_{\perp}<760\,\mathrm{km/s}$; the grid is non-uniform because $\Delta_{v}$
depends on $v_{\perp}$). For a fixed value of speed
$v_{\perp}=209\,\mathrm{km/s}$, 25 points for the polar angle and 50 points
for the azimuthal angle are required for resolution-limited angular coverage.
Covering the entire 3D cube in the template parameter space requires $\sim
2\times 10^{4}$ points. This is still computationally prohibitive in our
experience.

Since the most important criteria in choosing the bank is reducing the
probability of false negatives, the discussed resolution-limited bank
construction can be further augmented by the importance sampling technique
[52], where the parameter priors determine the density of sampling the
template parameter space. For example, in the GPS.DM search, the normal
velocity distribution, shown in Fig. 3, is a peaked distribution. Thus we need
more points on the $v_{\perp}$ grid near the maximum of the distribution and
less points in the tails to maintain the probability of false negatives
uniform throughout the bank. This can be accomplished by using a non-linear
grid or using Monte-Carlo sampling. To optimize the resulting grid, we would
further limit the spacing between the velocity grid points to be greater than
the resolution $\Delta_{v}$. The angular grids can be optimized with symmetry-
adapted choice of template bank that takes advantage of the network symmetry,
see Sec. IV.1.

Ref. [53] addressed a similar question of bank optimization in the context of
gravitational wave searches. They proposed to optimize the bank using the
singular value decomposition technique. However, these authors dealt with a
single sensor (not a network) and it remains to be seen if their technique can
be generalized to the network search.

An issue related to the network resolution, is the degree of similarity
between the templates. This degree of similarity will be quantified in the
following Sec. IV with cross-template correlation coefficients. The generated
bank can be optimized further based on this consideration. We also remind the
reader that the SNR-max statistic is the quantitative measure governing the
results of the search. The choice of template bank affects the SNR-max PDF,
CDF, and threshold values for false positives. These characteristics of the
SNR-max detection statistic is the focus of the remainder of this paper.

##  IV Template bank covariance matrix

Consider a template bank comprised of $M$ templates. Since the same data,
$\mathbf{d}^{w}$, is used to compute the SNR (3) $\rho_{k}$ for each template,
$\mathbf{s}_{k}$, then in general, different template SNRs are correlated,
i.e. cross-template correlation coefficients
$\Sigma_{ij}\equiv\langle\rho_{i}\rho_{j}\rangle\neq 0$. Because each SNR is a
Gaussian random variable with zero mean, see Sec. II, their joint PDF
$f\left(\bm{\rho}|\mathbf{\Sigma}\right)$ is given by the multivariate normal
distribution uniquely specified by the $M$-dimensional template bank
covariance matrix $\mathbf{\Sigma}$,

$f\left(\bm{\rho}|\mathbf{\Sigma}\right)=\frac{1}{\sqrt{\det(2\pi\bm{\Sigma}})}\exp\left(-\frac{1}{2}\bm{\rho}^{T}\bm{\Sigma}^{-1}\bm{\rho}\right)\,.$
(26)

Here $\bm{\rho}=\left(\rho_{1},...,\rho_{M}\right)^{T}$ and $\bm{\Sigma}^{-1}$
is the inverse of the covariance matrix $\bm{\Sigma}$.

We compute matrix elements of $\mathbf{\Sigma}$ in Appendix A. The result
reads

$\Sigma_{ij}=\frac{\mathbf{s}_{i}^{T}\mathbf{E}^{-1}\mathbf{s}_{j}}{\sqrt{\left(\mathbf{s}_{i}^{T}\mathbf{E}^{-1}\mathbf{s}_{i}\right)\left(\mathbf{s}_{j}^{T}\mathbf{E}^{-1}\mathbf{s}_{j}\right)}}\,,$
(27)

where the subscripts $i$ and $j$ identify templates. This expression depends
only on the templates and the noise covariance matrix $\mathbf{E}$. The
diagonal matrix elements of $\mathbf{\Sigma}$ are equal to 1, in agreement
with Ref. [31]. The off-diagonal matrix element vanishes if the templates are
orthogonal in the generalized sense,
$\mathbf{s}_{i}^{T}\mathbf{E}^{-1}\mathbf{s}_{j}=0$. Eq. (27) is simply
$(\mathbf{s}_{i}\cdot\mathbf{s}_{j})$ divided by the generalized norms of
$\mathbf{s}_{i}$ and $\mathbf{s}_{j}$, where the generalized dot product is
understood as
$(\mathbf{a}\cdot\mathbf{b})\equiv\mathbf{a}^{T}\mathbf{E}^{-1}\mathbf{b}$.

Qualitatively, the cross-template correlations $\Sigma_{ij}$’s quantify the
degree of similarity between two templates. If the two templates are
completely dis-similar (fully orthogonal), $\Sigma_{ij}=0$. If the two
templates are identical (up to an overall sign), $\Sigma_{ij}=\pm 1$. To
refine statements of Sec. III.3, it is desired to generate banks with the
smallest possible values of cross-template correlations to avoid testing data
streams against nearly-identical templates. In this sense, the most optimal
choice of the bank templates would be a fully-orthogonal bank. However, as we
show in Secs. IV.1 and IV.2, constructing a fully-orthogonal bank is not
possible for at least some of the applications.

The result (26) holds for any network with stationary white or colored
Gaussian noise and with and without cross-node correlations. It also holds for
any template bank choice. Namely the properties of $\mathbf{\Sigma}$ determine
the SNR-max statistic PDF, CDF, and threshold values. In the rest of this
section, we illustrate computation and properties of the template bank
covariance matrix with our toy problem of Sec. III.1 and the GPS.DM templates
of Sec. III.2.

###  IV.1 Toy problem and symmetry-adapted bank construction

Now we return to the toy problem of searching for line sweeps through a
circular topology network of Sec. III.1. In this problem, a template
$\mathbf{s}_{i}$ is uniquely specified by $(\bm{v}_{\perp})_{i}$, i.e. by the
incident direction $\hat{u}_{i}$ and the normal speed $(v_{\perp})_{i}>0$ of
the sweeping line. In component form,

$(s_{i})_{l}^{a}=\delta_{l,l_{a}^{i}},$ (28)

where $l_{a}^{i}$ is the epoch when sensor $a$ is crossed by the line moving
with $(\bm{v}_{\perp})_{i}$, see Eq. (16).

We are interested in properties of $\mathbf{\Sigma}$ as these determine SNR-
max CDF, PDF and threshold values. The bank covariance matrix elements, Eq.
(27), simplify to

$\Sigma_{ij}=\frac{1}{N_{\rm S}}\sum_{a=1}^{N_{\rm
S}}\delta_{l_{a}^{i},l_{a}^{j}}.$ (29)

The sum counts the number of sensors that are affected in the same epoch in
both templates. Apparently, $\Sigma_{ij}\geq 0$, as a result of the signals
being always positive. In addition, $\Sigma_{ij}\leq 1$, with $\Sigma_{ij}$
reaching the value of 1 only when both templates are identical.

We fix $v_{\perp}$ and discuss selections of incident directions $\hat{u}$ of
sweeping lines. The $C_{N_{\rm S}}$ symmetry of our network suggests a
symmetric choice of $\hat{u}_{i}:\hat{u}_{i}||\hat{n}_{a}$, i.e. sweeps
towards sensors. An example for an $M=5$ template bank is shown in Fig. 1
(right panel). Because of the rotational symmetry, all templates in the bank
can be easily generated by cyclic permutation of sensor labels. For the Fig. 1
bank, template 2 can be generated by relabeling sensors in the template 1,
Fig. 2 : $1\rightarrow 3,2\rightarrow 4,...,9\rightarrow 1,10\rightarrow 2$
and then swapping the rows of the resulting template, so that the sensors are
normally ordered.

We remind the reader that $\mathbf{\Sigma}$ is a symmetric matrix, with
diagonal elements $\Sigma_{ii}=1$. Based on the rotational symmetry arguments
for the constructed bank, we expect all matrix elements $\Sigma_{i,i\pm 1}$ to
be equal; the same applies to groups of matrix elements $\Sigma_{i,i\pm 2}$,
$\Sigma_{i,i\pm 3}$, etc. Here $i\pm n$ is understood as mod($i\pm n,M$). In
general, $\Sigma_{i,i\pm 1}\neq\Sigma_{i,i\pm 2}\neq\Sigma_{i,i\pm
3}\neq\ldots$. Such matrices belong to a class of symmetric Toeplitz matrices
that offer several computational advantages. This property is a reflection of
the underlying symmetry and would not hold for randomly-selected
${\hat{u}_{i}}$.

Our explicit computations for the $M=5$ bank of Fig. 1 with templates of Fig.
2 show that all off-diagonal matrix elements $\Sigma_{i,j\neq i}=1/5$, i.e.
all pairwise template overlaps have two simultaneously perturbed sensors. The
CDF and PDF of such template bank covariance matrices with $\Sigma_{i,j\neq
i}=r$ can be treated analytically, see Sec. V.4. We refer to such banks as
perfectly-squeezed. The CDF $C_{M=5}(Z|\mathbf{\Sigma}_{r=1/5})$ for such a
bank is given by Eq. (67).

Now we expand our bank of Fig. 1 to include all $N_{\rm S}=10$ possible
incident directions towards sensors. The discussed $M=5$ bank in Fig. 1 is
just a subset of this $M=10$ bank. We find that all off-diagonal matrix
elements of the bank covariance matrix between the Fig. 1 bank and newly added
templates vanish. The bank breaks into two mutually orthogonal sub-banks. Thus
$\mathbf{\Sigma}$ assumes a block diagonal form with $\Sigma_{i,j\neq i}=1/5$
for each block. Thereby, the joint probability distribution (26) factorizes
and the resulting SNR-max CDF is the product

$C_{M=10}(Z|\mathbf{\Sigma})=\\\
C_{M=5}(Z|\mathbf{\Sigma}_{r=1/5})C_{M=5}(Z|\mathbf{\Sigma}_{r=1/5}).$

We may continue the process of building the symmetry-adapted bank by
augmenting the $M=10$ bank with incidence directions $\hat{u}_{i}$ that
subdivide the arks between the adjacent $M=10$ bank directions into equal
parts, until we reach the network angular resolution $\Delta_{\theta}$,
introduced in Sec. III.1. This procedure relies on the intrinsic $C_{N_{\rm
S}}$ symmetry of our toy network. The presented ideas of symmetry-adapted bank
construction can be extended to 3D network with icosahedral point group
symmetry $I_{h}$. The GPS constellation, discussed in the following section,
unfortunately, lacks point symmetries and a different approach is required.

###  IV.2 GPS.DM search

Generation of GPS.DM templates for planar transients was discussed in Sec.
III.2. Compared to our toy network, GPS constellation has a reference clock,
it has a 3D geometry, lacks point symmetries, and positions of satellites
evolve in time.

We start with one of the un-optimized GPS.DM template banks used in the GPS.DM
search. The bank has $M=1024$ templates and was generated for the GPS
constellation orbital configuration for Feb. 14, 2005 08:20 UTC time. (We
specify time because the generated templates depend on the orbital positions
of the satellites, which vary throughout the day. Thus the template bank needs
to be re-generated every window.) The network had $N_{\rm S}=21$ sensors and
relative reference-sensor sensitivity parameter (12) $\xi=0.6$. Identical
white frequency noise was assumed for all the clocks. Templates were generated
by importance sampling using Monte-Carlo technique described in Sec. III, i.e.
the template event parameters were randomly drawn from the parameter space
weighed by priors. This bank generation did not take into account network
angular and normal speed resolutions.

Fig. 6 shows our computed distribution of the bank covariance matrix elements,
$\Sigma_{ij}$ for $i\neq j$. We observe that $\Sigma_{ij}$’s are centered
about the most probable value of $\Sigma_{ij}\approx 0.33$. All values are
positive and populate a range $0.0\leq\Sigma_{ij}\leq 1.0$. In this particular
template bank, there are 7 templates that are orthogonal to all other
templates, while some of the templates in the bank are nearly identical (
$\Sigma_{ij}\approx 1$), and 7 templates are identical, with $\Sigma_{ij}=1$.
Clearly, this is not an optimal choice for a template bank, as such aligned
(identical) templates would result in wasted computational time. A simple
solution to this problem would be to check for $\Sigma_{ij}\approx 1$ and if
so, to draw another template from the Monte-Carlo generator.

Figure 6: Distribution of cross-template correlation coefficients
$\Sigma_{i\neq j}$ of the template bank covariance matrix $\mathbf{\Sigma}$
for one of the GPS.DM template banks. Template bank contains $M=1024$
templates.

To quantify the distribution of $\Sigma_{ij}$ matrix elements and to develop
an analytical approach, in Appendix B we derive several properties of the
GPS.DM template bank covariance matrix for an idealized network of white noise
sensors with network noise covariance matrix $\mathbf{E}$ of Eq. (10). The
derivation is similar to that of Sec. IV.1 with an addition of the reference
sensor perturbation, so that the template elements are given by Eq. (21). The
derived template covariance matrix depends on a set of three parameters
$K_{1}$, $K_{2}$, and $R_{1}$:

$\displaystyle K_{1}^{s,s^{\prime}}$ $\displaystyle\equiv\sum_{a}^{N_{\rm
S}}\delta_{l_{a},l_{a}^{{}^{\prime}}},$ (30) $\displaystyle
K_{2}^{s,s^{\prime}}$ $\displaystyle\equiv\sum_{a}^{N_{\rm S}}\sum_{b\neq
a}^{N_{\rm S}}\delta_{l_{a},l_{b}^{{}^{\prime}}},$ $\displaystyle
R_{1}^{s,s^{\prime}}$ $\displaystyle\equiv N_{\rm
S}-\sum_{a}\left(\delta_{l_{R},l_{a}}+\delta_{l_{R},l_{a}^{{}^{\prime}}}\right)\,,$

where $l_{a}$ are epochs when sensor $a$ was perturbed in template
$\mathbf{s}$ and primed quantities refer to template $\mathbf{s}^{\prime}$.
Here $K_{1}$ is the number of sensors with a perturbation epoch identical to
both templates, $K_{2}$ is the cross-sensor degeneracy between the templates,
and $R_{1}$ quantifies the reference sensor and network sensor degeneracy.
These expressions apply to both same-template
($\mathbf{s}\equiv\mathbf{s}^{\prime}$) and the cross-template cases.

With these definitions, the cross-template correlation coefficients between
two templates $\mathbf{s}$ and $\mathbf{s}^{\prime}$ read

$\displaystyle\Sigma_{ss^{\prime}}=$ (31)
$\displaystyle\frac{aK_{1}^{s,s^{\prime}}+R_{1}^{s,s^{\prime}}+bK_{2}^{s,s^{\prime}}}{\sqrt{\left(aK_{1}^{s,s}+R_{1}^{s,s}+bK_{2}^{s,s}\right)\left(aK_{1}^{s^{\prime},s^{\prime}}+R_{1}^{s^{\prime},s^{\prime}}+bK_{2}^{s^{\prime},s^{\prime}}\right)}},$

where template-independent parameters are

$\displaystyle a$ $\displaystyle=1+(1-1/N_{\rm S})\xi\,,$ (32) $\displaystyle
b$ $\displaystyle=-\xi/N_{\rm S}\,.$ (33)

Here $\xi$ is the relative reference-sensor sensitivity parameter (12).

With these definitions, in Appendix B we determine “averaged” cross-template
correlation value $\overline{\Sigma}_{i\neq j}$ for our GPS network data. Here
“averaged” means that the degeneracy parameters (30) are equated to the
typical values using the most probable epoch spread $\overline{\Delta l}$, Eq.
(17) and the uniform occupancy approximation of Sec. III.2. The average
degeneracy parameters, computed in Appendix B, are
$\bar{K}_{1}^{s,s}=\bar{K}_{1}^{s^{\prime},s^{\prime}}=N_{\rm S}$,
$\bar{K}_{1}^{s,s^{\prime}}={N_{\rm S}}/\overline{\Delta l}$,
$\bar{K}_{2}^{s,s}=\bar{K}_{2}^{s^{\prime},s^{\prime}}=\bar{K}_{2}^{s,s^{\prime}}=(N_{\rm
S}-1){N_{\rm S}}/\overline{\Delta l}$, and
$\bar{R}_{1}^{s,s}=\bar{R}_{1}^{s^{\prime},s^{\prime}}=\bar{R}_{1}^{s,s^{\prime}}=\left(1-{2}/{\overline{\Delta
l}}\right)N_{\rm S}$. These values result in the “averaged” cross-template
correlation coefficient

$\overline{\Sigma}_{i\neq j}\approx\frac{1}{2+(1-1/N_{\rm S})\xi}\,.$ (34)

To access the accuracy of approximation (34), we consider the computed
distribution of cross-template correlation coefficients $\Sigma_{ij}$ of Fig.
6. This histogram was computed for the network of $N_{\rm S}=21$ Rb clocks
with sensitivity parameter $\xi=0.6$. Then Eq. (34) results in
$\overline{\Sigma}_{i\neq j}=0.388$. It is in $\sim 5\%$ agreement with the
most probable value inferred from the histogram, $\overline{\Sigma}_{i\neq
j}=0.398$. Over decades of the GPS constellation operation, satellites are
replaced as they reach their end of life. The most typical value of
$\xi\approx 0.3$, and then for $N_{\rm S}=30$ the Eq. (34) gives
$\overline{\Sigma}_{i\neq j}=0.44$.

We discussed template bank optimization strategies in Sec. III.3. One of the
refinements of these strategies is to “squeeze” the distribution of off-
diagonal template covariance matrix elements around the average value of
$\Sigma_{i\neq j}$, $r$, with some small spread $\delta\mathbf{\Sigma}$, so
that $\mathbf{\Sigma}=\mathbf{\Sigma}_{r}+\delta\mathbf{\Sigma}$, where off-
diagonal matrix elements of $\mathbf{\Sigma}_{r}$ are set to $r$ and
$\left(\Sigma_{r}\right)_{ii}=1$. One “squeezing” technique that can be used
with success (although with large computational cost) is by performing Monte-
Carlo and eliminating templates whose correlation element is above or below a
specified spread in correlation. In Sec. V.4 we provide analytical results for
the SNR-max CDF and PDF of such $\mathbf{\Sigma}_{r}$; the bank “squeezing”
improves the accuracy of approximating $\mathbf{\Sigma}$ with
$\mathbf{\Sigma}_{r}$. Fig. 7 shows an example of “squeezed” bank distribution
of template covariance matrix elements (not including diagonal 1’s) for the
most probable covariance matrix element $r=0.39$ and spread of
$\delta{\Sigma}_{ij}=\pm 0.2$ of the perturbation matrix.

Figure 7: (Color online) Distribution of off-diagonal matrix elements
$\Sigma_{ij}$ of the template bank covariance matrix for a squeezed template
bank with the most probable covariance matrix element $r=0.39$ and maximum
spread in matrix elements $\delta\Sigma_{ij}=\pm 0.2$. Template bank contains
$M=1024$ templates. All diagonal matrix elements $\Sigma_{ii}=1$ and are not
accounted for here.

##  V SNR-max distribution for a template bank

To reiterate, at this point we have discussed the properties of templates and
the template bank covariance matrix $\mathbf{\Sigma}$. Now we return to the
discussion of the SNR-max statistic and derive its distributions. The starting
point is the joint probability distribution (26) for SNRs in the template
bank. The joint PDF is fully characterized by the template bank covariance
matrix $\mathbf{\Sigma}$, which is a symmetric and positive semi-definite
matrix with diagonal matrix elements $\Sigma_{ii}=1$. Off-diagonal matrix
elements of $\mathbf{\Sigma}$ are given by Eq. (27). In this section, we
determine an analytic form for the CDF and PDF of the SNR-max statistic (6),
$z=\max(|\rho_{1}|,..,|\rho_{M}|)$. Apparently, $z\geq 0$.

We start the derivation by computing the CDF of the SNR-max statistic. CDF
$C_{M}\left(Z|\bm{\Sigma}\right)$ is the probability that the SNR-max
statistic $z$ takes a value less than or equal to $Z$. The condition $z\leq Z$
implies that all the individual template SNRs lie in the interval
$-Z\leq\rho_{k}\leq Z$. Thereby, we can express the CDF as an integral of the
joint probability distribution (26),

$\displaystyle C_{M}\left(Z|\bm{\Sigma}\right)=$ (35)
$\displaystyle\frac{1}{\sqrt{\det(2\pi\bm{\Sigma})}}\int_{-Z}^{+Z}d\rho_{1}...\int_{-Z}^{+Z}d\rho_{M}\exp\left(-\frac{1}{2}\bm{\rho}^{T}\bm{\Sigma}^{-1}\bm{\rho}\right)\,.$

This $M$-dimensional integral can be visualized as an integral of a
multivariate normal distribution over an $M$-dimensional hyper-cube centered
at $0$ and having a side of $2Z$.

To date there is no closed form for the CDF of a multivariate normal
distribution, so the $M$-dimensional integral requires numerical evaluation
(typically using Monte-Carlo techniques [54, 55]). For dynamically evolving
networks, such as the GPS constellation where satellites are moving, this
approach is computationally prohibitive, as templates change with the network
geometry, and have to be generated anew for every window. Every newly
generated template bank would have a different template bank covariance matrix
$\bm{\Sigma}$, new CDF, and a new threshold for false-positives. Partially,
this problem can be mitigated by observing that there is a 12-hour periodicity
in the satellite orbits, i.e. network geometry repeats every 12 hours. Then
the computations can be reduced by storing template banks and threshold values
for the 12 hour period. Even with this shortcut, one still has to deal with a
collection of 1440 template banks for individual $\Delta_{t}=30\,\mathrm{s}$
epochs.

The cumulative distribution can be expressed using the underlying probability
density $p_{M}\left(\rho|\bm{\Sigma}\right)$,
$C_{M}\left(Z|\bm{\Sigma}\right)=\int_{0}^{Z}p_{M}\left(z|\mathbf{\Sigma}\right)dz$,
or

$p_{M}\left(z|\mathbf{\Sigma}\right)=\left(\frac{d}{dZ}C_{M}\left(Z|\mathbf{\Sigma}\right)\right)_{Z\rightarrow
z}.$ (36)

We carry out the differentiation of the SNR-max CDF (37) in Appendix C. The
resulting SNR-max PDF reads

$\displaystyle p_{M}\left(z|\bm{\Sigma}\right)=$
$\displaystyle\frac{2}{\sqrt{\det\left(2\pi\bm{\Sigma}\right)}}\sum_{k=1}^{M}\left({\displaystyle\prod\limits_{n=1,n\neq
k}^{M}}\int_{-z}^{+z}d\rho_{n}\right)\times$ (37)
$\displaystyle\times\left[\exp\left(-\frac{1}{2}\bm{\rho}^{T}\bm{\Sigma}^{-1}\bm{\rho}\right)\right]_{\rho_{k}=z}.$

In general, this SNR-max PDF is sufficient, yet it still involves $M$
evaluations of $(M-1)$-dimensional integrals. Thus, for large banks this PDF
still requires computationally-expensive numerical evaluation.

These considerations force us to explore approximate techniques for evaluating
SNR-max CDF and PDF. In Sec. V.1 we review the SNR-max CDF and PDF for
uncorrelated templates or fully-orthogonal bank ($\bm{\Sigma}=\bm{I}$). In
Sec. V.2, we consider a simple case of two correlated templates. In Sec. V.3
we derive an approximate SNR-max PDF and CDF for a nearly-orthogonal template
bank, i.e. when off-diagonal $|\Sigma_{ij}|\ll 1$. Lastly, in Sec. V.4 we
consider the SNR-max PDF and CDF for a “squeezed” template bank, when off-
diagonal matrix elements $\Sigma_{ij}\approx r$ for all templates.

Finally, for completeness, we mention an apparent simplification when
$\bm{\Sigma}$ has a block-diagonal form, a case encountered in the toy
problem, Sec. IV.1. This happens when the template bank can be decomposed into
mutually orthogonal sub-banks of dimensions $M^{\prime}$ and
$M^{\prime\prime}$, $M^{\prime}+M^{\prime\prime}=M$. If $\bm{\Sigma}$ contains
two blocks $\bm{\Sigma}^{\prime}$ and $\bm{\Sigma}^{\prime\prime}$, spanning
SNR sub-spaces $\bm{\rho}^{\prime}$ and $\bm{\rho}^{\prime\prime}$, the joint
PDF (26) factorizes
$f\left(\bm{\rho}|\mathbf{\Sigma}\right)=f\left(\bm{\rho}^{\prime}|\mathbf{\Sigma}^{\prime}\right)f\left(\bm{\rho}^{\prime\prime}|\mathbf{\Sigma}^{\prime\prime}\right)$.
Then the CDF factorizes as well,

$C_{M}\left(Z|\bm{\Sigma}\right)=C_{M^{\prime}}\left(Z|\bm{\Sigma}^{\prime}\right)C_{M^{\prime\prime}}\left(Z|\bm{\Sigma}^{\prime\prime}\right)\,.$

###  V.1 Fully-orthogonal bank (independent templates)

For a fully-orthogonal bank, all the off-diagonal matrix elements of the
template bank covariance matrix $\bm{\Sigma}$ vanish and $\bm{\Sigma}$ reduces
to an identity matrix $\bm{I}$. Then the general SNR-max CDF (35) is simply

$\displaystyle C_{M}\left(Z|\bm{I}\right)$
$\displaystyle=\frac{1}{\left(2\pi\right)^{M/2}}\left(\int_{-Z}^{+Z}d\rho\exp\left(-\frac{1}{2}\rho^{2}\right)\right)^{M}$
(38)
$\displaystyle=\left[\text{erf}\left(\frac{Z}{\sqrt{2}}\right)\right]^{M}\,,$

where Gaussian error function
$\mathrm{erf}\left(x\right)=(2/\sqrt{\pi})\int_{0}^{x}du\,e^{-u^{2}}$.
Further, differentiating CDF (38) with respect to $Z$, we arrive at the SNR-
max PDF for a bank of $M$ independent templates

$p_{M}\left(z|\bm{I}\right)=M\sqrt{\frac{2}{\pi}}e^{-\frac{z^{2}}{2}}\left(\text{erf}\left(\frac{z}{\sqrt{2}}\right)\right)^{M-1}.$
(39)

These results subsume the limiting case of the $M=1$ bank. In this case, we
are dealing with a single template, so that the SNR-max statistic (6),
$z=|\rho|$, is an absolute value of the SNR for that one template. From the
SNR distribution (5) we conclude that $z$ is a half-normally distributed
variable with probability distribution

$p_{1}(z)=\sqrt{\frac{2}{\pi}}e^{-\frac{z^{2}}{2}}\,.$ (40)

Integrating it, we obtain CDF $C_{1}(Z)=\mathrm{erf}\left(Z/\sqrt{2}\right)$.
These results are in agreement with Eqs. (38,39) for $M=1$.

###  V.2 Bank of $M=2$ correlated templates

With the results for a fully-orthogonal bank established, we start quantifying
the role of cross-template correlations with an analytically tractable case of
a bank consisting of two templates. For the $M=2$ case, the cross-template
covariance matrix depends only on one off-diagonal matrix element
$\Sigma_{12}=\Sigma_{21}\equiv r$. Explicitly,

$\displaystyle\bm{\Sigma}$
$\displaystyle\bm{=}\left(\begin{array}[c]{cc}1&r\\\ r&1\end{array}\right),$
(43) $\displaystyle\bm{\Sigma}^{-1}$
$\displaystyle=\frac{1}{1-r^{2}}\left(\begin{array}[c]{cc}1&-r\\\
-r&1\end{array}\right),$ (46)

and $\det\mathbf{\Sigma}=1-r^{2}$. Then the joint probability distribution
(26) reads

$f\left(\rho_{1},\rho_{2}|r\right)=\frac{1}{2\pi\sqrt{1-r^{2}}}\exp\left(-\frac{\rho_{1}^{2}+\rho_{2}^{2}-2r\rho_{1}\rho_{2}}{2(1-r^{2})}\right).$

Applying Eq. (37), we arrive at the SNR-max PDF for the $M=2$ bank

$\displaystyle p_{2}\left(z|\,r\right)$
$\displaystyle=\sqrt{\frac{2}{\pi}}e^{-\frac{z^{2}}{2}}\times$ (47)
$\displaystyle\left(\text{erf}\left(\frac{1}{\sqrt{2}}\sqrt{\frac{1+r}{1-r}}z\right)+\text{erf}\left(\frac{1}{\sqrt{2}}\sqrt{\frac{1-r}{1+r}}z\right)\right).$

To support our analytical result, in Fig. 8 we compare results of a Monte-
Carlo simulation with the derived SNR-max PDF (47), and we find the histogram
and the PDF in agreement. We carried out these simulations using a single
($N_{\rm S}=1$) white noise (Gaussian) data stream with zero mean and
$\sigma^{2}=1$, and with the two templates having a correlation coefficient
$r=1/2$. The positive skew observed in the SNR-max PDF is due to the effect of
templates being correlated. False-positive thresholds determined for
correlated template banks are affected by the positive skew observed in this
SNR-max PDF, see Sec. VI.

Figure 8: (Color online) Comparison of the Monte-Carlo simulation (blue) with
the closed-form SNR-max PDF (47) (red curve) for a bank of $M=2$ templates
with correlation coefficient $r=1/2$.

We observe the following properties of the PDF (47). The distribution
$p_{2}\left(z|r\right)$ is symmetric with respect to flipping the sign of $r$,
i.e. $p_{2}\left(z|-r\right)=p_{2}\left(z|r\right)$. If there is no
correlation between the two SNR variables ($r=0$), Eq. (47) reduces to

$p_{2}\left(z|r=0\right)=2\sqrt{\frac{2}{\pi}}e^{-\frac{z^{2}}{2}}\text{erf}\left(\frac{z}{\sqrt{2}}\right)\,,$
(48)

which is identical to the general result (39) for the $M=2$ orthogonal
template bank. For fully correlated SNRs, $r=\pm 1$, the probability
distribution tends to the half-normal distribution (40)

$p_{2}\left(z|\,r\rightarrow\pm
1\right)=\sqrt{\frac{2}{\pi}}e^{-\frac{z^{2}}{2}}=p_{1}\left(z\right).$ (49)

In this case, the two templates are identical (or differ by an overall sign)
and the PDF reduces to that of a single template case Eq. (49). These two
cases of $r=0$ and $r=\pm 1$ agree with previous literature [56].

To further illustrate the effect of correlation on the SNR-max PDF, we vary
the correlation coefficient $r$ from 0 to 1, see Fig. 9. Due to the
$r\rightarrow-r$ symmetry of the distribution, positive and negative template
correlations have identical PDFs. Fig. 9 shows that when the correlation
$|r|>0$, the SNR-max PDF becomes positively skewed (skewed toward zero). When
the templates are fully correlated, $|r|=1$, the distribution $p_{2}$ reduces
to Eq. (49), i.e. the half-normal distribution $p_{1}$ characteristic of a
single template bank. This is expected qualitatively as the two random
variables $\rho_{1}$ and $\rho_{2}$ become identical, so that
$z=\max(|\rho_{1}|,|\rho_{2}|)=|\rho_{1}|=|\rho_{2}|$.

Figure 9:  (Color online) Probability distribution for SNR max statistic for
$M=2$ bank with varying cross-template correlation, $0\leq r\leq 1$ ( $r=0$
blue curve, $r=0.75$ red curve, and $r=1$ orange curve, respectively). As the
magnitude of the correlation coefficient increases from $r=0$ to larger $|r|$,
the distribution becomes increasingly positively skewed, and when $r=\pm 1$
(orange curve) we recover the half-normal SNR-max distribution for $M=1$ bank,
see Eq. (49).

For a template bank composed of two templates, $M=2$, the SNR-max CDF
$C_{2}(Z|r)$ can be determined by numerically integrating the two-template
SNR-max PDF in Eq. (47)

$\displaystyle C_{2}(Z|r)=$
$\displaystyle\int_{0}^{Z}dz\sqrt{\frac{2}{\pi}}e^{-\frac{z^{2}}{2}}\times$
(50)
$\displaystyle\left(\text{erf}\left(\frac{1}{\sqrt{2}}\sqrt{\frac{1+r}{1-r}}z\right)+\text{erf}\left(\frac{1}{\sqrt{2}}\sqrt{\frac{1-r}{1+r}}z\right)\right).$

An equivalent representation,

$\displaystyle
C_{2}(Z|r)=\frac{1}{4}\int_{-\infty}^{\infty}du\,\frac{e^{-u^{2}/2}}{\sqrt{2\pi}}\times$
(51)
$\displaystyle\left[\mathrm{erf}\left(\frac{Z+\sqrt{r}u}{\sqrt{2(1-r)}}\right)-\mathrm{erf}\left(\frac{-Z+\sqrt{r}u}{\sqrt{2(1-r)}}\right)\right]^{2}\,,$

is derived in Sec. V.4.

Fig. 10 shows the CDF $C_{2}(Z|r=1/2)$ resulting from evaluating, Eq. (51)
numerically. From this CDF, thresholds for specific false positive rates will
be inferred in Sec. VI.

Figure 10: (Color online) Cumulative distribution for a template bank composed
of $M=2$ templates with a correlation coefficient $r=1/2$.

###  V.3 SNR-max distribution for a nearly-orthogonal bank

Suppose that the correlations between various templates are small, i.e.
$|\Sigma_{ij}|\ll 1$ for all $i\neq j$, such that the template covariance
matrix $\mathbf{\Sigma}$ can be represented as a sum of the identity matrix
$\bm{I}$ (fully-orthogonal template bank) and a perturbation matrix $\bm{A}$.
In this case, the bulk of the SNR-max CDF and PDF is determined by the
distributions for fully-orthogonal bank, Eqs. (38, 39), with corrections that
depend on $\bm{A}$:

$\displaystyle C_{M}\left(Z|\bm{\Sigma}\right)\approx$ $\displaystyle
C_{M}\left(Z|\bm{\Sigma}=\bm{I}\right)+\delta_{I}C_{M}\left(Z|\bm{\Sigma}\right)\,,$
(52) $\displaystyle p_{M}\left(z|\bm{\Sigma}\right)\approx$ $\displaystyle
p_{M}\left(z|\bm{\Sigma=I}\right)+\delta_{I}p_{M}\left(z|\bm{\Sigma}\right)\,.$
(53)

Here we use the symbol $\delta_{I}$ to emphasize that the expansion is about
$\bm{\Sigma}=\bm{I}$. In Appendix D, we demonstrate that these corrections are
proportional to the trace of ${\bm{A}^{2}}$. Below we summarize the derivation
and results from Appendix D.

The general expression for the SNR-max CDF (35) depends on the determinant and
the inverse of the template bank covariance matrix. The determinant can be
approximated as

$\det\left(\bm{I}+\bm{A}\right)\approx
1-\frac{1}{2}\mathrm{Tr}{\left(\bm{A}^{2}\right)}.$ (54)

Here we omitted terms containing $\mathrm{Tr}{\left(\bm{A}\right)}=0$ as, by
construction, the diagonal matrix elements of the perturbation matrix $\bm{A}$
vanish. Notice that
$\mathrm{Tr}{\left(\bm{A}^{2}\right)}=\mathrm{Tr}(\bm{\Sigma}^{2})-M$.
Similarly, the inverse of $\bm{\Sigma}$ can be approximated as

$\bm{\Sigma}^{-1}=\left(\bm{I}+\bm{A}\right)^{-1}\approx\bm{I}-\bm{A}+\bm{A}^{2}.$
(55)

Based on these approximations and Taylor expansion of the exponential in the
SNR-max CDF (35), in Appendix D we derived the correction to the SNR-max CDF,

$\displaystyle\delta_{I}C_{M}\left(Z|\bm{\Sigma}\right)$ $\displaystyle=$ (56)
$\displaystyle\left(\operatorname{Tr}(\bm{\Sigma}^{2})-M\right)\frac{Z^{2}e^{-Z^{2}}}{2\pi}\left[\text{erf}\left(\frac{Z}{\sqrt{2}}\right)\right]^{M-2}\,.$

Similarly, the correction to the SNR-max PDF is found by taking the derivative
of the above correction to the CDF, leading to

$\displaystyle\delta_{I}p_{M}\left(z|\bm{\Sigma}\right)=\left(\mathrm{Tr}(\bm{\Sigma}^{2})-M\right)\frac{ze^{-z^{2}}}{\pi}\left(\operatorname{erf}\left(\frac{z}{\sqrt{2}}\right)\right)^{M-3}$
$\displaystyle\times\quad\left[\frac{(M-2)ze^{-z^{2}/2}}{\sqrt{2\pi}}+\left(1-z^{2}\right)\operatorname{erf}\left(\frac{z}{\sqrt{2}}\right)\right]\,.$
(57)

The advantage of these approximate formulas is that the multi-dimensional
integrations in the general SNR-max CDF and PDF expressions (35, 37) have been
carried out explicitly. If the bank changes due to time evolution of network
geometry, e.g., for orbiting satellites, only $\mathrm{Tr}(\bm{\Sigma}^{2})$
needs to be reevaluated. There is, however, a region of validity of these
approximations (see Appendix D),

$z\ll 1/\sqrt{\max\left|A_{ij}\right|}\,,$ (58)

which limits practical applicability of these approximations.

As a test of these approximate formulae, we consider a limiting case of two-
template bank, see Sec. V.2, with the closed-form SNR-max PDF (47). In this
case, $\mathrm{Tr}(\bm{\Sigma}^{2})-M=2r^{2}$ and Eq. (V.3) becomes

$\delta_{I}p_{2}\left(z|r\right)=2r^{2}\left(1-z^{2}\right)\frac{ze^{-z^{2}}}{\pi}\,.$
(59)

This result is in agreement with the small-$r$ expansion of Eq. (47).

###  V.4 SNR-max distribution for a squeezed template bank

In Sec. IV we encountered a practical search scenario when matrix elements of
the bank covariance matrix $\mathbf{\Sigma}$ are centered about a non-zero
value $r$ (see Fig. 6). Motivated by this scenario, in this section we
consider an approximation for the SNR-max CDF when the off-diagonal matrix
elements of $\mathbf{\Sigma}$ exhibit small deviations from some value $r$.
(As follows from the derivation of this section, the optimal choice of $r$
corresponds to the average of off-diagonal matrix elements of
$\mathbf{\Sigma}$.) This is the case of squeezed banks discussed in Sec. IV.
Then $\mathbf{\Sigma}=\mathbf{\Sigma}_{r}+\delta\mathbf{\Sigma}$, where the
perturbation matrix $\delta\mathbf{\Sigma}$ is “small” compared to the highly-
structured reference covariance matrix $\mathbf{\Sigma}_{r}$. To the leading
order in $\delta\mathbf{\Sigma}$, the SNR-max CDF
$C_{M}(Z|\mathbf{\Sigma}_{r}+\delta\mathbf{\Sigma})$ is determined by
$C_{M}(Z|\mathbf{\Sigma}_{r})$. Similar arguments apply to PDF
$p_{M}(Z|\mathbf{\Sigma}_{r})$. Below, we evaluate these leading order
contributions, leaving the cumbersome derivation of the
$\delta\mathbf{\Sigma}$-correction for future work.

Explicitly, we assume that matrix elements of the reference matrix
$\mathbf{\Sigma}_{r}$ are given by

$\left(\mathbf{\Sigma}_{r}\right)_{ij}=\delta_{ij}+(1-\delta_{ij})r\,,$ (60)

i.e., all off-diagonal matrix elements of $\mathbf{\Sigma}_{r}$ are identical
and equal to $r$. For example, $\mathbf{\Sigma}_{r}$ for $M=3$ reads

$\mathbf{\Sigma}_{r}=\left(\begin{array}[c]{ccc}1&r&r\\\ r&1&r\\\
r&r&1\end{array}\right)\,.$ (61)

In Sec. V.2 we considered a special case of such matrices for $M=2$. It is the
highly-structured nature of $\mathbf{\Sigma}_{r}$ that allows us to derive a
closed-form expression for the SNR-max CDF $C_{M}(Z|\mathbf{\Sigma}_{r})$
below.

We derived the inverse of $\mathbf{\Sigma}_{r}$ entering the multivariate
normal distribution. The inverse retains the general structure of
$\mathbf{\Sigma}_{r}$ with matrix elements

$\left(\mathbf{\Sigma}_{r}^{-1}\right)_{ij}=X\delta_{ij}+Y\,,$ (62)

where

$\displaystyle X$ $\displaystyle=1/(1-r)\,,$ (63) $\displaystyle Y$
$\displaystyle=-\frac{r}{1+(M-1)r}X\,.$ (64)

The derived inverse $\mathbf{\Sigma}_{r}^{-1}$ reproduces that for the $M=2$
case, Eq. (46). With the matrix inverse (62), the joint PDF (26) reads

$\displaystyle
f\left(\bm{\rho}|\mathbf{\Sigma}_{r}\right)=\frac{1}{\sqrt{\det(2\pi\mathbf{\Sigma}_{r})}}\times$
(65)
$\displaystyle\exp{\left\\{-\frac{1}{2(1-r)}\left(\sum_{i=1}^{M}\rho^{2}_{i}-\frac{r}{1+(M-1)r}\left(\sum_{i=1}^{M}\rho_{i}\right)^{2}\right)\right\\}}\,,$

with $\det(\mathbf{\Sigma}_{r})=\left(1-r\right)^{M-1}\left(1+(M-1)r\right)$.
This joint PDF is equivalent to

$\displaystyle
f\left(\bm{\rho}|\mathbf{\Sigma}_{r}\right)=\frac{1}{(2\pi(1-r))^{M/2}}\times$
(66)
$\displaystyle\int_{-\infty}^{+\infty}du\quad\frac{e^{-u^{2}/2}}{\sqrt{2\pi}}\prod_{i=1}^{M}\exp\left(-\frac{\left(\rho_{i}-\sqrt{r}u\right)^{2}}{2(1-r)}\right)\,,$

which can be verified by direct integration. Remarkably, this integral
representation has the advantage of decoupling the SNRs $\rho_{i}$.

The properties of multi-variate normal distribution for the class of
covariance matrices represented by Eq. (60) have been considered in the
literature (see Ref. [57] and references therein). We warn the reader,
however, that the literature contains misprints (for example, Eq. (65) has
misprints in Ref. [57]), uses different notation, and proceeds through the
less transparent (at least to our taste) “exchangeable-normal-variables”
proofs.

The CDF $C_{M}(Z|\mathbf{\Sigma}_{r})$ is still an $M$-dimensional integral
(35) of the multivariate normal distribution (65). With the alternative
representation (66) of the joint PDF, this integral can be reduced to a one-
dimensional integral,

$\displaystyle
C_{M}(Z|\mathbf{\Sigma}_{r})=\frac{1}{2^{M}}\int_{-\infty}^{\infty}du\,\frac{e^{-u^{2}/2}}{\sqrt{2\pi}}\times$
(67)
$\displaystyle\left[\mathrm{erf}\left(\frac{Z+\sqrt{r}u}{\sqrt{2(1-r)}}\right)-\mathrm{erf}\left(\frac{-Z+\sqrt{r}u}{\sqrt{2(1-r)}}\right)\right]^{M}.$

This CDF is equivalent to the result of Ref. [57] expressed in a different
notation. For $r=0$, we recover the CDF for independent templates, Eq. (38).

Differentiating this CDF with respect to $Z$, we determine the SNR-max PDF

$\displaystyle
p_{M}(z|\mathbf{\Sigma}_{r})=\frac{1}{2^{M}}\frac{M}{\pi}\frac{1}{\sqrt{1-r^{2}}}\int_{-\infty}^{\infty}du\,e^{-u^{2}/2}\times$
(68)
$\displaystyle\left[\mathrm{erf}\left(\frac{z+\sqrt{r}u}{\sqrt{2(1-r)}}\right)-\mathrm{erf}\left(\frac{-z+\sqrt{r}u}{\sqrt{2(1-r)}}\right)\right]^{M-1}$
$\displaystyle\times\left[\exp\left(-\frac{(z+\sqrt{r}u)^{2}}{2(1-r)}\right)+\exp\left(-\frac{(z-\sqrt{r}u)^{2}}{2(1-r)}\right)\right]\,.$

Given the bank covariance matrix $\bm{\Sigma}$, what is the optimal choice of
$r$? The answer to this question is not obvious - there are many possibilities
like taking $r$ as the most probable value, or the median, or the average of
the $\Sigma_{i\neq j}$ matrix elements distribution. Below we show that the
optimal choice of $r$ is the mean value of off-diagonal matrix elements of the
bank covariance matrix $\bm{\Sigma}$. Indeed, consider the
small-$\delta\mathbf{\Sigma}$ expansion of the determinant normalizing the
joint probability distribution,

$\displaystyle\det(\mathbf{\Sigma}_{r}+\delta\mathbf{\Sigma})\approx$
$\displaystyle\left(1+\mathrm{Tr}{(\mathbf{\Sigma}_{r}^{-1}\delta\mathbf{\Sigma})}\right)\det\mathbf{\Sigma}_{r}$
$\displaystyle=$
$\displaystyle\left(1+Y\sum_{ij}(\delta\mathbf{\Sigma})_{ij}\right)\det\mathbf{\Sigma}_{r}\,,$

where $Y$ is the off-diagonal matrix element of $\mathbf{\Sigma}^{-1}_{r}$
given in Eq. (64). To suppress the $\delta\mathbf{\Sigma}$ correction, we have
to require that $r$ is chosen in such a way that
$\sum_{ij}(\delta\mathbf{\Sigma})_{ij}=0$. This condition leads to the optimal
choice

$r=\frac{1}{M(M-1)}\sum_{i>j}\Sigma_{ij}\,,$ (69)

corresponding to the mean value of off-diagonal matrix elements of the bank
covariance matrix.

Threshold values for specific values of false positive rates for squeezed
template banks are discussed in Sec. VI.

##  VI Threshold values for the SNR-max statistic

Practical applications require certain level of confidence in the detection
statistic. Typically, confidence levels such as $90\%,95\%$, and $99\%$ are
suitable for many applications requiring less sensitivity far into the tails
of a distribution. Threshold values are determined by false positive rates,
$q$, and for the enumerated confidence levels, $q$ values are $0.1,0.05$, and
$0.01$, respectively. The meaning of a false positive rate is the probability
of falsely accepting the hypothesis in question, when the data should have
supported the null hypothesis. In the context of searches for transients, a
false positive happens when the SNR-max statistic $z$ exceeds a certain
threshold $Z^{\ast}$, but this in fact, is due to statistical randomness
mimicking the transient and not due to the sought exotic physics. Suppose for
a desired false positive rate, $q$, we want to determine the threshold value
$Z^{\ast}$ associated with that false positive rate. We may simply relate the
false positive rate to the threshold value by using the CDF

$1-q=C_{M}\left(Z^{\ast}|\bm{\Sigma}\right).$ (70)

As $C_{M}\left(Z|\bm{\Sigma}\right)$ is a monotonically increasing function of
$Z$, $Z^{\ast}$ increases for decreasing rate $q$.

Even without explicit computations, we qualitatively expect that for a fixed
false-positive rate $q$,

  1. (i)

Threshold value $Z^{\ast}$ increases with the increasing size $M$ of template
bank.

  2. (ii)

Threshold value $Z^{\ast}$ decreases with the growing cross-template
correlations.

Indeed, the random data have more chance to match one of the templates in a
larger template bank leading to a higher threshold value, supporting statement
(i). As to the statement (ii), suppose we keep the bank size $M$ the same, but
increase one of the cross-template correlation coefficients $\Sigma_{ij}$ all
the way to 1. Then the two templates $i$ and $j$ become identical, effectively
reducing the bank to $M-1$ templates. Then we can apply the statement (i)
resulting in $Z^{\ast}$ being lowered. In particular, this means that for a
fixed rate of false positives, a bank of correlated templates has a lower SNR-
max threshold than a fully-orthogonal bank (uncorrelated or independent
templates.)

Evaluation of threshold values for a desired false positive rates for a large
bank with arbitrary bank covariance matrix $\bm{\Sigma}$ requires Monte-Carlo
integration of an $M$-dimensional multivariate distribution, see Sec. V. While
Monte-Carlo integration can yield the desired threshold value for any
$\bm{\Sigma}$ and dimensionality $M$, the computational cost can be
prohibitive when $M$ is large and when sampling far into the tail of the
distribution, as needed for small false positive rates. Therefore, in this
section we present several analytical results that avoid or greatly reduce
numerical costs.

For example, if the $M$ templates in the bank are fully orthogonal, i.e.
$\bm{\Sigma}=\bm{I}$, the analytic CDF is given by Eq. (38) and thresholds are

$Z^{*}=\sqrt{2}\,\mathrm{erf}^{-1}\left((1-q)^{1/M}\right)\,.$ (71)

This formula illustrates our qualitative statement (i). For $q\ll 1$, $Z^{*}$
is given by the conventional threshold formula for a Gaussian distribution,
$Z^{*}=\sqrt{2}\,\mathrm{erf}^{-1}\left(1-q_{\mathrm{eff}}\right)$ with
$q_{\mathrm{eff}}=q/M$. As $M$ is increased, $q_{\mathrm{eff}}$ decreases,
resulting in larger values of $Z^{*}$.

###  VI.1 Threshold for two-template bank

The discussion for the case of the two correlated template SNR-max CDF in Sec.
V.2 suggests that for a fixed value of $q$, the threshold value $Z^{\ast}$
moves to lower values the more correlated the templates are,

$Z^{\ast}\left(r=1\right)<Z^{\ast}\left(r\right)<Z^{\ast}\left(r=0\right).$
(72)

Fig. 11 shows the effect of correlations on the SNR-max CDF for the case of
two templates. Here we vary the correlation coefficient in the range $0\leq
r\leq 0.9$. It is evident from comparing the $r=0.75$ (red) and the $r=0.9$
(blue) curves that as the correlation increases, the CDF approaches its
limiting value of 1 quicker. Then the thresholds $Z^{*}$ for specific false
positive rates $q$ decrease with increasing cross-template correlation $r$,
which is in agreement with the qualitative statement (ii).

Figure 11: (Color online) SNR-max CDF for a bank of two templates with varying
correlation coefficient: $r=0$ (yellow curve), $r=0.75$ (red curve), and
$r=0.9$ (blue curve). The inflection points correspond to the maxima of PDFs.

Table 1 compiles the SNR-max threshold values for two templates with varying
correlations $r$: $0$, $0.25$, $0.5$, and $0.9$, and for varying false
positive rates $q$: $10^{-2}$, $10^{-4}$, $10^{-6}$, and $10^{-8}$.
Examination of the table shows that the highest threshold values occur for the
case of uncorrelated templates regardless of the false positive rate $q$. For
very small false positive rates $q$, the effect of correlation is nearly
negligible on the threshold value unless the correlation is very strong,
$r\approx 1$. Whereas increasing the desired false positive rate, the effect
of correlation on the threshold value becomes more pronounced even for mild
correlations $r\approx 0.25$. This effect of cross-template correlations on
thresholds becomes even more pronounced as we increase the number of templates
in the bank, see Sec. VI.3.

$q$ | $r=0$ | $r=0.25$ | $r=0.5$ | $r=0.9$
---|---|---|---|---
$10^{-8}$ | 5.848 | 5.848 | 5.848 | 5.833
$10^{-6}$ | 5.027 | 5.027 | 5.027 | 5.003
$10^{-4}$ | 4.056 | 4.056 | 4.054 | 4.014
$10^{-2}$ | 2.807 | 2.805 | 2.795 | 2.716
Table 1: Threshold values for a template bank consisting of $M=2$ templates.
Threshold values are tabulated for varying cross-template correlation
coefficients, $0\leq r\leq 0.9$ and false positive rates, $10^{-8}\leq q\leq
10^{-2}$.

###  VI.2 Threshold for a nearly-orthogonal template bank

Suppose we can represent the bank covariance matrix as a sum of the reference
matrix $\mathbf{\Sigma}_{0}$ and a small correction $\delta\mathbf{\Sigma}$.
This is the case for a nearly-orthogonal template bank of Sec. V.3, where
$\mathbf{\Sigma}_{0}=\mathbf{I}$, and for a “squeezed” bank of Sec. V.3,
$\mathbf{\Sigma}_{0}=\mathbf{\Sigma}_{r}$. Then the thresholds for specific
false positive rates can be calculated in a perturbative fashion. The dominant
term in calculating the threshold value is determined by the CDF
$C_{M}(Z^{*}|\mathbf{\Sigma}_{0})$ with a corrective term $\delta
C_{M}(Z^{*}|\mathbf{\Sigma}_{0})$,

$C_{M}(Z|\mathbf{\Sigma})=C_{M}(Z|\mathbf{\Sigma}_{0})+\delta
C_{M}(Z|\mathbf{\Sigma}_{0})\,.$ (73)

The threshold is determined from Eq. (70), $1-q=C_{M}(Z^{*}|\mathbf{\Sigma})$,
with $Z^{*}=Z_{0}^{*}+\delta Z^{*}$. Here $Z_{0}^{*}$ is the reference
threshold and $\delta Z^{*}$ is the correction to the threshold. The reference
threshold satisfies the equation $1-q=C_{M}(Z^{*}_{0}|\mathbf{\Sigma}_{0})$.
We would like to find the correction to the threshold $\delta Z^{*}$, so that

$1-q=C_{M}(Z^{*}_{0}+\delta Z^{*}|\mathbf{\Sigma}_{0})+\delta
C_{M}(Z^{*}_{0}|\mathbf{\Sigma}_{0})\,.$ (74)

Explicitly,

$\displaystyle C_{M}(Z^{*}_{0}+\delta Z^{*}|\mathbf{\Sigma}_{0})$
$\displaystyle\approx C_{M}(Z^{*}_{0}|\mathbf{\Sigma}_{0})+\delta
Z^{*}\times\frac{\partial C_{M}}{\partial Z}$
$\displaystyle=C_{M}(Z^{*}_{0}|\mathbf{\Sigma}_{0})+\delta Z^{*}\times
p_{M}(Z^{*}_{0}|\mathbf{\Sigma}_{0}).$

From here, the correction to the threshold is

$\delta Z^{*}\approx-\frac{\delta
C_{M}(Z^{*}_{0}|\mathbf{\Sigma}_{0})}{p_{M}(Z^{*}_{0}|\mathbf{\Sigma}_{0})}\,.$
(75)

Now we apply this formalism to the determination of the threshold for a
nearly-orthogonal template bank, i.e. $\bm{\Sigma}_{0}=\bm{I}$. Then
$Z^{*}_{0}=\sqrt{2}\,\mathrm{erf}^{-1}\left((1-q)^{1/M}\right)$, Eq. (71). The
corrective CDF term and the SNR-max PDF in Eq. (75) are per Sec. V.3,

$\displaystyle\delta C_{M}(Z_{0}^{*}|\bm{I})$
$\displaystyle=\mathrm{Tr}{\left(\bm{A}^{2}\right)}\frac{Z_{0}^{\ast
2}e^{-Z_{0}^{\ast
2}}}{2\pi}\left[\text{erf}\left(\frac{Z^{\ast}_{0}}{\sqrt{2}}\right)\right]^{M-2},$
$\displaystyle p_{M}(Z^{*}_{0}|\bm{I})$
$\displaystyle=M\sqrt{\frac{2}{\pi}}e^{-\frac{Z_{0}^{\ast
2}}{2}}\text{erf}\left(\frac{Z^{\ast}_{0}}{\sqrt{2}}\right)^{M-1}\,.$

Plugging these into Eq. (75), we arrive at the correction to the threshold
value,

$\delta
Z^{*}=-\frac{\mathrm{Tr}{\left(\bm{A}^{2}\right)}}{M}\left[\text{erf}\left(\frac{Z^{\ast}_{0}}{\sqrt{2}}\right)\right]^{-1}Z_{0}^{\ast
2}e^{-Z_{0}^{\ast 2}}\frac{1}{2\sqrt{2\pi}}.$ (76)

Since $\mathrm{Tr}{\left(\bm{A}^{2}\right)}>0$, we see that $\delta
Z^{\ast}<0$, so that $Z^{\ast}<Z_{0}^{\ast}$. This is consistent with our
qualitative expectation (ii) that cross-template correlations reduce the
threshold values for a given rate of false positives.

We remind the reader that for a nearly-orthogonal bank,
$\mathrm{Tr}{\left(\bm{A}^{2}\right)}=\mathrm{Tr}(\bm{\Sigma}^{2})-M$. The
advantage of our approximation (76) is that if the bank changes due to time
evolution of network geometry, e.g., for orbiting satellites, only
$\mathrm{Tr}(\bm{\Sigma}^{2})$ needs to be reevaluated to readjust the
threshold for the new bank.

###  VI.3 Threshold for squeezed template bank

Squeezed template banks introduced in Sec. V.4 have a special structure of the
bank covariance matrix, $\mathbf{\Sigma}\approx\mathbf{\Sigma}_{r}$, in which
all off-diagonal matrix elements have the same value of $r$ (mean value of
cross-template correlation coefficients in the bank). This is directly
relevant to the GPS.DM search for DM bubble walls. Analytic SNR-max CDF (67)
$C_{M}(Z|\mathbf{\Sigma}_{r})$ was developed in Sec. V.4. This representation
sharply reduces computational cost in determining the thresholds for specific
false positive rates, as the $M$-dimensional integral is recast into a one-
dimensional integral.

We used the SNR-max CDF (67) for squeezed banks and Eq. (70) to determine SNR-
max threshold values for desired rates of false positives. The results are
presented in Fig. 12 and Table 2. Fig. 12 shows the SNR-max CDFs for a bank of
$M=1024$ templates with cross-template correlation coefficient varying in the
range $0\leq r\leq 0.9$. It is clear that the effect of cross-template
correlation strongly affects the underlying CDF: for larger $r$, the CDF is
quickly shifted toward lower values of $Z$. The values of SNR-max threshold
are determined by the intersection of the computed CDFs with the horizontal
purple dashed line corresponding to a fixed false positive rate. Increasing
cross-template correlations causes decreasing thresholds $Z^{*}$ in agreement
with our qualitative expectation (ii) and discussion for $M=2$ banks. We can
understand this shift in the CDF and thresholds, by thinking about how the
inflection points (located at the most probable value) vary between the
different PDF’s associated with the correlation value $r$.

Numerical values for the SNR-max thresholds are compiled in Table 2. Once
again we observe that the threshold $Z^{*}$ decrease with increasing
correlation value $r$ for all the tabulated false positive rates. For $r=0.25$
and $r=0.5$ the decrease in the threshold is nearly negligible, while for
larger $r=0.9$ the decrease in the threshold becomes much more apparent,
similar to the case of $M=2$ template bank. We also see that as the false
positive rate decreases, the effect of correlation becomes more drastic on the
threshold values.

Figure 12: (Color online) SNR-max CDFs $C_{M}(Z|\mathbf{\Sigma}_{r})$ for a perfectly squeezed bank of $M=1024$ templates for various values of cross-template correlation coefficient: $r=0$ (blue curve), $r=0.5$ (red curve), and $r=0.9$ (yellow curve). The horizontal purple dashed line is drawn at a constant value of 0.9 corresponding to the rate of false positives $q=0.1$. Intersections of the CDF curves with the horizontal line determine the values of the threshold $Z^{\ast}$. Increasing cross-template correlations causes a decrease in thresholds values $Z^{*}$.  $q$ | $r=0$ | $r=0.25$ | $r=0.5$ | $r=0.75$ | $r=0.9$
---|---|---|---|---|---
$10^{-8}$ | 6.807 | 6.807 | 6.804 | 6.733 | 6.462
$10^{-6}$ | 6.109 | 6.109 | 6.097 | 5.972 | 5.705
$10^{-4}$ | 5.327 | 5.325 | 5.283 | 5.076 | 4.745
$10^{-2}$ | 4.417 | 4.390 | 4.247 | 3.903 | 3.487
Table 2: Compilation of SNR-max statistic threshold values $Z^{\ast}$ for for
a perfectly squeezed bank of $M=1024$ templates. $r$ cross-template
correlation coefficient and $q$ is false positive rate.

For completeness, we investigate a more practical case of how the thresholds
vary for template banks with template covariance matrix off-diagonal elements
not all equal to $r$, but there is a spread in correlation elements values. As
a specific example, we consider a template bank composed of $M=3$ templates
with bank covariance matrix

$\mathbf{\Sigma}=\left(\begin{array}[c]{ccc}1&r_{12}&r_{13}\\\
r_{12}&1&r_{23}\\\ r_{13}&r_{23}&1\end{array}\right)\,.$ (77)

We fix $r_{12}=0.33$, $r_{13}=0.23$, and $r_{23}=0.43$. This bank has a 0.2
spread in cross-template correlation coefficients about mean value of
$r=0.33$. In Table 3, we present results of two threshold computations: the
Markov-Chain Monte-Carlo (MCMC) method using publicly-available package [58]
for integration of the exact CDF (35), and approximate squeezed-bank CDF (67)
with $r=0.33$. For this bank, the two thresholds agree to three significant
figures.

$q$ | Exact, Eq. (35) | Approx., Eq. (67)
---|---|---
$10^{-8}$ | 5.913 | 5.915
$10^{-6}$ | 5.103 | 5.104
$10^{-4}$ | 4.149 | 4.150
Table 3: Comparison of threshold values from Eq. (35) and Eq. (67) for the
case of $M=3$ template bank (see Eq. (77) and text) and varying false positive
rates, $10^{-8}\leq q\leq 10^{-4}$.

As a next numerical experiment, we increase the number of templates to
$M=1024$ for a squeezed template bank, with a distribution of cross-template
correlations $\Sigma_{ij}$ shown in Fig. 7. In this bank, the most probable
value is 0.390 with a range of total correlation values being
$0.25\leq\Sigma_{ij}\leq 0.65$. The median value of the $\Sigma_{ij}$
distribution is 0.395, and an average value $r=0.398$. First, we compute the
SNR-max threshold values using the MCMC [58]. MCMC technique evaluates the
exact $M$-dimensional integral CDF (35); the results are presented in Table 4.
Second, we use the approximate CDF $C_{M}(Z|\mathbf{\Sigma}_{r})$ with the
most probable 0.390 and an average value of 0.398. The thresholds are
evaluated for the false positive rates $q=10^{-4}$, $q=10^{-6}$, and
$q=10^{-8}$.

Table 4 demonstrates that choosing the reference cross-template correlation
$r$ as the mean value of $\Sigma_{i\neq j}$ (versus most probable value)
results in closer agreement with the exact result. This is consistent with our
arguments presented in Sec. V.4. The approximation
$\mathbf{\Sigma}\approx\mathbf{\Sigma}_{r}$ works better for smaller values of
false positives. Even for $q=10^{-8}$, the approximation is in $2\%$ agreement
with the MCMC result, demonstrating the utility of our approximation.

$q$ | Exact | $r=0.390$ | $r=0.398$
---|---|---|---
$10^{-8}$ | 6.67 | 6.8063 | 6.8062
$10^{-6}$ | 6.04 | 6.1073 | 6.1070
$10^{-4}$ | 5.31 | 5.3122 | 5.3108
Table 4: Comparison of threshold values from the exact Eq. (35) (computed with
MCMC, column marked “exact”) and Eq. (67) (both for the most probable (middle
column) and average correlation matrix element value (last column)) for the
case of $M=1024$ templates (typical template bank vs. perfectly squeezed
template bank) and varying false positive rates, $10^{-8}\leq q\leq 10^{-4}$.

###  VI.4 Threshold for varying number of templates $M$ in template bank

SNR-max thresholds for specific false positive rates depend on the number of
templates in the template bank. Fig. 13 shows the effect of the number of
templates $M$ on the thresholds for varying size of perfectly-squeezed banks,
i.e. with bank covariance matrix (60). In all three panels of that figure, the
rate of false positives is fixed (horizontal line at 0.9). In all panels, we
plot the same CDFs $C_{M}(Z|\mathbf{\Sigma}_{r})$ for $M=100$ (blue curve) to
$M=500$ (red curve) to $M=1000$ (yellow curve). The value of the cross-
template coefficient $r$ decreases from $r=0.9$ (panel (a)) to 0.5 in panel
(b) and to 0.1 in panel (c). For all three panels, we observe that as the
number of templates increases from $M=100$ to $M=500$ and to $M=1000$ the
threshold (intersection of the horizontal purple dotted line with the CDF
curves) shifts to higher $Z$ values. As the correlation coefficient decreases,
panels (a) $\rightarrow$ (c), the differences (spread) in the threshold values
increase. These observations further support our qualitative statements (i)
and (ii): generically, we expect a larger threshold value with increasing
number of templates in the template bank, as well as an increasing threshold
value with decreasing cross-template correlations.

Figure 13: (Color online) SNR-max CDFs $C_{M}(Z|\mathbf{\Sigma}_{r})$ for a
varying number of templates in a perfectly-squeezed bank: $M=100$ (blue
curve), $M=500$ (red curve), and $M=1000$ (yellow curve). The false positive
rate is fixed at $q=0.1$ (horizontal purple dotted line is drawn at the
$1-q=0.9$ value). The thresholds are determined by the value of $Z$ at the
intersection of the CDF curves with the false positive rate horizontal line.
Three panels correspond to different values of the cross-template correlation
coefficient $r$: (a) $r=0.9$, (b) $r=0.5$, and (c) $r=0.1$.

##  VII Conclusions and the utility of entanglement

The signal-to-noise ratio (SNR) statistic and template bank searches have seen
a variety of applications in modern physics and other fields. Here we focused
on the max-SNR detection statistic in the context of searches for transient
signals with networks of precision quantum sensors. For large template banks,
especially for dynamically evolving networks, computation of threshold rates
is computationally prohibitive. We developed several approximate methods for
computing properties of the SNR-max statistic that forgo much of the
computational costs. We developed the probability and cumulative distribution
for the SNR-max statistic correlated from $M$ templates in a template bank.
Generically the cross-template correlations give rise to diminished thresholds
for desired false positive rates. While our paper was motivated by dark matter
and other exotic physics searches, the results of Secs. III—VI have a wide
range of applicability. We also anticipate that our detailed exposè on
template bank construction and optimization in the cases of the pedagogical
“toy” network and the GPS network (Sec. III) will be useful in the nascent
applications of precision quantum sensors and their networks in fundamental
physics.

What is the utility of entanglement in the searches for transients? For
individual sensors, the entanglement or spin-squeezing is a useful resource,
as it improves the accuracy of a single-shot measurement typical of a search
for short transients. Cross-node or geographically-distributed entanglement
[13] is not useful and is, in fact, detrimental to the network searches for
transients. Indeed, projective measurement on a single node collapses the
distributed wave-function, effectively rendering all other nodes deaf to the
transient. Then the network loses both velocity and angular resolution.

## Acknowledgments

We would like to thank M. Pospelov, J. Singh, and T. Tscherbul for
discussions. This work was supported in part by the U.S. National Science
Foundation grant PHY-1806672.

##  Appendix A Derivation of matrix elements of template bank covariance
matrix $\bm{\Sigma}$, Eq. (27)

Matrix elements of the template bank covariance matrix are defined as

$\Sigma_{ij}=\langle\rho_{i}^{w}\rho_{j}^{w}\rangle\,,$ (78)

where SNR statistics $\rho_{i}^{w}$ and $\rho_{j}^{w}$ are given by Eq. (3)
and the averaging is over signal-free data. From here on we drop the window
superscript $w$. Explicitly,

$\langle\rho_{i}\rho_{j}\rangle=\frac{\langle(\mathbf{d}^{T}\mathbf{E}^{-1}\mathbf{s}_{i})(\mathbf{d}^{T}\mathbf{E}^{-1}\mathbf{s}_{j})\rangle}{\sqrt{(\mathbf{s}_{i}^{T}\mathbf{E}^{-1}\mathbf{s}_{i})(\mathbf{s}_{j}^{T}\mathbf{E}^{-1}\mathbf{s}_{j})}}\,.$
(79)

To carry out the averaging, we expand the products in component form

$\mathbf{d}^{T}\mathbf{E}^{-1}\mathbf{s}_{i}=\sum_{\alpha,\beta}d_{\alpha}E_{\alpha\beta}^{-1}(s_{\beta})_{i},$
(80)

where we used Greek letters for compound indices that encode both the sensor
and the epoch indexes, see Sec. II. Then the numerator of Eq. (79) can be
rearranged into

$\displaystyle\left\langle\sum_{\alpha,\beta,\gamma,\delta}d_{\alpha}E_{\alpha\beta}^{-1}(s_{\beta})_{i}d_{\gamma}E_{\gamma\delta}^{-1}(s_{\delta})_{j}\right\rangle=$
$\displaystyle=\sum_{\alpha,\beta,\gamma,\delta}\left\langle
d_{\alpha}d_{\gamma}\right\rangle
E_{\alpha\beta}^{-1}E_{\gamma\delta}^{-1}(s_{\beta})_{i}(s_{\delta})_{j}\,.$
(81)

For signal-free data, the combination $\left\langle
d_{\alpha}d_{\gamma}\right\rangle$ is the matrix element of the noise
covariance matrix, $\left\langle
d_{\alpha}d_{\gamma}\right\rangle=E_{\alpha\gamma}\,.$ This simplifies Eq.
(81) to

$\displaystyle\sum_{\beta,\delta}(s_{\beta})_{i}(s_{\delta})_{j}\sum_{\alpha,\gamma}E_{\alpha\gamma}\left(E^{-1}\right)_{\gamma\delta}\left(E^{-1}\right)_{\alpha\beta}=$
(82)
$\displaystyle\sum_{\alpha,\beta}(s_{\alpha})_{j}\left(E^{-1}\right)_{\alpha\beta}(s_{\beta})_{i}=\mathbf{s}_{j}^{T}\mathbf{E}^{-1}\mathbf{s}_{i}.$

Here we used the identity,
$\sum_{\gamma}E_{\alpha\gamma}(E^{-1})_{\gamma\delta}=\bm{I}_{\alpha\delta}=\delta_{\alpha\delta}.$
With this simplification, the template bank covariance matrix element (79)
becomes

$\Sigma_{ij}=\frac{\mathbf{s}_{j}^{T}\mathbf{E}^{-1}\mathbf{s}_{i}}{\sqrt{\left(\mathbf{s}_{i}^{T}\mathbf{E}^{-1}\mathbf{s}_{i}\right)\left(\mathbf{s}_{j}^{T}\mathbf{E}^{-1}\mathbf{s}_{j}\right)}}\,.$
(83)

From here, the variances of individual SNR statistics are
$\sigma^{2}_{\rho_{i}}=\langle(\rho_{i})^{2}\rangle=\Sigma_{ii}=1$, in
agreement with Ref. [31]. This concludes our proof of Eq. (27).

##  Appendix B Template bank covariance matrix element construction

We are interested in the behavior of the template covariance matrix elements
$\Sigma_{ij}$ for our generalized signal template construction (Sec. III). To
streamline notation, we rename $\mathbf{s}_{i}$ and $\mathbf{s}_{j}$ to
$\mathbf{s}$ and $\mathbf{s}^{\prime}$ respectively. The template covariance
matrix element (27) in this notation reads

$\Sigma_{ss^{\prime}}=\frac{\mathbf{s}^{T}\mathbf{E}^{-1}\mathbf{s^{\prime}}}{\sqrt{\left(\mathbf{s}^{T}\mathbf{E}^{-1}\mathbf{s}\right)\left(\mathbf{s^{\prime}}^{T}\mathbf{E}^{-1}\mathbf{s^{\prime}}\right)}}.$
(84)

For concreteness, we focus on a white noise network of identical sensors $a$
referenced to a common reference sensor $R$. In this case, the inverse
covariance matrix is given by Eq. (11), which we reproduce below

$\left(E^{-1}\right)_{lm}^{ab}=A\,\delta_{lm}\delta^{ab}+B\,(1-\delta^{ab})\delta_{lm},$
(85)

where coefficients

$\displaystyle A$ $\displaystyle=\frac{1}{\sigma^{2}}\frac{1+(1-1/N_{\rm
S})\xi}{1+\xi}\,,$ (86) $\displaystyle B$
$\displaystyle=-\frac{1}{\sigma^{2}}\frac{1}{N_{\rm S}}\frac{\xi}{1+\xi}$ (87)

are expressed in terms of $\xi=N_{\rm S}\sigma_{R}^{2}/\sigma^{2}$. In our
index convention, the letters at the beginning of alphabet $a,b,\ldots$
enumerate the sensors (range $\overline{1,N_{\rm S}}$) , while letters in the
middle of alphabet $l,m,\ldots$ — the epochs (index range $\overline{1,J}$ in
a data window).

Since we consider the signal amplitude $h$ to be the same for the network and
reference sensors as described in Sec. III, we may simply take out $h$ from
the templates. Then the GPS.DM templates (thin DM walls) can be represented
using Kronecker symbols,

$s_{l}^{a}=\delta_{l,l_{a}}+(-1)\delta_{l,l_{R}}\,.$ (88)

The signal for each sensor $a$ consists of a positive spike at epoch $l_{a}$,
and a negative spike at the reference epoch $l_{R}$ which is identical for all
the templates, $l_{R}=(J+1)/2$. Then, the individual templates are
characterized by the set of network sensor perturbation epochs

$\\{l\\}={l_{1},l_{2},...,l_{N_{\rm S}}}\,.$ (89)

for the template $\bm{s}$ and a similar (but not identical) set
$\\{l^{\prime}\\}$ for the template $\bm{s}^{\prime}$.

In general, $l_{a}$ can be equal to $l_{R}$, meaning that the sensor $a$ and
the reference sensor are affected by the transient in the same epoch; in this
case the signal reduces to zero, $s_{l}^{a}=0$. In addition, several sensors
can be affected by the transient signal in the same epoch. In this case, we
refer to this subset of sensors as “degenerate”. When sensor $a$ is degenerate
with sensor $b$ for template $\bm{s}$, $l_{a}=l_{b}$ in the set $\\{l\\}$.

In our derivations below, we will use the fact that $l_{a}$ and $l_{R}$ are
guaranteed to be within the same template window, so that

$\sum_{l}\delta_{l,l_{a}}=\sum_{l}\delta_{l,l_{R}}=1.$ (90)

First, consider the numerator of Eq. (84),
$\mathbf{s}^{T}\mathbf{E}^{-1}\mathbf{s}^{\prime}$, which can be expanded into

$\displaystyle\mathbf{s}^{T}\mathbf{E}^{-1}\mathbf{s}^{\prime}$
$\displaystyle=$ (91) $\displaystyle A\sum_{a,b}^{N_{\rm
S}}\sum_{l,m}^{J}s^{a}_{l}\delta^{ab}\delta_{lm}\,s^{\prime
b}_{m}+B\sum_{a}^{N_{\rm S}}\sum_{b\neq a}^{N_{\rm
S}}\sum_{l,m}^{J}s^{a}_{l}\delta_{lm}s^{\prime b}_{m}.$

The first sum can be represented as

$\displaystyle\sum_{l,a}s^{a}_{l}s^{\prime a}_{l}$
$\displaystyle=\sum_{l,a}\left(\delta_{l,l_{a}}-\delta_{l,l_{R}}\right)\left(\delta_{l,l_{a}^{{}^{\prime}}}-\delta_{l,l_{R}}\right)=$
$\displaystyle\left(\sum_{l,a}\delta_{l,l_{a}}\delta_{l,l_{a}^{{}^{\prime}}}\right)+$
$\displaystyle\left(\sum_{l,a}\delta_{l,l_{R}}-\sum_{l,a}\delta_{l,l_{R}}\delta_{l,l_{a}}-\sum_{l,a}\delta_{l,l_{R}}\delta_{l,l_{a}^{\prime}}\right)$
$\displaystyle\equiv K_{1}^{s,s^{\prime}}+R_{1}^{s,s^{\prime}}.$

Here we introduced

$K_{1}^{s,s^{\prime}}\equiv\sum_{a}\delta_{l_{a},l_{a}^{{}^{\prime}}},$ (92)

which counts the number of sensors with a perturbation epoch identical to both
templates. If the two templates are identical,
$\mathbf{s}=\mathbf{s}^{\prime}$, then $K_{1}^{s,s}=N_{\rm S}$. Similarly,

$R_{1}^{s,s^{\prime}}\equiv N_{\rm
S}-\sum_{a}\left(\delta_{l_{R},l_{a}}+\delta_{l_{R},l_{a}^{{}^{\prime}}}\right)$
(93)

characterizes the network-reference sensor degeneracy. If both templates have
no degeneracies between network and the reference sensors, then
$R_{1}^{s,s^{\prime}}=N_{\rm S}$.

The second sum in Eq. (91) simplifies to

$\displaystyle\sum_{a}^{N_{\rm S}}\sum_{b\neq a}^{N_{\rm
S}}\sum_{l}^{J}s^{a}_{l}s^{\prime b}_{l}=$ (94) $\displaystyle\sum_{a}^{N_{\rm
S}}\sum_{b\neq a}^{N_{\rm
S}}\sum_{l}^{J}\left(\delta_{l,l_{a}}-\delta_{l,l_{R}}\right)\left(\delta_{l,l_{b}^{\prime}}-\delta_{l,l_{R}}\right)=$
$\displaystyle\sum_{a}^{N_{\rm S}}\sum_{b\neq a}^{N_{\rm
S}}\delta_{l_{a},l_{b}^{\prime}}+\sum_{a}^{N_{\rm S}}\sum_{b\neq a}^{N_{\rm
S}}(1)-\sum_{a}^{N_{\rm S}}\sum_{b\neq a}^{N_{\rm
S}}\delta_{l_{R},l_{a}}-\sum_{a}^{N_{\rm S}}\sum_{b\neq a}^{N_{\rm
S}}\delta_{l_{R},l_{b}^{\prime}}$ $\displaystyle\equiv
K_{2}^{s,s^{\prime}}+R_{2}^{s,s^{\prime}}\,,$

with

$\displaystyle K_{2}^{s,s^{\prime}}$ $\displaystyle\equiv\sum_{a}^{N_{\rm
S}}\sum_{b\neq a}^{N_{\rm S}}\delta_{l_{a},l_{b}^{{}^{\prime}}}\,,$ (95)
$\displaystyle R_{2}^{s,s^{\prime}}$ $\displaystyle\equiv N_{\rm S}(N_{\rm
S}-1)-\sum_{a}^{N_{\rm S}}\sum_{b\neq a}^{N_{\rm
S}}\left(\delta_{l_{R},l_{a}}+\delta_{l_{R},l_{b}^{{}^{\prime}}}\right)$ (96)
$\displaystyle=(N_{\rm S}-1)R_{1}^{ss^{\prime}}\,.$

$K_{2}^{s,s^{\prime}}$ defines the cross-template network degeneracy between
network sensors $a,b\neq a$. $R_{2}^{s,s^{\prime}}$ quantifies the cross-
template degeneracy between the reference and the network sensors. If there
are no degeneracies between the network sensors $a,b\neq a$ and the reference
sensor $R$, $R_{2}^{s,s^{\prime}}=N_{\rm S}(N_{\rm S}-1)$. For two identical
templates, $K_{2}^{s,s}\equiv\sum_{a}^{N_{\rm S}}\sum_{b\neq a}^{N_{\rm
S}}\delta_{l_{a},l_{b}}$. This definition of $K_{2}^{s,s}$ depends on the
degeneracy of the network with respect to the template. If, for instance, all
the sensors are affected by the transient at different epochs, $l_{a}\neq
l_{b}$ and $a\neq b$, $K_{2}^{s,s}=0$. Thus $K_{2}^{s,s}$ quantifies the
degeneracy of the network with respect to template $\mathbf{s}$. Further, if
all of the sensors are degenerate ($l_{a}=l_{b}$ for all $a,b$) then
$K_{2}^{s,s}=N_{\rm S}(N_{\rm S}-1)$.

Combining these results, we arrive at the expression for the
$\mathbf{s}^{T}\mathbf{E}^{-1}\mathbf{s^{\prime}}$ product in the numerator of
Eq. (84)

$\displaystyle\mathbf{s}^{T}\mathbf{E}^{-1}\mathbf{s^{\prime}}=AK_{1}^{s,s^{\prime}}+CR_{1}^{s,s^{\prime}}+BK_{2}^{s,s^{\prime}}\,,$
(97)

where $C=A+(N_{\rm S}-1)B=\frac{1}{\sigma^{2}}\left(\frac{1}{1+\xi}\right)$.
The products in the denominator of Eq. (84), such as
$\mathbf{s}^{T}\mathbf{E}^{-1}\mathbf{s}$, can be obtained from Eq. (97) with
substitution $s^{{}^{\prime}}\rightarrow s$.

With Eq. (97) we can express the elements of the bank covariance matrix (84)
as

$\displaystyle\Sigma_{ss^{\prime}}=$ (98)
$\displaystyle\frac{AK_{1}^{s,s^{\prime}}+CR_{1}^{s,s^{\prime}}+BK_{2}^{s,s^{\prime}}}{\sqrt{\left(AK_{1}^{s,s}+CR_{1}^{s,s}+BK_{2}^{s,s}\right)(AK_{1}^{s^{\prime},s^{\prime}}+CR_{1}^{s^{\prime},s^{\prime}}+BK_{2}^{s^{\prime},s^{\prime}})}}.$

This expression can be simplified further,

$\displaystyle\Sigma_{ss^{\prime}}=$ (99)
$\displaystyle\frac{aK_{1}^{s,s^{\prime}}+R_{1}^{s,s^{\prime}}+bK_{2}^{s,s^{\prime}}}{\sqrt{\left(aK_{1}^{s,s}+R_{1}^{s,s}+bK_{2}^{s,s}\right)\left(aK_{1}^{s^{\prime},s^{\prime}}+R_{1}^{s^{\prime},s^{\prime}}+bK_{2}^{s^{\prime},s^{\prime}}\right)}},$

with

$\displaystyle a$ $\displaystyle=1+(1-1/N_{\rm S})\xi\,,$ (100) $\displaystyle
b$ $\displaystyle=-\xi/N_{\rm S}\,.$ (101)

Eq. (B) is an exact expression for the thin planar transients, such as thin DM
walls.

Now we introduce a concept of bank-averaged cross-template correlation
coefficient $\overline{\Sigma}_{ss^{\prime}}$. This quantity is useful in
quantifying properties of the SNR-max statistic for squeezed banks of Sec.
V.4. We assume that the number of templates is large so that the bank covers
the accessible parameter space for the sweep parameters, discussed in Sec.
III. In addition, we assume that we can average over spatial positions of the
nodes, see the uniform-occupancy approximation for a spherically-distributed
network of Sec. III.2.

With these approximations, we determine the typical values
$\bar{K}_{1},\bar{K}_{2}$, and $\bar{R}_{1}$ of the degeneracy coefficients.
The key to computing these quantities is the most probable template spread,
Eq. (25)

$\overline{\Delta l}=2R/(v_{\perp}^{p}\Delta_{t})\,,$ (102)

where $v^{p}_{\perp}=209\,\mathrm{km/s}$ the most probable speed of the
encounter. For GPS.DM data $\overline{\Delta l}=8$, see Sec. III.2. Then, a
typical template extends over $\overline{\Delta l}$ epochs and perturbs
$N_{\rm S}/\overline{\Delta l}$ sensors each epoch. Leading to

$\displaystyle\bar{K}_{1}^{s,s}$
$\displaystyle=\bar{K}_{1}^{s^{\prime},s^{\prime}}=N_{\rm S},$
$\displaystyle\bar{K}_{1}^{s,s^{\prime}}$ $\displaystyle=N_{\rm
S}/\overline{\Delta l},$ (103) $\displaystyle\bar{K}_{2}^{s,s}$
$\displaystyle=\bar{K}_{2}^{s^{\prime},s^{\prime}}=\bar{K}_{2}^{s,s^{\prime}}={N_{\rm
S}(N_{\rm S}-1)}/{\overline{\Delta l}},$ $\displaystyle\bar{R}_{1}^{s,s}$
$\displaystyle=\bar{R}_{1}^{s^{\prime},s^{\prime}}=\bar{R}_{1}^{s,s^{\prime}}=\left(1-{2}/{\overline{\Delta
l}}\right)N_{\rm S}\,.$

With these values, the “bank-averaged” bank covariance matrix element reads

$\displaystyle\overline{\Sigma}_{ss^{\prime}}$
$\displaystyle=\frac{1}{2+\xi(1-1/N_{\rm S})}\,.$ (104)

For a network of white noise sensors without a reference sensor, such as the
GNOME network of magnetometers [39], the noise covariance matrix is determined
by Eq. (8). For this network, $a=\frac{1}{\sigma^{2}}$, $b=0$, and
$R_{1}^{s,s^{\prime}}=0$. Then the bank-averaged cross-template correlation
coefficient simplifies to

$\displaystyle\overline{\Sigma}_{ss^{\prime}}={1}/{\overline{\Delta l}}\,.$
(105)

##  Appendix C Derivation of SNR-max PDF, Eq. (35)

SNR-max CDF is given by an $M$-dimensional integral over a hyper-cube,

$C_{M}\left(Z|\bm{\Sigma}\right)=\int_{-Z}^{+Z}d\rho_{1}...\int_{-Z}^{+Z}d\rho_{M}f\left(\bm{\rho}|\mathbf{\Sigma}\right)\,,$
(106)

with joint probability distribution (26),

$f\left(\bm{\rho}|\mathbf{\Sigma}\right)=\frac{1}{\sqrt{\det(2\pi\bm{\Sigma}})}\exp\left(-\frac{1}{2}\bm{\rho}^{T}\bm{\Sigma}^{-1}\bm{\rho}\right)\,.$
(107)

SNR-max PDF, the focus of our derivation, is a derivative of CDF,

$p_{M}\left(Z|\mathbf{\Sigma}\right)=\frac{d}{dZ}C_{M}\left(Z|\mathbf{\Sigma}\right)\,.$
(108)

To compute this derivative, we apply the Leibniz rule,

$\displaystyle\frac{d}{dZ}\int_{a(Z)}^{b(Z)}g(u,Z)du=$ (109)
$\displaystyle=\int_{a(Z)}^{b(Z)}\frac{\partial g(u,Z)}{\partial
Z}du+g(b(Z),Z)\frac{db}{dZ}-g(a(Z),Z)\frac{da}{dZ}.$

As a preliminary step, we introduce partial integrals of the joint PDF

$I_{n}\left(\rho_{1},\ldots\rho_{n};Z\right)\equiv\left(\prod_{k=n+1}^{M}\int_{-Z}^{Z}d\rho_{k}\right)f\left(\bm{\rho}|\mathbf{\Sigma}\right)$
(110)

with boundary values $I_{0}(;Z)=C_{M}\left(Z|\mathbf{\Sigma}\right)$ and
$I_{M}(\bm{\rho};)=f\left(\bm{\rho}|\mathbf{\Sigma}\right)$. Apparently,

$I_{n}\left(\rho_{1},\ldots\rho_{n};Z\right)=\int_{-Z}^{Z}d\rho_{n+1}I_{n+1}\left(\rho_{1},\ldots\rho_{n+1};Z\right)\,.$
(111)

As we repeatedly apply the Leibniz rule to the differentiating the CDF (106),
we encounter the general structure

$\displaystyle\left({\prod\limits_{m=1}^{n}{\int\limits_{-Z}^{Z}{d{\rho_{m}}}}}\right)\frac{d}{{dZ}}\left({\prod\limits_{l=n+1}^{M}{\int\limits_{-Z}^{Z}{d{\rho_{l}}}}}\right)f\left(\bm{\rho}|\mathbf{\Sigma}\right)=$
$\displaystyle=\left({\prod\limits_{m=1}^{k}{\int\limits_{-Z}^{Z}{d{\rho_{m}}}}}\right)\frac{d}{{dZ}}{I_{n}\left(\rho_{1},\ldots\rho_{n};Z\right)},$

which can be evaluated as

$\begin{array}[]{l}\frac{d}{{dZ}}{I_{n}}\left({{\rho_{1}},\cdots,{\rho_{n}};Z}\right)=\\\
\frac{d}{{dZ}}\int_{-Z}^{Z}d{\rho_{n+1}}{I_{n+1}}\left({{\rho_{1}},\cdots,{\rho_{n+1}};Z}\right)\\\
=\int_{-Z}^{Z}d{\rho_{n+1}}\frac{d}{{dZ}}{I_{n+1}}\left({{\rho_{1}},\cdots,{\rho_{n+1}};Z}\right)+\\\
+{\left.{{I_{n+1}}\left({{\rho_{1}},\cdots,{\rho_{n+1}};Z}\right)}\right|_{{\rho_{n+1}}=Z}}+\\\
+{\left.{{I_{n+1}}\left({{\rho_{1}},\cdots,{\rho_{n+1}};Z}\right)}\right|_{{\rho_{n+1}}=-Z}}.\end{array}$

Applying this recursion relation to
$I_{0}(;Z)=C_{M}\left(Z|\mathbf{\Sigma}\right)$ we eventually reach
$I_{M}(\bm{\rho};)=f\left(\bm{\rho}|\mathbf{\Sigma}\right)$ for which the
derivative with respect to $Z$ vanishes. The result reads

$\displaystyle
p_{M}\left(Z|\bm{\Sigma}\right)=\frac{1}{\sqrt{\det\left(2\pi\bm{\Sigma}\right)}}\times$
(112)
$\displaystyle\sum_{k=1}^{M}[\left(\int_{-Z}^{+Z}{\prod\limits_{n=1,n\neq
k}^{M}}d\rho_{n}\right)\left[\exp\left(-\frac{1}{2}\bm{\rho}^{T}\bm{\Sigma}^{-1}\bm{\rho}\right)\right]_{\rho_{k}=+Z}$
$\displaystyle+\left(\rho_{k}=-Z\right)].$

The quadratic form $\bm{\rho}^{T}\bm{\Sigma}^{-1}\bm{\rho}$ is invariant under
$\bm{\rho}\rightarrow-\bm{\rho}$. Thus by changing the variables
$\bm{\rho}\rightarrow-\bm{\rho}$ in the second term we can show that it is
equal to the first term. Thereby, the SNR-max PDF reduces to

$\displaystyle
p_{M}\left(z|\bm{\Sigma}\right)=\frac{2}{\sqrt{\det\left(2\pi\bm{\Sigma}\right)}}\times$
(113) $\displaystyle\sum_{k=1}^{M}{\prod\limits_{n=1,n\neq
k}^{M}}\left(\int_{-z}^{+z}d\rho_{n}\right)\left[\exp\left(-\frac{1}{2}\bm{\rho}^{T}\bm{\Sigma}^{-1}\bm{\rho}\right)\right]_{\rho_{k}=+z}.$

This concludes our proof of Eq. (35) of the main text.

##  Appendix D Derivation of SNR-max CDF and PDF for a nearly-orthogonal bank,
Eqs. (56,V.3)

Suppose that the correlations between the various templates were small, i.e.
$|\Sigma_{ij}|\ll 1$ for $i\neq j$, such that the template covariance matrix
$\mathbf{\Sigma}$ can be thought of as the identity matrix $\bm{I}$ plus some
perturbation matrix $\bm{A}$: $\mathbf{\Sigma}=\bm{I}+\bm{A}$. In this case,
the bulk of the SNR-max CDF and PDF is determined by the distributions for
fully-orthogonal bank with corrections that depend on $\bm{A}$:

$\displaystyle C_{M}\left(Z|\bm{\Sigma}\right)\approx$ $\displaystyle
C_{M}\left(Z|\bm{I}\right)+\delta_{I}C_{M}\left(Z|\bm{\Sigma}\right)\,,$ (114)
$\displaystyle p_{M}\left(z|\bm{\Sigma}\right)\approx$ $\displaystyle
p_{M}\left(z|\bm{I}\right)+\delta_{I}p_{M}\left(z|\bm{\Sigma}\right)\,.$ (115)

Here we use the symbol $\delta_{I}$ to emphasize that the expansion is about
$\bm{\Sigma}=\bm{I}$. The leading terms are given by Eqs. (38, 39), reproduced
below,

$\displaystyle C_{M}\left(Z|\bm{I}\right)$
$\displaystyle=\left[\text{erf}\left(\frac{Z}{\sqrt{2}}\right)\right]^{M}\,,$
(116) $\displaystyle p_{M}\left(z|\bm{I}\right)$
$\displaystyle=M\sqrt{\frac{2}{\pi}}e^{-\frac{z^{2}}{2}}\left(\text{erf}\left(\frac{z}{\sqrt{2}}\right)\right)^{M-1}.$
(117)

The goal of this Appendix is to derive the corrections
$\delta_{I}C_{M}\left(Z|\bm{\Sigma}\right)$ and
$\delta_{I}p_{M}\left(z|\bm{\Sigma}\right)$.

The determinant of the template covariance matrix can be approximated as

$\det\left(\bm{I}+\bm{A}\right)\approx
1+\mathrm{Tr}{\left(\bm{A}\right)}+\frac{1}{2}\left(\mathrm{Tr}{\left(\bm{A}\right)}\right)^{2}-\frac{1}{2}\mathrm{Tr}{\left(\bm{A}^{2}\right)}+...$
(118)

Here the elements of the perturbative matrix $\bm{A}$ are

$A_{ij}=\left(1-\delta_{ij}\right)\Sigma_{ij}.$ (119)

Since the perturbation matrix $\bm{A}$ has zeros along the diagonal,
$\mathrm{Tr}{\left(\bm{A}\right)}=0$, and the determinant simplifies to

$\det\left(\bm{\Sigma}\right)=\det\left(\bm{I}+\bm{A}\right)\approx
1-\frac{1}{2}\mathrm{Tr}{\left(\bm{A}^{2}\right)}.$ (120)

Notice that $\mathrm{Tr}{\left(\bm{A}^{2}\right)}$ can be related to
$\mathrm{Tr}{\left(\bm{\Sigma}^{2}\right)}$. Indeed,
$\bm{\Sigma}^{2}=\bm{I}+2\bm{A}+\bm{A}^{2}$, so that

$\mathrm{Tr}{\left(\bm{A}^{2}\right)}=\mathrm{Tr}{\left(\bm{\Sigma}^{2}\right)}-M\,.$
(121)

Also since, by construction, $\bm{A}$ is a symmetric matrix,
$\mathrm{Tr}{\left(\bm{A}^{2}\right)}=\sum_{ij}\left(A_{ij}\right)^{2}>0$.

As for the inverse of the template covariance matrix, this can also be
approximated as

$\bm{\Sigma}^{-1}=\left(\bm{I}+\bm{A}\right)^{-1}\approx\bm{I}-\bm{A+A}^{2}.$
(122)

Now that we have our approximate inverse template covariance matrix and
determinant, we can simplify the SNR-max CDF (35),

$\displaystyle
C_{M}\left(Z|\bm{\Sigma}\right)=\frac{1}{\left(2\pi\right)^{M/2}\sqrt{\det\bm{\Sigma}}}$
(123)
$\displaystyle\int_{-Z}^{+Z}d\rho_{1}...\int_{-Z}^{+Z}d\rho_{M}\exp\left(-\frac{1}{2}\bm{\rho}^{T}\bm{\Sigma}^{-1}\bm{\rho}\right).$

Using expansion (122), we represent the exponential in $C_{M}$ as

$\displaystyle\exp\left(-\frac{1}{2}\bm{\rho}^{T}\bm{\Sigma}^{-1}\bm{\rho}\right)$
(124)
$\displaystyle\approx\exp\left(-\frac{1}{2}\bm{\rho}^{T}\bm{\rho}\right)\exp\left(\frac{1}{2}\bm{\rho}^{T}\bm{A\rho}\right)\exp\left(-\frac{1}{2}\bm{\rho}^{T}\bm{A}^{2}\bm{\rho}\right).$

Expanding out the second and the third exponentials up to $O(\bm{A}^{2})$ we
obtain

$\displaystyle\exp\left(-\frac{1}{2}\bm{\rho}^{T}\bm{\Sigma}^{-1}\bm{\rho}\right)$
(125)
$\displaystyle\approx\exp\left(-\frac{1}{2}\bm{\rho}^{T}\bm{\rho}\right)\times$
$\displaystyle\left(1+\frac{1}{2}\bm{\rho}^{T}\bm{A\rho+}\frac{1}{8}\left(\bm{\rho}^{T}\bm{A\rho}\right)^{2}-\frac{1}{2}\bm{\rho}^{T}\bm{A}^{2}\bm{\rho}\right).$

This expansion converges as long as $|\rho|\ll
1/\sqrt{\max\left|A_{ij}\right|}$.

Using these expansions, the correction to the SNR-max CDF reads

$\displaystyle\delta_{I}C_{M}\left(Z|\bm{\Sigma}\right)\approx$
$\displaystyle\frac{1}{4}\mathrm{Tr}{\left(\bm{A}^{2}\right)}C_{M}\left(Z|\bm{I}\right)$
(126) $\displaystyle+\frac{1}{2}\sum_{i\neq
j}A_{ij}\langle\rho_{i}\rho_{j}\rangle_{Z}$
$\displaystyle+\frac{1}{8}\sum_{i\neq j}\sum_{k\neq
l}A_{ij}A_{kl}\langle\left(\rho_{i}\rho_{j}\right)\left(\rho_{k}\rho_{l}\right)\rangle_{Z}$
$\displaystyle-\frac{1}{2}\sum_{ij}\left(\bm{A}^{2}\right)_{ij}\langle\rho_{i}\rho_{j}\rangle_{Z}\,.$

Here we kept all the terms up to $O(\bm{A}^{2})$ and introduced “averaging”
notation

$\langle
g(\bm{\rho})\rangle_{Z}\equiv\int_{-Z}^{+Z}d\rho_{1}...\int_{-Z}^{+Z}d\rho_{M}g(\bm{\rho})f\left(\bm{\rho}|\bm{I}\right)\,,$
(127)

where subscript $Z$ emphasizes that the $M$-dimensional integral is evaluated
over hyper-cube centered at $\bm{\rho}=0$ and of side $2Z$.