Upper Limits on the Isotropic Gravitational-Wave Background from the first part of LIGO, Virgo, and KAGRA's fourth Observing Run

TL;DR

This paper reports upper limits on the isotropic gravitational-wave background using data from LIGO, Virgo, and KAGRA's first part of the fourth observing run, constraining theoretical models and merger rates.

Contribution

It provides the most stringent constraints to date on the gravitational-wave background from compact binary coalescences, improving sensitivity and combining background limits with GW event data.

Findings

No statistically significant background detected.

Constraints on energy density $ ext{Ω}_{ m GW}$ improved by factors of 1.7 and 2.1.

Estimated background from binary black holes is $ ext{Ω}_{ m BBH} imes 10^{-9}$.

Abstract

We present results from the search for an isotropic gravitational-wave background using Advanced LIGO and Advanced Virgo data from O1 through O4a, the first part of the fourth observing run. This background is the accumulated signal from unresolved sources throughout cosmic history and encodes information about the merger history of compact binaries throughout the Universe, as well as exotic physics and potentially primordial processes from the early cosmos. Our cross-correlation analysis reveals no statistically significant background signal, enabling us to constrain several theoretical scenarios. For compact binary coalescences which approximately follow a 2/3 power-law spectrum, we constrain the fractional energy density to $\Omega_{\rm GW}(25{\rm Hz})\leq 2.0\times 10^{-9}$ (95% cred.), a factor of 1.7 improvement over previous results. Scale-invariant backgrounds are constrained to $\Omega_{\rm GW}(25{\rm Hz})\leq 2.8\times 10^{-9}$, representing a 2.1x sensitivity gain. We also place new limits on gravity theories predicting non-standard polarization modes and confirm that terrestrial magnetic noise sources remain below detection threshold. Combining these spectral limits with population models for GWTC-4, the latest gravitational-wave event catalog, we find our constraints remain above predicted merger backgrounds but are approaching detectability. The joint analysis combining the background limits shown here with the GWTC-4 catalog enables improved inference of the binary black hole merger rate evolution across cosmic time. Employing GWTC-4 inference results and standard modeling choices, we estimate that the total background arising from compact binary coalescences is $\Omega_{\rm CBC}(25{\rm Hz})={0.9^{+1.1}_{-0.5}\times 10^{-9}}$ at 90% confidence, where the largest contribution is due to binary black holes only, $\Omega_{\rm BBH}(25{\rm Hz})=0.8^{+1.1}_{-0.5}\times 10^{-9}$.

Tables

Table 1. Table 1: Maximum calibration uncertainties on the magnitude and phase of detector strain at 68% confidence level (CL), for the LIGO Hanford (LHO) and LIGO Livingston (LLO) detectors.

Interferometer	LHO	LLO
Magnitude	$6.9\%$	$4.1\%$
Phase	${4.3}^{\circ }$	${5.1}^{\circ }$

Table 2. Table 2: Search results for an isotropic GWB using cross-correlation optimal filtering for power-law GWBs with α = 0 , 2 / 3 , 3 \alpha={0,2/3,3} . For each case, we report the point estimate and associated 1 ​ σ 1\sigma uncertainty of the amplitude Ω ^ α \hat{\Omega}_{\alpha} at 25 ​ Hz 25\,{\rm Hz} , along with the frequency range from 20 ​ Hz 20\,{\rm Hz} to f 99 % f_{99\%} , which contains 99% of the total sensitivity. These estimates are derived by combining the individual values from independent frequency bins following Eq. ( 11 ). The final two columns show results obtained by combining data from all observing runs (O1 − - O4a) following the prescription shown in  Abbott  et al. ( 2021b ) ( Eq. (7)). Note that Virgo data contributed to the O3 measurements, implying three separate baseline contributions to the O3 results.

Power law	$f_{99\%}^{\mathrm{O4a}}\left[\mathrm{Hz}\right]$	${\hat{\Omega }}^{\mathrm{O4a}}/{10}^{-9}$	$f_{99\%}^{\mathrm{O1}-\mathrm{O4a}}$ [Hz]	${\hat{\Omega }}^{\mathrm{O1}-\mathrm{O4a}}/{10}^{-9}$
0	54.8	$-2.2\pm 5.8$	58.2	$-1.3\pm 4.7$
2/3	83.6	$-1.6\pm 4.5$	86.8	$-1.2\pm 3.5$
3	376.7	$-0.1\pm 0.8$	336.6	$-0.3\pm 0.6$

Table 3. Table 3: Upper limits at 95% CL on the amplitude Ω ref \Omega_{\rm ref} of the GWB assuming uniform (left) and log-uniform (right) prior. We quote the upper limits obtained employing data up to the O4a data run (O1 − - O4a), comparing these with those obtained previously  Abbott  et al. ( 2021b ) (O1 − - O3), at individual values of α \alpha (which corresponds to imposing a delta function prior on α \alpha at the quoted value), and marginalized over the α \alpha posterior (assuming a Gaussian prior on α \alpha ). We observe an overall improvement between a factor 1.2 to 2.3. The improvement in upper limits compared to O3 is more significant for smaller spectral indices. This is because the greatest enhancement in detector sensitivity occurred at low frequencies. Furthermore, the difference between the upper limit for α = 0 \alpha=0 and the marginalized case is minimal, due to statistical fluctuations from the MCMC.

	Uniform prior			Log-uniform prior
$\alpha$	O1-O4a	O1-O3	Improvement	O1-O4a	O1-O3	Improvement
0	$8.6\times {10}^{-9}$	$1.7\times {10}^{-8}$	2.0	$2.8\times {10}^{-9}$	$5.8\times {10}^{-9}$	2.1
2/3	$6.3\times {10}^{-9}$	$1.2\times {10}^{-8}$	1.9	$2.0\times {10}^{-9}$	$3.4\times {10}^{-9}$	1.7
3	$1.0\times {10}^{-9}$	$1.3\times {10}^{-9}$	1.3	$3.2\times {10}^{-10}$	$3.9\times {10}^{-10}$	1.2
Marginalized	$1.5\times {10}^{-8}$	$2.7\times {10}^{-8}$	1.8	$2.9\times {10}^{-9}$	$6.6\times {10}^{-9}$	2.3

Table 4. Table 4: Upper limits (ULs) at 95% confidence level (CL) on the amplitude of each individual polarization assuming a mixed-polarized background. The are improved by over a factor of 2 compared to ULs obtained from O1 − - O3 data  Abbott  et al. ( 2021b ) . We observe the UL for the scalar-polarized case is less stringent as the ORF for scalar polarization is smaller in amplitude than those for other polarizations  Callister  et al. ( 2017 ) .

Polarization	O1$-$O3	O1$-$O4a	Improvement
Tensor	$6.4\cdot {10}^{-9}$	$2.6\cdot {10}^{-9}$	$\sim 2.5$
Vector	$7.9\cdot {10}^{-9}$	$2.9\cdot {10}^{-9}$	$\sim 2.7$
Scalar	$2.1\cdot {10}^{-8}$	$7.4\cdot {10}^{-9}$	$\sim 2.8$

Equations

$${\Omega }_{\mathrm{GW}}(f)=\frac{1}{{\rho }_{\mathrm{c}}}\frac{\text{d}{\rho }_{\mathrm{GW}}(f)}{\text{d}\ln f},$$

$${\hat{\Omega }}_0(f)=\frac{\text{Re}[C_{IJ}(f)]}{{\gamma }_{IJ}(f)S_0(f)},$$

$$C_{IJ}(f)=\frac{2}{T}{\tilde{s}}_I^{\ast }(f){\tilde{s}}_J(f).$$

$$S_0(f)=\frac{3H_0^2{F}_0^2(f)}{10{\pi }^2}\frac{1}{f^3}.$$

$${\Omega }_{\mathrm{GW}}(f)={\Omega }_{\alpha }{F}_{\alpha }(f),F_{\alpha }(f)=\text{quantity}(\frac{f}{25\mathrm{Hz}}{)}^{\alpha }.$$

$${\sigma }_0^2(f)=\frac{1}{2T\Delta f}\frac{P_I(f)P_J(f)}{{\gamma }_{IJ}^2(f)S_0^2(f)},$$

$${\hat{\Omega }}_0(f)$$

$${\sigma }_0(f)$$

$${\hat{\Omega }}_{{\alpha }_1}(f)={\hat{\Omega }}_{{\alpha }_2}(f)\frac{F_{{\alpha }_1}(f)}{F_{{\alpha }_2}(f)}.$$

$${\hat{\Omega }}_{\alpha }{|}_{25\mathrm{Hz}}$$

$${\sigma }_{\alpha }{|}_{25\mathrm{Hz}}$$

$$p({\hat{\Omega }}_0|\Theta )\propto \exp \left[-\frac{1}{2}\sum _f\frac{({\hat{\Omega }}_0(f)-{\Omega }_{\mathrm{GW}}(f;\Theta ){)}^2}{{\sigma }_0^2(f)}\right],$$

$${\Omega }_{\mathrm{GW}}(f;\Theta )=\sum _{(\mathrm{p})}{\beta }^{(\mathrm{p})}(f){\Omega }_{\mathrm{ref}}^{(\mathrm{p})}{F}_{{\alpha }_{\mathrm{p}}}(f),$$

$${\hat{\Omega }}_{\text{mag, }IJ}(f)=\frac{2}{T}\frac{\left|T_I(f)\right|\left|T_J(f)\right|\sqrt{{\sum }_{ab}|{\tilde{m}}_{I,a}^{\star }(f){\tilde{m}}_{J,b}(f){|}^2}}{{\gamma }_{IJ}(f)S_0(f)},$$

$$p({\hat{\Omega }}_0|\Theta ,\Lambda )\propto \exp \left[-\frac{1}{2}\sum _f\frac{({\hat{\Omega }}_0(f)-\Lambda {\Omega }_{\mathrm{GW}}(f;\Theta ){)}^2}{{\sigma }_0^2(f)}\right],$$

$$p(\Lambda )=\frac{1}{\sqrt{2\pi {\sigma }_{\Lambda }^2}}\exp \left[-\frac{(\Lambda -1{)}^2}{2{\sigma }_{\Lambda }^2}\right].$$

$$\frac{{\log }_{10}{\Omega }_{\mathrm{ref},\max }-{\log }_{10}{\Omega }_{\mathrm{ref},\min }}{2}=\mathrm{3.5.}$$

$$\Theta =\{{\Omega }_{\mathrm{ref},\mathrm{T}},{\Omega }_{\mathrm{ref},\mathrm{V}},{\Omega }_{\mathrm{ref},\mathrm{S}},{\alpha }_{\mathrm{T}},{\alpha }_{\mathrm{V}},{\alpha }_{\mathrm{S}}\}.$$

$${\Omega }_{\mathrm{GW}}(f|\Theta )={\Omega }_{\mathrm{ref}}{F}_{\alpha }(f)+{\Omega }_{\mathrm{MAG}}(f|{\Theta }_{\mathrm{MAG}}).$$

$$|T_{I,\mathrm{ITS}}(f)|={\kappa }_I{\left(\frac{f}{10\mathrm{Hz}}\right)}^{-{\beta }_I},$$

$${\Omega }_{\mathrm{GW}}(f;{\theta }_k)=\frac{f}{{\rho }_c{H}_0}{\int }_0^{z_{\max }}dz\frac{R_m(z;{\theta }_k)⟨\frac{d{E}_{\mathrm{GW}}}{d{f}_s}(f_s;{\theta }_k)⟩}{(1+z)E({\Omega }_{\mathrm{m},0},{\Omega }_{\Lambda ,0},z)},$$

$$R_{\mathrm{BBH}}(z)=C({\alpha }_z,{\beta }_z,z_{\mathrm{p}})\frac{R_0(1+z{)}^{{\alpha }_z}}{1+{\left(\frac{1+z}{1+z_{\mathrm{p}}}\right)}^{{\alpha }_z+{\beta }_z}},$$

$$p({\hat{\Omega }}_0,\{d_i\}|{\Lambda }_{\mathrm{BBH}})=p_{\mathrm{BBH}}(\{d_i\}|{\Lambda }_{\mathrm{BBH}})\times p_{\mathrm{stoch}}({\hat{\Omega }}_0|{\Lambda }_{\mathrm{BBH}}).$$

Full text

I Introduction

The superposition of individually unresolved gravitational-wave (GW) signals accumulated throughout cosmic history gives rise to a gravitational-wave background (GWB) Romano and Cornish (2017); Cornish and Romano (2015) which produces a persistent incoherent signal in GW detectors. The GWB is commonly referred to as “stochastic” as it can be characterized statistically as a mean zero stochastic field Romano and Cornish (2017).

GWB sources can be categorized as cosmological or astrophysical based on their origin. The former provides a unique probe of the early Universe, as GWs propagate directly from the inflationary epoch, while the latter results from an ensemble of unresolved astrophysical events, offering insights into source populations at high redshift. Potential contributors to the GWB include distant unresolved compact binary mergers Marassi et al. (2011); Zhu et al. (2011a); Rosado (2011a); Wu et al. (2012); Marassi et al. (2011); Zhu et al. (2012), core-collapse supernovae Sandick et al. (2006); Buonanno et al. (2005); Marassi et al. (2011); Zhu et al. (2012), rotating neutron stars Zhu et al. (2010); Buonanno et al. (2005); Sandick et al. (2006); Ferrari et al. (1999); Lasky et al. (2013); Rosado (2012); Zhu et al. (2011b), stellar core collapse Crocker et al. (2017, 2015) and superradiance of axion clouds around black holes Brito et al. (2015); Takahashi et al. (2024); Sasaki et al. (2016); Arvanitaki and Dubovsky (2011); Brito et al. (2017). Some of the cosmological sources are cosmic strings Kibble (1976a); Hindmarsh and Kibble (1995); Vachaspati et al. (2015); Schmitz and Schröder (2024); Abbott et al. (2018a, 2021a); Auclair et al. (2020); Damour and Vilenkin (2000, 2005a); Siemens et al. (2007a); Olmez et al. (2010); Vachaspati and Vilenkin (1985); Blanco-Pillado and Olum (2017); Kawasaki et al. (2010), first-order phase transitions in the early Universe Witten (1984); Hogan (1986); Mazumdar and White (2019); Hindmarsh et al. (2021); Turner and Wilczek (1990); Kosowsky et al. (1992); Kamionkowski et al. (1994); Romero et al. (2021); Badger et al. (2023), primordial black holes Carr (1975); Carr et al. (2021); Matarrese et al. (1998); Baumann et al. (2007); Kohri and Terada (2018); Romero-Rodríguez et al. (2022); Inui et al. (2024); Boybeyi et al. (2025), a stiff equation of state in the early Universe Giovannini (1998); Peebles and Vilenkin (1999); Giovannini (1999); Figueroa and Tanin (2019a, b); Duval et al. (2024), and gravitational waves generated during inflation Barnaby et al. (2012, 2011); Maleknejad (2016); Thorne et al. (2018); Bartolo et al. (2016); Starobinsky (1979); Badger et al. (2024); Mukhanov and Chibisov (1981); Grishchuk (1974); Rubakov et al. (1982); Abbott and Wise (1984).

The Laser Interferometer Gravitational-wave Observatory (LIGO) Collaboration et al. (2015), Virgo Acernese et al. (2014), and KAGRA Akutsu et al. (2021) collaborations (LVK collaborations Abbott et al. (2018b)) have previously placed upper limits on both isotropic Abbott et al. (2021b) and anisotropic Abbott et al. (2021c) GWBs using data from the first three observing runs collected by the LIGO Hanford, LIGO Livingston, and Virgo interferometers. These searches relied on cross-correlating data between detector pairs Romano and Cornish (2017). The implications of the upper limits placed on the GWB energy density spectrum with LVK data have been far-reaching, setting competitive bounds on early Universe physics Abbott et al. (2018a); Romero et al. (2021); Romero-Rodríguez et al. (2022); Inui et al. (2024); Abbott et al. (2021a); Duval et al. (2024); Badger et al. (2024, 2023), and expanding our knowledge of the astrophysical binary population properties complementary to what is obtained from individual GW sources Callister et al. (2020); Abbott et al. (2021b); Turbang et al. (2024); Lalleman et al. (2025).

The most compelling astrophysical sources of a GWB in our detector range are compact binary coalescences (CBCs). The CBC background is expected to be an approximately isotropic Cusin et al. (2018); Jenkins et al. (2018); Capurri et al. (2021) spectral power-law signal, peaking around a few hundred Hz (depending on the features of the CBC population Regimbau (2011a); Capurri et al. (2021); Renzini and Golomb (2024); Ebersold and Regimbau (2025)). The signal is not expected to be Gaussian due to the intermittent nature of the background, in the frequency range probed Drasco and Flanagan (2003). In this work, we update our models to estimate the expected contributions from binary black holes (BBHs), binary neutron stars (BNSs) and neutron star-black hole (NSBH) coalescences, incorporating the most recent observational data from O4a Abac et al. (2025a, b); LIGO Scientific, VIRGO, KAGRA Collaborations (2025a, b). For the BBH population in particular, we adopt the methodology of Callister et al. (2020) to constrain the redshift-dependent merger rate by combining individually resolved events with upper limits from the cross-correlation analysis presented in this paper. This approach allows us to construct a full posterior distribution for the resulting energy density spectrum, accounting for current uncertainties in the mass, spin and redshift distributions of the BBH population. Additionally, by combining the contributions from all classes of CBCs and comparing the resulting GWB spectrum with the projected sensitivity curves of the upgraded Advanced LIGO detectors (A+) Barsotti et al. and Virgo detector (V+) Rocchi et al. (2025), we demonstrate that one year of time-coincident data between detectors in the network at design sensitivity would enable us to substantially probe and potentially detect the expected signal arising from the CBC population.

In this work we follow the same methodology as previous observing runs, cross-correlating data from LIGO Hanford and Livingston detectors during the first portion of the fourth observing run, hereinafter referred to as O4a. During this first part of the run, Virgo was not in science mode, hence why it is not included in our analysis. In this paper, we focus exclusively on the isotropic GWB, while complementary results from an anisotropic search are presented in a companion paper LIGO Scientific, VIRGO, KAGRA Collaborations (prepa). The cross-correlation search method employed is optimal for a continuous Gaussian signal, given appropriate caveats Romano and Cornish (2017); Matas and Romano (2021), while it is sub-optimal but still sensitive to an intermittent signal Lawrence et al. (2023). Several alternative search methods are currently being explored tailored to intermittent signals Smith and Thrane (2018); Buscicchio et al. (2023); Lawrence et al. (2023); Sah and Mukherjee (2023).

As we do not detect a GWB, we set upper limits on its amplitude. We include inference on a set of possible GWB models, including an astrophysical signal sourced by CBCs and alternative GWB polarization modes Nishizawa et al. (2009) arising in theories of gravity beyond General Relativity (GR). As in O3, loud glitches are excluded from the analysis K. Riles and J. Zweizig (2021); A. Matas, I. Dvorkin, T. Regimbau, and A. Romero (2021). Additionally, as correlated magnetic noise is a potential limiting factor in stochastic background searches especially as detector sensitivity improves, we construct a magnetic noise budget, and find it lies well below the current GWB sensitivity. Alongside this paper, the LVK collaboration is also releasing a companion paper exploring the cosmological implications of our results LIGO Scientific, VIRGO, KAGRA Collaborations (prepb), focusing on the constraints these place on potential early-Universe sources of the GWB.

This paper is organized as follows. In Section II, we outline the analysis methods employed throughout. Section III details the procedures implemented to remove identifiable sources of detector artifacts that can bias the results and implications of the search. Section IV presents the results of our analysis. Finally, in Section V, we discuss the astrophysical implications of our findings and provide constraints on the BBH merger rate.

II Methods

The isotropic stochastic analysis presented here targets the spectral properties of the GWB, modeling the signal as isotropic, unpolarized, stationary, and Gaussian in the limit of long observing time Romano and Cornish (2017); Renzini et al. (2022). The signal is also assumed to be weak with respect to the detector sensitivity Matas and Romano (2021).

The GWB spectrum may be expressed in terms of the fractional GW energy density spectrum $\Omega_{\rm GW}(f)$, defined as

[TABLE]

where $\text{d}\rho_{\rm GW}$ is the GW energy density in the frequency band $[f,f+\text{d}f]$, and $\rho_{\rm c}=3H_{0}^{2}c^{2}/(8\pi G)$ is the critical energy density in the Universe, with $c$ the speed of light and $G$ Newton’s constant. Integrating $\Omega_{\rm GW}(f)$ over $\text{d}\ln f$ returns the total dimensionless GW energy density. In what follows, we adopt the Planck 2015 cosmological parameters, fixing the Hubble constant to $H_{0}=67.9\,\mathrm{km}\,\mathrm{s}^{-1}\,\mathrm{Mpc}^{-1}$ and the present-day matter density parameter to $\Omega_{\mathrm{m,0}}=0.3065$ Ade et al. (2016).

The analysis presented is carried out with the pygwb analysis package Renzini et al. (2023, 2024). In the following we summarize the methodological approach and specific choices made in the analysis of the O4a dataset. Methods employed here do not vary significantly from those of previous observing runs O1 Abbott et al. (2017), O2 Abbott et al. (2019), and O3 Abbott et al. (2021b), so that the resulting spectra may be ultimately combined into a single result. Changes with respect to previous analyses are highlighted in the text.

II.1 Cross-correlation spectra

Following Allen and Romano (1999), we employ cross-correlation spectra calculated from pairs of independent, time-coincident detector strain datastreams to construct an estimator for $\Omega_{\rm GW}(f)$:

[TABLE]

where $C_{IJ}$ is the cross-spectral density (CSD) calculated from the Fourier-transformed strain data $\tilde{s}(f)$ of duration $T$ from two non-coincident detectors $I$ and $J$:

[TABLE]

The geometric quantity $\gamma_{IJ}(f)$ is referred to as the overlap reduction function (ORF) Finn et al. (2009); Romano and Cornish (2017), which quantifies the frequency-dependent decoherence that affects the cross-correlation of data from two geographically separated interferometers. This reduction depends on the individual detector responses as well as their relative separation and orientation Romano and Cornish (2017). The function $S_{0}(f)$ converts units of GW strain power into fractional energy density, whilst encoding also a frequency weighting specific to a chosen spectral model $F_{0}(f)$:

[TABLE]

As discussed in Matas and Romano (2021), in the case of noise-dominated data streams the cross-correlation statistic described above is a quasi-optimal 111Sub-optimality is due to corrections to the variance Eq. (6) of the same order of magnitude of the signal, which become negligible when the signal is weak compared to the detector noise. estimator for the background spectrum.

In line with previous analyses Abbott et al. (2017, 2019, 2021b), here we assume a power-law spectral model expressed as an amplitude $\Omega_{\alpha}$ with respect to the fixed reference frequency $f_{\rm ref}=25\,{\rm Hz}$ scaling with a power-law index $\alpha$,

[TABLE]

In Eq. (2) specifically, the chosen spectral model is fixed to a power law with spectral index $\alpha=0$, as indicated by the subscript of $S_{0}$ and that of the estimator in Eq. (2). The corresponding variance estimate is Romano and Cornish (2017)

[TABLE]

where $P_{I}(f)$ is the one-sided power spectral density for detector $I$, and $\Delta f$ is the frequency resolution employed in the analysis.

In practice, due to noise non-stationarities and computational efficiency the ${\hat{\Omega}}_{0}(f)$ spectrum is estimated from a weighted average over a large number of data segments, where for each time segment $i$ there is a corresponding CSD estimate $C_{IJ,\,i}$ and PSD estimates $P_{I,\,i}$, $P_{J,\,i}$. Further details about this cross-correlation computation may be found in the pygwb publication Renzini et al. (2023).

II.2 Frequency estimator

We construct a narrow-band estimator ${\hat{\Omega}}_{0}(f)$ with associated variance $\sigma_{0}^{2}(f)$:

[TABLE]

where $i$ labels data from $t_{i}$ to $t_{i}+\Delta t$. In practice, the cross-spectral densities are calculated from the data with a 50% overlap to make up for the impact of Hann-windowing each segment, hence each independent time estimate is achieved by handling even and odd segments separately as described in Renzini et al. (2023).

This calculation is performed employing the $S_{0}(f)$ weighting function. However, note that it can be performed by choosing any other power-law weighting $\alpha$, and it is possible to re-weight from one index, $\alpha_{1}$, to another, $\alpha_{2}$, as

[TABLE]

Broadband estimates for a spectral index $\alpha$ and at reference frequency $25\,{\rm Hz}$ are obtained by optimally weighting the narrowband estimator, as described in Abbott et al. (2021b); Renzini et al. (2023):

[TABLE]

II.3 Parameter estimation

Parameter estimation is performed here using a hybrid approach that combines frequentist and Bayesian analysis techniques Mandic et al. (2012). Specifically, the frequentist estimators from the previous sections, ${\hat{\Omega}}_{0}(f)$ and $\sigma_{0}(f)$, combined over the entire observing period, are used as input data to calculate posterior probability distributions of GWB model parameters. These in turn are used to compute upper bounds on the spectral amplitude of the GWB. In Matas and Romano (2021) the authors show that in the weak-signal approximation this hybrid approach does not lose information compared to a fully Bayesian search, which produces posterior distributions from the full time series data rather than from the time-averaged frequency estimators.

We assume that, in each frequency bin, the un-weighted spectral estimate ${\hat{\Omega}}_{0}$ is Gaussian-distributed with variance $\sigma^{2}_{0}$ Abbott et al. (2021b):

[TABLE]

where $\Omega_{\rm GW}(f;\Theta)$ is an assumed model for the GWB signal, parametrized via a set of model parameters $\Theta$. Posteriors for the model parameters $\Theta$ are obtained by employing the pygwb library Renzini et al. (2023) and the Bilby library Ashton et al. (2019) for detector modeling and likelihood evaluation. As in previous work, e.g. Abbott et al. (2021b), this analysis is performed using data from multiple detector pairs, $I,J$ (i.e., baselines). Data from different baselines are assumed to be independent, such that the combined inference assumes likelihood (12) holds for data from each individual baseline, ${\hat{\Omega}}_{0}^{IJ}(f),\,\sigma_{0}^{IJ}(f)$, and the total likelihood is the result of the product of the likelihood calculated for each baseline.

We consider different models for the GWB:

• •

Noise: $\Omega_{\rm GW}(f;\Theta)=0$. We implicitly include uncorrelated Gaussian noise as part of every model that follows in the prior choices for the model parameters. This noise model is not strictly correct as there are potential sources of correlated noise, such as Schumann resonances Janssens et al. (2023a). However, these effects may be neglected in the present search due to detector sensitivity, as shown in Sec. III.6.

• •

Power law: $\Omega_{\rm GW}(f;\Theta)=\Omega_{\rm ref}F_{\alpha}(f)$, as per Eq. (5). In this analysis, the two parameters $\Theta=\{\Omega_{\rm ref},\alpha\}$ are inferred independently, and hence we choose this notation to avoid confusion with the cross-correlation estimators described above. The reference frequency is set to $f_{\rm ref}=25\,{\rm Hz}$ as in past analyses Abbott et al. (2021b) due to the location of peak stochastic sensitivity, for lower values of $\alpha$. Note that the precise sensitivity depends on the spectral index $\alpha$.

• •

Scalar-vector-tensor (SVT) power law: this model assumes that the signal may include vector and scalar polarizations as well as the standard tensor polarization predicted by General Relativity. This background is modeled as a superposition of individual contributions, labeled $\rm(p)=\{T,V,S\}$ where each component is individually modeled by a power law Callister et al. (2017):

[TABLE]

where $\beta^{\rm(p)}(f)\equiv\gamma_{IJ}^{\rm(p)}(f)/\gamma_{IJ}(f)$ is the ratio between the ORF for polarization (p) and the standard (tensor) polarization Nishizawa et al. (2009).

• •

Magnetic: describes correlations between pairs of detectors induced by large-scale coherent magnetic fields that can mimic a GWB signal. These are modeled in terms of magnetometer correlations and a transfer function between the local magnetic field and the strain channel of the detectors. Further discussion on this model is given in section IV.4.

• •

Compact binary coalescences (CBC): the predominant contribution to the GWB is expected to be the superposition of all merging compact binaries in the observable Universe Regimbau (2011a). The total GW fractional energy density due to compact binary mergers, $\Omega_{\rm CBC}$, can be parametrized as a function of the mass distribution of compact binaries and their merger rate as a function of redshift. This model will be discussed in Section V. Assuming $\Omega_{\rm CBC}$ is the only contributor to the stochastic signal in the analyzed frequency band, we can infer mass and merger rate population parameters.

To test these different models against the data, we calculate Bayes factors for each model. Unless otherwise specified, we compute the Bayes factor between the hypothesis of a specific signal model with respect to the noise model, $\Omega_{\rm GW}(f;\Theta)\equiv 0$.

III Data Quality

This section describes the data quality and handling specifics for the strain data collected during the first portion of the fourth observing run (O4a) by the LIGO Hanford and Livingston detectors. Note that results presented in the following sections include the combination of this data with previous datasets O1, O2, O3, whose relevant instrument performances and data quality features are described in Abbott et al. (2016a); Martynov et al. (2017); Nuttall et al. (2015); Abbott et al. (2018c), Davis et al. (2021); Brooks et al. (2021); Driggers et al. (2019), Buikema et al. (2020); Tse et al. (2019); Nguyen et al. (2021); Soni et al. (2020) respectively. The O3 analysis also includes data from Virgo Fiori et al. (2020); Acernese et al. (2023a, 2019, 2020, 2022, b). For further details on the performance of the LIGO detectors and the quality of their data during O4a Capote et al. (2025); Jia et al. (2024), see also Soni et al. (2025).

III.1 Data

The O4a observing period ran from 15:00 UTC on May 24, 2023 to 16:00 UTC on January 16, 2024. The total coincident time for the Hanford-Livingston (HL) baseline in this period was 126.57 days, before applying data quality cuts. The strain data used in this analysis were calibrated strain data with narrow-band features from photon calibrator and actuation injections and power mains noise, subtracted Viets and Wade (2021), followed by the application of non-stationary noise subtraction Vajente et al. (2020). The correlated magnetic noise was monitored using low-noise LEMI120 magnetometers located at each site lem . The same data processing was carried out for both strain and magnetometer data to allow a comparison between the magnetic and GW analyses.

Both strain and magnetic datastreams are initially downsampled from the original sample rate of 16384 Hz to 4096 Hz. This downsampling is used because the GW sources targeted in this search populate frequencies below $2$ kHz. The data are high-pass filtered using a 16th-order Butterworth filter with a knee frequency of 11 Hz, which is constructed using second-order sections. As in analyses described in Abbott et al. (2021b); Renzini et al. (2023), each data stream is divided into 50$\%$ overlapping, Hann-windowed 192-second segments, coarse-grained to a frequency resolution of 1/32 Hz. The analyzed frequency range spans $20-1726$ Hz; the minimum frequency is chosen based on the detector sensitivities Soni et al. (2025), while the maximum frequency is chosen to be sufficiently below the Nyquist frequency to avoid aliasing effects.

In the following text, tables, and figures, we label the LIGO Hanford detector as LHO and the LIGO Livingston as LLO.

III.2 Time domain cuts: vetoed times

Data quality for analysis purposes was assessed jointly with the Detector Characterization groups of LIGO, Virgo, and KAGRA. Data quality flags are employed to define observing time periods to be vetoed from analysis. Most significant to this search are the “Category 1” vetoes. These vetoes correspond to times when the detector was known to be operating outside of its nominal condition. For more details on the O4a data quality and the specific time-domain cuts made for stochastic data quality purposes, please refer to Soni et al. (2025), Sec. 4.3. After applying vetoes, the remaining viable time-coincident data is 126.47 days. The relevant segment information used to exclude vetoed data from this analysis will be available on GWOSC LIG (2025).

III.3 Time domain cuts: gating

In the third LVK observing run (O3), a large population of loud glitches Abbott et al. (2021b) led to very frequent flagging of data segments by non-stationarity checks (see Sec. III.4). To mitigate this effect, time-domain gating of the data was introduced Riles et al. (2020); Matas et al. (2021): loud glitches are removed from the data stream by multiplying the data with an inverse Planck-taper window with parameters tailored to the specific non-stationary feature (e.g., duration) McKechan et al. (2010).

In O4a, data are preprocessed with an auto-gating procedure as described in Renzini et al. (2024). The parameters chosen for the auto-gating were optimized for this specific search, minimizing the amount of data lost whilst retaining stationarity. A total of $0.093\%$ of analyzed Hanford data and $0.068\%$ of Livingston data were gated; more details can be found in Sec. 4.3.2. of Soni et al. (2025).

III.4 Time domain cuts: stationarity checks

As in previous runs, we apply a non-stationarity cut (often referred to as the $\Delta\sigma$ cut Renzini et al. (2023)) by removing times where the square root of the variance in Eq. (6) is found to vary by more than $20\%$ between neighbouring segments. We cut the union of segments that do not satisfy the stationarity test assuming a set of different spectral weights: $\alpha=\{-5,0,3,5\}$, implemented as in Eq. (2). As each power law accumulates SNR in a different frequency band, these may flag different segments depending on their specific spectral narrowband/broadband features. This cut additionally removes 8.1$\%$ of the analysed segments, leaving 108.41 days of viable data.

III.5 Frequency domain cuts: narrowband features

In addition to broadband features that determine the exclusion of entire data segments from the analysis, also narrowband deviations from Gaussianity are investigated, and data cuts are applied to exclude problematic frequency bins. These frequencies are determined by monitoring the coherence between the GW channels at different sites, and between the GW channel and auxiliary channels at the same site. Bins that show evidence of instrumental noise or that manifest strong non-Gaussian behavior are excluded from the analysis: for more details, see Soni et al. (2025), Sec. 2.1. The list of problematic frequencies, i.e., the notch list, is determined at the end of the observing run and then applied identically throughout the entire data analysis process (with a small exception in parameter estimation, as discussed in Sec. IV.2). This list is included in the appended data release.

The notches include the following components:

• •

Calibration lines: these are stationary lines injected into each detector’s differential arm feedback loop so as to continually monitor the detector response and the corresponding calibration amplitude factor. These include the stationary calibration lines injected throughout O4a Karki et al. (2016), and one calibration line at $24.5$ Hz that was turned on between July 25, 2023, and August 9, 2023.

• •

Pulsar injections: each injection corresponds to a simulation of a quasi-monochromatic, persistent GW signal as expected from a rapidly rotating, isolated neutron star. These were injected throughout the entire O4a run at both LIGO sites and emulate how canonical continuous-wave (CW) signals would appear in the detectors Biwer et al. (2017).

• •

Quadruple suspension violin modes: resonances of the detector mirror suspension fibers Collaboration et al. (2015).

• •

Powerline harmonics: integer harmonics of the 60 Hz power mains present at both detectors.

• •

Other artifacts: DuoTone signals (960Hz, 961Hz) and their 1 Hz side band present in both LIGO detectors in O4a from the timing system, used to synchronize data collection across the global detector network and in each interferometer Sullivan et al. (2023).

Overall, the notch list cuts 8.4% of the spectrum; however, note that in the frequency range $[20,300]$ Hz where the detectors are most sensitive, only 2.7% of the spectrum is cut.

III.6 Magnetic noise budget

Correlated magnetic noise present across multiple detectors is a probable hindrance for stochastic analyses Thrane et al. (2013, 2014); Coughlin et al. (2018). Examples of such correlated magnetic noise are Schumann resonances Schumann (1952a, b), standing waves in the cavity between the Earth’s surface and the ionosphere that may be excited by lightning strikes; and at higher frequencies, $f>100$ Hz, correlated fields from lightning strikes. Magnetic fields of any source couple to the detector via electric cables in the detector buildings or couple to the actuation magnets controlling the detector test mass mirrors Fiori et al. (2020).

Magnetic fields are monitored with sensitive probes placed outside the detector buildings to produce a magnetic “noise budget”, quantifying the strength of correlated magnetic noise present across detectors, which is then compared to the sensitivity of the detector network to a stochastic signal Janssens et al. (2023b). If this budget becomes large enough, correlated magnetic fields can contribute to correlated noise which may mimic a signal in the stochastic analysis presented here.

We compute the magnetic correlated noise budget via

[TABLE]

where $T_{I}(f)$ and $T_{J}(f)$ are the magnetic coupling functions of detectors $I$ and $J$, and $\tilde{m}_{I,a}(f)$ is the Fourier-transformed magnetic data from the outside probes of detector $I$ pointing in the direction $a$, which can be either along the X or Y arms of a detector, and its sum runs over the four possible combinations of $a$ and $b$. Magnetic coupling functions consist of the product of two terms: an outside-to-inside coupling, $T_{\rm OTI}$, and an inside-to-strain coupling, $T_{\rm ITS}$. The $T_{\rm OTI}$ term is given by the ratio of the strength of a magnetic source measured outside the detector buildings to the strength of the same source measured inside the buildings, and reduces to an overall constant averaged over several independent measurements, while $T_{\rm ITS}$ is estimated via hardware injections as explained below. Uncertainties on these measurements include the uncertainty on the outside-to-inside coupling as measured in O3 at LLO, and, to be conservative, an intrinsic uncertainty factor of 2 as the magnetometers in the building (used for the measurements) are not located at the exact location where the magnetic field couples to the strain.

The $T_{\rm ITS}$ functions are estimated by each detector site team by injecting far-field magnetic fields within the detector building and measuring the response in the strain channel Nguyen et al. (2021); Janssens et al. (2023a); Fiori et al. (2020). A value is recorded as a measurement when both the magnetometer sensor and strain channel amplitude spectral densities (ASDs) exceed their respective set thresholds; otherwise, if only the magnetometer ASD exceeds the threshold, the value is treated as an upper limit. The coupling function is interpolated to the frequencies of the computed magnetic CSD, $\operatorname{Re}\left[\tilde{m}_{I}^{\star}(f)\tilde{m}_{J}(f)\right]$. This does not create any artifacts in the coupling functions as these have a finer frequency resolution than the CSD. We include a calibration factor in the magnetic CSDs222The calibration factor is an intrinsic conversion factor between the measured current in the magnetometers and the actual magnetic field. and take the quadrature sum of the magnetic CSDs of four baselines, one for every combination of the two outside magnetometers at each detector used for computing the budget, as a most conservative estimate. These quantities are then combined in Eq. (14) to produce the magnetic budget shown in Fig. 1. Here the red bands show the uncertainties on the budget, compared to the $2\sigma$ power-law integrated sensitivity (PI) curve Thrane and Romano (2013), showing the sensitivity of our search to power-law backgrounds.

The final budget indicates that no correlated (magnetic) noise observation is expected. Some narrowband features are expected, such as the harmonics of the US power lines at 60 Hz, and the ORF zero-crossings (related to the fact that a single baseline was considered for this analysis).

While the magnetic budget computation in O3 was performed differently333In O3, weekly magnetic couplings were used and, due to fluctuations, computing a budget every single week produces a wide uncertainty band instead of a single line Nguyen et al. (2021)., we can compare the two budgets at an order-of-magnitude level. At lower frequencies $f<60$ Hz, the O4a budget is consistently within the band of the O3 budget, but in a narrow frequency band $60\,\mathrm{Hz}<f<65$ Hz, the O4a budget becomes slightly larger than the maximum value in O3. This is related to a zero-crossing of the ORF virtually enhancing the magnetic budget due to its presence in the denominator of Eq. (14), and does not accurately represent a magnetic response of the detector and thus the actual magnetic budget in O4a. At higher frequencies $f>65$ Hz, the O4a budget is broadly consistent with the O3 magnetic budget.

III.7 Calibration uncertainties

Calibration uncertainties account for both statistical uncertainties and systematic errors arising from the calibration process of LVK data444Calibration is the procedure used to convert the digital output of a gravitational-wave detector into the relative displacement of the test masses within the detector. Sun et al. (2020, 2021). These are incorporated in parameter estimation Vitale et al. (2021); Renzini et al. (2023) by introducing a Gaussian-distributed, positive-definite calibration factor $\Lambda$. As a result, the likelihood presented in Eq. (12) is modified as (see also App. B of Renzini et al. (2023))

[TABLE]

where $\Lambda$ is marginalized over, assuming the prior

[TABLE]

The calibration variance $\sigma_{\Lambda}$ is estimated in collaboration with the LIGO Calibration Group. Hourly calibration uncertainties provided by the Calibration Group are combined to obtain uncertainties in both the magnitude and phase of the strain at the 68% confidence level (CL). These uncertainties are frequency-dependent; however, we adopt a conservative approach and select the largest uncertainty across the entire frequency range for both magnitude and phase. The corresponding maximum uncertainties are summarized in Table 1.

Provided it remains approximately below 5 degrees, the phase uncertainty does not introduce a significant bias in stochastic background estimates Yousuf et al. (2023). Hence we restrict the calculation of calibration uncertainty to the magnitude term only, taking the quadrature sum of the magnitude uncertainties for each detector strain. This yields a final value of $\sigma_{\Lambda}=8.1\%$.

IV Results

We present results obtained with the analysis methods described above using data from the O4a observing run combined with data from the O1$-$O3 observing runs. We do not find evidence of a background signal in our analyses, as outlined below, hence we set upper limits on several models. We perform three different analyses: (i) we calculate the point estimate spectrum and variance of the GWB signal using cross-correlation optimal filtering (Sec. IV.1) and set upper limits (Sec. IV.2) considering a power-law spectral model for the overall background, both in the case of a fixed known spectral index $\alpha$, and marginalizing over $\alpha$; (ii) we set upper limits on SVT model parameters (Sec. IV.3); and (iii) we perform joint inference in the case where a correlated magnetic signal is present in the data together with a GWB (Sec. IV.4).

IV.1 Cross-correlation estimates

We produce broadband estimates ${\hat{\Omega}}_{\alpha}$ and associated variance using an optimal filter for three different power-law models motivated by different potential GWB sources:

• •

$\alpha=0$ corresponds to a scale-invariant signal, which approximately describes a GWB arising from a cosmic string network Kibble (1976b); Sarangi and Tye (2002); Damour and Vilenkin (2005b); Siemens et al. (2007b); Abbott et al. (2018a) or slow-roll inflation Starobinsky (1979); Turner (1997, 1997); Bar-Kana (1994) in the LVK frequency band.

• •

$\alpha=2/3$ characterizes the CBC GWB when the inspiral phase dominates, which is a good approximation in the LVK frequency band Regimbau (2011b). However, this approximation may break down for mergers of binaries originating from Population III stars Périgois et al. (2021) or for heavy BBH mergers with masses exceeding the pair-instability mass gap Ezquiaga and Holz (2021).

• •

$\alpha=3$ is a fiducial choice which provides an approximate description of a subset of astrophysical sources, such as supernovae Marassi et al. (2009); Sandick et al. (2006). It also corresponds to a GWB with a flat strain power spectral density, following $S_{h}(f)\propto f^{-3}\Omega_{\rm GW}(f)$ Allen and Romano (1999). This is the spectral shape detectors are most sensitive to by construction, as the optimal filter depends solely on the noise PSD, hence frequencies are weighted according to the detector sensitivity thereby maximizing the cross-correlation search sensitivity.

The point estimate spectrum $\hat{\Omega}_{0}(f)$ computed as defined in Eq. (7) is shown in Fig. 2. The 1$\sigma$ contours calculated from O4a data as well as the entire O1$-$O4a dataset are compatible with 0, as expected in the presence of uncorrelated noise and in the absence of a correlated gravitational-wave signal. This behavior is consistent with the expectation that, after applying the data quality cuts described in Section III, the data remain dominated by Gaussian noise. This is further demonstrated by the Kolmogorov-Smirnov (KS) test, which returns a p-value of $p=0.99$ in favour of Gaussianity.

Point estimates ${\hat{\Omega}}_{\alpha}$ and $1\sigma$ uncertainty using the O4a data only as well as combined results from all the advanced observing runs are presented in Table 2. Point estimates for all spectral models are consistent with 0 within $1\sigma$, which implies no detection of a GWB. We also report the upper edge of the frequency band that contributes 99% of the sensitivity to each spectral-weighted search, $f_{99\%}$.

We observe a decrease in variance with increasing spectral index, which is expected since, as mentioned above, our search is most sensitive to a flat strain spectrum ($\alpha=3$). Notably, the sensitivity improvement compared to O1$-$O3 results is a factor of $1.6$ for $\alpha=0$, $1.6$ for $\alpha=2/3$, and $1.4$ for $\alpha=3$. This is driven by the fact that the most significant sensitivity enhancement in our detectors occurred at lower frequencies Capote et al. (2025), below 40 Hz, as discussed below.

Altogether, these observations indicate that there is no evidence for a detectable GWB in the analyzed data. We also note that O4a contributes primarily at frequencies below 60 Hz, as highlighted in the right-hand-side panel of Fig. 2 which shows the portion of the spectrum in the $20-60\,{\rm Hz}$ band, where the $1\sigma$ contours from O4a and O1$-$O4a data are almost overlapping. The frequency dependence of the sensitivity change between O4a and O3 is described in Fig. 3, which shows the ratio of the PSDs calculated in the two runs. We observe that this ratio lies consistently above 1 in the range $10-60\,\rm Hz$, which matches the range highlighted in Fig. 3 where the O4a data contributes the most to the measurement of a stochastic spectrum. The effect of this frequency-dependent improvement in sensitivity to the GWB between O4a and past runs is further quantified below.

IV.2 Upper Limits on the GWB

Given the absence of evidence for a signal (see Table 2), we set upper limits on the amplitude of the GWB, assuming a power-law model as in Eq. (5). The narrowband estimators ${\hat{\Omega}}_{0}(f)$ and $\sigma_{0}(f)$ are employed in inference as described in Sec. II.2, where the frequency lines described in Sec. III.5 are notched from the spectra. We explore different priors for both the amplitude and the spectral index. Specifically, we consider uniform and log-uniform priors for the amplitude, and uniform and Gaussian priors on the spectral index. A log-uniform prior is more sensitive to signals of small amplitude, yielding tighter constraints, while a uniform prior provides more conservative limits. In both cases, we set the prior range for the $\Omega_{\rm ref}$ parameter between $\Omega_{\rm ref,min}=10^{-13}$ and $\Omega_{\rm ref,max}=10^{-6}$.

We perform a fixed model analysis, fixing the spectral index to $\{0,2/3,3\}$, where each value corresponds to a specific source as indicated above. We also perform joint inference on the amplitude and spectral index, where we test both a Gaussian prior and uniform prior on the latter. The uniform prior range is set to $\alpha\in[-7,7]$ 555This prior has been widened with respect to past analyses Abbott et al. (2021b), which employed $\alpha\in[-5,5]$, to investigate a small mode in the posterior in $\alpha$ which appeared at the edge of the prior range, as discussed below.. In the case of a Gaussian prior on $\alpha$, we impose zero mean and a standard deviation of

[TABLE]

The choice of the prior over $\alpha$ can be understood as follows. The log-uniform prior over $\Omega_{\rm ref}$ induces some implicit prior over $\alpha$ that can be shown Callister et al. (2017) to be a triangular prior centered on $\alpha=0$ and non-zero for $|\alpha|\leq(\log_{10}\Omega_{\rm ref,max}-\log_{10}\Omega_{\rm ref,min})$. To avoid a vanishing prior outside of this range, we choose a Gaussian prior for $\alpha$ with a standard deviation comparable with the triangular prior Callister et al. (2017).

We set spectral-independent constraints on the signal amplitude by marginalizing over the spectral index posteriors obtained with the Gaussian prior. We take the results obtained with the log-uniform prior on the amplitude as our main result.

Table 3 summarizes results for the fixed spectral index and marginalized analyses. We set 95% confidence level upper limits on $\Omega_{\rm ref}$ at $f_{\rm ref}=25\,{\rm Hz}$ at individual values of $\alpha$, $\Omega_{\rm ref}^{(\alpha)}$, fixing the $\alpha$ prior to a delta function at the chosen values; we also set upper limits marginalized over the Gaussian $\alpha$ posterior. For comparison, we report the corresponding upper limits from O1$-$O3 analyses. As expected, log-uniform priors on the background amplitude yield tighter constraints compared to uniform priors. Broadband constraints at 25 Hz are tightest for $\alpha=3$ as this corresponds to a flat strain power spectrum, yielding the most broadband detectable signal, while constraints are weakest for a scale-invariant spectrum as this signal corresponds to a strain power spectrum with the most negative spectral index out of the set. Furthermore, marginalized results over the spectral index yield overall weaker constraints as these admit higher uncertainty on the signal model. Overall, we observe improvements in the upper limits ranging from a factor of 1.2 to 2.3.

Posterior distributions for $\Omega_{\rm ref}$ and $\alpha$ obtained in these analyses are shown in Fig. 4 assuming a log-uniform prior on $\Omega_{\rm ref}$ and either a Gaussian or uniform prior on $\alpha$. As may be observed in Fig. 4, $\alpha$ is not constrained, and the prior on $\alpha$ has little impact on the $\Omega_{\rm ref}$ posteriors. The Bayes factor comparing the signal-versus-noise hypothesis to the noise-only hypothesis is $\log_{10}\mathcal{B}=-0.19\pm 0.03$ ($\log_{10}\mathcal{B}=-0.11\pm 0.03$) at 68% confidence level assuming a Gaussian (uniform) prior on $\alpha$, further confirming no evidence of a signal in the data.

We point out a small narrow mode that appears in the posterior plot at $\alpha=5$, present for both prior assumptions and more evident in the case of a uniform prior on $\alpha$. Note that $\alpha=5$ is not associated with any known physical GW signals and usually coincides with the edge of the broad prior range assumed for this analysis, however in this case the prior was broadened to confirm no railing was occurring. This mode has been traced back to weak, broad-band residual noise coherence between the LHO and LLO detectors during O4a in the $650-850\,\rm Hz$ band, which corresponds to the frequency range most sensitive to an $\alpha=5$ power-law signal. The ORF introduced in Eq. (2) for the LHO–LLO baseline is highly suppressed at these frequencies (to $\sim 0.2-0.1\%$ relative to its magnitude at 1 Hz), strongly suggesting this weak coherence is detector noise and not GW signal. Excluding this frequency range from the analysis removes this feature completely. The four 60 Hz power-line harmonics in this range, [660, 720, 780, 840] Hz, appear significantly wider and noisier than in past observing runs, hence notches around these lines have been widened accordingly for this analysis. Numerous other weak unidentified noise features such as broadened lines and bumps have been found in both detectors. While these are not loud enough to significantly impact the analysis presented here, as demonstrated by broad-band signal point estimates (Table 2) as well as the Bayes’ factors above, the significance of this feature will be revisited in the next portion of the O4 data.

Finally, Fig. 5 shows upper limits on $\Omega_{\rm ref}$ at 95% CL as a function of $\alpha$. Note that $\alpha$ values in this plot are not directly motivated by specific sources, but rather are chosen to highlight the overall trend: upper limits decrease as the spectral index increases, which is expected due to the extrapolation of the low-frequency sensitivity at higher frequencies, and as a consequence of $\Omega_{\rm GW}\propto f^{3}S_{h}(f)$. Fig. 5 also includes upper limits obtained up to and including the O3 data run. The relative improvement between these and the limits obtained including O4a decreases as $\alpha$ increases, indicating that the sensitivity enhancement between runs is more significant at lower frequencies, as discussed above.

IV.3 Non-GR polarizations

We employ the approach first presented in Callister et al. (2017) to constrain potential deviations from GR. The detection of scalar or vector polarization modes would provide direct evidence of a violation of GR, whereas their absence allows us to place stringent constraints on extended theories of gravity.

Multiple runs with differently polarized backgrounds were conducted. First, we consider a mixed-polarized background, as introduced in Sec. II.3. In this scenario, the parameters of our model are

[TABLE]

We set a log-uniform prior on the amplitudes $\Omega_{\rm ref,(p)}$, for $\rm(p)=\{T,V,S\}$, over the range $10^{-13}$ to $10^{-6}$ and a Gaussian prior on the spectral indices, centered at 0 with a standard deviation of 3.5. We set 95% CL upper limits on $\Omega_{\rm ref,(p)}$, presented in the third column of Table 4. We find an improvement in the upper limits of over a factor of 2 with respect to O3 results Abbott et al. (2021b). It is worth noting that the upper limit for the scalar-polarized case is less stringent as the ORF for scalar polarization is approximately three times smaller in amplitude than those for other polarizations Callister et al. (2017). As a result, our detector network is inherently less sensitive to a scalar-polarized background.

Additionally, we set upper limits on the amplitudes $\Omega_{\rm ref,T}(25\rm Hz)\leq 2.9\cdot 10^{-9}$, $\Omega_{\rm ref,V}(25\rm Hz)\leq 2.7\cdot 10^{-9}$ and $\Omega_{\rm ref,S}(25\rm Hz)\leq 2.9\cdot 10^{-9}$, assuming single-polarized backgrounds. The upper limits are similar across polarizations due to the reduced detector network in O4a, since only LHO and LLO were operational, unlike O3 which included all three detectors. This limits our ability to distinguish polarization modes.

The signal versus noise hypothesis Bayes factor yields $\log_{10}\mathcal{B}=-0.87\pm 0.01$ at 68% confidence level, indicating no evidence for the presence of a signal in our data. In this case, the signal is defined as the equal-weighted combination of the individual polarization components, $\{\rm T,V,S\}$, as well as their combinations: tensor–vector $(\rm TV)$, tensor–scalar $(\rm TS)$, vector–scalar $(\rm VS)$, and tensor–vector–scalar $(\rm TVS)$. This approach treats all signal configurations equally.

Additionally, the Bayes factor for non-GR versus GR-supported polarized backgrounds yields $\log_{10}\mathcal{B}=-0.44\pm 0.03$ at 68% confidence level. This result suggests that a tensor-polarized background (in accordance with GR) is favoured over vector, scalar, and mixed-polarized alternatives. However, in the absence of a detection we expect this result to be driven by the Occam penalty factor introduced by overfitting un-informed parameters in the search.

IV.4 Joint fit for GWB and magnetic noise

Correlated magnetic noise can provide a significant limitation to the search for an isotropic stochastic signal Thrane et al. (2013, 2014); Coughlin et al. (2018); Janssens et al. (2023a). To mitigate this, the likelihood presented in Eq. (12) may be extended to jointly model correlated magnetic noise appearing in the GW detectors and a GWB power-law signal, as described in Meyers et al. (2020):

[TABLE]

In Section III.6, we showed that correlated magnetic noise is not expected to contaminate O4a data. Nevertheless, this analysis provides a complementary check by simultaneously fitting both the GWB and the correlated magnetic noise. Instead of utilizing the measured coupling functions at the detectors from Section III.6, we now model the magnetic term $\Omega_{\rm MAG}$ with the parameters $\Theta_{\rm MAG}=\{\kappa_{I},\kappa_{J},\beta_{I},\beta_{J}\}$, approximating the inside-to-strain magnetic coupling functions $T_{I,\rm ITS}(f)$ as a simple power law,

[TABLE]

where $\kappa_{I}$ is the magnitude of the amplitude at $10\,\rm Hz$ and $\beta_{I}$ is the spectral index of the power law. The prior distribution for $\kappa_{I}$ is log-uniform in the range $10^{-25}$ to $10^{-22}\,\rm pT^{-1}$ for both detectors and the prior on the spectral index $\beta_{I}$ is uniform in the range $0-12$ and $1-10$ for HLO and LLO respectively, related to measured minimum and maximum values of the spectral index during magnetic injections. As before, the GWB model parameters are $\Omega_{\rm ref}$ and $\alpha$. This allows us to estimate either the evidence for a simultaneous measurement of a background and correlated magnetic noise, or estimate the evidence for correlated magnetic noise only.

We set joint constraints on the power-law GWB and the magnetic correlated noise, employing the same priors for the magnetic parameters as in the O3 analysis Abbott et al. (2021b). We find the log Bayes factor comparing a magnetic noise model to a (Gaussian) noise-only model is $\log_{10}\mathcal{B}_{N}^{\rm MAG}=-0.048$, indicating no evidence in favour of the presence of magnetic noise. We also find that for a joint model with a power-law background and magnetic noise modeled as described above that $\log_{10}\mathcal{B}_{N}^{\rm PL+MAG}=-0.306$, showing no preference for such a model compared to a noise-only model.

V Astrophysical Implications

Leveraging the upper limits on the GWB presented in Section IV.2, we examine the implications of these results for the GWB of astrophysical origin, namely from the coalescing population of compact binaries including BBHs, BNSs, and NSBHs. We compare our upper limits to the most recent predictions for the energy density spectrum associated with different classes of CBC sources. In particular for BBHs, following Ref. Callister et al. (2020), we combine BBH observations in the local Universe with stochastic upper limits to constrain their merger rate at high redshifts, and predict the uncertainty in the BBH background using a fully GW data-driven approach.

V.1 CBC Model

To model the GWB from binary systems we follow in part the approach of LVK isotropic search papers from the O1 to O3 observing runs Abbott et al. (2016b, 2017, 2019, 2021b), and additionally we introduce advances in our knowledge of the CBC population Abbott et al. (2021d, 2022); LIGO Scientific, VIRGO, KAGRA Collaborations (2025b). We categorize the population into classes, each labeled by $k$, to differentiate between BBHs, BNSs and NSBHs. Each class is characterized by distinct values of source parameters, such as masses and spins, which we denote as $\theta_{k}$. The total astrophysical background is then the sum of contributions from all classes. The contribution of class $k$ to the background can be expressed as an integral over the redshift $z$, given by Phinney (2001); Regimbau (2011a); Zhu et al. (2011c); Rosado (2011b)

[TABLE]

where $R_{m}(z;\theta_{k})$ denotes the binary merger rate in the source frame per unit comoving volume, $\langle dE_{\rm GW}/df_{s}(f_{s},\theta_{k})\rangle$ is the energy spectrum emitted by a single binary evaluated in terms of the source frequency $f_{s}=(1+z)f$ and averaged over the ensemble properties of source class $k$, and $E(\Omega_{\mathrm{m,0}},\Omega_{\Lambda,0},z)=\sqrt{\Omega_{\mathrm{m,0}}(1+z)^{3}+\Omega_{\Lambda,0}}$ accounts for the dependence of comoving volume on cosmology, assuming $\Omega_{\Lambda,0}=1-\Omega_{\mathrm{m,0}}$. We impose an upper redshift cutoff at $z_{\rm max}=10$, beyond which the compact binary merger rate in conventional formation scenarios is expected to be negligible Zhu et al. (2011d). We do not consider potential contributions to the GWB from formation channels involving Population III stars Martinovic et al. (2022) or primordial black holes Mandic et al. (2016); Mukherjee and Silk (2021).

The energy spectrum $\langle dE_{\rm GW}/df_{s}\rangle$ is determined from the strain waveform of the binary system. For computing the backgrounds of BBH and NSBH we use IMRPhenomXP, a phenomenological inspiral-merger-ringdown frequency-domain model for gravitational-wave signals emitted by quasi-circular precessing BBHs Pratten et al. (2020, 2021). For BNS background computation, we employ the IMRPhenomXP_NRTidalv3 waveform model, which tapers the merger-ringdown part from the baseline IMRPhenomXP model and has the option to include tidal effects Abac et al. (2024a); Colleoni et al. (2025); Dietrich et al. (2019). To compute the GWB from CBCs throughout the Universe according to Eq. (21), we use methods developed for the GWB calculation in Abbott et al. (2021b); Ebersold and Regimbau (2025), as well as popstock, a Python-based module that predicts the background spectrum and its uncertainty given a CBC population and merger rate evolution Renzini and Golomb (2024).

V.2 Binary black holes

To model the BBH background we generally follow the approach in Abbott et al. (2021b) with some adjustments in the presumed mass, redshift, and spin distributions of BBHs.

In Abbott et al. (2021b) it was assumed that the BBH primary mass distribution is approximated by a mixture between a Gaussian peak near $35\,M_{\odot}$ and a power-law continuum, that mass ratios are power-law distributed, and, for simplicity, that spins are negligible. We update these assumptions to match the latest understanding of the black hole population in GWTC-4 LIGO Scientific, VIRGO, KAGRA Collaborations (2025a). Following LIGO Scientific, VIRGO, KAGRA Collaborations (2025b), we assume that black hole primary masses follow a mixture between two Gaussians, centered near $10\,M_{\odot}$ and $35\,M_{\odot}$, and a broken power-law continuum extending up to $200\,M_{\odot}$. Whereas we previously neglected spins, we now include a realistic model for the black hole spin distribution; we follow LIGO Scientific, VIRGO, KAGRA Collaborations (2025b) and assume spin magnitudes are described by a truncated normal distribution and spin orientations by a mixture between isotropic and preferentially-oriented subpopulations. Black hole mass ratios are again assumed to be power-law distributed. As in Abbott et al. (2021b), we do not fix the parameters of these models, but instead show a projection for $\Omega_{\rm BBH}(f)$ that self-consistently includes observational uncertainties in the BBH mass, spin, and mass ratio distributions.

Previous forecasts of the BBH background included uncertainties in the ensemble of intrinsic black hole properties and the local BBH merger rate, while imposing a fixed fiducial model describing how this rate evolved with redshift. In Abbott et al. (2021b), for example, we assumed a BBH birth rate proportional to low-metallicity star formation, with a power-law time delay distribution between binary formation and merger. Here, we revise our projections to additionally reflect systematic uncertainties in the BBH merger rate. We assume the BBH merger rate density follows the broken power law Fishbach et al. (2018)

[TABLE]

where parameters $\alpha_{z}$ and $\beta_{z}$ represent the low and high redshift power laws to which the merger rate history asymptotes, respectively, while $z_{\rm p}$ is the redshift at which the merger rate history peaks. Specifically, Eq. (22) implies that $R_{\rm BBH}(z)\propto(1+z)^{\alpha_{z}}$ for redshifts $z\ll z_{\rm p}$ and $R_{\rm BBH}(z)\propto(1+z)^{-\beta_{z}}$ for $z\gg z_{\rm p}$, and the normalization constant $C(\alpha_{z},\beta_{z},z_{\rm p})$ is defined such that $R_{\rm BBH}(z)=R_{0}$ at $z=0$. This parametrization allows us to cover a broad range of theoretical predictions for the redshift evolution of the BBH merger rate. It accommodates scenarios involving isolated binary evolution, where the merger rate closely follows the metallicity-dependent cosmic star formation rate with a distribution of time delays between formation and merger Dominik et al. (2013); Mapelli et al. (2017); Santoliquido et al. (2020). At the same time, it can represent dynamically assembled mergers in dense star clusters, which are less directly tied to star formation and often exhibit a merger rate that peaks at higher redshifts Mapelli et al. (2022); Kou et al. (2024).

We combine all direct detections of BBHs in GWTC-4 with our updated constraints on the astrophysical GWB to measure the underlying distributions of BBH primary masses, mass ratios, spins, and redshifts. We assume these distributions are well-described by the models listed above, and we follow the methods presented in LIGO Scientific, VIRGO, KAGRA Collaborations (2025b); Callister et al. (2020); Turbang et al. (2024); Lalleman et al. (2025). In order to self-consistently combine data from direct BBH detections with constraints on the GWB, we adopt a factorized likelihood given by

[TABLE]

Here, $p_{\rm stoch}(\hat{\Omega}_{0}|\Lambda_{\rm BBH})$ is the likelihood of our cross-correlation measurements [see Eq. (12)], given a BBH population model denoted $\Lambda_{\rm BBH}$. Meanwhile, $p_{\rm BBH}(\{d_{i}\}|\Lambda_{\rm BBH})$ is the likelihood of the direct compact binary detections comprising GWTC-4 LIGO Scientific, VIRGO, KAGRA Collaborations (2025b). In practice, rather than evaluating Eq. (23) directly, we sample the likelihood in a two step process, first sampling from the posterior $p_{\Lambda_{\rm BBH}|\rm BBH}(\{d_{i}\})$ and then rejection sampling these draws using $p_{\rm stoch}(\hat{\Omega}_{0}|\Lambda_{\rm BBH})$. The resulting $\Omega_{\rm BBH}(f)$ spectra calculated from each posterior sample $\Lambda$ after rejection sampling are shown in the left-hand side of Fig. 6. At a reference frequency of $25\,\mathrm{Hz}$, we predict $\Omega_{\mathrm{BBH}}(25\,\mathrm{Hz})=0.8^{+1.3}_{-0.5}\times 10^{-9}$ before rejecting according to the GWB upper limits. After rejection sampling, the predicted BBH GWB amplitude is $\Omega_{\mathrm{BBH}}(25\,\mathrm{Hz})=0.8^{+1.1}_{-0.5}\times 10^{-9}$. Here and henceforth, all uncertainties are quoted as 90% credible intervals.

The inferred BBH merger rate $R_{\mathrm{BBH}}(z)$ is shown in the right side of Fig. 6. Black and red curves again show 90% credible bounds and individual posterior samples, respectively. The present-day merger rate is inferred to be $R_{0}=15.4^{+6.0}_{-4.5}\,\mathrm{Gpc}^{-3}\,\mathrm{yr}^{-1}$. For comparison, dotted black curves indicate 90% credible bounds obtained previously when analyzing BBH mergers in GWTC-3 Abbott et al. (2021e, 2022). New results are consistent with previous estimates, although it is now found that the rate with which $R_{\mathrm{BBH}}(z)$ increases ($\alpha_{z}=3.2^{+1.1}_{-1.0}$) is near the upper end of previous bounds ($\alpha_{z}=2.7^{+2.2}_{-1.9}$). The high redshift parameters $\beta_{z}$ and $z_{\rm p}$ remain poorly constrained ($\beta_{z}=5.1^{+4.5}_{-4.6}$, $z_{\rm p}=2.6^{+1.2}_{-1.3}$).

The forecasted $\Omega_{\mathrm{BBH}}(25\,\mathrm{Hz})$ reported above is approximately a factor of two larger than previously estimated in Abbott et al. (2021b) and Abbott et al. (2022). This increase is due primarily to two effects. The increased sample of direct compact binary detections indicates that the BBH merger rate likely increases more quickly with redshift than previously measured. Secondly, whereas past forecasts used simple parametric estimates for the energy spectrum of a given binary Ajith et al. (2011), our updated forecasts adopt the significantly more mature IMRPhenomXP waveform model to compute radiated energy.

The non-detection of a stochastic background provides an additional constraint on the population history of BBHs, effectively ruling out scenarios with high merger rates at large redshifts that would otherwise yield a detectable background signal Callister et al. (2020). While the background upper limits effectively do not impact the inferred merger rate, these imply a modest tightening of the upper uncertainty bound on the BBH background amplitude as described above. This result implies that, when modeling the BBH merger history via Eq. (22), the direct detections already require $\alpha_{z}$ to be sufficiently small that most possible astrophysical backgrounds lie below the current detection threshold. In future observing runs, the exact degree to which stochastic search results are or are not informative will depend on the achieved sensitivity and run duration, as well as updated inference of $\alpha_{z}$.

V.3 Binary neutron stars

Since neither the local merger rate nor the mass and spin distribution of BNS mergers are well constrained from current GW catalogs, we cannot adopt detailed data-driven models of the BNS population. Instead, we forecast the astrophysical background from BNSs using the methodology of Abbott et al. (2021b), which adopts a fixed star formation history convoluted with a time delay distribution and simplified assumptions for mass and spin. We assume that BNS progenitor formation follows the star formation rate described in Madau and Fragos (2017) (Eq.(22) with $\alpha_{z}=2.6$, $\beta_{z}=3.6$ and $z_{p}=2.2$), with a distribution $p(t_{d})\sim t_{d}^{-1}$ of time delays $t_{d}$ between binary formation and merger, where 20 Myr $<t_{d}<$ 13.5 Gyr. Following LIGO Scientific, VIRGO, KAGRA Collaborations (2025b), we assume a simple fixed population model, with a uniform mass distribution between $1-2.5\,M_{\odot}$ and isotropic spin orientations with uniformly distributed spin magnitudes below 0.4. Under these assumptions, the local merger rate is estimated to be $61^{+113}_{-48}$ $\mathrm{Gpc}^{-3}\mathrm{yr}^{-1}$. When modeling $\Omega_{\mathrm{BNS}}(f)$, we do not consider any tidal effects nor contributions from a potential post-merger signal Lehoucq et al. (2025). These modeling choices lead to an estimated BNS GWB shown in Fig. 7, specifically at the reference frequency of 25 Hz, we predict $\Omega_{\mathrm{BNS}}(25\mathrm{Hz})={3.6^{+6.7}_{-2.8}\times 10^{-11}}$.

V.4 Neutron star-black hole binaries

Similar to BNSs, few compact binaries with masses consistent with NSBH systems have been observed Abbott et al. (2020, 2021f); Abac et al. (2024b), and so the details of the NSBH population remain largely unknown. As in the BNS case above, we again follow the simplified modeling approach undertaken previously in Abbott et al. (2022). The black hole masses are drawn from a log-uniform distribution between 3 and $50\,M_{\odot}$, and neutron star masses from a uniform distribution between $1-2.5\,M_{\odot}$. Black hole spins are assumed to be isotropically distributed with uniform spin magnitude below 0.99, while the neutron star spin follows the same distribution as for BNS. The energy emission from NSBH binaries is calculated using the BBH waveform model IMRPhenomXP, which provides a reliable approximation, especially for systems with large mass ratios Foucart et al. (2013). However, the model does not include tidal deformability or the possibility of tidal disruption of the neutron star, effects that can affect the emitted GW signal from binaries with more comparable component masses Pannarale et al. (2015). For the given population model, the local merger rate of NSBH binaries is estimated to be $30^{+34}_{-19}$ $\mathrm{Gpc}^{-3}\mathrm{yr}^{-1}$ LIGO Scientific, VIRGO, KAGRA Collaborations (2025b). Moreover, we assume that the rate evolves in redshift in a manner identical to BNS systems as described above.

With these assumptions we find $\Omega_{\mathrm{NSBH}}(25\mathrm{Hz})={5.0^{+5.6}_{-3.1}\times 10^{-11}}$. The 90% credible band over the whole frequency range is displayed in Fig. 7. We note that $\Omega_{\mathrm{NSBH}}(f)$ is predicted to be significantly flatter than in past estimates, deviating noticeably from a canonical $f^{2/3}$ power law above $\sim 30\,\mathrm{Hz}$. This change is primarily due to our adoption of an updated waveform model, in which the energy spectra of highly unequal-mass binaries are suppressed at high frequencies.

V.5 Total compact binary background

In Fig. 7 (right hand side), we present an updated estimate of the combined CBC GWB. Our model predicts this background at $\Omega_{\mathrm{BBH+BNS+NSBH}}(25\mathrm{Hz})={0.9^{+1.1}_{-0.5}\times 10^{-9}}$. Additionally, we show the estimated sensitivity of the search described in this paper by presenting the $2\sigma$ power-law integrated sensitivity curve including data from O1$-$O4a. Although our estimate for the background amplitude is below current limits, it may become accessible with the A+ configuration of the LIGO detectors. This may be seen from the comparison shown in Fig. 7 with the target $2\sigma$ power-law integrated sensitivity curve for the O5 data run (dashed curve), which includes substantial planned updates for the LIGO sites Capote et al. (2025); The LIGO Scientific Collaboration, Virgo Collaboration, and KAGRA Collaboration (2022) and the Virgo detector Rocchi et al. (2025), and assumes one year of uninterrupted data. It is important to note that these are target sensitivities that are expected to be reached by the end of the observing run, hence the O5 target sensitivity curve in Fig. 7 should be interpreted as an optimistic limit for the observing scenario.

VI Conclusions

In this paper, we presented an isotropic search conducted on data spanning from the first LVK observing run (O1) to the first part of the fourth observing run (O4a), incorporating 115.51 days of new coincident LIGO-Hanford and LIGO-Livingston observation. No evidence for a GWB is found, and we limit the net energy density of the GWB to $\Omega_{\rm ref}(25\,\mathrm{Hz})\leq 2.9\times 10^{-9}$ under a log-uniform prior ($\leq 1.5\times 10^{-8}$ under a uniform prior), an improvement of approximately a factor of 2 relative to previous bounds Abbott et al. (2021b).

We additionally constrain alternative theories of gravity by searching for non-standard gravitational-wave polarizations, including scalar and vector modes beyond general relativity’s tensor predictions. Our analysis improves existing bounds on these exotic polarization states by factors of two compared to previous datasets. However, the absence of any detectable background prevents us from distinguishing between standard tensor-only scenarios and models incorporating additional polarization modes.

As GWB searches improve in sensitivity, previously negligible terrestrial noise sources, including magnetic coupling from Schumann resonances and lightning activity, might contaminate our results and require characterization and possible mitigation. To quantify the potential impact of magnetic couplings on our result, we combine measurements of the magnetic environments surrounding LHO and LLO with measured transfer functions between external magnetic fields and measured strain, estimating a net magnetic noise budget. We find that, given current search sensitivities, contamination due to correlated magnetic environments does not significantly impact our results.

We conclude by examining how close our sensitivity has come to detecting the expected astrophysical background from the population of merging compact binaries across the Universe. We present updated estimates of energy-density spectra arising from distant binary black holes, binary neutron stars, and neutron-star black hole binaries, incorporating updated measurements of these sources’ merger rates and demographics. Combining these three source classes, we estimate a total astrophysical background amplitude $\Omega_{\mathrm{BBH+BNS+NSBH}}(25\mathrm{Hz})={0.9^{+1.1}_{-0.5}\times 10^{-9}}$, consistent with but marginally smaller than previous estimates Abbott et al. (2022). This decrease is primarily driven by lower inferred BNS and NSBH merger rates, given the lack of such systems detected during O4a. This estimate lies at least a factor of three below the sensitivity of our search, and is thus consistent with a current non-detection of a GWB.

Several companion papers further extend the methodology presented here and elaborate on additional implications of a GWB non-detection. In this work, we have limited ourselves to exploring a purely isotropic stochastic background. Alternate analyses targeting anisotropic backgrounds are presented in LIGO Scientific, VIRGO, KAGRA Collaborations (prepa). Additionally, although we have commented on the astrophysical background expected from merging compact binaries, there exist a host of predicted mechanisms that may yield GWBs of cosmological origin. Implications of our updated constraints on the GWB on such cosmological backgrounds are discussed in LIGO Scientific, VIRGO, KAGRA Collaborations (prepb).

VII Acknowledgements

This material is based upon work supported by NSF’s LIGO Laboratory, which is a major facility fully funded by the National Science Foundation. The authors also gratefully acknowledge the support of the Science and Technology Facilities Council (STFC) of the United Kingdom, the Max-Planck-Society (MPS), and the State of Niedersachsen/Germany for support of the construction of Advanced LIGO and construction and operation of the GEO 600 detector. Additional support for Advanced LIGO was provided by the Australian Research Council. The authors gratefully acknowledge the Italian Istituto Nazionale di Fisica Nucleare (INFN), the French Centre National de la Recherche Scientifique (CNRS) and the Netherlands Organization for Scientific Research (NWO) for the construction and operation of the Virgo detector and the creation and support of the EGO consortium. The authors also gratefully acknowledge research support from these agencies as well as by the Council of Scientific and Industrial Research of India, the Department of Science and Technology, India, the Science & Engineering Research Board (SERB), India, the Ministry of Human Resource Development, India, the Spanish Agencia Estatal de Investigación (AEI), the Spanish Ministerio de Ciencia, Innovación y Universidades, the European Union NextGenerationEU/PRTR (PRTR-C17.I1), the ICSC - CentroNazionale di Ricerca in High Performance Computing, Big Data and Quantum Computing, funded by the European Union NextGenerationEU, the Comunitat Autonòma de les Illes Balears through the Conselleria d’Educació i Universitats, the Conselleria d’Innovació, Universitats, Ciència i Societat Digital de la Generalitat Valenciana and the CERCA Programme Generalitat de Catalunya, Spain, the Polish National Agency for Academic Exchange, the National Science Centre of Poland and the European Union - European Regional Development Fund; the Foundation for Polish Science (FNP), the Polish Ministry of Science and Higher Education, the Swiss National Science Foundation (SNSF), the Russian Science Foundation, the European Commission, the European Social Funds (ESF), the European Regional Development Funds (ERDF), the Royal Society, the Scottish Funding Council, the Scottish Universities Physics Alliance, the Hungarian Scientific Research Fund (OTKA), the French Lyon Institute of Origins (LIO), the Belgian Fonds de la Recherche Scientifique (FRS-FNRS), Actions de Recherche Concertées (ARC) and Fonds Wetenschappelijk Onderzoek - Vlaanderen (FWO), Belgium, the Paris Île-de-France Region, the National Research, Development and Innovation Office of Hungary (NKFIH), the National Research Foundation of Korea, the Natural Sciences and Engineering Research Council of Canada (NSERC), the Canadian Foundation for Innovation (CFI), the Brazilian Ministry of Science, Technology, and Innovations, the International Center for Theoretical Physics South American Institute for Fundamental Research (ICTP-SAIFR), the Research Grants Council of Hong Kong, the National Natural Science Foundation of China (NSFC), the Israel Science Foundation (ISF), the US-Israel Binational Science Fund (BSF), the Leverhulme Trust, the Research Corporation, the National Science and Technology Council (NSTC), Taiwan, the United States Department of Energy, and the Kavli Foundation. The authors gratefully acknowledge the support of the NSF, STFC, INFN and CNRS for provision of computational resources. This work was supported by MEXT, the JSPS Leading-edge Research Infrastructure Program, JSPS Grant-in-Aid for Specially Promoted Research 26000005, JSPS Grant-in-Aid for Scientific Research on Innovative Areas 2402: 24103006, 24103005, and 2905: JP17H06358, JP17H06361 and JP17H06364, JSPS Core-to-Core Program A. Advanced Research Networks, JSPS Grants-in-Aid for Scientific Research (S) 17H06133 and 20H05639, JSPS Grant-in-Aid for Transformative Research Areas (A) 20A203: JP20H05854, the joint research program of the Institute for Cosmic Ray Research, University of Tokyo, the National Research Foundation (NRF), the Computing Infrastructure Project of the Global Science experimental Data hub Center (GSDC) at KISTI, the Korea Astronomy and Space Science Institute (KASI), the Ministry of Science and ICT (MSIT) in Korea, Academia Sinica (AS), the AS Grid Center (ASGC) and the National Science and Technology Council (NSTC) in Taiwan under grants including the Science Vanguard Research Program, the Advanced Technology Center (ATC) of NAOJ, and the Mechanical Engineering Center of KEK.