Waveform systematics in gravitational-wave inference of signals from binary neutron star merger models incorporating higher order modes information
A. B. Yelikar R. O' Shaughnessy J. Lange A. Z. Jan



















Accurate information from gravitational wave signals from coalescing binary neutron stars provides essential input to downstream
interpretations, including inference of the neutron star population and equation of state. 
However, even adopting the currently most accurate and physically motivated models available for parameter estimation
(PE) of BNSs, these models remain subject to waveform modeling uncertainty: differences between these models may
introduce biases in recovered source properties.  In this work, we describe injection studies investigating these systematic differences between the two best waveform models available for BNS currently, and . We demonstrate that for BNS sources observable by current second-generation detectors, differences for low-amplitude signals are significant for certain sources. 





Waveform systematics in gravitational-wave inference of signals from binary neutron star merger models incorporating higher order modes information
    A. Z. Jan
    April 29, 2024
===================================================================================================================================================






§ INTRODUCTION










Since the discovery of gravitational waves from  GW150914 <cit.>, the Advanced Laser Interferometer Gravitational-Wave Observatory (LIGO) 
 <cit.> and  Virgo
<cit.>
detectors continue to discover gravitational waves (GW) from coalescing binary black holes (BBHs) and neutron stars.
The properties of each source are inferred by comparing each observation to some estimate(s), commonly called an approximant, 
for the GWs emitted when a BBH merges.  
As illustrated recently with GW190521 <cit.>,
GW190814 <cit.>, GW190412 <cit.>, and the discussion in GWTC-3 <cit.>, these
approximations have enough differences with respect to to each other to produce noticeable differences in inferred posterior distributions, consistent with prior work 
<cit.>.

Despite ongoing generation of new waveforms with increased accuracy <cit.>, these previous investigations suggest that waveform model systematics
can remain a limiting factor in inferences about individual events <cit.> and populations 
<cit.>.

Waveform systematics could be particularly pernicious for detailed analyses to infer the nuclear equation of
state from GW observations. For analyses not involving postmerger physics, these approaches look for the subtle impact
of matter on the pre-merger inspiral radiation, due to tidal deformations and altered inspiral rate <cit.>. 
Even though the GW signal from the early inspiral is well understood because tidal effects are small and accumulate only at the very end of the
inspiral, they're embedded deep within the most challenging strong field component of the GW signal.  

One known limitation of most previous investigations of waveform systematics for BNS is the neglect of  higher-order
modes (HOM). Current state-of-the-art BNS models  <cit.> and  <cit.> incorporate higher-order-modes enabling the exploration of these effects.  



For example, GW190412 <cit.>, which was a merger of two black holes that were highly asymmetric in 
masses, 30M_⊙ and 8M_⊙, demonstrated the existence and importance of HOM in parameter inference of GW from binary mergers <cit.>.  

Using models that incorporate HOM can significantly impact the
inferred parameters of sources identified with current-generation instruments for GW170817 or GW170817-like signals, as demonstrated in  <cit.>. Despite their expected significance to parameter inference, most studies of BNS systematics omit them and rarely 
perform large-scale parameter inference studies to fully assess the impact of systematics, although a similar study was done in  <cit.>.

For example, several mismatch studies are mostly done
for models having only leading-order (2,2) mode <cit.>. 
A study done with fiducial BNS signals with HOMs 
argued that biases in inferring the reduced tidal parameter could be larger than the statistical 90% only for very high 
SNR signals ∼80 <cit.> in the LIGO-Virgo band. 
Recent work by Narikawa <cit.> looked at the effects of multipoles by comparing MultipoleTidal model to PNTidal
and NRTidalv2 waveforms, showing that mismatches and phases do differ between them for systems with higher mass and 
large tidal deformabilities. 





This paper is organized as follows.
In Section <ref>, we review the use of RIFT for parameter inference; the waveform models used in this
work; and the techniques used in <cit.> to assess systematic error. 
We describe one fiducial ensembles of synthetic sources,  targeted at the most common  (low) amplitude sources.
In Section <ref>, we use two well-studied waveform models to demonstrate the impact of contemporary model
systematics. We show that model systematics will be
important, at a level which must impact population results and consistency tests like PP plots.   
In Section <ref>, we summarize our results and discuss their potential applications to future GW sources and
population inference.





§ METHODS






 §.§ RIFT review



A coalescing compact binary in a quasi-circular orbit can be completely characterized by its intrinsic
and extrinsic parameters.  By intrinsic parameters, we refer to the binary's detector-frame masses m_i, spins χ_i, and any quantities
characterizing matter in the system, Λ_i. By extrinsic parameters, we refer to the seven numbers needed to 
characterize its spacetime location and orientation; luminosity distance (d_L), right ascension (α), declination (δ), inclination (ι), polarization (ψ), coalescence phase (ϕ_c), and time (t_c).
  We will express masses in solar mass units and
 dimensionless nonprecessing spins in terms of cartesian components aligned with the orbital angular momentum
 χ_i,z, as we use waveform models that do not account for precession. We will use λ, θ to
refer to intrinsic and extrinsic parameters, respectively. 


    λ : (ℳ,q,χ_1,z,χ_2,z,Λ_1,Λ_2)


    θ : (d_L,α,δ,ι,ψ,ϕ_c,t_c)



RIFT <cit.>
consists of a two-stage iterative process to interpret gravitational wave data d via comparison to
predicted gravitational wave signals h(λ, θ). In one stage, for each  λ_β from some proposed
“grid” β=1,2,… N of candidate parameters, RIFT computes a marginal likelihood 

    L_ marg≡∫ L(λ ,θ )p(θ )dθ

from the likelihood L(λ ,θ ) of the gravitational wave signal in the multi-detector network,
accounting for detector response; see the RIFT paper for a more detailed specification <cit.>.  

In the second stage,  RIFT performs two tasks. First, it generates an approximation to L(λ) based on its
accumulated archived knowledge of marginal likelihood evaluations 
(λ_β, L_β). This approximation can be generated by Gaussian processes, random forests, or other
suitable approximation techniques. Second, using this approximation, it generates the (detector-frame) posterior distribution

    p_ post= L_ marg(λ )p(λ)/∫ dλ L_ marg(λ ) p(λ ).

where prior p(λ) is prior on intrinsic parameters like mass and spin.    The posterior is produced by
performing a Monte Carlo integral:  the evaluation points and weights in that integral are weighted posterior samples,
which are fairly resampled to generate conventional independent, identically distributed “posterior samples.”
For further details on RIFT's technical underpinnings and performance,   see
<cit.>.


  




 §.§ Waveform models



The tidal waveform models used in this study are , ,
and . NRTidalv2 models <cit.> 
are improved versions of NRTidal <cit.> models, which 
are closed-form tidal approximants for binary neutron star coalescence and have been analytically added to selected 
binary black hole GW model to obtain a binary neutron star waveform, either in the time or in the frequency domain. 
The  <cit.> tidal model is based on the binary black hole hybrid model NRHybSur3dq8, 
which is constructed via an interpolation of NR waveforms. It includes all modes ℓ≤ 4, (5,±5) but not (4,±1) and (4,0) 
and models tidal effect up to Λ_1,2 < 5000. This model combines the accuracy of surrogate waveforms with the efficiency of 
PN models.  <cit.>  is another but unique time-domain EOB formalism that includes tidal effects for all modes ℓ≤ 4, but no 𝑚=0 and 
models tidal effect up to Λ_1,2 < 5000 and for spins up to 0.5.





 §.§ Fiducial synthetic sources and PP tests






We consider one universe of 100 synthetic signals for a
3-detector network (HLV), with masses
drawn uniformly in m_i in the region bounded by /M_⊙∈ [1.2,1.4], η∈ [0.2, 0.25] and Λ 
for each object uniformly distributed up to 1000. These bounds are expressed in terms of =(m_1m_2)^3/5/(m_1+m_2)^1/5 
and η=m_1m_2/(m_1+m_2)^2, and encompass the detector-frame parameters of neutron stars observed till date 
<cit.>. 
The extrinsic parameters are drawn uniformly in sky position and isotropically in Euler angles, with source luminosity
distances drawn proportional to d_L^2  between 90Mpc and 240Mpc for low SNR injections.

All our sources have non-precessing spins, with each component assumed to be uniform between [-0.05,0.05], 
this is due to limitations of the  model. 


For complete reproducibility, we use  , and , starting the 
signal evolution and likelihood integration at 30Hz, performing all analysis with 4096Hz time series in Gaussian noise with
known advanced LIGO design PSDs <cit.>.  The BNS signal is generated for 
300 seconds but analysis was performed only on 128 seconds of data. 
For each synthetic event and interferometer, we use the same noise realization for all waveform 
approximations. Therefore, the differences between them arise solely due to waveform systematics.
The  model is utilized with two settings: a) ℓ=5 and b) ℓ=2, which includes only the dominant quadrupole mode. 
and approximants are used with ℓ=4 and ℓ=2 settings respectively. 



Fig.  <ref> shows the cumulative SNR distribution (under a “zero-noise" assumption) 
of the specific synthetic population generated from this distribution. Compared to GW170817's confident detection, 
which was a BNS merger that occurred at 40Mpc detected by LIGO-Virgo with a SNR of 32.4, 
the majority of the signals in this fiducial population have SNRs below or near the typical detection criteria for a BNS merger, 
with some having high enough amplitudes. 

By using a very modest-amplitude population to assess the impact of
waveform systematics, we demonstrate their immediate impact on the kinds of analyses currently being performed on real
observations, let alone future studies.

One way to assess the performance of parameter inference is a probability-probability plots (usually denoted PP-plot) <cit.>.
  Using  on each source k, with true parameters λ_k, we estimate
the fraction of the posterior distributions which is below the true source value λ_k,β   [P̂_k,β(<λ_k,β)] 
for each intrinsic parameter β.  After reindexing the sources so P̂_k,β(λ_k,β) increases with k for some
fixed β, the panels of Figure <ref> for example, show a
plot of k/N versus P̂_k(λ_k,β) for all binary parameters for different scenarios.














 §.§ JS test



To more sharply identify subtle differences introduced by waveform systematics, we will directly compare pairs of inferred posterior probability distributions
deduced with different waveforms but from the same set of data to each other.   Many pairwise error diagnostics have
been used in the literature in general and with RIFT, in particular, <cit.>.   In this study,
motivated by previous work <cit.>, we use the one-dimensional Jensen-Shannon (JS) divergence
J(p,q) = (D(q|m) + D(p|m)) where D(a|b) = ∫ dx a(x) log_2 a(x)/b(x)  and m=(a+b)/2.  The JS divergence is symmetric, ranges between 0 (for
identical distributions) and 1.   For the multidimensional problems described
here, we adopt the median JS divergence over all parameters.

Analyses of O3 using multiple waveforms suggest that binary black holes analyzed with different contemporary waveforms
will produce answers differing by O(0.02)  <cit.>.




§ RESULTS



Our investigations corroborate our central expectation: inferences computed with different waveforms are frequently
substantially different, even for BNS and even for near-threshold events.  The most extreme contrast appeared between
TEOBResumS and other waveform models, where for our near-threshold synthetic events we found ubiquitous qualitatively different
inferences.




 §.§ Anecdotal examples

As an illustrative example of the systematics explored more comprehensively in the population studies below, Figure <ref>
shows the results of parameter inference using multiple recovery waveforms applied to the same synthetic data source,
here a low-amplitude NRHybSur3dq8Tidal-lmax5  injection.   

Despite its low signal amplitude, this example shows that posterior distributions derived from the same synthetic data
will differ, depending on the GW signal model used to interpret it.





 §.§ JS divergences: Demonstrating and quantifying waveform systematics





Figure <ref> shows the cumulative distribution of combined JS divergences of parameters , q, χ_eff 
and Λ, between analyses performed with
NRHybSur3dq8Tidal ℓ_max=5 low amplitude injections.   Specifically, the JS divergence is calculated between an
analysis performed using precisely the same model used for injections on the one hand, and the alternative model listed
in the legend on the other.   Inferences performed with all state-of-the-art models that
include tidal physics often produce qualitatively similar inferences, with JS divergences typically less than 10^-2.
Some relatively modest disagreement expected  between (a) different waveform models and (b) the expected modest impact of
higher-order modes for low-mass sources. 
By contrast,  more than 10% of inferences have JS divergences larger than 10^-1 (mean over all parameters), in all
cases when using models that include similar mode content (but a different waveform model). These calculations suggest that frequently, both waveform
systematics and higher order modes produce noticably different results.




The green line in Figure <ref> shows that tides must be included to avoid significantly biasing the interpretation of a
typical low-amplitude source.  This JS divergence corresponds to inferences that neglect tides entirely (via a
point-particle IMRPhenomD model), even though the true full model includes tides. In this case, the JS divergence is
frequently larger than 0.1, indicating substantial disagreement with the best possible interpretation. 





To investigate the impact of  waveform systematics specifically,  we computed
mismatches between the injected and inference waveforms' (2,2) modes. Mismatch is a simple inner-product-based estimate of waveform similarity between two model predictions <cit.>
h_1(λ) and h_2(λ) at identical model parameters λ:

    M(λ) = 1 - max_t_c,ϕ_c|h_1e^i(2π f t_c +ϕ_c) h_2|/|h_1||h_2|

In this expression, the inner product
⟨ a|b⟩_k≡∫_-∞^∞2dfã(f)^*b̃(f)/S_h,k(|f|) is  implied by the k^th detector's
noise power spectrum S_h,k(f), which for the purposes of waveform similarity is assumed to be the
advanced LIGO instrument, H1. In practice, we adopt a low-frequency cutoff f_ min so all inner products are modified to

    ⟨ a|b⟩_k≡ 2 ∫_|f|>f_ mindf[ã(f)]^*b̃(f)/S_h,k(|f|).


The left panel of Figure
<ref> shows the results for injections recovered with , with the mismatch shown as a
color scale on top of the injected source SNR and cumulative JS divergence (summed over four one-dimensional JS
divergences for parameters , q, χ_eff 
and Λ).  While the mismatches are within the waveform accuracy requirements <cit.> for most of the injections(>10^-2), higher mismatches don't correlate well with extreme JS divergences.  Rather, below some modest SNR, the random noise realization seems to interact adversely with these large mismatches to produce nearly unconstrained posteriors, such that the similarity
between inferences becomes stochastic and diverges at low amplitude. 
The right panel of Figure
<ref> shows qualitatively similar behavior, using comparisons of and .








 §.§ PP plots

The differences between waveforms are significant enough that their imprint can even impact bulk diagnostics such as a PP
plots  <cit.>, which average the impact of waveform systematics over a large population of randomly chosen events.


Figures  <ref> and  <ref> provide another representation of the analyses presented above in the context of JS
divergence: synthetic sources generated with the  ℓ_max=5 and (ℓ_max=4) models. In each panel, colored dots show the empirical cumulative distribution of the posterior quantiles of the injections – the PP plot for each parameter, with colors corresponding to the parameters indicated in the legend. Fig.  <ref> in which the same model was used for both injection
and recovery for a particular panel, we see PP plots for every parameter are consistent with  P(<p)=p, as expected. However, Fig.  <ref> where analyses used for recovery, shows that omitting tidal physics entirely can bring in distinct inconsistencies with P=p for injections where tides are significant and important. 












§ CONCLUSIONS



In this paper, we demonstrated that the interpretation of typical low-amplitude BNS sources will frequently exhibit noteworthy
differences, depending on the adopted model for analysis. Specifically, we showed that a JS divergence between
inferences constructed between two different state-of-the-art waveforms would be larger than 10^-1 for a significant population of mergers.
These large differences persist even though the mismatch between the dominant (2,± 2) mode of this state-of-the-art
waveforms are small, and even though the SNR of our test sources is low.

Additionally, corroborating previous work  <cit.>, we demonstrate that tidal effects are
essential to include even in interpreting our
population of modest-SNR sources.  Specifically, we showed that neglecting tidal physics in parameter inference causes a PP plot to
deviate significantly from the expected diagonal behavior, indicating a biased recovery of mass and/or spin parameters. 

Our study stands in contrast with the expectations of several previous studies, which have argued that
the effects of waveform systematics for these low-mass, low-amplitude sources will be small.   For example, investigations done in  <cit.> suggest that systematic differences will supersede statistical differences
for sources only for high SNR for the current GW detectors. 


Our investigation only demonstrated notable differences in the conclusions derived from different waveform models.
Further study is required to assess to what extent these differences propagate into conclusions derived from a
population of sources or if they average out over the population.




The authors would like to thank Carl-Johan Haster for helpful comments.
ROS  gratefully acknowledges support from NSF awards NSF PHY-1912632, PHY-2012057, and AST-1909534.
AY acknowledges support from NSF PHY-2012057 grant. This material is based upon work supported by NSF’s LIGO Laboratory 
    which is a major facility fully funded by the
    National Science Foundation.
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.
LIGO Laboratory and Advanced LIGO are funded by the 
    United States National Science Foundation (NSF) as
    well as the Science and Technology Facilities Council (STFC) 
    of the United Kingdom,
    the Max-Planck-Society (MPS), 
    and the State of Niedersachsen/Germany 
    for support of the construction of Advanced LIGO 
    and construction and operation of the GEO600 detector.
Additional support for Advanced LIGO was provided by the 
    Australian Research Council.
Virgo is funded through the European Gravitational Observatory (EGO),
    by the French Centre National de Recherche Scientifique (CNRS),
    the Italian Istituto Nazionale di Fisica Nucleare (INFN),
    and the Dutch Nikhef,
    with contributions by institutions from Belgium, Germany, Greece, Hungary,
    Ireland, Japan, Monaco, Poland, Portugal, Spain.
The authors are grateful for computational resources provided by the 
    LIGO Laboratory and supported by National Science Foundation Grants
    PHY-0757058, PHY-0823459 and PHY-1626190, and IUCAA LDG cluster Sarathi.

















§ MARGINALIZING OVER WAVEFORM SYSTEMATICS



Recently, several studies have investigated techniques to mitigate waveform systematics by marginalizing over waveform
model uncertainty, using an ensemble of current waveform models  <cit.>. 

 In these procedures, the waveform-marginalized posterior is the  weighted
average of the posteriors p_k(θ) derived from each
 waveform model k alone, weighted by the evidence Z_k  for (and prior p_k for) each model k: p(θ) = [∑_k
   p_k(θ)p_kZ_k]/∑_qp_q Z_q. Extending on this approach, Ashton and Diectrich <cit.> presented
   a hypermodel approach for analysing BNS signals with a suite of models to identify waveform systematics. 

In this appendix, we outline the approach outlined in  <cit.>: we use the RIFT parameter inference engine to construct the marginalized likelihoods and
posterior distributions.  Particularly for analysis of BNS, RIFT's low cost provides distinct advantages when attempting
these or other strategies to  mitigate waveform uncertainty by incorporating results from  costly models. 


Suppose we have two models A and B for GW strain, and use them to interpret a particular GW source.  We have prior
probabilities p(A|λ) and p(B|λ), characterizing our relative confidence in these two models for a source
with parameters λ.[For simplicity we will assume there are no internal model hyperparameters, but
  the method is easily generalized to include them.]   Suppose we have produced a RIFT analysis with each model for this
event, and have marginal likelihood functions L_A(λ) and L_B(λ) evaluated at a
single point λ.  
We can therefore construct the marginal likelihood for λ by averaging over both models:

    L_av(λ) = p(A|λ)  L_A(λ) + p(B|λ)  L_B(λ)

For simplicity the calculations in this work always adopt p(A|λ)=p(B|λ)=1/2.
We can therefore transparently integrate multi-model inference into RIFT as follows.  We assume we have a single grid
of points λ_k such that both (λ_k,  L_A(λ_k) and (λ_k,  L_B(λ_k)  can be
interpolated to produce reliable likelihoods and thus posterior distributions p_A(λ) and p_B(λ),
respectively.  At each point λ_k we therefore construct L_av(λ_k) by the above procedure. 
We then interpolate to approximate L̂(λ) versus the continuous parameters λ.


Operationally speaking, we construct model-averaged marginal likelihoods by the following procedure.  First, we
construct a fiducial grid for models A and B, for example by joining the grids used to independently analyze A and B.
We use an algorithm to integrate the extrinsic likelihood (ILE), a process where each candidate GW signal is compared to
a regular grid of candidate source parameters to produce an array of candidate likelihood values, to evaluate L_A(λ_k) and L_B(λ_k) on this grid <cit.>.  We construct L_av(λ_k) as above.  We use the combinations (λ_k,  L_av) with an algorithm to construct the
intrinsic posterior (CIP) from this sampled data, to construct a
model-averaged posterior distribution <cit.>. 


Our procedure resembles the approach suggested by Ashton and Khan, but we have organized the
calculation differently.  In that approach, AK used the evidences Z_A = ∫ L_A
p(λ) dλ and Z_B for the two
waveform models.  While we can compute both quantities with very high accuracy, we prefer to directly average between
waveform models at the same choice of intrinsic parameters (i.e., via Eq. (<ref>)) , to insure that marginalization over waveform models is completely decoupled from the interpolation
techniques used to construct L̂ from the sampled data.






