Null-stream-based Bayesian Unmodeled Framework to Probe Generic Gravitational-wave Polarizations

The additive noise observation model of GW with all possible polarization modes in a $D$-detector network writes

$$\boldsymbol{d}(t;\Delta\boldsymbol{t})=\boldsymbol{F}(\alpha,\delta,\psi,t)\boldsymbol{h}(t)+\boldsymbol{n}(t;\Delta\boldsymbol{t})$$

(1)

where

$$\boldsymbol{d}(t;\Delta\boldsymbol{t})=\begin{bmatrix}d_{1}(t+\Delta t_{1})\\ d_{2}(t+\Delta t_{2})\\ \vdots\\ d_{D}(t+\Delta t_{D})\end{bmatrix}$$

(2)

is the observed strain outputs at shifted times,

$$\begin{split}\boldsymbol{F}(\alpha,\delta,\psi,t)&=\begin{bmatrix}\boldsymbol{f}_{+}&\boldsymbol{f}_{\times}&\boldsymbol{f}_{b}&\boldsymbol{f}_{l}&\boldsymbol{f}_{x}&\boldsymbol{f}_{y}\end{bmatrix}\\ &=\begin{bmatrix}F_{1}^{+}&F_{1}^{\times}&F_{1}^{b}&F_{1}^{l}&F_{1}^{x}&F_{1}^{y}\\ F_{2}^{+}&F_{2}^{\times}&F_{2}^{b}&F_{2}^{l}&F_{2}^{x}&F_{2}^{y}\\ \vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\ F_{D}^{+}&F_{D}^{\times}&F_{D}^{b}&F_{D}^{l}&F_{D}^{x}&F_{D}^{y}\end{bmatrix}\end{split}$$

(3)

is the beam pattern matrix,

$$\boldsymbol{h}(t)=\begin{bmatrix}h_{+}(t)&h_{\times}(t)&h_{b}(t)&h_{l}(t)&h_{x}(t)&h_{y}(t)\end{bmatrix}^{T}$$

(4)

are the polarization modes where $T$ denotes transpose and the plus mode, cross mode, scalar breathing mode, scalar longitudinal mode, vector x mode and vector y mode are labeled with the subscripts $+$, $\times$, $b$, $l$, $x$ and $y$ respectively,

$$\boldsymbol{n}(t;\Delta\boldsymbol{t})=\begin{bmatrix}n_{1}(t+\Delta t_{1})\\ n_{2}(t+\Delta t_{2})\\ \vdots\\ n_{D}(t+\Delta t_{D})\end{bmatrix}$$

(5)

is the detector noise at shifted times, $\alpha$, $\delta$ are the right ascension and declination of the source location respectively, $\psi$ is the polarization angle, $\Delta\boldsymbol{t}=\{\Delta t_{1},\Delta t_{2},...,\Delta t_{D}\}$ where

$$\Delta t_{j}=-\frac{\boldsymbol{r}_{j}\cdot\hat{\boldsymbol{N}}}{c}$$

(6)

is the time delay of the signal arrival at each detector with reference to the Earth center, $\boldsymbol{r}_{j}$ is the coordinates of the $j^{\text{th}}$ detector in the Earth-centered coordinate system,

$$\hat{\boldsymbol{N}}=\begin{bmatrix}\cos(\delta)\cos(t_{\text{gmst}}-\alpha)\\ -\cos(\delta)\sin(t_{\text{gmst}}-\alpha)\\ \sin(\delta)\end{bmatrix}$$

(7)

where $t_{\text{gmst}}$ is the Greenwich mean sidereal time, and $c$ is the speed of gravity. For a two-arm interferometer, the beam pattern functions take the following forms when the two arms are directed to the x-axis and y-axis respectively Nishizawa et al. (2009):

$$\begin{split}F_{+}(\theta,\phi,\psi)=&\frac{1}{2}(1+\cos^{2}\theta)\cos{2\phi}\cos{2\psi}\\ &-\cos\theta\sin{2\phi}\sin{2\psi}\end{split}$$

(8)

$$\begin{split}F_{\times}(\theta,\phi,\psi)=&-\frac{1}{2}(1+\cos^{2}\theta)\cos{2\phi}\sin{2\psi}\\ &-\cos{\theta}\sin{2\phi}\cos{2\psi}\end{split}$$

(9)

$$F_{x}(\theta,\phi,\psi)=\sin{\theta}(\cos{\theta}\cos{2\phi}\cos{\psi}-\sin{2\phi}\sin{\psi})$$

(10)

$$F_{y}(\theta,\phi,\psi)=-\sin{\theta}(\cos\theta\cos{2\phi}\sin{\psi}+\sin{2\phi}\cos{\psi})$$

(11)

$$F_{b}(\theta,\phi)=-\frac{1}{2}\sin^{2}\theta\cos{2\phi}$$

(12)

$$F_{l}(\theta,\phi)=\frac{1}{\sqrt{2}}\sin^{2}\theta\cos{2\phi}$$

(13)

where $\theta$ and $\phi$ are the polar angle and the azimuthal angle of the source location in the Earth-centered coordinate system respectively, and $\psi$ is the polarization angle. Fig. 1 demonstrates the effect of different polarization modes on a ring of free-falling test particles.

In frequency domain, the observation model writes

$$\tilde{\boldsymbol{d}}[k;\Delta\boldsymbol{t}]=\boldsymbol{F}(\alpha,\delta,\psi,t_{\text{event}})\tilde{\boldsymbol{h}}[k]+\tilde{\boldsymbol{n}}[k;\Delta\boldsymbol{t}]$$

(14)

where

$$\tilde{\boldsymbol{d}}[k;\Delta\boldsymbol{t}]=\begin{bmatrix}\tilde{d}_{1}[k]e^{i2\pi k\Delta t_{1}/T}\\ \tilde{d}_{2}[k]e^{i2\pi k\Delta t_{2}/T}\\ \vdots\\ \tilde{d}_{D}[k]e^{i2\pi k\Delta t_{D}/T}\end{bmatrix}$$

(15)

is the Fourier transform (see Appendix A) of the time-shifted strain outputs, $T$ is the duration of data in second, $\boldsymbol{F}(\alpha,\delta,\psi,t)\approx\boldsymbol{F}(\alpha,\delta,\psi,t_{\text{event}})$ is approximated to be constant over the data segment which is valid for a transient signal and is evaluated at the event time $t_{\text{event}}$,

$$\tilde{\boldsymbol{h}}[k]=\begin{bmatrix}\tilde{h}_{+}[k]&\tilde{h}_{\times}[k]&\tilde{h}_{b}[k]&\tilde{h}_{l}[k]&\tilde{h}_{x}[k]&\tilde{h}_{y}[k]\end{bmatrix}^{T}$$

(16)

is the Fourier transform of the polarization modes, and

$$\tilde{\boldsymbol{n}}[k;\Delta\boldsymbol{t}]=\begin{bmatrix}\tilde{n}_{1}[k]e^{i2\pi k\Delta t_{1}/T}\\ \tilde{n}_{2}[k]e^{i2\pi k\Delta t_{2}/T}\\ \vdots\\ \tilde{n}_{D}[k]e^{i2\pi k\Delta t_{D}/T}\end{bmatrix}$$

(17)

is the Fourier transform of the time-shifted noise in each detector. Assume the detector noise follows a stationary Gaussian distribution and the noise is independent between detectors, we have

$$\braket{\tilde{n}_{j^{\prime}}^{*}[k^{\prime}]\tilde{n}_{j}[k]}=\frac{1}{2\Delta f}S_{j}[k]\delta_{k^{\prime}k}\delta_{j^{\prime}j}$$

(18)

where $\braket{\cdot}$ denotes expectation over the noise realizations, $S_{j}[k]$ is the one-sided power spectral density (PSD) of the noise in detector $j$, $\delta_{jk}$ denotes Kronecker delta, and $\Delta f$ is the frequency resolution of the discrete Fourier transform.

For brevity, in the following sections we write the time-shifted data as the time-shift operator $\boldsymbol{\mathcal{T}}\left(\cdot;\Delta\boldsymbol{t}\right):\mathbb{C}^{D\times K}\rightarrow\mathbb{C}^{D\times K}$ applied on the unshifted data defined by

$$\boldsymbol{\mathcal{T}}(\tilde{\boldsymbol{x}};\Delta\boldsymbol{t})=\tilde{\boldsymbol{x}}\odot\mathcal{S}(\Delta\boldsymbol{t})$$

(19)

where

$$\tilde{\boldsymbol{x}}=\begin{bmatrix}\tilde{x}_{1}[1]&\tilde{x}_{1}[2]&\ldots&\tilde{x}_{1}[K]\\ \tilde{x}_{2}[1]&\tilde{x}_{2}[2]&\ldots&\tilde{x}_{2}[K]\\ \vdots&\vdots&\ddots&\vdots\\ \tilde{x}_{D}[1]&\tilde{x}_{D}[2]&\ldots&\tilde{x}_{D}[K]\end{bmatrix}$$

(20)

is the unshifted data, $(\mathcal{S}(\Delta\boldsymbol{t}))_{jk}=e^{i2\pi(k-1)\Delta t_{j}/T}$ where $j=1,2,...,D$ and $k=1,2,...,K$ or more explicitly

$$\begin{split}&\mathcal{S}(\Delta\boldsymbol{t})\\ &=\begin{bmatrix}e^{i2\pi(0)\Delta t_{1}/T}&e^{i2\pi(1)\Delta t_{1}/T}&\ldots&e^{i2\pi(K-1)\Delta t_{1}/T}\\ e^{i2\pi(0)\Delta t_{2}/T}&e^{i2\pi(1)\Delta t_{2}/T}&\ldots&e^{i2\pi(K-1)\Delta t_{2}/T}\\ \vdots&\vdots&\ddots&\vdots\\ e^{i2\pi(0)\Delta t_{D}/T}&e^{i2\pi(1)\Delta t_{D}/T}&\ldots&e^{i2\pi(K-1)\Delta t_{D}/T}\end{bmatrix}\end{split}\,,$$

(21)

$\odot$ denotes the elementwise product, and $K$ is the number of frequency bins.

Null stream Gürsel and Tinto (1989) is a specific linear combination of strain outputs from a network of detectors such that the linear combination only contains noise regardless of the GW waveform. Let us consider a tensorial signal as follows

$\displaystyle\tilde{\boldsymbol{s}}(f)=\begin{bmatrix}\boldsymbol{f}_{+}&\boldsymbol{f}_{\times}\end{bmatrix}\begin{bmatrix}\tilde{h}_{+}(f)\\ \tilde{h}_{\times}(f)\end{bmatrix}\,,$

(22)

the signal $\tilde{\boldsymbol{s}}(f)$ could be regarded as living in the subspace spanned by $\boldsymbol{f}_{+}$ and $\boldsymbol{f}_{\times}$ as shown in Fig. 2. The hyperplane spanned by $\boldsymbol{f}_{+}$ and $\boldsymbol{f}_{\times}$ in the figure defines all possible strain measurement of a tensorial signal in each detector. We could project away the hyperplane to remove any tensorial signal regardless of the waveform of $\tilde{h}_{+}(f)$ and $\tilde{h}_{\times}(f)$. The residual remains in the null space of $\boldsymbol{f}_{+}$ and $\boldsymbol{f}_{\times}$ is the so-called null stream which contains noise only.

One application of null stream is to distinguish between coincident GW transients and coincident noise glitches. If the observed coherent excess power is a genuine GW transient, one should obtain a null stream consistent with noise. For noise glitches, in general, the resultant null stream would still contain excess power. Null stream is therefore used in conducting vetos Wen and Schutz (2005); Ajith et al. (2006); Chatterji et al. (2006) and GW searches associated with gamma-ray bursts and other astrophysical triggers Sutton et al. (2010). In this section, we review the method to construct a null stream. The matrix representation follows from Ref. Sutton et al. (2010). The only difference is that we include the time-shift operation explicitly in the construction of null stream since we do not know the source location a priori in contrast to the follow-up analysis of gamma-ray bursts in Ref. Sutton et al. (2010).

Suppose we observe a signal at time $t_{\text{event}}$, the frequency domain observation model writes

$$\tilde{\boldsymbol{d}}[k]=\boldsymbol{\mathcal{T}}\left(\boldsymbol{F}(\alpha,\delta,\psi,t_{\text{event}})\tilde{\boldsymbol{h}};-\Delta\boldsymbol{t}\right)[k]+\tilde{\boldsymbol{n}}[k]\,.$$

(23)

Before performing the null projection, we first whiten the data since this would make it easier to handle the statistics of the residual noise. Dividing both sides of Eq. (23) by $\sqrt{\frac{1}{2\Delta f}S[k]}$, we have

$$\tilde{\boldsymbol{d}}_{w}[k]=\boldsymbol{\mathcal{T}}\left(\boldsymbol{F}_{w}(\alpha,\delta,\psi,t_{\text{event}})\tilde{\boldsymbol{h}};-\Delta\boldsymbol{t}\right)[k]+\tilde{\boldsymbol{n}}_{w}[k]$$

(24)

where

$$\tilde{d}_{w,j}[k]=\frac{\tilde{d}_{j}[k]}{\sqrt{\frac{1}{2\Delta f}S_{j}[k]}}\,,$$

(25)

$$F_{w,j}[k]=\frac{F_{j}}{\sqrt{\frac{1}{2\Delta f}S_{j}[k]}}$$

(26)

and

$$\tilde{n}_{w,j}[k]=\frac{\tilde{n}_{j}[k]}{\sqrt{\frac{1}{2\Delta f}S_{j}[k]}}$$

(27)

where $F_{j}$ represents the $j^{\text{th}}$ row of $\boldsymbol{F}$, $\boldsymbol{F}_{w}\in\mathbb{C}^{D\times M\times K}$ is the noise-weighed beam pattern matrix, and $M$ is the number of column vectors in $\boldsymbol{F}$. The construction of null stream consists of two steps: 1) time-shifting the data with the time delay implied by the source location and event time i.e. $\Delta\boldsymbol{t}=\Delta\boldsymbol{t}(\alpha,\delta,t_{\text{event}})$, and 2) performing the null projection. Denote the null operator as $\boldsymbol{P}:\mathbb{C}^{D\times D\times K}\rightarrow\mathbb{C}^{D\times D\times K}$ defined by

$$\boldsymbol{P}\left({}\cdot{};\alpha,\delta,\psi,t_{\text{event}}\right)=(\boldsymbol{I}-\boldsymbol{F}_{w}(\boldsymbol{F}_{w}^{\dagger}\boldsymbol{F}_{w})^{-1}\boldsymbol{F}_{w}^{\dagger})\boldsymbol{\mathcal{T}}\left(\cdot;\Delta\boldsymbol{t}\right)$$

(28)

where $\boldsymbol{I}\in\mathbb{C}^{D\times D\times K}$ is a stack of $K$ $D\times D$ identity matrices and $\dagger$ denotes conjugate transpose. The matrix operation in Eq. (28) is performed frequency-wise. One should notice that the null operator is completely determined by the constituent $\boldsymbol{f}_{m}$ of $\boldsymbol{F}$ and the set of parameters $\{\alpha,\delta,\psi,t_{\text{event}}\}$ is independent of the waveform $\tilde{\boldsymbol{h}}$. When the set of parameters are correctly specified, one could verify that the null operator projects away $\tilde{\boldsymbol{h}}$ without requiring any knowledge of $\tilde{\boldsymbol{h}}$ itself, and the residual $\tilde{\boldsymbol{z}}\in\mathbb{C}^{D\times K}$ is

$$\begin{split}\tilde{\boldsymbol{z}}&\coloneqq\boldsymbol{P}(\tilde{\boldsymbol{d}}_{w};\alpha,\delta,\psi,t_{\text{event}})\\ &=(\boldsymbol{I}-\boldsymbol{F}_{w}(\boldsymbol{F}_{w}^{\dagger}\boldsymbol{F}_{w})^{-1}\boldsymbol{F}_{w}^{\dagger})\boldsymbol{\mathcal{T}}\left(\tilde{\boldsymbol{d}}_{w};\Delta\boldsymbol{t}\right)\\ &=(\boldsymbol{I}-\boldsymbol{F}_{w}(\boldsymbol{F}_{w}^{\dagger}\boldsymbol{F}_{w})^{-1}\boldsymbol{F}_{w}^{\dagger})\left(\boldsymbol{F}_{w}\tilde{\boldsymbol{h}}+\boldsymbol{\mathcal{T}}(\tilde{\boldsymbol{n}};\Delta\boldsymbol{t})\right)\\ &=(\boldsymbol{I}-\boldsymbol{F}_{w}(\boldsymbol{F}_{w}^{\dagger}\boldsymbol{F}_{w})^{-1}\boldsymbol{F}_{w}^{\dagger})\boldsymbol{\mathcal{T}}\left(\tilde{\boldsymbol{n}}_{w};\Delta\boldsymbol{t}\right)\end{split}$$

(29)

which is the residual noise in the lower dimensional space. To understand the reduced dimensionality, one should be reminded that null projection removes the data on the hyperplane spanned by the beam pattern vectors $\boldsymbol{f}_{m}$, and we may call the hyperplane as the GW space. Only the data in the null space that is orthogonal to the GW space survive and it is the so-called null stream. One could perform a rotation $\mathcal{R}$ to rotate the residuals to align along the principal axes to visualize the reduced dimensionality (see Appendix B.1). Although we write the null operator in Eq. (28) as a function of the polarization angle $\psi$, the null operator is in fact independent of the polarization angle if we consider non-degenerate polarization modes (see Appendix B.2).

In this section, we present the basis formulation to test GW polarizations. Suppose we would like to compare two polarization hypotheses $\mathcal{H}_{1}$ and $\mathcal{H}_{2}$ which beam pattern matrices are $\boldsymbol{F}_{1}\in\mathbb{R}^{D\times m_{1}}$ and $\boldsymbol{F}_{2}\in\mathbb{R}^{D\times m_{2}}$ respectively where $m_{1}<m_{2}$ without loss of generality. One should notice the resultant null streams have unequal dimensionalities. To have a fair comparison between different models, the residuals have to be of the same dimensionality. Instead of performing the most generic null projection with $\boldsymbol{F}_{2}\in\mathbb{R}^{D\times m_{2}}$, one could perform the null projection to remove a smaller subspace to match the dimensionality of $\boldsymbol{F}_{1}\in\mathbb{R}^{D\times m_{1}}$. We propose the basis formulation to restrict the null projection to remove a smaller subspace constrained by the polarization basis modes. To illustrate the idea, we may consider a scalar-tensorial signal model with the scalar breathing mode in addition to the two tensorial modes

$$\tilde{\boldsymbol{s}}[k]=\boldsymbol{f}_{+}\tilde{h}_{+}[k]+\boldsymbol{f}_{\times}\tilde{h}_{\times}[k]+\boldsymbol{f}_{b}\tilde{h}_{b}[k]\,.$$

(30)

The most generic null projection would remove three dimensions from the data since we have three linearly independent beam pattern vectors. Without loss of generality, we choose the $+$ mode as the basis, and we can always decompose any polarization mode $\tilde{h}_{m}[k]$ into the parallel component and the orthogonal component with respect to the basis i.e.

$$\tilde{h}_{m}[k]=C_{\parallel}^{m}\tilde{h}_{\parallel}[k]+C_{\perp}^{m}\tilde{h}_{\perp}[k]$$

(31)

where $\tilde{h}_{\parallel}[k]=\tilde{h}_{+}[k]$ and $C_{\parallel},C_{\perp}\in\mathbb{C}$. With this formulation, we can rewrite Eq. (30) as

$$\begin{split}&\tilde{\boldsymbol{s}}[k]\\ &=(\boldsymbol{f}_{+}+C_{\parallel}^{\times}\boldsymbol{f}_{\times}+C_{\parallel}^{b}\boldsymbol{f}_{b})\tilde{h}_{\parallel}[k]+(C_{\perp}^{\times}\boldsymbol{f}_{\times}+C_{\perp}^{b}\boldsymbol{f}_{b})\tilde{h}_{\perp}[k]\\ &=\boldsymbol{f}_{\parallel}\tilde{h}_{\parallel}[k]+\boldsymbol{f}_{\perp}\tilde{h}_{\perp}[k]\end{split}$$

(32)

where

$$\begin{split}&\tilde{h}_{+}[k]=\tilde{h}_{\parallel}[k]\\ &\tilde{h}_{\times}[k]=C_{\parallel}^{\times}\tilde{h}_{\parallel}[k]+C_{\perp}^{\times}\tilde{h}_{\perp}[k]\\ &\tilde{h}_{b}[k]=C_{\parallel}^{b}\tilde{h}_{\parallel}[k]+C_{\perp}^{b}\tilde{h}_{\perp}[k]\end{split}$$

(33)

and

$$\begin{split}&\boldsymbol{f}_{\parallel}=\boldsymbol{f}_{+}+C_{\parallel}^{\times}\boldsymbol{f}_{\times}+C_{\parallel}^{b}\boldsymbol{f}_{b}\\ &\boldsymbol{f}_{\perp}=C_{\perp}^{\times}\boldsymbol{f}_{\times}+C_{\perp}^{b}\boldsymbol{f}_{b}\,.\end{split}$$

(34)

We could then perform the (partial) null projection to remove the smaller subspace spanned by $\boldsymbol{f}_{\parallel}$ only. Qualitatively, we have reduced the most generic scalar-tensor GW space spanned by $\boldsymbol{f}_{+}$, $\boldsymbol{f}_{\times}$ and $\boldsymbol{f}_{b}$ to the scalar-tensor GW subspace spanned by $\boldsymbol{f}_{\parallel}$ where the polarization modes are described by a single basis. With this formulation, the dimensionality of the GW space is then determined by the number of basis modes in contrast to that in the generic null projection determined by the number of beam pattern vectors. In other words, this formulation allows us to capture the non-tensorial components parallel to the basis mode(s) (in general we can have more than one basis mode). One attractive feature of this formulation is that we only need to choose the basis modes of some polarization types, but we do not need to assume the waveform of the basis modes.

In the basis formulation, we have the freedom to choose the basis mode(s) without explicitly modeling them. Also, we have the freedom to choose the number of basis mode(s). In terms of testing GR, we want to compare the likelihood between the tensor hypothesis and the non-tensor hypotheses. Since there are two polarization modes in tensorial polarizations, we could choose either one or two basis mode(s) to parameterize the tensor hypothesis. The competing hypotheses are then formulated with the same number of basis modes accordingly. For brevity, we may call the one-basis-mode analysis as the $L=1$ analysis, and the two-basis-mode analysis as the $L=2$ analysis where $L$ denotes the number of basis mode(s). In the $L=2$ analysis of the tensor hypothesis, the signal space contains all possible tensorial polarizations, and the waveforms of the plus and cross polarization can take any independent shape. In the $L=1$ analysis, an additional structure is imposed on the plus and cross polarizations i.e. they are linearly dependent in the Fourier domain, and therefore the signal space is more restricted compared to that in the $L=2$ analysis.

In general, a higher number of basis modes gives a greater degree of freedom to fit the data but it also reduces the distinguishing power between different polarization hypotheses. The major reason is that with a higher number of basis modes we remove a larger subspace from the data, and therefore the residual has a lower dimensionality. One should be reminded that the polarization subspaces spanned by the beam pattern vectors are in general not orthogonal to each other. For example, removing the tensorial subspace could also remove part of the non-tensorial component (if any) of the signal. The polarization subspaces, in general, have a greater overlap in the $L=2$ analysis than that in the $L=1$ analysis, and one should expect the analysis with a higher number of basis modes has a weaker distinguishing power between the polarization hypotheses.

Although the $L=1$ analysis, in general, has a stronger distinguishing power than the $L=2$ analysis, the $L=2$ analysis helps to capture the orthogonal polarization mode(s) that is/are possibly missed in the $L=1$ analysis. In practice, we can always perform both analyses and examine the consistency of the results.

In this section, we present the parameterization of the $L=1$ analysis and the $L=2$ analysis. Each polarization hypothesis is characterized by the beam pattern vector $\boldsymbol{f}_{m}$ involved in the construction of the null operator. To facilitate the mathematical expressions, we may define $\mathbb{F}_{p}$ as the set of labels of polarizations corresponding to the polarization hypothesis $\mathcal{H}_{p}$ as follows

$$\begin{split}&\mathbb{F}_{T}=\{+,\times\}\\ &\mathbb{F}_{V}=\{x,y\}\\ &\mathbb{F}_{S}=\{b,l\}\\ &\mathbb{F}_{TV}=\mathbb{F}_{T}\cup\mathbb{F}_{V}\\ &\mathbb{F}_{TS}=\mathbb{F}_{T}\cup\mathbb{F}_{S}\\ &\mathbb{F}_{VS}=\mathbb{F}_{V}\cup\mathbb{F}_{S}\\ &\mathbb{F}_{TVS}=\mathbb{F}_{T}\cup\mathbb{F}_{V}\cup\mathbb{F}_{S}\end{split}$$

(35)

where $+$, $\times$, $b$, $l$, $x$ and $y$ denotes the labels for the plus polarization, cross polarization, scalar breathing polarization, scalar longitudinal polarization, vector $x$ polarization, and vector $y$ polarization respectively.

Since the model parameters $\boldsymbol{\theta}=\{\alpha$, $\delta,\psi,\boldsymbol{\mathcal{A}},\boldsymbol{\varphi}\}$ are not known (we take the event time $t_{\text{event}}$ from the search pipelines as known), we adopt a Bayesian analysis to marginalize the unknown parameters to obtain the Bayesian evidence of the competing hypotheses.

To improve the sensitivity, we analyze the data in the time-frequency domain and target the region spanned by the candidate signal through time-frequency clustering. We use the Wilson-Daubechies-Meyer (WDM) time-frequency transform Necula et al. (2012) because of its superior spectral leakage control. Given a multi-detector discrete time series $\boldsymbol{x}\in\mathbb{R}^{D\times N}$ where $D$ is the number of detectors and $N$ is the number of time bins, denote the time-frequency transform as $\textbf{TF}:\mathbb{R}^{D\times N}\rightarrow\mathbb{R}^{D\times J\times K}$ where $J$ is the number of time bins and $K$ is the number of frequency bins which performs the WDM time-frequency transform on the time series of each detector independently. The corresponding time-frequency representation is denoted with a subscript as $\boldsymbol{x}_{\text{TF}}$.

One should notice that the construction of the null operator involves the noise-weighed beam pattern function matrix $\boldsymbol{F}_{w}$ in Eq. (26). The noise PSD of LIGO-Virgo interferometers typically exhibits sharp spectral lines due to the resonance of wires in suspension, electrical supply power line, and injected calibration lines Littenberg and Cornish (2015). To reduce the spectral leakage, we compute the residual $\tilde{\boldsymbol{z}}$ in Eq. (29) in the frequency domain which has a higher frequency resolution than the time-frequency representation. Then, we perform the inverse Fourier transform $\mathcal{F}^{-1}$ to transform the residual back to the time domain and then transform the time domain residual to the time-frequency domain by TF. As shown in Eq. (29), when the null operator $\boldsymbol{P}$ is constructed with the correct parameters, the residual $\tilde{\boldsymbol{z}}$ is only noise in the lower dimensional space regardless of $\tilde{\boldsymbol{h}}$. The time-frequency representation of the residual $\boldsymbol{z}_{\text{TF}}$ is also only noise with the null energy $E_{\text{null}}$ Sutton et al. (2010) defined as

$$E_{\text{null}}(\boldsymbol{z}_{\text{TF}})=\sum_{j=1}^{D}\sum_{\{k,l\}\in\mathcal{C}}\left\lVert z_{\text{TF},j}[k,l]\right\rVert^{2}$$

(45)

where $z_{\text{TF},j}[k,l]$ denotes the coefficient of the time-frequency representation $\boldsymbol{z}_{\text{TF}}$ corresponding to the $k^{\text{th}}$ time index and $l^{\text{th}}$ frequency index of the residual of the $j^{\text{th}}$ detector, and $\mathcal{C}$ denotes the set of time-frequency indices that the candidate signal occupies on the time-frequency plane follows the $\chi^{2}$ distribution with a degree of freedom $(D-L)N$ where $L$ is the number of basis modes and $N$ is the number of time-frequency pixels that the candidate signal spans when the null operator is correctly constructed. The degree of freedom is half of that in Ref. Sutton et al. (2010) since the WDM coefficients are real but in Ref. Sutton et al. (2010) Fourier transform is used and the real and imaginary components double the degree of freedom.

The likelihood is then given by

$$p(\boldsymbol{d}|\alpha,\delta,\psi,\boldsymbol{\mathcal{A}},\boldsymbol{\varphi},\mathcal{H})=\chi^{2}_{\text{DoF}}(E_{\text{null}}(\boldsymbol{z}_{\text{TF}}))$$

(46)

where $\boldsymbol{z}_{\text{TF}}=\textbf{TF}(\mathcal{F}^{-1}(\boldsymbol{P}_{\mathcal{H}}(\tilde{\boldsymbol{d}};\alpha,\delta,\psi,\boldsymbol{\mathcal{A}},\boldsymbol{\varphi})))$, $\boldsymbol{P}_{\mathcal{H}}$ denotes the null operator constructed with the $\boldsymbol{f}_{m}$ implied by the hypothesis $\mathcal{H}$ and the number of basis modes $L$, and $\chi_{\text{DoF}}^{2}(\cdot)$ is the $\chi^{2}$ probability density function with the degree of freedom DoF implied by the hypothesis $\mathcal{H}$. The Bayesian evidence is then evaluated by marginalizing all parameters

$$\begin{split}&p(\boldsymbol{d}|\mathcal{H})\\ &=\int p(\boldsymbol{d}|\alpha,\delta,\psi,\boldsymbol{\mathcal{A}},\boldsymbol{\varphi},\mathcal{H})p(\alpha,\delta,\psi,\boldsymbol{\mathcal{A}},\boldsymbol{\varphi}|\mathcal{H})d\alpha d\delta d\psi d\boldsymbol{\mathcal{A}}d\boldsymbol{\varphi}\end{split}$$

(47)

where $p(\alpha,\delta,\psi,\boldsymbol{A},\boldsymbol{\varphi}|\mathcal{H})$ is the prior distribution of the parameters given $\mathcal{H}$. The posterior odds between $\mathcal{H}_{1}$ and $\mathcal{H}_{2}$ defined as

$$\begin{split}\mathcal{O}_{\mathcal{H}_{2}}^{\mathcal{H}_{1}}&=\frac{p(\mathcal{H}_{1}|\boldsymbol{d})}{p(\mathcal{H}_{2}|\boldsymbol{d})}\\ &=\frac{p(\boldsymbol{d}|\mathcal{H}_{1})}{p(\boldsymbol{d}|\mathcal{H}_{2})}\times\frac{p(\mathcal{H}_{1})}{p(\mathcal{H}_{2})}\end{split}$$

(48)

that decribes the ratio of probabilities between the two hypotheses $\mathcal{H}_{1}$ and $\mathcal{H}_{2}$ given the observed data $\boldsymbol{d}$, $\frac{p(\boldsymbol{d}|\mathcal{H}_{1})}{p(\boldsymbol{d}|\mathcal{H}_{2})}$ is the ratio of model evidences and is also called the Bayes factor, and $\frac{p(\mathcal{H}_{1})}{p(\mathcal{H}_{2})}$ is the prior odds between the hypotheses $\mathcal{H}_{1}$ and $\mathcal{H}_{2}$ which describes the a priori belief of the ratio of probabilities of the two hypotheses. When we are ignorant about the relative probability between the competing hypothesis, we take the prior odds to be $1$, and then the posterior odds would be equal to the Bayes factor. In the following sections, we use the Bayes factor as the detection statistic. It should be understood that Eq. (47) is presented for generality, $\{\psi,\boldsymbol{A},\boldsymbol{\varphi}\}$ are not present in the integral when the hypothesis $\mathcal{H}$ does not have the extra parameterization.

As a follow-up analysis, in the analysis pipeline, we require the time-frequency cluster found to span across the event time reported by the search pipeline. This could fail when the signal is too weak, and there are loud noise transients around the event. If the cluster does not span across the event time, the analysis is performed in the frequency domain, and the whole data chunk is being analyzed. The null energy is then computed with the frequency domain residual $\tilde{\boldsymbol{z}}$

$$E_{\text{null}}(\tilde{\boldsymbol{z}})=\sum_{j=1}^{D}\sum_{k=1}^{K}\left\lVert\tilde{z}_{j}[k]\right\rVert^{2}$$

(49)

where $\tilde{z}_{j}[k]$ denotes the residual at frequency bin $k$ of the $j^{\text{th}}$ detector, $D$ denotes the number of detectors and $K$ denotes the number of frequency bins. $2E_{\text{null}}$ then follows the $\chi^{2}$ distribution with degree of freedom $2(D-M)K$ when the null operator is correctly constructed.

One would still need to pay attention to the orthogonal polarization component $\tilde{h}_{\perp}$ with the partial null projection. In terms of testing GR, the orthogonal polarization component must not be significant in the tensor hypothesis, and otherwise, this could lead to false GR violation. For CBC signals, the plus- and the cross-polarization modes in the dominant 22-mode only differ by the amplitude (due to the inclination angle) and the phase by $\pi/2$ which are frequency-independent Creighton and Anderson (2011). The addition of the higher-order harmonics only adds subdominant corrections. The plus mode and the cross mode are still well described by a single basis mode, and the orthogonal component is vanishingly small. In the example of the $L=1$ analysis presented in Sec. II.4.1, the orthogonal component would be significant if the dipole radiation only enters the non-tensorial polarizations but not the tensorial polarizations like the Brans-Dicke theory Chatziioannou et al. (2012b). We would argue that detecting the non-tensorial component parallel to the basis modes is sufficient for us to detect a GR violation. However, we also want to know how much uncaptured $\tilde{h}_{\perp}$ are in the data to understand how well the model can explain the data. This is important for us to understand the systematics of the ranking of the non-tensor hypotheses if we observe a deviation from GR.

With the same spirit of the residuals test Abbott et al. (2019b); Collaboration et al. (2020), we analyze whether the null energy defined in Eq. (45) with the maximum-likelihood parameters are consistent with noise or not. The consistency could be quantified using the plug-in p-value Demortier (2007) defined by

$$p_{\text{plug-in}}(E_{\text{null}})=\int_{E_{\text{null}}}^{\infty}\chi^{2}_{\text{DoF}}(E)dE$$

(50)

where $E_{\text{null}}$ is the null energy with the maximum-likelihood parameters. One should notice that the plug-in p-value does not distribute uniformly between $0$ and $1$. Instead, it should always be around $p_{\text{plug-in}}(\text{DoF}-2)$ (assume $\text{DoF}\geq 2$) since the mode of a $\chi^{2}$ distribution with a degree of freedom DoF is $\text{DoF}-2$ that is the $E_{\text{null}}$ with the highest possible likelihood. A significantly smaller $p_{\text{plug-in}}(E_{\text{null}})$ than the $p_{\text{plug-in}}(\text{DoF}-2)$ indicates the presence of an orthogonal polarization component. The $p_{\text{plug-in}}$ could therefore be served as a diagnostic tool, but one shoule be reminded that it cannot be interpreted as the usual p-value since it does not distriute uniformly between $0$ and $1$ under the null hypothesis. A more interpretable statistic could be obtained using the double bootstrap method to calculate the adjusted plug-in p-value Demortier (2007); Davison and Hinkley (1997) that distributes uniformly between $0$ and $1$.

Calibration errors are present in instruments Cahillane et al. (2017); Accadia et al. (2010). Although a previous study Vitale et al. (2012) shows that calibration-induced errors of Advanced LIGO and Advanced Virgo are not a significant detriment to accurate parameter estimation, including the effects of calibration errors has been a standard practice in the LIGO Scientific Collaboration and the Virgo Collaboration Abbott et al. (2019a, 2020a) to improve the accuracy of results and for completeness. The details of including effects of calibration errors into the framework are discussed in Appendix C.

The Bayesian model evidences are computed with nested sampling using MultiNest Feroz and Hobson (2008); Feroz et al. (2009, 2019a). The prior of source sky position $\hat{\Omega}=(\alpha$, $\delta)$ is taken to be uniform over the sky sphere. The prior of polarization angle $\psi$ is taken to be uniform over $[0,\pi]$. The prior of each relative amplitude in $\boldsymbol{\mathcal{A}}$ is taken to be uniform over $[0,2]$. The prior of each relative phase in $\boldsymbol{\varphi}$ is taken to be uniform over $[0,2\pi]$. 1024 live points are used. The sampling efficiency is set to be $0.3$ as recommended in the GitHub repository of MultiNest Feroz et al. (2019b).