pygwb: A Python-based Library for Gravitational-wave Background Searches

The collection of gravitational waves (GWs) that are either too weak or too numerous to be individually resolved is commonly referred to as the gravitational-wave background (GWB). A confident detection and model-driven characterization of such a signal will provide invaluable information about the evolution of the universe and the population of GW sources within it. We present a new, user-friendly, Python-based package for GW data analysis to search for an isotropic GWB in ground-based interferometer data. We employ cross-correlation spectra of GW detector pairs to construct an optimal estimator of the Gaussian and isotropic GWB, and Bayesian parameter estimation to constrain GWB models. The modularity and clarity of the code allow for both a shallow learning curve and flexibility in adjusting the analysis to one’s own needs. We describe the individual modules that make up pygwb, following the traditional steps of stochastic analyses carried out within the LIGO, Virgo, and KAGRA Collaboration. We then describe the built-in pipeline that combines the different modules and validate it with both mock data and real GW data from the O3 Advanced LIGO and Virgo observing run. We successfully recover all mock data injections and reproduce published results.


INTRODUCTION
Since the first direct gravitational wave (GW) detection Abbott et al. (2016), the field of GW astrophysics has exploded, now encompassing a wide range of instrumental and observational campaigns across the globe.These detection efforts monitor a vast range of frequencies, from the nanohertz to the kilohertz, and are sensitive to a multitude of GW sources emitting therein.While the GW sources in each band may present extremely different characteristics, a potential candidate for all GW measurements is a gravitational-wave background (GWB), given by the collection of all GWs too faint to be individually resolved, or by the incoherent overlap of a large number of signals in the same band Regimbau (2011); Christensen (2019); Renzini et al. (2022).This sort of signal has been targeted in several different datasets Abbott et al. (2004Abbott et al. ( , 2007Abbott et al. ( , 2009Abbott et al. ( , 2017Abbott et al. ( , 2019Abbott et al. ( , 2021c) ) using search methods which estimate the GW strain signal power modelling the signal as stochastic, frequently resorting to cross-correlation of multiple independent observations Allen & Romano (1999).These searches are often referred to as stochastic searches by the GW detection community, and these backgrounds are often referred to as stochastic gravitational-wave backgrounds (SGWBs), even though, in practice, not all target background signals are fully described by stochastic variables1 , and this definition may imply an approximation.So far, no confident detection of a GWB has been claimed.
With this paper we present pygwb Renzini et al. (2023a), a new Python-based package tailored to searches for isotropic GWBs with current ground-based interferometers, namely the Laser Interferometer Gravitational-wave Observatory (LIGO) Aasi et al. (2015a), the Virgo observatory Acernese et al. (2014), and the KAGRA detector Akutsu et al. (2020), and with the potential to be expanded and adapted to several other detection efforts.The core analysis tools, described in detail in what follows, are heavily inspired by the LIGO, Virgo, and KAGRA Collaboration (LVK) stochastic analysis code, stochastic.m.The latter consists of a set of MATLAB scripts easily parallelizable on a high-throughput computing cluster, and has been used in LVK data analysis for the past data acquisition runs Abbott et al. (2007Abbott et al. ( , 2009Abbott et al. ( , 2017Abbott et al. ( , 2019Abbott et al. ( , 2021c)).These include the 3 observing runs: O1 (September 2015 to January 2016), O2 (November 2016 to August 2017), and O3 (April 2019 to March 2020), performed with Advanced LIGO Hanford and Livingston, and Advanced Virgo Acernese et al. (2014) for part of O2 and O3.Data from Virgo has been included in stochastic analyses as of the latest observing run.The analysis consists in the calculation of an optimal statistic Allen & Romano (1999) from the data of multiple interferometers, which is directly related to the amplitude of the GWB signal.
A notable change throughout the years of stochastic GW analyses has been the constant shift towards Bayesian parameter estimation Mandic et al. (2012); Abbott et al. (2021c).To date, there is no preferred stochastic parameter estimation software, and different groups have employed private scripts.To extend the scope of the stochastic search beyond the optimal statistic, we include a parameter estimation module in pygwb based on the Bilby package Ashton et al. (2019) which allows the user to test both predefined and user-defined models and obtain posterior distributions on the parameters of interest.
The steady inflow of ever-improving GW data open for analysis et al (LIGO Scientific Collaboration & Collaboration) has been a catalyst for open-source GW data analysis codebase development.By adopting the Python language and focusing on user-friendliness, flexibility, and portability, we intend to introduce stochastic searches to the wider GW community.Detecting a GWB with ground-based interferometers will be a community effort, and we expect search pipelines to evolve along the way.The format and structure of pygwb facilitates this evolution, and conversely, the package is suitable for beginners approaching GWB data analysis for the first time.
This paper is structured as follows.In Sec. 2, concepts related to the characterization and detection methods of a GWB are reviewed.A detailed overview of the individual modules that make up the pygwb package follows in Sec. 3, outlining the steps of LVK stochastic analyses.Several manager objects which store relevant data and handle the analysis internally are described in Sec. 4. The built-in pygwb pipeline which combines individual modules and performs the search for an isotropic GWB is presented in Sec. 5. To test the capabilities of the pipeline, mock datasets with a variety of simulated signals are analyzed in Sec.6.1.To conclude, results from the analysis of the third LVK collaboration observing run, O3, are presented and compared with collaboration results in Sec.6.2.

THE ISOTROPIC STOCHASTIC ANALYSIS
A SGWB is characterized by its spectral emission, which is the target of stochastic GW searches.The spectrum is typically parametrized by the GW fractional energy density spectrum Ω GW (f ), such that where dρ GW is the energy density of GWs in the frequency band f to f + df , and ρ c is the critical energy density in the Universe.When integrated over d log f , Ω GW (f ) gives the total dimensionless GW energy density.The Ω GW (f ) spectrum is thus directly related to the intensity of GWs.Specifically, from Eq. ( 1) it follows that Allen & Romano (1999) where the strain spectral density S h (f ) is defined as the polarization-averaged second moment of the stochastic GW strain field, decomposed into its polarization components h + and h × , assuming statistical homogeneity.The unit vectors n, n span the 2-sphere, while f ∈ R. In the plane wave formalism, h + and h × in Eq. ( 3) are the Fourier coefficients of the time-domain strain fields.If these are stochastically distributed, these give rise to a SGWB which we describe solely through the statistical moments of the distribution.In particular, a Gaussian SGWB is fully described by its second moments, hence the spectral density in Eq. ( 3) is the primary target of a search which assumes the signal to be both stochastic and Gaussian.More details on these quantities can be found for example in Romano & Cornish (2017).Laser interferometers such as LIGO and Virgo are sensitive to the strain field in the time domain coming from all directions, h(t).These detectors measure the GW strain filtered through a linear response function F (see definition in Romano & Cornish (2017)) plus a detector noise component n, which we may write in shorthand as where " " indicates a convolution operation.Given that the SGWB signal is weak and hard to distinguish from instrumental noise, cross-correlating two independent, time-coincident datastreams with uncorrelated noise is an effective way to construct an estimator for Ω GW (f ).We assume our target stochastic GW signal is stationary, Gaussian, and isotropic.We further assume the detector noise is Gaussian and uncorrelated between detectors, which is a fair assumption in the case of ground-based interferometers at current detector sensitivity2 (after specific mitigation) Abbott et al. (2021c); Janssens et al. (2023), and that the noise amplitude is much larger than the signal amplitude.Under these assumptions3 , it has been shown Aasi et al. (2015b); Romano & Cornish (2017) that the cross-correlation-based minimum-variance unbiased estimator (MVUE) of Ω GW at a frequency bin f and the corresponding variance is given by, and where C IJ,f is the one-sided cross-spectral density (CSD) and P I,f is the one-sided (auto-)power spectral density (PSD) of strain data d t from two detectors (I, J), as defined below in Sec.3.24 .Note that throughout this work we will denote continuous functions of the frequency with the notation (f ), whereas discrete functions of the frequency will be denoted with a subscript f .Typically, in this paper, discrete functions of frequency are estimators for continuous functions, and in Equations such as Eqs.( 5) and ( 6) which mix discrete and continuous functions our notation implies that continuous functions are evaluated at the discrete set of frequencies for which we know the value of the discrete functions.In the above, T is the duration of data used to produce the above spectral densities, and γ IJ (f ) is the cross-correlated GW response, or overlap reduction function (ORF), which is the polarization-and sky-averaged cross-correlation of the individual detector responses, F I .The ORF normalized for a pair of perpendicular-arm interferometers is given by Allen & Romano (1999) γ where n is the unit vector on the sky, in an arbitrary basis5 , x I − x J is the difference between the position vectors of the two detectors I and J respectively, and A spans the polarization basis.The ORF quantifies the reduction in sensitivity of the cross-correlation stochastic search due to the detectors not being co-aligned and co-located, and having different non-trivial responses.The function S 0 is defined as Romano & Cornish (2017); Renzini et al. ( 2022) and converts a GW strain power spectrum into a fractional energy density.The derivation of S 0 is shown in Allen & Romano (1999), and note that it includes the normalization factor of the ORF, 5/8π, which ensures γ IJ (f ) ≡ 1 for co-aligned, co-located detectors.
There are two important considerations to make regarding the estimator in Eqs. ( 5) and ( 6).Firstly, the implementation of a discrete Fourier transform (DFT) over a finite time T in the estimator of the continuous non-periodic quantity Ω GW (f ) may create spectral artifacts, as seen in Press et al. (2007); Whelan (2004).We outline how this is handled in Sec.3.2.Secondly, as the estimator is initially derived as a minimal variance estimator in the time domain Allen & Romano (1999), the narrow-band frequency estimator in Eq. ( 5) is actually obtained from a broad-band one, as will be clarified in Sec.3.3.In the rest of this paper, we refer to ΩGW,f as the optimal estimator of the signal spectrum Ω GW (f ).The optimality of the estimator can either be justified by the proof that this is an MVUE, or equivalently by showing that it maximizes a reasonable likelihood for the data.When performing parameter estimation as outlined in Sec.3.6, we in fact employ a Gaussian likelihood which is maximized by ΩGW,f .
In stochastic analyses with current interferometers, we take advantage of long observing times to improve detection statistics.In practice, the data are segmented into smaller chunks and analyzed individually before they are optimally combined to produce an estimate.This is convenient due to potential non-stationarities in the detector noise over both short time-scales, such as the length of an individual data segment, and long time-scales, such as the total observation time, as well as reducing computational costs.Assuming each time segment is independent, we perform a weighted average over all segments to calculate ΩGW,f for long observations.This average can be thought of as an approximation to the ensemble averages in Eq. ( 3).Hence the more independent observations are averaged over, the better the measurement.The averaging procedure is described in full in Sec.3.3.
The narrow-band statistic of Eqs. ( 5) and ( 6) assumes each frequency bin is independent.The information from each bin can be combined under the assumption of a known GW spectral density distribution.In GWB analyses, it is most common to assume a power-law spectral shape for Ω GW , where α is the spectral index of the signal, and f ref is a reference frequency, and Ω ref is defined as Under this assumption, the rescaling can be used to re-weight the estimate of the spectrum ΩGW,f , obtained for α = 0, to optimize the statistic for a specific spectral index α at a chosen reference frequency f ref , reducing the search to the estimation of a single number, Ω ref .
This procedure is referred to as re-weighting and is clarified in Sec.3.3.Alternatively, it is also possible to keep α as a free parameter in the analysis, and estimate both Ω ref and α from the data.This is described in Sec.3.6.

INDIVIDUAL MODULES
What follows is a detailed step-by-step presentation of the stochastic analysis pipeline.We follow the natural structure of the code for clarity as we introduce each module individually.To start, we present the preprocessing module which pre-conditions the time-domain strain from GW detectors for spectral analysis.In spectral, we explain the power spectrum and cross-spectrum calculations, which produce the P I,f and C IJ,f spectra in Eqs. ( 5) and ( 6).We then describe postprocessing, which includes the averaging procedures employed over large datasets to obtain an optimal estimate of the signal amplitude starting from the quantities in Eqs. ( 5) and ( 6), and knowledge of the expected spectral shape.In delta-sigma cut and notch, we present modules which focus on data quality checks, and the implementation of relevant time-domain and frequency-domain data cuts.We then describe the built-in parameter estimation module pe, based on Bilby Ashton et al. ( 2019), a Python-based Bayesian inference library widely used in GW data analysis.Finally, we present the simulator module, which includes different mock-data injection techniques for GWB study and detection validation.
A schematic of the pygwb package is presented in Fig. 1.This includes the manager objects Interferometer, Baseline, and Network, presented in Sec. 4.

preprocessing
Pre-processing is the first step of stochastic GW data analysis, in which data are read, downsampled and high-pass filtered.The pipeline can use public data available from the Gravitational-wave Open Science Center (GWOSC) et al (LIGO Scientific Collaboration & Collaboration), private data (data stored on the LVK servers restricted to members of the collaboration), or local data.Data are read using existing gwpy Macleod et al. (2021) TimeSeries methods.we denote the raw data measured at detector I over the time period T by s I (t k ) in what follows, where t k are discrete times given by t k ≡ kδt.The values of k are positive integers between 0 and T /δt − 1 and δt is the sampling period, which in LIGO, Virgo and KAGRA interferometers is 1/(16384 Hz).The raw strain data from the two interferometers, s I (t k ) and s J (t k ), are downsampled to a user-defined sampling frequency f samp , using a user-defined re-sampling window (a Hamming window by default).The downsampling is performed to reduce the memory and computational requirements of the analysis.This is achieved using an existing gwpy TimeSeries filtering method for strain data.Note that selecting an f samp implies fixing a Nyquist frequency of f Nyquist = f samp /2 for the analysis.The Nyquist frequency is the highest frequency included in the Fourier expansion at a given sampling rate.Hence, frequencies above it cannot be probed.To avoid this becoming a limitation, f samp should be chosen high enough to contain the full spectrum of the signal of interest, within reasonable sensitivity of the detector.
The low-frequency content of ground-based interferometer data (in particular below 10Hz) is dominated by seismic and control noise Buikema et al. (2020).For this reason, frequencies below a given (user-defined) cutoff frequency are high-pass filtered, i.e., excluded from the analysis.In previous isotropic GWB searches Abbott et al. (2009Abbott et al. ( , 2017Abbott et al. ( , 2019Abbott et al. ( , 2021c)), the input data are high-pass filtered using a 16th-order Butterworth filter with a specific knee frequency.A 16th-order Butterworth filter is built by first computing its transfer function (in zero-pole-gain form) using the scipy library and then filtering the data with the relevant gwpy TimeSeries method.The design of the high-pass filter is fixed in the module, only allowing the user to specify the knee frequency.The default value of the knee frequency is 11 Hz, which was chosen to avoid the spectral leakage from the noise power spectrum below 20 Hz Abbott et al. (2021c).See Fig. 2 for an example of data before and after pre-processing.
At this point, the data may also be screened for large bursts of power in the detector data with high SNR, or glitches, due to instrumental or environmental disturbances, which are known to bias estimates of stochastic analyses Usman et al. (2016); Pankow et al. (2018); Davis et al. (2021); Acernese et al. (2022b); Davis & Walker (2022).Historically, segments with loud glitches were flagged and excluded from analysis by non-stationarity cuts (see Sec. 3.4).In O3, a In blue squares, we show the manager objects of the code that handle the analysis internally.These manager objects query (red arrows) different modules for specific objects, calculations, or quantities (rounded bubbles), imported (grey arrows) by either internal (i.e., within pygwb) or external modules (i.e., outside of pygwb).Internal modules are indicated in red, while external modules are indicated in green.
series of exceptionally loud glitches appeared in the data that led to large fractions of data being removed by previously employed non-stationarity cuts Abbott et al. (2021c).Hence, an alternative technique called gating was employed to address these loud glitches, and drastically reduce the amount of data removed Matas et al. (2021).Gating is performed internally by pygwb by multiplying the data by an inverse Planck-taper window McKechan et al. (2010).Time periods around samples in the whitened data that have an absolute value above a chosen threshold are marked for gating independently for each interferometer.The width of the gate must be sufficiently large to remove the entirety of the relevant glitch.The required width may change based on the data quality of the specific data in the analysis and hence must be empirically determined.The tapering length of the window must also be sufficiently long to minimize the addition of artifacts by the gating; 0.25 seconds is found to be sufficient Davis & Walker (2022).This technique is generically beneficial for the analysis of data that are non-Gaussian, such as real gravitational-wave detector data.
Gating implemented in pygwb is highly customizable to the specific needs of the analysis; default gating parameters are shown below in Table 1.For more details on gating and parameter choices see Davis & Walker (2022).Finally, the module also allows to perform a time-shifted analysis in which one of the two timeseries is shifted in time by an integer number of seconds before the cross-correlation is performed.This technique is employed as a detector noise characterization tool, since it removes the potential correlation due to a broadband GWB, while preserving instrumental correlations with coherence times greater than the applied time shift, like nearly-sinusoidal spectral artifacts from, e.g.electronics Covas et al. (2018) 6 .The time shift is a user defined parameter which should always be greater than the light travel time between detectors (i.e., 10 ms for the LIGO Hanford and Livingston detectors) and smaller than the segment duration.Typically, a time shift of 1s is used.

spectral
The role of the spectral module is to compute, for each time segment of duration T , the discrete frequency domain quantities C IJ,f , P I,f and P J,f used in Eqs. ( 5) and ( 6).The one-sided cross-and auto-power spectral densities C IJ and P I , respectively, of a single segment are defined as where sf are DFTs of s(t k ) defined by sf ≡ T −δt where f = mδf , with m a natural number between 0 and 1/(2 δt δf ), and δf the desired frequency resolution, chosen such that 1/(2δtδf ) is an integer.The segmented data are windowed before calculating Fourier transforms to avoid spectral leakage due to discontinuities at the ends of the segments.The user may define their own choice of window, which defaults to the Hann window if none is selected.The spectral module uses methods from scipy.signal to calculate spectrograms st f of the given data, which are then used to calculate the list of C t IJ,f , P t I,f , and P t J,f quantities, corresponding to different time segments labelled by t in the dataset.By default, these are calculated with a 50% time overlap to account for the impact of the windowing.However, the user may redefine the overlap between consecutive segments to be used throughout the analysis to better suit any choice of window.
Different averaging procedures are employed to reduce the fluctuations in the spectra estimates and compress the data.The procedures we employ are selected to minimize sensitivity loss.In the estimates of P t I,f and P t J,f we employ Welch's estimation method of PSDs Welch (1967a), which is known to produce minimum variance estimates of the PSD, implemented as follows.Each segment is divided into sub-segments of duration 1/δf which are DFTed individually.The auto-correlated power |s I,f | 2 is then averaged over the sub-segments to obtain estimates of P t I,f and P t J,f for time t.This procedure returns spectra at the desired frequency resolution δf , which is typically much larger than the original resolution 1/T Hz.
As the power varies slowly with frequency7 , we can average over neighboring frequencies using a process known as coarse-graining Talbot et al. (2021).This is the default procedure employed in the CSD estimation.The resulting spectra are returned at the desired frequency resolution δf .Note that the data are zero-padded before calculating Fourier transforms for C t IJ,f to avoid wrap-around problems arising from finite data Abbott et al. (2004); Whelan (2004); Press et al. (2007), and hence coarse-graining is required to achieve the desired frequency resolution.Zeropadding simply entails appending a vector of zeros equal to the length of the segment before taking the Fourier transform.
To further reduce fluctuations in the PSD estimates, the P t I,f quantities are averaged over neighboring segments to obtain the final estimate P t I,f of the PSD at a given time t.This is appropriate as the noise in GW detectors is (most often) approximately stationary over periods of a few minutes.We often refer to the initial (un-barred) quantities as "naive" and the final (barred) quantities as "average" estimates in the rest of this paper, to avoid confusion.By default, only nearest neighbors are used for the calculation, such that the PSD at time t is an average of the naive PSDs calculated for times t − T and t + T .The user may define any even number D of segments to be used to perform this average, which are taken before and after the reference time t such that the PSD is averaged over naive PSDs at times t − DT /2 and t + DT /2.

postprocessing
Once a set of data, comprised of an uninterrupted stretch of timeseries data, has been pre-processed and average PSD and CSD estimates have been calculated for each segment of data within the set, one can combine those separate time segments to construct a final, time-averaged estimate of the GWB amplitude.
Due to the aggressive windowing choice we typically make, and the subsequent overlapping of time segments, we must be careful in combining time segments together.The overlapping and windowing cause overlapping time segments to be correlated with one another.Within each processed set, individual time segments must be combined while accounting for this covariance.A detailed calculation and discussion of this covariance can be found in Lazzarini & Romano (2004), while effective approximations to that full calculation can also be used (see, e.g.Sec.IIIB of Ain et al. (2015)).
To start, we construct the estimate of the GWB in a single segment t.As detailed in Sec. 2, the GWB search is often framed in terms of constructing a point estimate for Ω GW (f ref ), the energy density of the GWB at the specific frequency f ref , assuming a power-law for the GWB with spectral index α.We refer to the estimator of this quantity as Ωα ref , in general, and for a single time segment of data, it can be constructed using a weighted average over the individual frequency bin estimators Ωf and σ f described in Eqs. ( 5) and ( 6) calculated per segment t, as where the rescaling H ref,α (f ) is defined in Eq. ( 10).The average variance spectrum per segment, σ2 t,f , is calculated using average PSDs described in Sec.3.2.These broadband quantities can be calculated for each time segment t, and then this set of estimators at each time can be combined to account for the overlap between time segments discussed above.We first lay out how to perform this combination assuming we have calculated the quantities above for each individual time segment.Then, we discuss how to alternatively average the estimators in each frequency bin over time independently, before combining them into an integrated quantity at the end.The latter calculation is normalized such that it gives the same result as the former.To avoid heavy notation we drop the bars that indicate average quantities in the rest of this section -all variances used for the following calculations are average variances as defined above.
To construct an estimator for the GWB using a set of measurements in short, overlapping time segments, we first combine the segments that are non-overlapping.If the overlap between segments is 50% or less, then this amounts to separately performing inverse-noise-weighted averaging over the even-and odd-indexed segments: where the quantities Ω t ≡ Ωα ref,t and σ t ≡ σ α ref,t for each time segment t. Analogous expressions are calculated for Ω even and σ even .Subscripts refer to even/odd time segments, and we drop here the subscripts GW , ref , and α used to construct the integrated quantities to lighten the notation.We refer to the final, frequency-and time-averaged estimate as Ωref for now.
Next, we calculate the cross-covariance between point estimates in the odd and even segment combinations Lazzarini & Romano (2004), where M is the number of independent segments and so 2M − 1 is the total number of overlapping segments, with the window factors w4 ovl and w4 as defined in App. A. For the sake of compactness, we rewrite this as where k = w4 ovl / w4 .The covariance matrix between even/odd segment sets is then defined as which we use to construct the optimal combination of segments to obtain the point estimate Ωref and its variance σ 2 ref .
These are given by: with where i, j indices label odd/even quantities.The bias factor b avg which arises due to harsh windowing of the data has been included in Eq. ( 23).The derivation of the bias factor is described in App. A. If combining over non-overlapping segments, then σ 2 oe = 0, and this method reduces to the typical inverse-noise-weighted average that one would expect.The above expressions are for a broadband estimator, but in practice the postprocessing module combines over time segments before combining over frequency bins.We refer to the estimated narrowband quantities as Ωref,f and σ ref,f .This notation indicates that, once a power-law spectral model is applied, the estimate in a frequency bin represents an estimate of the GWB at the reference frequency of the power law, assuming the chosen spectral shape.
We normalize Ωref,f and σ ref,f such that, when performing a weighted average over frequency bins after combining overlapping time segments we get the same result as Eqs.( 22) and ( 23) (which assume construction of a broadband estimator before combining overlapping time segments).This results in the following expression for σ and a corresponding expression for Ωref,f , The even and odd estimators for each frequency bin are defined as in Eqs. ( 15) and ( 16), except applied to individual bin-by-bin estimators calculated at each time segment.As discussed above, these expressions have been normalized such that The postprocessing module implements the above expressions to estimate Ωref,f and σ ref,f at a fixed α, and returns them in form of an OmegaSpectrum object, which sub-classes the classic gwpy.FrequencySeries and adds two key attributes: the spectral index α and the reference frequency f ref at which the spectrum is calculated.By default, pygwb assumes a power-law spectral index α = 0 and a reference frequency f ref = 25 Hz when constructing the above estimators.To explicitly include the α dependence in our results, we refer to the final postprocessed spectra as Ωα derived using Eq. ( 9), which implies the following relation between amplitudes at different reference frequencies, This allows to quickly calculate time-and frequency-averaged estimates of the GWB amplitude associated with a specific power-law model.The default Hubble constant H 0 , required in the scaling S 0 (f ) in Eq. ( 8), is chosen to be H 0 = 67.7 km/(Mpc•s), drawn from the Planck 2018 observations Aghanim et al. (2020) and imported directly from the astropy package.This is an attribute of the OmegaSpectrum and may be re-set by the user.

delta-sigma cut
In general, the noise level in ground-based detectors changes slowly on time-scales of tens of minutes to hours.The variance σ 2 GW (see Eq. ( 6)) associated to each segment is an indicator of that level of noise, which typically changes at roughly the percent level from one data segment to the next.However, there are occasional very loud disturbances to the detectors, such as glitches, which violate the Gaussianity of the noise.Auto-gating procedures are in place, as explained in Sec.3.1, to remove loud glitches from the data; however the procedure does not remove all nonstationarities.To avoid biases due to these noise events, an automated technique to exclude them from the analysis has been developed Abbott et al. (2007).To this end, the pygwb package includes the delta-sigma cut module, which flags specific segments to be cut from the analyzed set.Note that inverse-noise-weighting, as explained in Sec.
3.3, also reduces the effect of non-Gaussian noise artifacts.
The "∆σ cut" calculation consists in comparing the σ GW of a segment t, σ t , to that of its nearest neighbors and flagging it for removal in case their values differ by more than a chosen threshold.Conceptually, the calculation is based on the simple inequality, where i is a segment index.However, in practice we perform an analogous, more sophisticated calculation, which compares the naive and average segment variances σ t,α and σt,α .These are derived from the unweighted naive and average segment variances computed with Eq. ( 6) using naive and average PSDs per segment (see Sec. 3.2 for details), respectively, which are then reweighted by the index α, as shown in Eq. ( 14).The final expression employed in the calculation is which also takes into account the bias factors that arise due to the different impacts of windowing on naive and average quantities (see App.A for details).Past analyses have used a threshold of 0.2, as this has been shown to yield a Gaussian distribution for the remaining (un-cut) segment variances Abbott et al. (2009).For more details on this choice see Meyers (2018).The ∆σ cut calculation is performed assuming different spectral indices α as each power law is sensitive to a different frequency band (see Fig. 5).The union of all the segments flagged for each α is taken, leading to a full list of segments to discard from the analysis.The default choice of α values in the delta-sigma cut module is α = {−5, 0, 3}, as this adequately covers most of the frequency band of LVK searches, from 20-1726 Hz Abbott et al. (2021c), at current sensitivity.These may be easily modified by the user.This would be especially recommended if the search were carried out over a different set of frequencies, or for data from detectors with a spectral sensitivity different than that for Advanced LIGO, Advanced Virgo, or KAGRA.Often, the value of α = 5 is also considered, and was employed in the most recent LVK isotropic search Abbott et al. (2021d).The analysis performed at a spectral index α = −5 is mostly sensitive to non-stationary effects in the ∼ 15 − 50 Hz range, while in the case of α = 0 the analysis is sensitive to effects between ∼ 40 − 80 Hz, for α = 3 from ∼ 90 − 500 Hz, and finally α = 5 is most sensitive to fluctuations at the higher frequencies, above ∼ 500 Hz.These higher frequencies are not always included in this sort of analysis due to reduced sensitivity in this range, hence α = 5 is not a default value used for the cut.
As the ∆σ cut only compares neighboring segments, long stretches of loud noise-contaminated data can pass the test and be included in the analysis.We are currently working to improve this by monitoring and flagging longer stretches of non-stationary noise and prolonged loud noise conditions.

notch
Ground-based laser interferometers present many narrow-frequency noise artifacts which are typically persistent in time, and are generally referred to as noise lines.Some examples are calibration lines and mechanical resonances Davis et al. (2021); Acernese et al. (2022a);van Remortel et al. (2022).The notch module provides the framework to properly deal with these noise lines in the case of the search for an isotropic GWB.The solution is to "notch out" these noise lines, i.e., set the values of the spectra at the affected frequency bins to zero.Note that the notch module is not built to identify these lines, as this is typically done by detector characterization experts working closely with instrumentalists running the detectors.Rather, the final product of the notch module is a frequency mask which may be applied to the relevant spectra in the analysis.
The key object of the notch module is the StochNotchList, which is a list of StochNotch objects.A StochNotch object represents a physical noise line which has been identified and needs to be removed from the data analysis.The object has a minimum and maximum frequency indicating the contaminated frequency region.Furthermore, it also comes with a descriptive string which allows the user to keep track of the reason why the line was notched.All the different StochNotch objects for a certain analysis are then stored in the StochNotchList which contains the entire list of lines to be notched from the analysis.
The notch mask used to apply a set of notches within the analysis is constructed conservatively, such that any frequency that has overlapping frequency content with the noise lines defined in the StochNotchList will be removed when applying the notch mask.To maintain generality, we discuss here a generic estimated spectrum Ωf , where its value at frequency f estimates Ω GW (f ) in the frequency range [f − δf /2, f + δf /2], where δf is the chosen frequency resolution, as defined in Sec.3.2.If a noise line has any overlap with the interval [f − δf /2, f + δf /2], the f frequency bin is excluded.This implies that a hypothetical delta-peak noise line at f + δf /2, leads to notching both f as well as f + δf .We present the creation of a notch mask with an example in Fig. 6, which illustrates how our conservative notching strategy excludes frequency bins based on different scenarios of noise lines.
The current code is set up to apply the same notches to an entire stretch of data, which can be considered "timeindependent" notching.To allow for time-dependent notching we could either use the current notch module and split the analysis in different segments, each having their own notch list.Alternatively, one could extend the current module with an additional parameter which keeps track of which times have to be notched.Since typically the majority of the notched lines in the search for an isotropic GWB with data from the LIGO and Virgo detectors are present during the entire dataset, the possible gain of implementing time-dependent notching is expected to be limited.

pe
Starting from an estimate of the GWB spectrum ΩGW,f , with variance σ 2 GW,f , it is possible to place stringent constraints on the GWB amplitude using a hybrid frequentist-Bayesian approach.We consider the general case where we have a set of GWB measurements ΩIJ GW,f from different detector pairs, or baselines, IJ.We define a Gaussian likelihood for B pairs of detectors, where Ω M (f |Θ) is the GWB model and Θ are its parameters.Bayes' theorem is used to obtain posterior distributions on the model parameters, where the priors p(Θ) are employed.In practice, when performing parameter estimation on a large dataset, we take the post-processed, unweighted (i.e., α = 0) estimate Ω0,IJ ref,f to be the measured GWB spectrum in each frequency bin, and plug it into Eq.( 33).Note that it is necessary for the input spectra used in parameter estimation to be unweighted as any other value would constitute a model choice and bias results.In this example there are five noise lines, from left to right: a noise line ending at the lowest frequency bin, a noise line entirely contained in one frequency bin, a noise line spread across two frequency bins, a noise line spread across multiple frequency bins, and a noise line from bin-edge to bin-edge.After our notching procedure, the data is reduced to the bins marked by the red dots.For visual convenience we have changed the amplitude in these remaining frequency bins by a factor 0.9.
Within the pygwb package, we include the pe module to perform parameter estimation as an integral part of the analysis, which naturally follows the computation of the optimal estimate of the GWB.This is a notable improvement compared to previous LVK analyses, where data products and parameter estimation were handled independently by packages in different programming languages.Furthermore, the pe module is a simple and user-friendly toolkit for any model builder to constrain their physical models with GW data.
The pe module is built on class inheritance, with GWBModel as the parent class.The methods of the parent class are functions shared between different GWB models, e.g., the likelihood formulation in Eq. ( 33), as well as the noise likelihood, given by Eq. ( 33) with Ω M (f |Θ) ≡ 0. It is possible to include calibration uncertainty by modifying the calibration epsilon parameter, which defaults to 0. For details on the marginalization over calibration uncertainty, see App.B and Whelan et al. (2014).The GW polarization used for analysis is user-defined, and defaults to standard General Relativity (GR) polarization (i.e., tensor).More details on possible polarization choices can be found in Sec.4.2.In our implementation of pe, we rely on the Bilby package Ashton et al. (2019) to perform parameter space exploration, and employ the sampler dynesty by default Speagle (2020).The user has flexibility in choosing the sampler as well as the sampler settings.
Child classes in the pe module inherit attributes and methods from the GWBModel class.Each child class represents a single GWB model, and combined they form a catalog of available GWB models that may be probed with GW data.The inheritance structure of the module makes it straightforward to expand the catalog, allowing users of the pygwb package to add their own Ω M (f |Θ) models.The flexibility of the pe module allows the user to combine several GWB models defined within the module.A particularly useful application of this is the modelling of a GWB in the presence of correlated magnetic noise, as discussed in Meyers et al. (2020), or the simultaneous estimation of astrophysical and cosmological GWBs Martinovic et al. (2021).The pygwb documentation Renzini et al. (2023b) contains information on the existing models in the catalog, with a description of the GWB models and their parameters.

simulator
To both design optimized stochastic analyses and understand our sensitivity to different categories of signals, it is essential to be able to readily simulate realistic interferometer data.To this end, the simulator module is primarily designed to generate data that corresponds to an isotropic SGWB with a given PSD.
The GWB data in a network of interferometers satisfy a specific correlation matrix, which includes the set of ORFs of the entire detector network to account for the spatial separation and relative orientation of the detectors.Given a generic signal PSD, S h , the correlation matrix C(f ) is given by Here γ IJ (f ) is the normalized ORF of the baseline IJ as shown in Eq. ( 7), hence γ II (f ) ≡ 1, and P I is the noise PSD of interferometer I.We have introduced a boldface notation which indicates matrices and vectors which span the detector space.The fact that the cross-correlation between detectors for I = J only depends on the signal PSD assumes the noise is uncorrelated across all detectors.The simulation of data correlated according to C(f ) proceeds as follows.First, a vector of white, uncorrelated frequency-domain data are generated, v f , with a certain frequency resolution ∆f .Then, the data are linearly transformed into the correlated C space by, where Λ and E are the eigenvalue and eigenvector matrices of C, respectively, calculated in each frequency bin.This transformation results in data x f that presents the correct correlation, and has been colored with the injected noise and signal power spectra, where appropriate.Finally, the frequency-domain data vector is inverse-discrete Fourier transformed (IFTed) to obtain a data vector in the time domain.Data generation in the frequency domain, followed by the IFT to the time domain, can introduce edge-effects in the simulated data segments.These may be avoided by splicing multiple data segments Rabiner & Gold (1975).The splicing procedure combines neighboring data segments by windowing and overlapping them, and thus requires more data segments than the actual desired number of segments.
Concretely, we consider the example where N seg = 1 and detail the splicing procedure below.As the number of desired segments is 1, 2N seg +1 = 3 data segments are simulated.Assuming these are simulated following the procedure outlined above, we denote these three time-domain data segments by x 0 , x 1 and x 2 .A sine window, defined as for 0 ≤ j < N , where N is the number of samples per segment, is used to window the data, which are then combined as and finally we obtain a single segment of time domain data z, In the above expressions, 0 represents an array of zeros with length N/2, used to pad the segments, whereas the bracket notation stands for python array slicing.The splicing procedure can introduce a bias in the simulated power spectrum due to the spectral properties of the window that is applied.However, the simulator module was tested for several values of the spectral index of a power-law signal PSD, ranging between −3 and 3, yielding a correct injection for these spectral indices.The user should nevertheless exercise caution when using the simulator module and be aware of the possible introduction of a bias outside the range of tested spectral indices due to splicing.
We show an example of simulated data in Fig. 7.The injected signal and noise PSDs are plotted together with the calculated PSD of a simulated data segment.A thorough testing of the simulator module is performed in Sec.6.1.1.

MANAGER OBJECTS
pygwb counts three manager objects the user can interface with: Interferometer, Baseline, and Network, which are defined in the detector, baseline, and network modules, respectively.Each object is in charge of storing and saving relevant data, and handles data analysis internally.The manager objects are designed such that the user need never call a method from a module directly, but rather will invoke the manager which queries the relevant module to perform the calculation.For details on how to use these objects, see the complete set of tutorials in the pygwb online documentation Renzini et al. (2023b).

detector
The detector module is designed to organize the data products related to a GW detector and provides functions to create and process its internal data.It is formally defined as a subclass of the Bilby Interferometer class Ashton et al. (2019).In what follows, we describe the additional features we have developed, and refer the reader to the Bilby documentation for the built-in properties inherited from the parent class.
By default, the Interferometer object is initialized by taking geometrical information of a GW detector such as latitude, longitude and elevation as arguments.
from pygwb import detector my_detector = detector.Interferometer(*args, **kwargs) While the above method allows the user to customize the detector's specification, one can initialize the object based on existing GW detectors as done in bilby's Interferometer class by calling the get empty interferometer method.
Once initialized, this object provides several ways to read in and process timeseries data, all of which internally call the preprocessing module, using information such as a channel name to query the data, or pointing to a numpy array or a gwpy object directly.Additionally, the Interferometer object includes processing methods relying on the spectral module described in Sec.3.2 to calculate naive and averaged spectrograms from the stored timeseries data.
A pair of Interferometer objects can be used to initialize a Baseline object, as described below.These are then imported as attributes of the Baseline object and store data products specific to each detector.

baseline
The Baseline module is by design the core of the pygwb stochastic analysis.Its main role is to manage the crosscorrelation between Interferometer data products, combine these into a single cross-spectrum, which represents the point estimate of the analysis, and calculate the associated error, as introduced in Sec. 2.
The standard initialization of a Baseline object simply requires a pair of Interferometer objects.
from pygwb import baseline H2H2_baseline = baseline.Baseline("H1-H2",H1, H2) Here H1 and H2 are Interferometer objects.It is also possible to load a previously stored Baseline object in pickle format by calling the relevant class method.
The data loaded into the Interferometer objects are automatically imported into the Baseline object upon initialization.The Baseline object relies on the spectral module to calculate cross-correlations between the data streams, following the methodology shown in Sec.3.2.Similarly, it relies on the postprocessing module to obtain the point estimate Ωα ref and its variance σ α ref , as described in Eqs.(22)(23).The user may choose to calculate point estimate and sigma spectra or point values; in the latter case, the spectra are automatically stored to facilitate subsequent analyses.
Calculating Ωα ref , as well as performing parameter estimation on the GWB spectrum, requires the two-detector ORF, γ IJ , shown in Eq. ( 7).The ORF is calculated at Baseline object initialization, then stored as an attribute.By default, we assume GR, which presents two independent degrees of freedom for the strain field, typically A = {+, ×} in the transverse-traceless gauge.For a precise derivation of this function and detector response definitions, see for example Romano & Cornish (2017).
The Baseline object is also equipped to probe circularly polarized backgrounds Seto & Taruya (2007), and non-GR polarizations in the GWB, such as scalar and vector backgrounds Callister et al. (2017).This requires selecting a different choice of A, according to the chosen polarization type, which can be declared when calculating Ωα ref or the ORF directly.Details on the expressions for non-GR γ IJ functions may be found in the appendix of Callister et al. (2017).

network
The network module is designed to handle two different tasks.Its primary purpose is to combine results from different Baseline objects.Similarly to the Baseline object, the Network object imports Baseline objects as attributes which may be invoked through the Network.In addition to this functionality, it can also be used to simulate cross-correlated data across a network of detectors.Both signal-only and signal and noise data can be simulated using a Network object.The network module handles all data generation by querying the simulator module.
The Network object can be initialized in two ways.By default it is initialized through a list of Interferometer objects.The combined point estimate and sigma spectra are stored as attributes of the Network.These are combined by performing an inverse-noise-weighted average over the individual Baseline final spectra, assuming the data are uncorrelated between baselines, i.e., assuming each baseline provides independent information.This is a valid approximation when working in the large noise limit.Further details can be found in Allen & Romano (1999).
The Network is also designed to produce appropriately correlated simulated data for a network of interferometers.The Network can either simulate data from scratch, or add simulated data to pre-existing data, if the interferometers used to initialize the object contain strain data.The latter is simply done as strain adds coherently in the time domain.This functionality relies on the simulator module which performs the data simulation, as discussed in Sec.3.7.

ANALYSIS PIPELINE
The previous sections contain a detailed description of each of the modules of the pygwb package.We now present an overview of the package analysis scripts, which combine the various modules into a GWB analysis pipeline.The pipeline has several default values which may be changed according to the user's requirements.However, we note that thanks to the flexibility of the pygwb package, one can also easily construct an ad-hoc pipeline.

pygwb pipe
The core script of our analysis suite, pygwb pipe, is designed to carry out the bulk of the stochastic analysis.It combines the pygwb modules in order to go from the unprocessed data to the optimally averaged Ωα spectra for a single baseline.To read in the analysis parameters, pygwb pipe interfaces with the parameters module, specifically designed to handle the analysis parameters, either passed through an initialization file (param file) or declared in the command line.The module includes the Parameters dataclass which stores the chosen parameters.The pipeline may be run from the command line as follows.
pygwb_pipe --param_file {path_to_param_file} All param file parameters may be alternatively passed from the command line directly.If a mixture of parameter file and command line parameters are passed, the latter will override their corresponding values stored in the parameter file.Additionally, a set of pipeline-specific parameters may be passed from the command line for ease of use, such as whether to apply data quality cuts.A full list of parameters and their description may be found in Table 1.
After reading in the parameters, two Interferometer objects are created accordingly, and data are loaded in and preprocessed using the preprocessing module.Depending on the value of the gate data parameter in the initialization file, the gating outlined in Sec.3.1 also takes place at this stage.Subsequently, a baseline object is created using the pair of interferometer objects.Recall that the baseline module plays a central role in the pipeline and handles the computation of the various quantities interest, including the (average) PSDs and CSDs of the baseline, relying on the spectral module.This is described in more detail in Sec.s 3.2 and 4.2.
The delta-sigma cut is then performed, and optimally averaged spectra and overall point estimate are calculated with the relevant Baseline methods.The delta-sigma cut is applied by default, but may also be calculated and applied at a later stage.Finally, the spectra, the overall point estimate, and the pickled baselines (if requested), are saved as output.By default, the output is in numpy binary file format, npz.
In realistic scenarios, we analyze year-long datasets and running pygwb pipe in series is sub-optimal.However, a long dataset can be split into smaller jobs and parallelized on a cluster.The output of each job is then combined into a single set of result spectra Ωα ref,f and σ ref,f using the pygwb combine script.The latter simply takes a weighted average over all jobs, assuming each job is an independent measurement of the signal.At this stage it is possible to implement final post-processing choices, such as re-weighting the spectra to a desired α and f ref , as well as change the default Hubble constant H 0 at which results are reported.
Details on running the pipeline and combination scripts may be found in our online documentation Renzini et al. (2023b).

statistical checks
With the statistical checks module, we provide a tool to perform initial statistical analyses of a pygwb run result set, and visualize them in pre-formatted plots.We identify five broad categories of checks.
The first set calculates the running point estimate for Ωα ref and σ α ref quantities as a function of time, as more data segments are added to the analysis.The values of α and f ref are those used in the analysis and may not be changed at this point.The running averages are cumulative weighted averages of time-ordered segments, and do not take segment-by-segment correlation into account.In case of detection, these converge to a biased point estimate and σ, as proper postprocessing is not applied (see Sec. 3.3).However, the visualization of running quantities is extremely useful to identify trends in the data, and ultimately will flag a possible detection.The module also provides a linear trend analysis, fitting the evolution of the parameters described above as a function of time.
The second set focuses on the signal-to-noise ratio (SNR) spectrum as a function of frequency, defined as The absolute value, real, and imaginary part of the SNR are calculated, as well as the cumulative SNR.An example of these plots using the first sub-set of O3 data further described in Sec.6.2 is given in Fig. 8.These plots are a faithful representation of the "noisiness" of each frequency bin and how much each bin contributes to the analysis.The third set of checks produces the IFT of the point estimate spectrum, which should peak around zero seconds in case of a detection.Time-shifting the data in two detectors by more than the coherence time between the two detectors breaks the coherence between the two data streams, removing any evidence of a GWB signal.Note that the coherence time is determined by the bandwidth of our signal, which is of order 100 Hz, resulting in a coherence time of 10 ms.Hence, a GWB signal will only peak around zero time lag between the output of the detectors.
The fourth set studies the effect of the ∆σ data quality cut described in Sec.3.4 on the analysis run.To this end, we display several quantities before and after the cut is applied to the data, including the segment values of Ωα ref,i , σ α ref,i , and ∆σ α ref,i , and the deviations in SNR,  as a function of time.Here angle brackets indicate an arithmetic mean over all segments i.We also plot a histogram of the values of ∆SNR before and after the cut.This distribution should be centred around 0, with a smooth narrower distribution after the application of the ∆σ cut.We additionally plot the ∆SNR as a function of individual σ α ref,i .Finally, we plot the distribution of the ratios σ 2 ref,i / σ 2 ref,i , which should peak around 1. Some representative plots are shown as an example in Fig. 9.
The last set of checks concerns a Kolmogorov-Smirnov (KS) test that is used to verify that the ∆SNR i are consistent with a Gaussian distribution.The KS test implementation of this module returns the KS test statistic, which is the maximal deviation from the Gaussian cumulative distribution function, as well as the p-value.These values can be used to make statements about the Gaussianity of the data Dodge (2008).In addition, the cumulative distribution function is plotted for the data as well as for a Gaussian distribution.

TESTING
To comprehensively test the pygwb analysis suite, we employ an efficient workflow to analyze datasets of increasing complexity.The datasets considered in this paper are: 1. Continuous SGWB: A loud stationary and continuous stochastic signal generated with the simulator module, injected in Advanced LIGO Hanford and LIGO Livingston assuming design A+ sensitivity Barsotti et al. (2018).
2. Realistic compact binary coalescence (CBC) GWB: A realistic background of merging binary black holes (BBHs) and binary neutron stars (BNSs), injected in Advanced LIGO Hanford and LIGO Livingston assuming design A+ sensitivity.

O3 dataset:
The full Advanced LIGO Hanford and LIGO Livingston dataset from the third LVK observing run Abbott et al. (2021c).
The continuous SGWB (dataset 1) is an idealized observing scenario, as our stochastic model matches the target signal perfectly by design, and as the signal is stationary and continuous our approach is optimal Drasco & Flanagan (2003); Lawrence et al. (2023).The CBC background (dataset 2) is a realistic scenario where the target signal is generated according to astrophysical models, informed by GW detections.In this case the signal is non-Gaussian, and we expect our approach to be un-biased Meacher et al. (2015); Regimbau et al. (2012Regimbau et al. ( , 2014) ) but sub-optimal Drasco & Flanagan (2003); Lawrence et al. (2023), due to the intermittent nature of the signal which is not taken into account in the search method.For more details on the time-domain characteristics of these two types of signals and the detection challenges these present, see for example Regimbau (2022).Finally, the O3 Advanced LIGO dataset (dataset 3) presents all the complexity of analyzing real GW detector data, which includes non-stationary noise, a large data volume, and expensive computational requirements.We handle large datasets by splitting the data into smaller pygwb pipe jobs, assuming each job is independent; these are then combined using the pygwb combine script (see Sec. 5 for details).We then employ a parameter estimation script, pygwb pe, based on the pe module described in Sec.3.6, to perform parameter estimation on specific models.For more details on how to run this sort of analysis, we refer users to the online documentation for the most up-to-date workflow instructions Renzini et al. (2023b).In the following, we present the different datasets and summarize our analysis results.

Stationary and continuous stochastic gravitational-wave background
We employ the Network (Sec.4.3) to generate a stationary and continuous SGWB signal with a fixed PSD, S h (f ).This allows us to simultaneously test the module and the whole analysis pipeline.The injected SGWB is scale-invariant, i.e., Ω GW (f ) is constant over frequencies, This is converted to S h (f ) using the relation in Eq. ( 2).The noise P n (f ) is taken to be Gaussian, colored using the the Advanced LIGO noise PSD Aasi et al. (2015a).One hundred days of consecutive data are simulated at a sampling rate of 1024 Hz.Each of the one hundred days is analyzed separately, and we recover a distribution of Ω0 25 point estimates, shown in Fig. 10  construct a distribution of recovered point estimates, which is useful to assess the ability of the simulator module to inject a stochastic stationary signal.We then perform parameter estimation on the combined one hundred days, presented in Fig. 10 (right).We assume a log-uniform prior from 10 −11 − 10 −6 for Ω ref and Gaussian prior with mean 0 and standard deviation 1.5 for α.This shows a recovery within 1σ for Ω ref = 1.06 × 10 −7 and within 2σ for the spectral index α inj = 0.
The tests above illustrate that the simulator module is able to successfully inject a stochastic stationary signal and that the pygwb pipeline is able to recover this injection.

Gravitational-wave background from a coalescing compact binary population
The inspiral and merger of two compact objects emit a characteristic GW signal.We generate datasets containing a GWB signal resulting from the superposition of GW signals from a set of CBC populations including BBHs and BNSs.To simulate the signals, we employ the code used in the past for the Einstein Telescope (ET) mock data and science challenges (MDSCs) (Regimbau et al. 2012(Regimbau et al. , 2014;;Meacher et al. 2016) and for the Advanced LIGO and Advanced Virgo MDSC (Meacher et al. 2015).The Monte Carlo algorithm that we use for the generation of a compact binary population up to redshift z = 10 is extensively described in Regimbau et al. (2012) and Regimbau et al. (2015).We summarize below the main steps of the simulations.
To generate a CBC population we assume a merger rate per unit redshift (Belczynski et al. 2006;Berger et al. 2007;Belczynski & Kalogera 2001;Bulik et al. 2004), where dV c /dz is the co-moving volume element and r c the coalescence rate as a function of redshift Regimbau (2011).
The element of co-moving volume assumes a ΛCDM cosmology from Planck 2018 Aghanim et al. (2020) (Hubble parameter H 0 = 67.7 km/s/Mpc, Ω m = 0.31 and Ω Λ = 1 − Ω m ).We assume a coalescence rate normalized to a local rate r c (0) = 1 Mpc −3 Myr −1 for BNS coalescences and r c (0) = 3 Mpc −3 Myr −1 for BBH coalescences, assuming the star formation rate from Hopkins & Beacom (2006) and a minimum delay time between binary formation and merger of 20 Myr for BNSs and 50 Myr for BBHs; see (Dominik et al. 2012;Neijssel et al. 2019) for more details.These choices give rise to a dataset composed by 87% of BNSs and 13% of BBHs.
The time intervals τ between consecutive CBC events in our population are obtained by sampling an exponential distribution P (τ ) = exp(−τ /τ ), where τ is the average time between consecutive events.This is consistent with the assumption that the coalescence times t c of the events behave as a Poisson process Regimbau et al. (2015).The coalescence redshift is drawn from the normalized coalescence rate p(z) = τ dR/dz(z) within z ∈ [0, 10].The sky position n of each source is generated isotropically on the sky.The GW polarization angle ψ, the phase angle φ 0 at the coalescence time, and the cosine of the inclination angle of the orbital plane to the line of sight ι are all drawn from uniform distributions.The mass function of the components in the BBHs is chosen to be a power-law-plus-peak (PLPP) from the preferred case presented in the LVK collaboration CBC population inference paper (Abbott et al. 2021b) or a simple power-law (PL) (Abbott et al. 2021a), while the BNS masses are drawn from uniform distribution between 1 and 3 M .The BBH mass functions are used to label the two datasets presented below.Spins are neglected in both cases.
For each source, the signal waveform is generated in the time domain.For BNSs, we use the TaylorT4 time-domain waveform (Buonanno et al. 2003).For BBH signals, we use the EOBNRv2 (Buonanno et al. 2009) time-domain waveform from numerical relativity.These are then injected into the LIGO Hanford and Livingston detectors, with the addition of colored Gaussian noise generated from the LIGO A+ Design O'Reilly et al. (2020a); Barsotti et al. (2018) expected sensitivity curve, to produce the final datasets.
Following the above prescriptions, we generate two six-months datasets with sampling frequency 1024 Hz, labelled PLPP and PL, formed by two different CBC populations.The two populations differ by the mass distributions of the BBHs and the average time between consecutive events, as seen in Table 2.The latter is chosen such that the GWB amplitudes of the two datasets match, for ease of comparison.Furthermore, to obtain a SNR large enough to confidently detect the injected GWB, a small amplification of the signal is required.To this end, the amplitude of the CBC waveforms is multiplied by 1.5 and 1.7 for the PLPP and the PL datasets, respectively, resulting in an injected value of Ω ref = 2.05 × 10 −9 .
The Ω GW (f ) spectrum relative to the each dataset is obtained by summing the contributions from individual coalescences (Meacher et al. 2015), and is illustrated in Fig. 11.As may be observed, in the case of CBC signals Ω GW (f ) increases as f 2/3 from the inspiral phase (and then as f 5/3 from the BBH merger phase) before reaching a peak and steeply decreasing Marassi et al. (2011).This motivates fixing the spectral index parameter to α = 2/3 in our searches.Fig. 11 also shows the power-law integrated sensitivity (PI) curve (Thrane & Romano 2013) for the Hanford-Livingston baseline, assuming the design A+ sensitivity for the two detectors (Barsotti et al. 2018), an observation time T obs = 6 months, and a desired sensitivity of SNR=5.Given that the PI curve is almost tangent to Ω ref of the two datasets, we expect to observe the GWB signals with SNR ∼ 5.
We analyze the two datasets in the frequency band 20 − 500 Hz, using a frequency resolution of 1/32 Hz and a segment duration of 192 s.We choose α = 2/3, f ref = 25 Hz, and H0 = 67.7 km/s/Mpc for this analysis.The results of the analysis are summarized in Table 2.We recover the PLPP injection within 1 σ, and observe it with SNR = 5.4, while recovering the PL injection within 1 σ, with SNR = 5.0.We attribute the differences in the recoveries to the specific data and noise realizations within the datasets.We then proceed with estimating the parameters α and Ω α ref modelling Ω GW (f ) as a simple power-law in frequency as given by Eq. ( 9).We assume a log-uniform prior over Ω 0 ref in the range [Ω 0 min , Ω 0 max ] = [10 −11 , 10 −8 ], and a Gaussian prior on α with mean 2/3 and standard deviation (log 10 Ω 0 min − log 10 Ω 0 min )/2 = 1.5.Note that the priors in Ω α ref are defined for α = 0.The choice of the prior over α can be understood as follows.The log-uniform prior over Ω 0 ref induces some implicit prior over α that can be shown to be a triangular prior centred on α = 0 and non-zero for |α| ≤ (log 10 Ω 0 max − log 10 Ω 0 min ).To avoid a vanishing prior outside of this range, we choose a Gaussian prior for α with standard deviation comparable with the triangular prior, centered on α = 2/3 to better match the injected GWB.
Parameter estimation corner plots are shown in Fig. 12.For both datasets, Ω ref and α are recovered within 1 σ.

O3
In this section we present results from the application of the pygwb analysis suite to the full LIGO Hanford and LIGO Livingston O3 dataset et al (LIGO Scientific Collaboration & Collaboration).We set upper limits on the signal from a SGWB and confirm these are consistent with previously published collaboration results Abbott et al. (2021c).
The O3 data run collected between April 1, 2019 and March 27, 2020, divided into two sub-sets with an interruption between October 1 and November 1, 2019, with a total coincident livetime of 205.4 days between LIGO Hanford  2021c) within 1%, with the previous cut excluding an extra 0.06% of the time.We believe this small variation to be due to a different window bias factor used in the two analyses (the bias factor calculation used here is outlined in App.A).
We calculate broad-band integrated estimates between 20 − 1726 Hz of Ω GW (f ref ) for different power-law spectral models, applying the released O3 notchlist LIGO Scientific Collaboration, Virgo Collaboration and KAGRA Collaboration (2021) to exclude known problematic frequencies Covas et al. (2018).A summary of the values for the point estimate and uncertainty for these is presented in Table 3.The uncertainties σ α ref agree within 1% with previously published LVK results, presented in Abbott et al. (2021c).The point estimates for Ω α ref fluctuate notably more than the uncertainties.We believe this to be due to small differences in the analyses, to which the point estimates are more sensitive, such as individual start and end time of each pipeline job, and the differences in the non-stationarity cuts described above.
Finally, we perform parameter estimation to constrain Ω GW (f ref = 25Hz) ≡ Ω 25 and the spectral index α with O3 data.We employ a log-uniform prior on Ω 25 spanning [10 −13 , 10 −6 ], and present results for two different priors on α: a uniform prior between [−4, 4] and a Gaussian prior centered around 0 with norm 3.5 (the latter matches the choice in Abbott et al. (2021c)).To account for calibration uncertainty, we marginalize over the uncertainty parameter λ as described in App.B, with combined calibration error for Hanford and Livingston of 1.48%, as in Abbott et al. (2021c).Parameter estimation confirms Ω 25 is consistent with 0 and α remains unconstrained, as may be seen in Fig. 14.These results agree with the previous parameter estimation carried out in Abbott et al. (2021c).
Our results are quoted at the value of the Hubble parameter H 0 = 67.9km/s/Mpc, in line with published results.This is not the built-in value of H 0 , defined in Sec.3.3; however, rescaling is straightforward as it is an overall multiplication factor, which may be changed when post-processing the run with the pygwb combine script or manually using the built-in functions of the OmegaSpectrum object, as explained in Sec.s 3.3 and 5.
We would like to note that this entire analysis was carried out on a large computing cluster and completed in less than five hours of human time.This is an example of the computational efficiency of our package.3 −1.5 ± 0.9 −1.3 ± 0.9 Table 3. Summary of pygwb search results on O3 dataset.On the left, three columns summarising point estimates from the weighted optimal statistic, at different spectral indices α.On the right, three columns summarising Bayesian upper limits (UL) with log-uniform prior on Ω 0 25 and either uniform or Gaussian prior on α.These results are consistent with no detection of the amplitude of the background (Ω α ref is consistent with 0), nor its spectral shape (α remains unconstrained).

CONCLUSIONS
We present a new Python-based package tailored to GWB searches with current ground-based interferometers.We opt for a modular code, where each module performs specific tasks of the GWB data analysis.The modularity of the code results in large flexibility and offers the possibility to customize the pipeline according to one's own needs.With the use of Python language, the user-friendliness and flexibility of the code, we aim to bring GWB searches to the wider GW community, as the detection of a GWB with ground-based interferometers draws potentially closer.With increasing amounts of GW data, pygwb also answers the need for an open-source GWB data analysis tool.
In this paper, we show the application of the pygwb package to mock datasets, illustrating how the various modules can be assembled to form a search pipeline, and showing what a GWB detection could look like with our analysis approach.To conclude, we run the pygwb pipeline on real GW data from the third observing run (O3) of the LVK collaboration, and recover results in agreement with published results.Both analyses serve as a validation of the software.
The pygwb package is designed to evolve along the way and address new analysis needs as they arise.This is facilitated by the structure and the format of pygwb, and the management of the online Git repository.The pygwb team invites input from the broader community, under the form of Git issues and pull requests.New contributors to the code are always welcome.Official updates and releases of the code will be handled and reviewed internally by the software and review teams, which are due to evolve.
We are aware other analysis methodologies exist which accommodate different features of specific GWBs, such as potential anisotropy Ain et al. (2018), and the intermittency of the BBH background Smith & Thrane (2018); Lawrence et al. (2023).We look forward to interfacing with these methods and, where useful and appropriate, improving the current codebase to support and encompass more analysis schemes.
Finally, we are particularly excited at the prospect of broadening the scope of the package to include support for next generation detectors such as ET Maggiore et al. (2020) and Cosmic Explorer (CE) Reitze et al. (2019).While the science cases and design properties of these detectors are still under development, there is evident interest in targeting GWBs with these detectors within the community Regimbau et al. (2012Regimbau et al. ( , 2014)); Sathyaprakash et al. (2011), and a notable increase in sensitivity compared to present-day interferometers is expected.
where we apply a Hann window with amplitude {w i } at each sample i, as well as the overlapping of our chunks of data.The correction factor is given by Welch (1967b) In practice, we ignore the term (K − 1)/K, as it leads to extra corrections that are O(K −2 ) that are quite small.We can now define a bias correction factor based on the windowing we choose and the number of averages used in constructing PI,f .Defining N eff = κ −1 K, we have where we have used simplified notation again where the hat indicates our estimator for Eq. ( 6) and the unhatted indicates the true value.
Taking the square root of both sides and inverting it gives us where the bias factor, b(N eff ), is given by assuming N eff is large.In Sec.3.4, two different bias factors are discussed.In one case, the "naive" σ is estimated using one segment of length T , which results in fewer effective averages, and a larger bias correction than our typical estimate of σ which uses two adjacent segments of length T and there twice as many averages.

B. MARGINALIZING OVER CALIBRATION UNCERTAINTY
Given measurements { Ωi } with uncertainties σ 2 i , as shown in Sec.3.6 the following likelihood function can be used to perform parameter estimation on the GWB: Here, the { Ωf } are a set of estimators for the GW energy density at discrete frequencies f , Ω M (f |Θ) is a model for the energy density with parameters Θ, and N is a normalization constant.We will consider only a single baseline and neglect the sum over detector pairs IJ appearing in Eq. ( 33); if multiple detector pairs exist, the derivation below can be replicated for each pair.Eq. (B10) assumes that our estimators { Ωf } are direct, unbiased measurements of the underlying energy-density spectrum.In general, however, the imperfect amplitude and phase calibration of GW detectors will break this assumption.We can account for calibration uncertainty by amending our likelihood to introduce a new parameter λ: The parameter λ is an unknown multiplicative factor that encapsulates potential calibration inaccuracy.In the case of perfect amplitude calibration (λ = 1), then { Ωf } are direct measurements of the underlying (unknown) energy spectrum.But if our calibration is imperfect (λ = 1), then { Ωf } are instead measurements of some multiple λΩ(f ) of the GWB spectrum.Although we do not know λ, it is possible to estimate the uncertainty on instrumental calibration.We will therefore model λ itself as an unknown variable drawn from a normal distribution centered at 1 (corresponding to perfect calibration) but with a variance 2 : where is a known amplitude calibration uncertainty.Additionally, we impose the constraint that λ be positive: we expect errors in the amplitude of strain measurements but not their sign.In this case, the probability distribution for λ becomes normalized to unity on the interval λ ∈ (0, ∞).Eq. ( B13) is our prior on λ.
We can now use Eq.(B13) to marginalize our likelihood (Eq.( B11)) over the unknown calibration factor λ. The marginalized likelihood is given by p({ Ωf }|Θ) = p({ Ωf }|Θ, λ) p(λ)dλ (B14) If we define and ) Marginalization of calibration uncertainty is built into the pygwb pe module, and this calculation is automatically triggered when passing a calibration error = 0. Additional information on the treatment of calibration uncertainties can be found in Whelan et al. (2014).

Figure 1 .
Figure1.Schematic overview of pygwb analysis flow.In blue squares, we show the manager objects of the code that handle the analysis internally.These manager objects query (red arrows) different modules for specific objects, calculations, or quantities (rounded bubbles), imported (grey arrows) by either internal (i.e., within pygwb) or external modules (i.e., outside of pygwb).Internal modules are indicated in red, while external modules are indicated in green.

Figure 2 .
Figure 2. Comparison between the amplitude spectral density (ASD) of a raw (blue solid line) and pre-processed (orange solid line) 192 s segment of LIGO Livingston O3 data.Pre-processing consists of downsampling the data to 4096 Hz and then removing the low frequency content below 10 Hz.

Figure 3 .
Figure 3.An example of cross-and auto-power spectral densities of the Hanford and Livingston detector data during O3.

Fig. 3
shows the cross-and auto-power PSDs of 192 s of data from the Hanford and Livingston detectors during O3, while Fig.4shows a two-hour spectrogram of Hanford data during O3, produced with the spectral module.

Figure 4 .
Figure 4.An example spectrogram showing two hours of LIGO Hanford data during O3.The visible vertical columns correspond to noisy segments, which are usually removed from the analysis (see Sec. 3.4).
ref,f and σ α ref,f .One of the advantages of averaging over time before averaging over frequency is that one can reweight Ωα ref,f and σ α ref,f to be estimators for different choices of α without needing to average over all time segments again for a new choice of α.The OmegaSpectrum object has a built-in method to perform a reweighting to change either f ref or α used to calculate Ω ref , employing the relation

Figure 5 .
Figure 5.In this plot, power-law spectra with different spectral indices are compared to the O3 sensitivity curve of LIGO-Livingston.Each power law is sensitive to a different frequency band.This makes it necessary to repeat the ∆σ cut assuming different α, since this allows to check for noise fluctuations in the whole range of frequencies analyzed.The O3 sensitivity curve for LIGO-Livingston was retrieved from O'Reilly et al. (2020b).

Figure 6 .
Figure6.Example of how the notching of noise lines (orange curve) applied to the discrete measurements of the spectrum ΩGW,f (blue stars) leads to a final set of measurements (red dots).The vertical shaded regions indicate the bins, where even bins are white and odd bins are light blue.The orange line traces out the noise lines such that a noise line is present where the orange curve is zero.The analyzed data spans [5.0, 6.875] Hz, in the un-shaded region.In this example there are five noise lines, from left to right: a noise line ending at the lowest frequency bin, a noise line entirely contained in one frequency bin, a noise line spread across two frequency bins, a noise line spread across multiple frequency bins, and a noise line from bin-edge to bin-edge.After our notching procedure, the data is reduced to the bins marked by the red dots.For visual convenience we have changed the amplitude in these remaining frequency bins by a factor 0.9.

Figure 7 .
Figure 7. Example of injection using the simulator module.The noise PSD (orange) and the injected signal (red) are clearly discernible as a part of the PSD of the simulated data (blue).

Figure 8 .
Figure 8. Left: Absolute value of the SNR spectrum as a function of frequencies.Right: Sigma spectrum as a function of frequency.

Figure 9 .
Figure 9. Left: Point estimate, sigma and deviates ∆SNRi as a function of time before the delta-sigma cut (red) and after the cut (blue).Right: Distribution of the deviates ∆SNRi as a function of time before the delta-sigma cut (red) and after the cut (blue).

Figure 10 .
Figure 10.Left: Distribution of the recovered point estimate for each day in the dataset.The injected value is denoted by the red line.Ω ref and σ ref denote the mean and the standard deviation of the one hundred point estimates.Right: Parameter estimation performed on the one hundred days, obtained assuming a log-uniform prior from 10 −11 − 10 −6 for Ω ref and Gaussian prior with mean 0 and standard deviation 1.5 for α.The injected values are denoted by the black lines, while the contours represent the 1σ, 2σ, and 3σ contours.The vertical dashed lines represent the 1σ confidence interval.
Figure11.ΩGW(f ) for the dataset corresponding to the PLPP (blue) and the PL (red) models for the mass function.The black line is the power-law integrated sensitivity (PI) curve for an observation time of six months and an expected SNR = 5, assuming the HL baseline with Advanced LIGO plus design sensitivity.The simulated signals intersect the PI curve, hence they are expected to be detected with an SNR of at least 5.

Figure 12 .
Figure 12.PE Left: Corner plot obtained from running the parameter estimation over the PLPP dataset.Right: Corner plot obtained from running the parameter estimation over the PL dataset.Each plot shows the posteriors on Ω ref and α obtained assuming a log-uniform prior on Ω 0 ref from 10 −11 -10 −8 and a Gaussian prior on α with mean 2/3 and standard deviation of 1.5, respectively, denoted by the gray dashed lines.The injected values are represented by the black lines, indicating a recovery of both the amplitude of the signal and α within 1 σ.The vertical blue dashed lines represent the 2σ confidence interval.

Figure 13 .
Figure 13.Estimated cross-correlation spectrum Ω025 ± σ0 25 from O3 data.By eye, it is possible to spot several narrowband artifacts (lines) which are subsequently excluded from our analysis.

Figure 14 .
Figure14.Parameter estimation results with pygwb pe on LVK O3 data, using a log-uniform prior on Ω25, and a uniform prior (left) or a Gaussian prior on α (right), as described in the text.The priors are denoted by the gray dashed lines.

Table 1 .
Default parameters for the pygwb pipe script as well as the Parameters dataclass.Most of these choices reflect defaults chosen in the past when analysing LIGO and Virgo data.Notably, the default start and end times for the analysis are not meaningful and represent placeholders for the user-defined times.A default initialization file is included in the package with meaningful start and end times present in the O3 open dataset.

Table 2 .
Kass & Raftery (1995) B GW noise are 11.1 and 9.2 for the PLPP and PL datasets, respectively, indicating strong evidenceKass & Raftery (1995)for the presence of signal over noise only.Parameters and results of each dataset.The first row refers to the PLPP dataset, while the second row to the PL one.The second and third columns display the average time between two successive binary mergers, τ , and the waveform amplification factor, a.The last three columns illustrate the recovered point estimate with 1 σ uncertainty on the quantity Ωα ref = 25 Hz, α = 2/3), the corresponding SNR, and the log-Bayes factor B GW noise .
ref (f