Enhancing IMR Consistency Tests for Gravitational Wave Transients

Time-Domain Analysis

The time domain likelihood function for gravitational wave time series data $\textbf{d} = [d_i]$ and waveform model $\textbf{h} = [h_i]$ with Gaussian noise is expressed as $\mathcal{L}(\textbf{d} | \textbf{h}) = \frac{1}{\sqrt{2\pi\det \textbf{C}}} \exp{\left(-\frac12(\textbf{d} - \textbf{h})^\dag \cdot \textbf{C}^{-1} \cdot (\textbf{d} - \textbf{h})\right)}$ where $\textbf{C}$ is the noise correlation matrix defined by $C_{ij} = \big\langle n_in_j \big\rangle$, where $n_k=d_k-h_k$, as the ensemble average.

We analyze the noise covariance matrix, $C_{ij}=\lim\limits_{n\rightarrow\infty}\frac{1}{n}\displaystyle\sum\limits_{\ell=0}^{n-1}s^{\ell}_is^{\ell}_{\Delta_{ji}+i}$, where $\Delta_{ji}=j-i$, which usually is dependent on time $t_i=i\Delta t$ and displacement $\tau_{ij}=\Delta_{ji}\Delta t$. However, one can make the assumption that the noise is wide-sense stationary, where the mean and variance are both constant in time. Under this assumption, any constant can be added to the indices in the above equation to obtain the same result. This makes $\textbf{C}$ symmetric since the factors commute, and Toeplitz since the elements along the diagonals are equal. Additionally, since $C_{0\ \Delta_{ji}}=C_{-\Delta_{ji}\ 0}=C_{0\ -\Delta_{ji}}$, the elements of $\textbf{C}$ are even functions of $\Delta_{ji}$.

Additionally, one can assume that the data is ergodic, meaning that new realizations of $\textbf{n}$ are obtained via time. Under this assumption and the properties of the elements of $\textbf{C}$, the ensemble averages can be replaced with time averages to give, $C_{ij}=\lim\limits_{n\rightarrow\infty}\frac{1}{n}\displaystyle\sum\limits_{\ell=0}^{n-1}s^{\ell}_{\ell}s^{\ell}_{\Delta_{ji}+\ell}=\lim\limits_{n\rightarrow\infty}\frac{1}{2n}\displaystyle\sum\limits_{\ell=-n}^{n-1}s^{\ell}_is^{\ell}_{\Delta_{ji}+i}=\frac{1}{2}\rho((j-i)\Delta t)$, where $\rho$ describes the correlation between points in the time series.

For a stationary random process, the covariance takes a symmetric Toeplitz form given by $C_{ij} =\frac{1}{2}\rho((j-i)\Delta t)$ where $\rho$ is the autocovariance function. This is estimated empirically by autocorrelating a long stretch of noise-only data, where $\rho(k) \approx \frac1{N_{\rho}}\sum\limits_{i=0}^{N_{\rho}-1}n_in_{i+k} \quad (N_{\rho}\gg N)$ for $|k|<N_{\rho}$ and zero otherwise. If, in addition to stationarity, we impose periodic boundary conditions, then $\rho(k) = \rho(N - k)$, and $C$ will be circulant. Circulant matrices are diagonalized by the discrete Fourier transform, implying the Fourier amplitudes are given by $\tilde{C}_{ij} = \big\langle \tilde{n}_i\tilde{n}_j \big\rangle = \frac{1}{2} (N\Delta t) S(|f_{ij}|) \delta_{ij}$ where the PSD is derived from the cyclic autocorrelation function through a discrete Fourier transform $S(|f_j|) = 2\Delta t \sum\limits_{j=0}^{N-1} \rho(k)\exp{-\frac{2\pi ijk}{N}}$.

Starting from a given PSD estimate $S(|f|)$ sampled at $N$ frequencies $f_j$ , we can invert the above equation to obtain the corresponding autocorrelation function estimate is given by $\rho(k) = \frac{1}{2T} \sum\limits_{k=0}^{N_{\text{FFT}}-1} S(|f_j|)\exp{\frac{2\pi ijk}{N}}$ as long as the PSD was estimated from data segments of length $N_{\text{FFT}}$ much longer than the analysis segment length $N$. A long segment allows us to truncate the resulting $N_{\text{FFT}}-$long $\rho(k)$ to length $N$ before constructing the covariance matrix. Thereby the covariance matrix can be estimated as a symmetric Toeplitz form.

Once we have constructed a covariance matrix, we can analyze a noisy data stream with the likelihood is evaluated as $\ln\mathcal{L}(\textbf{d}|\textbf{h}) = -\frac{1}{2}\sum\limits_{i,j=0}^{N-1}\left(\tilde{d}_i-\tilde{h}_i\right)C^{-1}_{ij}\left(\tilde{d}_j-\tilde{h}_j\right) = \sum\limits_{i,j=0}^{N-1}\tilde{d}_iC^{-1}_{ij}\tilde{h}_j-\frac12\sum\limits_{i,j=0}^{N-1}\tilde{d}_iC^{-1}_{ij}\tilde{d}_j=\mathrm{SNR}^{\mathrm{inner}}-\frac12\mathrm{SNR}^{\mathrm{opt}}$

Gating and In-Painting analyses

If we assume that a detector’s noise is wide-sense stationary and ergodic, only then its covariance $\textbf{C}$ is a symmetric Toeplitz matrix with elements given by the autocorrelation function of the data. If the autocorrelation function goes to zero in some finite amount of time that is less than $T/2$, then the covariance matrix is asymptotically equivalent to a circulant matrix. The inverse of the covariance matrix can then be well-approximated by that of the equivalent circulant matrix.

Instead of considering the full set of time samples, if we wish to only evaluate the truncated set $\textbf{d}_{\rm tr} = \{d_0, ..., d_a, d_{a+M}, ..., d_{N-1}\}$, with $M > 1$. The data between the time steps $[a,a+M)$ is said to be gated.

The probability density function of the truncated noise is still a multivariate normal distribution (excising dimensions from a multivariate normal is equivalent to marginalizing over those dimensions). The challenge is that the covariance matrix of the truncated noise $\textbf{C}_{\rm tr}$ is no longer Toeplitz. Its eigenvectors can no longer be approximated by that of a circulant matrix, and so the above expression for the likelihood is no longer valid. The inverse of the covariance matrix needs to be found by other means. We use gating and in-painting to find the likelihood of the truncated time series. This was applied to the problem of matched filtering; we apply it to parameter estimation.

Define $\textbf{n}' = \textbf{n}_g + \textbf{x}$, where $\textbf{n}_g$ is the noise with the gated times $t\in [a, a+M)\Delta t$ zeroed out, and $\textbf{x}$ is a vector that is zero everywhere except in the gated times. If the non-zero elements of $\textbf{x}$ are such that $(\textbf{C}^{-1} \textbf{n}')[k] = 0$ for all $k \in [a, a+M)$, then $\textbf{n}'^\dag \textbf{C}^{-1} \textbf{n}'$ will be the same as the truncated version $\textbf{n}_{\rm tr}^\dag \textbf{C}_{\rm tr}^{-1} \textbf{n}_{\rm tr}$. Our aim is to solve the equation $\textbf{C}^{-1}(\textbf{n}_g + \textbf{x}) = 0$ in the gated region. Since $\textbf{x}$ is zero outside of the gated region, $\textbf{C}^{-1}x$ only involves the $[a, a+M)$ rows and columns of $\textbf{C}^{-1}$, which form an $M \times M$ Toeplitz matrix. We therefore solve for $\textbf{x}$ such that $\overline{\textbf{C}^{-1}}\,\overline{\textbf{x}} = -\overline{\textbf{C}^{-1}\textbf{n}_g}$, where the overbar indicates the $[a, a+M)$ rows (and columns) of the given vector (matrix). This can be solved numerically using a Toeplitz solver. Adding $\textbf{x}$ to the gated noise (in painting) will then yield the same result as if we had truncated the noise and the covariance matrix.

Note that if the gate spans the entire beginning of the data segment, the truncated covariance matrix $\textbf{C}_{\rm tr}$ is Toeplitz, and so could be inverted numerically using a Toeplitz solver. The advantage of using in-painting is that it involves solving for an $M\times M$ matrix rather than an $(N-M)\times (N-M)$ matrix. In the standard in-painting approach, one first solves for the in-painting vector $\textbf{x}$ using a Toeplitz solver and then augments the gated noise $\textbf{n}_g$ by setting $\textbf{n}' = \textbf{n}_g + \textbf{x}$.

Drafting is currently in progress. Sorry for any inconvenience. Please check back soon!