Probing Dynamical Gravity with the Polarization of Continuous Gravitational Waves

TL;DR

This paper introduces new methods to detect and analyze the polarization content of continuous gravitational waves, enabling tests of general relativity and its alternatives using data from advanced detectors.

Contribution

It presents a novel framework for identifying and constraining polarization modes in continuous gravitational wave signals, independent of specific models.

Findings

Detectability threshold for polarization signals at strain $h\gtrsim 3\times10^{-27}$

Detection of non-tensorial modes depends on their power, not tensorial strain strength

Framework tested successfully on simulated data from advanced detectors

Abstract

The direct detection of gravitational waves provides the opportunity to measure fundamental aspects of gravity which have never been directly probed before, including the polarization of gravitational waves. In the context of searches for continuous waves from known pulsars, we present novel methods to detect signals of any polarization content, measure the modes present and place upper-limits on the amplitude of non-tensorial components. This will allow us to obtain new model-independent, dynamical constraints on deviations from general relativity. We test this framework on multiple potential sources using simulated data from three advanced-era detectors at design sensitivity. We find that signals of any polarization will become detectable and distinguishable for characteristic strains $h\gtrsim 3\times10^{-27} \sqrt{1~{\rm yr}/T}$, for an observation time $T$. We also find that our ability to detect non-tensorial components depends only on the power present in those modes, irrespective of the strength of the tensorial strain.

Equations

$$F_{+}=\frac{1}{2}\left[({\mathbf{w}}_x\cdot {\mathbf{d}}_x{)}^2-({\mathbf{w}}_x\cdot {\mathbf{d}}_y{)}^2-({\mathbf{w}}_y\cdot {\mathbf{d}}_x{)}^2+({\mathbf{w}}_y\cdot {\mathbf{d}}_y{)}^2\right],$$

$$F_{\times }=({\mathbf{w}}_x\cdot {\mathbf{d}}_x)({\mathbf{w}}_y\cdot {\mathbf{d}}_x)-({\mathbf{w}}_x\cdot {\mathbf{d}}_y)({\mathbf{w}}_y\cdot {\mathbf{d}}_y),$$

$$F_{\mathrm{x}}=({\mathbf{w}}_x\cdot {\mathbf{d}}_x)({\mathbf{w}}_z\cdot {\mathbf{d}}_x)-({\mathbf{w}}_x\cdot {\mathbf{d}}_y)({\mathbf{w}}_z\cdot {\mathbf{d}}_y),$$

$$F_{\mathrm{y}}=({\mathbf{w}}_y\cdot {\mathbf{d}}_x)({\mathbf{w}}_z\cdot {\mathbf{d}}_x)-({\mathbf{w}}_y\cdot {\mathbf{d}}_y)({\mathbf{w}}_z\cdot {\mathbf{d}}_y),$$

$$F_{\mathrm{b}}=\frac{1}{2}\left[({\mathbf{w}}_x\cdot {\mathbf{d}}_x{)}^2-({\mathbf{w}}_x\cdot {\mathbf{d}}_y{)}^2+({\mathbf{w}}_y\cdot {\mathbf{d}}_x{)}^2-({\mathbf{w}}_y\cdot {\mathbf{d}}_y{)}^2\right],$$

$$F_{\mathrm{l}}=\frac{1}{2}\left[({\mathbf{w}}_z\cdot {\mathbf{d}}_x{)}^2-({\mathbf{w}}_z\cdot {\mathbf{d}}_y{)}^2\right].$$

$$F_{+}(\vartheta ,\phi ;\psi =0)=-\frac{1}{2}\left(1+{\cos }^2\vartheta \right)\cos 2\phi ,$$

$$F_{\times }(\vartheta ,\phi ;\psi =0)=-\cos \vartheta \sin 2\phi ,$$

$$F_{\mathrm{x}}(\vartheta ,\phi ;\psi =0)=-\sin \vartheta \sin 2\phi ,$$

$$F_{\mathrm{y}}(\vartheta ,\phi ;\psi =0)=\sin \vartheta \cos \vartheta \cos 2\phi ,$$

$$F_{b/l}(\vartheta ,\phi ;\psi =0)=\mp \frac{1}{2}{\sin }^2\vartheta \cos 2\phi ,$$

$$h(t)=\sum _p{F}_p(t)h_p(t),$$

$$h_p(t)=a_p\cos \left(\varphi (t)+{\varphi }_p\right),$$

$$h_{+}(t)=h_0\frac{1}{2}(1+{\cos }^2\iota )\cos \varphi (t),$$

$$h_{\times }(t)=h_0\cos \iota \sin \varphi (t),$$

$$h_{\mathrm{x}}=h_{\mathrm{y}}=h_{\mathrm{s}}=h_{\mathrm{l}}=0,$$

$$h_0=\frac{16{\pi }^{2}G}{c^4}\frac{ϵI_{zz}{f}_{\mathrm{rot}}^2}{r},$$

$$\varphi (t)=2{\varphi }_{\mathrm{rot}}(t)+{\varphi }_{(\mathrm{GW}-\mathrm{EM})},$$

$$h(t)=\Lambda (t)e^{i\varphi (t)}+{\Lambda }^{\ast }(t)e^{-i\varphi (t)},$$

$$\Lambda (t)\equiv \frac{1}{2}\sum _{p=1}^5{a}_p{e}^{i{\varphi }_p}{F}_p(t_k;\psi ,\alpha ,\delta ),$$

$$\varphi (t)=2\pi \left(2f_{\mathrm{rot}}\tau +{\dot{f}}_{\mathrm{rot}}{\tau }^2+...\right),$$

$$\tau (t)=t+{\Delta }_{\mathrm{R}}+{\Delta }_{\mathrm{E}}+{\Delta }_{\mathrm{S}}+{\Delta }_{\mathrm{binary}}.$$

$$h^{\prime }(t)\equiv h(t)e^{-i\varphi (t)}=\Lambda (t)+{\Lambda }^{\ast }(t)e^{-i2\varphi (t)}$$

$$B_{\mathrm{expected}}(t_k)=\Lambda (t_k)+n(t_k),$$

$$p(B_k\mid \sigma ,{\mathcal{H}}_{\mathrm{N}})=\frac{1}{2\pi {\sigma }^2}\exp \left(-\frac{|B_k{|}^2}{2{\sigma }^2}\right).$$

$$P(\mathbf{B}\mid {\mathcal{H}}_{\mathrm{N}})=\prod _{j=1}^{N_{\mathrm{S}}}{A}_j{\left(\sum _{k={\kappa }_j}^{K_j}{\left|B_k\right|}^2\right)}^{-s_j},$$

$$P(\mathbf{B}\mid {\mathcal{H}}_{\mathrm{N}})=\prod _{i=1}^{N_{\mathrm{D}}}\prod _{j=1}^{{N_{\mathrm{S}}}_i}{A}_{i,j}{\left(\sum _{k={\kappa }_{i,j}}^{K_{i,j}}{\left|B_{i,k}\right|}^2\right)}^{-s_{i,j}},$$

$${\Lambda }_{\mathrm{t}}(t)=\frac{1}{2}\left[a_{+}{e}^{i{\varphi }_{+}}{F}_{+}(t;\psi =0)+a_{\times }{e}^{i{\varphi }_{\times }}{F}_{\times }(t;\psi =0)\right],$$

$${\Lambda }_{\mathrm{GR}}(t)=\frac{1}{2}{h}_0{e}^{i{\varphi }_0}$$

$${\Lambda }_{\mathrm{st}}(t)=\frac{1}{2}$$

Full text

††preprint: LIGO-P1600305

I Introduction

The recent detection of gravitational waves (GWs) by the advanced Laser Interferometer Gravitational-Wave Observatory (aLIGO) heralds the beginning of the long-awaited era of GW astronomy Abbott et al. et al. (2016a, b). One of the main goals of this field is to use GWs as a probe of fundamental physics in the highly dynamical and strong-field regimes of gravity, as predicted by the general theory of relativity (GR). The first few GW detections have already been used to place some of the most stringent constraints on deviations from GR in this domain, which is inaccessible to laboratory, Solar System or cosmological tests of gravity.

However, it has not been possible to use LIGO signals to learn about the polarization content of GWs Abbott et al. et al. (2016c), a measurement highly relevant when comparing GR to many of its alternatives Will (1993, 2014). The reason for this is that the relative orientation of the two LIGO detectors makes it nearly impossible to unequivocally characterize the polarizations of transient GW signals like the compact-binary coalescences (CBCs) observed so far. In fact, at least five noncoaligned quadrupolar detectors would be needed to break the degeneracies of all five nondegenerate polarizations allowed by generic metric theories of gravity Eardley et al. (1973a, b).

Existing observations that are usually taken to constrain the amount of allowed non-GR polarizations can do so only in an indirect and strongly model-dependent manner. For example, measurements of the orbital decay of binary systems are sensitive to the total radiated GW power, but do not probe the waves directly (see e.g. Weisberg et al. (2010); Freire et al. (2012), or Stairs (2003); Wex (2014) for reviews). In the context of specific alternative theories (e.g. scalar-tensor) such observations can indeed constrain the power in extra polarizations; however, they provide no direct, model-independent information on the actual polarization content of the gravitational radiation. Thus, there may be multiple theories, with different polarization content, that still predict the correct observed GW emitted power. Because other traditional tests of GR (like Solar System tests) have no bearing on GWs, there currently exist no direct measurements of GW polarizations.

Unlike CBC transients, continuous gravitational waves (CWs) are, by definition, long-lasting narrow-band signals. Although they have not yet been observed Abbott et al. (2016a); Aasi et al. (2015); Abbott et al. (2016b); Aasi et al. (2014, 2016); Abbott et al. et al. (2017), CWs are expected to be emitted by stable systems, like spinning neutron stars with an asymmetric moment of inertia Thorne (1987). If detected, such signals would allow for tests of gravity complementary to those achievable with transients, including the study of GW polarizations Isi et al. (2015).

In Isi et al. (2015) we showed that it is possible to search for CWs in a polarization-agnostic way and to disentangle the polarization content if a signal is present. However, the data analysis methods proposed were based on a frequentist approach to statistics and suffered from the associated limitations. In this paper, we reframe the ideas of Isi et al. (2015) in a more sophisticated Bayesian framework that allows us to achieve the following novel goals:

Model-independent detection: determine whether a set of GW detector data, prepared for any given known pulsar and from one or multiple detectors, provides evidence for the presence of an astrophysical signal of any polarization content. 2. 2.

Model selection: in the presence of a signal, determine whether the data favor GR or a generic non-GR model, as well as comparing specific alternative theories among themselves and to GR; combine data for multiple sources into a single statement about the validity of GR. 3. 3.

Inference: if the data favor the presence of a GR signal, place constraints on specific alternative theories using the tools of Bayesian parameter estimation.

Furthermore, while Isi et al. (2015) treated only the case of a single detector, we are now able to consider the generic case of a network of detectors.

We present Bayesian methods to achieve the three goals above in the context of searches targeted to known pulsars and present sensitivity estimates for the advanced detector era, including the first generic estimates of sensitivity to nontensorial CW polarizations ever published. In Sec. II, we review the basics of beyond-Einstein polarizations and the targeted pulsar CW search. In Sec. III, we phrase our problem in the language of model selection and explain the construction of hypotheses that will allow us to distinguish GR from non-GR signals. In Sec. IV we specify the details of our analysis, and we explain our results in Sec. V. Finally, we summarize our findings and explain caveats in section VI.

II Background

II.1 Polarizations

GWs can be decomposed into different polarizations, which arise from the linearly independent components of the three-dimensional, rank-two tensor representing the spatial metric perturbation Eardley et al. (1973b). A generic metric theory of gravity may thus allow any combination of up to six independent modes: plus ($+$), cross ($\times$), vector x (x), vector y (y), breathing (b) and longitudinal (l). The effect of each of these modes is represented in Fig. 1. The rotational properties of the fields underlying any given theory determines which polarizations the theory supports: $+$ and $\times$ correspond to tensor fields (helicity $\pm 2$), x and y to vector fields (helicity $\pm 1$), and b and l to scalar fields (helicity [math]).

The components of the tensor and vector pairs are not separable, in the sense that a signal model that includes one element of the group must also include the other (e.g. it is not possible to have a model that allows plus $+$ but not $\times$), because the distinction between $+$ and $\times$, or x and y, is contingent on the frame of reference (e.g. relative orientation of source and detector).

Einstein’s theory only allows the existence of the $+$ and $\times$ polarizations. On the other hand, scalar-tensor and massive-graviton theories may also predict the presence of some b and/or l component associated to the theory’s extra scalar field Will (2014). On top of tensor and scalar modes, bimetric theories, like Rosen or Lightman-Lee theories, also predict vector modes Chatziioannou et al. (2012). Furthermore, less conventional theories might, in principle, predict the existence of vector or scalar modes only, while still possibly being in agreement with all other non-GW tests of GR (see Mead (2015) for an example). Although all these different theoretical frameworks serve as motivation for our study, our approach to the measurement of GW polarizations is phenomenological and, thus, theory-agnostic (Sec. III). It is important to underscore that the detection of a GW signal with a non-GR polarization, no matter how small, is sufficient to falsify GR (note the converse is not true, however).

Because different polarizations have geometrically distinct effects, GW detectors will react differently to each mode. This is manifested in the detector response function $F_{p}$ for each polarization $p$, which encodes the effect of a linearly $p$-polarized GW with unit amplitude, $h_{p}=1$. Ground-based GW detectors, like LIGO and Virgo are quadrupolar antennas that perform low-noise measurements of the strain associated with the differential motion of two orthogonal arms. Their detector response function can thus be written as Nishizawa et al. (2009); Błaut (2012); Isi et al. (2015); Poisson and Will (2014):

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Here, the spatial vectors ${\bf d}_{x}$, ${\bf d}_{y}$ have unit norm and point along the detector arms such that ${\bf d}_{z}={\bf d}_{x}\times{\bf d}_{y}$ is the local zenith; the direction of propagation of the wave from a source at known sky location (specified by right ascension $\alpha$, and declination $\delta$) is given by ${\bf w}_{z}$, and ${\bf w}_{x}$, ${\bf w}_{y}$ are such that ${\bf w}_{z}={\bf w}_{x}\times{\bf w}_{y}$. We choose ${\bf w}_{x}$ to lie along the intersection of the equatorial plane of the source with the plane of the sky, and let the angle between ${\bf w}_{y}$ and the celestial north be $\psi$, the polarization angle.

Because of their symmetries, the breathing and longitudinal modes are fully degenerate to networks of quadrupolar antennas (see e.g. Sec. VI of Chatziioannou et al. (2012)). This means that no model-independent measurement with such a network can possibly distinguish between the two, so it is enough for us to consider just one of them explicitly. We will refer to the scalar modes jointly by the subscript “s”.

The response of gravitational detectors to signals of a given polarization and direction of propagation can be represented, as in Fig. 2, by a spherical polar plot in which the radial coordinate corresponds to the sensitivity given by Eqs. (1–6). In the frame of a given detector, this can be written as [see e.g. Eqs. (13.98) in Poisson and Will (2014) with $\psi=-\pi/2$, to account for the different wave-frame definition]:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $\vartheta$ and $\varphi$ are the polar an azimuthal coordinates of the source with respect to the antenna at any given time (with detector arms along the $x$ and $y$-axes), and we have fixed the wave frame so that $\psi=0$. The representation of Fig. 2 makes it clear that quadrupolar detectors will generally be more sensitive to some polarizations than others, although this will vary with the sky location of the source. For example, for all but a few sky locations, quadrupolar antennas will respond significantly less to a breathing signal than a plus or cross signal.

For a given detector, polarization angle and sky location, the antenna patterns of Eqs. (1–6) become simple, distinct functions of time determined by the rotation of the Earth. This can be pictured by noting that, as the Earth spins on its axis, the angular location of the source with respect to detector will change, tracing an arc on the surfaces of Fig. 2 with varying radial distance. As we explain in Sec. III.1.1, the $F_{p}$’s of polarizations with different rotational properties can be distinguished even in the absence of information on the source orientation; for the minority of cases in which such information exists, it can be taken into account to better distinguish among specific signal models (see Sec. IV).

Because their characteristic period (a sidereal day) is much longer than the CBC timescale (order of minutes or less), the $F_{p}$’s are treated as constants in transient searches; however, this simplification is not valid for CW searches, since their coherent-integration time can be of the order of months to years. As we have pointed out before, this can be used to distinguish the polarization content of a signal Isi et al. (2015). Assuming wave frequency and speed are the same for all modes, the only differences between CWs of different polarizations arise from the sidereal-period amplitude modulations caused by each antenna pattern.

II.2 Continuous waves

II.2.1 Signal

A CW is an almost-monochromatic gravitational perturbation with constant intrinsic amplitude and phase evolution $\phi(t)$. For arbitrary polarization content, such a GW will induce a strain in a quadrupolar detector which can be written as:

[TABLE]

where the sum is over the five independent polarizations, $p\in\{+,~{}\times,~{}{\rm x},~{}{\rm y},~{}{\rm s}\}$; the $F_{p}$’s are those of Eqs. (1–5), and thus implicitly depend on the relative location and orientation of source and detector by means of $\psi$, $\alpha$ and $\delta$; the $h_{p}$ term encodes the amplitude and phase of the wave before being projected onto the frame of the detector:

[TABLE]

where $a_{p}$ is a time-independent amplitude with a functional dependence on source parameters determined by each particular theory of gravity; $\phi(t)$ the phase evolution, a consequence of the dynamics of the source in that theory; and $\phi_{p}$ a phase offset for each polarization. The polarization amplitudes $a_{p}$ and phases $\phi_{p}$ may take arbitrary values depending on the specific theory of gravity and emission mechanism.

In GR, there are several ways in which a neutron star could emit CWs, but the most likely is the presence of a nonaxisymmetry in the star’s moment of inertia Jones and Andersson (2002). For this type of triaxial, nonprecessing source, GR predicts:

[TABLE]

[TABLE]

[TABLE]

where $\iota$ is the inclination angle between the spin axis of the source and the observing line-of-sight, and $h_{0}$ is an overall amplitude given by:

[TABLE]

where $r$ is the source distance, $f_{\rm rot}$ its rotation frequency around the principal axis $z$, $\mathbf{I}$ the moment-of-inertia tensor and $\epsilon\equiv(I_{xx}-I_{yy})/I_{zz}$ the equatorial ellipticity. For the triaxial case, the GW frequency $f$ is twice the rotational value $f_{\rm rot}$, so that we can write:

[TABLE]

where $\phi_{\rm rot}$ is the rotational phase as measured via electromagnetic (EM) observations and $\phi_{({\rm GW}-{\rm EM})}$ is a potential, constant phase offset between the GW and EM signals that can be absorbed into the definition of the $\phi_{p}$’s in Eq. (13).

Note that other emission mechanisms may result in GW radiation at $f=f_{\rm rot}$ Zimmermann and Szedenits (1979), or even noninteger powers $f_{\rm rot}$ Owen et al. (1998); Bondarescu et al. (2009); Jones (2015). Furthermore, alternative theories of gravity may (and, in general, will) support signals at any harmonic. Although in this paper we only consider the case in which only the second rotational harmonic appears in the GW phase, the analysis can be easily generalized to also include contributions from the fundamental and other multiples of $f_{\rm rot}$ (see Sec. VI).

II.2.2 Targeted search

We would like to search a given set of data (from one or more detectors) for CW signals coming from a specific candidate pulsar which has already been observed and timed electromagnetically. Timing solutions are obtained through the pulsar timing package TEMPO2 Edwards et al. (2006). We want to achieve this regardless of polarization content, and to reliably distinguish between the different modes present.

If we assume all polarizations share the same phase evolution, then detector response is the only factor distinguishing CW polarizations and, thus, all the relevant information is encoded in the sidereal-day-period amplitude modulation of the signal. This allows us to focus on a narrow frequency band around the expected GW frequency by processing the data following the complex-heterodyne method presented in Niebauer et al. (1993) and Dupuis and Woan (2005). This procedure is summarized below.

A signal like Eqs. (12, 13) can be rewritten in the form:

[TABLE]

[TABLE]

with $*$ indicating complex conjugation and $\phi(t)$ given by a Taylor expansion around $f=2f_{\rm rot}$:

[TABLE]

where $\tau$ is itself a function of time given by:

[TABLE]

In the above, $\tau$ is the time measured by a clock inertial with respect to the pulsar; $t$ is the time as measured at a given detector; $\Delta_{\rm R}$ is the Roemer delay; $\Delta_{\rm E}$ is the Solar-System Einstein delay; $\Delta_{\rm S}$ is the Solar-System Shapiro delay; $\Delta_{\rm binary}$ is the delay originating from the motion of the pulsar in its binary (a term that vanishes for isolated sources) Dupuis and Woan (2005).

It is important to remember that, the $F_{p}$’s are functions of the source orientation and sky location relative to the detector, so we have made this dependence explicit in Eq. (20) by writing $F_{p}(t_{k})$ as $F_{p}(t_{k};\psi,\alpha,\delta)$. Also, recall that these functions have a characteristic period of a sidereal day (${\sim}10^{-5}$ Hz).

Because the phase evolution $\phi(t)$, including all corrections from Eq. (22), is known (with known uncertainties) from electromagnetic observations, we can digitally heterodyne the data by multiplying by $\exp{\left[-i\phi(t)\right]}$ so that the signal therein becomes:

[TABLE]

and the frequency modulation of the first term is removed, while that of the second term is doubled. A series of low-pass filters can then be used to remove the quickly varying term, which enables the down-sampling of the data by averaging over minute-long time bins. As a result, $\Lambda(t)$ is the only contribution from the original signal left in our data, and hence we can use Eq. (20) as the template for our search. Note that, although we started with real-valued data, after this process the data are now complex.

From Eq. (23) we see that, in the presence of a signal, the heterodyned and down-sampled noisy detector strain data ${B_{k}}$ for the $k$th time bin (which can also be labeled by the Earth-frame GPS time-of-arrival at the detector, $t_{k}$) are expected to be of the form:

[TABLE]

where $n(t_{k})$ is the heterodyned, filtered and downsampled noise in bin $k$, which carries no information about the GW signal. Note then that ${B_{k}}(t_{k})-\Lambda(t_{k})$ should be expected to have the statistical properties of noise, a fact that will be used below in defining likelihoods.

III Method

III.1 Model selection

We use the tools of Bayesian model selection (also known as second-level inference) to determine whether the data contain a signal and, if so, whether that signal agrees with the GR prediction or not. Our procedure is hierarchical and consists of the following stages:

detection: select between signal and noise models; 2. 2.

test of GR: if a signal is present, select between GR and non-GR models; 3. 3.

upper limits: if GR is favored, place upper limits on nontensorial strain amplitudes, in the context of specific alternative polarization models.

This subsection covers only the first two items in this list, since the placement of upper limits belongs in the section on parameter estimation. We treat the case of a single data set in III.1.1 and III.1.2, and we show how to combine results from multiple analyses in III.1.3; we offer some considerations about how to approach the problem of non-Gaussian noise in III.1.4.

III.1.1 Hypotheses

For any given pulsar, we would first like to use reduced (i.e. heterodyned, filtered and downsampled) GW data to decide between the following two logically disjoint hypotheses:

noise (${{\cal H}_{\rm N}}$): no signal, the data are drawn from a Gaussian distribution of zero mean and some (possibly slowly varying) standard deviation; 2. 2.

signal (${{\cal H}_{\rm S}}$): the data contain noise drawn from a Gaussian distribution and a signal with the assumed phase evolution and any polarization content.

In order to perform model selection, we need to translate these hypotheses into the corresponding Bayesian models; this means setting a likelihood function derived from the expected noise properties and picking a multidimensional prior distribution over all parameters. It is important to underscore that a Bayesian model is defined by the choice of these two probability distributions.

For ${{\cal H}_{\rm N}}$, the construction of the likelihood is straightforward. First, let $\sigma$ be the standard deviation of the detector noise at or near the expected GW frequency; then, for each complex-valued data point ${B_{k}}$, Gaussianity implies:

[TABLE]

Here, and throughout this document, a lower-case $p$ is used for probability densities, while an uppercase $P$ is used for discrete probabilities.

If the data are split into ${N_{\rm S}}$ segments of lengths $s_{j}$ ($j=1,\dots,{N_{\rm S}}$) over which the standard deviation $\sigma_{j}$ is assumed to remain constant, we can analytically marginalize over this parameter to obtain a likelihood for the entire data set ${\bf B}$ in the form of a Student’s $t$-distribution Dupuis and Woan (2005); Pitkin et al. (2017):

[TABLE]

with $A_{j}=(s_{j}-1)!/2\pi^{s_{j}}$, $\kappa_{j}=1+\sum_{n=1}^{j}s_{n-1}$, $K_{j}=\kappa_{j}+s_{j}-1$ and $s_{0}=0$. Data streams from ${N_{\rm D}}$ detectors can be analyzed coherently by generalizing this to:

[TABLE]

where $i$ indexes detectors, ${B_{i,k}}\equiv B_{i}(t_{k})$ is the datum corresponding to the $i$th detector at time $t_{k}$, and $A_{i,j}$, $\kappa_{i,j}$ and $K_{i,j}$ are defined analogously to $s_{j}$, $\kappa_{j}$ above. The splitting of the data into segments of constant standard deviation may be achieved with a strategy similar to the Bayesian-blocks algorithm of Scargle (1998), and explained in detail in Pitkin et al. (2017).

Note that the likelihood $p({\bf B}\mid\vec{\theta},{{\cal H}})$ of some hypothesis ${{\cal H}}$, is the probability of observing the data ${\bf B}$ assuming ${{\cal H}}$ is true and given a specific choice of free parameters $\vec{\theta}$ from the model’s parameter space $\Theta$. However, in the case of the noise (“null”) hypothesis, as defined by the Student’s $t$ likelihood above, there are no free parameters. Consequently, $\Theta=\varnothing$ and $p({\bf B}\mid\vec{\theta},{{\cal H}_{\rm N}})=P({\bf B}\mid{{\cal H}_{\rm N}})$.

The case of ${{\cal H}_{\rm S}}$ requires more careful attention. One could be tempted to use Eq. (24) to define a likelihood like Eq. (27) with the substitution $|{B_{k}}|\rightarrow|{B_{k}}-\Lambda_{k}|$, for $\Lambda_{k}\equiv\Lambda(t_{k})$ including all polarizations like in Eq. (20); the priors would reflect uncertainties in measured source parameters and extend over reasonable ranges for $a_{p}$ and $\phi_{p}$. However, for most realistic prior choices, that would correspond to a hypothesis that assigns most of the prior probability to regions of parameter space for which $a_{p}\neq 0$ for all $p$, thus downweighting more conservative models (including GR) that we would like to prioritize. This is simply because the subspace in parameter space corresponding to any of these smaller subhypotheses (which, for example, fix one of the $a_{p}$’s to be zero) has infinitely less volume (i.e. it offers infinitesimally less support) than its complement; hence any practical choice of prior probability density will also assign this subspace infinitely less weight, and so the prior for the corresponding subhypothesis will be vanishingly small.

Formally, the inadequacy of the naive construction of ${{\cal H}_{\rm S}}$ as proposed in the previous paragraph is related to the logical independence of nested hypotheses. We refer to this important point multiple times in the following sections; in particular, we discuss it in the context of odds computations in the text surrounding Eq. (38). We refer readers not familiar with this line of reasoning to a similar discussion in Li et al. (2012), or, more generally, to Ch. 4 in Sivia and Skilling (2006) or Ch. 28 in MacKay (2003).

Instead, we will construct ${{\cal H}_{\rm S}}$ from two logically disjoint component hypotheses:

GR signal (${{\cal H}_{\rm GR}}$ or ${{\cal H}_{\rm t}}$): the data contain Gaussian noise and a tensorial signal with the assumed $\phi(t)$; 2. 2.

non-GR signal (${{\cal H}_{\rm nGR}}$): the data contain Gaussian noise and a signal with non-GR polarization content, but with the assumed $\phi(t)$.

The tensorial hypothesis is embodied most generally by a signal model such that

[TABLE]

where $a_{+}$, $a_{\times}$, $\phi_{+}$ and $\phi_{\times}$ are free parameters, and we pick a specific polarization frame by setting $\psi=0$ (we are allowed to do this because of a degeneracy between $\psi$ and $a_{+}$, $a_{\times}$ explained in Appendix A). An alternative parametrization can be derived from the triaxial emission model of Eqs. (14–16), namely

[TABLE]

where the free parameters are now $h_{0}$, $\phi_{0}$, $\iota$ and $\psi$ [in the notation of Eqs. (13, 20), $\phi_{+}=\phi_{0}$ and $\phi_{\times}=\phi_{0}-\pi/2$]. This is the parametrization used in most traditional GR-only searches (see e.g. Abbott et al. et al. (2017); Pitkin et al. (2017)).

The templates of Eq. (28) and Eq. (29) span the same signal space; therefore, if we pick parameter priors properly related by their Jacobian, the respective hypotheses (${{\cal H}_{\rm t}}$ and ${{\cal H}_{\rm GR}}$) will be logically equivalent (i.e. ${{\cal H}_{\rm t}}\equiv{{\cal H}_{\rm GR}}$). However, we will sometimes want to restrict $\psi$ or $\iota$ in Eq. (29) to incorporate measurements of the source orientation (see Table 3 in Aasi et al. et al. (2014)), and compare those results to the unconstrained model of Eq. (28). In such cases, ${{\cal H}_{\rm t}}$ and ${{\cal H}_{\rm GR}}$ are no longer equivalent: the former corresponds to a free-tensor signal, while the latter now corresponds to a GR triaxial signal for some given source orientation [i.e. a signal with the functional dependence on $\iota$ and $\psi$ of Eq. (29)]. Because of lack of any orientation information, this is a distinction without a difference for most pulsars. (See Appendix A for more details.)

The non-GR hypothesis, ${{\cal H}_{\rm nGR}}$, can itself be seen as a composite hypothesis encompassing all the signal models that depart from GR in some way, i.e. models that include polarizations other than $+$ and $\times$. We denote such subhypotheses with a subscript listing the polarizations included in the signal. For example, “st” (meaning “scalar plus tensor”) corresponds to a model with unrestricted scalar and tensor contributions:

[TABLE]

With this notation extended to the names of the relevant hypotheses, we may then write ${{\cal H}_{\rm nGR}}$ as the logical union (“or” junction, $\lor$)

[TABLE]

where, for convenience, we have defined the non-GR subscript set $\tilde{M}$:

[TABLE]

Just as before, we may equivalently use the triaxial parametrization, Eq. (29), for the tensor modes in the non-GR hypotheses by instead defining $\tilde{M}$ as

[TABLE]

where, for example, GR+s denotes a signal template like

[TABLE]

and similarly for GR+v and GR+sv, with the added vector modes. Again, the two definitions of $\tilde{M}$, Eqs. (32, 33), are equivalent unless orientation information is incorporated in the way explained above.

By the same token, the signal hypothesis can be built from the logical union of ${{\cal H}_{\rm GR}}$ or ${{\cal H}_{\rm t}}$, and ${{\cal H}_{\rm nGR}}$:

[TABLE]

with $M$ defined similarly to $\tilde{M}$, but also including the tensor-only hypothesis, ${{\cal H}_{\rm GR}}$ or ${{\cal H}_{\rm t}}$:

[TABLE]

The validity of Eqs. (III.1.1, 35) is contingent on the mutual logical independence of all the ${\cal H}_{m}$’s. This requirement is satisfied by construction, since each of the ${\cal H}_{m}$’s is defined to exclude regions of parameter space that would correspond to other hypotheses nested within it (e.g. ${{\cal H}_{\rm GR+s}}$ is defined over all values of the scalar amplitude except $a_{\rm s}=0$, to avoid including ${{\cal H}_{\rm GR}}$). In practice, however, it is not necessary to explicitly exclude these infinitesimal regions of parameter space, as will be explained in the following section.

III.1.2 Odds

We can construct a Bayesian model for ${{\cal H}_{\rm S}}$ starting from its components: for each subhypothesis ${\cal H}_{m}$ for $m\in M$, we use a likelihood function like Eq. (27) with the substitution $|{B_{i,k}}|\rightarrow|{B_{i,k}}-\Lambda_{m,i,k}|$, i.e.

[TABLE]

(where $\Lambda_{m,i,k}$ is the template corresponding to model $m$, for detector $i$ and time-bin $k$), and suitable priors on the model parameters $\vec{\theta}_{m}\in\Theta_{m}$; then, we combine the posteriors with priors on the models themselves to obtain the posterior for ${{\cal H}_{\rm S}}$. This last step allows us to incorporate our a priori beliefs about the validity of each of the components. This procedure is represented schematically in Fig. 3 and fleshed out below.

The choice of model priors can be made clearer by considering the posterior probability for the signal model. Given some set of detector data ${\bf B}$ and underlying assumptions $I$ (suppressed from the following expressions), the posterior probability for ${{\cal H}_{\rm S}}$ is:

[TABLE]

by Eq. (35) and because the components are all logically independent [i.e. ${\cal H}_{m_{1}}\land{\cal H}_{m_{2}}={\rm False}$, hence $P({\cal H}_{m_{1}}\land{\cal H}_{m_{2}}\mid{\bf B})=0$ for any $m_{1},~{}m_{2}\in M$ such that $m_{1}\neq m_{2}$]. Note that this is true even for hypotheses that may contain each other as special cases. For instance, even though the GR template can be obtained from GR+s by setting the scalar amplitude to $a_{\rm s}=0$, the points in the GR+s parameter space satisfying this condition define an infinitesimally-thin slice in parameter space that offers no support to the prior distribution and is thus assigned no weight (see similar discussion in Li et al. (2012)).

We can expand each term on the RHS of Eq. (38) using Bayes’ theorem:

[TABLE]

Each of the terms on the RHS is made up of three factors: a marginalized likelihood $P({\bf B}\mid{\cal H}_{m})$, a prior on the model $P({\cal H}_{m})$, and a normalization constant $P({\bf B})$.

The marginalized likelihood (also known as evidence) is computed from the data:

[TABLE]

where $p({\bf B}\mid\vec{\theta}_{m},{\cal H}_{m})$ is itself the likelihood. The evaluation of the multidimensional integral of Eq. (40) is the most computationally intensive part of our analysis (see Sec. IV for details).

We are free to choose the model priors (discussed in Sec. IV), as long as we satisfy the constraint:

[TABLE]

This is a statement about the exhaustiveness and disjointedness of the hypotheses we are considering: we assume that reality will agree with one and only one of the hypotheses at hand. (As we will see in Sec. VI, this assumption might not hold; for example, the noise may not be Gaussian.) The particular choice of prior for each model will encode our expectations about the corresponding theory (before seeing the data), and thus allow for some degree of subjectivity.

Note that we cannot directly compute $P({\bf B})$ in a straightforward manner and without assuming that our hypothesis set is indeed exhaustive (which is not the case for non-Gaussian detector noise, see Sec. III.1.4). However, the need for this computation can be avoided by looking at relative probabilities, i.e. odds. The odds for ${{\cal H}_{\rm S}}$ versus ${{\cal H}_{\rm N}}$ is defined as:

[TABLE]

Using Bayes’ theorem again and canceling the $P({\bf B})$ factors, this simplifies to:

[TABLE]

where, in the second equality, we have used the definition of the Bayes factor:

[TABLE]

for any two hypotheses ${\cal H}_{i}$, ${\cal H}_{j}$.

The odds in Eq. (43) can be used as a detection statistic to determine whether it is likely that the data contain a signal (of any polarization) or not. Once the presence of a signal has been established, a similar ratio can be constructed to assess agreement with GR:

[TABLE]

This ratio encodes the relative probability that there is a GR violation. Because it is now assumed that there is a signal in the data, $P({{\cal H}_{\rm N}})=0$ and the model priors must instead satisfy:

[TABLE]

We can reduce the number of computations needed to obtain ${\cal O}^{\rm S}_{\rm N}$ and ${\cal O}^{\rm nGR}_{\rm GR}$ by using the fact that:

[TABLE]

This means that we need to evaluate an integral like Eq. (40) seven times per set of data, to compute ${\cal B}^{m}_{\rm N}$ for each $m$ in $M$. Those seven numbers, together with the evidence for ${{\cal H}_{\rm N}}$, are enough to compute all the quantities of interest.

Instead of asking about a generic deviation from GR, we may also compare GR to a particular alternative theory. For such purpose, we will usually assign equal prior weight to GR and its alternative to compute:

[TABLE]

where ${\cal H}_{j}$ may be any of the hypotheses in $\tilde{M}$ or an even more specific hypothesis. (The latter case demands an extra execution of the inference code.)

III.1.3 Multiple data sets

So far we have assumed that the data ${\bf B}$, corresponding to one or more GW detectors, can be analyzed coherently; however, there are cases in which we would like to combine results from sets of data analyzed incoherently. Examples are data sets corresponding to different sources or observation periods. Our Bayesian framework makes it possible to combine the respective odds in order to make an overall model selection statement (in our case, about the presence of signal or the validity of GR).

For instance, we may analyze data for ${N_{\rm P}}$ pulsars and ask about the probability that any of them contain a signal; treating each as an independent observation, the combined probability can be constructed from the odds above. Letting ${\cal H}_{{\rm S}_{i}},~{}{\cal H}_{{\rm N}_{i}}$ respectively denote signal and noise hypotheses for the $i$th source, while ${{\cal H}_{\rm S_{any}}}$ corresponds to a signal being present in any of the sources and ${{\cal H}_{\rm N_{all}}}$ corresponds to Gaussian noise in data for all sources:

[TABLE]

where we have used the exclusivity and exhaustiveness of the signal and noise hypotheses, i.e.

[TABLE]

[TABLE]

with $i$ indexing data sets. Note that the data sets for different sources (${\bf B}_{i}$’s) are not conditionally independent under ${{\cal H}_{\rm S_{any}}}$ or ${{\cal H}_{\rm N_{all}}}$. Also, Eq. (III.1.3) does not enforce the requirement that, if signals are present in multiple sources, they all correspond to the same model from Eq. (36); such a constraint could be implemented at this stage, but is more easily enforced by examining individual values of ${\cal O}^{m}_{\rm N}$ when necessary.

The construction of Eq. (III.1.3) implicitly assigns model priors to each of the meta-hypotheses ${{\cal H}_{\rm S_{any}}}$ and ${{\cal H}_{\rm N_{all}}}$ such that:

[TABLE]

where we have assumed the priors for signal vs noise are equal for all sources, i.e. $P({{\cal H}_{{\rm S}_{i}}})=P({{\cal H}_{\rm S}})$ and $P({{\cal H}_{{\rm N}_{i}}})=P({{\cal H}_{\rm N}})$ for all $i$. When making combined statements for multiple sources, we may wish to choose $P({{\cal H}_{\rm S}})/P({{\cal H}_{\rm N}})$ such as to produce any desired value of $P({{\cal H}_{\rm S_{any}}})/P({{\cal H}_{\rm N_{all}}})$, say $P({{\cal H}_{\rm S_{any}}})=P({{\cal H}_{\rm N_{all}}})$. Furthermore, one may wish to weight each pulsar differently within ${{\cal H}_{\rm S_{any}}}$ by incorporating information about the source distance (or other parameters) into the priors via a parametrization like Eq. (17); this may improve the sensitivity of the ensemble odds to weak signals in the set, as suggested in Fan et al. (2016). However, using such a parametrization generally implies committing to a specific gravitational theory (or family of theories). We choose not to take such approach in this study.

Besides combining data for multiple pulsars, for a given source, we could also (incoherently) combine the results of analyses using data from different observation periods. Since the astrophysical CWs we are considering should either be present in all ${N_{\rm R}}$ observation runs or in none of them, the relevant odds, generalizing Eq. (43), are:

[TABLE]

where we have again used ${\bf B}=\{{\bf B}_{j}\}_{j=1}^{N_{\rm R}}$ to refer to the totality of data, with $j$ indexing observation runs. The independence of the ${\bf B}_{j}$’s, conditional on ${\cal H}_{m}$ and ${{\cal H}_{\rm N}}$, is applied on the last line to write the result in terms of the individual Bayes factors for each run, $\left({\cal B}^{m}_{{\rm N}}\right)_{j}$.

Similarly, we can use multiple data sets to make a single statement about deviations from GR. Once we have made ${N_{\rm P}}$ detections from different sources, the odds for a GR violation is:

[TABLE]

where, again, $i$ indexes sources; this is a generalization of Eq. (45). (See Sec. IIID of Li et al. (2012) for an analogous derivation.)

III.1.4 Non-Gaussian noise

Up to this point, like most other CW studies, we have assumed that the detector noise is Gaussian. However, although previous work has indicated that this is generally a very good approximation Isi et al. (2015); Abbott et al. et al. (2017), it is not exactly true for actual detector noise (for some frequencies more so than others). Happily, most of the model selection statements expounded so far are valid also in the presence of non-Gaussian instrumental noise, after some light reinterpretation.

If the assumption of Gaussianity does not hold, the hypotheses constructed in Sec. III.1.1 are no longer exhaustive: the data may not only be explained by Gaussian noise or a signal (GR or otherwise), but also by non-Gaussian artifacts that are impossible to satisfactorily model. Nevertheless, the computation and interpretation of evidences and odds remain unchanged for all the hypotheses under consideration.

Because “noise” no longer just means “Gaussian noise”, ${\cal O}^{\rm S}_{\rm N}$ (which compares the signal model vs Gaussian noise) has to be treated more carefully for detection purposes. Indeed, instrumental features that are clearly non-Gaussian (e.g. a loud, narrow-band artifact wandering across the frequency of interest) will generally result in a relatively large value of ${\cal O}^{\rm S}_{\rm N}$, even if there is no detectable astrophysical signal in the data. This issue affects the standard GR searches as well Abbott et al. et al. (2017), although perhaps to a lesser degree due to the reduced signal parameter space.

It is possible to mitigate this problem by constructing a hypothesis that captures some key characteristic of instrumental features and helps discriminate those from real astrophysical signals. Perhaps the best way to do this is to take advantage of the fact that an astrophysical CW must manifest itself coherently across detectors, while the same is not true for detector artifacts Keitel et al. (2014). We can thus define an instrumental feature hypothesis (${{\cal H}_{\rm I}}$) to encompass the cases in which the data are composed of Gaussian noise, or features that look like astrophysical signals but are not coherent across detectors (viz. they do not have a consistent phase evolution and they are best described by different waveform parameters).

Formally, we define ${{\cal H}_{\rm I}}$ by:

[TABLE]

where the subscript $d$ identifies detectors, and $\land$ is the logical “and” junction. This definition does not explicitly encompass instrumental features that are coherent across some subset of the detectors. Also, note that Eq. (55) implicitly contains a term equivalent to the usual noise hypothesis ${{\cal H}_{\rm N}}=\bigwedge_{d}{{\cal H}_{{\rm N}_{d}}}$. Similarly, it also contains a term corresponding to the presence of signals in all detectors ($\bigwedge_{d}{{\cal H}_{\rm S}}_{d}$). Importantly, such an incoherent term is not equivalent to the coherent signal hypothesis ${{\cal H}_{\rm S}}$, as given by the multidetector likelihood of Eq. (37):

[TABLE]

While the evidence integral of Eq. (40) factorizes into single-detector terms for ${{\cal H}_{\rm N}}$ (due to the null parameter space), the same is not true for ${{\cal H}_{\rm S}}$. Furthermore, because it does not demand detector coherence, the RHS of Eq. (56) is associated with a considerably larger parameter space than the LHS. Thus, in the presence of an astrophysical signal, model selection will favor ${{\cal H}_{\rm S}}$ due to its smaller Occam’s penalty. The same is true, of course, when comparing ${{\cal H}_{\rm S}}$ to ${{\cal H}_{\rm I}}$ as a whole.

From Eq. (55), it is straightforward to write the evidence for ${{\cal H}_{\rm I}}$ as

[TABLE]

and use this to construct the odds comparing against ${{\cal H}_{\rm S}}$:

[TABLE]

Here we have used Eq. (III.1.4), together with the fact that $P({{\cal H}_{{\rm S}_{d}}}\mid{{\cal H}_{\rm I}})+P({{\cal H}_{{\rm N}_{d}}}\mid{{\cal H}_{\rm I}})=1$ and $P({\bf B}\mid{{\cal H}_{\rm N}})=\prod_{d}P({\bf B}\mid{{\cal H}_{\rm N}}_{d})$, to write ${\cal O}^{\rm S}_{\rm I}$ as a function of the detector-coherent signal vs noise Bayes factor ${\cal B}^{\rm S}_{\rm N}$, the single-detector signal vs noise Bayes factors ${\cal B}^{{\rm S}_{d}}_{\rm N}$, and model priors $P({{\cal H}_{\rm S}})$, $P({{\cal H}_{\rm I}})$ and $P({{\cal H}_{{\rm S}_{d}}}\mid{{\cal H}_{\rm I}})$.

As usual, we are free to choose the model priors to give more or less weight to different hypotheses. For example, we recover the choice of Abbott et al. et al. (2017) (Appendix A3) by setting $P({{\cal H}_{{\rm S}_{d}}}\mid{{\cal H}_{\rm I}})=0.5$ for all $d$ and $P({{\cal H}_{\rm S}})=P({{\cal H}_{\rm I}})\times 0.5^{N_{\rm D}}$ such that:

[TABLE]

(When comparing to Appendix A3 of Abbott et al. et al. (2017), however, note that in that work “I” is used to denote both the background information and the “incoherent-signal-or-noise” hypothesis, which can be identified with our ${{\cal H}_{\rm I}}$.)

There is reason to believe that $\ln{\cal O}^{\rm S}_{\rm I}$, with model priors as in Eq. (59), is quite good at picking out instrumental features, even for data from just two instruments Abbott et al. et al. (2017). (Note that we would expect the discriminatory power of $\ln{\cal O}^{\rm S}_{\rm I}$ to grow with the number of detectors available.) However, at the end of the day, we can never be fully confident that ${{\cal H}_{\rm I}}$ will indeed capture all nonastrophysical disturbances. To address this, we may always treat $\ln{\cal O}^{\rm S}_{\rm N}$ and $\ln{\cal O}^{\rm S}_{\rm I}$ as any generic detection statistic and use estimates of the background distribution to establish significance.

III.2 Parameter estimation

Besides choosing between different models, we can use Bayesian statistics to obtain posterior probability density functions (PDFs) on the parameters of a given template (first–level inference). In the absence of a loud signal, this can be used to obtain credible intervals that yield upper–limits on the amplitudes of GR deviations.

For a model ${{\cal H}}{}$ with $N$ parameters, an $N$-dimensional posterior PDF covering the parameter space $\Theta$ can be obtained from Bayes’ theorem:

[TABLE]

for $\vec{\theta}$ in $\Theta$, and with $p(\vec{\theta}\mid{{\cal H}})$ the prior over $\Theta$. To obtain a one-dimensional PDF for a single parameter (call it $\theta_{i}$), the $N$-dimensional distribution must be marginalized over all nuisance parameters (viz. all parameters except $\theta_{i}$):

[TABLE]

where $0<j\leq N$, such that $j\neq i$, and $\Theta^{\prime}$ denotes the parameter space $\Theta$ with the $i$th dimension removed. Note that the equality has been replaced by a relation of proportionality because we have excluded the evidence $P({\bf B}\mid{{\cal H}})$ from the expression. (Although of great importance for model selection, this quantity is uninteresting for the purposes of parameter estimation and can be treated as a simple normalization constant.) As discussed in Sec. IV, we evaluate Eq. (III.2) with the same algorithm used to compute the evidence.

Eq. (III.2) can be used to place upper limits on model parameters; in particular, we will use it to place limits on the amplitude of GR deviations. Consider, for instance, the case of a scalar-tensor theory that can be encapsulated by our GR+s model as described in the previous section; the $95\%$-credible upper limit on the strength of the breathing mode is $h_{\rm s}^{95\%}$, defined by:

[TABLE]

where $\min(h_{\rm s})$ is the minimum value of $h_{\rm s}$ allowed by the prior.

Note that there may be reasons to compute posteriors under different priors than when computing Bayes factors. In particular, it is conventional to present upper limits obtained using a uniform prior over some broad range of the amplitude parameters. With a uniform prior, the posterior is trivially related to the likelihood. This approach produces a more conservative upper limit than other choices, e.g. a Jeffreys prior (see Appendix B).

IV Analysis

We quantify our ability to use Bayesian model selection to detect CW signals and determine their polarization content as described above. To do this, we use one year of simulated data from three advanced interferometric detectors at design sensitivity: LIGO Hanford (H1), LIGO Livingston (L1) and Virgo (V1). Detector noise is simulated by drawing from a Gaussian distribution with zero mean and variance corresponding to the power spectral density (PSD) of each detector at the GW frequency of the pulsar. (Previous work has shown that these are good assumptions for actual reduced detector data Isi et al. (2015); Abbott et al. et al. (2017).)

As mentioned in the previous section, the key step in our analysis is the computation of the evidence integral of Eq. (40) for the hypotheses under consideration (one noise model, plus seven signal submodels; see Sec. III.1.1). We carry this out using a version of the Bayesian inference code used for the targeted pulsar search in Pitkin et al. (2012, 2017), which we modified to handle signals from theories other than GR. This inference code is itself built on the implementation of Skilling’s nested-sampling algorithm Skilling (2006) in the LALInference package Veitch et al. (2015), part of the LIGO Algorithm Library Suite lal . This is the same inference software used for CBC analyses, including GW150914 Abbott et al. (2016c).

In computing likelihoods, we take source location, frequency and frequency derivatives as known quantities (relevant uncertainties are negligible for this analysis). Unless otherwise stated, priors uniform in the logarithm are used for amplitude parameters ($h_{0}$ or $h_{p}$’s), since these are the least informative priors for scaling coefficients (also known as “Jeffreys priors”) Jaynes (1968); we make the somewhat arbitrary choice of restricting the strain amplitudes to the $10^{-28}$–$10^{-24}$ range (this is of little consequence for model selection, as explained in Appendix B). Flat priors are placed over all phase offsets ($\phi_{0}$ and all the $\phi_{p}$’s).

All plots for the Crab pulsar (PSR J0534+2200) in Sec. V are produced using known values of its orientation parameters, $\cos\iota$ and $\psi$, and with the triaxial parametrization of tensor modes; for other pulsars, however, the free-tensor parametrization is used instead. (See Sec. III.1.1 and Appendix A.)

We follow common practice by adopting the principle of indifference (see e.g. Ch. 5 of Sivia and Skilling (2006)) in assigning equal prior probability to the signal and noise models, i.e. we let

[TABLE]

We must also decide how to split the prior among the different ${\cal H}_{m}$’s when computing ${\cal O}^{\rm S}_{\rm N}$ and ${\cal O}^{\rm nGR}_{\rm GR}$. In the former case we choose to distribute the prior weight uniformly among all signal models, so that:

[TABLE]

with $|M|=7$ the cardinality of $M$ [i.e. the number of signal models that go into the construction of ${{\cal H}_{\rm S}}$, see Eq. (36)]. In the latter, however, we prioritize GR by setting:

[TABLE]

[TABLE]

This distribution is illustrated schematically in Fig. 4. Note that these are not the only justifiable options; for example, we might want to prioritize ${{\cal H}_{\rm GR}}$ when constructing ${{\cal H}_{\rm S}}$ in order to better handle a noise background that does not conform to our assumption of Gaussianity. (Other strategies to tackle non-Gaussian noise are discussed in Sec. III.1.4.) In any case, the code is sufficiently flexible to make different choices for the model priors if desired.

To study our method in the presence of signal, we perform several injections of scalar, vector and tensor polarizations (and combinations thereof) for all the 200 pulsars analyzed in Abbott et al. et al. (2017). The simulated signals have a range of signal-to-noise ratios (SNRs), which we proxy below by their effective strain amplitudes. We define these in terms of the $a_{p}$’s from Eq. (13) by:

[TABLE]

[TABLE]

[TABLE]

for tensor, vector and scalar signals respectively. Each simulated signal is generated with a random value of the nuisance phase parameters ($\phi_{0}$ or $\phi_{p}$’s). GR injections are always carried out using the triaxial template of Eq. (29), with random orientation parameters ($\psi$ and $\iota$) when those are not known. Location is always taken to be fixed at the known value for each pulsar.

V Results

V.1 Model selection

V.1.1 Signal versus noise

We first show that ${\cal O}^{\rm S}_{\rm N}$, as defined in Eq. (43), can be used to discriminate signals of any polarization from Gaussian noise, without significant loss of sensitivity to GR signals. The black histogram in Fig. 5 shows the distribution of the natural logarithm of this quantity (henceforth, “log-odds”), obtained from the analysis of an ensemble of noise instantiations corresponding to a single source—in this case, the Crab pulsar. For comparison, the gray (unhatched) histogram in Fig. 5 is the analogous distribution for $\ln{\cal B}^{\rm GR}_{\rm N}$ [note that ${\cal B}^{\rm GR}_{\rm N}={\cal O}^{\rm GR}_{\rm N}$ if we assign equal priors to the GR and Gaussian noise models, cf. Eq. (43) with $m={\rm GR}$]; this is the value computed in regular, GR-only targeted pulsar searches, although with different signal amplitude priors Abbott et al. et al. (2017). Note that odds carry an intrinsic probabilistic meaning in terms of gambling probabilities, and a background histogram like this is not required to interpret their value (see e.g. Sivia and Skilling (2006)).

For both quantities shown in Fig. 5, a negative value marks a preference for the noise model (${{\cal H}_{\rm N}}$, as defined at the beginning of section III.1.1). However, note that a conservative (as determined by the priors) analysis should not be expected to strongly favor ${{\cal H}_{\rm N}}$, since the presence of a weak signal below the noise threshold cannot be discarded; this explains why the ranges in the plots of Fig. 5 do not extend to more negative values. Generally speaking, the magnitude of the signal prior volume (viz. the volume of parameter space allowed by the signal model, weighted by the prior function) will determine the mean of background distributions like Fig. 5, which will be more negative the greater the signal volume. This is a manifestation of an implicit Occam’s penalty.

The relationship between the Bayes factors for different signal hypotheses vs noise is illustrated in Fig. 6, which shows violin plots representing the noise-ensemble distributions of $\ln{\cal B}^{m}_{\rm N}$ for all models discussed in III.1.1. The values for $m\in\{\rm s,v,sv,GR,GR+s,GR+t,GR+sv\}$ are combined to produce $\ln{\cal O}^{\rm S}_{\rm N}$ in Fig. 5. As explained above, the “GR” label indicates that the tensor modes have been parametrized using the triaxial model of Eq. (29), with orientation parameters fixed at the known values for the Crab pulsar; on the other hand, the “t” label corresponds to the free-tensor template of Eq. (28). We include both parametrizations to demonstrate the effect of assuming a triaxial emission mechanism and restricting the orientation of the source (see also Appendix A).

Interestingly, Fig. 6 reveals the relationship between $\ln{\cal B}^{m}_{\rm N}$ and the number of degrees of freedom (a proxy for the prior volume) of model $m$: models with more degrees of freedom have a greater prior volume and are correspondingly downweighted, resulting in more negative values of $\ln{\cal B}^{m}_{\rm N}$; this is a manifestation of the Occam’s penalty automatically applied by the Bayesian analysis (see e.g. Ch. 28 in MacKay (2003)). We underscore that this feature arises naturally from the computation of the evidence integral, and not from manually downweighting either model a priori.

If the data contain a sufficiently loud signal of any polarization, the evidence for ${{\cal H}_{\rm S}}$ will surpass that for ${{\cal H}_{\rm N}}$, and this can be used to establish a detection. Fig. 7 shows the response of $\ln{\cal O}^{\rm S}_{\rm N}$ and $\ln{\cal B}^{\rm GR}_{\rm N}$ to the presence of GR and non-GR signals. In particular, the second panel in Fig. 7 shows results for injected signals of the vector-only model of Mead (2015), but the behavior would be the same for scalar-only signals. The general features of these plots confirm our expectations that for weak, subthreshold signals, the analysis should not be able to distinguish between the signal and noise models, yielding a Bayes factor close to unity (more precisely, a value of $\ln{\cal O}^{\rm S}_{\rm N}$ consistent with the background distributions of Fig. 5). Note that, in agreement with Fig. 5, the noise baseline for $\ln{\cal O}^{\rm S}_{\rm N}$ lies below that of $\ln{\cal B}^{\rm GR}_{\rm N}$, due to its greater prior volume.

For stronger (detectable) signals, the basic form of our likelihood functions, Eq. (25), leads us to expect $\ln{\cal O}^{\rm S}_{\rm N}$ to scale linearly with the square of the signal-to-noise ratio (SNR):

[TABLE]

where the variance $\sigma^{2}$ proxies the PSD and we let ${\boldsymbol{\Lambda}}_{\rm MP}$ be the time-series vector corresponding to the maximum probability template; for a stationary PSD, this implies $\ln{\cal B}^{m}_{\rm N}\propto h_{\rm inj}^{2}$, as observed in Fig. 7. The spread around the trendline is due to the individual features of each noise instantiation and (much less so) to numerical errors in the computation of the evidence, Eq. (40). For details on numerical uncertainty, see Appendix C.

From the left panel of Fig. 7, we conclude that $\ln{\cal O}^{\rm S}_{\rm N}$ can be as good an indicator of the presence of GR signals as $\ln{\cal B}^{\rm GR}_{\rm N}$ itself; this implies that we may include non-GR polarizations in our search without significantly sacrificing sensitivity to GR signals. However, the power of $\ln{\cal O}^{\rm S}_{\rm N}$ lies in responding also to non-GR signals. As an example of this, the right panel in Fig. 7 shows $\ln{\cal O}^{\rm S}_{\rm N}$ and $\ln{\cal B}^{\rm GR}_{\rm N}$ as a function of the amplitude of a fully non-GR injection. Here, we have chosen to inject a particular model of vector signal developed in Mead (2015), but the results are generic.

Note that, for sufficiently loud signals, ${{\cal H}_{\rm GR}}$ becomes preferable over ${{\cal H}_{\rm N}}$ (hence $\ln{\cal B}^{\rm GR}_{\rm N}>0$), even when the injection model does not match the search; this is because the noise evidence drops faster than GR’s and becomes very small (i.e. the data do not look at all like Gaussian noise, although they do not match the expected GR signal well either). The particular SNR at which this occurs will depend on the overlap between the antenna patterns of the injection and those of GR, and will consequently vary among sources.

For the interesting case of scalar-tensor theories (here, templates composed of GR plus an extra breathing component, and denoted “GR+s”), the behavior is slightly different. This is both because GR+s has an extra amplitude degree of freedom ($a_{\rm s}$) and, as discussed in Sec. III.1.1, because ${{\cal H}_{\rm GR}}$ can be recovered as a special case of ${{\cal H}_{\rm GR+s}}$ (namely, when $a_{\rm s}\rightarrow 0$). In Fig. 8, we present the log-odds of signal versus noise hypotheses as a function of injected GR ($x$-axis) and scalar ($y$-axis) strengths. These plots divide the ${h_{\rm s}}$–${h_{\rm t}}$ plane in roughly two regions where the associated signal model (${{\cal H}_{\rm S}}$, ${{\cal H}_{\rm GR}}$ or ${{\cal H}_{\rm GR+s}}$) is preferred (black) and where it is not (red). The latter corresponds to the area of parameter space associated with subthreshold signals that cannot be detected.

As expected, the best coverage is obtained when analyzing the data using the model matching the injection, GR+s, (rightmost plot) or the all-signal model (leftmost plot). In both these cases, the results improve with either scalar or tensor SNR. In contrast, the GR analysis (center plot) is sensitive to tensor strain, but, as evidenced by the extended red region in the central plot, it misidentifies strong scalar signals as noise. Nevertheless, if the scalar component is larger than ${\sim}5\times 10^{-26}$, the GR analysis will disfavor the noise hypothesis, even for a small tensor component, as in the right panel of Fig. 7; this is the same behavior observed in Fig. 7. In contrast, the any-signal analysis is sensitive to the total power of the injected signal, regardless of polarization.

We have produced distributions of background $\ln{\cal O}^{\rm S}_{\rm N}$, like those of Fig. 5, for all 200 known pulsars in the sensitive band of the three detectors under consideration (same set analyzed in Abbott et al. et al. (2017). In Fig. 9, these are represented by their respective means and one-sided standard deviations as a function of the pulsar’s GW frequency. The frequency dependence is explained by variations in the instrumental noise spectra. This is explained by the fact that, for a particular prior choice, more information is gained from the data if the noise floor is lower: with less noise it is possible to discard the presence of weaker signals, so the value of $\ln{\cal O}^{\rm S}_{\rm N}$ decreases.

V.1.2 GR vs non-GR

In the presence of a signal, ${\cal O}^{\rm nGR}_{\rm GR}$, as defined by Eq. (45), indicates whether there is reason to believe there is a GR violation or not. Because there could always be an unresolvably small departure from GR, we do not expect our analysis (with priors as chosen) to ever strongly favor the GR hypothesis; rather, in the presence of a GR signal we will find that $\ln{\cal O}^{\rm nGR}_{\rm GR}$ remains relatively close to zero, simply meaning that there is no strong evidence for or against non-GR features. This is indeed the behavior observed in the left panel of Fig. 10, where $\ln{\cal O}^{\rm nGR}_{\rm GR}$ is shown to be roughly insensitive to tensor injection amplitude. For values of ${h_{\rm t}}$ below certain threshold (which, in this case, is around $3\times 10^{-27}$), the search does not detect a signal and, consequently, no information is gained for or against ${{\cal H}_{\rm GR}}$, i.e. $\ln{\cal O}^{\rm nGR}_{\rm GR}\sim 0$. The difference between the two populations (below and above threshold) is determined mainly by the choice of amplitude priors.

The behavior of ${\cal O}^{\rm nGR}_{\rm GR}$ is less ambiguous in the presence of a non-GR signal. For instance, if the data contain a detectable signal that completely lacks tensor components, then ${\cal O}^{\rm nGR}_{\rm GR}$ will unequivocally reflect this. This is evidenced by the growth of $\ln{\cal O}^{\rm nGR}_{\rm GR}$ with injected nontensorial SNR in the right panel of Fig. 10. In other words, while the analysis is inconclusive for GR injections because it cannot discard the presence of subthreshold non-GR components hidden by the noise, vector signals are are clearly identified as not conforming to GR. This is a reflection of the fact that, as mentioned in the introduction, any evidence of a nontensorial component is fatal for GR, but absence of non-GR components does not mean Einstein’s theory is necessarily correct.

As might be expected, ${\cal O}^{\rm nGR}_{\rm GR}$ responds to non-GR signals that include a tensor component with a combination of features from both panels of Fig. 10. As an example, the left plot of Fig. 11 shows ${\cal O}^{\rm nGR}_{\rm GR}$ in the presence of GR+s injections, as a function of injected tensor and scalar amplitudes. This plot can be split into three clearly demarcated regions: one in which the signal is not detected (light red, bottom left), one in which the signal is detected and the non-GR model is preferred (black, top), and one which the signal is detected but where the evidence for a deviation from GR is not clear due to the predominance of the tensorial component (darker red, bottom right). The first corresponds to the subthreshold population on either side of Fig. 10, while the second and third correspond to the above-threshold populations on the right and left sides of Fig. 10 respectively; indeed, note that a horizontal slice taken over the red region of the left plot produces a series of points like those in the left panel of Fig. 10. For reference, Fig. 11 also includes the direct comparison of GR+s and GR on the right.

We can make a stronger statement about the agreement of the data with GR by making use of signals from multiple sources, as discussed in Sec. III.1.3. The power of combining multiple signals is illustrated in Fig. 12, where $\ln{\cal O}^{\rm nGR}_{\rm GR}$, as defined in Eq. (54), is plotted vs number of GR (left) and GR+s (right) signals detected. Note that this presumes that, for each source, the presence of a signal has already been established from the value of $\ln{\cal O}^{\rm S}_{\rm N}$. Computing the ensemble $\ln{\cal O}^{\rm nGR}_{\rm GR}$, as done here, is a good way of summarizing the information contained in the data about the relative likelihoods between the two models, but it provides no information not already present in the set of individual single-source odds.

V.2 Parameter estimation

When no conclusive evidence for a CW is found in the data, we are still interested in placing upper limits on the strength of possible signals (up to some credibility), and this is done as explained in Sec. III.2. By the same token, if a signal consistent with GR is detected, we can always place an upper limit on the amplitude of non-GR modes, even if the odds indicate there is no clear sign of a GR violation.

For instance, we can get a quantitative estimate of our sensitivity to scalar modes from a given source by looking at the distribution of $h_{\rm s}^{95\%}$, defined in Eq. (62), computed for a set of noise-only data instantiations. Such distribution for the Crab pulsar is presented in the left panel of Fig. 13. Similarly, the right panel presents estimates for the sensitivity to vector modes coming from the Crab pulsar, assuming a vector-tensor model. In this case, however, the quantity plotted is the upper limit on total, effective vector strain amplitude ${h_{\rm v}}$, Eq. (68). These plots include distributions produced using the same log-uniform prior used to obtain Bayes factors, as well as more conservative ones obtained using uniform amplitude priors (see Appendix B). In either case, the magnitude of $h_{\rm v}^{95\%}$ is comparable to that of $h_{\rm s}^{95\%}$.

Interestingly, our ability to measure scalar and vector amplitudes is unaffected by the presence of other modes. We illustrate this for the Crab pulsar in Fig. 14, which results from analyzing data with GR+s (left) and GR+v (right) injections. There we plot $h_{\rm s}^{95\%}$ as a function of scalar and tensor injection amplitudes on the left, and $h_{\rm v}^{95\%}$ as a function of vector and tensor injection amplitudes on the right. From these plots, one can conclude that $h_{\rm s}^{95\%}$ and $h_{\rm v}^{95\%}$ are sensitive only to the corresponding scalar and vector components, and not by ${h_{\rm t}}$. (It is worth emphasized that the upper limits, $h_{\rm s}^{95\%}$ and $h_{\rm v}^{95\%}$, are well-defined even when the non-GR component is strong enough to be detected, as is the case for the darker-colored regions.)

As shown previously in the literature, the mean of distributions like those of Fig. 13 will scale with $\sqrt{S_{\rm n}(f)/T}$, where $S_{\rm n}(f)$ is the effective PSD of the detector noise at the expected GW frequency $f$, and $T$ is the integration time (cf. Eq. (26) and Fig. 1 of Dupuis and Woan (2005)). Because of this, the mean of this distribution will vary with the source’s expected GW frequency, as shown in Fig. 15. Following convention, these upper limits are computed using uniform amplitude priors, which means they are a factor of a few less stringent than those obtained with a log-uniform prior (see Fig. 13 and Appendix B). Also, for completeness, Fig. 15 also includes the expected tensor upper limits, $h_{\rm t}^{95\%}$. Note that those values are not the same as would be obtained by the standard GR-only search, because that analysis looks at the triaxial $h_{0}$ of Eq. (29), rather than ${h_{\rm t}}$.

In order to compare our sensitivity to the different polarizations, in Fig. 16 we histogram the the tensor and vector upper limits as a ratio of the scalar upper limits—this includes the $h_{\rm t}^{95\%}$ and $h_{\rm v}^{95\%}$ values shown in Fig. 15, as well as the limits on the individual amplitudes from which they are constructed ($h_{\rm+}^{95\%}$, $h_{\rm\times}^{95\%}$, $h_{\rm x}^{95\%}$ and $h_{\rm y}^{95\%}$). The mean of these distributions (vertical dashed lines in Fig. 16) indicate that, for most pulsars, the scalar upper limit is slightly larger in magnitude than those for the $+$, $\times$, x or y modes; this systematic effect is a manifestation of the decreased sensitivity of quadrupolar detectors to scalar waves, which was discussed in Sec. II.1 (see, in particular, Fig. 2). The fact that the difference between $h_{\rm s}^{95\%}$ and $h_{\rm t}^{95\%}$, or $h_{\rm v}^{95\%}$, is less pronounced can be easily be explained as a statistical factor arising from the definitions of ${h_{\rm t}}$ and ${h_{\rm v}}$ as square-roots of sums of squares, Eqs. (67, 68). Both these scalings are discussed in more detail in Appendix D.

VI Conclusion

We have developed a Bayesian framework to detect CW signals from known sources regardless of polarization content, to disentangle the modes present in a given signal, and to constrain the amplitudes of extra polarizations that may be hiding under the noise. We have implemented this as an extension of LIGO’s Bayesian targeted CW search pipeline Pitkin et al. (2017), and thus benefit from the power of the nested sampling algorithm on which it is based.

We have tested our methods on one year of simulated noise for three advanced-era detectors at design sensitivity (H1, L1, V1), and prepared for a set of multiple known sources in their frequency band. This allows us to estimate our future sensitivity to CW polarizations, in this most optimistic case. Under these conditions and for the Crab pulsar in particular, we expect signals of any polarization to become detectable for characteristic strain amplitudes $h\gtrsim 3\times 10^{-27}$ (Figs. 7 and 8); this threshold will vary among sources, due to differences in position (sky location and orientation) and detector PSD at the expected GW frequency (cf. e.g. Fig. 9). Furthermore, the value of this threshold will decrease linearly with the square-root of the observation time Dupuis and Woan (2005).

A signal louder than the detection threshold will allow us to determine whether its polarization content is consistent with GR or not, and the strength of this statement will depend almost exclusively on the power of the non-GR component (Figs. 10 and 11). In other words, from a model-selection standpoint, the non-GR hypothesis will only be unequivocally favored if the total power in non-GR modes is greater than the threshold value, regardless of the strength of the GR modes. However, for signals that do not satisfy this, we may always place upper limits on nontensorial amplitudes and thus constrain deviations from GR; for instance, Fig. 15 presents the most optimistic projections for 95%-credible upper limits for scalar and vector amplitudes of CW signals from all pulsars in our set ($h_{\rm s}^{95\%}\sim 4\times 10^{-27}$ and $h_{\rm v}^{95\%}\sim 3\times 10^{-27}$, in the best case). As far as we are aware, these are the first generic estimates of sensitivity to scalar and vector polarizations ever published 111Note that sensitivity estimates presented in Isi et al. (2015) were restricted to the specific vector-only model of Mead (2015).

From our projected upper limits, we have found that, at design sensitivity, the LIGO-Virgo network will be generally less sensitive to continuous scalar signals than to the individual vector or tensor modes by factors of 0.45–0.7, depending on the location of the source (Fig. 16); this diminished sensitivity to scalar modes stems from the quadrupolar nature of the detector antenna patterns (Fig. 2 and Appendix D). Also, our injection studies indicate that the upper limits on the amplitudes of nontensorial modes will be roughly unaffected by the presence or absence of a tensor signal in the data (Fig. 14).

Although the results presented here made use of simulated Gaussian noise, the procedure is identical for actual detector data. Furthermore, the assumption of Gaussianity has been shown to hold relatively well for real CW data Isi et al. (2015), so the actual sensitivity limits should not be far from those presented here. If the data are strongly non-Gaussian, however, one must be careful in using $\ln{\cal O}^{\rm S}_{\rm N}$ for detection purposes and may instead wish to adopt one of the strategies suggested in Sec. III.1.4.

Another important limitation of our results is that here we only consider CW signals emitted at $f=2f_{\rm rot}$, while it is to be expected that other mechanisms (within GR or not) allow emission at other harmonics, $f=f_{\rm rot}$ in particular. Yet, the only change required to account for this is to modify the template in Eq. (12) to include terms at different harmonics; the ability to do this already exists within our current infrastructure. We also assume that other aspects of the waves, like their speed, remain in agreement with the GR prediction, an assumption that will be relaxed in a future study.

Acknowledgements.

The authors would like to thank Ian Jones and Walter Del Pozzo for carefully reading this manuscript and providing insightful suggestions; we also thank Tjonnie Li and Carver Mead, as well as many colleagues in the LIGO Scientific Collaboration Continuous Waves group, for many useful comments. LIGO was constructed by the California Institute of Technology and Massachusetts Institute of Technology with funding from the National Science Foundation and operates under cooperative agreement PHY-0757058. We are grateful for computational resources provided by Cardiff University, and funded by an STFC grant supporting UK Involvement in the Operation of Advanced LIGO. MP is funded by the STFC under grant number ST/N005422/1. Plots produced using Matplotlib Hunter (2007). This paper carries LIGO Document Number LIGO-P1600305.

Appendix A Tensor models

A conceptual distinction can be drawn between the triaxial GR model and a free-tensor model that includes $+$ and $\times$ but does not restrict their relative amplitudes (denoted “t”). The former has four free parameters (overall amplitude, $h_{0}$; overall phase, $\phi_{0}$; inclination, $\iota$; polarization, $\psi$) and corresponds to a signal template of the form [same as Eq. (29)]:

[TABLE]

This is a reparametrization of the free-tensor model, which also has four parameters (plus amplitude, $a_{+}$; cross amplitude, $a_{\times}$; plus phase, $\phi_{+}$; cross phase, $\phi_{\times}$) and whose template is [same as Eq. (28)]:

[TABLE]

If $\psi$ and $\iota$ are known, it is clear that the two models are different, since ${{\cal H}_{\rm GR}}$ has two free parameters ($h_{0}$, $\phi_{0}$) and ${{\cal H}_{\rm t}}$ has four ($a_{+}$, $a_{\times}$, $\phi_{+}$, $\phi_{\times}$). If the orientation is not fixed, however, the two models span the same signal space. This is because there is a degeneracy between $\psi$ and $a_{+}$, $a_{\times}$ due to the way the antenna patterns transform under changes in $\psi$:

[TABLE]

[TABLE]

with $\psi^{\prime}=\psi+\Delta\psi$. Eqs. (73) and (74) can be derived from Eqs. (1) and (2) respectively, as in Błaut (2012) (or see, e.g., Sec. 9.2.2 in Thorne (1987)). Consequently, changing $\psi\rightarrow\psi^{\prime}$ in Eq. (72) is equivalent to leaving $\psi$ fixed [at, say, $\psi=0$ as in Eq. (72)] while replacing the plus and cross complex amplitudes by:

[TABLE]

[TABLE]

This is illustrated in Fig. 17.

These rotational properties are easily understood by recalling that GW polarizations can be defined in any frame, although a given signal might look more or less simple given the choice of frame. Eqs. (73) and (74) provide the transformation between frames that are coaligned except for a rotation of $\Delta\psi$ around their $z$-axes. Because waveform predictions, e.g. Eq. (71), are made in specific frames, it is important to orient the wave frame appropriately when working with a given theory and emission mechanism. However, if the signal parametrization is independent of any theory, e.g. Eq. (72), one is free to pick any frame (i.e. any $\psi$).

The relationship between the different tensor model parameters is reflected in the posterior probability plots of Fig. 18. For fixed orientation, both the triaxial (a) and free-tensor (b) analyses accurately determine the amplitude and phase of the injected signal. In panel (b), $a_{+}$ and $a_{\times}$ are constrained to lie within a region consistent with $h^{2}_{\rm inj}=a_{+}^{2}+a_{\times}^{2}$ and $a_{+}/a_{\times}=(1+\cos^{2}\iota)/(2\cos\iota)$, for an effective injection amplitude given by $h_{\rm inj}^{2}=h_{0}^{2}(1+\cos^{2}\iota)^{2}/4+h_{0}^{2}\cos^{2}\iota$, as in Eq. (67). When the orientation is allowed to vary, we observe the expected correlations between the recovered triaxial amplitude and the orientation parameters in panel (c); in this case, $\psi$ and $\cos\iota$ will also become correlated, as better shown in Fig. 19.

The degeneracy between $\psi$ and $a_{+}$, $a_{\times}$ is particularly evident in Fig. 18d, where the one-dimensional PDF for $\psi$ shows that this parameter cannot be constrained, even for a loud signal. Furthermore, joint posteriors between $\psi$ and $a_{+}$ & $a_{\times}$ confirm that this is due to the degeneracy from Eqs. (75) and (76), as seen by comparing these two-dimensional PDFs to Fig. 17. Physically, this is a consequence of the fact that we are free to orient the polarization frame as we wish.

Because their signal templates are degenerate when $\psi$ and $\iota$ are allowed to vary, the distinction between ${{\cal H}_{\rm GR}}$ and ${{\cal H}_{\rm t}}$ is not really meaningful for unfixed orientation. This can be seen from the values of $\ln{\cal B}^{\rm t}_{\rm GR}$ in the cases of known and unknown orientations, as in Fig. 20. On the left panel, ${{\cal H}_{\rm GR}}$ is defined with specific values of $\psi$ and $\cos\iota$ that match the injections; on the right, the ${{\cal H}_{\rm GR}}$ priors allow $\psi$ and $\cos\iota$ to range over their full ranges, and the injections are performed with random values of both. When the orientation is fixed, ${{\cal H}_{\rm GR}}$ will always be preferred to ${{\cal H}_{\rm t}}$ for resolvable signals because of its lower Occam’s penalty; however, that is not true for unfixed orientation. Note that, in the strictest sense, the two hypotheses are not logically equivalent unless their parameter priors are related by the Jacobian of the coordinate transformation between the two parametrizations, Eqs. (71, 72); only in that case all regions of signal space are treated equally by ${{\cal H}_{\rm GR}}$ and ${{\cal H}_{\rm t}}$. This explains the variation in values of $\ln{\cal B}^{\rm t}_{\rm GR}$ on the right of Fig. 20.

If one knew the source orientation and one believed that the only viable mechanism for producing CWs at the assumed frequency in GR is the triaxial model embodied by Eq. (71), then one could include the free-tensor hypothesis and all of its derivatives (i.e. t, st, vt, stv) in the non-GR set $\tilde{M}$, on top of $\{\rm s,v,sv,GR+s,GR+v,GR+sv\}$. Doing so would mean treating a tensor-only signal that does not conform to Eq. (71) as evidence of a GR violation, rather than of a different emission mechanism within GR. Given the many simplifications intrinsic to the triaxial model, however, having that much confidence in its validity seems unwarranted; hence we choose to not take that approach.

Appendix B Amplitude priors

Previous CW Bayesian searches targeted to known pulsars have always applied a flat prior on the signal amplitude parameter Abbott et al. et al. (2017). This is because flat priors, if wide enough, cause the posterior to be only determined by the likelihood (up to normalization), yielding more conservative upper limits on the signal strength. Furthermore, unlike with priors uniform in the logarithm of the quantity, upper limits derived with flat priors will generally not depend on the limits set by the prior (again, assuming the range allowed extends from zero amplitude to some large value that does not truncate the likelihood).

Upper limits obtained using log-uniform priors (uniform in the logarithm of the quantity) will, generally, be dependent on the range of the prior, although not strongly. For example, consider a one-dimensional problem on some positive parameter $x$. For simplicity, further assume we have a flat likelihood between $x=0$ and an upper cutoff at $x=x_{\rm max}$; then, $x_{\rm max}$ will necessarily also be an upper bound for the posterior. Because the likelihood is uniform, below the cutoff the posterior will be determined, up to normalization, by the prior only, i.e. for $x<x_{\rm max}$,

[TABLE]

Now consider a log-uniform prior $p(x\mid{{\cal H}})\propto{\rm d}(\log{x})\propto 1/x$, with a lower bound $x_{\rm min}$, such that $0<x_{\rm min}<x_{\rm max}$. Because such prior is uniform in the $\log x$, this implies that the 95%-credible upper limit on $x$ will be given by:

[TABLE]

Since $x_{\rm max}$ is set by the likelihood (by construction), if the prior is changed by rescaling $x_{\rm min}$ by a factor $\alpha$,

[TABLE]

then, for a given set of data, the upper limit becomes $x^{95\%}_{\alpha}$, satisfying:

[TABLE]

Thus, the dependence of the upper limit on the range defined by the log-uniform prior is quite weak, as illustrated in Fig. 21. This explains why upper limits obtained with a log-uniform prior differ only by a factor of a few from those obtained with a flat one, as seen in Fig. 13.

However, the flat priors do not properly represent our ignorance of the scale of the signal amplitude. This problem manifests itself in negative Bayes factors that too quickly favor the noise hypothesis if no loud signal is clearly present, rather than reflecting our expectation that a signal might be hiding under the noise. This can be seen in Fig. 22, where we show the distributions of $\ln{\cal B}^{\rm GR}_{\rm N}$, obtained for several noise-only data instantiations for the Crab pulsar, corresponding to flat and log-uniform priors in the GR amplitude parameter, $h_{0}$; a uniform prior results in lower values of $\ln{\cal B}^{\rm GR}_{\rm N}$ that strongly favor ${{\cal H}_{\rm N}}$. This behavior is not specific to the GR model.

For most of our analysis, we choose to apply priors uniform in the logarithm of all amplitude quantities. However, for the sake of consistency with previous searches and in order to make our limits more conservative, we also present upper limits produced using flat amplitude priors, as shown in Fig. 13.

Appendix C Numerical error

The fractional numerical error in the computation of the natural logarithm of the evidence by nested sampling is usually estimated by:

[TABLE]

where $N_{\rm live}$ is the number of of live points and $H$ is the information gained in the analysis:

[TABLE]

a quantity that is easy to estimate from the output of the nested sampling code Skilling (2006); Keeton (2011).

An example of the actual statistical error as function of SNR is presented in Figs. 23 & 24, where the injected GR signal amplitude serves as proxy for $\rho$ (for fixed PSD). From these plots it becomes apparent that, although the actual error might exceed the estimator of Eq. (81), its absolute magnitude is quite small and should not affect our results. In any case, Eq. (81) indicates that any level of accuracy may be achieved by increasing the number of live points (at the cost of increased computational burden). For more details on the numerical error of the nested sampling algorithm in LALInference, we refer the reader to Sec. IVB of Veitch and Vecchio (2010).

Appendix D Upper-limit ratios

When comparing upper limits for the different modes, as in Fig. 16, two scalings become apparent: first, the $+$, $\times$, x, and y upper limits are, on average, more stringent than those for the scalar polarization by a factor of ${\sim}1.8$; second, the upper limits on ${h_{\rm t}}$ (${h_{\rm t}}$) are a factor of ${\sim}1.3$ larger than those on the individual $+$ and $\times$ (x and y) amplitudes.

The scaling between the scalar upper limit and those for the other individual strain amplitudes can be accounted for by the decreased sensitivity of quadrupolar GW detectors to scalar waves. For a single instrument (that is, not a network), this can be appreciated visually from Fig. 2, by noting that for most sky locations the magnitude of the scalar response is considerably less than for the other modes.

To properly evaluate the effect of the detector geometry on the analysis, however, it is necessary to look at the relative SNRs of unit-amplitude scalar, vector and tensor GWs from a given source, as they are received by the detector network under consideration (H1, L1, V1) after some fixed observation time. Assuming all detectors have comparably noise levels, the network SNR can be proxied by the root-mean-square (RMS) amplitude of the effective network antenna patterns, defined by

[TABLE]

for each polarization $p$, some long observation time $T$, and where the sum is over detectors $d$. [Here we have fixed the source and detector parameters so that the $F_{p}$’s of Eqs. (1–5) are now just simple functions of time.] We may then compute this for all five polarizations and for multiple sources to obtain a sky-average of the ratio of the scalar RMS antenna pattern to those of the other polarizations. We find this ratio to be roughly ${\sim}0.55$ for all polarizations, in agreement with Fig. 16, since we should expect

[TABLE]

where the average $\langle\cdot\rangle$ is taken over multiple sources distributed across the sky.

The relation between the ${h_{\rm t}}$ (${h_{\rm v}}$) upper limits and those for their component amplitudes, $+$ and $\times$ (x and y), can be easily understood by noting that, if using flat priors and in the absence of signal, the marginalized posteriors for each of the component amplitudes ($h_{+}$, $h_{\times}$, $h_{\rm x}$, $h_{\rm y}$) will roughly be described by a one-sided normal distribution. Consequently, it can be shown that posterior for the square-root of the sum of the squares of two of these quantities will be given by a chi distribution with two degrees of freedom. Considering the definitions of Eqs. (67, 68). It is straightforward to show (numerically or analytically) that this explains the observed factor of ${\sim}1.3$ difference between $h_{\rm t}^{95\%}$ ($h_{\rm v}^{95\%}$) and $h_{\rm+}^{95\%}$ or $h_{\rm\times}^{95\%}$ ($h_{\rm x}^{95\%}$ or $h_{\rm y}^{95\%}$).