Towards a first design of a Newtonian-noise cancellation system for Advanced LIGO
M Coughlin,
1
N Mukund,
2
J Harms,
3, 4
J Driggers,
5
R Adhikari,
6
and S Mitra
2
1
Department of Physics, Harvard University, Cambridge, MA 02138, USA
2
Inter-University Centre for Astronomy and Astrophysics ,
Ganeshkhind, Pune University Campus Pune 411 007, India
3
Universit`a degli Studi di Urbino “Carlo Bo”, I-61029 Urbino, Italy
4
INFN, Sezione di Firenze, Firenze 50019, Italy
5
LIGO Hanford Observatory, Richland, WA, 99352, USA
6
Division of Physics, Math, and Astronomy, California Institute of Technology, MS 100-36, Pasadena, CA, 91125, USA
Newtonian gravitational noise from seismic fields is predicted to be a limiting noise source at
low frequency for second generation gravitational-wave detectors. Mitigation of this noise will be
achieved by Wiener filtering using arrays of seismometers deployed in the vicinity of all test masses.
In this work, we present optimized configurations of seismometer arrays using a variety of simplified
models of the seismic field based on seismic observations at LIGO Hanford. The model that best
fits the seismic measurements leads to noise reduction limited predominantly by seismometer self-
noise. A first simplified design of seismic arrays for Newtonian-noise cancellation at the LIGO sites
is presented, which suggests that it will be sufficient to monitor surface displacement inside the
buildings.
PACS numbers: 95.75.-z,04.30.-w
I. INTRODUCTION
With the recent detection of gravitational waves pro-
duced by a binary black-hole system [1], efforts to de-
velop novel technology to maximize the sensitivity of the
existing gravitational-wave detectors such as Advanced
LIGO [2] and Advanced Virgo [3] are being reinforced.
A possible upgrade of the advanced detectors with min-
imal impact on the detector infrastructures is the can-
cellation of so-called Newtonian noise (NN) produced by
terrestrial gravity fluctuations [4–6].
The idea of NN cancellation is to monitor the sources
of gravity perturbations, which are generally associated
with fluctuating mass density in the vicinity of the test
masses. For example, microphones can be used to re-
trieve information about certain density perturbations
in the atmosphere, and seismometers provide informa-
tion about density perturbations of the ground. Predic-
tions based on a detailed characterization of the LIGO
sites show that seismic surface fields give the dominant
contribution to NN [7]. In this case, a NN cancellation
scheme can be realized using an array of seismometers
deployed at the surface near the test masses [6].
In the absence of NN observations, studies of NN can-
cellation schemes rely on precise modelling. Models of
seismic NN have been gradually refined over the past
decades [6, 8–10]. Atmospheric NN was identified as a
limiting noise source in superconducting gravimeters at
mHz frequencies, and noise cancellation successfully im-
plemented using pressure sensors [11]. However, it is to
be expected that atmospheric NN cancellation in future
large-scale GW detectors will be significantly more com-
plicated due to a greater variety of atmospheric phenom-
ena that can cause NN [12].
The design of a seismic NN cancellation scheme for
the advanced detectors needs to address the optimization
of the filter used to calculate a NN estimate from seis-
mic data, and also the optimal placement of seismome-
ters around test masses [7]. The traditional cancella-
tion scheme is based on Wiener filters as was already
implemented in feed-forward configuration with respect
to other GW detector noise [13–15]. Wiener filters are
typically calculated from observed correlations between
sensors and target channels, but it is also possible to
employ models of correlations informed by observations.
The latter is advantageous if prior knowledge of the seis-
mic field can be used to constrain the correlation models
(such as seismic speeds, or the location of known seismic
sources).
Optimizing the shape of seismic arrays is important for
maximizing the efficacy of the noise cancellation. Drig-
gers et al. [7] first explored schemes to mitigate NN
by estimating and subtracting it from the interferometer
data stream. They simulated seismic fields with local and
non-stationary components, achieving subtraction levels
of about a factor of 10 with a finite-impulse response
implementation of the Wiener filters. In this paper, we
explore using measured seismic fields from an array pre-
viously stationed at LIGO Hanford to inform the models
for the seismic fields we use to perform array optimiza-
tion. We include four correlation models fit to the ob-
served data. For all of the models, we use the LHO ar-
ray data to inform parameters in these models, including
seismic speeds. At that point, we use high-dimensional
samplers computing optimal arrays for these models by
minimizing the expected noise residuals.
The structure of the paper is as follows. We review the
formalism of NN and the models used in this analysis in
section II. We discuss how previous measurements from
an array of accelerometers at the LHO site inform our
seismic field models in section III. In section IV, we de-
scribe how we perform the array optimization. In section
arXiv:1606.01716v1 [gr-qc] 6 Jun 2016
2
V, we provide the results of analytic and sampler opti-
mizations for reference seismic fields. Our conclusions
are summarized in section VI.
II. NEWTONIAN NOISE CANCELLATION
A. Plane-wave Rayleigh Newtonian noise
We briefly review Rayleigh-wave NN (see section 3.4.2
in [6] for details). Rayleigh waves are characterized by
an evanescent displacement field, which means that dis-
placement amplitude decreases exponentially with dis-
tance to the surface. Density perturbations caused by
Rayleigh waves can be produced in two ways. First, one
needs to consider density changes due to normal surface
displacement
ξ
z
(
~%,ω
) at frequency
ω
, where
~%
= (
x,y
) is
the horizontal coordinate vector. Second, Rayleigh waves
produce density changes inside the medium. We start
with gravity perturbations from pure surface displace-
ment, and focus on homogeneous media in a half-space
for simplicity.
If surface displacement is associated with a single
Rayleigh wave, then the corresponding gravity acceler-
ation along the horizontal direction
x
of a test mass at
location
~%
0
can be written
δa
x
(
~%
0
,ω
) =
−
2
π
i
Gρ
0
ξ
z
(
~%
0
,ω
) exp(
−
kh
)
k
x
/k,
(1)
where
G
is the gravitational constant,
ρ
0
the density
of the homogeneous medium,
k
the wavenumber,
h
the
height of the test mass above ground, and the wave vector
is
~
k
= (
k
x
,k
y
) with
|
~
k
|
=
k
. Displacement has the same
phase at points inside the medium sharing the same hor-
izontal coordinates. This means that additional contri-
butions to the gravity perturbation from density changes
inside the medium can be accounted for by multiplying
the last equation with a number
γ
(
ν
). It depends on the
Poisson’s ratio
ν
of the medium. For Rayleigh waves, the
value of this factor is about
γ
≈
0
.
8 [6]. This means that
NN from density perturbations of the medium partially
cancel NN from surface displacement.
B. Wiener filtering
Noise cancellation using Wiener filters exploits corre-
lations between reference sensors and a target channel to
provide a coherent estimate of certain noise contributions
to the target channel [16]. Wiener filters are the opti-
mal, linear noise-cancellation filters provided that data
from all channels are described by stationary, random
processes. In the following, we assume this to be the
case.
For seismic NN cancellation, the reference sensors are
the
N
seismometers, and their mutual correlations will
be denoted
C
SS
in this paper, which is understood to be
a
N
×
N
matrix containing the cross spectral densities of
normal surface displacement
C
(
ξ
z
;
~%
1
,~%
2
,ω
) between two
seismometers located at
~%
1
, ~%
2
. The matrix contains seis-
mometer instrumental noise on its diagonal. The target
channel is the output data of the GW detector, but for
simplicity, we only consider NN from a specific test mass.
In this case, the
N
seismometers are located in the vicin-
ity of this test mass. The NN spectral density of a test-
mass located at
~%
0
is denoted by
C
TT
≡
C
(
δa
x
;
~ρ
0
,ω
). It
is straight-forward to extend the analysis to a cancella-
tion of the total NN from all four test masses.
The correlation between NN and seismometers is de-
noted by
~
C
ST
, which is a
N
-component vector containing
the cross spectral densities
C
(
ξ
z
,δa
x
;
~%,~%
0
,ω
) between
normal surface displacement observed at
~%
and NN of a
test mass at
~%
0
. The Wiener filter in frequency space can
be written
~w
=
~
C
>
ST
·
C
−
1
SS
, where
C
−
1
SS
is the inverse ma-
trix of
C
SS
. The coherent NN estimate is now obtained by
multiplying the Wiener filter with the amplitudes of seis-
mic displacement observed by the
N
seismometers. The
matrix
C
SS
and the vector
~
C
ST
depend on seismometer
locations.
Subtracting the coherent NN estimate from the ac-
tual test-mass NN, an average relative noise residual
R
is achieved, which is determined by
R
= 1
−
~
C
>
ST
·
C
−
1
SS
·
~
C
ST
C
TT
.
(2)
The cancellation performance can be limited by seis-
mometer instrumental noise, or by the amount of infor-
mation extracted from the seismic field. We call the first
limitation
instrumental
, and the second
geometrical
. For
example, a seismometer array with seismometers having
distances of order kilometer to each other cannot be used
to cancel NN if the seismic wavelength is of order 100 m
independent of the seismometer instrumental noise (un-
less the seismic field is extremely simple, e. g. associated
with a single plane Rayleigh wave).
C. Array optimization
Given a fixed number of seismometers, the optimal ar-
ray is found by changing seismometer locations and min-
imizing the noise residual
R
. From eq. (2), one finds
that the lowest possible residual with
N
seismometers is
R
min
(
ω
) = 1
/
(
Nσ
(
ω
)
2
) with
σ
(
ω
) being the signal-to-
noise ratio of the seismometers assumed to be equal for
all sensors. This residual can only be achieved when ge-
ometrical limitations are overcome, which is not always
possible even for an infinite number of seismometers (see
fig. 6).
We could now construct specific plane-wave composi-
tions of the seismic field and use eq. (1) to evaluate and
minimize the noise residual in eq. (2). While this has
some merit for analytical understanding of the optimiza-
tion, the goal of this paper is to present a method to use
observable quantities of the seismic field to model NN
3
and its cancellation. Seismic fields can show great com-
plexity for example due to seismic scattering, or presence
of local sources, in which case the plane-wave model is
impractical. The key here is to realize that all correla-
tion functions required for calculating the noise residual
can be expressed in terms of the two-point spatial corre-
lations
C
(
ξ
z
;
~%
1
,~%
2
,ω
) of the Rayleigh field:
C
(
δa
x
;
~ρ
0
,ω
) = (2
πGρ
0
γ
(
ν
))
2
∫∫
d
2
k
(2
π
)
2
d
2
k
′
(2
π
)
2
S
(
ξ
z
;
~
k,
~
k
′
,ω
)
k
x
k
k
′
x
k
′
e
−
hk
e
−
hk
′
e
i
~%
0
·
(
~
k
−
~
k
′
)
,
C
(
δa
x
;
~%
0
,ω
) = (
Gρ
0
γ
(
ν
))
2
∫∫
d
2
%
d
2
%
′
C
(
ξ
z
;
~%,~%
′
,ω
)
K
(
~%,~%
0
)
K
(
~%
′
,~%
0
)
,
C
(
ξ
z
,δa
x
;
~%
0
,~%,ω
) =
−
2
π
i
Gρ
0
γ
(
ν
)
∫∫
d
2
k
(2
π
)
2
d
2
k
′
(2
π
)
2
S
(
ξ
z
;
~
k,
~
k
′
,ω
)
k
x
k
e
−
hk
e
i(
~%
0
·
~
k
−
~%
·
~
k
′
)
,
C
(
ξ
z
,δa
x
;
~%
0
,~%,ω
) =
Gρ
0
γ
(
ν
)
∫
d
2
%
′
C
(
ξ
z
;
~%,~%
′
,ω
)
K
(
~%
′
,~%
0
)
.
(3)
with
S
(
ξ
z
;
~
k,
~
k
′
,ω
)
≡
∫∫
d
2
%
d
2
%
′
C
(
ξ
z
;
~%,~%
′
,ω
)e
−
i(
~%
·
~
k
−
~%
′
·
~
k
′
)
K
(
~%
1
,~%
2
)
≡
x
1
−
x
2
(
h
2
+
|
~%
1
−
~%
2
|
2
)
3
/
2
.
(4)
Simplified versions of these equations for a homogeneous
field,
C
(
ξ
z
;
~%,~%
′
,ω
) =
C
(
ξ
z
;
~%
−
~%
′
,ω
), or a homogeneous
and isotropic field,
C
(
ξ
z
;
~%,~%
′
,ω
) =
C
(
ξ
z
;
|
~%
−
~%
′
|
,ω
), are
presented in [6].
The first and third line in eq. (3) follow from eq. (1)
by calculating expectation values of products of displace-
ment amplitudes at different wavenumbers (note that we
do not impose the condition
k
=
ω/c
, which holds for
waves propagating at speed
c
). The second and fourth
line provide the desired link between seismic observa-
tion and NN modelling. It is important that these equa-
tions only require seismic correlations and no other prior
knowledge of the seismic field such as seismic speeds. Of
course, the simplest analytical models of
C
(
ξ
z
;
~%,~%
′
,ω
)
that we can construct will have additional parameters,
and seismic wavelength or equivalently seismic speed may
be one of them. We are now prepared to either di-
rectly use measured correlations between seismometers,
or models fit to observed correlations, to optimize seis-
mometer arrays for NN cancellation. Consequently, the
path towards a NN cancellation scheme includes a de-
tailed site characterization to measure two-point spatial
correlations of the seismic field near all test masses.
III. LIGO HANFORD ARRAY
In 2012, from April through November, an array of
44 Wilcoxon Research 731-207 accelerometers [17] was
deployed at an end station building at LIGO Hanford.
Figure 1 shows the locations of accelerometers in the
vicinity of the vacuum enclosure. Each accelerometer’s
signal was conditioned, including amplification of a fac-
tor of 100, before being acquired digitally and saved at
1024 Hz sampling rate. The array configuration followed
as closely as possible the shape of a spiral since it has
favorable properties for analyzing seismic fields [18].
23.6 m
10.8 m
FIG. 1: Layout of the instrument floor at the Y-end station of
the LIGO Hanford Observatory. Green circles indicate place-
ment of accelerometers in the vicinity of the vacuum enclo-
sure. Dotted line indicates the spiral shape of the array. Due
to ongoing work at the time, the array is not centered around
the vacuum chamber that holds the test mass mirror.
Two analyses of the array data are presented in the fol-
lowing, which are directly relevant to NN modelling. The
first analysis is a measurement of Rayleigh-wave speed in
the frequency range 10 – 20 Hz. The method used here
is to decompose the seismic field into plane harmonics
(see section 3.6.3 of [6]), and collect the phase speeds as-
sociated with the maximum-amplitude component over a
period of a week. Figure 2 shows that mean seismic speed
at 10 Hz and 15 Hz is around 300 m/s consistent with pre-
vious site-characterization measurements [19], while at
20 Hz measured speeds around 380 m/s are higher than
observed previously. This anomalous dispersion can be
4
200
300
400
500
600
700
0
100
200
300
400
500
600
Speed [m/s]
10Hz
15Hz
20Hz
FIG. 2: Histograms of seismic speed measurements using
LHO array data at 10, 15, and 20 Hz.
explained by the concrete slab of the laboratory building,
which has greater effect at short seismic wavelengths. A
similar result can be expected for the LIGO Livingston
site, where Rayleigh-wave speeds in this frequency range
estimated from field measurements are around 200 m/s
[20]. While the Rayleigh-wave speed is not required to
evaluate the integrals in eq. (3), it can be used to define
simple correlation models as shown in section V A.
The speed histograms also indicate that seismic noise
is dominated by Rayleigh waves. The high-speed tail
of the distribution, which is potentially associated with
seismic body waves, has very small values. This is a
very important conclusion for NN modelling and greatly
impact plans for NN cancellation. It remains to be seen
of course, if the situation is similar at Livingston and
other stations of the Hanford detector.
The second analysis is a measurement of seismic corre-
lations
C
(
ξ
z
;
~%
1
,~%
2
,ω
) between all 29 accelerometers used
for this study. Only 29 out of 44 accelerometers were
used due to the others being mislabeled in the data aqui-
sition. Figure 3 shows values of
C
(
ξ
z
;
~%,ω
) at 15 Hz with
~%
=
~%
2
−
~%
1
and
~%
=
~%
1
−
~%
2
. If the seismic field were
homogeneous, then
C
(
ξ
z
;
~%
1
,~%
2
,ω
) =
C
(
ξ
z
;
~%
2
−
~%
1
,ω
),
and the function outlined by the scatter plot would be
smooth. However, there are parts in this plot where high
correlation values exist very close to low correlation val-
ues. So the observed
C
(
ξ
z
;
~%
2
−
~%
1
,ω
) is not smooth, and
therefore the seismic field is inhomogeneous. The most
likely cause of inhomogeneities are local sources, which
affect some seismometers of the array more strongly
than others altering correlations throughout the array.
Nonetheless, in this paper, we will only fit homogeneous
models to the scatter plot, assuming that we can under-
stand the result as a smoothly changing correlation func-
tion
C
(
ξ
z
;
~%
2
−
~%
1
,ω
) perturbed by a sufficiently small
number of outliers.
If inhomogeneities were more pronounced, then NN
modelling and array optimization must be directly based
on eq. (3) valid for inhomogeneous fields. Since it is very
−15
−10
−5
0
5
10
15
−10
−5
0
5
10
Position [m]
Position [m]
Coherence
−1
−0.5
0
0.5
1
FIG. 3: Measured correlation function
C
(
ξ
z
;
~%,ω
) for the
LIGO Hanford site at 15 Hz.
hard to conceive a model of
C
(
ξ
z
;
~%
1
,~%
2
,ω
) that could
be fit to observed inhomogeneous correlations, the bet-
ter strategy may be to calculate the integrals numeri-
cally over observed values, possibly using interpolation
techniques. These schemes were tested in preparation
of this article, but it was found that array optimization
based on numerical integration requires a larger num-
ber of seismometers. The crucial difference to optimiza-
tion in homogeneous fields is that with
N
seismometers,
for each seismometer, the function
C
(
ξ
z
;
~%
1
,~%
k
,ω
) with
k
= 1
,...,N
is only sampled
N
times, while an esti-
mate of the homogeneous correlation
C
(
ξ
z
;
~%
2
−
~%
1
,ω
) is
calculated from
N
×
N
samples. So given a certain pre-
cision goal, which is determined by the targeted noise-
suppression factor, a much larger number of seismome-
ters is required in the case of inhomogeneous fields com-
pared to homogeneous fields. Consequently, the level of
homogeneity of the seismic field is one of the most im-
portant properties concerning NN cancellation.
IV. ARRAY OPTIMIZATION
There are different methods to find the optimal array,
which are all based on minimizing eq. (2). The method
used in [6] first applies analytic transformations of the
equation to determine the optimal array as zeros of a new
set of equations. Alternatively, one can directly find the
minimal residual by applying generic high-dimensional
sampling algorithms such as nested sampling, Metropo-
lis Hastings MCMC, or particle swarm optimization (the
last was used in [7]).
The results in this paper were calculated with a
Metropolis Hastings MCMC algorithm implementing
adaptive simulated annealing (SA), which statistically
guarantees obtaining solutions close to global minima
[21, 22]. The cost function involves a parameter space
with dimensionality 2
N
(for
N
sensors) with numer-
ous local minima. The algorithm should optimize all
seismometer locations simultaneously to obtain near to
global solutions within a reasonable execution time. SA
5
is a probabilistic method whose name draws inspiration
from a metallurgic annealing process where heating fol-
lowed by slow cooling results in achieving minimum en-
ergy configuration in a crystal lattice . The success of the
optimization depends on careful selection of the probabil-
ity distribution
g
T
(
x
) used to span the multi-dimensional
parameter space, acceptance criteria
h
(∆
E
) are used to
select the newer points and the temperature or the an-
nealing schedule, which determine the level of vigorous-
ness associated with the search. Algorithm is constructed
to begin as a random exploration of the parameter space
and then progress with time towards a greedy search.
Usually the acceptance function is given by
h
(∆
E
) =
1
/
(1 + ∆
E/T
) where ∆
E
is the difference in cost func-
tion between the new and the previous state. The state
space probability distribution is given by a Cauchy dis-
tribution
g
(∆
x
) =
T/
(∆
x
2
+
T
2
)
(
D
+1)
/
2
. The expo-
nential annealing schedule is used to calculate temper-
ature at future time steps
T
i
(
k
) =
T
0
i
exp(
c
i
k
1
/D
) with
c
i
=
m
i
exp(
−
n
i
/D
). Here
m
i
and
n
i
act as free param-
eters that are used to fine tune the optimization process.
The variation in sensitivities of different parameters
during the course of the search is addressed by reanneal-
ing or rescaling the time-steps after every user-defined
number of accepted states. This essentially decreases
or increases the range for those parameters, which show
more or less sensitivity at the current annealing sched-
ule. Accuracy of the search after each SA iteration is
further enhanced by performing a polling based pattern
search, which is well suited for carrying out local min-
imization when the problem at hand is ill-defined and
non-differentiable. In our implementation, we perform
parallel analysis at each iteration by executing a MAT-
LAB implementation of adaptive SA on different com-
puting cores and selecting the best configuration for use
in further iterations. Here the annealing schedule and
the process of newer samples generation get adaptively
modified so as to achieve accelerated convergence.
For optimizing the array locations, the algorithm is
as follows. We will use either the analytic or measured
seismometer correlations, denoted
C
SS
. Arrays with
N
seismometers are then generated by the algorithms. The
C
SS
, and thereby the
C
SN
using eq. (3), is then com-
puted for the various cases of seismic correlations. This
process is repeated until the algorithm’s stopping crite-
rion is reached, and thereafter the array with the smallest
residual is taken to be the optimal array.
We make a number of assumptions in the following
analysis. We assume that the seismic field is stationary,
or in other words, the calculated arrays ensure best possi-
ble cancellation of the stationary background of the NN,
while there may exist better arrays when only higher per-
centiles of the NN are to be cancelled. We also assume a
signal-to-noise ratio (SNR) of 100 for the seismic sensors.
This may be a challenge in low-noise environments, but
at the existing sites of the LIGO and Virgo detectors,
the most sensitive instruments in the NN band such as
GS-13, STS-2, or T240 have SNRs of a few 1000. We
choose for the test mass to be suspended 1.5 m above
ground, which is approximately the height of the LIGO
test masses. It should be noted though that the test-
mass height has no effect on optimal arrays for plane-
wave Rayleigh fields. In our simulation, we search for
optimal sensor locations within a 100 m
×
100 m surface
area with a test mass at its center. This is conservative,
as this is larger than the area from which interesting NN
contributions are expected [23].
V. RESULTS
The structure of this section is as follows. We present
the models we use to fit to the measured LHO array seis-
mic field from section III in subsection V A. We then use
those models in subsection V B to find optimal seismic
arrays for each of the models presented.
A. Models
We solve the integrals of eq. (3) for various models of
the seismic correlation
C
SS
, including an isotropic plane-
wave (IPW) model, which takes the form of a Bessel
function [7], an isotropic Gaussian model, an anisotropic
plane-wave (APW) model, and the case of a single plane
wave (SPW).
Explicit expressions for
C
SS
exist for the IPW model,
C
SS
=
J
0
(
k
|
~r
i
−
~r
j
|
), the Gaussian model,
C
SS
=
exp(
−|
~r
i
−
~r
j
|
2
/σ
2
), and the SPW model,
C
SS
= cos(
~
k
·
(
~r
i
−
~r
j
)). In the APW case, the seismic correlations
are calculated numerically by averaging plane-wave con-
tributions over propagation directions with a direction
dependent Gaussian weight
A
(
φ
) = exp(
−
(
φ
−
φ
0
)
2
/δ
2
),
where
φ
is the polar angle quantifying the direction of a
Rayleigh wave.
Using the LHO array data presented in section III, we
perform fits of the above models to the measured
C
SS
.
Taking 300 m/s as the seismic speed, we find that at 15 Hz
then, the length of seismic waves is
λ
= 20 m, which
directly determines the wavenumbers of the IPW, SPW,
and APW models. For the anisotropic model, we find
best fit parameters of
φ
0
= 110
◦
and
δ
= 52
◦
.
The SPW and Gaussian models are not based on in-
dependent fits to observed correlations, but are based on
the above best-fit parameter values. The reason for this
is that both models are poor representations of observed
correlations, but serve as lower (SPW) and upper (Gaus-
sian) bound of achievable noise residuals.
The APW value of
φ
0
is also used for the propagation
direction of the SPW model, which, together with the
wavelength parameter from above, determines the wave
vector
~
k
. The Gaussian model parameter is defined as
σ
=
λ/π
, which yields
σ
= 6
.
5 m. In this way, the Gaus-
sian model is at least a good fit to observed short-range
correlations, while neglecting deviations at larger sensor
6
−15
−10
−5
0
5
10
15
−10
−5
0
5
10
Position [m]
Position [m]
Coherence
−1
−0.5
0
0.5
1
−15
−10
−5
0
5
10
15
−10
−5
0
5
10
Position [m]
Position [m]
Coherence
−1
−0.5
0
0.5
1
−15
−10
−5
0
5
10
15
−10
−5
0
5
10
Position [m]
Position [m]
Coherence
−1
−0.5
0
0.5
1
FIG. 4: From the left to right are the IPW, Gaussian, and APW models for the seismic fields.
separations. The Gaussian model is an isotropic, homo-
geneous correlation model, but it implies a significant
impact from local seismic sources, which can also include
scattering. There is no underlying plane-wave represen-
tation of the seismic field since this would necessarily give
rise to a IPW or APW type model. The IPW, APW,
and Gaussian models evaluated for the sensor pairs of
the LHO array are shown in figure 4.
B. Optimal Arrays
The structure of the optimal array analysis is as fol-
lows. For each of the models presented in section V A, we
use the optimization scheme presented in section IV to
minimize the NN residuals, given by eq. (2), at 15 Hz. We
perform the optimization for each of these cases for ar-
rays containing from 1 to 20 sensors, which allows for the
exploration of how many sensors are necessary to achieve
a required subtraction.
As discussed above, the optimization analysis is per-
formed minimizing residuals at 15 Hz, corresponding to
a wavelength of
λ
= 20 m assuming a seismic speed of
300 m/s, and a seismometer SNR of 100. On the left
of figure 5, we show optimized seismometer arrays with
six sensors for all correlation models. The key features
of all arrays optimized for a single test mass are that
they are symmetric with respect to the
x
-axis pointing
in the direction of the interferometer arm, with their dis-
tance to the test mass located at the center of the plot
proportional to the wavelength of interest. This optimal
distance depends also on the sensor SNR.
The IPW and APW models yield similar optimal ar-
rays with distances between seismometers and the test
mass being slightly shorter for the APW model. This
greatly facilitates the design of NN cancellation systems,
since anisotropies of the field play a minor role. For the
SPW model, there are no distinct solutions for optimal
seismometer locations. Instead, seismometers can be de-
ployed anywhere on the two dashed lines and lines par-
allel to this at larger distances separated by half a wave-
length. Also the Gaussian model produces a qualitatively
different optimal array, which is further evidence of the
fact that there is no underlying plane-wave representa-
tion of the seismic field.
For all correlation models, seismometer are placed
away from the test mass by a significant fraction of the
seismic wavelength. This is because a sensor located di-
rectly under the test mass has vanishing correlation with
NN, which means that a sensor not directly under but
very close to the test mass must have extremely high SNR
to show significant correlation with NN. Therefore, given
a finite SNR, optimal array configurations with any num-
ber of sensors place seismometers sufficiently far away
from the test mass according to their SNR, but not so
far that correlation with test-mass NN becomes too small
for effective cancellation.
Interestingly, in the case of the IPW and APW models,
the 6-sensor optimal arrays also represent the optimal
locations of arrays with a higher number of sensors. In
other words, the optimization algorithm yields solutions
where seismometers are placed on top of each other. In
practice, this can be realized by placing seismometers
as close to each other as possible, or to conclude that
instead of adding new seismometers to a 6-sensor array,
seismometers with higher SNR should be used. For the
same reason, the 6-sensor configuration is independent
of the seismometer SNR, which is not true for optimal
arrays with less than 6 sensors where seismometers with
higher SNR are positioned closer to the test mass.
In the right plot of figure 5, optimized arrays are again
shown for the IPW model. The optimal array at the cor-
ner station has a different configuration since seismome-
ters are used to cancel NN from both test masses simul-
taneously. The effect is greatest on placement of seis-
mometers that are close to both test masses. According
to these results, seismic arrays for NN cancellation are
located completely inside the buildings.
Figure 6 shows subtraction residuals. It is assumed
that seismometers have a frequency-independent SNR of
100. In the left plot, residual noise spectra relative to
the unsuppressed NN spectra are shown using 6-sensor
arrays in the case of the single test-mass cancellation,
and 10 sensors for the 2-mass cancellation. The Gaussian
residuals are consistently high over the entire frequency
range supporting its role as worst-case scenario. In all
spectra, residuals increase towards lower frequencies. It
is a consequence of the fact that the ability of an array
to extract information from the seismic field as relevant
to NN cancellation decreases if the separation between
sensors becomes significantly smaller than the seismic
wavelength. At high frequencies, wavelengths become
7
−20
−10
0
10
20
−25
−20
−15
−10
−5
0
5
10
15
20
25
X [m]
Y [m]
APW
IPW
SPW
Gaussian
FIG. 5: On the left are the locations of sensors resulting from numerically minimizing the subtraction residual for the IPW,
APW, SPW, and Gaussian correlation functions. The dashed lines indicate that optimal locations of seismometers for the
SPW can lie anywhere on these lines. On the right is an example arrangement of seismometers at both end stations and the
corner stations, optimized for cancellation at 15 Hz. The end station placement is based on the IPW model discussed in the
text. The corner station arrangement is also based on the IPW model, but this time simultaneously minimizing noise residuals
from both test masses.
smaller than distances between sensors decreasing corre-
lations between seismometers and NN, again degrading
the cancellation performance.
For the APW and IPW models with one or two test
masses, increasing sensor SNR and using the same seis-
mometer locations, residuals at frequencies below the
minimum at 15 Hz decrease accordingly, while the residu-
als above 15 Hz do not change significantly. Adding sen-
sors at larger distances to the test mass would not signif-
icantly change residuals above 15 Hz, but lead to further
modest decrease of residuals below 15 Hz. Adding sensors
at smaller distances to the test mass would not have a
great impact on residuals below 15 Hz, but would extend
the frequency band with good cancellation performance
to higher frequencies. In summary, if it is necessary to
improve residuals at frequencies below the minimal resid-
ual, then the best option is to use seismometers with
higher SNR. If it is necessary to lower residuals above
the minimum, then the best option is to add further seis-
mometers at shorter distance to the test mass. With
respect to the SPW model, adding sensors at smaller
and larger distances to the test mass or increasing SNR
improve residuals at all frequencies, and potentially elim-
inate the oscillatory behavior of SPW residuals at higher
frequencies. For the Gaussian model, neither adding sen-
sors nor increasing the SNR has a significant effect on
noise residuals, unless the extended array is reoptimized
for the new SNR or number of seismometers.
Scaling the seismometer locations relative to the test
mass by a constant factor, the residual noise spectra of
the plane-wave models with a single test mass are sim-
ply shifted towards higher or lower frequencies otherwise
maintaining absolute values and shapes. This is not the
case for the Gaussian model where NN residuals depend
on the test-mass height (this dependence cancels out in
all plane-wave models), and for the 2-mass IPW model
where the distance between test masses introduces an-
other reference length.
The right plot in figure 6 shows noise residuals at
15 Hz as a function of number of sensors. The residu-
als can be compared with the absolute sensor-noise limit
1
/
(SNR
√
N
) shown as dashed line, which marks the min-
imal residuals that can possibly be achieved with
N
sen-
sors independent of the seismic correlation model. Al-
ready with 2 or 3 seismometers, residuals are strongly
reduced for all plane-wave models. Beyond
N
= 4 sen-
sors, the IPW, APW, and SPW curves approximately
follow an asymptote proportional to 1
/
√
N
. In other
words, 4 sensors are already able to overcome all geomet-
ric limitations so that residuals are purely sensor-noise
limited. Nonetheless, the absolute sensor noise limit of
1
/
(SNR
√
N
) is not reached by the IPW and APW mod-
els independent of the number of sensors. In the case of
the IPW, residuals for large numbers of sensors lie about
a factor
√
2 above the absolute sensor-noise limit.
In summary, depending on the properties of the seis-
mic field, potential subtraction factors using 6-sensor ar-
rays span 3.5 – 250, with the APW model representing
the best fit to the observed data and corresponding to a
suppression by 140. Since inhomogeneities of the seismic