# Parameter estimation with non stationary noise in gravitational waves data

Sumit Kumar Max-Planck-Institut für Gravitationsphysik (Albert-Einstein-
Institut), D-30167 Hannover, Germany  Leibniz Universität Hannover, D-30167
Hannover, Germany  Alexander H. Nitz Max-Planck-Institut für
Gravitationsphysik (Albert-Einstein-Institut), D-30167 Hannover, Germany
Leibniz Universität Hannover, D-30167 Hannover, Germany  Xisco Jiménez Forteza
Max-Planck-Institut für Gravitationsphysik (Albert-Einstein-Institut), D-30167
Hannover, Germany  Leibniz Universität Hannover, D-30167 Hannover, Germany

###### Abstract

The sensitivity of gravitational-waves detectors is characterized by their
noise curves which determine the detector’s reach and the ability to
accurately measure the parameters of astrophysical sources. The detector noise
is typically modelled as stationary and Gaussian for many practical purposes.
However, physical changes in the state of detectors due to environmental and
instrumental factors, including extreme cases where a detector discontinues
observing for some time, introduce non-stationarity into the noise. Even slow
evolution of the detector sensitivity will affect long duration signals such
as binary neutron star (BNS) mergers. Mis-estimation of the noise behavior
directly impacts the posterior width of the signal parameters. This becomes an
issue for studies which depend on accurate localization volumes such as i)
probing cosmological parameters (such as Hubble constant, clustering bias)
using cross-correlation methods with galaxies, ii) doing electromagnetic
follow-up using localization information from parameter estimation done from
pre-merger data. We study the effects of dynamical noise on the parameter
estimation of the GW events. We develop a new method to correct dynamical
noise by estimating a locally-valid pseudo PSD which is normalized along the
time-frequency track of a potential signal. We do simulations by injecting the
BNS signal in various scenarios where the detector goes through a period of
non-stationarity with reference noise curve of third generation detectors
(Cosmic explorer, Einstein telescope). As an example, for a source where mis-
modelling of the noise biases the signal-to-noise estimate by even $10\%$, one
would expect the estimated localization volume to be either under or over
reported by $\sim 30\%$; errors like this, especially in low-latency, could
potentially cause follow-up campaigns to miss the true source location.

gravitational waves — power spectral density —binary neutron stars — third
generation detectors

##  1 Introduction

Since the first detection of gravitational waves (GW) from the merger of two
compact objects in 2015 (Abbott et al., 2016), GW science has come a long way
and GW detections have become a routine (Abbott et al., 2021a; Nitz et al.,
2021b). LIGO observatories in Livingston and Hanford (Aasi et al., 2015) have
already finished three science runs, namely O1, O2, and O3 (Abbott et al.,
2021a). The GW observatories: Virgo (Acernese et al., 2015) in Italy and KAGRA
in Japan (Akutsu et al., 2021) joined the LIGO-Hanford and LIGO-Livingston
observatories in O2 and O3 respectively. Currently, there are more than 90
confident GW detections from the merger of compact objects as identified by
various data analysis pipelines summed up as GWTC-3 catalog (Abbott et al.,
2021a) and 4-OGC catalog (Nitz et al., 2021b). In the upcoming decade, the
network of ground base GW detectors is set to be expanded with the addition of
LIGO-India (Saleem et al., 2022). The current ground base detectors are also
set to undergo phases of improvement in sensitivity (Abbott et al., 2018).
Additionally, the proposed third generation of ground based detectors
(Shoemaker et al., 2021) such as the Einstein telescope (ET) (Sathyaprakash et
al., 2012), and Cosmic Explorer (CE) (Reitze et al., 2019; Evans et al., 2021)
are expected to have an order of magnitude better sensitivity and will be able
to probe lower gravitational-wave frequencies ($f_{low}\sim 3-5$ Hz). This
will help in rapidly expanding the gravitational-wave catalog.

The ever growing catalog enables us to answer important questions about the
intrinsic population properties of the sources; these include the mass
function, spin distribution, and redshift distribution (Abbott et al., 2021b).
By inferring the localization volumes of individual GW merger events with the
help of parameter estimation (PE) methods, and by cross correlating the
localization volumes with galaxy catalogs, one can also estimate the
cosmological parameters such as the Hubble constant (Abbott et al., 2021c),
and large scale distribution of matter (Vijaykumar et al., 2020; Mukherjee et
al., 2021a; Libanore et al., 2021; Mukherjee et al., 2021b; Cañas-Herrera et
al., 2021). The reach of third-generation detectors also implies a large
population of loud BNS mergers which will be prime candidates for
electromagnetic-followup campaigns. It will also feasible to infer the sky
location from the early part of the signal and alert the electromagnetic (EM)
telescope before a merger and the subsequent electromagnetic emission occurs
(Chan et al., 2018; Nitz & Dal Canton, 2021). The correct estimation of the
localization volumes requires accurately understanding and modelling the noise
of the detector along with the signals.

Data from the GW detectors is analysed using the time series strain data
$s(t)$ which might contain a signal $h(t)$ and the detector noise $n(t)$. The
noise is typically assumed to be Gaussian, stationary, and colored. This makes
it easy to model the noise.

In GW searches, the data is analyze in chunks and the the power spectral
density (PSD) of the noise is typically estimated for a segment of data $\sim
256-1024$s using Welch’s method (Abbott et al., 2020; Usman et al., 2016). It
is assumed that the PSD is constant during this time. This PSD is then used to
calculate the signal-to-noise ratio (SNR) $\rho$ by applying matched filtering
to the data for a bank of template waveforms (Allen et al., 2012). Whenever
$\rho$ is obtained above a certain threshold, it is called a single detector
trigger. The local trigger density depends on the the PSD and if a local
deviation from the average PSD is not taken into account, it can affect the
estimation of the local trigger density (Mozzon et al., 2020).

For parameter estimation, the estimated PSD is used to calculate the
likelihood function and the PSD is assumed to be constant during the time of
the signal in the detector frequency band. Again, any significant local time
variation from the estimated PSD will affect the likelihood function and
therefore affect the width of the posterior and hence the reported
uncertainties in the marginalized distribution of individual parameters. For
the current generation detectors, we do not expect this bias to be
statistically significant for loud signals because most of the signals are
short lived ($\sim\mathcal{O}(1-10)$ seconds for BBH and
$\sim\mathcal{O}(100)$ seconds for BNS) in the detector sensitivity band. But
for third generation detectors, because of the better sensitivity at low
frequencies, we expect the signals to be lasting for
$~{}\sim\mathcal{O}(100-60000)$ seconds. Due to very long signal durations,
any change in noise characteristics due to evolution in the state of the
detectors will cause the noise to be non-stationary; modelling the non-
stationarity is needed to ensure the most accurate source localization and
parameter estimates. An important aspect of long duration signals is that a
fraction of them can be identified from few minutes to few hours before the
merger can be identified and located in sky with localization information
(Nitz & Dal Canton, 2021). This can be used to alert various telescopes to
point in right direction to be ready for observing the electromagnetic follow-
up observations. A mis-estimation of localization volume for pre-merger alert
will end up in either wasting telescope time and resources (in case of over
estimated localization volume) or in worst case scenario missing out on
observing the event (under estimated localization volume).

In this paper we explore the modeling of non-stationary noise and it’s
implications on parameter estimation. We introduce an efficient method to
account for non-stationarity which is slowly varying compared to the signal
evolution. This is accomplished by calculating a corrected PSD along the track
of the signal by dividing the time-frequency track of a signal into shorter
duration segments and estimating a local PSD for each time bin separately. The
corrected PSD can be constructed by selecting the frequency range from each
local PSD estimate corresponding to the expected signal contribution in the
same time range. This method of estimating the PSD depends on getting a
preliminary constraint from searches for trigger parameters to identify the
track of the signal.

##  2 Gravitational waves data analysis

The signal from the merger of two compact objects: binary black holes (BBHs),
binary neutron stars (BNS), or neutron star-Black hole (NSBH), as seen in
detector, consists of a time series $h(t)$ buried in the noisy data of the
detector $n(t)$.

$s(t)=\begin{cases}n(t)+h(t),&\text{When signal is present}\\\
n(t),&\text{otherwise}\end{cases}$ (1)

The shape of the signal $h(t)$ depends on the intrinsic parameters of the
source such as the component masses ($m_{1,2}$), the component spins
($s_{1x,2x},s_{1y,2y},s_{1z,2z}$) as well as on the extrinsic parameters such
as the luminosity distance to the source ($D_{L}$), the inclination angle of
the binary to line of sight ($\iota$), the polarization angle ($\psi$), the
coalescence phase ($\phi$), the sky localization (right ascension: $RA$,
declination: $dec$). The noise $n(t)$ has contributions from many sources such
as quantum noise, thermal noise, Newtonian noise, and transient noise of
unknown origins (Abbott et al., 2020). The overall noise is typically
described as a stochastic process modelled as a Gaussian distribution:

$p(\mathbf{n})=\frac{1}{\sqrt{2\pi|C_{ij}|}}\exp(\frac{1}{2}\sum_{i,j}(n_{i}-\mu)C_{ij}^{-1}(n_{j}-\mu))$
(2)

Where $\mathbf{n}$ is a time series vector of a certain noise realization with
a given sample rate. $\mu$ is the expectation value of the noise and $C_{ij}$
is the covariance between $i^{th}$ and $j^{th}$ bin in the time series. Both
these quantities can be estimated from the data. The noise is said to be
stationary if the covariance matrix $C_{ij}$ only depend on the lag between
bins $i$ and $j$ i.e. $|t_{i}-t_{j}|$. A useful property of stationary noise
is that the covariance matrix takes a positive diagonal form in frequency
domain i.e. $C_{ij}=\delta_{ij}S_{n}(f_{i})$. Here $S_{n}(f_{i})$ is defined
as power spectral density which is the Fourier transform of the correlation
function $C(\tau)$ where $\tau=|t_{i}-t_{j}|$ is the time lag. The validity of
the assumptions of the noise being Gaussian and stationary can be tested by
various techniques, for example, with large number of injections in simulated
Gaussian noise and detector data, the percentile-percentile test is done with
posterior samples to see if the recovered credible intervals matches with the
statistical expectations (Kulkarni & Capano, 2021). One can also use the
distribution of PSD variation statistic $v_{s}$ to identify the stretches of
the data which do not agree with the expected statistical distribution of
$v_{s}$ for the Gaussian and stationary noise (Mozzon et al., 2020), or by
measuring the noise covariance of the data (Edy et al., 2021).

Figure 1: Different types of local variation in the PSD are considered here.
To track the local variation in PSD, we use the PSD variation statistic
$v_{s}$ described in (Mozzon et al., 2020). In all the figures, the lower
panel shows the strain data for 1000 seconds. Top panel shows the local PSD
variation statistic with respect to the PSD estimated for the entire duration
(1000 seconds). In the inset, the distribution of all the PSD variation values
is shown. Different scenarios of non-stationarities considered in these
figures are: (a) Baseline noise model with Gaussian and stationary
assumptions. The local variations in the PSD are due to the Gaussian random
fluctuations. (b) The baseline noise model is acted upon by a filter with a
bump in the middle. We use $B(t)=B_{0}sin(\frac{\pi}{\Delta T}t)$ with
$[B_{0},\Delta T]$ = $[0.2,200sec]$ so that it produces a bump in the
amplitude of positive half sine wave between t=400 to t=600 seconds. (c) The
floor of the baseline PSD is lifted in the middle (at t=500 seconds) with
$B(t)=0.2~{}\text{for}~{}t>500$sec (d) The floor of the PSD is lifted at t=400
second and drops down to the base scale slowly in due time by t=600 seconds.

###  2.1 Effects of estimated PSD on the infered signal properties

The signal to noise ratio (SNR) of a signal template $h$ is defined as
$\rho\equiv\frac{(s|h)}{\sqrt{(h|h)}}$ where the inner product $(a|b)$ between
two time series $a(t)$ and $b(t)$ (and their counterpart in Fourier domain
$a(f)$ and $b(f)$) is defined as:

$(a|b)\equiv 2\int_{-\infty}^{\infty}\frac{a(f)b^{\star}(f)}{S_{n}(|f|)},$ (3)

where, $S_{n}(f)$ is the estimated power spectrum of the noise over a chosen
time window. The GW searches construct the SNR time series to identify a
signal in the data. Whenever a template is found to have SNR greater than a
pre-defined threshold, it is stored as a single detector trigger. These single
detector triggers are then passed through a ‘coincidence test’ where only
those triggers, which are within the light travel time (plus few milliseconds
to account for errors in estimating trigger time) between the two detectors,
survive. As the estimation of SNR depends on estimated PSD $S_{n}(f)$, mis-
estimation of PSD, will have an effect on the SNR and hence on the trigger
density. There are methods available to correct for these effects at the level
of the searches, such as applying a psd drift correction (Zackay et al., 2021)
and applying dynamic normalization using the PSD variation statistic (Mozzon
et al., 2020). However, note that both methods assume a constant change in the
overall amplitude of the noise and that this is constant over the duration of
the signal. We may expect this to be a reasonable assumption for short-
duration signals.

For the parameter estimation, we describe the signal $h(\theta)$ in terms of
parameters $\theta$ of the waveform model and perform Bayesian analysis to
obtain the posterior distribution $p(\theta|d)$ which represents the
probability distribution over the parameters $\theta$ given the observed data
$d$ and model,

$p(\theta|d)=\frac{\mathcal{L}(d|\theta)\pi(\theta)}{p(d)},$ (4)

where $\mathcal{L}(d|\theta)$ is the likelihood function: the probability of
obtaining the data given the parameter $\theta$ describing model.
$\pi(\theta)$ is the prior probability distribution for the parameters, and
$p(d)$ is known as Bayesian evidence or marginalized likelihood which acts as
normalization factor. For the Gaussian noise model, the likelihood function in
the frequency domain takes the form:

$\mathcal{L}\propto\exp(-\frac{1}{2}n^{\star}(f)S_{n}^{-1}(f)n(f))\,.$ (5)

As shown in (Edy et al., 2021), the effects of mis-estimating the PSD can be
calculated by doing a linear series expansion of the waveform template $h$
around the true signal value $h_{0}$. The mis-estimation of the $S_{n}(f)$
only results in over (or under) estimating the posterior width. In other
words, if we use two different estimates of $S_{n}(f)$, the posteriors will
still be centred at the same values but the width of the posteriors will be
different. Specifically, at linear order in the signal $h_{0}$, the covariance
matrix of the waveform parameters – the inverse of the Fisher matrix
(Vallisneri, 2008) – scales linearly with $S_{n}(f)$, while the SNR of the
signal as $S_{n}^{-1}(f)$. Thus, as a rule of thumb, a PSD variation of type
$S_{n}(f)\rightarrow\alpha S_{n}(f)$ – a variation of the noise floor – will
broaden/narrow the width on the posterior distributions of the physical
parameters by the same factor $\alpha$, which results in increasing their
inferred uncertainty. This may become very relevant for the studies where the
localization volumes are very important, such as i) pre-merger localization,
ii) estimating Hubble parameter using cross correlation methods with galaxy
catalogs, and iii) probing large scale structures using two point correlation
function of localization volumes with next generation GW detectors. For all
the studies mentioned above, the estimation of correct localization volume is
important for unbiased estimators. In a study by Mozzon et al. (2021), through
simulation of BNS signals, it is shown that though non-stationary noise will
not be a limiting factor for resolving Hubble tension for current generation
of detectors, it will play important roll for 3G detectors. On the other hand,
the non-stationary models considered are not expected to introduce any bias on
the estimated parameters at least at linear order (Edy et al., 2021).

Figure 2: We show a frequency-time track of a GW signal from the merger of
binary neutron star. The track of the signal can be used to identify the times
corresponding to the frequency bins to estimate the PSD of the noise.

If the data is non-stationary, the covariance matrix in the frequency domain
also contains the off diagonal elements and the estimation of the likelihood
function 5 becomes computationally expensive. So it is difficult to just model
the non-stationarities in a straight forward way without compromising the
computational costs. In one of the studies, (Edy et al., 2021) shows that it
is important to account for the non-stationarity in the signals that lasts
longer than few minutes such as BNS signals in current detectors and BBH/BNS
signals in third generation GW detectors.

Table 1: The description of various scenarios of non stationary noise considered. We use noise sensitivity curve for Cosmic Explorer design sensitivity (Hall et al., 2021) with low frequency cutoff of $f_{low}=5Hz$. The non-stationarity in the noise is introduced by applying a filter the baseline noise model with a time varying amplitude for a fixed time window. To estimate the SNR, we estimate the PSD by two methods: i) Using a long segment of 8000 seconds surrounding the trigger time $t_{c}$, using Welch’s method by dividing the whole segment into smaller chunks of 8 seconds. And ii) Estimated along the frequency-time track of the reference BNS signal. The PSD estimated along the track of the signal is used as reference to estimate fractional difference in SNR $=\absolutevalue{\delta\rho}/\rho_{\text{ref}}$.  S. |  |  |
---|---|---|---
No. | Amplitude variation |  Description | SNR difference (%)
1. | $B(t)=\begin{cases}B_{0},&:t>t_{shift}\\\ 0,&:\text{otherwise}\end{cases}$ |  A positive shift, described by a step function, in the amplitude of the Gaussian noise at $t_{shift}=t_{c}-190$ seconds. $B_{0}=0.8$, and $t_{c}$ is the trigger time of the signal. Therefore the signal spend last 190 seconds (before merger) buried in the noise characterized by, where noise floor is raised according to the equation. | 43.6
2. | $B(t)=\begin{cases}-\frac{B_{0}}{\Delta T}(t-t_{1})+B_{0},&:t_{1}<t<t_{1}+\Delta T\\\ 0,&:\text{otherwise}\end{cases}$ |  A positive shift in the amplitude of the noise followed by a linear decay to original amplitude (base scale). We use $B_{0}=0.8$, $\Delta T=600$ seconds, and $t_{1}=t_{c}-650$. The signal spent 600 seconds of time in non stationary noise where the noise floor was above the base scale. | 11.47
3. | $B(t)=\begin{cases}B_{0}\sin{(\frac{\pi(t-t_{1})}{\Delta T})},&:t_{1}<t<t_{1}+\Delta T\\\ 0,&:\text{otherwise}\end{cases}$ |  A bump described by a half positive cycle of sine function with $B_{0}=0.8$, $\Delta T=600$ seconds, and $t_{1}=t_{c}-650$. The signal spend 600 seconds of time in non-stationary noise where noise floor was lifted above the base scale. | 21.19
4. | $B(t)=\begin{cases}B_{0},&:t<t_{shift}\\\ 0,&:\text{otherwise}\end{cases}$ |  A drop in the amplitude of the Gaussian noise at $t_{shift}=t_{c}-10$ seconds. $B_{0}=0.8$, and $t_{c}$ is the trigger time of the signal. So the signal spend all the time until 10 seconds prior to the merger in the noise characterized by where noise floor is raised with respect to base scale. | 9.11

##  3 Estimating local variation in the PSD

To check for the stationarity of the data, one can estimate the noise
covariance. One can also use other methods such as PSD variation statistics
described in (Mozzon et al., 2020). The non-stationarity of the data will lead
to difference in actual power spectrum $S_{A}(f)$ and estimated one $S_{E}(f)$
modelled with the help of PSD variation statistics $v_{s}$,

$S_{A}(f)=v_{s}S_{E}(f)$ (6)

where $S_{E}(f)$ is the estimated power spectral density over long stretch of
data using the Welch method with smaller overlapping segments. Using a filter
$\mathcal{F}(f)=\frac{h(f)}{S_{E}(f)}$, with approximating $|h(f)|\propto
f^{-7/6}$ (dominant amplitude behaviour of the CBC signal during the inspiral
phase), $v_{s}$ can be estimated by computing the convolution between the data
$s(f)$ and filter $\mathcal{F}(f)$.

The non-stationarity in the detector noise can be attributed to various
factors of known or unknown origin. There are different non-stationary noise
models considered in the literature (Vajente et al., 2020; Hebbal et al.,
2019; Edy et al., 2021). A few examples are:

  * •

local fluctuations in the detector noise due to non stationary component: It
might be the case that there is some non-stationary component in the overall
noise. This can be modelled as follows: $n(t)=An_{1}(t)+B(t)n_{2}(t)$ (Edy et
al., 2021). Here the first component $n_{1}(t)$ is stationary Gaussian noise
and second component $n_{2}(t)$ represents noise of unknown origin with a
time-dependent component $B(t)$. A special case of above is when the noise is
Gaussian but the overall amplitude varies as a function of time such that
$n(t)=(1+B(t))A_{0}n_{1}(t)$.

  * •

A shift in the overall PSD: There are times when detector goes out of lock and
comes back online but with overall noise floor raised above the reference PSD
and slowly comes back to the more stable configuration.

Figure 3: The estimated amplitude spectral density (ASD) as a function of
frequency. The blue curve represents the the ASD estimated for long segment
($~{}8000$ seconds) around the signal and red curve represent the ASD
estimated along the track of the signal as described in the text. We inject
the BNS signal in the non stationary noise as described in the Table 1. The
shaded region represents the frequency range along the track of the signal
when detector went through the period of non stationary noise. During this
time, the noise amplitude is increased as a function of time. This is
reflected in the raised (estimated) noise floor. The figures (a),(b),(c), and
(d) corresponds to the scenarios (1),(2),(3), and (4) in the table 1
respectively

In figure 3, we show some examples of the types of non-stationarities
described above. We use a baseline noise model (Gaussian, stationary colored
noise) and apply a filter upon it to scale the amplitude of Gaussian noise as
function of time according to desired non stationary nature. This may not be
the most general way to model the non-stationarity but we consider it to test
our correction method in parameter estimation. We use PSD variation statistic
$v_{s}$ to track the non-stationary nature of the noise. For the examples
shown in the figure, we generated $\sim 1000$ seconds of data with Gaussian
noise with one of the available noise model ’aLIGOZeroDetHighPower’ in the
publicly available software package LALSuite (LIGO Scientific Collaboration,
2018). We call this baseline noise model. For the baseline model, the local
PSD variation is due to purely Gaussian random fluctuations so the PSD
variation statistic $v_{s}$ have mean around $<v_{s}>=1$ with standard
deviation consistent with that of the Gaussian fluctuations. We generate the
non stationary noise by applying a filter on baseline noise such that the
amplitude of the Gaussian noise varies as desired function. This type of PSD
variation can be modelled as,

$n(t)=(1+B(t))A_{0}n_{G}(t)$ (7)

where $A_{0}n_{G}(t)$ is the stationary, Gaussian, colored noise, $B(t)$
models the variation in the amplitude. $B(t)=0$ corresponds to the baseline
noise model. As shown in figure 3, we introduce some arbitrary non-
stationarity in the noise at various times and use PSD variation statistics
$v_{s}$ to track it.

##  4 Modeling the non-stationarity for parameter estimation

In most of the standard parameter estimation (PE) pipelines, to obtain the
posterior distribution 4, a sampling method (such as MCMC, nested sampling) is
used with a likelihood model and prior probability distribution over the model
parameters. The PSD is estimated for a segment of the data around the signal
using the Welch’s method with smaller overlapping segments or spline-modelling
such as implemented by the Bayesline algorithm (Littenberg & Cornish, 2015).
For example, to estimate the PSD for PE in 4OGC and 3OGC analysis (Nitz et
al., 2021b, a), 512 seconds of data centered around trigger is used to
estimate the PSD with overlapping segments of 8 seconds (unless the data is
corrupted in this window, in which case, this window is shifted and/or
shrunk). The estimated PSD is then used in the calculation of likelihood
function 5. The averaged PSD over the long segment surrounding the signal wash
away any local fluctuations (or non-stationarities) that might be present
along the track of the signal. In order to incorporate the non-stationarity,
one has re-define the likelihood function by incorporating non diagonal terms
in covariance matrix estimated in frequency domain.

In this work, we propose a new method to model certain types of non-
stationarities where the strain data can be divided into smaller chunks of
”stationary” or ”semi-stationary” parts: such as slowly varying Gaussian noise
amplitude, sudden lift or drop in the noise floor, even missing some duration
of data (for example detector being offline) for very long duration signals.
For 3G GW detectors such as Einsteins Telescope (ET) or cosmic explorer (CE),
due to their low frequency reach, the BBH signals can be
$\sim\mathcal{O}(100)-\mathcal{O}(1000)$ seconds long and BNS signals can be
$\sim\mathcal{O}(1000)-\mathcal{O}(60000)$ seconds long depending on the value
of low frequency cutoff. During this time, it is possible that detector goes
through periods of non-stationarities. Even for the current generation
detectors, a BNS signal in the detector band is usually few hundred seconds
long. When the track of the signal in detector band is long enough, we can
estimate the PSD along the track of the signal. To do so, we can use reference
signal parameters obtained from the searches as trigger with maximum SNR. The
track of the reference signal can be drawn in the frequency-time graph and it
can be broken in frequency bins (and corresponding time bin). We can then
estimate the PSD for each frequency bins by using strain data from
corresponding time bin and stack the estimated PSD in each frequency bin to
estimate a corrected PSD along the track of the signal. In figure 2, we show
frequency time track of a reference BNS signal. We divide the track of the
signal in frequency-time graph into smaller bins. We assume that the
variations in the detector noise across each time (frequency) bins are
negligible so that we can model the detector noise in each bin to be Gaussian
and stationary.

In Table 1, we show various scenarios of arbitrary non stationary noise we
considered in simulations. We inject a BNS signal with parameters
($m_{1,source}$ = $1.5M_{\odot}$, $m_{2,source}$ = $1.3M_{\odot}$, $D_{L}$ =
$100$ Mpc) in the data such that the detector goes through a period of non-
stationarity during the course of the signal. We estimate the SNR using two
estimated PSDs: i) assuming the detector noise is constant during the signal
and estimating average PSD from long stretch of data around signal, and ii)
estimating PSD along the track of the signal as described above. We estimate
the fractional SNR difference $(\frac{\delta_{\rho}}{\rho})$ arising from two
PSDs and found that there can be significant fractional SNR difference (upto
$\sim 40\%$) for the arbitrary non stationary noise we consider here. In
reality, these differences might not be as extreme but even a few percent
difference in fractional SNR will have an effect on localization volumes. We
use the baseline noise model from the design sensitivity of 3G detector cosmic
explorer. We use 3G detector noise curve only for the reason that BNS signals
are expected to last for very long time in the detector frequency band. The
current method can be applied to current detectors as well when we see the
non-stationarity (similar to ones considered in the examples) along the track
of a long BNS signal lasting few hundred seconds.

In figure 3, we show the estimated PSD in both ways for all the injection
scenarios described in table 1. The frequency and time bins are chosen such
that the local variation in the PSD can be tracked.

Table 2: The specifications of 3G detector network (location, low frequency cutoff $f_{\mathrm{low}}$) considered in this study. We use the design sensitivity noise curves for ET and CE. These detectors configuration for CE and ET are also used in previous work (Hall & Evans, 2019; Nitz & Dal Canton, 2021; Kumar et al., 2021) Abbreviation | Observatory | $f_{\textrm{low}}$ | Latitude | Longitude
---|---|---|---|---
$CE^{U}$ | Cosmic Explorer USA | 5.2 | 40.8 | -113.8
$CE^{A}$ | Cosmic Explorer Australia | 5.2 | -31.5 | 118.0
$E$ | Einstein Telescope | 5 | 43.6 | 10.5

###  4.1 Simulations and Parameter Estimation

To put this method to test, we inject a GW signal from the merger of a BNS
($m_{1,source}$ = $1.5M_{\odot}$, $m_{2,source}$ = $1.3M_{\odot}$, $D_{L}$ =
$100$ Mpc) in a network of third generation gravitational wave detector. We
consider two cosmic explorer (CE) detectors (one located in USA and other in
Australia), and Einstein telescope (ET) (Europe). We summarize the properties
of the detector network in Table 2. The CE is proposed to be built in two
phases, known as CE1 and CE2 (Hall et al., 2021), with second phase having
better sensitivity. We consider the design sensitivity for both the detectors.
$E_{1}$, $E_{2}$, and $E_{3}$ in our notation denote the three effetively
independent detectors that make up the ET observatory. We choose 3G detector
network to demonstrate the differences solely because it is very likely for
signals to be present when detector undergoes the non stationary times. The
methods described here will be applicable to any scenario with long duration
signals where the data can be divided into smaller chunks of stationary (or
semi stationary) times to estimate the PSD along the track of the signal.

In the 3G detector network we considered, we introduce arbitrary non-
stationarities modelled in 7 as follows:

  * •

$CE^{U}$: We introduce non-stationarity described by a step function
$B(t)=B_{0}$ for $t>t_{c}-180$ seconds,i.e. noise floor is raised around 180
seconds before the trigger time $t_{c}$.

  * •

$CE^{A}$: Here we introduce the non-stationarity in the detector noise by
dropping the noise floor of the detector at base level 10 seconds prior to the
trigger time. Mathematically, $B(t)=B_{0}$ for $t<t_{c}-10$.

  * •

$E_{1}$: In this detector we again used positive step function $B(t)=B_{0}$
for $t>t_{c}-50$ seconds.

  * •

$E_{2}$: Noise floor drops to base scale 50 seconds prior to trigger time:
$B(t)=B_{0}$ for $t<t_{c}-50$.

We do not introduce any non-stationarity in the detector $E_{3}$. We use
$B_{0}=0.8$ for this simulation. Notice that all these non-stationarities will
have different effect in the fractional SNR difference of the detector
network.

We do Bayesian parameter estimation, using publicly available code PyCBC
Inference (Biwer et al., 2019). We use uniform priors in detector frame chirp
mass $\mathcal{M}$, mass ratio $(q=m_{1}/m_{2});m_{1}>=m_{2}$, comoving
volume, and time of arrival $t_{c}$. We use the isotropic prior for RA, dec,
inclination angle, and polarization. We use TaylorF2 waveform model (Blanchet
et al., 1995; Faye et al., 2012) implemented in LALSuite (LIGO Scientific
Collaboration, 2018) to simulate the signal to be injected in the noise, and
for recovery while doing parameter estimation. To estimate the likelihood
function, we make use of the heterodyne likelihood model described in
(Cornish, 2010; Finstad & Brown, 2020; Zackay et al., 2018). To reduce the
computational costs, we use the sample rate of 1024 to generate the data. We
use low frequency cutoff $f_{low}=5.1$ Hz for $CE^{U,A}$ and $f_{low}=5$Hz for
$ET$, and high frequency cutoff $f_{high}=512$ Hz for all detectors to
estimate likelihood function.

In figure 4, we show the difference in the localization volume inferred from
two methods of estimating the PSD. It shows the difference localization
volumes for two cases. The difference in network SNR corresponding to maximum
likelihood obtained for this example is $4.76\%$.

Figure 4: This figure shows difference in localization volume for right
ascension and declination angle arising for two cases. Blue contours
represents the $50\%$ and $90\%$ credible intervals when an constant PSD is
used for parameter estimation. The orange contours represents same when a PSD
estimated along the track of the signal is used The red lines represents the
injected values. For this simulation, we used a loud BNS injection
($m_{1,source}:1.5~{}M_{\odot}$, $m_{2,source}:1.3~{}M_{\odot}$,
$d_{L}:100$Mpc) in a third generation detector network of two cosmic explorer
and Einstein telescope. See text for the details of detector network and types
of non-stationarity considered.

##  5 Summary and Discussion

The current methods of parameter estimation (PE) using GW data assume
stationary noise. Estimation of the PSD over a large chunk of data washes away
any local fluctuations due to any possible non-stationarity in the data. The
stationary noise has a useful property that the correlation function of time
series between two times $t_{1}$ and $t_{2}$ depends only on the time
difference $|t_{1}-t_{2}|$. In the frequency domain, it means that the
covariance matrix is diagonal. These assumptions help in simplifying the
likelihood function and make it faster to calculate. If any non-stationarity
is present in the data and not accounted for, it will not primarily shift
estimates of source parameters, but will most directly affect width of the
posterior distributions and hence reported errors on parameters. A simple
method to incorporate (any) non stationarity present in the data is to
estimate the off diagonal terms in the likelihood function 5. While exact,
this makes the estimation of the likelihood function impractically expensive
for most purposes.

In this work, propose a fast method to incorporate the non stationarity
specially for the long signals such as BNS. We propose to estimate the PSD
along the track of the signal in the frequency-time graph. We assume that the
data can be broken into smaller chunks where it can be assumed to be
stationary or semi stationary (such as slowly varying amplitude of the
Gaussian noise). We do simulations by injecting a long BNS signal in various
scenarios of arbitrary non stationary noise with 3G detectors and show that
there there will be difference in the SNR estimated if these effects ts are
not accounted for. We also show that it will lead to the difference in
estimated localization volume obtained from parameter estimation.

This method is reliable for the types of non-stationarities where the strain
data can be broken into piece-wise smaller segments of semi-stationary data
such as when amplitude of the Gaussian noise is changing slowly in the smaller
time bins. It also relies on identifying nearly the right track of the signal
in frequency-time graph. In practice, identifying this track should not be
problematic; several bootstrapping schemes are practical, including using the
initial point estimates for the coalescence time and chirp mass from standard
compact-binary searches.

We note that mis-characterization of the noise also has an effect not-only on
the estimation of signal parameters but also on the detection of distance
sources due to the mis-estimation of a signal’s SNR. The method for
calculating a corrected PSD can also be efficiently applied to searches by a
simple two-step procedure. In the first stage, an average PSD is used to
calculated an initial SNR estimate and identify potential candidates. Once the
initial candidates are identified, the SNR can be re-weighted using the
corrected PSD. This is natural extension of the dynamic normalization outlined
in Mozzon et al. (2020).

We present this as method to incorporate certain types of non-stationarities
in the data for parameter estimation. For future work, we plan to modify this
method to incorporate more general types of non-stationarities by estimating
the PSD along the track of the signal using PSD variation statistic. We also
plan to apply this method to study non-stationarities present in the data from
the current detectors by identifying those stretches and doing PE by
performing BNS injections.

We acknowledge the Max Planck Gesellschaft. We thank the computing team from
AEI Hannover for their significant technical support. SK thanks members of
Observational relativity and cosmology group at AEI for feedback and useful
comments. SK also thanks Aditya Vijaykumar for feedback.

## References

  * Aasi et al. (2015) Aasi, J., et al. 2015, Class. Quant. Grav., 32, 074001, doi: 10.1088/0264-9381/32/7/074001
  * Abbott et al. (2016) Abbott, B., Abbott, R., Abbott, T., et al. 2016, Physical Review D, 93, doi: 10.1103/physrevd.93.122003
  * Abbott et al. (2018) Abbott, B. P., et al. 2018, Living Rev. Rel., 21, 3, doi: 10.1007/s41114-020-00026-9
  * Abbott et al. (2020) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2020, Classical and Quantum Gravity, 37, 055002, doi: 10.1088/1361-6382/ab685e
  * Abbott et al. (2021a) Abbott, R., et al. 2021a.  https://arxiv.org/abs/2111.03606
  * Abbott et al. (2021b) —. 2021b, Astrophys. J. Lett., 913, L7, doi: 10.3847/2041-8213/abe949
  * Abbott et al. (2021c) —. 2021c.  https://arxiv.org/abs/2111.03604
  * Acernese et al. (2015) Acernese, F., et al. 2015, Class. Quant. Grav., 32, 024001, doi: 10.1088/0264-9381/32/2/024001
  * Akutsu et al. (2021) Akutsu, T., et al. 2021, PTEP, 2021, 05A101, doi: 10.1093/ptep/ptaa125
  * Allen et al. (2012) Allen, B., Anderson, W. G., Brady, P. R., Brown, D. A., & Creighton, J. D. E. 2012, Phys. Rev. D, 85, 122006, doi: 10.1103/PhysRevD.85.122006
  * Biwer et al. (2019) Biwer, C. M., Capano, C. D., De, S., et al. 2019, Publ. Astron. Soc. Pac., 131, 024503, doi: 10.1088/1538-3873/aaef0b
  * Blanchet et al. (1995) Blanchet, L., Damour, T., Iyer, B. R., Will, C. M., & Wiseman, A. G. 1995, Phys. Rev. Lett., 74, 3515, doi: 10.1103/PhysRevLett.74.3515
  * Cañas-Herrera et al. (2021) Cañas-Herrera, G., Contigiani, O., & Vardanyan, V. 2021, The Astrophysical Journal, 918, 20, doi: 10.3847/1538-4357/ac09e3
  * Chan et al. (2018) Chan, M. L., Messenger, C., Heng, I. S., & Hendry, M. 2018, Phys. Rev. D, 97, 123014, doi: 10.1103/PhysRevD.97.123014
  * Cornish (2010) Cornish, N. J. 2010.  https://arxiv.org/abs/1007.4820
  * Edy et al. (2021) Edy, O., Lundgren, A., & Nuttall, L. K. 2021, Physical Review D, 103, doi: 10.1103/physrevd.103.124061
  * Evans et al. (2021) Evans, M., et al. 2021.  https://arxiv.org/abs/2109.09882
  * Faye et al. (2012) Faye, G., Marsat, S., Blanchet, L., & Iyer, B. R. 2012, Class. Quant. Grav., 29, 175004, doi: 10.1088/0264-9381/29/17/175004
  * Finstad & Brown (2020) Finstad, D., & Brown, D. A. 2020, Astrophys. J. Lett., 905, L9, doi: 10.3847/2041-8213/abca9e
  * Hall & Evans (2019) Hall, E. D., & Evans, M. 2019, Classical and Quantum Gravity, 36, 225002, doi: 10.1088/1361-6382/ab41d6
  * Hall et al. (2021) Hall, E. D., et al. 2021, Phys. Rev. D, 103, 122004, doi: 10.1103/PhysRevD.103.122004
  * Hebbal et al. (2019) Hebbal, A., Brevault, L., Balesdent, M., Talbi, E.-G., & Melab, N. 2019, Bayesian Optimization using Deep Gaussian Processes.  https://arxiv.org/abs/1905.03350
  * Kulkarni & Capano (2021) Kulkarni, S., & Capano, C. D. 2021, On the reliability of parameter estimates in the first observing run of Advanced LIGO.  https://arxiv.org/abs/2011.13764
  * Kumar et al. (2021) Kumar, S., Vijaykumar, A., & Nitz, A. H. 2021.  https://arxiv.org/abs/2110.06152
  * Libanore et al. (2021) Libanore, S., Artale, M. C., Karagiannis, D., et al. 2021, JCAP, 02, 035, doi: 10.1088/1475-7516/2021/02/035
  * LIGO Scientific Collaboration (2018) LIGO Scientific Collaboration. 2018, LIGO Algorithm Library - LALSuite, free software (GPL), doi: 10.7935/GT1W-FZ16
  * Littenberg & Cornish (2015) Littenberg, T. B., & Cornish, N. J. 2015, Phys. Rev. D, 91, 084034, doi: 10.1103/PhysRevD.91.084034
  * Mozzon et al. (2021) Mozzon, S., Ashton, G., Nuttall, L. K., & Williamson, A. R. 2021, Does non-stationary noise in LIGO and Virgo affect the estimation of $H_{0}$?  https://arxiv.org/abs/2110.11731
  * Mozzon et al. (2020) Mozzon, S., Nuttall, L. K., Lundgren, A., et al. 2020, Classical and Quantum Gravity, 37, 215014, doi: 10.1088/1361-6382/abac6c
  * Mukherjee et al. (2021a) Mukherjee, S., Wandelt, B. D., Nissanke, S. M., & Silvestri, A. 2021a, Phys. Rev. D, 103, 043520, doi: 10.1103/PhysRevD.103.043520
  * Mukherjee et al. (2021b) Mukherjee, S., Wandelt, B. D., & Silk, J. 2021b, Mon. Not. Roy. Astron. Soc., 502, 1136, doi: 10.1093/mnras/stab001
  * Nitz et al. (2021a) Nitz, A. H., Capano, C. D., Kumar, S., et al. 2021a, The Astrophysical Journal, 922, 76, doi: 10.3847/1538-4357/ac1c03
  * Nitz & Dal Canton (2021) Nitz, A. H., & Dal Canton, T. 2021, The Astrophysical Journal Letters, 917, L27, doi: 10.3847/2041-8213/ac1a75
  * Nitz et al. (2021b) Nitz, A. H., Kumar, S., Wang, Y.-F., et al. 2021b, 4-OGC: Catalog of gravitational waves from compact-binary mergers.  https://arxiv.org/abs/2112.06878
  * Reitze et al. (2019) Reitze, D., et al. 2019, Bull. Am. Astron. Soc., 51, 035.  https://arxiv.org/abs/1907.04833
  * Saleem et al. (2022) Saleem, M., et al. 2022, Class. Quant. Grav., 39, 025004, doi: 10.1088/1361-6382/ac3b99
  * Sathyaprakash et al. (2012) Sathyaprakash, B., et al. 2012, Class. Quant. Grav., 29, 124013, doi: 10.1088/0264-9381/29/12/124013
  * Shoemaker et al. (2021) Shoemaker, D. H., et al. 2021.  https://arxiv.org/abs/2112.12718
  * Usman et al. (2016) Usman, S. A., et al. 2016, Class. Quant. Grav., 33, 215004, doi: 10.1088/0264-9381/33/21/215004
  * Vajente et al. (2020) Vajente, G., Huang, Y., Isi, M., et al. 2020, Phys. Rev. D, 101, 042003, doi: 10.1103/PhysRevD.101.042003
  * Vallisneri (2008) Vallisneri, M. 2008, Phys. Rev. D, 77, 042001, doi: 10.1103/PhysRevD.77.042001
  * Vijaykumar et al. (2020) Vijaykumar, A., Saketh, M. V. S., Kumar, S., Ajith, P., & Choudhury, T. R. 2020\.  https://arxiv.org/abs/2005.01111
  * Zackay et al. (2018) Zackay, B., Dai, L., & Venumadhav, T. 2018.  https://arxiv.org/abs/1806.08792
  * Zackay et al. (2021) Zackay, B., Venumadhav, T., Roulet, J., Dai, L., & Zaldarriaga, M. 2021, Physical Review D, 104, doi: 10.1103/physrevd.104.063034