††thanks: NHFP Einstein fellow

# A morphology-independent test of the mixed polarization content of transient
gravitational wave signals

Katerina Chatziioannou  Department of Physics, California Institute of
Technology, Pasadena, California 91125, USA  LIGO Laboratory, California
Institute of Technology, Pasadena, CA 91125, USA  Center for Computational
Astrophysics, Flatiron Institute, 162 5th Ave, New York, NY 10010, USA
Maximiliano Isi  LIGO Laboratory, Massachusetts Institute of Technology, 185
Albany St, Cambridge, MA 02139, USA  Department of Physics and Kavli Institute
for Astrophysics and Space Research, Massachusetts Institute of Technology, 77
Massachusetts Ave, Cambridge, MA 02139, USA  Carl-Johan Haster  LIGO
Laboratory, Massachusetts Institute of Technology, 185 Albany St, Cambridge,
MA 02139, USA  Department of Physics and Kavli Institute for Astrophysics and
Space Research, Massachusetts Institute of Technology, 77 Massachusetts Ave,
Cambridge, MA 02139, USA  Tyson B. Littenberg  NASA Marshall Space Flight
Center, Huntsville, AL 35812, USA

###### Abstract

Gravitational waves in general relativity contain two polarization degrees of
freedom, commonly labeled plus and cross. Besides those two _tensor_ modes,
generic theories of gravity predict up to four additional polarization modes:
two _scalar_ and two _vector_. Detection of nontensorial modes in
gravitational wave data would constitute a clean signature of physics beyond
general relativity. Previous measurements have pointed to the unambiguous
presence of tensor modes in gravitational waves, but the presence of
additional generic nontensorial modes has not been directly tested. We propose
a model-independent analysis capable of detecting and characterizing mixed
tensor and nontensor components in transient gravitational wave signals,
including those from compact binary coalescences. This infrastructure can
constrain the presence of scalar or vector polarization modes on top of the
tensor modes predicted by general relativity. Our analysis is morphology-
independent (as it does not rely on a waveform templates), phase-coherent, and
agnostic about the source sky location. We apply our analysis to data from
GW190521 and simulated data and demonstrate that it is capable of placing
upper limits on the strength of nontensorial modes when none are present, or
characterizing their morphology in the case of a positive detection. Tests of
the polarization content of a transient gravitational wave signal hinge on an
extended detector network, wherein each detector observes a different linear
combination of polarization modes. We therefore anticipate that our analysis
will yield precise polarization constraints in the coming years, as the
current ground-based detectors LIGO Hanford, LIGO Livingston, and Virgo are
joined by KAGRA and LIGO India.

##  I Introduction

Gravitational wave (GW) detections by LIGO Aasi  _et al._ (2015) and Virgo
Acernese  _et al._ (2015) have made it possible to test general relativity
(GR) in the dynamically extreme and strong-field regimes Abbott  _et al._
(2016a, b); Yunes  _et al._ (2016); Abbott  _et al._ (2017a); Ezquiaga and
Zumalacárregui (2017); Sakstein and Jain (2017); Creminelli and Vernizzi
(2017); Baker  _et al._ (2017); Abbott  _et al._ (2017b, 2019a, 2019b); Isi
_et al._ (2019); Abbott  _et al._ (2020a). The three instruments in the US and
Europe will soon be joined by KAGRA in Japan Akutsu  _et al._ (2020), while
this decade will also see the addition of LIGO India Iyer  _et al._ (2011). A
growing detector network will bring not only improvements in duty cycle and
sky-localization capabilities, but also enable qualitatively new tests of
gravity: a larger network will allow us to thoroughly probe the polarization
structure of GWs, thus testing one of the key predictions of GR.

In vacuum GR, GWs possess only two polarization states, corresponding to the
two helicities of a massless spin-2 (“tensor”) particle; in the linear basis,
these are conventionally denoted plus ($+$) and cross ($\times$). However,
Einstein’s theory is special in this respect: theories beyond GR generally
introduce additional gravitational degrees of freedom that tend to manifest in
nontensorial polarizations Will (2014). In fact, generic metric theories of
gravity may allow up to six independent GW polarization states: besides the
two tensor modes (helicity $\pm 2$), there may also exist two vector modes
(helicity $\pm 1$), and two scalar modes (helicity 0) Eardley  _et al._
(1973a, b). The two vector modes are often denoted $v_{1}$ and $v_{2}$ (or $x$
and $y$), and the two scalar modes are called breathing ($b$) and longitudinal
($l$). Detecting any of these alternative polarizations, no matter how small
in amplitude, would constitute incontrovertible evidence of nonstandard
physics.

The different polarization states are distinguished by their local geometric
effects. Each mode can be represented as a spatial basis tensor for the metric
perturbation in a local Cartesian frame in which $\hat{z}$ is aligned with the
direction of propagation and $\hat{x}$ and $\hat{y}$ are oriented in some
predetermined way (e.g., rotated by an arbitrary polarization angle $\psi$
relative to some reference), e.g., Eq. (8) in Isi and Weinstein (2017). In the
small-antenna limit (i.e., for detectors such as LIGO and Virgo whose size is
much smaller than the GW wavelength, sometimes referred to as the “long-
wavelength approximation”), the effect of different polarization modes on a
detector can be summarized through antenna pattern functions $F^{p}$ that
encode the sensitivity of the instrument to GWs of polarization $p$ as a
function of its direction of propagation. For an astrophysical source seen by
a global detector network, it is convenient to write these factors as
$F^{p}(\Omega,\psi;t)$, for the source sky location
$\Omega\equiv(\alpha,\delta)$, in terms of the right ascension $\alpha$ and
declination $\delta$, plus the polarization angle $\psi$ defining the linear-
polarization frame within the plane of the sky (conventionally tied to the
orientation of the source in some specific way); we adopt the usual
normalization for the $F^{p}$’s as given by, e.g., Eqs. (12–17) in Isi and
Weinstein (2017). Since the relative orientation of the detector and source
will vary with Earth’s rotation, the $F^{p}$ factors are implicit functions of
time $t$; however, for transient signals (much shorter than a sidereal day),
we can set $t=t_{0}$ for some time of arrival $t_{0}$ and treat the antenna
patterns as constants.

We can take advantage of the antenna patterns to disentangle the polarizations
present in a given signal. The output $h_{I}$ of any given detector $I$ to a
GW with generic polarizations can be written in the time domain as

$h_{I}(t)=\sum_{p}F^{p}_{I}(\Omega,\psi)\,h^{p}(t+\Delta t_{I}(\Omega))\,,$
(1)

where the sum is over polarization states $p=\\{+,\times,v_{1},v_{2},b,l\\}$,
and the $\Delta t_{I}(\Omega)$ represent time of flight delays between
detectors, which also depend on the source sky location. Uncertainty in the
source location will therefore produce uncertainty in the polarization content
due to both the $F^{p}_{I}$ and the $\Delta t_{I}$ quantities. We have
suppressed the time dependence of the antenna patterns by assuming the signal
duration is short, so the only time dependence is introduced by the waveforms
$h^{p}(t)$ for each polarization, which depend on the intrinsic parameters of
the source in a way specific to each theory of gravity. Unfortunately,
differential-arm detectors such LIGO and Virgo are only sensitive to a
specific (traceless) linear combination of the two scalar modes, so networks
of such detectors can distinguish at most _five_ , not six, polarization
states; it is thus sufficient to include only one of the scalar modes ($b$,
without loss of generality) in the sum of Eq. (1).

Crucially, up to a time delay, the $h^{p}$’s are the same for all detectors in
Eq. (1), so a sufficiently large network would allow us to reconstruct the
polarization content of arbitrary signals, and thus detect (or exclude) the
presence of nontensorial modes. In particular, five non-cooriented
differential-arm detectors would be required to invert Eq. (1) and uniquely
determine the five polarizations distinguishable by such instruments, even if
the sky location was perfectly known. On the other hand, if only one or two
detectors are available, Eq. (1) will not be invertible and we will be unable
to break all the (breakable) polarization degeneracies of a transient signal
in a model-independent way. In that case, one could resort to theory-specific
predictions for the waveform phasing to look for smoking-gun features of any
given nontensor mode, but the results would be restricted to that theory.
Alternatively, one could rely on _persistent_ GW signals (such as continuous
GWs from nonaxisymmetric neutron stars or stochastic GW backgrounds) to
distinguish polarizations through the time dependence of the antenna patterns
over a long observation period Isi  _et al._ (2015, 2017); Callister  _et al._
(2017); Abbott  _et al._ (2018a, b); unfortunately, such signals have yet to
be detected.

The small size of the existing global detector network has so far precluded
full-fledged GW polarization studies. This is worsened by the fact that the
two LIGO instruments are nearly coaligned, and thus do not measure fully
independent combinations of polarizations; as a consequence, their joint power
to differentiate GW polarizations is less than that of two independent
detectors Takeda  _et al._ (2018). For this reason, initial detections made
with LIGO alone carried no implications for polarizations Abbott  _et al._
(2016a). It was only once Virgo joined the network that some limited
polarization studies became viable. Although the LIGO-Virgo network is
incapable of fully inverting Eq. (1), it can be sufficient to distinguish
between some extreme polarization alternatives, such as scenarios in which the
polarizations of a given GW are purely tensor, versus purely vector, versus
purely scalar Isi and Weinstein (2017). Such a simplified analysis was carried
out first for GW170814 Abbott  _et al._ (2017b), and then, most notably, for
GW170817 Abbott  _et al._ (2019a), which produced evidence vastly favoring the
full-tensor model over the two nontensor alternatives (full vector, or full
scalar). Different flavors of this test have since been applied to a number of
other events Abbott  _et al._ (2019b); Takeda  _et al._ (2021); Abbott  _et
al._ (2020a); Haster (2020). However, studies targeting spin-weight mixtures
(e.g., combinations of scalar and tensor, or scalar and vector), which are
more interesting theoretically Chatziioannou  _et al._ (2012), have remained
out of reach.111A model-dependent test of mixed tensor-scalar polarizations
was recently proposed Hiroki  _et al._ (2021). However, this study uses an _ad
hoc_ model that does not correspond to the phenomenology of nontensorial modes
in beyond-GR theories—cf. Eq. (5) in Hiroki  _et al._ (2021) with Eqs. (47),
(63), and (79) in Chatziioannou  _et al._ (2012). Beyond issues with the
specific model, a model-dependent approach makes it unclear how much is
learned from the (lack of) observation of nontensorial polarizations rather
than from constraints on the tensor phase evolution. Indeed in Hiroki  _et
al._ (2021) a single parameter was used to simultaneously constrain both the
tensor phase and amplitude, and the scalar amplitude. Since the GW phasing is
measured more accurately than the amplitude, we expect any constraints in
Hiroki  _et al._ (2021) to be dominated by the tensor phase evolution, similar
to Abbott  _et al._ (2019a, b), which is not a direct probe of scalar modes.
Indeed, after a requisite reparametrization, the GW170814 result in Hiroki
_et al._ (2021) is identical to the one in Abbott  _et al._ (2019b), which is
based purely on the phase evolution of tensor modes.

With the advent of KAGRA and, further in the future, LIGO India, there is now
incentive to develop and apply more powerful polarization probes. One possible
avenue is to extend model-independent parametrized tests such as the
parametrized post-Einsteinian framework Yunes and Pretorius (2009), commonly
applied to GW data Abbott  _et al._ (2020a), to allow for nontensorial
polarizations Chatziioannou  _et al._ (2012). A constrained maximum-likelihood
direct inversion of Eq. (1) was explored in Hayama and Nishizawa (2013).
Another possibility related to this is to exploit Eq. (1) to construct a
linear combination of detector outputs that is known a priori to be blind to
tensor GWs from a given sky location Guersel and Tinto (1989); an upper limit
on the signal present in such a _null-stream_ can be translated into a
constraint on nontensor polarizations Chatziioannou  _et al._ (2012); Hagihara
_et al._ (2019); Pang  _et al._ (2020). This approach, used most recently in
Abbott  _et al._ (2020a), has the advantage of not requiring a waveform model
for the nontensor polarizations; however, it also has some important
limitations. First, the original null-stream construct is premised on a known
sky location; although it can be extended to ease this restriction Abbott  _et
al._ (2020a); Wong  _et al._ (2020), that has only been done at the price of
limiting scope to full-tensor vs full-nontensor hypotheses, as for previous
methods. Second, existing implementations rely on phase-incoherent methods to
collect signal power from the null stream, i.e., the power is computed
independently in different time-frequency pixels Sutton  _et al._ (2010). This
means that a signal cannot be tracked coherently over time, and all
information about its phasing is lost. While this might be desirable when
targeting a stochastic signal, such as a white-noise burst, it represents a
severe handicap when it comes to coherent sources, notably compact binaries,
requiring higher signal to noise ratios (SNRs).

In this paper, we present a new approach to study of GW polarizations from
compact binaries that allows us to constrain any combination of helicities,
whether the sky location is known or not, and without sacrificing information
contained in the phase coherence of the signal. We achieve this by
generalizing the BayesWave algorithm, which can reconstruct GW signals of
generic morphology through sine-Gaussian wavelets Cornish and Littenberg
(2015); Cornish  _et al._ (2021), extending it to accommodate nontensor
polarizations. Building upon the robust and flexible BayesWave infrastructure
enables us to overcome all major limitations to previous polarization studies:
we can simultaneously fit for the sky location and polarization content of a
signal to obtain waveform reconstructions for any arbitrary set of helicities,
and, consequently, place bounds on the amplitude of nontensor components even
in the presence of a dominant GR contribution.

In the following sections, we outline the details of the method and
demonstrate its efficacy on a number of simulated GR and non-GR signals as
well as real data from GW190521. The framework can be applied to model any of
the seven possible combinations of tensor (T), vector (V), or scalar (S)
helicities, while maintaining the flexibility to exclude the presence of any
given mode by assigning it zero power. However, since we already have evidence
of the existence of tensor modes in nature (disfavoring the pure V and S
cases), all examples presented here include tensor modes (either TV, TS or
TVS). Given that LIGO Hanford, LIGO Livingston, Virgo and KAGRA (HLVK) are
expected to operate in the upcoming O4 observing run Abbott  _et al._ (2018c),
we will be one instrument short of five detectors; thus, we anticipate results
to be more constraining for the TV and TS than TVS. However, we find that HLVK
will already be sufficient to discern between the three possible spin-weights
within the TVS model, if not to perfectly reconstruct the relative
contributions of the five individual modes—that is, it can successfully
categorize the signal power as tensor, vector or scalar, without fully
breaking all degeneracies between $\\{+,\times,v_{1},v_{2},b,l\\}$. This
situation will further improve once LIGO India comes online later in the
decade Abbott  _et al._ (2018c), forming an HLVKI network that will enable
superior constraints on the full TVS model, as we demonstrate below.

The rest of the paper is organized as follows. In section II we discuss our
methodology in detail, including the tensor and nontensor signal models. In
section III we apply our method to an example GW signal: GW190521. In section
IV we focus on an HLVK network and apply our analysis to signals obeying GR
and estimate the upper limits on beyong-GR modes we can place for TV and TS
analyses. In section V we tackle example signals with nontensorial mode
content again in an HLVK network and demonstrate that our analysis can recover
them and characterize their morphology for TV and TS analyses. In section VI
we turn to full TVS analyses and demonstrate that they are feasible with a
four-detector network, but are considerably strengthened with a full five-
detector HLVKI network. In section VII we conclude.

##  II Methodology and signal model

The nontensor part of a GW waveform may take different forms for different
beyond-GR theories and for different sources, depending on the specific way in
which the theory modifies the GW generation and propagation processes.
Furthermore, modified gravity theories are not only expected to excite
nontensor GW modes, but also to modify the tensor part of the signal compared
to the GR prediction. For this reason, GR templates for the tensor part are
not strictly applicable in beyond-GR theories. In order for our analysis to
remain generic, we use BayesWave Cornish and Littenberg (2015); Cornish  _et
al._ (2021), and model both the tensor and the nontensor polarizations,
$h^{p}$ in Eq. (1), in a generic way. BayesWave employs a sum of sine-Gaussian
wavelets to model astrophysical signals as observed by the detectors. The
number of wavelets and the parameters of each are marginalized over with a
transdimensional sampler, offering the analysis the flexibility to recover
signals for which no accurate models exist.

###  II.0.1 Elliptical Polarization

The simplest signal model in BayesWave assumes that the GW signal is composed
of only the two GR polarizations with relative amplitudes and phases such that
the overall signal is elliptically polarized Cornish and Littenberg (2015). In
that case, the response of a detector $I$ to an impinging GW signal is a sum
of its response to each polarization, expressed in the frequency domain as

$h_{I}(f)=\left[F^{+}_{I}(\Omega,\psi)h^{+}(f)+F^{\times}_{I}(\Omega,\psi)h^{\times}(f)\right]e^{2\pi
if\Delta t_{I}(\Omega)}.$ (2)

In the above equation, $h_{I}(f)$ is the Fourier transform of the output of
detector $I$, and other quantities are the same as in Eq. (1). The signal
itself, $h^{+}(f)$ and $h^{\times}(f)$ is defined at geocenter and then
projected onto each detector; the amplitude of this projection is controlled
by the antenna pattern functions, while the term $e^{2\pi if\Delta
t_{I}(\Omega)}$ expresses the frequency-dependent phase shift of the signal
due to the time delay $\Delta t_{I}(\Omega)$ from geocenter to the detector
location.

In compact-binary-coalescence (CBC) analyses, $h^{+}(f)$ and $h^{\times}(f)$
are most commonly obtained from templates derived under the GR assumption. In
BayesWave’s elliptical polarization model, however, they are given by

$\displaystyle h^{+}(f)$
$\displaystyle=\sum_{n}\Psi(f;t_{0}^{n},f_{0}^{n},Q^{n},A^{n},\phi_{0}^{n}),$
(3) $\displaystyle h^{\times}(f)$ $\displaystyle=\epsilon h^{+}(f)e^{i\pi/2}.$
(4)

where $\Psi(f;t_{0}^{n},f_{0}^{n},Q^{n},A^{n},\phi_{0}^{n})$ is a sine-
Gaussian wavelet with five parameters: a central time $t_{0}^{n}$, a central
frequency $f_{0}^{n}$, a quality factor $Q^{n}$, and amplitude $A^{n}$, and a
phase $\phi_{0}^{n}$; its explicit form is given in Cornish and Littenberg
(2015). The superscript $n$ counts the discrete sum of wavelets.

Equations (3) and (4) imply that the plus waveform is expressed as a sum of
wavelets, while the cross waveform is derived from the plus one under the
assumption that the overall signal is elliptically polarized: $h^{\times}(f)$
is proportional to $h^{+}(f)$, with a phase shift of $\pi/2$ and an amplitude
scale given by the ellipticity $\epsilon$ (which is a real number). In other
words, the cross mode is given by the same wavelets as the plus mode, modulo a
constant phase shift and a constant amplitude scale that is the same for all
wavelets as they share the same ellipticity.

Although the signal model in Eqs. (3) and (4) is morphologically generic, Eq.
(2) is predicated on three assumptions:

  1. 1.

the detectors respond to each polarization mode as prescribed by the
corresponding antenna pattern;

  2. 2.

the signals propagate at the speed of light between detectors; and

  3. 3.

the signals do not disperse as they travel between detectors.

The first point is a restatement of our restriction to metric theories of
gravity, as well as the small-antenna assumption; the latter could be relaxed
if necessary. The second point is imposed by the $e^{2\pi if\Delta
t_{I}(\Omega)}$ factor in Eq. (2), where we compute $\Delta t_{I}(\Omega)$
assuming GWs travel the distance between each detector and the geocenter at
the speed of light. This does not preclude scenarios in which GWs propagate at
speeds $c_{g}\neq c$ over astrophysical distances, as long as the difference
is too small (or fully screened) to be detected over the comparatively short
distances between instruments; existing constraints on $c_{g}/c$ already
validate this assumption Abbott  _et al._ (2017c).222The fact that there exist
six polarizations in generic metric theories of gravity holds for GW
propagation speeds close, but not necessarily equal, to the speed of light
Will (2018). Lastly, the third point establishes that the wave packet should
not measurably change shape during flight between detectors, so that the only
change in the observed morphologies across instruments is explained by the
antenna patterns. This does not exclude beyond-GR models in which GWs disperse
on their long journey from the source.

BayesWave’s elliptical polarization signal model is very flexible and has been
shown to accurately recover GW signals as morphologically distinct as binary
black holes (BBHs) in GR and beyond Ghonge  _et al._ (2020), post merger
signals from BNS coalescences Chatziioannou  _et al._ (2017a); Torres-Rivas
_et al._ (2019), eccentric BBHs Dálya  _et al._ (2021), antichirping signals
Millhouse  _et al._ (2018), and white noise bursts Bécsy  _et al._ (2017).

###  II.0.2 Generic Polarization

The first step to generalize the elliptical polarization model in BayesWave
was taken in Cornish  _et al._ (2021), where generic (nonelliptical) tensor
polarizations were considered. The generic model again expresses the detector
response through a sum of the two GR tensor polarizations as

$h_{I}(f)=\left[F^{+}_{I}(\Omega,0)h^{+}(f)+F^{\times}_{I}(\Omega,0)h^{\times}(f)\right]e^{2\pi
if\Delta_{I}t(\Omega)},$ (5)

where

$\displaystyle h^{+}(f)$
$\displaystyle=\sum_{n}\Psi(f;t_{0}^{n},f_{0}^{n},Q^{n},A^{n}_{+},\phi^{n}_{0,+}),$
(6) $\displaystyle h^{\times}(f)$
$\displaystyle=\sum_{n}\Psi(f;t_{0}^{n},f_{0}^{n},Q^{n},A^{n}_{\times},\phi^{n}_{0,\times}).$
(7)

There are two critical differences between the elliptical and the relaxed
polarization model. First, the ellipticity parameter is absent and the
polarization angle has been set to zero throughout.333Setting $\psi=0$ does
not lead to loss of generality because the only purpose of $\psi$ is to define
the arbitrary orientation of the polarization frame, and its effect can be
fully absorbed by
$(A^{n}_{+},\phi^{n}_{0,+},A^{n}_{\times},\phi^{n}_{0,\times})$. See, e.g.,
App. A in Isi  _et al._ (2017). Second, each of the plus and cross
polarization is expressed as a sum of wavelets. The two sums are not
completely independent, but consist of the same number of wavelets that also
share time, frequency, and quality factor parameters. The amplitude and phase
of the wavelets is different for each mode. Effectively, the generic
polarization model replaces the ellipticity parameter and polarization angle
of the elliptical polarization model with a series of amplitudes and phases
for each wavelet. This could be equivalently parametrized as in Eqs. (3) and
(4) if the phase shift and ellipticity were made wavelet-specific parameters.

The wavelets in Eqs. (6,7) have the same quality factor and central frequency
for the plus and cross components. One way to interpret this is to take each
set of $\\{t_{0}^{n},f_{0}^{n},Q^{n}\\}$ to define a single, elliptically
polarized wavelet with plus and cross components of arbitrary amplitude and
phase, as given by the summands in Eqs. (6,7). This parametrization is
convenient because we do not expect the spectral content of the two
polarizations to be totally independent. This is true because (1) the same
physical processes that generate plus also generate cross, and (2) even if, in
principle, there could exist a source for which plus and cross were totally
morphologically independent in _some_ frame (i.e., for some unknown choice of
$\psi$ specific to each source) the plus and cross polarizations measured by
generic observers would be linear combinations of the two independent
components and, therefore, would look spectrally similar. In other words, even
if there existed a frame orientation in which the morphology of plus and cross
looked totally independent for some source, this will not be the case in the
BayesWave frame. Nevertheless, if somehow we detected a signal in which plus
and cross were indeed truly independent, the model of Eqs. (6,7) can still fit
that by increasing the number of wavelets.

Compared to the elliptical polarization, the generic polarization model has
more parameters and hence more flexibility to model morphologically complex
signals. The elliptical polarization model has $5n+4$ parameters, where $n$ is
the number of wavelets: five parameters intrinsic to each wavelet
($t_{0}^{n},f_{0}^{n},Q^{n},A^{n}_{+},\phi^{n}_{0})$, and four shared
parameters $(\alpha,\delta,\psi,\epsilon)$. In contrast, the relaxed
polarization model has $7n+2$ parameters: seven intrinsic wavelet parameters
($t_{0}^{n},f_{0}^{n},Q^{n},A^{n}_{+},\phi^{n}_{0,+},A^{n}_{\times},\phi^{n}_{0,\times})$,
and only two shared parameters ($\alpha,\delta$). For any given $n>1$, the
generic polarization model will always have more parameters than the
elliptical polarization one. Although this means the generic model is able to
fit a greater range of signal morphologies, this flexibility comes at the
price of an increased Ockham penalty per wavelet. For given data the
$(n+1)^{th}$ wavelet will result in a larger reduction in the prior (i.e. and
Ockham penalty) compared to the $n^{th}$ wavelet when using the generic
polarization model than it would with the elliptical model. This drop in the
prior needs to be compensated by an increase in the likelihood as the
$(n+1)^{th}$ wavelet presumably models some non-Gaussian feature in the data.
The result is that a wavelet in the generic polarization model needs on
average to result in a larger increase in the likelihood than a wavelet in the
elliptical polarization model in order to contribute to the posterior.

These considerations suggest that the generic polarization model has the
flexibility to model complicated signals, but might underperform the
elliptical polarization model for simple signals that are indeed elliptically
polarized. Tests of the generic polarization model were presented in Cornish
_et al._ (2021) where it was shown to more accurately reproduce spin-
precessing compact binary signals (which are not elliptically polarized) than
the elliptical polarization model. It was further demonstrated that the signal
tensor polarization content can be extracted in terms of the signal Stokes
parameters.

###  II.0.3 Beyond-GR Polarization

We extend the generic tensor model above to include nontensor polarizations
modes in the detector response. The resulting beyond-GR polarization model is
expressed as

$\displaystyle h_{I}(f)=$
$\displaystyle\left[F_{I}^{+}(\Omega,0)h^{+}(f)+F_{I}^{\times}(\Omega,0)h^{\times}(f)\right.$
$\displaystyle\left.+F_{I}^{v1}(\Omega,0)h^{v1}(f)+F_{I}^{v2}(\Omega,0)h^{v2}(f)\right.$
$\displaystyle\left.+F_{I}^{b}(\Omega,0)h^{b}(f)\right]e^{2\pi
if\Delta_{I}t(\Omega)},$ (8)

where the first line is the usual tensor part from the generic polarization
model. The second and third lines correspond to the vector and scalar modes
respectively. The detector response is characterized by the usual antenna
pattern functions for the two vector modes and the breathing mode,
$F^{v1}(\Omega,0),F^{v2}(\Omega,0),F^{b}(\Omega,0)$. As mentioned above, while
in general there are two possible scalar modes -the breathing and the
longitudinal- the response of interferometric detectors to them is degenerate,
so we follow common practice and consider only one; we will refer to it as the
scalar mode from now on. As before, all antenna pattern functions are
evaluated at zero polarization angle (the scalar antenna pattern function is
already independent of the polarization angle).

Our model for the geocenter signal for each polarization mode is

$\displaystyle h^{+}(f)$
$\displaystyle=\sum_{n_{T}}\Psi(f;t_{T}^{n},f_{T}^{n},Q_{T}^{n},A^{n}_{T+},\phi^{n}_{T+}),$
(9) $\displaystyle h^{\times}(f)$
$\displaystyle=\sum_{n_{T}}\Psi(f;t_{T}^{n},f_{T}^{n},Q_{T}^{n},A^{n}_{T\times},\phi^{n}_{T\times}),$
(10) $\displaystyle h^{v1}(f)$
$\displaystyle=\sum_{n_{V}}\Psi(f;t_{V}^{n},f_{V}^{n},Q_{V}^{n},A^{n}_{Vv1},\phi^{n}_{Vv1}),$
(11) $\displaystyle h^{v2}(f)$
$\displaystyle=\sum_{n_{V}}\Psi(f;t_{V}^{n},f_{V}^{n},Q_{V}^{n},A^{n}_{Vv2},\phi^{n}_{Vv2}),$
(12) $\displaystyle h^{b}(f)$
$\displaystyle=\sum_{n_{S}}\Psi(f;t_{S}^{n},f_{S}^{n},Q_{S}^{n},A^{n}_{Sb},\phi^{n}_{Sb}).$
(13)

The first two lines are the two tensor modes and they are identical to the
generic polarization case (cf. Eqs. ((6), (7))). The third and fourth lines
correspond to the two vector modes and are constructed similarly to the two
tensor modes: each mode is expressed as a sum of the same number of wavelets
who share their time, frequency, and quality factor parameters, but have
independent amplitudes and phases. The last line corresponds to the sole
scalar mode which is given by yet another sum of wavelets whose number and
parameters are independent from those of the tensor and the vector modes. To
remain generic, we also allow for independent time, frequency, and quality
factor between the tensor, the vector, and the scalar modes. Indeed existing
examples of polarization states in beyond-GR theories predict different time-
frequency content for each polarization mode Chatziioannou  _et al._ (2012).

Though we choose to construct our beyond-GR model based on the generic
polarization model, an alternative choice would be to make use of the
elliptical polarization model. In this case, the two vector modes would also
be elliptically polarized, but their ellipticity would be independent from the
ellipticity of the tensor modes and would thus be expressed as a new model
parameter. While in this study we choose to remain generic and not assume
elliptical polarization, we will explore the implications of this choice in
the future.

The beyond-GR polarization model employs the same priors for the wavelet
parameters as the generic polarization model: uniform in time, frequency,
phase, and quality factor, while the amplitude of each wavelet has a prior
determined by the SNR of the wavelet, as described in Cornish  _et al._
(2021). The prior on the total number of wavelets in all modes is flat between
$[1,\,50]$. This prior is based on the _total_ number of wavelets which means
that any individual mode (tensor, vector, scalar) can have zero wavelets, as
long as the total number in all modes is greater than one.

In the discussion below, it will be useful to quantify the signal power
contributed by each spin-weight, i.e. in modes of a common helicity up to
sign. Respectively for tensor, vector and scalar, this is at detector $I$

$\mathrm{SNR}^{T}_{I}=\left|\left|F^{+}_{I}h^{+}+F^{\times}_{I}h^{\times}\right|\right|\,,$
(14)
$\mathrm{SNR}^{V}_{I}=\left|\left|F^{v_{1}}_{I}h^{v_{1}}+F^{v_{2}}_{I}h^{v_{2}}\right|\right|\,,$
(15) $\mathrm{SNR}^{S}_{I}=\left|\left|F^{b}_{I}h^{b}\right|\right|\,,$ (16)

where the norm $\left|\left|x\right|\right|\equiv\sqrt{\left\langle
x|x\right\rangle}$ is defined in terms of the usual frequency-domain inner
product weighted by the one-sided noise power spectral density $S_{n}(f)$,

$\left\langle x\mid
y\right\rangle=4\Re\int_{0}^{\infty}\frac{x(f)y^{*}(f)}{S_{n}(f)}\mathrm{d}f\,,$
(17)

for any two Fourier domain functions $x(f)$ and $y(f)$. As with the total
signal SNR, we define the network spin-weight SNR as the quadrature sum of
spin-weight SNRs at each detector.

The SNRs defined above are a measure of the signal power _projected_ onto each
detector by each set of polarizations. The measurement of this “projected”
power ratio can then be used to estimate the intrinsic power ratio (i.e., the
ratio of scalar or vector to tensor GW power emitted by the source), once the
inferred sky location is taken into account.

In the presence of a signal, we may ask which (if any) of the polarization
models is favored. This is a model selection question that can formally be
answered by evaluating odds or marginalized-likelihood ratios (Bayes factors),
quantities natively computed by BayesWave. However, interpreting marginalized
likelihoods is more challenging than interpreting features of the posterior.
The fundamental reason is that Bayes factors depend on prior choices in a way
that is highly sensitive even to regions of the parameter space where the
likelihood offers no support.444In contrast, the posterior is affected by the
prior only in regions of nonnegligible likelihood support, and therefore the
effect of the prior is less pronounced and easier to interpret. For instance,
for any set of data, one can arbitrarily increase the preference for tensor
waves by sufficiently extending the prior range of the nontensorial
amplitudes. Therefore, Bayes factors are generally only meaningful when every
aspect of the prior can be set robustly, based either on first principles or
large scale simulations with realistic populations of signals. In our case,
carrying out such simulations would be impossible without committing to very
narrow, likely unrealistic, models of GR deviations and source populations. In
any case, Bayes factors do not provide an answer to the main question in which
we are interested: how strong could a nontensorial component of the signal be
and still be consistent with the data? For that reason, we focus on upper
limits to the SNR and do not present Bayes factors below.

##  III Gravitational Wave Events: GW190521

As a first application of the beyond-GR polarization model, we analyze the GW
signal GW190521 Abbott  _et al._ (2020b, c), the heaviest BBH detected to date
and the only one with a proposed accurate sky localization thanks to a
candidate optical electromagnetic counterpart Graham  _et al._ (2020). Below
we assume GW190521 indeed originates from this sky location and discuss the
results in detail as an example of the typical expected outcome of our
analysis on GW data. The quantitative results, however, are contingent on
GW190521 truly being associated with the candidate counterpart Graham  _et
al._ (2020), something that is still under debate Ashton  _et al._ (2020);
Palmese  _et al._ (2021). The physical applicability of the results of this
section is therefore limited and we approach them mainly as an instructive
exercise showcasing the methodology.

We analyze $4\,$s of data from the LIGO Hanford, LIGO Livingston, and Virgo
detectors around the time of the event from $16\,$Hz to $512\,$Hz Abbott  _et
al._ (2019c). The power spectral density of the Gaussian detector noise is
marginalized over using the BayesLine algorithm Littenberg and Cornish (2015);
Chatziioannou  _et al._ (2019). We assume that the sky location of the source
is inferred from its electromagnetic counterpart candidate Graham  _et al._
(2020): the right ascension $\alpha$ is 3.36 radians and the declination
$\delta$ is 0.61 radians. This sky location was selected in Graham  _et al._
(2020) to be consistent with the 3D source localization volume estimated from
GW data alone while assuming GR Abbott  _et al._ (2020b); Chen  _et al._
(2020); Ashton  _et al._ (2020); Haster (2020). Though our method does not
require a known sky location (see for example the simulated signals in Sec.
IV), that additional information is particularly helpful here given that only
three detectors were operational at the time of GW190521.

We analyze the data with three possible combinations of polarization modes:
pure tensor content corresponding to the GR prediction without assuming an
elliptically polarized signal (T), a mixed tensor and scalar content (TS), and
a mixed tensor and vector content (TV). Our results are presented in Figs. 1
and 2. As expected, we find that the T model is sufficient to describe the
signal. Indeed, the TS and TV analyses result in 71% and 77% of the posterior
samples having exactly zero wavelets for the beyond-GR polarization mode
respectively. The remaining fraction of samples contain some weak nontensorial
contribution which is consistent with statistical noise, as we show below.

Figure 1: Whitened signal reconstructions for GW190521 in each detector
relative to GPS time 1242442967.45. We plot the whitened data in grey dashed
lines, while shaded regions denote the 90% credible intervals using the tensor
(T; pink), tensor+scalar (TS; green), and tensor+vector (TV; blue) model. All
reconstructions are consistent with each other suggesting no deviation from
the GR polarization prediction.

Figure 1 shows the whitened signal reconstruction as a function of time in
each detector. The shaded regions give the 90% credible intervals for the
signal reconstruction under each polarization content model, with the T model
corresponding to the GR prediction. Compared to the analysis presented in
Figure 1 of Abbott  _et al._ (2020b), the T model here allows for generic plus
and cross polarizations without requiring an elliptical polarization. The TS
and TV reconstructions are consistent with the T one, with most of the SNR
placed in the tensor wavelets in both cases. We compute the network overlap
Ghonge  _et al._ (2020) between the median TS (TV) signal reconstruction and
the median T reconstruction and find 0.998 (0.995), showing high agreement
between the two. We also compute network overlaps between random draws from
the reconstruction posterior for the TS or TV models with draws from the T
model and find them to be statistically indistinguishable from overlaps
between random draws from the T distribution alone. This suggests that the
data do not offer any significant evidence for a deviation from GR in the
signal polarization content and tensor modes alone are sufficient to model the
data. We also compute the scalar (vector) part of the reconstruction in the TS
(TV) analysis and find it to be consistent with zero at the 90% level. If
plotted in Fig. 1, both scalar and vector reconstructions would be a straight
horizontal line at $0$, we thus omit them for clarity.

Figure 2: Network matched-filter SNR posterior for beyond-GR polarization
modes in GW190521, as defined in Eqs. (15) and (16). The inset shows the
inherent SNR prior and posterior for the vector case for the analysis. In the
main figure the posterior has been reweighted to a flat SNR prior. Numbers in
the legend of the main plot give the 90% upper limit on the SNR.

We quantify the upper limits on the beyond-GR polarization content through the
SNR of the scalar and vector modes in the TS and TV analyses respectively in
Fig. 2. We plot the network matched-filter SNR posterior and prior for each
analysis. In both cases we find that the SNR in beyond-GR polarizations is
consistent with zero, again indicating no evidence for a beyond-GR
polarization. The main figure shows the SNR posterior under a flat prior. The
corresponding 90% upper limits on the SNR of beyond-GR polarizations in
GW190521 is 4.1 and 4.7 for scalar and vector modes respectively. The
corresponding SNR in tensor modes is $15.3^{+1.5}_{-1.5}$
($15.3^{+1.5}_{-1.5}$) while the total SNR in all modes is
$15.4^{+1.6}_{-1.5}$ ($15.4^{+1.6}_{-1.5}$) for the TS(TV) analysis, in
agreement with Abbott  _et al._ (2020b, c) and the expectation that the signal
can be fully explained by tensor modes.

The inset in Fig. 2 shows the default SNR prior used by BayesWave during
sampling, and the corresponding posterior for the vector analysis; we obtain
the SNR distribution shown in the main figure through reweighing of the
posterior by the prior from the inset. The inherent prior in BayesWave imposes
there to be at least one wavelet present across all the included polarization
states, though this allows for $0$ vector wavelets provided that there is at
least $1$ wavelet in the tensor model. As a result, the SNR prior (pink) has a
sharp peak at $0$, corresponding to no vector wavelets. In addition, the SNR
prior has a second, smooth mode with a peak near SNR 10 and broad support up
to comparatively high SNR values, which corresponds to at least one vector
wavelet. The SNR posterior (purple) is shifted toward lower SNRs: the zero-
wavelets peak is enhanced and the smooth part is concentrated at lower SNR
values. Combined, these observations suggest that BayesWave indeed infers a
reduced posterior support for vector modes relative to the default prior, with
the data disfavoring the presence of vector power in the observed signal.

Though the above results and upper limits present a constraint on the presence
of a mixed polarization content in an observed GW signal, they still depend on
the sky location inferred from the counterpart candidate. We therefore treat
them mainly as an example of the constraints on beyond-GR theories possible
with our analysis and only consider them as an astrophysical measurement under
the caveat that the association of the GW and the counterpart signal is still
uncertain.

##  IV HLVK network and simulated signals in GR

After demonstrating the potential for mixed polarization constraints using
GW190521 as an example in Sec. III, in this section we quantify the expected
constraints from a 4-detector network using simulated signals. We consider a
detector network consisting of LIGO Hanford, LIGO Livingston, Virgo, and
KAGRA, all of which are expected to be operational during the upcoming fourth
observing run Abbott  _et al._ (2018c). For simplicity we assume that each
detector is operating with its design sensitivity Evans  _et al._ (2020), an
assumption that is almost certainly too optimistic at least for KAGRA. We do
not assume a priori knowledge of the sky location of any signal in this
section, instead the sky location is marginalized over during sampling.

Figure 3: Matched-filter SNR 90% upper limits for the scalar (left) and vector
(right) polarization modes for GW150914-like (top) and GW151226-like (bottom)
signals, as function of right ascension $\alpha$ and declination $\delta$ of
the source. The simulated signals share the same parameters apart from the sky
location and the distance, which is changed so as to keep the SNR constant at
50. Each signal is injected in separate and independent random Gaussian noise
realizations. The nontensorial SNR upper limit does not have a strong
dependence on the sky localization, confirming that a 4-detector HLVK network
has a nearly uniform sky coverage.

We begin by exploring the upper limits we can place on the presence of a
beyond-GR polarization content in the case of signals that obey GR. We
simulate signals using the IMRPhenomXPHM model Pratten  _et al._ (2020a, b);
García-Quirós  _et al._ (2020, 2021) with sky locations randomly chosen from a
distribution that is uniform on the sphere. The intrinsic parameters are
chosen to be consistent either with a GW150914- Abbott  _et al._ (2016c) or
with a GW151226-like555The detector-frame masses of the GW151226-like
simulation were multiplied by a factor $1.3$ relative to the original GW151226
signal in order for the simulation to be contained within a data segment of 8s
duration. The GPS time of the GW150914-like signals is $1126259462.42$ and
that of the GW151226-like ones is $1135136350.63$ binary Abbott  _et al._
(2016d), in both cases with BH spins restricted to be aligned with the orbital
angular momentum. The two signals were chosen as representative examples of
typical high-mass and low-mass BBH systems respectively, while the SNR for all
simulated signals is fixed to 50 across the detector network. Since we are
interested in upper limits, our simulated data include a realization of
Gaussian noise colored by the assumed power spectral density of each detector
Evans  _et al._ (2020). For reference, we expect that the exact noise
realization affects the inferred posterior by ${\cal{O}}(\mathrm{SNR})$ in the
large SNR regime. We analyze all signals with a tensor+scalar (TS) and a
tensor+vector (TV) signal polarization model. Figure 3 shows the SNR 90% upper
limits as a function of the sky location for GW150914-like signals (top) and
GW151226-like signals (bottom) and for scalar part of the TS (left) and the
vector part of the TV (right) analysis.

In all cases we find that the expected SNR upper limits are $\sim(0,10)$,
suggesting that the upper limit obtained for GW190521 in the previous section
is typical even if the simulated signals here have a higher SNR (50) than
GW190521 (15), as we also expand on below.666The simulated signals in this
section were analyzed assuming a network of $4$ detectors, while only $3$
detectors were operational at the time of GW190521. However, unlike the
simulated signals, the GW190521 analysis assumed a known sky location.
Additionally, we find that upper limits on nontensorial modes are
approximately independent of an injection’s sky location. The sensitivity of a
network of two detectors has a strong directionality at a given time due to
the shape of the antenna pattern functions. But here we find that four
detectors across the globe achieve more uniform sky coverage, in agreement
with previous studies Fairhurst (2011); Pankow  _et al._ (2018, 2020). All SNR
upper limits in Fig. 3 are computed assuming the intrinsic prior BayesWave
uses on the SNR, and are not reweighted to a flat-in-SNR prior. They would
therefore correspond to the inset of Fig. 2. Additionally, some points have an
upper limit of exactly $0$; this is because at the 90% credible level the
beyond-GR polarization mode has $0$ wavelets, and thus SNR identically zero.

Figure 4: Matched-filter SNR 90% upper limits for the scalar (green) and
vector (red) polarization modes for a GW150914-like signal. The simulated
signals share the same parameters as the GW150914-like signals shown in Fig.
3, but with a single sky location and a varied distance. Each signal has been
injected into $20$ independent noise realizations. The solid line represents
the median upper limit across the different noise realizations, with the
shaded regions showing the 50% interval.

The upper limits presented in Fig. 3 were calculated with simulated signals at
an SNR of 50 across the detector network. Since these signals were simulated
assuming GR, this SNR is entirely due to tensor modes. To explore the
dependence of upper limits on the SNR of the tensor signal, we select a sky
location and simulate a GW150914-like signal at varying SNRs. The upper limit
on the SNR depends sensitively on the exact realization of Gaussian detector
noise. In the high SNR regime, the exact noise realization affects the SNR by
order ${\cal{N}}(0,1)$. However, in the low scalar/vector SNR regime that we
explore here, the dependence is more severe. Upper limits can thus be
particularly sensitive on the noise realization, as shown in the inset of Fig.
2. Indeed, the scalar or vector model can switch between zero and one
wavelets, causing the SNR itself to jump between zero and some value close to
five. As expected, we find that different noise realizations can result in a
small preference for zero or one wavelets, with a corresponding change in the
SNR upper limit. For this reason we repeat our analysis for $20$ noise
realizations for each injected tensor SNR.

Figure 4 shows the upper limit on the scalar and vector SNR as a function of
the SNR of the simulated GR signal, which is solely due to tensor modes.
Similar to Fig. 3, the SNR here has not been reweighted to a flat prior. The
upper limits we can place on beyond-GR polarization modes have a weak
dependence on the SNR of the tensor signal. For tensor SNRs below ${\sim}40$,
we find that most noise realizations result in the 90% upper limit on the SNR
of scalar and vector modes being zero, since the model favors no nontensor
wavelets. For tensor SNRs above ${\sim}40$, we find an upper limit on the
beyond-GR SNR of about ${\sim}5$, consistent with the addition of a single
wavelet at the 90% level. Additionally, the scalar SNR upper limit is
generally higher than the corresponding vector upper limit. The reason is that
a scalar wavelet has five parameters, while a vector (or a tensor) wavelet has
seven parameters. By parsimony, therefore, it is easier to add a scalar
wavelet than a tensor or vector one: when the injected GR signal is quite
loud, requiring more wavelets overall, the preference toward scalar wavelets
can lead to the addition of a scalar wavelet even if the fit would have been
slightly better with a tensor one. On the other hand, any leakage of tensor
SNR into the vector modes is driven mainly by random fluctuations.

Modulo the behavior intrinsic to the way that BayesWave places wavelets
(including discreteness in the SNR prior, and parsimony considerations), these
results are in agreement with the general expectation that, with enough
detectors, the power in tensor and nontensor modes should roughly be
orthogonal, regardless of waveform specifics. That was also the conclusion of
previous studies of continuous waves Isi  _et al._ (2017) and stochastic
backgrounds Callister  _et al._ (2017), where upper limits on the amplitudes
of nontensor modes were found to be independent of the tensor power, and
constraints on the ratio of amplitudes of the different polarization models
depended only on the total SNR of the signal.

##  V HLVK network and simulated signals beyond GR

Besides placing upper limits on beyond-GR physics, any test of the nature of
gravity should be capable of detecting and characterizing potential
deviations. Our infrastructure models the beyond-GR polarization modes in a
morphology-independent way, and should thus have the flexibility to recover
the waveform of any potential nontensor mode in the data, whatever its
specific phase evolution. To demonstrate this, we simulate signals that
violate GR through the presence of nontensorial polarization states. Since no
fully consistent calculation of a GW signal from a binary coalescence exists
in a theory with nontensor modes, the choice of simulated non-GR signals will
inevitably be somewhat ad hoc. We explore two phenomenologies for the
nontensor polarizations: (i) a “burst” morphology, where bursts of scalar or
vector radiation are emitted at various points during the coalescence; and
(ii) a “chirp” morphology, where the scalar or vector signal resembles a
typical GW chirp.

We select five random sky locations and assume the tensor part of the signal
to be given by the NRSur7dq2 Blackman  _et al._ (2017) model, a surrogate to
numerical relativity simulations in GR. We simulate signals with parameters
consistent with GW150914,777Component masses $m_{1}=38\,M_{\odot}$ and
$m_{2}=31\,M_{\odot}$; spin magnitudes $\chi_{1}=0.88$ and $\chi=0.49$; spin
tilts $\theta_{1}=1.67$ and $\theta_{2}=1.66$; azimuthal inter-spin angle
$\phi_{12}=4.19$; angle between orbital and total angular momenta
$\phi_{JL}=6.24$; inclination $\theta_{JN}=2.67$; reference phase
$\phi_{c}=2.41$; reference frequency $f_{\rm ref}=20\,\mathrm{Hz}$. and total
SNRs of 20 and 70 in an HLVK network. These SNR values correspond to a signal
of moderate strength where nontensorial content is barely detectable and a
strong signal where nontensorial modes are clearly identifiable, respectively.
For this demonstration, the tensor part of the simulated signals will obey GR;
however, the recovery is not informed or limited by this fact, remaining fully
general per Eqs. (9) and (10). We choose a tensor signal that obeys GR for the
same as the reason we choose ad hoc “burst” and “chirp” morphologies for the
beyond-GR polarizations: self-consistent calculations that can predict the
full inspiral-merger-ringdown signal in a theory beyond-GR are still in their
infancy, though recent progress has been made Okounkova  _et al._ (2019,
2020); Okounkova (2020); Witek  _et al._ (2020). Furthermore, we generally
expect the tensor part of the signal to remain close to GR - at least
qualitatively- even if additional helicities are present Chatziioannou  _et
al._ (2012); Abbott  _et al._ (2020a).

We simulate the non-GR polarization content of the signal as either scalar of
vector modes with a burst or a chirp morphology. The burst case corresponds to
two sine-Gaussian bursts of non-GR power: one centered 0.5 s before merger,
with central frequency 50 Hz and damping time 0.1 s; and another centered at
merger, twice as loud, with central frequency 100 Hz and the same damping
time. On the other hand, for the chirp case we use NRSur7dq2 also for the
nontensor part, but evaluate it for component masses twice as heavy as those
of the tensor signal—this results in a nontensor signal that mimics the tensor
part but with frequencies twice as low. This approach is inspired by the
calculations of Chatziioannou  _et al._ (2012) that show that the frequency
evolution of the scalar and vector signals is also given by harmonics of the
orbital phase, at least during the inspiral. Similar to the tensor part, the
specific choice of morphology for the nontensor part has no bearing on the
wavelet-based recovery, which is fully flexible.

Figure 5: Signal reconstruction for a beyond-GR injection in a 4-detector
network that consists of tensor and scalar modes. The signal is analyzed with
the TS model and we plot in shaded regions the 90% credible interval for the
full reconstruction (TS; pink), the tensor part of the reconstruction (T;
orange), and the scalar part of the reconstruction (S; green) in the Virgo and
the KAGRA detectors. The scalar mode is absent in Virgo but very prominent in
KAGRA thanks to their different sky responses, while the tensor mode has a
similar strength in both. This allows us to distinguish between the different
polarization modes of the GW signal.

We begin with a detailed discussion of a specific simulated signal in order to
lay out the main features of how polarization modes are extracted from GW
data. Figure 5 shows the reconstruction for one of the simulated signals with
tensor and scalar modes of the burst morphology, with a total SNR of 70. We
analyze the signal with the TS polarization model in a 4-detector network,
without assuming a known sky location, and plot the whitened reconstructions
for Virgo and KAGRA. The observed signal is morphologically different in each
detector due to the fact that the instruments respond differently to each
polarization mode from this sky location, as encoded in the respective antenna
patterns. In the case at hand, the Virgo detector has a negligible response to
scalar polarizations from the direction of this source, resulting in an
observed signal that is dominated by the tensor polarization and looks similar
to a typical merging BBH. On the other hand, KAGRA has a strong scalar
response for this source location, due to its different relative orientation;
as a result, the signal observed by KAGRA is a combination of tensor and
scalar polarizations and no longer resembles a merging BBH as described by GR.
We chose this extreme example, in which one of the detectors is almost totally
blind to one of the targeted helicities, to demonstrate the general fact that
the variation in the signal as seen by each detector, combined with the
arrival times at different instruments, allows us to simultaneously infer the
sky location of the source and its polarization content.

The colored shaded regions in Fig. 5 show the 90% credible intervals for the
signal reconstruction. Colors correspond to different helicity components
derived from the TS reconstruction. In pink, we plot the full reconstruction,
including both tensor and scalar modes: since the TS model matches the
polarization content of the simulated signal, this reconstruction agrees with
the data throughout, as expected. The green shaded region shows the scalar
part of the TS reconstruction: we recover minimal scalar signal in Virgo, but
a strong scalar signal in KAGRA. Finally, the orange region shows the tensor
part of the reconstruction: in Virgo, where there is no sensitivity to the
scalar mode, the tensor part is coincident with the full TS reconstruction and
resembles a merging BBH; although, the chirping morphology is also present in
KAGRA, it is obfuscated by the strong scalar component. Indeed, in KAGRA
neither the T nor S parts of the model fully reproduce the simulated data
alone, whereas TS does.

Figure 6: Breakdown of the signal from Fig. 5 in terms of the different
injected and recovered polarization modes. The full signal consists of
consists of tensor and scalar modes and it is analyzed with the TS model. We
plot in shaded regions the 90% credible interval for the full reconstruction
(TS; pink), the tensor part of the reconstruction (T; orange), and the scalar
part of the reconstruction (S; green) in the KAGRA detector. Same-color dashed
lines denote the corresponding polarization content of the injected signal. In
all combinations of polarization modes the injected and recovered signals
overlap.

To further demonstrate that this procedure leads to a reliable separation of
the different polarization modes, in Fig. 6 we break down the bottom panel of
Fig. 5 in terms of different combinations of polarization modes. We again plot
the full injected data and recovered reconstruction (TS; pink) and also the
injected and recovered part of the signal corresponding to T (orange) and S
(green) modes. We find that both the T and S parts of the reconstruction agree
with the corresponding modes of the injected signal throughout. This suggests
that our analysis can robustly recover not only the full signal, but also
separate its polarization components.

Figure 7: Whitened signal reconstructions for a TS (top) and a TV (bottom)
injection in each of the $4$ detectors relative to GPS time $1126259461$. The
top panel corresponds to an injected sky location of
$(\alpha,\delta)=(2.750,0.923)$ and the bottom to
$(\alpha,\delta)=(4.901,0.720)$. The left panels show a burst injection while
the right panels show a chirp nontensor injection. The solid grey line shows
the simulated data in each detector. Pink shaded regions show the 90% credible
interval for the reconstruction of the full signal containing both tensor and
scalar (left) or vector (right) modes. The scalar (S; left) or vector (V;
right) modes of this full signal reconstruction along are shown with green
shaded regions.

To discuss the features of the nontensor reconstruction in a bit more detail,
we also consider an SNR 70 TV burst injection from a different sky location
together with the scalar case in Fig. 5 and their corresponding chirp analog;
results for other sky locations are qualitatively similar. Figure 7 shows the
whitened reconstructions in each detector for the TS (top) and the TV (bottom)
injection of the burst (left) and chirp (right) type. In all cases, the full
reconstruction is consistent with the injected signal, regardless of whether
it exhibits a strong nontensor content. Comparing the beyond-GR
reconstructions across interferometers shows how the strength of each mode
differs per detector, allowing us to separate them. Comparing the burst and
chirp injections (left and right) panels shows that our analysis can also
separate tensor and nontensor modes even in the case where they are
morphologically similar and overlapping.

Figure 8: Spectra for a TS simulated signal. The left panel shows a burst
injection while the right panel shows a chirp injection. The black line is the
power spectral density of the noise, while the grey line is the spectrum of
the simulated data. The pink shaded region shows the spectrum of the full TS
model, while the green region shows only the S part of the full TS model.

Besides the time-domain reconstruction, we can also characterize the frequency
content of the signal. Figure 8 shows the spectrum of the TS injection for the
burst (left) and the chirp (right) injection; we find qualitatively similar
results for the TV injection. Again in all cases the full TS reconstruction
agrees with the full simulated data. For the burst case (left panel), the
nontensor reconstruction clearly recovers the two distinct bursts of scalar
radiation centered at $50\,$Hz and $100\,$Hz, picked up by the S part of the
model. In the chirp case (right panel), the interference between the tensor
and scalar signals results in a modulation of the overall signal spectrum that
is akin to the beating induced by spin-precession in GR, e.g., Apostolatos
_et al._ (1994); Schmidt  _et al._ (2012); Chatziioannou  _et al._ (2013);
Schmidt  _et al._ (2015); Chatziioannou  _et al._ (2017b); Fairhurst  _et al._
(2019). This modulation is absent in Virgo, since that instrument is
insensitive to scalar waves from this sky location.

Figure 9: Tensor part of the signal reconstruction from different runs using
TS (left) and TV (right) data and analyses. Grey lines show the tensor part of
the data. The orange shaded regions shows the T reconstruction from an
analysis of that exact tensor data. The purple shaded regions show the tensor
part of the TS (left) and TV (right) analysis of the burst injections of Fig.
7. All tensor reconstructions agree with each other, showing that our analysis
can accurately recover the tensor signal even when the simulated data contain
nontensor power.

The reconstruction and spectrum plots presented above suggest that the
scalar/vector part of the model successfully recovers the scalar/vector part
of the simulated signal, irrespective of the tensor modes. To further
demonstrate the absence of cross-helicity contamination, in Fig. 9 we simulate
signals that contain only the tensor part of the burst simulations in Fig. 7
(left panels), and analyze them with the T model. We then compare the pure-
tensor data (grey line), the T reconstruction (orange), and the tensor part of
the TS (left) or TV (right) reconstruction (purple) from the analyses in Fig.
7. If the TS (TV) reconstructions were successful, we should be able to
cleanly separate the T and S (V) parts, and find that the T part of the TS
(TV) reconstruction agrees with the T reconstruction obtained from T-only
data. Indeed, all reconstructions in Fig. 9 agree with each other, even though
the TS/TV models result in larger uncertainties, consistent with their higher
number of degrees of freedom. This shows that we can separate the tensor and
nontensor signal components and thus faithfully characterize the signal
morphology.

Figure 10: Posterior for the total SNR and the SNR per mode for simulated TS
(left) and TV (right) signals of the chirp morphology. Solid vertical lines
and histograms correspond to an SNR 70 injection, while dashed lines and
histograms correspond to SNR 20. The total and tensor SNR are given in green
and orange colors respectively. In pink we denote the SNR of the nontensorial
mode, namely scalar (left) or vector (right).

Another quantity that can be used to characterize the observed signal is the
recovered SNR in each spin weight, defined in Eqs. 14–16. We plot posteriors
for the SNRs of each mode of our simulated chirp signals in Fig. 10 for
injected total SNRs in all modes of 20 (dashed) and 70 (solid), again in an
HLVK detector network. We obtain qualitatively similar results for the burst
injections, so we omit them for clarity.

In both TS (left) and TV (right) cases we can recover the relative tensor-to-
nontensor SNR for loud signals, as demonstrated by the SNR 70 case.
Concretely, if we observed the SNR 70 signal shown on the left (right) panel,
we would be able to correctly infer the ratio of scalar (vector) to tensor
power of the signal as projected onto this specific detector network. For all
SNR 70 simulations, we also find that the SNR posteriors slightly
underestimate the true SNR value—an effect that becomes more pronounced as the
signal SNR decreases, as seen for the SNR 20 cases.

The reason BayesWave underestimates the signal SNR is twofold. Firstly, the
expectation value of the overlap between the true signal and the
reconstruction is $1-{\cal{O}}(\mathrm{D/SNR}^{-2})$ for sufficiently high SNR
values Chatziioannou  _et al._ (2017c) where $D$ is the dimensionality of the
signal model. This means that even though some posterior samples can achieve
overlaps with the true signal reaching up to $\sim 1$, we expect a statistical
spread in the overlap distribution that is inversely proportional to the
square of the SNR; this reduction in the overlap correspondingly affects the
recovered SNR. Secondly, for BayesWave the dimensionality $D$ of the signal
model is not fixed but rather depends on the injected SNR, as more wavelets
are required to fit louder signals with more resolvable structure. As an
outcome, BayesWave will recover on average a smaller fraction of the injected
signal than regular matched-filtering analyses (see, for example, Fig. 4 of
Ghonge  _et al._ (2020)).

Previous studies have shown that BayesWave can recover SNR 20 signals assuming
they are sufficiently short; however, such studies only involved tensor
signals Kanner  _et al._ (2016); Littenberg  _et al._ (2016); Bécsy  _et al._
(2017); Ghonge  _et al._ (2020); Lee  _et al._ (2021). In our case, the
injected signals contain both tensor and non tensor power, so the total SNR of
20 is spread over more wavelets. In other words, an SNR 20 tensor signal is
modeled by fewer wavelets than an SNR 20 signal with a mixed polarization
content. As a result, BayesWave will recover a smaller fraction of the TS or
TV signal power than it would for a T signal of the same SNR and duration. We
find that, for both TS and TV, the total SNR and the per-mode-SNR is
underestimated, suggesting that the analysis does not recover a totality of
the available signal power. Despite that, SNR 20 signals will still offer
clear indications of beyond-GR physics if nontensorial modes are present,
including an accurate estimate of the ratio of tensor to nontensor power in
the signal.

Figure 11: Sky localization posteriors for the signals from Fig. 7. We plot
90%-credible contours for the right ascension ($\alpha$) and sine of the
declination ($\sin\delta$), with black lines denoting the true values. Green
and orange lines show results from analyses on tensor-only data using the
generic and elliptical polarization tensor model respectively. Purple solid
and pink dashed lines show results from the TS analysis on TS data (left
panel) and the TV on TV data (right panel) for the chirp and burst simulations
respectively, as presented in Fig. 7.

Finally, Fig. 11 shows sky location posteriors for various analyses. We do not
assume the source sky location is known, but rather allow BayesWave to infer
it from the data. Extracting the sky location is necessary for converting the
observed SNR distributions from Fig. 10 to the intrinsic SNR of the source,
giving an estimate of the relative strength of each mode as emitted by the GW
source, as discussed below Eqs. 14–16. In the case of a nontensorial mode
detection, such information could be useful to study the properties of the
beyond-GR physics. Figure 11 shows the 90%-credible contours of the sky
posteriors for the same analyses as in Fig. 7. The two tensor posteriors
correspond to purely-tensor data analyzed either with the elliptical
polarization model of Sec. II.0.1 or the generic polarization model of Sec.
II.0.2. The TS (left) and TV (right) posteriors correspond to the analysis of
the mixed-polarization data from Fig. 7, for both the burst and chirp
simulations.

In the left panel, we find that all analyses recover the correct sky posterior
at the 90% level, although the uncertainties are larger for the TS analyses as
expected since it has more degrees of freedom. In both TS and TV cases the
signals with the chirp morphology results in a larger uncertainty over the
recovered sky location than the burst case. This is because in the former case
the tensor and nontensor signals are morphologically more similar, both being
consistent with chirping BBHs. As a result, the uncertainty on the recovered
parameters increases, as also noticed in Fig. 7. Additionally, the generic
tensor analysis results in a less precise sky localization relative to the
elliptical tensor one due to the additional degrees of freedom in Eq. (9). In
the right panel, we find qualitatively similar results; however, the generic
tensor-only analysis results in a small bias in the inferred sky location.
This shift is inherited by the mixed-polarization analyses, although it can be
less pronounced in that case due to the larger uncertainties. Even though we
demonstrate this using the TV panel on the right of Fig. 11, this behavior
arises from the tensor model, not from the TV model (as can be seen from the
green and orange curves).

Although this has no effect on our ability to test GR, we study the bias in
the inferred sky location in more detail and find it to be due to the large
flexibility of the generic (tensor) polarization model compared to the
elliptic (tensor) polarization model. It is therefore unrelated to the
implementation of the beyond-GR analyses presented here: the generic
polarization model of Sec. II.0.2 results in a prior that strongly disfavors
face-on or face-off systems Isi (2021). Therefore, if the system is close to
these two extremes, the inferred orientation could be biased due to the prior
and result in a corresponding bias to the sky location. We indeed study a
large number of simulated signals within GR (i.e. including only and analyzed
with solely tensor modes) and find that the amount of bias depends weakly on
the exact sky location and primarily on the inclination of the source. This is
an example of the effect discussed in Sec. II.0.2 where the extreme
flexibility of the generic polarization model can cause it to underperform
compared to the elliptical mode in cases where the signal is indeed
elliptically polarized. In future work we will explore potential mitigations
of this by using different parametrizations, imposing a prior that does not
disfavor face-on/face-off inclinations, or restricting to the elliptical
polarization model.

##  VI TVS constraints in four- and five-detector networks

Since only HLVK will be available in the near future, so far we have focused
on the study of two-helicity models (TV and TS), which we expect to be less
degenerate than their full, three-helicity counterpart (TVS). This expectation
stems from the fact that, as discussed in Sec. I, full inversion of the
antenna pattern matrix, Eq. (1), would demand access to five independent
detectors. Yet, this does not mean that meaningful statements about TVS cannot
be made with fewer instruments. Perfectly inverting Eq. (1) amounts to fully
categorizing the five polarization modes by measuring five amplitudes (or
SNRs); however, we might not need all this information if all we want is to
distinguish each of T, V and S, without concern for the specific distribution
of signal power _within_ each spin-weight (say, the precise relative amplitude
of plus versus cross, or $v_{1}$ versus $v_{2}$). Put differently, we are
interested in independently inferring T and V, but within -say- T the plus and
cross modes are not completely independent, as no theory predicts solely one
without the other, and similar for $v_{1}$ and $v_{2}$. In BayesWave’s signal
model, this is expressed in the fact that plus and cross (or $v_{1}$ and
$v_{2}$) share some common parameters, but T and V share none. In this section
we show that indeed the HLVK network will be sufficient to place interesting
constraints on TVS, although uncertainties will be reduced once India becomes
operational.

Figure 12: Posterior for the total SNR and the SNR per mode for simulated TS
(top) and TV (bottom) signals of the burst morphology from Fig. 7. Each signal
is injected in a 4-detector HLVK network and analyzed with a full TVS signal
model. Solid vertical lines and histograms of the same color correspond to the
injected and recovered SNR for each mode or combination of modes. The
total/tensor/vector/scalar SNR is given in green/orange/purple/pink
respectively. In the top (bottom) panel the recovered vector (scalar) SNR is
consistent with zero, consistent with the fact that the simulated signal had
no vector (scalar) power.

We revisit the two simulated signals from the left panels of Fig. 7 and
analyze them with the full TVS signal model of Eqs. (9–13); as before, we
simulate a measurement by the HLVK network at design sensitivity, without
assuming the sky location is known a priori. The signals contain only TS (TV)
modes, with the tensor signal corresponding to a merging BBH and the scalar
(vector) part has a burst morphology (see Sec. V). Figure 12 shows the
recovered and injected total SNR, as well as the SNR in each spin-weight as
defined in Eqs. (14–16). The per-spin-weight SNR posteriors reveal how well
the analysis can separate the different polarizations and reconstruct their
relative contribution to the overall signal.

The total SNR posteriors in Fig. 12 agree with the injected value as they did
in Fig. 10 for both the TS and TV signals. More interestingly, in both panels,
the SNR of the “absent” mode (vector for the top and scalar for the bottom
panel) is consistent with zero; this shows that even with only four detectors,
the analysis correctly infers which modes are present in the signal.
Additionally, we find that the vector upper limit is relatively more stringent
than the scalar one. This is reminiscent of Fig. 4, where we argued that it is
easier for BayesWave to constrain vector rather than scalar modes due to the
different number of wavelet parameters involved in each case. Despite the
success of the analysis in constraining the absent mode, we find that the
separation of modes is not ideal; this is revealed, for example, by the vector
and tensor posteriors in the bottom panel, which are respectively over and
under estimated as a consequence of unbroken degeneracies.

Figure 13: Whitened signal reconstructions for a burst injection containing
tensor and vector modes in a $5$-detector network as analyzed by the full TVS
signal model. The solid grey line shows the simulated data in each of the $5$
detectors. Pink shaded regions show the 90% credible interval for the
reconstruction of the full signal, while the scalar (S) and vector (V) part of
the full reconstruction is shown in orange and green bands respectively. As
expected given the simulated signal, the scalar reconstruction vanishes for
all times.

For a robust extraction of all five TVS modes, we turn to the five-detector
HLVKI network. For this purpose, we select the simulated signal from the
bottom panel of Fig. 12, which contains tensor and vector modes, and analyze
it in an HLVKI network with the full TVS model; we keep the SNR in HLVK the
same as before, so that the total network SNR increases due to LIGO India’s
contribution. We plot the 90%-credible interval for the signal reconstruction
in Fig. 13, for different combinations of polarization modes derived from the
full TVS analysis. The full reconstruction containing all TVS modes accurately
traces the injected data. The vector reconstruction is similar to the one in
Fig. 7, which was obtained by a TV analysis. On the other hand, the scalar
reconstruction is identically zero at the quoted credible level, which is what
we expect for an injection that has no scalar component, as was the case here.

Figure 14: Posterior for the total SNR and the SNR per mode for a burst
injection containing tensor and vector modes in a $5$-detector network as
analyzed by the full TVS signal model. Solid vertical lines and histograms of
the same color correspond to the injected and recovered SNR for each mode or
combination of modes. The injected scalar SNR is $0$, so the recovered
posterior is sharply peaked at this value.

The accurate extraction of each of the five polarization modes is also
demonstrated in Fig. 14, which shows the recovered and injected SNRs for
different combinations of polarization modes. Similar to the four-detector
case from Fig. 12, the analysis correctly concludes that the signal contains
no scalar power, recovering a scalar SNR that is consistent with zero. As
expected, however, the five-detector upper limit is more stringent than the
four-detector one. Moreover, the recovered vector, tensor, and total SNR are
now consistent with their injected values. Comparing to the bottom panel of
Fig. 12, we conclude that a five-detector network can more robustly separate
tensor, vector, and scalar modes.

Figure 15: SNR 90% upper limits for GW150914-like signals injected in a
$5$-detector HLVKI network as function of right ascension $\alpha$ and
declination $\delta$. Each signal is analyzed with the full TVS polarization
model that allows for all $5$ relevant polarization modes and we plot the
upper limit on scalar (top) and vector (bottom) polarization modes from this
single analysis. The simulated signals obey GR and have been injected in
separate and independent random Gaussian noise realizations. The upper limits
are more stringent that the ones obtained from a $4$-detector network in Fig.
3.

As a final step, we also consider signals that obey GR in the context of an
HLVKI network and full TVS analyses. We repeat the analysis of Sec. IV by
injecting the same $200$ GW150914-like signals, randomly distributed over the
sky, this time using a five-detector HLVKI network. We analyze the signals
with the full TVS model, and plot the resulting upper limits on the SNR of the
scalar (top) and vector (bottom) modes in Fig. 15. The main difference between
the analysis of Fig. 15 and that of Fig. 3, besides the addition of a
detector, is that now we obtain the SNR upper limits on the scalar and vector
modes from a single TVS analysis, as opposed to separate TS and TV analyses.
The HLVKI network can not only place constraints on all relevant polarization
modes simultaneously, but does so while also producing more stringent upper
limits (compare the color scales in Figs. 15 and 3).

##  VII Conclusions

We have presented an analysis capable of constraining the mixed polarization
content of GW signals. Our analysis is based on BayesWave and models the
tensor and nontensor polarization modes in terms of a sum of sine-Gaussian
wavelets. This flexible signal model makes no assumptions about the amplitude
or phase evolution of the signal, besides that the signal propagates at the
speed of light and does not suffer from dispersion during propagation between
detectors. This makes our analysis ideal for characterizing the morphology of
a nontensor polarization modes, if they are ever detected, as demonstrated
with example signals. In the most likely case of a null result and no detected
deviation from GR, we can place upper limits on the presence of nontensor
polarizations in the observed signals. Our results suggest that with a design-
sensitivity network of four detectors, as anticipated for O4, we will be able
to place constraints and/or characterize the presence of vector and/or scalar
modes in the detected signals, on top of the expected tensorial modes. A five-
detector network, as expected later in the decade, will strengthen constraints
and result in a more accurate separation of all five polarization modes to
which differential-arm detectors are sensitive. We tested our infrastructure
on CBC signals, since they are the only sources of GWs detected to date;
however, the analysis is generic and can handle other types of transient
signals, such as bursts from supernovae or other sources.

Our analysis extends and generalizes previous efforts to constrain the GW
polarization content by allowing for a mixed tensor/nontensor content, not
requiring a known source sky location, and taking advantage of the phase
coherence of the signal. Regarding the first point, results from GW170817,
GW170814, and other events place strong constraints against the complete
absence of tensor modes in the signal from both BNSs and BBHs. Therefore
upcoming analyses of the GW signal polarization content should account for the
presence of tensor modes. Indeed, in our analysis we test for a mixed content
of tensor and nontensor modes. It would be trivial to analyze signals assuming
a pure S or V content with our infrastructure if required.

Furthermore, the signal sky location does not need to be a priori known, but
is a model parameter that is marginalized over together with other model
parameters such as the wavelet number. This does away with the requirement of
an electromagnetic counterpart to the GW signal that could provide the sky
location, something expected to be rare especially for BBHs. Additionally,
using electromagnetic counterparts to identify the sky location of a signal in
order to perform tests of the polarization content contains an element of
circularity. Electromagnetic counterparts to GW signals are detected with
observing campaigns that either use the GW localization directly or compare
triggers to it. The initial GW localization is always computed assuming the
signal is composed solely of tensor modes that agree with GR. Therefore any
counterpart found in the localization area provided will have a sky location
that is by selection consistent with tensor modes. Though this does not
completely preclude the possibility of polarization constraints using these
sky localizations, it does mean that only signals whose potential nontensor
modes do not affect the sky localization enough to hinder a counterpart
identification can meet the analysis criteria, and thus is a strong selection
effect. Our analysis does not require an externally provided sky localization
and can instead be applied to all detected signals.

Finally, our methodology takes advantage of the fact that CBC signals with
nontensor polarization modes are phase-coherent similarly to their tensor
counterparts. We indeed model directly the phase evolution of each
polarization mode, an approach that is more sensitive than analyses that
detect excess power in the data caused by nontensor modes. Excess power
analyses are necessary for signals without phase coherence, such as a
stochastic background signal, but will be suboptimal for signals with a non
random phase evolution.

Given the number of detectors that will be operational in the upcoming O4 run,
we expect to obtain the most stringent constraints in the context of one
additional mode (scalar or vector), i.e., a TS or TV analysis. A full analysis
including all five possible modes with a four-detector network would contain
more degrees of freedom than detectors in the network and would thus not be
able to fully invert Eq. (1). However, this does not mean that a four detector
network cannot be used to make interesting statements about the full TVS
model, since we do not need to break _all_ possible degeneracies between the
five polarization modes to distinguish between the three spin-weights, T, V,
and S. Indeed, we find that HLVK will be able to place meaningful TVS
constraints, even when the sky location is unknown. The addition of LIGO India
to the detector network in O5 Abbott  _et al._ (2018c) will only serve to make
full TVS analysis more informative in the future. We demonstrate that the full
five-detector network will further facilitate such TVS analyses and robust
extraction of all five polarization modes available to differential-arm
detectors, even when the sky location of the source is unknown.

The signal model of Eq. (8) is designed to be as generic as possible, modeling
each of the five polarization modes in terms of sums of sine-Gaussian
wavelets. However, as more signals are detected we will be able to place
increasingly stringent constraints on their properties, and thus potentially
specialize the signal model. As an example, our analysis currently allows for
fully generic plus/cross and $v1/v2$ polarizations, but the majority of
signals observed to date are consistent with being elliptical polarized, since
they originate from nonprecessing compact binaries. One option is, therefore,
to restrict our signal model to elliptically polarized tensor and vector
components along the lines of Eq. (4). This would likely result in better
reconstructions of elliptical signals, as an elliptical model contains fewer
degrees of freedom than the generic polarization one.

Another possibility concerns the fully generic form of the tensor modes.
Indeed all current detections appear to be consistent with GR, we thus know
that GW signals contain tensor modes that resemble those of GR. A more
restricted signal model would then assume that the tensor modes are given by
regular CBC templates with some parametrized deviation from GR (such as the
parametrized post-Einsteinian framework Yunes and Pretorius (2009)), while the
nontensor modes are again modeled with sine-Gaussian wavelets. The possibility
of joint analyses using CBC templates and wavelets with BayesWave was recently
demonstrated in Chatziioannou  _et al._ (2021), so we plan to explore this
option in future work.

As the global detector network continues to grow, studies of the local
geometry of GWs will be significantly enriched by our enhanced ability to
break degeneracies between possible combinations of GW polarizations. The
constraints will also improve in quantity and quality thanks to the increased
number and variety of detections, as well as their greater average SNR. The
method presented here and its extensions will be crucial in fully taking
advantage of this, and will enable precision tests of general relativity with
GW polarizations.

###### Acknowledgements.

We thank Tom Callister for helpful comments on the manuscript. M.I. is
supported by NASA through the NASA Hubble Fellowship grant No. HST-
HF2-51410.001-A awarded by the Space Telescope Science Institute, which is
operated by the Association of Universities for Research in Astronomy, Inc.,
for NASA, under contract NAS5-26555. C.-J.H. acknowledge support of the
National Science Foundation, and the LIGO Laboratory. This research has made
use of data, software and/or web tools obtained from the Gravitational Wave
Open Science Center (https://www.gw-openscience.org), a service of LIGO
Laboratory, the LIGO Scientific Collaboration and the Virgo Collaboration.
This material is based upon work supported by NSF’s LIGO Laboratory which is a
major facility fully funded by the National Science Foundation. Virgo is
funded by the French Centre National de Recherche Scientifique (CNRS), the
Italian Istituto Nazionale della Fisica Nucleare (INFN) and the Dutch Nikhef,
with contributions by Polish and Hungarian institutes. The authors are
grateful for computational resources provided by the LIGO Laboratory and
supported by National Science Foundation Grants PHY-0757058 and PHY-0823459,
and for resources provided by the Open Science Grid Pordes  _et al._ (2007);
Sfiligoi  _et al._ (2009), which is supported by the National Science
Foundation award 1148698, and the U.S. Department of Energy’s Office of
Science. Software: matplotlib Hunter (2007), numpy Harris  _et al._ (2020),
scipy Virtanen  _et al._ (2020), gwpy Macleod  _et al._ (2020), seaborn Waskom
(2021).

## References

  * Aasi  _et al._ (2015) J. Aasi _et al._ (LIGO Scientific), Class. Quant. Grav. 32, 074001 (2015), arXiv:1411.4547 [gr-qc] .
  * Acernese  _et al._ (2015) F. Acernese _et al._ (Virgo Collaboration), Class. Quant. Grav. 32, 024001 (2015), arXiv:1408.3978 [gr-qc] .
  * Abbott  _et al._ (2016a) B. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 116, 221101 (2016a), [Erratum: Phys.Rev.Lett. 121, 129902 (2018)], arXiv:1602.03841 [gr-qc] .
  * Abbott  _et al._ (2016b) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. X6, 041015 (2016b), [erratum: Phys. Rev.X8,no.3,039903(2018)], arXiv:1606.04856 [gr-qc] .
  * Yunes  _et al._ (2016) N. Yunes, K. Yagi, and F. Pretorius, Phys. Rev. D 94, 084002 (2016), arXiv:1603.08955 [gr-qc] .
  * Abbott  _et al._ (2017a) B. P. Abbott _et al._ (LIGO Scientific and Virgo Collaboration), Phys. Rev. Lett. 118, 221101 (2017a), [Erratum: Phys. Rev. Lett. 121, 129901(E) (2018)], arXiv:1706.01812 [gr-qc] .
  * Ezquiaga and Zumalacárregui (2017) J. M. Ezquiaga and M. Zumalacárregui, Phys. Rev. Lett. 119, 251304 (2017), arXiv:1710.05901 [astro-ph.CO] .
  * Sakstein and Jain (2017) J. Sakstein and B. Jain, Phys. Rev. Lett. 119, 251303 (2017), arXiv:1710.05893 [astro-ph.CO] .
  * Creminelli and Vernizzi (2017) P. Creminelli and F. Vernizzi, Phys. Rev. Lett. 119, 251302 (2017), arXiv:1710.05877 [astro-ph.CO] .
  * Baker  _et al._ (2017) T. Baker, E. Bellini, P. G. Ferreira, M. Lagos, J. Noller, and I. Sawicki, Phys. Rev. Lett. 119, 251301 (2017), arXiv:1710.06394 [astro-ph.CO] .
  * Abbott  _et al._ (2017b) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 119, 141101 (2017b), arXiv:1709.09660 [gr-qc] .
  * Abbott  _et al._ (2019a) B. P. Abbott _et al._ (LIGO Scientific and Virgo Collaboration), Phys. Rev. Lett. 123, 011102 (2019a), arXiv:1811.00364 [gr-qc] .
  * Abbott  _et al._ (2019b) B. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. D 100, 104036 (2019b), arXiv:1903.04467 [gr-qc] .
  * Isi  _et al._ (2019) M. Isi, K. Chatziioannou, and W. M. Farr, Phys. Rev. Lett. 123, 121101 (2019), arXiv:1904.08011 [gr-qc] .
  * Abbott  _et al._ (2020a) R. Abbott _et al._ (LIGO Scientific, Virgo), “Tests of General Relativity with Binary Black Holes from the second LIGO-Virgo Gravitational-Wave Transient Catalog,” (2020a), arXiv:2010.14529 [gr-qc] .
  * Akutsu  _et al._ (2020) T. Akutsu _et al._ (KAGRA), (2020), arXiv:2005.05574 [physics.ins-det] .
  * Iyer  _et al._ (2011) B. Iyer _et al._ ,  _LIGO-India, Proposal of the Consortium for Indian Initiative in Gravitational-wave Observations (IndIGO)_ , Tech. Rep. LIGO-M1100296 (2011) https://dcc.ligo.org/LIGO-M1100296/public.
  * Will (2014) C. M. Will, Living Rev. Rel. 17, 4 (2014), arXiv:1403.7377 [gr-qc] .
  * Eardley  _et al._ (1973a) D. M. Eardley, D. L. Lee, A. P. Lightman, R. V. Wagoner, and C. M. Will, Phys. Rev. Lett. 30, 884 (1973a).
  * Eardley  _et al._ (1973b) D. M. Eardley, D. L. Lee, and A. P. Lightman, Phys. Rev. D 8, 3308 (1973b).
  * Isi and Weinstein (2017) M. Isi and A. J. Weinstein, “Probing gravitational wave polarizations with signals from compact binary coalescences,” (2017), arXiv:1710.03794 [gr-qc] .
  * Isi  _et al._ (2015) M. Isi, A. J. Weinstein, C. Mead, and M. Pitkin, Phys. Rev. D 91, 082002 (2015), arXiv:1502.00333 [gr-qc] .
  * Isi  _et al._ (2017) M. Isi, M. Pitkin, and A. J. Weinstein, Phys. Rev. D 96, 042001 (2017), arXiv:1703.07530 [gr-qc] .
  * Callister  _et al._ (2017) T. Callister, A. S. Biscoveanu, N. Christensen, M. Isi, A. Matas, O. Minazzoli, T. Regimbau, M. Sakellariadou, J. Tasson, and E. Thrane, Phys. Rev. X 7, 041058 (2017), arXiv:1704.08373 [gr-qc] .
  * Abbott  _et al._ (2018a) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 120, 031104 (2018a), arXiv:1709.09203 [gr-qc] .
  * Abbott  _et al._ (2018b) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 120, 201102 (2018b), arXiv:1802.10194 [gr-qc] .
  * Takeda  _et al._ (2018) H. Takeda, A. Nishizawa, Y. Michimura, K. Nagano, K. Komori, M. Ando, and K. Hayama, Phys. Rev. D 98, 022008 (2018), arXiv:1806.02182 [gr-qc] .
  * Takeda  _et al._ (2021) H. Takeda, S. Morisaki, and A. Nishizawa, Phys. Rev. D 103, 064037 (2021), arXiv:2010.14538 [gr-qc] .
  * Haster (2020) C.-J. Haster, Res. Notes AAS 4, 209 (2020).
  * Chatziioannou  _et al._ (2012) K. Chatziioannou, N. Yunes, and N. Cornish, Phys. Rev. D 86, 022004 (2012), [Erratum: Phys.Rev.D 95, 129901 (2017)], arXiv:1204.2585 [gr-qc] .
  * Hiroki  _et al._ (2021) T. Hiroki, S. Morisaki, and A. Nishizawa, (2021), arXiv:2105.00253 [gr-qc] .
  * Yunes and Pretorius (2009) N. Yunes and F. Pretorius, Phys. Rev. D 80, 122003 (2009), arXiv:0909.3328 [gr-qc] .
  * Hayama and Nishizawa (2013) K. Hayama and A. Nishizawa, Phys. Rev. D 87, 062003 (2013), arXiv:1208.4596 [gr-qc] .
  * Guersel and Tinto (1989) Y. Guersel and M. Tinto, Phys. Rev. D 40, 3884 (1989).
  * Hagihara  _et al._ (2019) Y. Hagihara, N. Era, D. Iikawa, A. Nishizawa, and H. Asada, Phys. Rev. D 100, 064010 (2019), arXiv:1904.02300 [gr-qc] .
  * Pang  _et al._ (2020) P. T. H. Pang, R. K. L. Lo, I. C. F. Wong, T. G. F. Li, and C. Van Den Broeck, Phys. Rev. D 101, 104055 (2020), arXiv:2003.07375 [gr-qc] .
  * Wong  _et al._ (2020) I. C. F. Wong, P. T. H. Pang, and R. K. L. Lo,  _Technical Notes on Null Stream Polarization Test_, Tech. Rep. LIGO-T2000405 (LIGO Scientific Collaboration, 2020).
  * Sutton  _et al._ (2010) P. J. Sutton _et al._ , New J. Phys. 12, 053034 (2010), arXiv:0908.3665 [gr-qc] .
  * Cornish and Littenberg (2015) N. J. Cornish and T. B. Littenberg, Class. Quant. Grav. 32, 135012 (2015), arXiv:1410.3835 [gr-qc] .
  * Cornish  _et al._ (2021) N. J. Cornish, T. B. Littenberg, B. Bécsy, K. Chatziioannou, J. A. Clark, S. Ghonge, and M. Millhouse, Phys. Rev. D 103, 044006 (2021), arXiv:2011.09494 [gr-qc] .
  * Abbott  _et al._ (2018c) B. P. Abbott _et al._ (KAGRA, LIGO Scientific, VIRGO), Living Rev. Rel. 21, 3 (2018c), arXiv:1304.0670 [gr-qc] .
  * Abbott  _et al._ (2017c) B. P. Abbott _et al._ (LIGO Scientific, Virgo, Fermi-GBM, INTEGRAL), Astrophys. J. Lett. 848, L13 (2017c), arXiv:1710.05834 [astro-ph.HE] .
  * Will (2018) C. M. Will,  _Theory and Experiment in Gravitational Physics_, 2nd ed. (Cambridge University Press, Cambridge, 2018).
  * Ghonge  _et al._ (2020) S. Ghonge, K. Chatziioannou, J. A. Clark, T. Littenberg, M. Millhouse, L. Cadonati, and N. Cornish, Phys. Rev. D 102, 064056 (2020), arXiv:2003.09456 [gr-qc] .
  * Chatziioannou  _et al._ (2017a) K. Chatziioannou, J. A. Clark, A. Bauswein, M. Millhouse, T. B. Littenberg, and N. Cornish, Phys. Rev. D 96, 124035 (2017a), arXiv:1711.00040 [gr-qc] .
  * Torres-Rivas  _et al._ (2019) A. Torres-Rivas, K. Chatziioannou, A. Bauswein, and J. A. Clark, Phys. Rev. D 99, 044014 (2019), arXiv:1811.08931 [gr-qc] .
  * Dálya  _et al._ (2021) G. Dálya, P. Raffai, and B. Bécsy, Class. Quant. Grav. 38, 065002 (2021), arXiv:2006.06256 [astro-ph.HE] .
  * Millhouse  _et al._ (2018) M. Millhouse, N. J. Cornish, and T. Littenberg, Phys. Rev. D97, 104057 (2018), arXiv:1804.03239 [gr-qc] .
  * Bécsy  _et al._ (2017) B. Bécsy, P. Raffai, N. J. Cornish, R. Essick, J. Kanner, E. Katsavounidis, T. B. Littenberg, M. Millhouse, and S. Vitale, Astrophys. J. 839, 15 (2017), arXiv:1612.02003 [astro-ph.HE] .
  * Abbott  _et al._ (2020b) R. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 125, 101102 (2020b), arXiv:2009.01075 [gr-qc] .
  * Abbott  _et al._ (2020c) R. Abbott _et al._ (LIGO Scientific, Virgo), Astrophys. J. Lett. 900, L13 (2020c), arXiv:2009.01190 [astro-ph.HE] .
  * Graham  _et al._ (2020) M. J. Graham _et al._ , Phys. Rev. Lett. 124, 251102 (2020), arXiv:2006.14122 [astro-ph.HE] .
  * Ashton  _et al._ (2020) G. Ashton, K. Ackley, I. M. Hernandez, and B. Piotrzkowski, (2020), arXiv:2009.12346 [astro-ph.HE] .
  * Palmese  _et al._ (2021) A. Palmese, M. Fishbach, C. J. Burke, J. T. Annis, and X. Liu, (2021), arXiv:2103.16069 [astro-ph.HE] .
  * Abbott  _et al._ (2019c) R. Abbott _et al._ (LIGO Scientific, Virgo), (2019c), arXiv:1912.11716 [gr-qc] .
  * Littenberg and Cornish (2015) T. B. Littenberg and N. J. Cornish, Phys. Rev. D91, 084034 (2015), arXiv:1410.3852 [gr-qc] .
  * Chatziioannou  _et al._ (2019) K. Chatziioannou, C.-J. Haster, T. B. Littenberg, W. M. Farr, S. Ghonge, M. Millhouse, J. A. Clark, and N. Cornish, Phys. Rev. D 100, 104004 (2019), arXiv:1907.06540 [gr-qc] .
  * Chen  _et al._ (2020) H.-Y. Chen, C.-J. Haster, S. Vitale, W. M. Farr, and M. Isi, (2020), arXiv:2009.14057 [astro-ph.CO] .
  * Evans  _et al._ (2020) M. Evans, R. Sturani, S. Vitale, and E. Hall,  _Unofficial sensitivity curves (ASD) for aLIGO, Kagra, Virgo, Voyager, Cosmic Explorer, and Einstein Telescope_, Tech. Rep. LIGO-T1500293 (2020).
  * Pratten  _et al._ (2020a) G. Pratten _et al._ , (2020a), arXiv:2004.06503 [gr-qc] .
  * Pratten  _et al._ (2020b) G. Pratten, S. Husa, C. García-Quirós, M. Colleoni, A. Ramos-Buades, H. Estellés, and R. Jaume, Phys. Rev. D 102, 064001 (2020b), arXiv:2001.11412 [gr-qc] .
  * García-Quirós  _et al._ (2020) C. García-Quirós, M. Colleoni, S. Husa, H. Estellés, G. Pratten, A. Ramos-Buades, M. Mateu-Lucena, and R. Jaume, Phys. Rev. D 102, 064002 (2020), arXiv:2001.10914 [gr-qc] .
  * García-Quirós  _et al._ (2021) C. García-Quirós, S. Husa, M. Mateu-Lucena, and A. Borchers, Class. Quant. Grav. 38, 015006 (2021), arXiv:2001.10897 [gr-qc] .
  * Abbott  _et al._ (2016c) B. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 116, 061102 (2016c), arXiv:1602.03837 [gr-qc] .
  * Abbott  _et al._ (2016d) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 116, 241103 (2016d), arXiv:1606.04855 [gr-qc] .
  * Fairhurst (2011) S. Fairhurst, Class. Quant. Grav. 28, 105021 (2011), arXiv:1010.6192 [gr-qc] .
  * Pankow  _et al._ (2018) C. Pankow, E. A. Chase, S. Coughlin, M. Zevin, and V. Kalogera, Astrophys. J. Lett. 854, L25 (2018), arXiv:1801.02674 [astro-ph.HE] .
  * Pankow  _et al._ (2020) C. Pankow, M. Rizzo, K. Rao, C. P. L. Berry, and V. Kalogera, Astrophys. J. 902, 71 (2020), arXiv:1909.12961 [astro-ph.HE] .
  * Blackman  _et al._ (2017) J. Blackman, S. E. Field, M. A. Scheel, C. R. Galley, C. D. Ott, M. Boyle, L. E. Kidder, H. P. Pfeiffer, and B. Szilágyi, Phys. Rev. D 96, 024058 (2017), arXiv:1705.07089 [gr-qc] .
  * Okounkova  _et al._ (2019) M. Okounkova, L. C. Stein, M. A. Scheel, and S. A. Teukolsky, Phys. Rev. D 100, 104026 (2019), arXiv:1906.08789 [gr-qc] .
  * Okounkova  _et al._ (2020) M. Okounkova, L. C. Stein, J. Moxon, M. A. Scheel, and S. A. Teukolsky, Phys. Rev. D 101, 104016 (2020), arXiv:1911.02588 [gr-qc] .
  * Okounkova (2020) M. Okounkova, Phys. Rev. D 102, 084046 (2020), arXiv:2001.03571 [gr-qc] .
  * Witek  _et al._ (2020) H. Witek, L. Gualtieri, and P. Pani, Phys. Rev. D 101, 124055 (2020), arXiv:2004.00009 [gr-qc] .
  * Apostolatos  _et al._ (1994) T. A. Apostolatos, C. Cutler, G. J. Sussman, and K. S. Thorne, Phys.Rev. D49, 6274 (1994).
  * Schmidt  _et al._ (2012) P. Schmidt, M. Hannam, and S. Husa, Phys. Rev. D 86, 104063 (2012), arXiv:1207.3088 [gr-qc] .
  * Chatziioannou  _et al._ (2013) K. Chatziioannou, A. Klein, N. Yunes, and N. Cornish, Phys. Rev. D 88, 063011 (2013), arXiv:1307.4418 [gr-qc] .
  * Schmidt  _et al._ (2015) P. Schmidt, F. Ohme, and M. Hannam, Phys. Rev. D91, 024043 (2015), arXiv:1408.1810 [gr-qc] .
  * Chatziioannou  _et al._ (2017b) K. Chatziioannou, A. Klein, N. Cornish, and N. Yunes, Phys. Rev. Lett. 118, 051101 (2017b), arXiv:1606.03117 [gr-qc] .
  * Fairhurst  _et al._ (2019) S. Fairhurst, R. Green, C. Hoy, M. Hannam, and A. Muir, (2019), arXiv:1908.05707 [gr-qc] .
  * Chatziioannou  _et al._ (2017c) K. Chatziioannou, A. Klein, N. Yunes, and N. Cornish, Phys. Rev. D 95, 104004 (2017c), arXiv:1703.03967 [gr-qc] .
  * Kanner  _et al._ (2016) J. B. Kanner, T. B. Littenberg, N. Cornish, M. Millhouse, E. Xhakaj, F. Salemi, M. Drago, G. Vedovato, and S. Klimenko, Phys. Rev. D 93, 022002 (2016).
  * Littenberg  _et al._ (2016) T. B. Littenberg, J. B. Kanner, N. J. Cornish, and M. Millhouse, Phys. Rev. D 94, 044050 (2016).
  * Bécsy  _et al._ (2017) B. Bécsy, P. Raffai, N. J. Cornish, R. Essick, J. Kanner, E. Katsavounidis, T. B. Littenberg, M. Millhouse, and S. Vitale, Astrophys. J. 839, 15 (2017), arXiv:1612.02003 [astro-ph.HE] .
  * Lee  _et al._ (2021) Y. S. C. Lee, M. Millhouse, and A. Melatos, Phys. Rev. D 103, 062002 (2021), arXiv:2102.10816 [gr-qc] .
  * Isi (2021) M. Isi, “Parametrizing gravitational-wave polarizations,”  (2021), LIGO-P2100158.
  * Chatziioannou  _et al._ (2021) K. Chatziioannou, N. Cornish, M. Wijngaarden, and T. B. Littenberg, Phys. Rev. D 103, 044013 (2021), arXiv:2101.01200 [gr-qc] .
  * Pordes  _et al._ (2007) R. Pordes _et al._ , J. Phys. Conf. Ser. 78, 012057 (2007).
  * Sfiligoi  _et al._ (2009) I. Sfiligoi, D. C. Bradley, B. Holzman, P. Mhashilkar, S. Padhi, and F. Wurthwrin, WRI World Congress 2, 428 (2009).
  * Hunter (2007) J. D. Hunter, Comput. Sci. Eng. 9, 90 (2007).
  * Harris  _et al._ (2020) C. R. Harris _et al._ , Nature 585, 357 (2020).
  * Virtanen  _et al._ (2020) P. Virtanen _et al._ , Nature Meth. 17, 261 (2020), arXiv:1907.10121 [cs.MS] .
  * Macleod  _et al._ (2020) D. Macleod _et al._ , “gwpy,”  (2020).
  * Waskom (2021) M. L. Waskom, Journal of Open Source Software 6, 3021 (2021).