Calibration Uncertainty’s Impact on Gravitational-Wave Observations

Modern interferometric Gravitational-Wave (GW) detectors [1, 2, 3] provide unprecedented access to fundamental physics on both the largest and smallest scales. Advanced detector technologies probe aspects of quantum interactions with macroscopic (classical) objects [4, 5], and the astrophysical signals detected probe gravity in the strong/dynamical-field regime [6, 7], the distribution of rare stellar remnants, including black holes (BHs) and neutron stars (NSs) [8, 9], the behavior of matter at densities up to several times that of stable nuclei (e.g., Refs. [10, 11]), and many other extreme environments that are otherwise difficult or impossible to access. Much of this success can be attributed to the detectors’ ability to faithfully report the astrophysical strain so that meaningful comparisons to theoretical expectations can be made. The process of predicting how astrophysical signals will appear within the detector is referred to as calibration. It typically proceeds by measuring the opto-mechanical response of the interferometer’s cavities and control loops.

Calibration has long been of general interest in the field. Although many calibration techniques exist, the basic physical measurement is the change in the detector output as a function of the change in the difference between the lengths of the interferometer’s arms. In current interferometric observatories (IFOs), one typically moves the test masses (mirrors) that bracket the Fabry-Pérot cavity in each arm by a known amount at a known frequency and then measures the detector response. Initial techniques involved driving the suspension system, composed of a cascade of pendula, that suspend each test mass and isolate it from many sources of terrestrial noise (see, e.g., Refs. [12, 13] and references therein). Recently, several authors have investigated Newtonian calibrators [14, 15], which drive the test masses through the Newtonian gravitational force exerted by an oscillating mass multipoles (rotating dumbbell). Additionally, several authors have studied the ability of astrophysical signals themselves to calibrate detectors. In this approach, one assumes that the correct theory of gravity is known and compares the detector output to the expected change in length induced by different astrophysical sources. Indeed, there is sustained interest in understanding both how to use astrophysical signals to calibrate our detectors and how uncertainty in calibration can affect our inference of the properties astrophysical signals [16, 17, 18, 19, 20].

Currently, the most successful calibration techniques rely on auxiliary laser systems (photon calibrators or PCALs [21, 22]) to exert a known radiation force (laser with known intensity) on the test masses. PCALs have met with great success, providing precise calibration measurements throughout the advanced detector era [23]. However, through the careful quantification of the detector response enabled by PCAL measurements, it is now understood that systematic modeling uncertainties in the detector response, including variations over time, are typically larger than the uncertainty in the directly-measured parameters of the model for the detector response [23, 13]. These systematic uncertainties are represented with a Gaussian process (GP) and uncertainty envelopes derived therefrom are available at a high cadence [24].

It has become common practice to marginalize over the calibration uncertainty at the time of individual events within inferences of those events’ astrophysical parameters. Several techniques exist, and they differ primarily in how they represent the calibration uncertainty. Traditionally, a cubic spline with knots spaced logarithmically in frequency is used [25]. The marginal prior distributions for the calibration error at each knot is chosen to match the measured calibration envelope at that frequency. However, it is not necessarily the case that this spline-prior is a good representation of the full GP, which includes correlations between frequencies and may not be a slowly varying function of frequency. Several authors have also investigated physical models of calibration uncertainty, even parametrizing the model of the detector response in addition to sampling from the full GP prior for the model systematics [19, 20]. Marginalization over calibration is then performed numerically. While this approach is computationally tractable, it does not offer an immediate interpretation of which aspects of the calibration uncertainty most affect the astrophysical inference or, in turn, which astrophysical signals could most improve our knowledge of the detectors’ calibration.

I show how to analyze GP calibration priors analytically, as the GP prior is conjugate to the GW likelihood under the common assumption of stationary Gaussian noise. This allows me to derive several expressions for the impact of calibration uncertainty on searches for, and the inference of, astrophysical signals. I also examine how astrophysical signals can inform our understanding of calibration in our detectors.

While these expressions should be of general interest, current calibration uncertainties are typically small ($O(\mathrm{few}\%)$ in amplitude and phase [23, 13]). That is, signals detected with current IFOs are quiet enough that the calibration uncertainty does not matter much [26, 27, 28, 29, 30]. This may not be the case for proposed third-generation (3G) detectors, in which the typical signal-to-noise ratio (S/N) may be larger. I clarify the S/N thresholds above which uncertainty in detector calibration will dominate over uncertainty from stationary Gaussian noise. What’s more, opto-mechanical effects in the detector response for 3G detectors will be much more complicated than in current detectors, as light’s round-trip travel time within the Fabry-Pérot cavities becomes comparable to the frequencies probed in the astrophysical strain [31]. Altogether, then, astrophysical inference of the loudest events detected with 3G detectors will almost certainly be limited by the understanding of the detector responses, and several common approximations for current detector responses will need to be revisited. This fact appears to not have been fully investigated within the 3G literature (see, e.g., discussion in Ref. [32]), and I hope that the analytic expressions presented in this manuscript will enable improved estimates of the scientific capabilities of 3G instruments.

I first summarize a few basic assumptions common in GW data analysis and introduce my notation in Sec. II. I then make use of the GP representation of the calibration uncertainty in Sec. III to obtain analytic expressions for the conditioned distribution for the calibration error given the observed data and an astrophysical signal as well as the marginal distribution for the observed data given an astrophysical signal integrated over uncertainty in detector calibration. Several limiting cases are explored in Sec. III.1, and Sec. III.2 explores extensions to networks of detectors and the joint inference over multiple signals. I discuss the possible impact on searches in Sec. IV. Sec. V discusses the measurability of astrophysical parameters based on the calibration-marginalized likelihood, with particular focus on how much information about the astrophysical parameters is lost due to uncertainty in the detector response. The similarity of calibration errors to waveform uncertainty and deviations from General Relativity (GR) is introduced within Sec. VI. I conclude with a discussion of calibration requirements in Sec. VII.

The likelihood of observing data $d$ in stationary Gaussian noise described by a one-sided power spectral density (PSD) $s$ within a single interferometer in the presence of a signal described by intrinsic parameters $\vartheta$ and extrinsic parameters $\varphi$ is typically written as

$$\ln p(d|\theta,\varphi,\mathcal{H}_{s})\supset-\frac{1}{2}\left(4\int\limits_{0}^{\infty}\mathrm{d}f\,\frac{|d-h(\vartheta,\varphi)|^{2}}{s}\right)$$

(1)

where $h=f_{\alpha}(\varphi)h_{\alpha}(\vartheta)$ is the strain projected in the interferometer, $f_{\alpha}$ is the detector response to the $\alpha$ polarization, $h_{\alpha}$ is the waveform in the $\alpha$ polarization, and repeated indices imply summation. In general, the detector response $f_{\alpha}$ is a function of the frequency [31]. The PSD is defined in terms of the noise autocorrelation function as

$$\left<n(t)n(t+\tau)\right>=\frac{1}{2}\int\mathrm{d}f\,e^{2\pi if\tau}s(f)$$

(2)

where $i=\sqrt{-1}$ and $\left<\cdot\right>$ denotes an average over noise realizations. $\mathcal{H}_{s}$ denotes the assumption of stationary Gaussian noise with a known PSD. See, e.g. Refs. [33, 34] for more discussion about PSD estimation. Under this assumption, the noise model is a stationary Gaussian process (GP) with (auto)correlations defined by Eq. 2.

However, as IFOs record only discrete data at regularly spaced intervals, it will be convenient to work explicitly with matrices in the following. In most cases, I suppress indices for simplicity. However, for notational clarity, I denote objects that have a single frequency index with lower case letters and objects that have two frequency indices with upper case letters. It will also be convenient to work with only real variables. As such, I explicitly include the real and imaginary parts of data separately so that

	$\displaystyle h_{2i}$	$\displaystyle=\mathcal{R}\{h(f_{i})\}$		(3)
	$\displaystyle h_{2i+1}$	$\displaystyle=\mathcal{I}\{h(f_{i})\}$		(4)

for $i\in\{1,\cdots,N_{\mathrm{frq}}$}. Therefore, for complex data recorded at $N_{\mathrm{frq}}$ discrete frequencies, vectors have $2N_{\mathrm{frq}}$ real entries and matrices have $4N_{\mathrm{frq}}^{2}$ real elements.

With this notation, I obtain

	$\displaystyle\ln$	$\displaystyle p(d|\theta,\varphi,\mathcal{H}_{s})$
		$\displaystyle=-2\sum_{i}^{N_{\mathrm{frq}}}\left[\Delta f\,\frac{\mathcal{R}\{d(f_{i})-h(f_{i})\}^{2}+\mathcal{I}\{d(f_{i})-h(f_{i})\}^{2}}{s(f_{i})}\right]$
		$\displaystyle\quad\quad\quad\quad\quad\quad-N_{\mathrm{frq}}\ln 2\pi-\sum\limits_{i}^{N_{\mathrm{frq}}}\ln\left(\frac{s(f_{i})}{4\Delta f}\right)$
		$\displaystyle=-\frac{1}{2}(d-h)^{\mathrm{T}}C^{-1}(d-h)$
		$\displaystyle\quad\quad\quad\quad\quad\quad-N_{\mathrm{frq}}\ln 2\pi-\frac{1}{2}\ln\mathrm{det}|C|$		(5)

where $C$ is the covariance matrix (expected to be diagonal) with

$$C_{2i\,2i}=C_{2i+1\,2i+1}\equiv\left(\frac{s(f_{i})}{4\Delta f}\right)$$

(6)

for stationary Gaussian noise. That is, there are neither correlations between frequencies nor between the real and imaginary parts of the data at each frequency. In general, nonstationary noise and/or windowing functions applied within discrete Fourier transforms may introduce off-diagonal terms within $C$ (see, e.g., Ref. [35]). I ignore these for the purposes of this study, although I do not assume $C$ is diagonal in what follows unless explicitly stated.

Now, if the modeled detector response is incorrect by a multiplicative factor such that $f^{(\mathrm{true})}_{\alpha}=f_{\alpha}(1+\delta)$, the actual projected strain in the detector will be $h^{(\mathrm{true})}=h(1+\delta)$. This assumes that the calibration uncertainty affects only the total projected strain in the interferometer and therefore impacts the detector response to each polarization in the same proportion.¹¹1This may not be true in general, though, as approximations to the interferometer’s response may break down differently for the response to each polarization. These errors may also depend on the source’s extrinsic parameters (i.e., relative orientation of the source and the detector). However, these effects are small compared to the GP uncertainty for current detectors, and the approximation that $\delta$ is independent of polarization and $\varphi$ should be reasonable. In general, $\delta$ can be complex and mixes the real and imaginary parts of the projected strain so that

	$\displaystyle\mathcal{R}\{h^{(\mathrm{true})}\}$	$\displaystyle=\mathcal{R}\{h\}+\mathcal{R}\{h\}\mathcal{R}\{\delta\}-\mathcal{I}\{h\}\mathcal{I}\{\delta\}$		(7)
	$\displaystyle\mathcal{I}\{h^{(\mathrm{true})}\}$	$\displaystyle=\mathcal{I}\{h\}+\mathcal{I}\{h\}\mathcal{R}\{\delta\}+\mathcal{R}\{h\}\mathcal{I}\{\delta\}$		(8)

which I express as

$$h^{(\mathrm{true})}=h+H\delta$$

(9)

where $H$ is a ($2\times 2$ block)-diagonal matrix with each block given by

$$H_{2\times 2}=\begin{bmatrix}\mathcal{R}\{h(f_{i})\}&-\mathcal{I}\{h(f_{i})\}\\ +\mathcal{I}\{h(f_{i})\}&\mathcal{R}\{h(f_{i})\}\end{bmatrix}$$

(10)

I can then replace $h$ in Eq. 5 with Eq. 9 to obtain the likelihood of observing $d$ in the presence of a true signal $h$ and calibration error $\delta$.

It is common to parametrize $\delta$ as

$$1+\delta=(1+\delta A)e^{i\delta\psi}$$

(11)

where both $\delta A$ and $\delta\psi$ are real. I instead use $\mathcal{R}\{\delta\}$ and $\mathcal{I}\{\delta\}$ because this allows me to marginalize over them analytically. Nonetheless, note that $\delta\approx\delta A+i\delta\psi$ in the limit $\delta A$, $\delta\psi\ll 1$, which is often the case. In this limit, then, uncertainty in the amplitude and phase, including correlations between the two, can be immediately translated into uncertainty on the real and imaginary parts of the calibration error.

Often, one has some prior knowledge of $\delta$ (both amplitude and phase) as a function of frequency, which allows one to place a reasonable prior on $\delta$. This is informed by direct measurements of the detector response and models of the interferometer’s cavities [23, 13]. With such a prior, one can construct a joint distribution for both $d$ and $\delta$ assuming $\delta$ is independent of the noise realization and source properties:

$$p(d,\delta|\vartheta,\varphi,\mathcal{H}_{s},\mathcal{H}_{\delta})=p(d|\delta,\vartheta,\varphi,\mathcal{H}_{s})p(\delta|\mathcal{H}_{\delta})$$

(12)

in which $\mathcal{H}_{\delta}$ denotes the prior assumptions about the calibration uncertainty. This as the starting point for Sec. III.

It has become common practice to model the calibration uncertainties with a GP [23] in order to capture possible systematic errors between the modeled and measured detector responses. These systematic errors typically dominate over the uncertainty in the detector response from the measured parameters of the opto-mechanical model [13], and, as such, I consider a GP prior for $\delta$ with mean $\gamma$ and covariance $\Gamma$ such that

	$\displaystyle\ln p(\delta|\mathcal{H}_{\delta})$	$\displaystyle=-\frac{1}{2}(\delta-\gamma)^{\mathrm{T}}\Gamma^{-1}(\delta-\gamma)$
		$\displaystyle\quad\quad-N_{\mathrm{frq}}\ln 2\pi-\frac{1}{2}\ln\mathrm{det}|\Gamma|$		(13)

This choice is particularly convenient because it allows me to analytically marginalize over calibration uncertainties in Eq. 12. That is, a GP prior for $\mathcal{R}\{\delta\}$ and $\mathcal{I}\{\delta\}$ is conjugate to the likelihood for stationary Gaussian noise.²²2One may be able to approximate more general priors for $\delta$ as sums of GPs. To wit, marginalizing over $\delta$ yields

	$\displaystyle\ln p(d|\vartheta,\varphi,\mathcal{H}_{s},\mathcal{H}_{\delta})=$	$\displaystyle\ln p(d|\vartheta,\varphi,\mathcal{H}_{s})$
		$\displaystyle-\frac{1}{2}\gamma^{\mathrm{T}}\Gamma^{-1}\gamma-\frac{1}{2}\ln\mathrm{det}|\Gamma|$
		$\displaystyle+\frac{1}{2}\nu^{\mathrm{T}}A^{-1}\nu+\frac{1}{2}\ln\mathrm{det}|A|$		(14)

where

$$\nu=A\left[H^{\mathrm{T}}C^{-1}(d-h)+\Gamma^{-1}\gamma\right]$$

(15)

is a vector that is linear in the data and

$\displaystyle A^{-1}=H^{\mathrm{T}}C^{-1}H+\Gamma^{-1}$

(16)

As will be seen in Sec. III.1, limiting cases of the prior knowledge about calibration errors determine the behavior of $A$. Indeed, a comparison of the two terms in $A$ is essentially a question of whether the typical fluctuation in projected strain from calibration errors ($\sim|h|^{2}\Gamma$) is large compared to the typical fluctuations from Gaussian noise ($\sim s/4\Delta f$).

Eq. III shows that the marginalization can be expressed as a multiplicative factor that modifies the likelihood obtained by assuming perfect calibration. In this way, one can quickly reweigh ($\vartheta,\varphi$)-samples drawn assuming perfect calibration under different assumptions about the calibration error post hoc. A similar strategy was investigated via numeric marginalization in Ref. [19].

Similarly, conditioning on $d$ yields

	$\displaystyle\ln p(\delta|d,\vartheta,\varphi,\mathcal{H}_{s},\mathcal{H}_{\delta})=$	$\displaystyle-\frac{1}{2}(\delta-\nu)^{\mathrm{T}}A^{-1}(\delta-\nu)$
		$\displaystyle\quad\quad-N_{\mathrm{frq}}\ln 2\pi-\frac{1}{2}\ln\mathrm{det}|A|$		(17)

One can approximate the marginal posterior on $\delta$ as a Monte Carlo sum of GPs such that

$$p(\delta|d,\mathcal{H}_{s},\mathcal{H}_{\delta})\approx\frac{1}{N_{\mathrm{smp}}}\sum\limits_{k}^{N_{\mathrm{smp}}}p(\delta|d,\vartheta^{(k)},\varphi^{(k)},\mathcal{H}_{s},\mathcal{H}_{\delta})$$

(18)

with $N_{\mathrm{smp}}$ $(\vartheta^{(k)},\varphi^{(k)})$-samples drawn from Eq. III. Alternatively, one could draw samples from Eq. 5 and then construct a sum with nontrivial weights using the correction from Eq. III.

At times, it will be more convenient to express Eq. III in terms of a single contraction involving $d$. This allows me to investigate the effective covariance matrix for $d$ induced by the marginalization over calibration uncertainty. I therefore write

	$\displaystyle\ln p(d|\vartheta,\varphi,\mathcal{H}_{s},\mathcal{H}_{\delta})$	$\displaystyle=-\frac{1}{2}(d-h-\mu)^{\mathrm{T}}B^{-1}(d-h-\mu)$
		$\displaystyle\quad\quad-N_{\mathrm{frq}}\ln 2\pi-\frac{1}{2}\ln\mathrm{det}|B|$		(19)

with effective (inverse) covariance matrix

$$B^{-1}=C^{-1}-C^{-1}HAH^{\mathrm{T}}C^{-1}$$

(21)

and

$$\mu=BC^{-1}HA\Gamma^{-1}\gamma$$

(22)

These are my basic results, and I explore their applications in what follows.

A search for a known signal morphology that assumes perfect calibration produces a maximum likelihood estimator for the signal amplitude (denoting the signal $h=a\tilde{h}$)

	$\displaystyle\hat{a}_{\sigma_{\Gamma}=0}$	$\displaystyle=\frac{d^{\mathrm{T}}C^{-1}\tilde{h}}{\tilde{h}^{\mathrm{T}}C^{-1}\tilde{h}}$
		$\displaystyle=\left(4\int df\frac{\mathcal{R}\{d^{\ast}\tilde{h}\}}{s}\right)\left(4\int df\,\frac{|\tilde{h}|^{2}}{s}\right)^{-1}$		(48)

which is an unbiased estimator with unit variance under the common convention $\tilde{h}^{\mathrm{T}}C^{-1}\tilde{h}=1$. When $\tilde{h}$ is normalized in this way, the estimator is often written as

$$\hat{a}_{\sigma_{\Gamma}=0}=\frac{d^{\mathrm{T}}C^{-1}\tilde{h}}{\sqrt{\tilde{h}^{\mathrm{T}}C^{-1}\tilde{h}}}$$

(49)

This is the normal definition of the matched-filter S/N. I examine the behavior of this statistic in the presence of calibration uncertainty.

Specifically, when $\gamma=0$,

	$\displaystyle E[\hat{a}_{\sigma_{\Gamma}=0}]$	$\displaystyle=a_{\mathrm{true}}$		(50)
	$\displaystyle\mathrm{Var}[\hat{a}_{\sigma_{\Gamma}=0}]$	$\displaystyle=E[(\hat{a}_{\sigma_{\Gamma}=0}-E[\hat{a}_{\sigma_{\Gamma}=0}])^{2}]$
		$\displaystyle=\frac{\tilde{h}^{\mathrm{T}}C^{-1}BC^{-1}\tilde{h}}{\tilde{h}^{\mathrm{T}}C^{-1}\tilde{h}}>1$		(51)

That is, the mean S/N is not affected for zero-mean calibration priors, but the variance is larger ($B\geq C$). That is, I assume $C$ is the covariance of $d$ in the absence of a signal. I assume this is known (measured) and is independnet of $\delta$ (Eq. 5).⁴⁴4Ref. [40], instead, suggests the measured $C$ will depend on $\delta$. This means one must include $C$’s dependence on $\delta$ when computing moments of $\hat{a}_{\sigma_{\Gamma}=0}$. They find a cancellation between the numerator and denominator of Eq. 49 at leading order in $\gamma$ when computing the expected value. In the context of their analysis, my result would hold if one performed the inference directly on (any deterministic function of) the raw detector output. Nonetheless, although not shown explicitly, their analysis suggests $\mathrm{Var}[\hat{a}_{\sigma_{\Gamma}=0}]\geq 1$ as well, which is my main result. This implies that, for a given observed S/N, the probability that noise alone (combined with calibration uncertainty) could have produced that statistic is higher. Therefore, the statistical significance of that S/N is lower. If one considers the the ratio of the estimator to its standard deviation as a measure of significance (a classical $z$-test), then $\hat{a}_{\mathrm{naive}}$ will produce strictly smaller $z$ than one would expect without calibration uncertainty. This makes sense intuitively; it should be harder to detect signal with a matched filter in the presence of calibration uncertainty.

Similar effects occur for unmodeled searches. The variance of the naive search statistic (maximized likelihood) is increased under the noise hypothesis in the presence of calibration uncertainty. Louder signals are needed to stand out against the background.

It is noteworthy that modern GW searches do not estimate the false alarm probability based on the assumption of stationary Gaussian noise. Instead, they employ a variety of bootstrap procedures (see, e.g., Refs. [41, 42, 43, 44, 45, 46]). De facto, then, there is some threshold above which the observed S/N must fall to be considered a confident detection. In this context, the broader variance of the naive S/N statistic in the presence of nonzero calibration uncertainty still hampers detectability. This is because the true S/N (sometimes called the optimal S/N, defined in Eq. 42) must be systematically larger than the threshold in order for the observed S/N to be above the threshold with high probability. That is, the optimal S/N must be shifted above the threshold by a larger amount when $\Gamma$ is nonzero because the variance of S/N is larger. Larger optimal S/N imply closer sources, and therefore the space-time volume to which detectors are sensitive is smaller.

Finally, one could derive a maximum likelihood estimator for the signal amplitude directly from $\ln p(d|\vartheta,\varphi,\mathcal{H}_{s},\mathcal{H}_{\delta})$. However, the corresponding expressions quickly become intractable and may be computationally prohibitive. As such, it is expected that they will be of little practical use.

In Sec. III, I obtained general expressions for the marginal likelihood given prior uncertainty on the detector calibration. I now investigate how the presence of imperfect calibration uncertainty can confound the inference of astrophysical properties in more detail. Specifically, I consider the Fisher information matrix for astrophysical parameters derived from the calibration-marginalized likelihood. In particular, I contrast this to what one might naively expect under the assumption that calibration errors were known to be absent.

The Fisher information matrix element between parameters $\vartheta_{a}$ and $\vartheta_{b}$ is defined as

$$I_{ab}\equiv\int\mathcal{D}d\,p(d|\vartheta)\left(\frac{\partial\ln p(d|\vartheta)}{\partial\vartheta_{a}}\right)\left(\frac{\partial\ln p(d|\vartheta)}{\partial\vartheta_{b}}\right)$$

(52)

and is inversely related to the lower bound on the covariance between these variables (Cramér-Rao bound). That is, larger $I_{ab}$ imply smaller covariances and vice versa. In the special case where $p(d|\vartheta)$ is Gaussian

$$\ln p(d|\vartheta)\supset-\frac{1}{2}(d-h-\mu)^{\mathrm{T}}B^{-1}(d-h-\mu)\\ -\frac{1}{2}\ln\mathrm{det}|B|$$

(53)

the Fisher information matrix becomes

$$I_{ab}=\left(\frac{\partial(h+\mu)}{\partial\vartheta_{a}}\right)^{\mathrm{T}}B^{-1}\left(\frac{\partial(h+\mu)}{\partial\vartheta_{b}}\right)\\ +\frac{1}{2}\mathrm{tr}\left[B^{-1}\frac{\partial B}{\partial\vartheta_{a}}B^{-1}\frac{\partial B}{\partial\vartheta_{b}}\right]$$

(54)

The general expression can be quite complicated, and I do not attempt to describe it in detail here. Instead, I again focus on limiting cases: when the calibration uncertainty is small and when it is large. In all cases, I assume zero-mean prior processes ($\gamma=\mu=0$) for simplicity.

For reference, in the naive case ($\Gamma=0$)

	$\displaystyle I^{(\mathrm{naive})}_{ab}$	$\displaystyle=(\partial_{a}h)C^{-1}(\partial_{b}h)$
		$\displaystyle=4\int df\frac{\mathcal{R}\{(\partial_{a}h)^{\ast}(\partial_{b}h)\}}{s}$		(55)

Finally, as has been pointed out previously (e.g., Refs. [40, 26, 47]), deviations in the waveform $h_{\alpha}\rightarrow h_{\alpha}(1+\delta)$ produce analogous terms in the likelihood: $h=f_{\alpha}h_{\alpha}\rightarrow h(1+\delta)$. If one can express the prior uncertainty on waveform errors or possible deviations from GR in terms of a GP, then the expressions derived in the context of calibration uncertainty, to a large part, hold unaltered. The interpretation of $\gamma$ and $\Gamma$ changes, but the mechanics of how to constrain them and how they might affect the inference remain the same. Just as one is able to analytically integrate out calibration uncertainties to obtain a marginal likelihood for the data, one can marginalize over uncertainty in the astrophysical signal to obtain the same. In this sense, the waveform constitutes a probabilistic map between a system’s properties (masses, spins, etc) and the expected signal from the true theory of gravity. Of particular interest may be hierarchical inference for the distribution of waveform errors assuming these errors are independently and identically distributed draws from the same GP prior for each event. This would be, in effect, a nonparametric extension to the parametrized hierarchical analysis proposed in Ref. [39] and used extensively within Refs. [6, 7].

A similar approach was explored in Ref. [48] using cubic splines to parametrize deviations from GR, although they did not perform hierarchical inference and only analyzed events separately. While splines are related to GPs (see Chapter 6 in Ref. [49]), a direct formulation of the problem in terms of GPs will likely allow for more control over prior assumptions and improved interpretability of the resulting constraints.

Although I leave a thorough investigation of the implications of GP priors for waveform errors and deviations from GR to future work, I note a few things. In this case, one expects perfect correlations in $\delta$ between IFOs (all detectors witness the same astrophysical signal) and the deviations to be drawn from the same distribution for all events (the signal of each event is determined by the same underlying theory of gravity). However, one may need to allow for separate deviations for each polarization ($\delta_{\alpha}$ may be different for each $\alpha$), mixing between polarizations ($\delta_{\alpha}$ may depend on $h_{\beta}$), and these deviations may not be best expressed as fractional changes in the waveform ($\delta_{\alpha}$ may be nonzero even when $h_{\alpha}=0\ \forall\ \alpha$).

As a final note, one may consider GP priors for both calibration ($\delta_{f}$) and waveform ($\delta_{h}$) uncertainty simultaneously. In this case, the projected strain in the detectors appears as $h=f_{\alpha}h_{\alpha}\rightarrow h(1+\delta_{f})(1+\delta_{h})$. The likelihood will contain a term like $|h\delta_{f}\delta_{h}|^{2}$, which may stymie analytic marginalization over both $\delta_{f}$ and $\delta_{h}$, although one should still be able to analytically marginalize over one of them. Again, I leave a full investigation to future work.

GW astrophysics has rapidly grown over the last few years with an ever increasing number of confidently detected compact binary coalescences [27, 28, 29, 30] and steadily improved upper limits on the rates of other sources (see, e.g., [53, 54]). This has been enabled in no small part by the remarkable understanding of the current interferometers’ responses to astrophysical strain. Indeed, in addition to providing accurate calibration for each interferometer in the network, analysts are able to carefully quantify the precision of their calibration estimates through GPs, and the incorporation of this information within astrophysical inference has become commonplace within the literature.

This manuscript exploits the fact that the dominant calibration uncertainties are described by GPs to analytically marginalize over them. Although several authors have already demonstrated numeric marginalization techniques, this analytic approach offers several advantages.

First, it may be computationally faster in practice. That is, it may be faster to compute a single inner product than to directly Monte Carlo sample over the GP calibration prior. The limiting step will likely be the inversion of $\Gamma$ and $A^{-1}=H^{\mathrm{T}}C^{-1}H+\Gamma^{-1}$ required to compute $B^{-1}$ for each signal model. This step may scale as $O(N_{\mathrm{frq}}^{3})$. Monte Carlo integration could be faster if very few Monte Carlo samples are needed for the integral to converge. This will likely be the case only when the calibration prior is already so precise that it is not changed significantly by the addition of astrophysical data (the current situation for all signals detected to date). That is, direct numeric integration may only be faster when the calibration uncertainties are too small to make much of a difference anyway. At the same time, one may also be able to approximate $A\approx\Gamma$ in this limit, removing the need to invert $\Gamma$ and $A$ and thereby reducing the computational cost from $O(N_{\mathrm{frq}}^{3})$ to $O(N_{\mathrm{frq}}^{2})$.

What’s more, analytic expressions allow for greater interpretability of the resulting conditioned and marginal distributions. For example, I am able to easily take several limits and explore the impact of large/small calibration uncertainties and correlations. I also show how to trivially extend the formalism to include networks of multiple detectors and multiple astrophysical signals, including possible correlations between the calibration errors witnessed in each detector/event pair. Equivalent investigations via numeric integration are possible, but may be more difficult to interpret and more costly to perform.

Specifically, I find that calibration uncertainties always reduce the measurability of astrophysical parameters. This can be seen in several ways. I explore the impact of calibration uncertainties on search sensitivity through the behavior of the standard matched-filter S/N, finding larger variance in the presence of calibration uncertainty and, therefore, reduced search sensitivity. Monte Carlo simulations of search sensitivity [55, 56, 57] may be able to exploit this analytic marginalization to better account for calibration uncertainties without recourse to more direct numerical treatments. However, perhaps the most striking example can be seen through the change in the Fisher information matrix elements.

The loss of information is a projection of the signal onto the calibration uncertainty’s prior covariance matrix. That is, the additional uncertainty in astrophysical parameters is directly related to whether the types of changes in the waveform induced by changing those parameters are common calibration errors. Additionally, it is not necessary to constrain the calibration uncertainty to be small at all frequencies for accurate inferences of individual signals. Instead, one need only constrain particular combinations of frequencies in order to minimize the information lost. This also suggests an algorithm to determine the frequencies at which one should measure the calibration more precisely in order to best reduce the uncertainty in astrophysical parameters. As a rule of thumb, these are (unsurprisingly) the frequencies at which the current calibration uncertainty most hurts the inference.

When calibration errors are large (and correlations between frequencies are not perfect), the uncertainty in astrophysical parameters becomes independent of the signal amplitude. The parameters of louder signals are not better constrained than the parameters of quieter signals. This sets a limit on the precision to which the exceptionally loud astrophysical signals can be measured. The relevant comparison is whether the S/N is greater or smaller than the inverse of the calibration error’s standard deviation ($\rho\lessgtr\sigma_{\Gamma}^{-1}$). Including an approximate proportionality constant, I find that $1.4+1.4\,M_{\odot}$ binaries with $\rho\gtrsim 10^{3}$ require calibration uncertainties $\lesssim 3\%$.