Global analysis for the LISA gravitational wave observatory

TL;DR

This paper presents a global data analysis approach for LISA, accounting for overlapping signals and noise, and estimates how confusion noise and parameter errors are affected by this comprehensive method.

Contribution

It introduces analytic estimates for signal subtraction residuals and noise impacts in LISA's global analysis, incorporating noise modeling for the first time.

Findings

Confusion noise is reduced in global analysis.

Waveform errors for sources increase with global analysis.

Parameter estimation errors are inflated due to overlapping signals.

Abstract

The Laser Interferometer Space Antenna (LISA) will explore the source-rich milli-Hertz band of the gravitational wave spectrum. In contrast to ground based detectors, where typical signals are short-lived and discrete, LISA signals are typically long-lived and over-lapping, thus requiring a global data analysis solution that is very different to the source-by-source analysis that has been developed for ground based gravitational wave astronomy. Across the LISA band, gravitational waves are both signals {\em and} noise. The dominant contribution to this so-called confusion noise (better termed unresolved signal noise) is expected to come from short period galactic white dwarf binaries, but all sources, including massive black hole binaries and extreme mass ratio captures will also contribute. Previous estimates for the galactic confusion noise have assumed perfect signal subtraction. Here we provide analytic estimates for the signal subtraction residuals and the impact they have on parameter estimation while for the first time incorporating the effects of noise modeling. The analytic estimates are found using a maximum likelihood approximation to the full global Bayesian analysis. We find that while the confusion noise is {\em lowered} in the global analysis, the waveform errors for individual sources are {\em increased} relative to estimates for isolated signals. We provide estimates for how parameter estimation errors are inflated from various parts of a global analysis.

Figures

Tables

Table 1. Table 1: Parameter Error Estimates

Parameter	${\sigma }_{\Gamma }$	${\sigma }_{\mathrm{MCMC}}$	$\Delta {\lambda }_{\mathrm{ML}}$	$\Delta {\lambda }_{\mathrm{MAP}}$	$\Delta {\lambda }_{\mathrm{mean}}$
$A$	0.141	0.135	-0.078	-0.048	-0.067
$f_0$	0.00055	0.00054	0.0010	0.0011	0.0010
$t_0$	0.125	0.123	0.195	0.214	0.217

Equations

$$S_p(f)$$

$$S_a(f)$$

$$A$$

$$E$$

$$T$$

$$p(\text{s}|\text{h})=e^{-G/2}{e}^{-(\text{s}-\text{h}(\vec{\lambda })|\text{s}-\text{h}(\vec{\lambda }))/2}=e^{-G/2}{e}^{-{\chi }^2/2},$$

$$\left(\text{g}|\text{k}\right)\text{vbox}..=2\sum _{I=\{A,E\}}{\int }_0^{\infty }\frac{{\tilde{g}}_I(f){\tilde{k}}_I^{\ast }(f)+{\tilde{g}}_I^{\ast }(f){\tilde{k}}_I(f)}{S_{n,I}(f)}df,$$

$$\text{mathbbm}E[{\tilde{n}}_I(f){\tilde{n}}_J^{\ast }(f^{\prime })]=\frac{1}{2}\delta (f-f^{\prime })S_{n,I}(f){\delta }_{IJ}.$$

$$G=\sum _{I=\{A,E\}}{\int }_0^{\infty }T\log \left[\pi T{S}_{n,I}(f)\right]df.$$

$$\text{h}(\vec{\lambda })={\text{h}}_T+{\partial }_i{\text{h}}_T\Delta {\lambda }^i+\mathcal{O}\left(\Delta {\lambda }^2\right),$$

$${\chi }^2=(\text{s}-\text{h}|\text{s}-\text{h})$$

$$=(\text{n}|\text{n})-2(\text{n}|{\partial }_i{\text{h}}_T)\Delta {\lambda }^i+({\partial }_i{\text{h}}_T|{\partial }_j{\text{h}}_T)\Delta {\lambda }^i\Delta {\lambda }^j+\mathcal{O}(\Delta {\lambda }^3).$$

$$\Delta {\lambda }^j=(\text{n}|{\partial }_i{\text{h}}_T){\left({\Gamma }^{-1}\right)}^{ij}+\dots$$

$${\Gamma }_{ij}=({\partial }_i{\text{h}}_T|{\partial }_j{\text{h}}_T)$$

$$C^{ij}=\text{mathbbm}E[\Delta {\lambda }^i\Delta {\lambda }^j]={\left({\Gamma }^{-1}\right)}^{ij}+\mathcal{O}{\left(\mathrm{SNR}\right)}^{-1},$$

$$\Delta \text{h}={\text{h}}_T-\text{h}\simeq -{\partial }_i{\text{h}}_T\Delta {\lambda }^i+\dots$$

$$\text{mathbbm}E[{\rho }_{\Delta h}^2]\text{vbox}..=\text{mathbbm}E[(\Delta \text{h}|\Delta \text{h})]=({\partial }_i{\text{h}}_T|{\partial }_j{\text{h}}_T)\text{mathbbm}E[\Delta {\lambda }^i\Delta {\lambda }^j]\approx {\Gamma }_{ij}{\left({\Gamma }^{-1}\right)}^{ij}=D.$$

$$\text{mathbbm}E[{\chi }^2]$$

$$\text{mathbbm}E[(\Delta \text{h}|\text{n})]=-\text{mathbbm}E[({\text{h}}_{,i}|\text{n})({\Gamma }^{-1}{)}^{ij}(\text{n}|{\text{h}}_{,j})]=-D,$$

$$\text{mathbbm}E[(\text{s}(f)-\text{h}(f))(\text{s}(f^{\prime })-\text{h}(f^{\prime }){)}^{\ast }]$$

$$\frac{1}{2}\delta (f-f^{{}^{\prime }})S_{\Delta \text{h}}(f)$$

$$=\text{mathbbm}E[\Delta \tilde{\text{h}}(f)\Delta {\tilde{\text{h}}}^{\ast }(f^{\prime })]$$

$$\alpha$$

$$\beta$$

$$p(\text{s}|\text{h})=\prod _{k=1}^N\frac{1}{\sqrt{2\pi {\sigma }^2}}{e}^{-(s_k-h_k{)}^2/(2{\sigma }^2)}.$$

$${\sigma }_0^2=\frac{1}{N}\sum _{k=1}^N{n}_k^2.$$

$$\log p(\text{s}|\text{h})$$

$$\Delta {\lambda }^j$$

$$\Delta \kappa$$

$$\text{mathbbm}E[{\sigma }^2]={\sigma }_0^2\left(1-\frac{D}{N}\right).$$

$$p(\text{s}|\Delta \vec{\eta })=\frac{1}{\sqrt{2\pi \text{mbox}det{\mathbf{\Gamma }}^{-1}}}{e}^{-\frac{1}{2}{\Gamma }_{\mu \nu }\Delta {\eta }^{\mu }\Delta {\eta }^{\nu }}$$

$$\boldsymbol{\Gamma}=\left(\begin{array}[]{cc}\Gamma_{0,ij}&(\textbf{n}|\partial_{j}\textbf{h}_{T})_{0}/\sigma_{0}^{2}\\ (\textbf{n}|\partial_{i}\textbf{h}_{T})_{0}/\sigma_{0}^{2}&2\sigma_{0}^{4}/N\end{array}\right)\,,$$

$${\left({\Gamma }^{-1}\right)}^{ii}\approx {\left({\Gamma }_0^{-1}\right)}^{ii}+\frac{2}{N}{\left({\Gamma }_0^{-1}\right)}^{il}{\left({\Gamma }_0^{-1}\right)}^{ki}(\text{n}|{\partial }_l{\text{h}}_T{)}_0(\text{n}|{\partial }_k{\text{h}}_T{)}_0.$$

$$S_n(f;\vec{\theta })=S_{n,0}(f)+\Delta {\theta }^a{\partial }_a{S}_n(f)+\frac{1}{2}\Delta {\theta }^a\Delta {\theta }^b{\partial }_a{\partial }_b{S}_n(f)+\dots ,$$

Full text

1 Introduction

The recent discovery of gravitational waves [1] and the outstanding success of the LISA pathfinder mission [2] have given new life to the LISA mission. Building on decades of study, an updated LISA mission concept [3] was recently submitted to address the European Space Agency’s “Gravitational Universe” science theme with a launch scheduled for the early 2030’s. The plan is to fly three identical spacecraft connected by six laser links forming a triangular detector with 2.5 million km long arms.

It has long been recognized that the LISA mission will suffer from “an embarrassment of riches”, delivering data sets so packed with signals that extracting information about individual sources will require the development of unique data analysis techniques. Significant attention was given to this problem through the 2000’s, culminating in a series of Mock LISA Data Challenges (MLDCs) [4, 5, 6, 7, 8, 9] that produced some promising proof-of-principle solutions. The demise of the original LISA project in 2010 halted this effort, but work is now resuming following the re-birth of the mission. In addition to finding an implementable solution to the global analysis problem, there is also interest in producing reliable estimates for the science that can be achieved, including the number of sources of each type that can be resolved, and how well they can be characterized. A key input to these studies are estimates for the confusion noise from unresolved sources, as this adds to the instrument noise, and reduces the signal-to-noise ratio (SNR) of the resolved systems. These confusion noise estimates use variants of the idealized iterative source subtraction scheme introduced by Timpano et al [10]. We recently applied this technique to produce confusion noise estimates [11] that were used in the design study for the new LISA mission concept [3]. There are, however, several deficiencies with the simple confusion noise estimates: it assume that the confusion noise is stationary when in fact it oscillates with a 6 month period; it neglects the parameter estimation errors for the subtracted signals and the waveform residuals; and it also neglects the impact the removal has on other resolvable signals such as massive black holes. Here we generalize the Timpano et al [10] method to account for the waveform residuals and the impact they have on the detection of other signals. Some of our results were derived previously by Cutler and Harms [12] in studies of foreground subtraction for the Big Bang Observer, but many results are new, including analytic estimates for power spectrum of the waveform residuals, incorporating the process of noise modeling, and the impact on parameter estimation for other sources. We find that the parameter estimation errors caused by other resolved signals are typically small compared to those due to instrument noise and unresolved signals. The exception to this rule is when two signals have very high overlap, such as sometimes occurs for galactic binaries with near identical orbital periods and sky locations [13].

Electromagnetic observations have identified $\sim 50$ galactic binaries with orbital periods that put their predicted gravitational signals in the LISA band [3]. Those that rise above the noise are refer to as “verification binaries”. Population synthesis models predict that there are far more detectable systems waiting to be discovered, though the estimates have been lowered in the past decade as on-going surveys have been used to re-calibrate the models [14, 15, 16]. It is estimated that there are hundreds of millions of galactic binaries GBs emitting gravitational waves in our galaxy. In the mid-band of the LISA sensitivity, between $\sim 0.5-3$ mHz, gravitational waves from these systems are expected to dominate over instrument noise, with the unresolved component producing what is termed “confusion noise”.

There have been several previous attempts at estimating the galactic confusion noise [17, 18, 10, 14, 11]. To characterize the confusion noise one must first determine how many galactic binaries are resolvable, but in order to figure out which binaries are resolvable, one must already have an estimate for the noise. The ideal solution is to perform a global fit, e.g. a full Bayesian analysis that fits both resolvable sources and noise at the same time as done by Littenberg [19]. Unfortunately, this procedure is extremely computationally intensive, and more efficient techniques are needed if we want to consider a range of population models and detector configurations for design studies. To this end, Timpano et al [10] developed an iterative subtraction scheme which starts with a simulated data set comprised of an instrument noise realization and the superposition of all gravitational waves produced by synthetic population of galactic binaries. The signal-to-noise ratio (SNR) of the GBs is calculated, and those above a specified threshold SNR are subtracted perfectly i.e. the true waveform is removed from the data stream using the simulated signal parameters. The noise estimate is updated after the bright signals are removed, and the SNRs of the remaining sources are re-computed. Those above the detection threshold are removed, and the whole process is repeated. It typically takes 5-6 iterations for the solution to converge.

It is this assumption of perfect signal recovery we wish to address in this paper. In reality the instrument plus confusion noise realization will randomly perturb the estimated parameters for the resolvable systems, resulting in an inaccurate signal recovery. Here we use the Maximum-Likelihood approximation and Fisher information matrix to estimate the parameter errors and waveform residuals.

The outline of this paper is as follows: In Section 2 we briefly review the galactic population model used to produce realizations of the LISA data used in the analysis. Next, in Section 3 we provide a review of relevant Maximum-Likelihood methods and how they can be used to estimation the noise-induced errors in signal recovery and parameter estimation of the resolved GBs. In Section 4 we extend the usual ML analysis to include noise spectral estimation, and in Section 5 we illustrate the relevance of the ML to a full Bayesian analysis using a simple model of a sinusoid in stationary Gaussian noise. In Section 6 we apply the Maximum-Likelihood approach to the global fitting of multiple signals and drive expressions for how the interaction between the signals impacts waveform and parameter estimation errors. We conclude in Section 7 by computing an improved estimate for the galactic confusion noise that takes into account parameter estimation errors in the bright source removal.

2 Instrument and Galactic Population Models

Our galaxy simulations use realizations of the the Toonen et al [20] population model provided by Valeriya Korol and Gijs Nelemans . The space density of interacting white dwarf binaries is reduced by a factor of ten relative to earlier models in response to the findings of recent observational studies[15, 16]. The population has $\sim 26$ million systems with gravitational wave frequencies above 0.1 mHz. The signals from these systems are simulated using an improved version of the fast waveform generation algorithm of Cornish and Littenberg [21]. The improved algorithm removes the need to sum over terms in the Fourier convolution by referencing the carrier frequency $f_{0}$ to the nearest integer multiple of the sample frequency, such that $f_{0}=m/T_{\rm obs}+\delta f$, and absorbing the factor of $e^{2\pi i\delta ft}$ into the slowly varying part of the signal. This removes the need for the sum in equation (A24) of Ref. [21], and significantly speeds up the waveform generation. In our analysis we use the full set of first-generation time-delay interferometry (TDI) [22] variables $\tilde{X}(f),\tilde{Y}(f),\tilde{Z}(f)$ given in Ref. [21], but when displaying results we show the more familiar Michelson-equivalent signals. To obtain the Michelson equivalent sensitivity we make use of the relation $S_{X}=4\sin^{2}\left(f/f_{*}\right)S_{M}$ where $S_{X}$ is the noise as seen in the TDI $X$ data channel and $S_{M}$ is the equivalent Michelson noise and $f_{*}=c/(2\pi L)$ is the transfer frequency. Our instrument noise model assumes white position noise $S_{p}(f)$ and colored acceleration noise $S_{a}(f)$ with spectral densities

[TABLE]

Under the assumption that the noise levels are the same in each link, we can form the noise-orthogonal $\{A,E,T\}$ channels [23]

[TABLE]

Below the transfer frequency $f_{*}$, where most signals are found, the $T$ channel is far less sensitive to gravitational waves, and does not contribute to our analysis.

Figure 1 shows a realization of instrument noise (assumed to be Gaussian and stationary), combined with the signal from all of the relevant GBs in the population synthesis realization plotted as the Michelson-equivalent sensitivity. We see that the galactic foreground rises above the instrument noise across the frequency range $0.3-20$ mHz.

3 Parameter estimation and waveform errors

In the high SNR regime, the likelihood is strongly peaked about the true model parameters, which allows for a Maximum-Likelihood analysis. Many useful results can be derived from a Taylor expansion of the likelihood about the true parameters. Here we provide a brief review of the Maximum-Likelihood (ML) approximation, deriving results for the parameter estimation errors and waveform residuals. We follow with a discussion of how the ML analysis relates to a global Bayesian analysis.

3.1 Maximum-Likelihood review

Consider the simple case of data s comprised of a gravitational wave signal $\textbf{h}_{T}=\textbf{h}(\vec{\lambda}_{T})$ and stationary, Gaussian noise n. The likelihood of observing s given the presence of a gravitational wave signal $\textbf{h}_{T}$ is then

[TABLE]

where

[TABLE]

defines the noise-weighted inner product taken across all independent data channels, and the one-sided noise spectral density in channel $I$ is given by the expectation value

[TABLE]

The noise $S_{n,I}(f)$ will include instrument noise and unresolved gravitational wave signals. The normalization factor $G$ is given by

[TABLE]

The traditional derivation of the maximum likelihood solution assumes that the noise model $S_{n,I}(f)$ is known, and that the signal model, $\textbf{h}(\vec{\lambda})$ is close to the true signal $\textbf{h}_{T}=\textbf{h}(\vec{\lambda}_{T})$. The signal model is then Taylor expanded about the true parameters:

[TABLE]

where $\Delta\vec{\lambda}=\vec{\lambda}-\vec{\lambda}_{T}$. The chi-squared in the likelihood can then be expanded as

[TABLE]

The maximum likelihood solution is found by setting $\partial_{i}\chi^{2}=0$, which yields

[TABLE]

where

[TABLE]

is the Fisher information matrix. Using the identity $\mathbbm{E}[(\textbf{n}|\textbf{g})(\textbf{n}|\textbf{k})]=(\textbf{g}|\textbf{k})$ we find that the error covariance matrix is given to leading order in the signal-to-noise ratio by the inverse of the Fisher information matrix:

[TABLE]

where the SNR $\rho$ is given by $\rho^{2}(\textbf{h})=(\textbf{h}|\textbf{h})$. See Vallisneri [24] for a more in depth presentation that discusses some of the potential pitfalls in using the Fisher Information matrix approximation parameter error estimation. It is important to note that there are higher order corrections to the signal parameters and covariance matrix which appear in Cutler and Flanagan [25] as equations (A31) and (A35).

3.2 Signal residuals

We can use the maximum likelihood approximation to study noise induced errors in the parameter recovery and signal subtraction for galactic binaries. A closely related analysis was performed by Cutler and Harms [12] in the context of subtracting the signals from neutron stars to allow for the detection of a primordial stochastic background for the Big Bang Observer mission concept. We extend their analysis to include noise modeling, and derive new expressions for the impact the foreground removal has on parameter estimation for other sources such as massive black hole mergers and extreme mass ratio insprials (EMRIs).

The noise-induced parameter estimation errors (9) result in waveform errors

[TABLE]

An example of the observed signal from a galactic binary and the noise-induced subtraction residual is shown in Figure 2.

Note that the residual is below the reference noise level as the waveform error is down-weighted by the level of overlap between noise and parameter derivatives of the signal.

The waveform error has zero mean, $\mathbbm{E}[\Delta\textbf{h}]=0$, and variance

[TABLE]

In the final step we have assumed that the error covariance matrix is approximated by the inverse of the Fisher matrix. The SNR of the residual depends only on the parameter dimension in the signal model and not upon the strength of the signal. Each term in the sum for $\Delta h$ is random walk induced by the noise realization, i.e. $\Delta h\sim n\sqrt{D}$ as there are $D$ terms. This means that $|\Delta\tilde{h}|^{2}\sim(n^{*}n)D$ and the inner product is weighted by the RMS noise resulting in a dependence only on model dimension. It can be shown that the variance of $\rho_{\Delta h}^{2}$ is $2D$ and the skew is $1/D^{2}$. The expectation value of the chi-squared can be written as

[TABLE]

where the last step follows from

[TABLE]

and $N$ is the number of data samples. We see that the signal residuals are anti-correlated with the noise, which results in a reduction in the chi-squared. Part of the noise gets absorbed by the signal model, which will ultimately result in a lowering of the confusion noise estimate relative to that found assuming perfect signal subtraction. Note that power spectrum of the residual $\textbf{s}-\textbf{h}$ has expectation value

[TABLE]

where

[TABLE]

is the power spectral density of $\Delta\textbf{h}$. Note that we made use of $\mathbbm{E}\left[\tilde{\textbf{n}}(f)(\textbf{n}|\partial_{i}\textbf{h}_{T})\right]=\partial_{i}\tilde{\textbf{h}}_{T}(f)$.

Our waveform model for a galactic binary has $D=9$ parameters: the sky location $(\theta,\phi)$; inclination and polarization angles $(\iota,\psi)$; amplitude $A$; reference phase $\phi_{0}$ and reference gravitational wave frequency $f_{0}$; and first and second frequency derivatives $\dot{f}_{0}$ and $\ddot{f}_{0}$. For most galactic binaries the frequency derivatives are poorly constrained and the effective model dimension is closer to $D=7$. The relevant quantity for estimating which systems have detectable frequency evolution are the number of frequency bins of evolution, $\alpha=\dot{f}_{0}T_{\rm obs}^{2}$ and $\beta=\ddot{f}_{0}T_{\rm obs}^{3}$, and the SNR. Roughly speaking, frequency evolution through $\sim(7/{\rm SNR})$ bins is detectable [26], and to leading post-Newtonian order we have

[TABLE]

The chirp mass ${\cal M}$ has been scaled to the mode of the population distribution [27]. From these expressions we see that only the loudest, most massive and highest-frequency systems will have a measurable second frequency derivative, and that most systems below 3 mHz will show no measurable frequency evolution at all. Including poorly constrained parameters in the model can lead to ill-conditioned Fisher matrices with badly behaved inverses. One solution is to reduce the model dimension by eliminating parameter $\lambda_{i}$ whenever the inner product $(\partial_{i}\textbf{h}_{T}|\partial_{i}\textbf{h}_{T})$ drops below some threshold. An alternative solution is to replace the Fisher matrix $\boldsymbol{\Gamma}$ with a matrix formed from the augmented Fisher matrix, ${\bf K}$, which includes derivatives of the priors (see Section 5 for details). We adopt the latter approach and include Gaussian priors on $\alpha,\beta$ centered on zero with width $\sigma=10$. The priors condition the matrix and when the parameters are poorly constrained by the data, have the effect of reducing the model dimension. Using the more stable approximation $\mathbbm{E}[\Delta\lambda^{i}\Delta\lambda^{j}]\approx\left(K^{-1}\right)^{ij}$, yields $\mathbbm{E}[\rho_{\Delta h}^{2}]\approx D_{\rm eff}$, where $D_{\rm eff}$ is the effective dimension of the model, defined by the number of parameters that have posterior distributions that differ measurably from their priors (a notion that can be made precise using the Kullback-Leibler divergence).

Figure 3 shows a histogram of square SNRs of the waveform residuals for the galactic binaries that are deemed detectable by the iterative subtraction scheme discussed in Section 7. The average value for there residual ${\rm SNR}^{2}$ of 7.6 is less than the full model dimension $D=9$, and consistent with our estimate for the effective dimension.

4 Maximum-Likelihood approximation with noise estimation

The standard treatment of the maximum likelihood expansion assumes that the noise spectrum is known. If the detectable gravitational wave signals are infrequent and short-lived, as is currently the case for compact binary mergers in LIGO, then “off-source” data from times where no detectable signals are present can be used to estimate the noise spectrum. These estimates will include instrument noise and unresolved gravitational wave signals. The option of making off-source estimates will not be available for LISA, and the noise spectrum will have to be inferred along with the signal model. Our derivation we will assume that we have a parameterized model for the noise, such as the cubic spline model used by the BayesLine algorithm [28].

To get an understanding for how noise modeling impacts the maximum likelihood calculation, consider a simple example with zero mean, additive, white Gaussian noise and $N$ data samples with likelihood

[TABLE]

The un-perturbed ($\textbf{h}=0$) noise level is given by sample variance

[TABLE]

We could expand $\sigma^{2}$ about the theoretical variance, $\sigma_{*}^{2}$, but it is simpler to expand $\sigma^{2}$ about the sample variance: $\sigma^{2}=\sigma_{0}^{2}+\Delta\kappa$ so that $\Delta\kappa=0$ when $\textbf{h}=0$. The typical difference between the sample variance and the theoretical variance will be by an amount that scales as the standard deviation of the sample variance, $\Delta\sigma_{0}^{2}=\sqrt{2}\sigma_{*}^{2}/\sqrt{N}$. The log likelihood can be expanded:

[TABLE]

where the notation $(a|b)_{0}$ denotes that the noise weighted inner product is taken with respect to the un-perturbed noise level $\sigma_{0}^{2}$. Setting $\partial_{\Delta\kappa}\log p(\textbf{s}|\textbf{h})=0$ and $\partial_{\Delta\lambda^{k}}\log p(\textbf{s}|\textbf{h})=0$, yields the ML solution

[TABLE]

We see that the leading order ML solution for the signal parameters is unchanged from the fixed noise case. The updated noise estimate $\sigma^{2}_{\rm ML}=\sigma^{2}_{0}+\Delta\kappa$ is lowered relative to its true value, as can be seen by taking the expectation value

[TABLE]

While the ML waveform removes some of the noise, this is now accounted for in the ML estimate for the noise, such that the expected value of the chi-squared is again just the dimension of the data: $\mathbbm{E}[\chi^{2}]=N$.

From expanding the likelihood around the signal and noise parameters $\vec{\eta}=\{\vec{\lambda},\kappa\}$ and maximizing the likelihood obtains the form

[TABLE]

where $\Gamma_{\mu\nu}=-\partial_{\mu}\partial_{\nu}\log p(\textbf{s}|\textbf{h})|_{\scriptsize\mbox{max}}$. Note that we are using Greek indicies to denote the entire collection of parameters. One can read off the Fisher matrix from the log-likelihood:

[TABLE]

where $\Gamma_{0,ij}$ is the Fisher matrix obtained from the signal-only ML analysis discussed in the previous section. The full Fisher matrix can be inverted by recognizing that the off-diagonal terms $(\textbf{n}|\partial_{j}\textbf{h}_{T})_{0}/\sigma_{0}^{2}$ are small compared to the block diagonal terms. We find that the variances for the signal parameters are inflated:

[TABLE]

On average $\left(\Gamma^{-1}\right)^{ij}\rightarrow\left(\Gamma^{-1}_{0}\right)^{ij}\left(1+2/N\right)$ where in the limit of large $N$ we obtain the original Fisher matrix. The signal model parameter variances are inflated by covariances with the noise model parameters as they both attempt to capture pieces of the signal. Note that covariance of the parameter shifts $\Delta\lambda^{j}$ and $\Delta\kappa$ from their true values, as computed in (23), does not equal the inverse of the Fisher matrix, $\mathbbm{E}\left[\Delta\eta^{\nu}\Delta\eta^{\mu}\right]\neq(\Gamma^{-1})^{\mu\nu}$. This is because the noise modeling changes the shape of the likelihood, and not just the location of the peak. However, we see from (25) that $\Gamma^{-1}$ does indeed describe the parameter covariances.

The ML expansion for a colored noise model is considerably more involved, and we relegate the details to Appendix A. To keep the notation simple we suppress the sum over data channels. Introducing the parameterized noise model

[TABLE]

(noise model derivatives are evaluated at the ML values after differentiation) where $S_{n,0}(f)$ is some smooth estimate of the instrument noise and unresolved signals (assuming perfect subtraction of the resolvable signals), we find the leading-order solution for the signal parameters has the same form as in (23), while the noise model parameters are given by

[TABLE]

The notation $(x|y)_{a}$ defines the inner product

[TABLE]

The integrals with the factor $(TS_{n,0}-2\tilde{n}^{*}\tilde{n})/S_{n,0}$ accounts for the difference between the theoretical noise model and fluctuation seen in a particular noise realization, i.e. the difference between $\sigma_{*}^{2}$ and $\sigma_{0}^{2}$ in the white noise toy model, as evidenced by its expectation value vanishing. Neglecting this difference and considering the expectation value of $\Delta\theta^{j}$ we obtain the simplification

[TABLE]

The white-noise case (23) can be recovered by setting $S_{n,0}=\frac{2}{T}\sigma_{0}^{2}$ and $\partial_{a}S_{n}(f)=\frac{2}{T}$, so that $\left[\int T\frac{S_{n,i}S_{n,j}}{S_{n,0}^{2}}df\right]^{-1}=\sigma_{0}^{4}/N$, and $(x|y)_{a}=(x|y)_{0}/\sigma_{0}^{2}$.

We can now compute the expectation value of the noise perturbation $\Delta S_{n}(f)=\Delta\theta^{a}\partial_{a}S_{n}(f)+\dots$

[TABLE]

Note that the perturbation to the noise model is negative, as it must be given that the signal model absorbs some of the noise. One would expect that $\Delta S_{n}(f)$ should be a smoothed representation of $-S_{\Delta h}(f)$ for an effective noise model, mopping up errors introduced by the signal ML.

Similar to the white noise model above we may obtain the signal model variances for a general noise model by making the same appeals to neglecting differences between the theoretical and sample variance and averaging over many noise realizations

[TABLE]

An effective noise model would minimize the factors $(\partial_{m}\textbf{h}_{T}|\partial_{n}\textbf{h}_{T})_{ab}$ such that $\left(\Gamma^{-1}\right)^{ii}\rightarrow\left(\Gamma_{0}^{-1}\right)^{ii}$. We can turn this into a more useful expression by taking advantage of the compact (in the frequency domain) nature of the GB signal and assume the noise PSD is roughly constant

[TABLE]

where $\Delta f$ is the bandwidth of the signal such that $T\Delta f$ is the number of frequency bins the GB spans. For sources that span many frequency bins such that the noise PSD cannot be assumed to be constant this serves as an upper limit for the increase in the parameter errors. Note that other terms exist when considering covariances of the signal model. Again, we see that when the source occupies a large bandwidth we recover the variances for when the noise is known. For a $3$ mHz source that experiences a $\Delta f\approx 0.6\times 10^{-6}$ Hz Doppler shift spreading due to LISA’s orbital motion, the parameter variances will be inflated by $10\%$ for a one year observation span, dropping to $3\%$ after four years.

5 Relating Bayesian Inference and Frequentist Maximum Likelihood Estimation

The LISA data will include overlapping signals from an unknown number of sources of different types. Bayesian inference provides a powerful and flexible framework for inferring the number and properties of the resolvable sources. In addition to the ontological differences between the Bayesian and Frequentist approach to statistical inference - Bayesian inference considers the data to be known and the signal parameters to be uncertain while Frequentist inference considers the signal parameters to be fixed and the data to be uncertain - Bayesian inference typically integrates over uncertainty (marginalization), while Frequentist analysis employs maximization. Despite these differences, the maximum likelihood analysis we have described can be used to estimate results from Bayesian inference by way of a Taylor expansion of the posterior distribution $p(\vec{\lambda}|\textbf{s})=p(\textbf{s}|\vec{\lambda})p(\vec{\lambda})/p(\textbf{s})$. Expanding about the mode of the posterior distribution (also known as the maximum a posteriori probability (MAP) estimate, $\vec{\lambda}_{\rm MAP}$), we have the quadratic approximation

[TABLE]

where $\Delta\vec{\lambda}=\vec{\lambda}-\vec{\lambda}_{\rm MAP}$ and

[TABLE]

When the likelihood is more informative than the prior, $\vec{\lambda}_{\rm MAP}\approx\vec{\lambda}_{\rm ML}$ and ${\bf K}$ is well approximated by the Fisher information matrix $\boldsymbol{\Gamma}$, though even small contribution from the derivatives of the log prior can have a an important stabilizing effect on ${\bf K}$.

To illustrate the relationship between the maximum likelihood analysis and Bayesian inference we produced simulated data consisting of stationary, Gaussian white noise with variance $\sigma^{2}$ and a sinusoidal signal $h(A,f_{0},t_{0},\phi_{0})=A\cos(2\pi f_{0}(t-t_{0})+\phi_{0})$. We held the phase parameter $\phi_{0}=\pi$ fixed in the analysis as otherwise there is a near perfect degeneracy between the time offset $t_{0}$ and the phase offset $\phi_{0}$, which significantly complicates the analysis. The simulated data consisted of $N=10^{4}$ evenly spaced samples spanning $T=100$ seconds, with $A=\sqrt{2}$, $f_{0}=0.25$ Hz, $t_{0}=1$ second and $\sigma^{2}=100$. The noise level was set to yield a signal-to-noise ratio of ${\rm SNR=10}$. A plot of the simulated data and signal are shown in the upper panel of Figure 4.

For the Bayesian analysis we assumed uniform priors on the signal and noise model parameters $(A,f_{0},t_{0},\sigma^{2})$ across a range that was much wider than the predicted statistical errors so that posterior distribution and the likelihood were effectively identical. A Markov Chain Monte Carlo (MCMC) simulation was used to compute the mean and variance of the signal parameters and the waveform error, while (23) and (25) were used to estimate the parameter shifts at maximum likelihood and the variances. The MCMC and ML derived values for the parameter shifts and standard deviations are listed in Table 1 for a particular noise realization. The marginalized posterior distributions for the parameters are compared to the predictions of the Gaussian approximation (35) in Figure 5. The agreement between the ML and MCMC seen in this example was typical of what we found when repeating the analysis for different noise realizations. That is not to say that the Gaussian approximation will be this accurate in more realistic settings where the noise is more complicated and the parameters are highly correlated [24], but it does provide useful order-of-magnitude estimates in most situations.

The displacement of the parameters from their true values is related to the waveform error $\Delta h$ shown in the lower panel of Figure 4 which displays the ML, MAP and mean waveform errors. Note that the mean waveform error is found by averaging the waveform errors, and not by using the mean parameter values to compute a waveform.

The good agreement between the maximum likelihood and mean waveform residual hides a key difference between the frequentist and Bayesian analyses: in the Bayesian global fit the waveform uncertainties are marginalized over, while the frequentist analysis uses point estimates. Rather than subtracting a particular point estimate for each signal from the data, the Bayesian approach subtracts a range of estimates for each signal such that the residual is consistent with the noise model. This procedure is illustrated in Figure 6 for the sinusoid signal model, where the waveform residuals from each iteration of the MCMC analysis are used to produce a histogram of the residual at each time sample. Also shown in Figure 6 is the MAP point estimate for the waveform residual. Notice that the full posterior distribution for the waveform residuals has significant spread about the point estimate.

6 Multiple Sources

The LISA data will contain many signals that partially overlap in both time and frequency. Extracting information about these signals necessitates finding a global solution for all signals that can be resolved - that is, signals that are both sufficiently loud and sufficiently distinct to be individually identified. The full LISA data stream can be written as

[TABLE]

where $\textbf{H}=\sum\textbf{h}_{T}$ denotes the sum total of all gravitational wave signals, which can be further separated into a resolved, $\textbf{H}_{R}$, and un-resolved $\textbf{H}_{U}$, component. The unresolved component is often referred to as “confusion noise”, the largest component of which is expected to come from white dwarf binaries in our galaxy. The ML analysis for single sources can be applied to multiple sources by replacing h with $\textbf{H}_{R}$, and by replacing n by $\textbf{n}+\textbf{H}_{U}$. The resolved signals will include a large number of bright galactic binaries, along with multiple supermassive black holes and EMRIs.

The parameter estimation for the resolved systems will be impacted by the unresolved signals, which add to the effective noise level, and by the other resolved signals due to signal overlap. To simplify the discussion imagine that the resolved signals consist of bright galactic binaries $\textbf{H}_{G}$ and a single massive black hole binary $\bar{\textbf{h}}$. The parameter vector $\vec{\lambda}$ runs over the full set of galactic binary parameters (denoted by indices early in the alphabet, $\lambda^{a},\lambda^{b}\dots$ etc) and the black hole parameters (denoted by indices later in the alphabet $\lambda^{i},\lambda^{j}\dots$ etc). The full set of signal parameters for the galactic sources and the black hole are indicated by Greek indices. The Fisher information matrix for the combined solution, $\Gamma_{\alpha\beta}$ can be broken into a block diagonal part $B_{\alpha\beta}$ formed from a galactic-binary block, $G_{ab}=(\partial_{a}\textbf{H}_{G}|\partial_{b}\textbf{H}_{G})$, and a black hole block $\bar{B}_{ij}=(\partial_{i}\bar{\textbf{h}}|\partial_{j}\bar{\textbf{h}})$, and a mixed block $M_{aj}=(\partial_{a}\textbf{H}_{G}|\partial_{j}\bar{\textbf{h}})$. The waveform error for the black hole signal is then

[TABLE]

The expectation value for the squared SNR of the black hole waveform residual is then

[TABLE]

Assuming the mixture terms $M_{aj}$ are small compared to the terms on the diagonal, the inverse of the full Fisher matrix can be expanded as

[TABLE]

The black-hole block of the inverse, $\left(\Gamma^{-1}\right)^{ij}$, lacks the second term since $\left(B^{-1}\right)^{i\mu}M_{\mu\nu}\left(B^{-1}\right)^{\nu j}=\left(\bar{B}^{-1}\right)^{ik}M_{kl}\left(\bar{B}^{-1}\right)^{lj}=0$. Therefore we have

[TABLE]

where $\bar{D}$ is the dimension of the black hole model. We see that the black hole waveform residuals are inflated from the isolated source result by the mixture term

[TABLE]

The second line is obtained using Isserilis’s theorem and dropping the cross terms $\mathbbm{E}[\Delta\lambda^{i}\Delta\lambda^{a}]$ which would produce higher order corrections. Using $\mathbbm{E}[\Delta\tilde{H}_{G}(f)\Delta\tilde{H}^{*}_{G}(f^{\prime})]=\frac{1}{2}S_{\Delta H_{G}}(f)\delta(f-f^{\prime})$ we have

[TABLE]

If we were to switch the roles of the black hole and the resolved galactic binaries we would find the the squared SNR of the galactic residuals were inflated by exactly the same amount:

[TABLE]

On the other hand, the squared SNR of the full residual is equal to the total parameter dimension:

[TABLE]

These results can be reconciled by noting that

[TABLE]

and using

[TABLE]

Thus we see that the extra residual for each source class is canceled by the cross-correlation of the residuals between the source classes. Note that the results for the SNR of the signal residuals are unchanged to the order we are considering when using the full noise model or the unperturbed noise model.

Next we consider the impact on the black hole parameter estimation errors caused by fitting the bright galactic binaries. The variance in the parameter estimation errors can be estimated from the diagonal entries of inverse of the full Fisher information matrix

[TABLE]

The second term in the expansion comes from correlations between the black hole signal and the resolved galactic binaries as is positive definite since $\textbf{x}^{T}\mathbbm{G}\textbf{x}\geq 0$ for a positive-definite matrix. Thus, the simultaneous fitting of the galactic binary signals and the black hole signal tends to inflate the parameter estimation errors. Expanding to leading order the second term is given by

[TABLE]

where $\partial_{a}H_{R}=\partial_{a}H_{G}$ as the derivatives are with respect to GB parameters. Using an estimate to the waveform errors of the resolved sources we can express this result in a more useful form:

[TABLE]

This implies that the covariance matrix inflates with the following upper bound:

[TABLE]

where $dN/db$ of sources resolved per frequency bin. In the next section we will obtain an estimate for the noise due to GB waveform errors which will allow us to obtain a more useful expression for this overlap term. A similar inflation of GB parameter variances results from the overlap with the BH signal.

7 Galactic Confusion Noise

To obtain an estimate of the Galactic confusion foreground we employ an iterative subtraction scheme. Previously, this scheme was performed with perfect removal of source waveform [10, 11], but clearly noise will lead to errors in the parameters estimation and signal subtraction, and a reduction of the confusion noise estimate.

The revised procedure is as follows: A simulated data set is produced that includes a realization of the instrument noise and the sum of the strains due to the galactic binaries $H$ from the galactic population model. A smooth fit to the power spectral density of the instrument noise and the signals is used as an initial estimate for the noise in each data channel. Next we identify sources which are loud (SNR $>7$) relative to this noise estimate and subtract the best-fit waveform $\textbf{h}(\vec{\lambda})_{{\rm bestfit}}=\textbf{h}_{T}-\partial_{i}\textbf{h}_{T}\Delta\lambda^{i}$ from the data. A smooth fit to the power spectral density of the remaining signals and noise is computed, and signals above the SNR threshold for the updated noise estimate are identified and subtracted. As can be seen in Figure 7, the subtraction procedure quickly converges.

The number of sources which can be resolved converges after just 5 or 6 iterations. The number density of sources, measured in terms of the number per a frequency bin, $dN/db$, is shown in Figure 8 for a single channel and dual channel analysis. We see that more sources can be resolved when multiple data channels are used in the analysis. Attempts have been made to deal with the identification and subtraction of signals which overlap [29, 30] and of how many sources per frequency band [31] can be resolved. Here we see that with two data channels and a 4-year mission the peak density is roughly one source resolved per ten frequency bins. In our simulations we made the simplifying assumption that the augment Fisher matrix for the galactic population is block diagonal. That is, we ignored correlations between galactic signals. This approximation is reasonable when $dN/db$ is small, but may be questionable in the highest density regions and for the occasional systems that happen to have high overlap. We will re-visit this complication in a future study, as the parameter estimation errors grow significantly for highly overlapping systems [13].

Figure 9 compares the Michelson-equivalent strain power spectral densities for the imperfect and perfect subtraction scheme. The dashed lines show galactic confusion noise and the solid lines show the combined instrument and confusion noise.

Note that the differences between the PSD arise where the most sources are resolved as one would expect (see Figure 8). In the dual $A,E$ channel the PSD is lower as indicated by the noise levels specified by the reference frequency $2$ mHz.

We can estimate the power spectral density of the combined waveform residual, $S_{\Delta H_{R}}(f)=\frac{T}{2}\mathbbm{E}\left[|\Delta\tilde{H}_{R}(f)|^{2}\right]$ by applying (13) to the full compliment of $N$ resolved binaries:

[TABLE]

Considering the contribution in a small frequency range $\Delta f$ centered at $f$ we find

[TABLE]

Figure 10 compares this analytic estimate to the numerical value found from the iterative subtraction scheme. The prediction lines up quite well with a small deviation at the frequencies where the number of resolved sources per bin peaks.

8 Discussion

We have used the maximum likelihood approximation to derive a number of analytic results pertaining to the LISA global analysis problem. A simple toy model was used to demonstrate the relevance of these estimates to a full Bayesian analysis, though we cannot guarantee that the approximations will be as reliable when applied to LISA data analysis. We extended the standard maximum likelihood analysis to include noise modeling, and found that the estimated noise level is lowered whenever signals are subtracted from the data. We applied our general results to the simultaneous fitting of a black hole binary and a collection of galactic binaries and found that the errors in the black hole waveform recovery are increased, as wells as the variances of the source parameters. It is important to note that it is not the overlap of the signals which cause the inflation of variances, but rather, it is due to the waveform subtraction errors. We concluded by incorporating parameter estimation errors in the estimation of the galactic confusion noise, and derived a useful expression that can be used to predict the reduction in the confusion noise in terms of the number density of resolved signals. Equations (34) and (52) provide quick estimates for how parameter estimation errors are inflated by noise fitting and source confusion.

Acknowledgments

We are grateful for the support provided by NASA grant NNX16AB98G.

Appendix A

Here we derive the more general signal plus noise model ML analysis results presented in section 4. Beginning with the log-likelihood obtained from equation (3)

[TABLE]

where the noise-weighted inner product is now parameterized by $\vec{\theta}$ as denoted by $(\cdot|\cdot)_{\vec{\theta}}$. We may expand about the true noise model $S_{n,0}(f)$ and the true signal model $\textbf{h}_{T}$ and maximized to obtain estimates of $\Delta\vec{\theta}$ and $\Delta\vec{\lambda}$. Expanding out the normalization constant and dropping terms which are constant with respect the maximization gives

[TABLE]

where $T$ is the observation period. An arbitrary noise-weighted inner product expanded out takes the form

[TABLE]

where

[TABLE]

With these pieces in hand the chi-squared piece of the log-likelihood, dropping constants with respect to the maximization, is

[TABLE]

Collecting terms, and maximizing results in the solution

[TABLE]

where $\Delta\textbf{h}$ need only be kept to leading order. The Fisher matrix can be obtained similarly to the toy white noise problem presented in section 4

[TABLE]

Inverting the Fisher matrix and considering the signal model variances

[TABLE]

9 References