text:
# Data conditioning for gravitational wave detectors: A Kalman filter for regressing suspension violin modes

## I Introduction

Within the next several years five kilometer-scale interferometric gravitational wave detectors (GEO , LIGO  and VIRGO ) will begin searching for astronomical sources of gravitational waves, something the smaller TAMA  prototype detector has already begun. The fundamental limitation on the sensitivity of these terrestrial detectors is seismic noise at low frequencies, thermal excitations of the mirror substrates and their suspensions at intermediate frequencies, and laser shot noise at higher frequencies. Here we focus on that part of the suspension noise associated with the violin modes of the suspension wires.

Gravitational waves interact with interferometric detectors by changing the separation, and thus the round-trip light travel time, between the mirrors in the detector arms. In order to be sensitive to passing gravitational radiation, the mirrors in terrestrial interferometeric detectors are suspended as pendula in order that they be approximately inertially free above the pendulumn frequency. Nevertheless, thermal and other excitations of the suspension wires move the mirrors and lead to (narrow-band) data stream artifacts, located at the violin-mode resonance frequencies of the suspension wires. These narrow-band features, which stand some 40–60 dB above the shot noise background in the LIGO I detector, are signals themselves that, unless properly treated, can confound efforts to detect actual gravitational wave signals.

In the waveband of the detectors now under development the transfer function, which determines how the detector output at a particular frequency is related to incident radiation at that frequency, is constant. The time evolution of the signal arising from any incident gravitational waves is thus determined solely by the waves themselves. The time evolution of the signal arising from suspension violin mode excitations, on the other hand, is strongly influenced by the the dynamics of the modes. We make use of that distinction to design a Kalman filter  that can, in real-time, detect violin mode excitations — thermal or otherwise — as distinct from actual gravitational wave signals. The measured excitation can be used to subtractively remove the modes from the detector output, significantly reducing the rms detector noise and simplifying the down-stream analysis. Additionally, these measured excitations can be used to diagnose the instrument’s behavior.

Other techniques have been described for the identification and removal of resonant line features in gravitational wave detector data. In particular, Sintes and Schutz  have discussed methods suitable for addressing nearly periodic artifacts that appear in a range of harmonics and whose phases, in each harmonic, are linearly related. These methods are not suitable for isolating the violin modes, however: though a full spectrum of violin mode harmonics do appear in the detector output, these modes are all independent and independently excited. Allen and Ottewill  have applied the multi-taper spectral method of Lees and Park to estimate the contributions of these resonant features to the power spectrum. This non-parametric method does not and cannot provide a time-domain estimate of the mode state or its instantaneous contribution to the detector output and so cannot be used for conditioning the data prior to analysis.

In section II we summarize the theory of Kalman Filtering, laying the ground for its application to identifying gravitational wave detector suspension violin modes. The construction of the Kalman filter for a particular problem requires a dynamical model for the underlying physical process giving rise to the signal. We discuss our model and how it comes to be represented as a set of state-space difference equations on the physical observables. In section V we investigate the performance of the filter using both simulated data sets, whose “signal” and “noise” content is known exactly, and experimental data taken at the LIGO 40 M prototype detector during November 1994.***Preliminary results are presented in . We summarize our conclusions in section VI.

## II The Kalman Filter

### A Introduction

A Kalman filter is a mechanism for predicting the multi-dimensional state of a dynamical system from a multi-dimensional observable. The observable is assumed to be linearly related to the system’s state, which is assumed to evolve linearly with time. The observable is assumed to be contaminated by an additive Gaussian measurement noise of known co-variance. The dynamical system described by the state may be driven by a known input in addition to Gaussian process noise of known co-variance.

For example, the system may be a pendulum, the state may be the pendulum bob’s instantaneous position and momentum, and the observable may be the bob’s position. The process noise is the thermal fluctuation in the wire and the measurement noise in this case may be the photon shot noise background.

When the process and measurement noise are Normal (i.e., Gaussian) the Kalman filter gives the optimal estimate of the system’s state, in the least-squares sense. When the noise is more complex than Gaussian the Kalman filter gives the best estimate of the state assuming that only the second order statistical characterization of the noise is known. Alternatively, when the evolution equations are not linear or the observables are not related linearly to the state, then the Kalman filter has a natural generalization that makes it the optimal means of predicting the state from the observables in a least-squares sense .

The derivation of the Kalman Filter equations is straightforward and instructive.For a more detailed discussion of Kalman Filters, see . In the particular application dealt with in this paper we assume that there is no driving force (aside from the process noise) and that the system state and observables are needed only at discrete times; correspondingly, we limit our derivation to that case. Denote the state of the system at sample time $`k`$ by $`𝝍[k]`$ and the observable at sample time $`k`$ by $`𝒛[k]`$:

$`𝝍[k]`$ $`:=`$ $`\left(\begin{array}{c}\text{The state vector of the dynamical}\hfill \\ \text{system at sample }k\text{; a }M\times 1\text{ vector}\hfill \end{array}\right)`$ (4)
$`𝒛[k]`$ $`:=`$ $`\left(\begin{array}{c}\text{The system observable at}\hfill \\ \text{sample }k\text{; a }N\times 1\text{ vector}\hfill \end{array}\right).`$ (7)
The state $`𝝍`$ evolves from one sample to the next according to the linear equation

$$𝝍[k+1]=𝑨𝝍[k]+𝒘[k],$$
(8)
where $`𝑨`$ is a linear operator and $`𝒘[k]`$ is an instance of process noise:

$`𝑨`$ $`:=`$ $`\left(\begin{array}{c}\text{Linear transformation that}\hfill \\ \text{takes }𝝍[k]\text{ to }𝝍[k+1]\text{;}\hfill \\ \text{}M\times M\text{ matrix}\hfill \end{array}\right)`$ (12)
$`𝒘`$ $`:=`$ $`\left(\begin{array}{c}\text{Process noise: a }M\text{-dimensional}\hfill \\ \text{Normally distributed random}\hfill \\ \text{variable; a }M\times 1\text{ vector.}\hfill \end{array}\right)`$ (16)
The process noise is assumed to be zero mean and white. It is then fully characterized by its co-variance matrix $`𝑾`$:

$$𝑾:=E[𝒘𝒘^T],$$
(17)
where $`E[]`$ represents an ensemble average across instances of its argument. The process noise co-variance is a symmetric $`M\times M`$ constant matrix.

The observable $`𝒛`$ at sample $`k`$ is related to the system state by the linear equation

$$𝒛[k]=𝑪𝝍[k]+𝒗[k],$$
(18)
where $`𝑪`$ is a linear operator and $`𝒗[k]`$ is an instance of measurement noise:

$`𝑪`$ $`:=`$ $`\left(\begin{array}{c}\text{The linear transformation}\hfill \\ \text{that relates the state }𝝍[k]\hfill \\ \text{to the observable }𝒛[k]\text{; an}\hfill \\ N\times M\text{ matrix}\hfill \end{array}\right)`$ (23)
$`𝒗[k]`$ $`:=`$ $`\left(\begin{array}{c}\text{Measurement noise: a }N\text{-dimensional}\hfill \\ \text{Normally distributed random variable;}\hfill \\ \text{}N\times 1\text{ vector}\hfill \end{array}\right).`$ (27)
The measurement noise is assumed to be zero mean, white,If the measurement noise is not white, it can be whitened by a suitable, linear redefinition of the observable. and have co-variance $`𝑽`$:

$$𝑽:=E[𝒗𝒗^T],$$
(28)
where $`E[]`$ represents an ensemble average across instances of $`𝒗[k]`$ of the measurement noise. The measurement noise co-variance is a symmetric $`N\times N`$ matrix.

Finally, we assume that the measurement noise and process noise are independent:

$$E[𝒗𝒘^T]=0.$$
(29)
### B Predicting the state from observations

Suppose that we have an estimate $`\widehat{𝝍}[k1]`$ of the process state at sample $`k1`$. Denote the error in this estimate by $`\widehat{𝒆}[k1]`$:

$$\widehat{𝒆}[k1]:=𝝍[k1]\widehat{𝝍}[k1].$$
(30)
This error is a random variable whose second moments we denote $`𝑷[k1]`$:

$$\widehat{𝑷}[k1]:=E[\widehat{𝒆}[k1]\widehat{𝒆}^T[k1]].$$
(31)
From the estimated state at sample $`k1`$, we can form an a priori estimate of the state at sample $`k`$ using the state transition operator $`𝑨`$:

$$\stackrel{~}{𝝍}[k]:=𝑨\widehat{𝝍}[k1].$$
(32)
Associated with $`\stackrel{~}{𝝍}[k]`$ is an error $`\stackrel{~}{𝒆}[k]`$ whose second moments we denote $`\stackrel{~}{𝑷}[k]`$:

$`\stackrel{~}{𝒆}[k]`$ $`:=`$ $`𝝍[k]\stackrel{~}{𝝍}[k].`$ (34)
$`\stackrel{~}{𝑷}[k]`$ $`:=`$ $`E[\stackrel{~}{𝒆}[k]\stackrel{~}{𝒆}^T[k]]`$ (35)
$`=`$ $`𝑨\widehat{P}[k1]𝑨^T+𝑾.`$ (36)
The final equality follows from equations 34, 8 and 17.

From $`\stackrel{~}{\psi }[k]`$, we can make a prediction $`\stackrel{~}{𝒛}`$ of the observable $`𝒛[k]`$:

$$\stackrel{~}{𝒛}[k]=𝑪\stackrel{~}{𝝍}[k].$$
(37)
Now compare the estimate $`\stackrel{~}{𝒛}[k]`$ to the actual measurement $`𝒛[k]`$. These will in general be different. The Kalman Filter uses the difference between the the actual and estimated measurement to generate a corrected estimate $`\widehat{𝝍}`$ of the state:

$$\widehat{𝝍}[k]=\stackrel{~}{𝝍}[k]+𝑲[k]\left(𝒛[k]\stackrel{~}{𝒛}[k]\right).$$
(38)
The linear operator $`𝑲[k]`$ is the Kalman gain, which we calculate in section II C. The refined state estimate $`\widehat{𝝍}[k]`$ has an associated error $`\widehat{𝒆}[k]`$ (cf. eq. 30) and covariance $`\widehat{𝑷}[k]`$ (cf. eq. 31):

$$\widehat{𝑷}=\left(𝑰𝑲𝑪\right)\stackrel{~}{𝑷}\left(𝑰𝑲𝑪\right)^T,$$
(39)
where $`𝑰`$ is the identity matrix and we have used 30, 34 and 38 to represent $`\widehat{𝑷}`$ in terms of the unknown $`𝑲`$ and the known $`𝑪`$ and $`\stackrel{~}{𝑷}`$, all at sample $`k`$.

### C Calculating the Kalman Gain

Our object is to find the Kalman gain $`𝑲`$ that minimizes an appropriate measure of the error $`\widehat{𝑷}`$ in the state estimate $`\widehat{\psi }`$ at each sample $`k`$.

Recall that $`\widehat{𝑷}`$ is a co-variance matrix of the error $`\widehat{𝒆}`$; consequently, its eigenvectors describe linear combination of the state vector elements whose errors are uncorrelated. Correspondingly, the sum of the eigenvalues — the trace of $`\widehat{𝑷}`$ — measures the variance of the total error in the state estimate. Adopting $`\text{tr}\widehat{𝑷}`$ as a measure of the error in the state estimate, we choose the Kalman gain $`𝑲`$ to minimize $`\text{tr}\widehat{𝑷}`$:

$`𝑲[k]`$ $`=`$ $`\stackrel{~}{𝑷}[k]𝑪^T/\left(𝑽+𝑪\stackrel{~}{𝑷}[k]𝑪^T\right),`$ (40)
where we have adopted the notation

$$𝑸=𝑺/𝑹$$
(42)
to denote $`𝑸`$ as a solution to the linear equation

$$𝑺=𝑸𝑹$$
(43)
for known $`𝑹`$, $`𝑺`$. Correspondingly, equation 39 for $`\widehat{𝑷}`$ can be simplified:

$`\widehat{𝑷}`$ $`=`$ $`\left(𝑰𝑲𝑪\right)\stackrel{~}{𝑷}\left(𝑰𝑲𝑪\right)^T`$ (44)
$`=`$ $`\left(𝑰𝑲𝑪\right)\stackrel{~}{𝑷}\left(𝑰𝑲𝑪\right)𝑲\left(𝑽+𝑪𝑲𝑪^T\right)𝑲^T+𝑲𝑽𝑲`$ (45)
$`=`$ $`\left(𝑰𝑲𝑪\right)\stackrel{~}{𝑷}.`$ (46)
### D Summary

To summarize, we begin with an estimate of the state $`\widehat{𝝍}`$ and the associated error $`\widehat{𝑷}[k]`$ at sample $`k1`$. From these, we

1. Form an a priori estimate of the state $`\stackrel{~}{𝝍}[k]`$, the error associated with that estimate $`\stackrel{~}{𝑷}[k]`$, the predicted measurement $`\stackrel{~}{z}[k]`$, and the Kalman gain $`K[k]`$, all at sample $`k`$ (eqs. 32, 36, 37, and 40);
2. Observe the system to determine $`𝒛[k]`$;
3. Form corrected estimates of the state $`\widehat{\psi }[k]`$ and the estimated error $`\widehat{𝑷}[k]`$ (cf. eq. 38, 46).

We end as we begin, ready to make a new estimate of the state at sample $`k+1`$. Once we are started, we can continue as long as observations can be made. To begin, we need only supply an initial estimate of the state and the error. These initial guesses are not critical as the algorithm continually corrects the estimated state and re-estimates the error with each iteration, converging on an optimum (in the sense of minimum $`\text{tr}\widehat{𝑷}`$) estimate of the state.

The cost per iteration — in either computation or memory — is constant. Only the present estimates of the state and the associated error need be retained from step to step: all of the history necessary to estimate the future values of the state and the associated error is encoded in these quantities.

## III Application: Suspension violin modes

### A Introduction

The cavity mirrors of terrestrial gravitational wave detector are suspended as pendula in order that they be inertially free at frequencies above the pendulum frequency. The suspension wires are under tension from the mirror mass and have a spectrum of “violin” modes. Each violin mode is, like the pendulum mode, endowed with $`k_BT`$ of energy from its contact with a thermal bath. The corresponding motion of the cavity mirrors leads to noise of thermal origin in the detector output. Other disturbances affect the suspension, increasing the overall noise power in the violin modes.

To minimize the impact of the thermal noise on the detector’s performance the suspensions are designed to have high $`Q`$. Laser shot noise dominates the detector noise power in a broad band about the violin mode resonant frequencies: only the peaks of the stochastically excited violin modes stand above this background. The peaks are strong, however: in the LIGO 40 M prototype they stand approximately 50 dB above the shot-noise floor and the contrast in the LIGO I detector will be approximately the same. In the LIGO 40 M prototype the fundamental violin mode resonant frequencies are in the (570, 606) Hz band with $`Q`$’s on order several times $`10^4`$, while in the LIGO I interferometer the fundamental violin modes will have frequencies near 330 Hz and $`Q`$’s on order $`10^5`$ . Figure 1 shows the power spectral density (PSD) of the differential-mode read-out (IFO\_DMRO) of the LIGO 40 M prototype interferometer in the neighborhood of the violin mode resonances for a typical stretch of data taken during November 1994. The units of the PSD are ADC-counts2/Hz. The line at 600 Hz is the ninth harmonic of the omnipresent 60 Hz power-main. Note how the violin mode resonances stand-out as sharp peaks some 40 dB above an approximately white (over this bandwidth) background, as described.

The combination of violin modes in an approximately white background maps naturally onto a model involving a stochastically driven dynamical process observed in the presence of white measurement noise:

* Each violin mode corresponds to a distinct dynamical processes driven by an independent noise source (thermal or otherwise);
* The observable is the motion of the cavity mirror as reflected in the gravity-wave channel of the interferometer;
* The measurement noise is the laser shot-noise, and other technical noises that affect the gravity-wave channel but do not drive the violin modes.

In the remainder of this section we describe how we restrict attention to the narrow band surrounding the violin mode, where the model associated with the use of a Kalman filter is particularly apt, and develop the state-space equations that describe the modes in that band.

### B The Observations

The violin mode resonances are very weakly damped; so, the noise power in the modes contributes significantly to the overall detector output only in a very narrow band about each mode’s resonant frequency. Correspondingly, to monitor one or more violin modes we focus attention on a narrow band in the the neighborhood of the mode resonant frequencies. For definiteness, denote the gravitational wave channel at full bandwidth (sample rate $`f_s`$) by $`g`$ and let the band of interest, which includes the mode frequency $`f_0`$, range from $`f_c\mathrm{\Delta }f/2`$ to $`f_c+\mathrm{\Delta }f/2`$, for $`\mathrm{\Delta }f`$ much less than $`f_c`$.

To exploit the narrow bandwidth of the violin modes we form the complex quantity

$$g^{}[k]=\mathrm{exp}\left(2\pi ik\frac{f_c}{f_s}\right)g[k].$$
(47)
The spectrum of $`g^{}`$ in the band $`[\mathrm{\Delta }f/2,\mathrm{\Delta }f/2]`$ is equal to the spectrum of $`g`$ in the band $`[f_c\mathrm{\Delta }f/2,f_c+\mathrm{\Delta }f/2]`$. Focus attention on just this narrow band by down-sampling $`g^{}`$ to the new sample rate $`2\mathrm{\Delta }f`$ using, e.g., a polyphase technique \[16, Chapter 4\].§§§The down-sampling operation also involves the appropriate low-pass anti-alias filtering. Denote the down-sampled $`g^{}`$ by $`z`$, which is complex-valued.

From $`z`$, it is convenient to construct the real vector sequence $`𝒛[k]`$:

$$𝒛[k]=\left(\begin{array}{c}\mathrm{}(z)[k]\\ \mathrm{}(z)[k]\end{array}\right).$$
(48)
The vector-valued time series $`𝒛[k]`$ plays the role of the observation in our Kalman filter implementation.

### C Dynamical Model

The gravitational wave channel $`g`$ is the sum of the violin modes (whose behavior we wish to determine) and other noise sources. In constructing $`𝒛`$ we focused attention on a narrow band about the mode resonant frequencies. Here we provide a model for the mode behavior as reflected in $`𝒛`$.

We first note that, near resonance, the details of the damping are unimportant; consequently, we can model each violin mode as an independent, viscously damped harmonic oscillator whose coordinate $`u`$ is the amplitude of the mode present in the detector output. The coordinate $`u`$ satisfies the differential equation

$$\ddot{u}+\frac{\omega _0}{Q}\dot{u}+\omega _0^2u=F(t),$$
(49)
where $`F`$ is the (stochastic) driving force. The “gravity-wave” channel is the sum of $`u`$ and other noises (and signals).

We can relate $`u`$ to the contribution of the violin modes to $`𝒛`$ by mixing $`u`$ with the phase factor $`\mathrm{exp}(i\omega _ct)`$ (where $`\omega _c`$ is equal to $`2\pi f_c`$):

$$\psi =u\mathrm{exp}(i\omega _ct).$$
(50)
The complex $`\psi `$ satisfies the differential equation

$$\ddot{\psi }+\left(\frac{\omega _0}{Q}+2i\omega _c\right)\dot{\psi }+\left(\omega _0^2\omega _c^2+i\frac{\omega _c\omega _0}{Q}\right)\psi =F\mathrm{exp}(i\omega _ct),$$
(52)
or, in the Laplace domain,

$$\left(sp_+\right)\left(sp_{}\right)\overline{\psi }(s)=\overline{F}(s+i\omega _c)$$
(53)
where

$`\overline{g}(s)`$ $`=`$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑t\mathrm{exp}(st)g(t)`$ (54)
$`p_\pm `$ $`=`$ $`{\displaystyle \frac{\omega _0}{2Q}}i\omega _c\pm i\omega _0\sqrt{1{\displaystyle \frac{1}{4Q^2}}}.`$ (55)
(Note that in equations 50 we have assumed that $`\psi `$ is causal — i.e., it vanishes at sufficiently earlier times — and used a two-sided Laplace transform. Written in this way the Laplace transform $`\overline{\psi }`$ of $`\psi `$ is related to its Fourier transform $`\stackrel{~}{\psi }`$ through the substitution of $`2\pi if`$ for $`s`$:

$$\stackrel{~}{g}(f):=_{\mathrm{}}^{\mathrm{}}𝑑t\mathrm{exp}(2\pi ift)g(t).$$
(56)
(Note that we have adopted the engineering convention for the Fourier transform.) After down-sampling, the real and imaginary parts of $`\psi `$ are the contributions of the violin modes to the real and imaginary parts of $`𝒛`$.

The two poles $`p_\pm `$ in the Laplace transform solution $`\overline{\psi }(s)`$ correspond to the positive and negative frequency resonances of the oscillator. When we choose $`\omega _c`$ close to $`\omega _0`$ we placed one of the resonances in $`\psi `$ near zero angular frequency and the other far away: correspondingly, one of the poles $`p_\pm `$ is near the origin and one is distant. Letting $`p_+`$ be the “near” pole and $`p_{}`$ be the “far” pole we write

$$\left(sp_+\right)\overline{\psi }(s)=\frac{\overline{F}(s+i\omega _c)}{sp_{}},$$
(57)
or, for $`|s||p_{}|`$,

$$\left(sp_+\right)\overline{\psi }(s)=p_{}^1\overline{F}(s+i\omega _c).$$
(58)
The spectrum of $`\psi `$ arising from this second equation, which differs from the first only by the source term, is identical to the spectrum of $`\psi `$ arising from the first as long as $`|\omega |`$ is much less than $`|p_{}|`$. Since down-sampling restricts the bandwidth of $`\psi `$ in just this way we adopt equation 58 as the evolution equation for $`\psi `$. Correspondingly, near zero angular frequency $`\psi `$ satisfies

$$\dot{\psi }p_+\psi =p_{}^1F(t)e^{i\omega _ct}.$$
(59)
Now write, instead of the complex $`\psi `$, the real vector $`𝝍`$ whose components are the real and imaginary parts of $`\psi `$:

$$𝝍=\left(\begin{array}{c}\mathrm{}(\psi )\\ \mathrm{}(\psi )\end{array}\right).$$
(60)
The first order differential equation 59 describing complex $`\psi `$ near zero frequency becomes a coupled set of first order differential equations

$$\dot{𝝍}=𝑨^{\mathbf{}}𝝍+𝑭$$
(62)
where

$`𝑨^{\mathbf{}}`$ $`=`$ $`\left(\begin{array}{cc}\mathrm{}(p_+)& \mathrm{}(p_+)\\ \mathrm{}(p_+)& \mathrm{}(p_+)\end{array}\right)`$ (65)
$`𝑭`$ $`=`$ $`F\left(\begin{array}{c}\mathrm{}(p_{}^1e^{i\omega _ct})\\ \mathrm{}(p_{}^1e^{i\omega _ct})\end{array}\right).`$ (68)
Our observation is not of continuous $`z(t)`$ but of sampled $`𝒛[k]`$. Suppose that the final sampling rate is $`\mathrm{\Delta }f`$ and approximate the evolution of $`𝝍`$ as free between samples, with the accumulated effect of the stochastic force $`𝑭`$ acting at the sample times. Then the sampled $`𝝍[k]`$ satisfies

$$𝝍[k]=𝑨𝝍[k1]+𝒘[k]$$
(70)
where

$$𝑨=\mathrm{exp}\left(\frac{𝑨^{}}{\mathrm{\Delta }f}\right)$$
(71)
and $`𝒘[k]`$ is a random vector related to $`𝑭`$. The quantity $`𝝍[k]`$ in this case is equal to the contribution of the violin mode to the observable $`𝒛[k]`$, or

$$𝒛[k]=𝑪𝝍[k]+𝒗[k],$$
(72)
where $`𝒗`$ represents the (approximately) white background in which the violin modes are immersed. Equations III C are the state-space form of the Kalman equations. To complete our description of the suspension modes in this form it remains to determine the properties of the process and measurement noises $`𝒘`$ and $`𝒗`$.

The statistical properties of $`𝒘`$ derive ultimately from the statistical properties of the force $`F`$ (cf. 49). In applying the Kalman filter we assume that $`𝒘`$, and consequently $`F`$, is white. Assuming that modes are driven principally by thermal forces we can use the Fluctuation-dissipation theorem  to find the power spectrum of the thermal force driving the mode:

$$S_F(f)=4k_BT\mathrm{}(Z)$$
(73)
where $`Z`$ is the system’s impedance. In our application,

$$Z=\frac{\omega ^2+i\gamma \omega \omega _0^2}{i\gamma \omega };$$
(74)
so, $`\mathrm{}(Z)`$ is independent of frequency, the thermal force acting on the modes is white, and the assumption made in applying the Kalman filter to this problem is satisfied exactly.

The quantity $`𝒗`$ is the measurement noise: i.e., the noise within which the violin mode signal is embedded. In our application the measurement noise is associated with the laser shot noise. In the current generation of large detectors the shot noise has the spectrum

$$S_h(f)=S_0\left[1+\left(\frac{f}{f_k}\right)^2\right],$$
(75)
with knee frequency $`f_k`$ approximately 100 Hz. the LIGO 40 M prototype the mode fundamental frequencies are confined to an approximately 40 Hz bandwidth about 590 Hz; in the LIGO I detector they will be in a similarly narrow band about 350 Hz. As long as we are able to confine attention to bandwidths small compared to the mode frequency the shot noise in that band will be very nearly white; correspondingly, $`𝒗`$ will be white and have a covariance $`𝑽`$ proportional to the unity matrix.

To summarize, we construct the observable $`𝒛`$ from the gravitational wave channel $`g`$ by mixing it with a local oscillator and down-sampling to the bandwidth $`\mathrm{\Delta }f`$ about zero frequency. We model the sampled observable as arising from the equations

$`𝝍[k]`$ $`=`$ $`𝑨𝝍[k1]+𝒗[k1]`$ (77)
$`𝒛[k]`$ $`=`$ $`𝑪𝝍[k]+𝒘[k]`$ (78)
where

$`𝑨`$ $`=`$ $`\mathrm{exp}\left({\displaystyle \frac{𝑨^{}}{\mathrm{\Delta }f}}\right)`$ (79)
$`𝑨^{}`$ $`=`$ $`\left(\begin{array}{cc}\omega _0/2Q& \omega _c\omega _0\sqrt{11/4Q^2}\\ \omega _c+\omega _0\sqrt{11/4Q^2}& \omega _0/2Q\end{array}\right)`$ (82)
### D Reconstruction

In applying the Kalman filter we determine, at each sample, the state estimate $`\widehat{𝝍}[k]`$ from which we form an estimate $`\widehat{𝒛}_v`$ of the contribution of the violin mode to the detector output:

$$\widehat{𝒛}_v=𝑪\widehat{𝝍}.$$
(83)
The violin mode contribution to the observation is an instrumental artifact that is unrelated to the incidence of gravitational waves on the detector; correspondingly, we may wish to use the prediction to subtractively remove the artifact from the data stream.In section IV we show that the Kalman filter estimate of the violin mode contribution to the observation is insensitive to the presence of a gravitational wave signal.

The prediction $`\widehat{𝒛}_v`$ for the violin mode contribution is formed at reduced bandwidth and offset frequency. From $`\widehat{𝒛}_v`$ we can construct the predicted contribution violin mode contribution to the full bandwidth observation ($`g[k]`$, cf. 47) by reversing the steps taken to form $`𝒛`$ from $`g`$ (cf. sec. III B). In particular, we first form the complex scalar sequence

$$\widehat{z}_v=\left(\begin{array}{cc}1& i\end{array}\right)\widehat{𝒛}_v$$
(84)
and up-sample it to the full sample rate, forming $`\widehat{G}_v^{}`$. The mixing operation undertaken in equation 47 can then be reversed,

$$\widehat{G}_v[k]=\widehat{G}^{}[k]\mathrm{exp}\left(2\pi ik\frac{f_c}{f_s}\right).$$
(85)
The result is still complex: it is missing the contributions from negative frequencies that were discarded in the low-pass filtering that accompanied the original down-sampling operation. Since the original observation was real, however, this contribution is just the complex conjugate of the positive frequency contribution, which we have; consequently, the estimated contribution of the violin mode to the data stream $`g[k]`$ is

$$\widehat{g}_v[k]=2\mathrm{}(\widehat{G}[k]).$$
(86)
The residual difference

$$g_r[k]=g[k]\widehat{g}_v[k]$$
(87)
is then free of violin mode artifact, to the accuracy of our estimation procedure.

### E Estimating the measurement and process noise

To use the Kalman filter we must have a model for the process dynamics ($`𝑨`$) and how the process contributes to the observable ($`𝑪`$). These we described for our particular application — the violin modes of an interferometer mirror suspensions (each mode an independent process) — in section III B. Additionally, we must have estimates for the process and measurement noise co-variance matrices ($`𝑾`$ and $`𝑽`$). Here we describe how we estimate these from a short segment of filter input, taken at the beginning of the input.

Consider the typical bandwidth of the observable $`𝒛`$, shown in figure 1. The process contributes significantly over only a small fraction of the bandwidth, where the resonance rises above the measurement noise floor. To estimate the measurement noise, we form the PSD of $`z`$ and take the average power far from the resonant peaks. This average power $`\overline{P}_m`$ is the total measurement noise (tr$`𝑽`$). We assume that there is no cross-covariance between the noise in $`\mathrm{}(z)`$ and $`\mathrm{}(z)`$, in which case

$$𝑽=\frac{1}{2}\overline{P}_m\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right).$$
(88)
To estimate the process noise we focus on that part of the spectrum near the resonant peak where the process dominates. We assume, in our model, that the process is driven by white noise. From the Kalman filter model, the power spectral density of the contribution to $`z`$ in this bandwidth is given by (cf. equation 59)

$$P_z(\omega )=\frac{4Q^2\text{tr}𝑾}{\omega _0^2+\left(2Q\omega _c\omega _0\left(4Q^21\right)^{1/2}\right)^2}.$$
(89)
With the exception of the total process noise (tr$`𝑾`$) the right hand side of equation 89 is known. The left hand side we measure. To estimate tr$`𝑾`$, we measure $`P_z`$ and sum it over the limited bandwidth dominated by the process. Setting that equal to the same sum executed on the right hand side of equation 89 we determine the proportionality constant tr$`𝑾`$. Assuming that the process noise is equally distributed between the real and imaginary parts of $`𝒛`$, with no cross-variance, we thus determine the estimated process noise covariance:

$$𝑾=\frac{1}{2}\left(\text{tr}𝑾\right)\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right).$$
(90)
## IV Characterization

Before providing results from the application of our filter on LIGO 40 M prototype data we characterize the filter’s behavior on “synthetic” or “mock” time series. These time series are constructed so that the violin mode state, process noise, measurement noise, and signal content are known exactly. By applying the filter to these data sets and investigating the difference between its estimation of these quantities and their actual, known values, we characterize the filter’s performance.

Four applications to synthetic data sets are described here:

1. A single isolated violin mode, driven by Gaussian white process noise and observed in the presence of Gaussian white measurement noise (cf.s̃ec. IV B);
2. Several closely spaced violin modes, driven by non-Gaussian white process noise and observed in the presence of non-Gaussian white measurement noise (cf. sec. IV C);
3. A single isolated violin mode, driven by Gaussian white process noise with a superposed impulsive excitation and observed in the presence of non-Gaussian white measurement noise (cf. sec. IV D); and
4. A single isolated violin mode, driven by non-Gaussian white process noise and observed in the presence of non-Gaussian white measurement noise with a superposed impulsive excitation (cf. sec. IV E).

The first case demonstrates the filter operating in the regime where the measurement noise and process noise satisfy the assumptions made in the derivation of the Kalman filter. The second case examines the filter’s ability to distinguish a violin mode from measurement noise, and to distinguish between two different modes, when the noise is not as simple as Gaussian and modes are closely spaced in frequency. In the third case we test the filter’s ability to respond to impulsive excitations of the suspension modes. In the final case we introduce an impulse into the measurement noise and observe the filter’s ability to discriminate between it and the contribution from the violin mode.

As they appear in the detector data stream, violin modes are characterized entirely by their resonant frequency $`f_0`$ and damping constant $`Q`$. Since we are dealing always with discrete-time signal processing, the scale of the frequency does not matter: the only relevant quantity (aside from $`Q`$) is the ratio of the resonant frequency to the sampling frequency $`f_s`$. Nevertheless, to foreshadow our application of the filter to data taken at the LIGO 40 M prototype, we assume that the sampling frequency is 9868.421 Hz, the sampling rate of the LIGO 40 M prototype, and we always choose the simulated violin modes to have resonant frequencies coincident with the frequencies of modes identified in that instrument. With this convention, the frequency and damping constants of the simulated modes used in the examples to follow are given in table I.

### A Generating synthetic data sets

Before discussing the results of our investigations into the Kalman filter’s behavior we describe how we generated the random time series used in those investigations.

In all cases the measurement noise is white: i.e., its power spectrum was constant over its bandwidth. For Gaussian noise measurement noise we used Matlab’s  randn random number generator. Non-Gaussian noise was modeled by a two-component mixture-Gaussian model:

$`P(x|\alpha ,\mu _1,\mu _2,\sigma _1,\sigma _2)`$ $`=`$ $`\alpha N(\mu _1,\sigma _1)+\left(1\alpha \right)N(\mu _2,\sigma _2)`$ (92)
$`\begin{array}{ccccccc}\hfill \alpha & =& 0.7,\hfill & & & & \\ \hfill \mu _1& =& 0,\hfill & & \hfill \mu _2& =& 0,\hfill \\ \hfill \sigma _1& =& 1,\hfill & & \hfill \sigma _2& =& 1.5.\hfill \end{array}`$ (96)
Here $`N(\mu ,\sigma )`$ represents a Normal process with mean $`\mu `$ and variance $`\sigma ^2`$, the parameters $`\alpha `$, $`\mu _k`$ and $`\sigma _k`$ are fixed with $`0\alpha 1`$. Figure 2 shows the distribution of this mixture Gaussian process together with the distribution of a Gaussian process with the same mean and variance. The mixture Gaussian distribution described by the parameters in equation IV A is strongly leptokurtic (i.e., it has a significant excess of high amplitude events).

To generate the contribution to the observation from the process we begin with the differential equation modeling the process (cf. 49):

$$\ddot{r}+\frac{\omega _0}{Q}\dot{r}+\omega _0^2r=F,$$
(97)
where $`r`$ is the position of the suspended mass. We make no assumptions here about driving force $`F`$. The detector gravitational wave channel consists of a superposition of measurement noise $`v`$ and the oscillator coordinate $`r`$. The solution to this equation can be expressed as an all-pole linear filter $`H`$ acting on the driving force $`F`$. In the Laplace domain the filter is given by

$`\overline{r}(s)`$ $`=`$ $`\overline{H}(s)\overline{F}(s)`$ (99)
$`\overline{H}(s)`$ $`=`$ $`{\displaystyle \frac{1}{s^2+s\omega _0/Q+\omega _0^2}}`$ (100)
$`=`$ $`{\displaystyle \frac{1}{\left(sp_+\right)\left(sp_{}\right)}}`$ (101)
where

$$p_\pm =\left(\omega _0/2Q\pm i\omega _0\sqrt{11/4Q^2}\right)$$
(102)
The frequency response of the filter at angular frequency $`\omega `$ is equal to $`\overline{H}(s)`$ evaluated on the imaginary axis:

$$\stackrel{~}{H}(\omega )=\overline{H}(i\omega ).$$
(103)
The two-pole filter $`H`$ operates in continuous time while the detector output we wish to simulate is discretely sampled at a rate $`f_s`$, which we assume is at least twice as great as the resonant frequency $`f_0`$. Using the analog filter $`\overline{H}(s)`$ as a prototype, we design a digital filter with a similar response at frequencies much less than the Nyquist frequency.

Suppose that the discretely sampled sequence $`r[k]`$, driven by the input $`F[k]`$, satisfies

$$u[k]=\underset{j=0}{\overset{N_b}{}}b[j]F[kj]\underset{j=1}{\overset{N_a}{}}a[j]u[kj],$$
(104)
where $`a[j]`$ and $`b[j]`$ are the constant filter coefficients. Define the $`z`$-transform of the sequence $`y[k]`$ by

$$\stackrel{ˇ}{y}\underset{k=\mathrm{}}{\overset{\mathrm{}}{}}y[k]z^k.$$
(105)
(The $`z`$-transform of a digital filter is analogous to the Laplace transform of a continuous filter.) In $`z`$-space, the linear filter described in equation 104 can be written

$`\stackrel{ˇ}{r}(z)`$ $`=`$ $`{\displaystyle \frac{\underset{j=0}{\overset{N_b}{}}b[j]z^j}{1+_{j=1}^{N_a}a[j]z^j}}\stackrel{ˇ}{F}(z)`$ (107)
$`=`$ $`\stackrel{ˇ}{G}(z)\stackrel{ˇ}{F}(z)`$ (108)
The frequency response of the digital filter $`G`$ is equal to $`\stackrel{ˇ}{G}(z)`$ evaluated on the unit circle in the complex $`z`$-plane:

$$\stackrel{ˇ}{G}\left(e^{2\pi if/f_s}\right)$$
(109)
where $`f_s`$ is the sample rate of the filtered sequence. If, by a conformal transformation, we map the right half $`s`$-plane on to the interior of the unit circle in the $`z`$-plane, then we will have mapped the analog prototype $`\overline{H}(s)`$ on to the digital filter $`\stackrel{ˇ}{G}(z)`$:

$$\stackrel{ˇ}{G}(z)=\overline{H}(s)|_{s=2\mathrm{\Omega }(z1)/(z+1)},$$
(110)
where $`\mathrm{\Omega }`$ is the so-called warping constant. This mapping will not be without distortion; nevertheless, at frequencies much less than the Nyquist frequency the distortion will be minimal. We can choose the warping constant $`\mathrm{\Omega }`$ to guarantee that the digital filter’s response matches exactly the analog filter’s response at exactly one frequency:

$$\stackrel{ˇ}{G}(e^{2\pi if^{}/f_s})=\overline{H}(2\pi if^{})$$
(111)
when

$$\mathrm{\Theta }=\frac{\pi f^{}}{\mathrm{tan}(\pi f^{}/f_s)}.$$
(112)
We choose the warping constant so that the response of the digital filter matches the response of the analog filter at the resonant frequency $`f_0`$. This completes the design of the digital filter $`\stackrel{ˇ}{G}(z)`$. Given a pseudo-random process $`F[k]`$ and the recursion relationship of equation 104 we evaluate $`r[k]`$, the simulated contribution of the violin mode to the detector output. Finally, to obtain the simulated gravitational-wave channel output, we add to $`r[k]`$ another pseudo-random process $`n[k]`$, which represents the measurement noise.

In our simulations we used the $`\mathrm{csc}Matlab`$  Normal random number generatorrandn wherever Normal deviates were called for: i.e., in forming either the measurement noise, the process noise or the stochastic driving force $`F`$.

### B Single mode: Normal, white process and measurement noise

Using the techniques described above we generated several hundred seconds of simulated data, consisting of a single mode with the frequency and damping constants given in the first column of table I. The process noise was white with variance of 2.38; the measurement noise was also white with variance 2.8.

After generating the data set we applied the Kalman filter, focusing on a 6 Hz bandwidth $`\mathrm{\Delta }f`$ centered about the mode’s resonant frequency. Following the discussion in section III E we estimated $`𝑾`$ and $`𝑽`$ using approximately 100 s of the simulated data.Since this data set is simulated the values of $`𝑾`$ and $`𝑽`$ are known a priori; however, in a real application they will not be known. We used the same procedure to estimate the filter as would be used in the field. We applied the Kalman filter to the subsequent simulated data, using these values of $`𝑾`$, $`𝑽`$, setting the initial state estimate $`\widehat{𝝍}`$ to zero and the initial estimated error $`𝑷`$ equal to $`𝑾`$$`𝑰`$ (where $`𝑰`$ is the identity matrix). Figure 3 summarizes, in three panels, the Kalman filter’s effectiveness in estimating the mode state. The upper panel shows the power spectrum of the simulated observations in a wide band about the mode frequency. The middle panel shows the power spectrum of the Kalman estimate of the mode contribution to the observable. Finally, the lower panel shows the power spectral density of the difference between the actual data and the estimated contribution from the mode.\**\**\**While the power spectral density in a relatively wide band about the mode resonant frequency is shown, subtraction of the estimate affects only a 6 Hz bandwidth about the mixing frequency $`f_c`$ (cf. sec. III B and eqns. 50, 85).. Comparing these three panels it is clear that the Kalman filter has correctly identified the modes contribution to the observable.

Since the data is simulated we know exactly the measurement noise and mode contributions to the observable. Correspondingly, we can evaluate the actual error in the filter’s estimate of the mode state. The mode (estimated) state is described by the complex ($`\widehat{\psi }[k]`$) $`\psi [k]`$, which we describe by its amplitude and phase

$`A`$ $`:=`$ $`|\psi |`$ (114)
$`\mathrm{tan}\varphi `$ $`:=`$ $`\mathrm{}(\psi )/\mathrm{}(\psi ).`$ (115)
The difference between the actual and estimated mode states we describe by the quantity $`\mathrm{\Delta }`$:

$$\mathrm{\Delta }:=\frac{2|\psi [k]\widehat{\psi }[k]|}{|\psi [k]+\widehat{\psi }[k]|}.$$
(116)
Figure 4 shows the instantaneous estimated mode amplitude, phase and error $`\mathrm{\Delta }`$. The top panel shows the estimated mode amplitude. The middle panel shows the estimate mode phase, less $`2\pi f_0t`$, where $`f_0`$ is the mode frequency. The bottom panel shows the error in the estimated mode state. Following an initial epoch, during which the filter is “locking-on” to the mode state, the relative error in the state estimate falls quickly to zero.

The middle and bottom panel of figure 3 show a slight depression in the power spectral density of the residual in the band that the Kalman filter operates. This slight depression is associated with the error in the state estimate, which is shown in the bottom panel in figure 4 to be on order several percent. The depression is uniform throughout the band, showing that the errors introduced by the subtraction of the estimated mode contribution from the data are white. For this reason it is unlikely that these errors can be reduced by, for example, a “better” estimate of $`𝑾`$, $`𝑽`$, $`f_0`$ or $`Q`$: changing any of these would lead to an error that was not white, but reflected either an over- or under-estimate of the power in the mode relative to the measurement error.

We expect that the continuously updated mode amplitude and phase estimates as a function of time, either analyzed directly or displayed in a strip chart of which the top two panels of figure 4 may be considered a snapshot, will be particularly useful for monitoring the state of the detector.

### C Multiple modes; non-Gaussian white measurement and process noise

In deriving the Kalman filter we assumed that the process and measurement noise were white; however, we never required that they be Gaussian. Here we investigate the effectiveness of the filter when these noises are non-Gaussian and when there are several modes present, at least two of which are very close together in frequency. This is a critical application, since in actual data there are multiple closely spaced violin modes resulting from several suspension wires of the mirrors, which cannot be estimated separately (cf. 1). The mode frequencies and damping constants for this example are given in the second column of table I.

For isolated modes (i.e., those whose resonant peaks do not overlap with each other above the level of the measurement noise), one can design and apply separately a Kalman filter for each mode. Closely spaced modes cannot be treated in isolation; however, since the mode dynamics are independent, the development of a Kalman filter that treats $`N`$ modes jointly is straightforward:

* Referring to section III B, fix a carrier frequency $`f_c`$ and a bandwidth $`\mathrm{\Delta }f`$ such that the modes of interest all have resonances in the band $`[f_c\mathrm{\Delta }f/2,f_c+\mathrm{\Delta }f/2]`$.
* Considered separately, the dynamics of the mode $`j`$ ($`1jN`$) and the contribution $`𝒛_j`$ of the mode to the measurement are described by the equations

$`𝝍_j[k]`$ $`=`$ $`𝑨_j𝝍_j[k1]+𝒘_j[k1]`$ (118)
$`𝒛_j[k]`$ $`=`$ $`𝑪_j𝝍_j[k].`$ (119)
* The Kalman equations describing the joint evolution of the $`N`$ modes, and their contribution to the observation $`𝒛`$ are then

$`𝚿[k]`$ $`=`$ $`𝑨𝚿[k1]+𝒘[k1]`$ (121)
$`𝒛[k]`$ $`=`$ $`𝑪𝚿[k]+𝒗[k],`$ (122)
where

$`\mathrm{\Psi }`$ $`:=`$ $`\left(\begin{array}{c}𝝍_1\\ \mathrm{}\\ 𝝍_N\end{array}\right),`$ (126)
$`=`$ $`\left(\begin{array}{c}\mathrm{}(\psi _1)\\ \mathrm{}(\psi _1)\\ \mathrm{}\\ \mathrm{}(\psi _N)\\ \mathrm{}(\psi _N)\end{array}\right),`$ (132)
$`𝑨`$ $`:=`$ $`\left(\begin{array}{ccc}𝑨_1& & \\ & \mathrm{}& \\ & & 𝑨_N\end{array}\right),`$ (136)
$`𝑾`$ $`:=`$ $`\left(\begin{array}{ccc}𝑾_1& & \\ & \mathrm{}& \\ & & 𝑾_N\end{array}\right),`$ (140)
$`𝑽`$ $`:=`$ $`\left(\begin{array}{ccc}𝑽_1& & \\ & \mathrm{}& \\ & & 𝑽_N\end{array}\right),`$ (144)
$`𝑪`$ $`:=`$ $`\left(\begin{array}{ccc}𝑪_1& \mathrm{}& 𝑪_N\end{array}\right)`$ (146)
$`=`$ $`\left(\begin{array}{ccccc}1& 0& \mathrm{}& 1& 0\\ 0& 1& \mathrm{}& 0& 1\end{array}\right),`$ (149)
and $`𝑾_j`$ and $`𝑽_j`$ are, as before, the process and measurement noise.

In constructing $`𝑾`$ we have approximated the individual modes as being driven by independent process noise sources. Correlated process noise leads to correlations in the states of the individual modes. Even in the absence of an explicit correlation in the assumed process noise covariance (i.e., off-diagonal terms in $`𝑾`$), the filter will, if warranted by the observations, correlate the state estimates. Lacking knowledge of the driving noise (beyond contact with independent thermal baths), assuming that the process noise is uncorrelated guarantees that there is no prejudice that the mode states are correlated without constraining the state estimates to be uncorrelated.

To evaluate the effectiveness of the Kalman filter in treating multiple modes driven by non-Normal noise, we applied a single Kalman filter to this three mode system, choosing $`f_c`$ equal to 578.6 Hz and $`\mathrm{\Delta }f`$ equal to 18 Hz. The top panel of figure 5 shows the power spectral density of the simulated gravity wave channel in a 60 Hz bandwidth about the carrier frequency $`f_c`$, while the bottom panel shows the power spectral density of the residual difference between the simulated observation and the filter’s prediction of the contribution from the violin modes. The residual broad feature that remains following the subtraction of the mode estimate is due to the non-Gaussian nature of the process noise that drives the mode. The filter has made a best-estimate of the mode state, given only knowledge of the driving force variance. In this example, however, the variance provides a very incomplete picture of the process noise. Still, the filter has clearly identified approximately 45 dB of the approximately 55 dB of the mode contribution to the measurement.

The state $`𝚿`$ is the direct sum of the state vector estimates $`𝝍_j`$ of the individual modes; consequently, the filter identifies not just the total contribution of each mode to the observation, but the state of each mode separately. The four panels of figure 6 show the estimated amplitude and phase (less $`2\pi f_jt`$ for each mode) of each of the two overlapping modes at 584.6 Hz ($`f_2`$) and 585.4 Hz ($`f_3`$). In this example the two modes are known to be uncorrelated; consequently, we expect that the state estimates should also be uncorrelated, which we have found to be the case.

### D Effect of a transient excitation of a mode

To evaluate the effectiveness of the Kalman filter in estimating the filter state in the presence of a transient excitation, we have simulated a single mode observed in the presence of the same mixture-Gaussian measurement noise described above (cf. equation IV A). The process noise is taken to be Gaussian with variance 2.38. At 20 seconds into the simulation we add to the process noise a single pulse of magnitude $`35\sigma `$, where $`\sigma `$ is the rms amplitude of the process noise. The top panel of figure 7 shows the estimated amplitude of the mode while the bottom panel shows error $`\mathrm{\Delta }`$ in the mode state estimate:

$$\mathrm{\Delta }:=\frac{2|\psi [k]\widehat{\psi }[k]|}{|\psi [k]+\widehat{\psi }[k]|}.$$
(150)
After the filter’s start-up transient the error in the estimated state is small until the transient excitation. There is, not surprisingly, a sudden increase in the error at the moment of the excitation; however, the filter quickly adjusts to the new state.

### E Effect of transients in the measurement noise

Sudden excitations in the observed time series may arise from sources other than excitations of the violin modes. To evaluate the effectiveness of the Kalman filter in discriminating between excitations of the mode and measurement noise transients we simulate a single line, driven by Gaussian noise. This mode is embedded in the same mixture-Gaussian measurement noise as described in section IV C, but with a pulse in the measurement noise 14 seconds into the simulation. The pulse lasts for a single sample and has an amplitude of $`35\sigma `$, where $`\sigma `$ rms amplitude of the measurement noise.

The lower panel of figure 8 shows $`\mathrm{\Delta }`$ (cf. 150), measuring the magnitude of the difference in the estimated and actual mode state; the upper panel shows the estimated state itself. The error in the mode state estimate always remains quite small, showing that the filter discriminates between excitations of the modes and excitations in the measurement noise.

### F Summary

The Kalman filter described here is able to correctly determine instantaneous state of the violin mode suspensions. This is true whether the mode is driven by, or observed in the presence of, Gaussian or non-Gaussian noise. The filter does not require that the measurement or process noise be stationary: it correctly tracks the mode after it has been impulsively excited and can distinguish between excitations of the mode and transients in the measurement noise.

## V Violin modes in LIGO 40 M prototype data

In this section we discuss the application of the Kalman filter to data taken at the LIGO 40 M prototype detector during November 1994. The general parameters of this run are described elsewhere . The particular results presented here all refer to the particular locked segment beginning at GPS time 469357019 s on 19 November 1994 and lasting for 2665.54496 s.

### A General character of the resonant features

At the time of the November 1994 data run each of the four cavity mirrors in the 40 M prototype was suspended by four cylindrical wires. Setting aside the coupling of the wires through their attachment to the suspended mirror there are two fundamental modes associated with each wire: one with motion nominally along the optical axis and one with motion perpendicular to that axis. Owing to imperfections in the wires, the attachments of the wires to the test masses, and the coupling of the modes through the motion of the mass all the modes will have some component along the optical axis. We thus expect 32 fundamental violin modes.

A high resolution power spectrum focused on the bandwidth (570 Hz, 606 Hz) and taken of the data from a single locked-segment (i.e., a single epoch during which the interferometer operated without interruption as a gravitational wave detector) of 44 m duration shows 28 distinct resonant features. One of these is clearly identified as the ninth harmonic of the 60 Hz power-main, leaving 27 different modes whose frequencies and quality factors are listed in table III.

Referring to table III note that, apart from modes 7, 16, 19, 24, and 25, the lines appear in pairs, with a high amplitude line closely associated with a lower amplitude line at a slightly higher frequency (less than or of order 0.2 Hz). This association of lines is clearly shown in figure 9. Of the five apparently unpaired lines, three show companions in a high resolution power spectrum taken in the band of the first harmonic, suggesting that all five of the “missing modes” may be so closely spaced in frequency with other, identified modes, that we are unable to resolve them.

### B Application

Using the techniques described in section III we have applied a Kalman filter designed to track the 27 violin modes identified in table III. The top panel of figure 10 shows the power spectrum of the IFO\_DMRO (interferometer differential mode read-out, or “gravity wave”) channel in the 45 Hz band beginning at 565.0 Hz, with the violin modes clearly present. The bottom panel shows the power spectrum of the residual after the modes identified by the filter have been subtractively removed. Very little power associated with these modes remains.††††††The 600 Hz line, which is the ninth harmonic of the power-main, was not modeled in our analysis.

Figure 10 also shows an overall approximately 2 dB reduction in the background noise level within the band, which is due to the “wings” of the predicted contribution of the violin modes to the total noise. Our experience with simulated data (cf. sec. IV BIV E) shows that the filter generally makes fractional prediction errors whose magnitude is on order 1% and that, as we have implemented the filter, this error tends to systematically over-estimate the contribution of the mode wings to the total noise. This is consistent with spectrum of the residual found here.

### C Preliminary characterization of mode motion and residual

The Kalman filter provides us with a decomposition of the observation into the contributions from the individual violin modes and a residual, which is, in a least-squares sense, as free from the influence of the modes as we can make it. The noise character of the violin modes, which are instrumental artifacts, can thus be studied separately from the noise character of the residual.

Figure 11 shows in histogram form the frequency with which samples of given amplitude (normalized by rms amplitude) appear in the $`(570,595)\text{Hz}`$ bandwidth of the original data (“raw”), in the contribution to this band from the modes as identified by the Kalman filter, and in the residual difference between the raw data stream and the estimated contribution from the modes. If the detector noise were strictly Gaussian the distribution of squared-amplitudes (normalized by the mean-squared amplitude) would follow an exponential distribution and the resulting curves would be, on the semi-log graph as shown, a straight lines with slope $`\mathrm{log}_{10}e`$. We have evaluated $`\chi ^2`$ for each histogram, using the left-most 47 bins shown in figure 11 and with the remainder of the data bins collapsed into a single additional bin, so that no bin has fewer than 40 events. The same binning was used for all three $`\chi ^2`$ statistic evaluations. Table IV gives the resulting $`\chi ^2`$, the number of degrees of freedom, and the corresponding $`p`$-value.

Referring to table IV we see that the Kalman filter has resolved the total noise in this band into two components that have very different noise character. While in all cases the statistical character of the noise is far from Gaussian, the contribution identified as arising from the violin modes is much better behaved (i.e., much closer to Gaussian) than the residual. The close correspondence between the $`\chi ^2`$ for the “raw” data and the identified violin mode contribution reflects the fact that the violin mode contribution, even while it is narrow-band, is of such great amplitude that it dominates the total detector noise in this band. Recalling our earlier results (cf. sec. IV D and IV E), which show how the Kalman filter is able to discriminate correctly between excitations of the modes and impulses in the measurement noise, we conclude that the difference in the noise character is real and that, in the LIGO 40 M, the violin mode contribution to the overall noise has many fewer outliers than the non-mode background. Since, in a real detector, it is the residual that will contain the gravitational wave signal, the importance of removing instrumental artifacts like the violin modes in order to discover the noise character of the residual cannot be overstated.

A still closer examination, focusing on the contributions from the individual modes, shows that certain modes have much greater excess noise than others. Figure 1214 shows, for twelve different violin modes, the frequency distribution of mode amplitudes. All of the modes shown here have a strong excess of high amplitude noise events: i.e., they are subject to strong, non-thermal forces.

## VI Conclusions

Violin modes of the test mass suspensions are a significant instrumental artifact in the signal band of modern interferometric gravitational wave detectors: at the mode resonant frequencies they stand between 40 and 50 dB above the ambient noise floor. These modes are not excited or otherwise disturbed by the passage of a gravitational wave; yet, they contribute almost all of the noise power in a wide bandwidth. Correspondingly, the ability to identify and remove their contribution to the detector noise and characterize the residual noise can significantly increase the detector’s sensitivity in the affected wave-band.

Treating the modes as a stochastic signal in a measurement noise background we have developed and demonstrated the use of a Kalman filter for identifying the state of the suspension wires and the contribution of the violin modes to the detector noise. A Kalman filter uses the known dynamics of the process that underlies a noise component (in this case the stochastic excitation of the suspension violin modes) to identify that component’s contribution to total system noise from measurements made on the total noise alone. When the stochastic driving force is Normal (i.e., Gaussian) the Kalman filter is optimal in the least squares sense. When the process noise is non-Gaussian, the Kalman filter estimate is optimal when only the second order statistics of the noise (i.e., its power spectrum) is known.

Using simulated data, whose mode and measurement noise contribution is known exactly, we have shown that the filter works well in the presence of both Gaussian and non-Gaussian measurement noise and excitations of the violin modes. We have also shown that the filter easily discriminates between transients in the measurement noise and transient excitations of the violin modes themselves, capturing the latter in its estimate of the mode state while rejecting the former as not arising from excitation of the mode. Using data taken at the LIGO 40 M we have demonstrated the use of the filter by isolating and subtractively removing the noise owing to 27 distinct violin modes, each standing between 40 and 50 dB above the ambient noise floor, over an approximately 40 Hz bandwidth.

The Kalman filter describes here operates in the time domain, providing instantaneous estimates of the suspension violin modes states. The computational cost is negligible, making it suitable for use as an on-line diagnostic as well as for on- or off-line data analysis.

Thermal and other excitations of the violin modes dominate the detector noise over a wide band. This additive contribution to the detector output will never include any gravitational wave signal, however. Separating the mode contribution from the total noise dramatically increases the detector’s sensitivity by significantly reducing the mean-square noise amplitude. It also, however, reveals a new “layer” of non-Gaussian and transient noise that must be characterized and, where possible, removed through its association with other instrumental and physical environment monitors. Using the Kalman filter to separate the total noise into a violin mode contribution and a residual we find that, in the LIGO 40 M prototype, the noise owing to the violin modes is much “cleaner” (i.e., much closer to Normal) than the residual noise, which is strongly leptokurtic. Since, in a real detector, it is the residual that will contain the gravitational wave signal, the importance of removing instrumental artifacts like the violin modes in order to discover the noise character of the residual cannot be overstated. As we dig deeper into the detector noise the challenges of identifying the origin of the newly uncovered noise components becomes greater; however, this is a challenge that must be faced and overcome if we are to make the most of the opportunities provided by the these new detectors.

## VII Acknowledgment

We are grateful to Albert Lazzarini for drawing our attention to Kalman Filtering and suggesting that it might be a useful tool for gravitational wave data analysis. We are also glad to thank the LIGO Laboratory at the California Institute of Technology for their hospitality during the 1997/8 academic year, when this work was begun. We also gratefully acknowledge the use of data taken at the LIGO 40 M interferometer during November 1994. SM wishes to thank Soumya D. Mohanty for many valuable insights and criticisms, and Constantin Brif for an extremely valuable set of notes on thermal noise in test mass suspensions of interferometric gravitational wave detectors. LSF is particularly glad to thank Peter Fritschel for valuable discussions. Finally, we are glad to thank Gabriela Gonzalez, Gary Sanders and David Tanner for their comments on this manuscript.

This work was supported by National Science Foundation award PHY 99-96213, PHY 98-00111 and their predecessors.


format: text

text:
# 1 Introduction

## 1 Introduction

In a discussion of generalized Born-Infeld theories  the construction of the Lagrangian was reduced to the solution of certain unilateral algebraic matrix equations and it was conjectured that the iterative solution of those equations is a sum of symmetric polynomials in the coefficients and of terms which are commutators. Equivalently, that the trace of the matrix solution is equal to the sum of traces of symmetric polynomials in the coefficients. The conjecture was later proven in  and by A. Schwarz in  using different methods and a slightly different form of the equation.

In the present note we combine Schwarz’ s idea of expressing the trace of the solution as a contour integral in the complex plane of the trace of the resolvent of the corresponding matrix with the idea used in  of using the basic property of the trace of the logarithm of matrices.

The two approaches can be easily combined by using the generalized Bezout theorem for matrix polynomials. The qualitative algebraic fact that the trace of the solution is given by a sum of traces of symmetrized polynomials in the coefficients emerges almost without any computation. The coefficients in the expansion are not hard to compute and we give an explicit expression for them. 111Our contour integral formulas are very similar to Schwarz’ s, but there seem to be some minor discrepancies between ours and his. The method used here applies equally to the two different equations considered in  and ; in the last section we clarify the relation between the two equations and their solutions.

Before closing this section we recall the statement of the generalized Bezout theorem (see e.g. ). Let

$$P(\lambda )=a_0+a_1\lambda +a_2\lambda ^2+\mathrm{}+a_n\lambda ^n,\lambda $$
(1.1)
be a polynomial with square matrix coefficients $`a_r`$ and $`x`$ a square matrix of the same order. Define

$$P(x)=a_0+a_1x+a_2x^2+\mathrm{}+a_nx^n$$
(1.2)
with the coefficients all on the left. It is easy to verify that

$$P(\lambda )P(x)=Q(\lambda ,x)(\lambda x),$$
(1.3)
where

$$Q(\lambda ,x)=\underset{l=0}{\overset{n1}{}}\lambda ^l\left(\underset{r=l+1}{\overset{n}{}}ϵa_rx^{rl1}\right).$$
(1.4)
In other words $`\lambda x`$ is a divisor of $`P(\lambda )P(x)`$ on the right (if we had taken all the coefficients on the right, then $`\lambda x`$ would have been a divisor on the left).

We shall study matrix equations of the type

$$P(x)=0,$$
(1.5)
which we shall call left-unilateral matrix equations (meaning that all the coefficients are on the left). If $`x`$ is a solution of (1.5), the characteristic polynomial $`P(\lambda )`$ of the equation can be factorized as

$$P(\lambda )=Q(\lambda ,x)(\lambda x).$$
(1.6)
It is clear that the Bezout theorem applies more abstractly if one considers $`a_r`$ and $`x`$ as elements of an associative, but possibly noncommutative, algebra. The same remark applies to the rest of this note, if one uses an appropriate algebraic definition of the trace as cyclic average (see ).

## 2 The Trace of the Solutions of Unilateral Matrix Equations

A. Schwarz  considers the unilateral matrix equation

$$x^n=1+ϵ\left(a_0+a_1x+\mathrm{}+a_{n1}x^{n1}\right).$$
(2.1)
where $`ϵ`$ is a small parameter. For $`ϵ=0`$ the equation has $`n`$ solutions in the complex plane, the $`n`$ roots of $`1`$. For $`ϵ0`$ each of these solutions admits a perturbative expansion as a formal power series in $`ϵ`$. A. Schwarz proves that $`\text{Tr}x^s`$ can be expressed as a power series whose coefficients are symmetrized products of the $`a_i`$ and gives an expression for the coefficients in terms of contour integrals in the complex plane.

We have found that a more explicit expression for $`\text{Tr}x^s`$ is provided by

$$\text{Tr}x^s=\text{Tr}\mathrm{\hspace{0.25em}1}+s\underset{k=1}{\overset{\mathrm{}}{}}\frac{ϵ^k}{n^k}\underset{n_0+\mathrm{}+n_{n1}=k}{}\frac{\text{Tr}𝒮(a_0^{n_0}\mathrm{}a_{n1}^{n_{n1}})}{n_0!\mathrm{}n_{n1}!}\underset{r=1}{\overset{k1}{}}\left(s+\underset{l=1}{\overset{n1}{}}ln_lrn\right).$$
(2.2)
This formula holds for positive as well as for negative values of $`s`$. Here the expansion is made around $`1`$. A similar expansion could be derived for all the $`n`$ iterative solutions of the equation. The normalization of the symmetrized product $`𝒮(a_0^{n_0}\mathrm{}a_n^{n_n})`$ is chosen in such a way as to give the ordinary product if the factors commute .

If we introduce the notation

$$\left(\begin{array}{c}\alpha \hfill \\ k\hfill \end{array}\right)=\{\begin{array}{cc}\underset{r=1}{\overset{k}{}}\frac{\alpha r+1}{r}\hfill & \text{ for }k=1,2,\mathrm{}\hfill \\ & \\ 1\hfill & \text{ for }k=0,\hfill \end{array}$$
(2.3)
which reduces to the usual definition of $`\left(\begin{array}{c}\alpha \hfill \\ k\hfill \end{array}\right)`$ for $`\alpha `$, then the result (2.2) can be rewritten

$`\text{Tr}x^s=`$ (2.8)
$`\text{Tr}\mathrm{\hspace{0.25em}1}+{\displaystyle \frac{s}{n}}{\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}ϵ^k(k1)!{\displaystyle \underset{n_0+\mathrm{}+n_{n1}=k}{}}{\displaystyle \frac{\text{Tr}𝒮(a_0^{n_0}\mathrm{}a_{n1}^{n_{n1}})}{n_0!\mathrm{}n_{n1}!}}\left(\begin{array}{c}\frac{1}{n}(sn+_{l=1}^{n1}ln_l)\hfill \\ \\ k1\hfill \end{array}\right).`$
As a first step to prove (2.2), we apply the generalized Bezout theorem. For (2.1) the characteristic polynomial is

$$P(\lambda )1\lambda ^n+ϵ\left(a_0+a_1\lambda +\mathrm{}+a_{n1}\lambda ^{n1}\right)$$
(2.9)
and, if $`P(x)=0`$,

$$P(\lambda )=Q(\lambda ,x)(\lambda x)$$
(2.10)
with

$$Q(\lambda ,x)=\underset{l=0}{\overset{n1}{}}\left(\underset{r=l+1}{\overset{n1}{}}ϵa_rx^{rl1}x^{nl1}\right)\lambda ^l.$$
(2.11)
Let $`x`$ be one of the solutions of (2.1), e.g. the one which reduces to $`1`$ for vanishing $`ϵ`$. (For the other solutions a similar procedure can be followed). We take the logarithm and then the trace of (2.10), and obtain

$$\text{Tr}\mathrm{log}P(\lambda )=\text{Tr}\mathrm{log}(Q(\lambda ,x))+\text{Tr}\mathrm{log}(\lambda x).$$
(2.12)
Differentiation with respect to $`\lambda `$ yields

$$\text{Tr}\frac{1}{P(\lambda )}P^{}(\lambda )=\text{Tr}\frac{1}{Q(\lambda ,x)}Q^{}(\lambda ,x)+\text{Tr}\frac{1}{\lambda x}.$$
(2.13)
More generally, we can multiply the above equation by any function $`f(\lambda )`$ of $`\lambda `$, which is regular in a neighbourhood of $`1`$, so as to get

$$\text{Tr}\frac{1}{P(\lambda )}P^{}(\lambda )f(\lambda )=\text{Tr}\frac{1}{Q(\lambda ,x)}Q^{}(\lambda ,x)f(\lambda )+\text{Tr}\frac{1}{\lambda x}f(\lambda ).$$
(2.14)
As a next step, we follow Schwarz’ s idea of using a contour integration to isolate the relevant part in the trace. We consider a small circle $`\mathrm{\Gamma }`$ around $`1`$ in the complex plane, or more generally a small closed curve winding once around $`1`$, and containing no other solution of (2.1) for $`ϵ=0`$. Then we compute the integral of equation (2.14) along it for small $`ϵ`$

$`(2\pi i)^1{\displaystyle _\mathrm{\Gamma }}𝑑\lambda \text{Tr}{\displaystyle \frac{1}{P(\lambda )}}P^{}(\lambda )f(\lambda )=`$
$`(2\pi i)^1{\displaystyle _\mathrm{\Gamma }}𝑑\lambda \text{Tr}{\displaystyle \frac{1}{Q(\lambda ,x)}}Q^{}(\lambda ,x)f(\lambda )+(2\pi i)^1{\displaystyle _\mathrm{\Gamma }}𝑑\lambda \text{Tr}{\displaystyle \frac{1}{\lambda x}}f(\lambda ).`$
The two integrals on the right hand-side can be evaluated through the Cauchy theorem in the following way. We are considering the case of small $`ϵ`$, and then $`Q`$ is a polynomial in $`x`$ and $`\lambda `$, which vanishes for $`\lambda `$ near $`e^{\frac{2\pi ik}{n}}`$, $`k=1,\mathrm{},n1`$, so that in this case its inverse $`Q^1`$ has no singularities near $`1`$. $`Q^{}`$ and $`f`$ are regular functions near $`1`$ and therefore the Cauchy theorem guarantees that the first term on the right hand-side vanishes:

$$(2\pi i)^1_\mathrm{\Gamma }𝑑\lambda \text{Tr}\frac{1}{Q(\lambda ,x)}Q^{}(\lambda ,x)f(\lambda )=0.$$
(2.16)
On the other hand $`\text{Tr}\frac{1}{\lambda x}f(\lambda )`$ has poles for $`\lambda `$ near $`1`$ with total residue $`\text{Tr}f(x)`$. As $`x`$ is close to $`1`$ for small $`ϵ`$, the Cauchy theorem yields

$$(2\pi i)^1_\mathrm{\Gamma }𝑑\lambda \text{Tr}\frac{1}{\lambda x}f(\lambda )=\text{Tr}f(x).$$
(2.17)
Finally, we obtain

$$\text{Tr}f(x)=(2\pi i)^1_\mathrm{\Gamma }𝑑\lambda \text{Tr}\frac{1}{P(\lambda )}P^{}(\lambda )f(\lambda ).$$
(2.18)
Therefore, the problem of computing $`\text{Tr}f(x)`$ amounts to evaluating the integral on the right hand-side of (2.18).

We can factorize

$$P(\lambda )=(1\lambda ^n)T(\lambda ),$$
(2.19)
where

$$T(\lambda )=1ϵ(\lambda ^n1)^1\underset{l=0}{\overset{n1}{}}a_l\lambda ^l.$$
(2.20)
Then

$$\frac{1}{P(\lambda )}P^{}(\lambda )=\frac{n\lambda ^{n1}}{1\lambda ^n}+\frac{1}{T(\lambda )}T^{}(\lambda )$$
(2.21)
and (2.18) becomes

$$\text{Tr}f(x)=\frac{1}{2\pi i}_\mathrm{\Gamma }𝑑\lambda \frac{n\lambda ^{n1}}{\lambda ^n1}\text{Tr}f(\lambda )+\frac{1}{2\pi i}_\mathrm{\Gamma }𝑑\lambda \text{Tr}\frac{1}{T(\lambda )}T^{}(\lambda )f(\lambda ).$$
(2.22)
The first integral on the right hand-side equals $`\text{Tr}f(1)`$. To evaluate the second term we integrate by parts

$$\frac{1}{2\pi i}_\mathrm{\Gamma }𝑑\lambda \text{Tr}\frac{1}{T(\lambda )}T^{}(\lambda )f(\lambda )=\frac{1}{2\pi i}_\mathrm{\Gamma }𝑑\lambda \text{Tr}\mathrm{log}(T(\lambda ))f^{}(\lambda ).$$
(2.23)
This is justified, because $`\text{Tr}\mathrm{log}T(\lambda )`$ and $`f(\lambda )`$ are regular on $`\mathrm{\Gamma }`$, so that

$$\frac{1}{2\pi i}_\mathrm{\Gamma }𝑑\lambda \text{Tr}\frac{d}{d\lambda }\left(\mathrm{log}(T(\lambda ))f(\lambda )\right)=0.$$
(2.24)
In this way we obtain

$$\text{Tr}f(x)=\text{Tr}f(1)\frac{1}{2\pi i}_\mathrm{\Gamma }𝑑\lambda \text{Tr}\mathrm{log}(T(\lambda ))f^{}(\lambda ).$$
(2.25)
We expand the logarithm

$$\mathrm{log}(1\alpha )=\underset{k=1}{\overset{\mathrm{}}{}}\frac{\alpha ^k}{k}\text{ for }\alpha \text{ small}$$
(2.26)
and then make the following change of variable

$$y=\lambda ^n.$$
(2.27)
This is possible, because if we start with a small closed curve $`\mathrm{\Gamma }`$ winding once around $`1`$, then after the variable transformation (2.27) we still have a closed curve winding once around $`1`$. We restrict ourselves to the case $`f(\lambda )=\lambda ^s`$.

$`{\displaystyle \frac{s}{2\pi i}}{\displaystyle _\mathrm{\Gamma }}𝑑\lambda \text{Tr}\mathrm{log}\left(1ϵ{\displaystyle \frac{_{l=0}^{n1}a_l\lambda ^l}{\lambda ^n1}}\right)\lambda ^{s1}=`$
$`s{\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{ϵ^k}{2\pi ik}}{\displaystyle _\mathrm{\Gamma }}𝑑\lambda \text{Tr}\left({\displaystyle \frac{_{l=0}^{n1}a_l\lambda ^l}{\lambda ^n1}}\right)^k\lambda ^{s1}=`$
$`{\displaystyle \frac{s}{n}}{\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{ϵ^k}{2\pi ik}}{\displaystyle _\mathrm{\Gamma }}𝑑y\text{Tr}\left({\displaystyle \frac{_{l=0}^{n1}a_ly^{\frac{l}{n}}}{y1}}\right)^ky^{\frac{s}{n}1}.`$
Notice that $`y^{\frac{l}{n}}`$ is regular in a neighbourhood of $`1`$. In this form we can already see that the result is symmetrized in the coefficients $`\{a_i\}`$, because they enter through expressions of the type

$$\left(\underset{l=0}{\overset{n1}{}}a_ly^{\frac{l}{n}}\right)^k=\underset{n_0+\mathrm{}+n_{n1}=k}{}\frac{k!}{n_0!\mathrm{}n_{n1}!}𝒮(a_0^{n_0}\mathrm{}a_{n1}^{n_{n1}})y^{{\displaystyle \frac{1}{n}}{\displaystyle \underset{l=1}{\overset{n1}{}}}ln_l}.$$
(2.29)
We obtain

$`\text{Tr}x^s`$
$`=\text{Tr}\mathrm{\hspace{0.25em}1}+{\displaystyle \frac{s}{n}}{\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}ϵ^k{\displaystyle \underset{n_0+\mathrm{}+n_{n1}=k}{}}{\displaystyle \frac{\text{Tr}𝒮(a_0^{n_0}\mathrm{}a_{n1}^{n_{n1}})}{n_0!\mathrm{}n_{n1}!}}{\displaystyle \frac{d^{k1}}{dy^{k1}}}y^{{\displaystyle \frac{1}{n}}\left({\displaystyle \underset{l=1}{\overset{n1}{}}}ln_l+sn\right)}|_{y=1}`$
$`=\text{Tr}\mathrm{\hspace{0.25em}1}+s{\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{ϵ^k}{n^k}}{\displaystyle \underset{n_0+\mathrm{}+n_{n1}=k}{}}{\displaystyle \frac{\text{Tr}𝒮(a_0^{n_0}\mathrm{}a_{n1}^{n_{n1}})}{n_0!\mathrm{}n_{n1}!}}{\displaystyle \underset{r=1}{\overset{k1}{}}}\left(s+{\displaystyle \underset{l=1}{\overset{n1}{}}}ln_lrn\right).`$
Here we have applied the Cauchy theorem in its more general form

$$(2\pi i)^1_C𝑑y\frac{f(y)}{(yy_0)^k}=\frac{1}{k1!}\frac{d^{k1}}{dy^{k1}}f(y)|_{y=y_0},$$
(2.31)
where $`C`$ a closed curve winding once around $`y_0`$, and $`f(y)`$ is a function which is regular inside $`C`$.

In the simplest case $`n=2`$ it is possible to solve the classical equation (i.e. the equation for commuting $`a_i`$) and this provides a closed expression for (2.2)

$$\text{Tr}x=\text{Tr}\left[ϵ\frac{a_1}{2}+𝒮\sqrt{1+ϵa_0+(ϵ\frac{a_1}{2})^2}\right],$$
(2.32)
where we again have chosen the solution which reduces to $`1`$ for $`ϵ=0`$. In the appendix we show how expanding the square root in (2.32) near $`1`$ we actually recover (2.2) in this case.

## 3 Relation between two unilateral matrix
equations

The same procedure which we have used to prove (2.2), namely applying the generalized Bezout theorem and then the Cauchy theorem, can be equally well followed for the equation

$$\mathrm{\Phi }=A_0+A_1\mathrm{\Phi }+\mathrm{}A_n\mathrm{\Phi }^n,$$
(3.1)
which is studied in . If we define

$$A(\lambda )=A_0+A_1\lambda +\mathrm{}A_n\lambda ^n$$
(3.2)
the characteristic polynomial of (3.1) is $`\lambda A(\lambda )`$ and the result which corresponds to (2.25) is

$$\text{Tr}f(\mathrm{\Phi })=\frac{1}{2\pi i}_C𝑑\lambda \text{Tr}\mathrm{log}(1\frac{A(\lambda )}{\lambda })f^{}(\lambda ).$$
(3.3)
Here $`C`$ is a closed curve winding once around $`0`$. We consider the particular case $`f(\lambda )=\lambda ^s`$ and expand the logarithm

$$\text{Tr}\mathrm{\Phi }^s=s\underset{k=1}{\overset{\mathrm{}}{}}\frac{(k1)!}{2\pi i}\underset{n_0+\mathrm{}n_n=k}{}\frac{\text{Tr}𝒮(A_0^{n_0}\mathrm{}A_n^{n_n})}{n_0!\mathrm{}n_n!}_C𝑑\lambda \lambda ^{_{l=0}^n(l1)n_l+s1}.$$
(3.4)
The Cauchy theorem has the effect of selecting the words of dimension $`_{l=0}^n(l1)n_l=s`$ in the expansion of the logarithm, i.e.

$$\text{Tr}\mathrm{\Phi }^s=s\text{Tr}\underset{k=1}{\overset{\mathrm{}}{}}\frac{1}{k}(A_0+\mathrm{}A_n)^k|_{_{l=0}^n(l1)n_l=s}$$
(3.5)
and we recover the result already stated in , where the concept of dimension of a word was introduced. It should be remarked that the solution to (3.1) for vanishing coefficients is $`\mathrm{\Phi }=0`$. Consistent with that, as already noted in in , the result (3.5) only holds for positive integers $`s`$: it has no sense to invert the solution in this case, but see below.

To study more closely the relation between equation (2.1) and (3.1) we make the Ansatz

$$x=1+\alpha \mathrm{\Phi }.$$
(3.6)
(Again, the same procedure could be followed for any of the roots of unity which are the solutions of (2.1) for $`ϵ=0`$.)

We choose the parameter $`\alpha `$ in such a way that

$$\alpha ^{n1}=n.$$
(3.7)
There is no reason for $`\alpha `$ to be real, since most of the solutions of (2.1) are not real, even for real coefficients.

Then the relation between the $`n`$ coefficients $`a_i`$ and the $`n+1`$ coefficients $`A_i`$ is easily found to be

$$A_l=\{\begin{array}{cc}\alpha ^{ln}ϵ\underset{r=l}{\overset{n1}{}}\left(\begin{array}{c}r\hfill \\ l\hfill \end{array}\right)a_r\hfill & \text{ for }l=0,1\hfill \\ \alpha ^{ln}\left(\begin{array}{c}n\hfill \\ l\hfill \end{array}\right)\alpha ^{ln}ϵ\underset{r=l}{\overset{n1}{}}\left(\begin{array}{c}r\hfill \\ l\hfill \end{array}\right)a_r\hfill & \text{ for }2ln.\hfill \end{array}$$
(3.8)
Some remarks can be made with respect to (3.8). First, observe that $`A_n=1`$ is fixed, but this had to be expected, because the equation (2.1) has one coefficient less than the equation (3.1).

Moreover, (3.8) is a linear relation, and it is invertible, so that once the result that the trace of the powers of a solution depends only on the symmetrized products of the coefficients is proven for one of the two equations, it immediately follows also for the other. Negative powers of $`x`$ can be expanded from (3.6) into a series of positive powers of $`\mathrm{\Phi }`$.

If $`ϵ`$ is small, only the first two coefficients $`A_0`$ and $`A_1`$ are automatically small, the other coefficients satisfy

$$A_l\stackrel{ϵ0}{}\alpha ^{ln}\left(\begin{array}{c}n\hfill \\ l\hfill \end{array}\right)\text{ for }2ln.$$
(3.9)
Therefore, the two expansions (2.2) and (3.5) do not necessarily hold for the same range of the coefficients.

## Appendix

In this appendix we check that (2.2) actually is the series expansion of (2.32) for $`s=1`$, $`n=2`$. We start by expanding the square root in (2.32) around $`1`$

$`\text{Tr}x`$ $`=`$ $`\text{Tr}\left[1+ϵ{\displaystyle \frac{a_1}{2}}+{\displaystyle \underset{r=1}{\overset{\mathrm{}}{}}}(1)^{r1}{\displaystyle \frac{1}{2^rr!}}(2r3)!!𝒮\left(ϵa_0+(ϵ{\displaystyle \frac{a_1}{2}})^2\right)^r\right]`$
$`=`$ $`\text{Tr}\mathrm{\hspace{0.25em}1}+ϵ{\displaystyle \frac{\text{Tr}a_1}{2}}`$
$`+{\displaystyle \underset{r=1}{\overset{\mathrm{}}{}}}(1)^{r1}{\displaystyle \frac{1}{2^r}}(2r3)!!{\displaystyle \underset{n_0=0}{\overset{r}{}}}{\displaystyle \frac{ϵ^{2rn_0}}{n_0!(rn_0)!}}\text{Tr}𝒮\left(a_0^{n_0}\left({\displaystyle \frac{a_1}{2}}\right)^{2(rn_0)}\right).`$
We introduce the new variables:

$$n_1=2(rn_0),k=n_0+n_1=2rn_0.$$
(A.2)
Then we can rewrite

$`\text{Tr}x`$ $`=`$ $`\text{Tr}\mathrm{\hspace{0.25em}1}+ϵ{\displaystyle \frac{\text{Tr}a_1}{2}}`$
$`+{\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{ϵ^k}{2^k}}{\displaystyle \underset{\stackrel{n_0+n_1=k}{n_1even}}{}}(1)^{\frac{n_1}{2}+n_01}{\displaystyle \frac{(n_11)!!}{n_0!n_1!}}(2n_0+n_13)!!\text{Tr}𝒮(a_0^{n_0}a_1^{n_1}),`$
where we have used the relation

$$\frac{n_1!}{\frac{n_1}{2}!2^{\frac{n_1}{2}}}=(n_11)!!\text{ for }n_1\text{ even.}$$
(A.4)
Due to the conditions $`n_0+n_1=k`$, $`n_1`$ even, the following relation holds

$$(1)^{\frac{n_1}{2}+n_01}(n_11)!!(2n_0+n_13)!!=\underset{r=1}{\overset{k1}{}}\left(1+n_12r\right)$$
(A.5)
and so formula (Appendix) can be brought into the more compact form

$$\text{Tr}x=\text{Tr}\mathrm{\hspace{0.25em}1}+\underset{k=1}{\overset{\mathrm{}}{}}\frac{ϵ^k}{2^k}\underset{n_0+n_1=k}{}\frac{1}{n_0!n_1!}\text{Tr}𝒮(a_0^{n_0}a_1^{n_1})\underset{r=1}{\overset{k1}{}}\left(1+n_12r\right),$$
(A.6)
which now coincides with (2.2) for $`s=1`$, $`n=2`$. It is no longer necessary to explicitly sum only over even values of $`n_1`$, because

$$\underset{r=1}{\overset{k1}{}}\left(1+n_12r\right)=0\text{for }n_1\text{ odd, }1n_1k,k>1.$$
(A.7)
## Acknowledgments

We would like to thank Dan Brace and Bogdan Morariu for many helpful discussions. This work was supported in part by the Director, Office of Science, Office of High Energy and Nuclear Physics, Division of High Energy Physics of the U.S. Department of Energy under Contract DE-AC03-76SF00098 and in part by the National Science Foundation under grant PHY-95-14797. B.L.C. is supported by the DFG (Deutsche Forschungsgemeinschaft) under grant CE 50/1-1.


text:
# HUDDLED MASSES YEARNING TO STREAM FREE: GLOBULAR CLUSTERS IN THE HEARTS OF GALAXY CLUSTERS

## 1 Introduction: Globular Clusters and their Uses

Globular clusters (GCs) provide key information for a wide variety of areas in astrophysics and cosmology. Absolute GC ages (see reviews by Salaris et al.  and Chaboyer et al. ) constrain the cosmological model by setting a lower bound on the age of the universe. The observed GC age range (e.g., Stetson et al. 1996; Sarajedini et al. 1997) constrains the time-scale of galaxy formation, and the distribution in metallicity can provide further details on the formation processes (e.g., Ashman & Zepf 1992; Ajhar et al. 1994; Forbes et al. 1997).

Globular clusters are also an important part of the distance scale. Their horizontal-branch magnitude is the basis for the “Population II” RR Lyra distances within our Galaxy (e.g., Carney et al. 1992; Carretta et al. 2000) and provides an important self-consistency check (e.g., Walker 1992; Ajhar et al. 1996; Fusi Pecci et al. 1996) of the Cepheid distance scale on which most extragalactic distances are based. In addition, the mean luminosity of the GCs in external galaxies has been used extensively as a standard-candle distance indicator reaching as far as the Coma cluster (e.g., Harris 1991; Baum et al. 1995, 1997; Whitmore et al. 1995; Kavelaars et al. 2000; Ferrarese et al. 2000).

Besides these other, more traditional, “uses,” globular clusters are now recognized as unique and important probes of the centers of galaxy clusters. The present paper summarizes some recent work on this topic.

## 2 Globular Clusters as Kinematical Tracers of the Cluster Potential

The existence of massive dark halos around spiral galaxies has been known for over 20 years (Rubin et al. 1978) from their flat rotation curves. The presence of a dominant dark matter component in galaxy clusters has been known even longer from the galaxy dynamics (Zwicky 1937), and repeatedly confirmed by X-ray, lensing, and other dynamical studies. Yet, massive dark halos around individual giant ellipticals have been notoriously difficult to demonstrate and measure. The stellar velocity dispersion can only be traced to $`\mathrm{\hspace{0.17em}10}`$ kpc, roughly 20% of the distance to which spiral rotation curves can be followed, and then it becomes necessary to resort to various, usually sparse, tracers of the halo potential. For the same reason, the mass distribution in cluster cores is poorly constrained at radii $`r10`$–50 kpc.

A major step forward came with the work of Cohen & Ryzhov (1997), who used the 10-m Keck telescope to measure velocities for 205 confirmed GCs around M87, the giant elliptical at the dynamical center of the Virgo cluster. Velocities for an additional 16 GCs were added by Cohen (2000). These authors found that the 1-d velocity dispersion of the M87 GC population rises from $`\sigma 300`$ km s-1 at $`r6`$ kpc to $`\sigma >450`$ km s-1 at $`r>30`$ kpc (assuming the Virgo Cepheid distance of 16 Mpc, so that $`1^{}=4.7`$ kpc). Interestingly, the rotation also increases outward from the center, reaching $`v_{\mathrm{rot}}300`$ km s-1 (Cohen 2000). Cohen & Ryzhov concluded that the total mass increases as $`M(<r)r^{1.7}`$, or $`\rho r^{1.3}`$, which agrees well with the mass profile inferred from X-ray data (Nulsen & Böhringer 1995) in the region of overlap. A similar result was found by Romanowsky & Kochanek (2000), who determined the mass profile using joint constraints from the GCs and halo stars and concluded that the potential of the Virgo cluster is dominant by $`r20`$ kpc. Thus, the GC velocity data trace the transition from the dynamics of the M87 stellar halo to that of surrounding Virgo cluster.

Kissler-Patig et al. (1999) studied a sample of 74 GC velocities around NGC 1399, the central cD galaxy in the Fornax cluster ($`d19`$ Mpc), using data drawn from three previous spectroscopic studies. They also found a rising velocity dispersion profile, bridging the “velocity gap” between the NGC 1399 stellar halo and the galaxy population of the Fornax cluster. These authors concluded that most of the extended population of GCs were stripped from other Fornax galaxies, as advocated by Côté et al. (1998).

Recently, a team led by R. Sharples at Durham and S. Zepf at Yale (Sharples et al. 1997; Zepf et al. 2000) and another from Caltech (J. Cohen, P. Côté, and myself) have obtained velocity data for large samples of GCs around NGC 4472. Although NGC 4472 is not a central cluster galaxy, it is the brightest Virgo member and the center of its own smaller subcluster $``$1 Mpc away from M87 in projection. It has also been the subject of a very detailed photometric study using Washington photometry (Geisler et al. 1996). The GC system contains two unusually well demarcated metallicity subpopulations, as compared to most other giant ellipticals (e.g., Gebhardt & Kissler-Patig 1999; Woodworth & Harris 2000; Brodie et al. 2000).

Figure 1 compares the velocity dispersion data for 272 NGC 4472 GCs, split into blue and red subpopulations following Geisler et al. (1996), to that of the stellar halo and the satellite galaxy population. The red GCs share the kinematics of, and have colors similar to, the stellar halo, bespeaking a common origin. The blue GCs have a significantly higher dispersion and show the transition to the dynamics of galaxy population. The blue population is also more extended in its spatial distribution. However, both subpopulations appear equally old, within the uncertainties (Beasley et al. 2000), and the same is true in M87 (Cohen et al. 1998).

Although it appears universally true that blue GC populations have more spatially extended distributions, it is unclear whether or not the striking difference in the dynamics of the two subpopulations is a common feature of the GCs around central cluster galaxies. There is no evidence for this effect in the available M87 GC data (Cohen et al. 1998), but the sample was not so clearly separated into blue and red components. In NGC 1399, the sample size is too small in to constitute a reliable test. Further spectroscopic observations in these two important systems are needed. Velocity measurements for the brightest GCs around the Hydra cluster cD NGC 3311, a giant galaxy with a very diffuse cD halo and a broad GC color distribution lacking any apparent multiple peaks (Brodie et al. 2000), would be extremely valuable and is now feasible with the VLT.

## 3 Globular Cluster Number as a Mass Tracer

At distances much beyond the Virgo and Fornax clusters, it becomes increasingly difficult to measure accurate velocities for individual GCs. For instance, at $`cz3500`$ km s-1, the very brightest GCs have an apparent magnitude $`V22.6`$, requiring very long integrations to obtain adequate signal in the absorption line spectra. However, abundant GC populations can be photometrically detected around bright galaxies to many times this distance, and their properties contain useful information on the host galaxies and galaxy clusters.

### 3.1 Measuring GC Specific Frequencies

The globular cluster “specific frequency” $`S_N`$ of a galaxy is the number of GCs per unit $`M_V=15`$ of galaxy luminosity (Harris & van den Bergh 1981). Thus, a dwarf galaxy with $`M_V=`$15 and a single GC would have $`S_N=\mathrm{\hspace{0.17em}1}`$. Spiral galaxies such as the Milky Way or M31 typically have $`S_N=0.5`$–1, while most ellipticals have $`S_N=2`$–4. However, cD galaxies such as M87 can have $`S_N12`$, which means a total GC population of $`N_{\mathrm{GC}}12`$,500 (see the review by Harris 1991). However, data on other such “high-$`S_N`$” systems were very sparse, as the total number of GCs is difficult to estimate from GC counts if the data do not approach the turnover magnitude $`m^0`$ of the globular cluster luminosity function (GCLF). At the distance of the Coma cluster, the GCLF has a $`V`$-band turnover $`m_V^027.5`$.

To address this problem, we (Blakeslee & Tonry 1995) developed a new method for simultaneously measuring specific frequencies and GCLF widths in relatively distant galaxies without the requirement of complete point source detection to the GCLF turnover. The method borrows from the surface brightness fluctuations (SBF) method (Tonry & Schneider 1988) for measuring galaxy distances. In the standard SBF technique, the fluctuations produced by the Poisson statistics of stars in an early-type galaxy are used to derive an average flux for the stellar population, yielding the distance (see the recent review by Blakeslee, Ajhar, & Tonry 1999). Here, we apply the same analysis methods to measure the fluctuations from unresolved GCs around more distant galaxies.

The usual way of determining $`S_N`$ in galaxies beyond Virgo is to count the brightest GCs and extrapolate, yielding a result which is strongly dependent upon the assumed GCLF width $`\sigma `$. In our technique, we also first identify and count the brightest GCs, determining $`S_N`$ as a function of $`\sigma `$, but we then excise these bright ones and measure the psf-convolved variance (“fluctuations”) in the image from the rest of the GC population. This variance is directly proportional to the number density of remaining GCs; thus, we again determine $`S_N`$ as function of $`\sigma `$. These two determinations will be consistent over some range in $`\sigma `$ which should encompass the true value. Figure 2 provides an illustration. Even if there is relatively little variation in the GCLF $`\sigma `$ among cD galaxies, the combined information from the counts of the brightest GCs and the variance from the fainter ones produces much tighter constraints on the total GC number.

### 3.2 Ground-based Suveys and Interpretations

We studied two different samples of GC populations in Abell cluster central galaxies using this new method. The first was a statistically complete sample of 19 northern clusters with early-type brightest members within $`cz<10,000`$ km s-1 observed with the 2.4-m telescope at MDM Observatory on Kitt Peak (Blakeslee et al. 1997). The second was a smaller sample of more massive clusters to twice this distance with Keck (Blakeslee 1999). In all, $`S_N`$ values were measured for 32 bright galaxies in 25 clusters. Figure 3 shows data for the most distant cluster studied, A2124 at $`z=0.066`$, which also happens to be the nearest cluster in which strong gravitational lensing of distant galaxy has been observed (Blakeslee & Metzger 1999)

The data show that $`S_N`$ for central cluster galaxies correlates well with overall cluster properties including velocity dispersion and X-ray luminosity (see Figure 4), as well as other measures of cluster richness. In fact, the luminosities of the cD galaxies show relatively little variation (e.g., Postman & Lauer 1995), so the correlations exhibited in Figure 4 are really between the total number of GCs in the cluster center and the mass of the cluster. Luminous galaxies displaced from the cluster dynamical centers do not seem to participate in these correlations (although see Woodworth & Harris 2000).

We have argued (Blakeslee et al. 1997; Blakeslee 1999) that the observations are most easily explained if the GCs surrounding cD galaxies formed at early times and with an approximately universal efficiency of roughly $`0.7`$GC per $`10^9M_{}`$. However, the present-day galaxy luminosity is fairly insensitive to this quantity, resulting in the observed correlations of GC frequency with cluster density. Thus, rather than having anomalously large GC populations, the galaxies are simply underluminous for their prominent positions in the centers of very rich clusters. In this “missing light” scenario, the gas is removed from the central galaxy into the intracluster medium at the cluster formation epoch, halting any subsequent star formation in the central galaxy that would serve to lower $`S_N`$. A similar model has been discussed by Harris et al. (1998), and McLaughlin (1999) has shown that the numbers work out in detail for the amount of gas surrounding the nearby galaxies M87, NGC 1399, and NGC 4472.

Blakeslee (1999) suggested that there may be evidence for later $`S_N`$ evolution in these systems: the central galaxies in the more spiral-rich, less evolved clusters had marginally higher GC specific frequencies than those in the more centrally concentrated, spiral-poor clusters. Again the difference in the $`S_N`$ values was because of the lower luminosities of the central galaxies in the less evolved clusters, not because of larger GC populations of in the more evolved one. Thus, continuing star formation and dynamical evolution in the spiral-rich clusters may allow the central galaxy to grow in luminosity, without adding large numbers of new GCs; the net result is a decrease in GC frequency. Following McLaughlin et al. (1994), who made a similar suggestion from a smaller, much less homogeneous, dataset, we call this process “$`S_N`$ dilution.” It will be interesting to see if this effect proves true; for now, the evidence is inconclusive.

### 3.3 Recent HST Data

Of course, with HST, it is possible to reach the GCLF turnover at the distance of the Coma cluster and obtain a more direct estimate of $`S_N`$. No one has yet done a systematic study of GC populations in a large sample of cD galaxies in rich clusters using HST, as the required integration times are substantial. However, Harris et al. (2000) have recently used HST to measure $`S_N4`$ for NGC 4874 in Coma. This is significantly lower than our value of $`S_N=9\pm 2`$ (Blakeslee et al. 1997). There are two basic reasons for this. First, they assumed a distance modulus 0.35 mag greater than ours, and this translates directly into a $`35`$% difference in $`S_N`$. Harris et al. point out that there is very close agreement between the observed GC counts when compared over the same radial and magnitude ranges. For instance, if the value of $`S_N`$ at $`r1^{}`$ shown in Figure 2 above is transformed to their distance scale, we obtain very close agreement with $`S_N4`$.

The more important difference stems from the radial profiles measured in the two studies. We found that the GC distribution was more extended than the halo light, so that $`S_N`$ increased with radius, yielding a significantly higher $`S_N`$. Harris et al. also found a very extended distribution in the inner parts, but then a sharp truncation of the GC system at $`r\mathrm{\hspace{0.17em}2}`$arcmin. It is most likely that problems with background source contamination or the $`I`$-band surface photometry at large radii (the rest of our data were all taken in the $`R`$-band, where the sky is much darker and CCD’s flatten better, to avoid this problem) have affected the ground-based $`S_N`$ measurement. Still, the sharp truncation found near the limit of the HST data is difficult to understand. It would be worth revisiting this system with the next generation of HST imaging cameras (ACS, WFC3) and their much wider fields of view.

## 4 Summary

The rich globular cluster populations at the centers of galaxy clusters constitute important probes of the central cluster potential. In nearby clusters, they can be used to directly trace the transition from the potential of the central galaxy to that of surrounding cluster. In the case of NGC 4472, the brightest Virgo elliptical, this transition is well traced by the more spatially-extended blue GC subpopulation, whereas the red GCs appear to share the dynamics of the stellar halo. By obtaining large samples of GC velocities around both central and non-central giant ellipticals, it will be possible to determine the mass profiles of the galaxies and constrain the slope of the dark matter profile near the centers of galaxy clusters. This will provide a useful comparison for high-resolution $`N`$-body simulations (e.g., Moore et al. 1998).

At larger distances where GCs become too faint for spectroscopic measurements, observations of GC populations around cD galaxies can give a reasonable indication of the cluster core mass. The reason for the elevated specific frequencies of the central galaxies is most likely that their gas was removed and their star formation halted sometime after the formation of the GCs. In this case, it should be possible to observe $`S_N`$ evolution in clusters of different evolutionary states, but the available data are inconclusive. As always, future observations of GC systems from the ground and with HST will be helpful in resolving the many remaining questions.

## Acknowledgments

I thank Mike Beasley for supplying Figure 1 and helpful comments. This work was supported at Caltech by a Fairchild Fellowship and at the University of Durham by a PPARC Rolling Grant in Extragalactic Astronomy and Cosmology.