4-OGC: Catalog of Gravitational Waves from Compact Binary Mergers

We present the fourth Open Gravitational-wave Catalog (4-OGC) of binary neutron star (BNS), binary black hole (BBH), and neutron star–black hole (NSBH) mergers. The catalog includes observations from 2015 to 2020 covering the first through third observing runs (O1, O2, O3a, and O3b) of Advanced LIGO and Advanced Virgo. The updated catalog includes seven BBH mergers that were not previously reported with high significance during O3b for a total of 94 observations: 90 BBHs, 2 NSBHs, and 2 BNSs. The most confident new detection, GW200318_191337, has component masses 49.1−12.0+16.4M⊙ and 31.6−11.6+12.0M⊙; its redshift of 0.84−0.35+0.4 (90% credible interval) may make it the most distant merger so far. We estimate the merger rate of BBH sources, assuming a power-law mass distribution containing an additive Gaussian peak, to be 16.5−6.2+10.4(25.0−8.0+12.6) Gpc−3 yr−1 at a redshift of z = 0 (0.2). For BNS and NSBH sources, we estimate a merger rate of 200−148+309 Gpc−3 yr−1 and 19−14+30 Gpc−3 yr−1, respectively, assuming the known sources are representative of the total population. We provide reference parameter estimates for each of these sources using an up-to-date model accounting for instrumental calibration uncertainty. The corresponding data release also includes our full set of subthreshold candidates.


INTRODUCTION
Gravitational-wave astronomy has entered an era of regular and routine observation of compact-binary mergers.This achievement was made possible by the secondgeneration gravitational-wave observatories, led by the twin Advanced LIGO (Hanford and Livingston) (Aasi et al. 2015) and Advanced Virgo (Acernese et al. 2015) observatories, which have been operating since 2015 and 2017, respectively.This period has seen continued improvement in their astrophysical reach over their three completed observing runs (O1-O3) (Abbott et al. 2020b), with the pace rapidly increasing from 3 merger observations in O1 to dozens in the first half of O3 (O3a) (Abbott et al. 2019a(Abbott et al. , 2020c(Abbott et al. , 2021a;;Venumadhav et al. 2019a;Nitz et al. 2019cNitz et al. , 2021a)); the vast majority of these are binary black hole mergers (BBHs).To date, there is only a single binary neutron star (BNS) observation, GW170817 (Abbott et al. 2017a), which has been corroborated by extensive electromagnetic observations (Abbott et al. 2017b).In addition, GW190425 is a potential heavy BNS merger (Abbott et al. 2020a) and recently two sources with masses compatible with merging neutron star-black hole binaries have been reported (Abbott et al. 2021b).The plethora of BBH observations, in addition to exceptional events such as GW190521 with total mass ∼ 150M (Abbott et al. 2020d;Capano et al. 2021), are beginning to constrain formation scenarios (Abbott et al. 2020e;Gerosa & Fishbach 2021;Edelman et al. 2021;Zevin et al. 2021) and deviations from general relativity (Abbott et al. 2021c;Wang et al. 2021).
We expect the current observatories to continue to improve in sensitivity over the next few years (Abbott et al. 2020b) due to ongoing active commissioning and the inclusion of new technologies into the existing sites (Aasi et al. 2013).The upcoming O4 observing run is scheduled to begin at the end of 2022 with a fiducial BNS range of 160-190 Mpc (Abbott et al. 2020b).In addition, we can expect that a fourth ground based gravitationalwave observatory, KAGRA, will join the O4 observing run (Akutsu et al. 2020).In the mid-to-late 2020's a fifth observatory will join the worldwide network with the construction of LIGO India (Unnikrishnan 2013;Saleem et al. 2021).Efforts are also advancing for thirdgeneration observatories, Cosmic Explorer and the Einstein Telescope, which promise an order of magnitude sensitivity improvement (Reitze et al. 2019;Evans et al. 2021;Punturo et al. 2010).
This work provides a comprehensive catalog of gravitational-wave observations from merging BNS, BBH, and NSBH sources.The analysis is based on a deep archival search using the public data from the LIGO and Virgo observatories (Vallisneri et al. 2015;Abbott et al. 2021d) which spans 2015-2020 and includes all existing observing runs (O1-O3).This catalog updates the results of our previous 3-OGC catalog (Nitz et al. 2021a) by including analysis of data from the second half of the third observing run (O3b) which very recently became public.The next expected public data release would be at the end of 2024 if the current release delay of 18 months with a 6-month cadence is maintained.Included in our companion data release is the complete set of sub-threshold candidates and parameter estimates for significant candidates (Nitz & Kumar 2021).We make available our sub-threshold candidates so that they may aid follow-up studies, including those which cross-correlate candidates with other archival datasets such as from gamma-ray bursts (Burns et al. 2019;Nitz et al. 2019b), high-energy neutrinos (Countryman et al. 2019), or optical transients (Andreoni et al. 2018;Setzer et al. 2018).Archival analyses have the potential to uncover distant or faint populations.
We find a total of 94 mergers which pass our significance threshold, P astro > 0.5 or false alarm rate (FAR) less than once per 100 years.The vast majority, 90, of these are BBH mergers; 7 are reported during O3b for the first time here with high significance, 3 additional were previously reported marginal candidates in the 3-OGC catalog (Nitz et al. 2021a).We find the previouslyreported BNS and NSBH mergers (Abbott et al. 2017a(Abbott et al. , 2020a(Abbott et al. , 2021b)), however, no new confident BNS or NSBH mergers are observed.Our results are broadly consistent with the recent GWTC-3 catalog produced by the LVK (LIGO-Virgo-KAGRA) collaborations (Abbott et al. 2021e,f).

SEARCH FOR COMPACT-BINARY MERGERS
Our catalog includes results from the analysis of the complete set of public LIGO and Virgo data collected over a period of five years (Vallisneri et al. 2015;Abbott et al. 2021d).To identify the signature of gravitational waves from compact-binary mergers, we use matched filtering to extract the signal-to-noise ratio (SNR) of a potential signal (Allen et al. 2012).Candidates are identified by looking for peaks in the single-detector SNR time series and correlating these triggers between observing detectors.Each candidate is assessed for consistency between the expected signal morphology and the data (Allen 2005;Nitz 2018) and ranked using additional factors such as the observed rate of triggers, multidetector coherence (Nitz et al. 2017;Davies et al. 2020) and the local data quality (Mozzon et al. 2020;Abbott et al. 2018).This procedure is the same as that used for the prior 3-OGC catalog (Nitz et al. 2021a) and is implemented using the open-source PyCBC toolkit (Nitz et al. 2021b;Usman et al. 2016).This toolkit has also been used to construct many matched-filter based gravitational-wave searches (Abbott et al. 2019a(Abbott et al. , 2020c, 2021e,g), 2021e,g), including those that account for eccentricity (Nitz et al. 2019a), are focused on intermediate-mass BBHs (Chandra et al. 2021), or target subsolar-mass mergers (Nitz & Wang 2021a,b,c).

LIGO and Virgo Observing Period
The public LIGO and Virgo data spans the time period from 2015-2020 and the three completed observing runs O1-O3 (Vallisneri et al. 2015;Abbott et al. 2021d).As in 3-OGC (Nitz et al. 2021a), we include results from data released outside the nominal observation runs, notably this includes data around GW170608 (Abbott et al. 2017c) and GW190814 (Abbott et al. 2020f).Table 1 shows the number of days that different network configurations were observing.Our updated catalog includes the recently released second half of the third observing run (O3b).The full O3 observing run contains 152 days of triple-detector observing time and ∼ 200 days when both LIGO observatories were observ-Table 1. Analyzed time in days for different global network observing scenarios.The abbreviations H, L, and V are used for the LIGO-Hanford, LIGO-Livingston, and Virgo observatories, respectively.Each time period is exclusive of the others.Some data is excluded from the full public dataset due to analysis requirements (O(1)%).However, data around GW170608 and GW190814 which was made available separately from the bulk data release is included (Vallisneri et al. 2015).ing.We analyze all data when multiple detectors are observing and also identify sources in data when only LIGO-Livingston or LIGO-Hanford are observing.Due to the decreased range and relatively large population of non-Gaussian transient noise, we do not consider data when only Virgo is observing.In Fig. 1 we show the sky-averaged fiducial BNS range of each instrument as a function of time.All instruments have made significant gains in sensitivity during O3 as compared to O2, with Virgo reaching 60 Mpc range and LIGO-Livingston periodically exceeding 140 Mpc.
The data has been calibrated by the LVK collaborations and made available through the Gravitationalwave Open Science Center (GWOSC) (Viets et al. 2018;Acernese et al. 2018;Bhattacharjee et al. 2021;Estevez et al. 2021;Vallisneri et al. 2015;Abbott et al. 2021d).Noise subtraction using auxiliary witness channels is ap-plied to the bulk data release (Davis et al. 2019;Vajente et al. 2020;Estevez et al. 2019;Rolland et al. 2019).Our analysis also makes use of the data quality information compiled by the LVK detector characterization teams (Davis et al. 2021) to exclude times around hardware signal injections and times affected by adverse instrumental behavior.

Search Space
Our search uses matched-filtering to identify signals within the data; matched filtering requires a model of the expected signal to filter the data.To identify sources within a broad target region of parameters (component masses m 1,2 and spins), we use a discrete bank of template waveforms.For this catalog, we use the same template bank, shown in Fig. 2, that we previously used in the 3-OGC analysis (Nitz et al. 2021a).
The templates are chosen to ensure that a minimum fraction of the optimal SNR of a potential signal is recovered; typically O(97%) for many searches (Dal Canton & Harry 2017).The bank is constructed in 4 parts using a brute-force stochastic algorithm (Harry et al. 2009;Ajith et al. 2014).These regions include a focused low-mass-ratio, low-spin, BBH region with a target SNR recovery > 99.5%.In addition, there are BNS, NSBH, and broad-parameter BBH banks which target SNR recovery > 97%.
The template bank is designed to recover gravitational-wave signals from non-precessing quasicircular sources.Our search accounts only for the effects of the dominant-mode gravitational-wave signal and does not include the effects of higher-order modes.Neglecting these effects will reduce the search sensitivity to sources which strongly exhibit these features, such as for highly-inclined, high-mass-ratio, or highly precessing binaries (Harry et al. 2016) or where there remains residual eccentricity (Ramos-Buades et al. 2020;Wang & Nitz 2021).Development of optimal search strategies for these sources is an ongoing endeavor (Harry et al. 2018(Harry et al. , 2016) ) and techniques that don't rely on matched filtering also target these sources (Klimenko et al. 2008(Klimenko et al. , 2016;;Tiwari et al. 2016).To model the gravitational-wave signal, we use a combination of Tay-lorF2 (for BNS sources) (Sathyaprakash & Dhurandhar 1991;Droz et al. 1999;Blanchet 2002;Faye et al. 2012), SEOBNRv4 ROM (for BBH and NSBH) (Taracchini et al. 2014;Bohéet al. 2016), and IMRPhenomD (focused BBH) (Husa et al. 2016a;Khan et al. 2016b).

Candidate Selection and Significance
Candidates are assigned a ranking statistic value according to their SNR, consistency with the expected signal morphology (Allen 2005;Nitz 2018), and coherence between observing detectors (Nitz et al. 2017;Davies et al. 2020).The statistical significance of multi-detector candidates is assessed by comparing to empirically estimated background from artificially produced analyses.By construction these cannot contain astrophysical sources as they are produced by time-shifting the data from one or more detectors by a constant greater than the light-travel-time between detectors (Babak et al. 2013;Usman et al. 2016).This technique has been used in many past analyses (Nitz et al. 2019c;Venumadhav et al. 2019b;Abadie et al. 2012;Abbott et al. 2009Abbott et al. , 2020c, 2021e), 2021e).For multi-detector candidates, the estimated background distribution is used to establish the false alarm rate (FAR) of the search as a function of ranking statistic value.Candidates detected in a sin- Signals often have multiple templates which will produce a candidate.Here, we only show the template which produced a candidate with the lowest false alarm rate; the parameters of the selected template can only be considered as crude point estimate of the true parameters, so they may differ significantly from the posterior estimates in Table 5.
gle detector are instead assessed against off-source observation time.In particular, we use time when both LIGO detectors are observing which allows for confident multi-detector observations to be excised from the background, minimizing potential signal contamination.
To limit the effects of non-stationary noise and non-Gaussian noise transients, we further restrict our singledetector analysis to candidates which arise from either the BNS region, focused-BBH with chirp mass M < 60, or NSBH with total mass M < 50.The highest mass templates of each region are the most difficult to distinguish from non-Gaussian transient noise due to their short duration.The probability of astrophysical origin P astro is calculated using the empirically measured background and comparing to the distribution of observed ranking statistic values produced by a simulated source population as part of a two-component mixture model (Farr et al. 2015b).For the single-detector analysis, the background distribution is extrapolated using the method of (Nitz et al. 2020).This is the same methodology as previously used the 3-OGC analysis.We limit our assignment of P astro to the single-detector analyses using the method of (Nitz et al. 2020) and multi-detector candidates from the focused-BBH region, where the vast majority of ob-servations are found.Due to the uncertain population distribution in other regions, we do not assign a probability of astrophysical origin.
In the 3-OGC analysis, we implicitly assumed a detection-prior that is flat in redshifted chirp-mass (Nitz et al. 2021a) within the focused-BBH region; we make a marginal improvement to this step by averaging the obtained astrophysical probability over this fiducial population scenario with the one obtained from the smoothed observed distribution of redshifted chirp-masses of the high confidence observations (FAR < 1 per 100 yrs).Due to the broad observed distribution in chirp masses, this produces only very modest changes to the estimated astrophysical probability compared to those previously reported in 3-OGC; the primary impact is to mildly increase support for sources with M <∼ 50 in lieu of higher mass sources.A future improvement for the assessment of marginal candidates would be to include a physical model that can simultaneously fit the population along with the foreground and background distributions.

OBSERVATIONAL RESULTS
All events identified by our search with a probability of astrophysical origin P astro > 0.5 or inverse false alarm rate (IFAR) > 100 years are listed in Table 3.We report 94 events, 90 of which are BBHs, 2 originate from BNS mergers, and 2 are the result of coalescing NSBHs.The table contains information on the observational status of the three different instruments LIGO Hanford, LIGO Livingston, and Virgo, as well as the recovered SNR for the detectors which triggered on the event.Finally, we also provide the Global Positioning System (GPS) times, estimates of the P astro for BBH signals, and the observed IFARs.
Our search finds 7 new BBH signals during O3b which pass our significance criteria and were not reported as marginal or confident observations by previous studies (Abbott et al. 2021e).These events are GW191224 043228, GW200106 134123, GW200129 114245, GW200210 005122, GW200214 223307, GW200305 084739, and GW200318 191337.In addition, 3 previously marginal BBH candidates from the previous 3-OGC catalog (Nitz et al. 2021a) now pass our criteria; this is due to the updated prescription for calculating P astro (see sec 2.3).No new BNS or NSBH were discovered at high confidence.
The most confident new detection is GW200318 191337 with a P astro = 0.97.
It is detected in both LIGO Hanford and LIGO Livingston with SNR ≤ 6.2.GW200318 191337 may be the most distant observed merger with redshift z∼ 0.84 +0.4  −0.35 .GW200129 114245 is a potential new mass ratio ∼ 2 BBH with primary mass 79.1 +40.2 −37.6 M .GW200129 114245 is detected by both LIGO Hanford and LIGO Livingston with SNR ≤ 6 and occurs less than 5 hours after the previous event GW200129 065458.
The results of our single-detector analysis is shown in Fig. 3 for BNS, BBH, and now NSBH sources.There have now been 8 mergers detected in a single detector: 6 BBHs, 1 BNS, 1 NSBH.The previously reported NSBH merger GW200105 162426 (Abbott et al. 2021b) is found as a near threshold detection in our analysis.It was observed only in LIGO Livingston with a SNR of 13.3.LIGO Hanford was not observing at the time.It has P astro ∼ 0.5, similar to GW190425, due to its wide separation from the background distribution.The assumed merger rate for similar NSBH sources is determined by the inclusion of GW200115 whose singledetector LIGO Hanford trigger had a ranking statistic marginally larger than all collected single-detector background.
We present 30 marginal sub-threshold events following the criterion P astro > 0.2 or IFAR > 0.5 years in Table 4.The majority of these are consistent with BBH mergers except for 170722 065503 and 191219 163120 which are consistent with a BNS and a NSBH merger, respectively.In total 12 marginal candidates were previously reported in our most recent 3-OGC catalog; these are 151011 192749, 170425 055334, 170704 202003, 190426 053949, 190509 004120, 190530 030659, 190630 135302, 190704 104834, 190707 071722, 190808 230535 and 190821 050019.The remainder of our sub-threshold candidates are included in the corresponding data release (Nitz & Kumar 2021).
We visually inspect spectrograms of the data around each marginal event to check for non-stationary behavior or an excess of non-Gaussian transient noise which could adversely affect the identified candidate.We observed excess power consistent with a blip glitch (Cabero et al. 2019) in several cases at off-source times.In a few other events we observe slight deviations in power from the expected stationary Gaussian noise in the frequency range 60-110 Hz, e.g. in the Hanford data a few seconds before the time of events (e.g.200102 095606 or 200310 090144).Further investigation would be required to determine if potential candidates could be adversely affected in these cases.
Our results are otherwise broadly consistent with those of the recently published GWTC-3 catalog (Abbott et al. 2021e).There are eight events reported in GWTC-3 which do not pass our candidate threshold.These are GW191103 012549, GW191113 071753, GW191219 163120, GW200208 222617, GW200210 092254,  2020) is used to extrapolate the background distribution and estimate the probability of astrophysical origin.The BBH analysis uses the statistic λs (Nitz et al. 2020) while the others use a re-weighted SNR statistic (Babak et al. 2013;Nitz et al. 2019c).GW190425 appears in both the BNS and NSBH analysis as candidates were produced from both regions.
GW200220 061928, GW200220 124850, and GW200322 091133.Both potential new NSBH candidates GW191219 163120 and GW200210 092254 from GWTC-3 are not assigned high confidence by our search.GW191219 163120 and GW200220 124850 are detected as near threshold triggers with a P astro > 0.2 or IFAR > 0.5 years.All other events are assigned a lower significance.We expect these small differences in the population of near-threshold events are consistent with differing analysis choices.We note that our candidate threshold is marginally more conservative; GWTC-3 includes candidates with P astro > 0.5 in any of several analyses (Abbott et al. 2021e).

Binary Source Parameters
In order to infer the properties of the observed compact-binary mergers, we use PyCBC inference (Biwer et al. 2019) to perform Bayesian parameter estimation with the standard likelihood assuming the detector noise to be stationary, Gaussian and uncorrelated between the detectors.The parameter estimation results for a total of 94 events are summarized in Fig. 4 and Table .5. for all events that pass our detection criteria in 4-OGC.The 5th and 95th quantile values are marked with a bar, respectively.Different colors are used to aid associating each event with its posterior estimates.
For estimating the source parameters of the BBHs, we use IMRPhenomXPHM (Pratten et al. 2020; LAL-Suite 2020), which models precessing binaries and includes higher-order modes.For the two NSBH observations, we use the IMRPhenomNSBH (Thompson et al. 2020) waveform that models a non-precessing NSBH with the neutron star mass, m NS ≤ 3M , tidal deformability Λ ∈ [0, 5000] and the mass ratio no higher than 100 (LALSuite 2020).For the BNS mergers, we use the IMRPhenomD NRTidal waveform model (Khan et al. 2016a;Husa et al. 2016b;Dietrich et al. 2017Dietrich et al. , 2019;;LALSuite 2020).This model includes the two tidal deformabilities, Λ 1 and Λ 2 of the two components, which characterize the neutron star equation-of-state.A dynamical nested sampler dynesty (Speagle 2020) is used to sample over the intrinsic parameters mass m 1,2 , spin s 1,2 (in the case of BNS, also tidal deformabilities Λ 1,2 , in the case of NSBH, we include the NS tidal deformability Λ 2 ), and extrinsic parameters luminosity distance d L , inclination angle ι, polarization angle Ψ, right ascension α, declination δ, coalescence time t c , and phase φ c .For BBH sources, the likelihood is marginalized over numerically over Ψ, while for BNS and NSBH sources, it is analytically marginalized over φ c .
We also take into account the effect of calibration uncertainties in amplitude and frequency on the parameter estimation of each event.We use the calibration uncertainty envelopes obtained from the GWOSC for previously known mergers (Vallisneri et al. 2015).For the new BBH events, we use the calibration uncertainty envelope from the nearest available time, assuming the calibration envelopes do not change rapidly.The calibration errors in amplitude and phase of the data in each detector are described by frequency-dependent splines (Farr et al. 2015a) whose parameters are directly sampled over in our analysis.We find that the estimation of intrinsic parameters is not significantly affected by the calibration uncertainties.
The priors for the intrinsic and extrinsic parameters, the power spectral density (PSD) estimation methods as well as the sampler settings are broadly consistent with those used in the 3-OGC analysis (Nitz et al. 2021a).For the BBH events, uniform priors on source-frame component masses and merger time are used.For the spins, we use uniform priors for the magnitude of the spin and isotropic for the orientation.The distance prior is assumed to be uniform in comoving volume assuming a flat ΛCDM cosmological model (Ade et al. 2015).For the rest of the extrinsic parameters, we consider an isotropic prior distribution in the sky localization and binary orientation.
For the two NSBH systems, the prior on the tidal deformability of the neutron star is considered to be uniform in the range [0, 5000].Other prior distributions are the same as BBH events except that the NS spin magnitude is uniform in [0, 0.05] and the BH is uniform in [0, 0.5].This latter range is chosen based on the parameter space where the waveform model is valid (Thompson et al. 2020;Abbott et al. 2021b).
In the case of BNS mergers, we use the same prior distributions for the component masses, comoving volume, merger time, and orientation, as used in estimating the BBH parameters.However, for GW170817 124104 we fix the sky location to that from the observed electromagnetic counterpart, α = 3.446 and δ = −0.408Abbott et al. (2017b).The spin magnitudes for the two components are uniform in the range [0, 0.05].We also vary the tidal deformability parameters of the NS components, Λ 1 and Λ 2 independently in the range [0, 5000].

Binary Black Holes
A total of 90 BBH mergers have now been observed which pass our significance threshold.The growing population will eventually enable distinguishing possible formation channels and the potential contribution from each (Talbot & Thrane 2018;Abbott et al. 2019b;Roulet et al. 2020;Abbott et al. 2020g, 2021f).Notable features expected of the BBH population were the lower (Gupta et al. 2020;Zevin et al. 2020) and upper (Fishbach & Holz 2020;Edelman et al. 2021) mass gaps, the former is an observational expectation from X-ray binaries (Farr et al. 2011;Özel et al. 2010;Bailyn et al. 1998), while the latter is expected due to pairinstability supernovae (Woosley 2017;Marchant et al. 2019;Stevenson et al. 2019;van Son et al. 2020).
Several events have components plausibly within the lower mass gap, notably GW190427 180650, GW190725 174728, GW190924 021846, GW190930 133541, GW200210 005122 and GW200316 215756 have posterior support for the secondary mass below 5 M .GW190814, however, has secondary mass well constrained to be within this region with secondary mass 2.6 +0.2 −0.1 M .The upper mass gap is expected to occur at ∼ 50-120 M (Yoshida et al. 2016;Woosley 2017;Belczynski et al. 2016;Marchant et al. 2019;Woosley 2019;Stevenson et al. 2019;van Son et al. 2020).Several detections, including GW190521 (Abbott et al. 2020d), have primary masses within this range.One explanation is that these may be second-generation mergers (Rodriguez et al. 2019;Kimball et al. 2020;Gerosa & Fishbach 2021).Other exceptional observations include the high-mass-ratio GW190814 merger (q ∼ 9); several proposals exist to explain its forma- The observed distribution of source parameters (e.g.masses and spins) and redshift distribution from a large catalog of compact-binary mergers can be used to constrain the population models that predict intrinsic source distribution as well as the rate of mergers (Abbott et al. 2019bRoulet et al. 2021;Zhu et al. 2021a).In figure 5, we show the distribution of component masses obtained from stacking all the posteriors marginalized over all other parameters, with and without taking into account the zeroth order selection effect that arises due to the loudness of the signal as a function of source parameters.To infer the source population, we use the same method as in the 3-OGC analysis Nitz & Capano (2021).We combine the posterior samples from each event to obtain large mass samples.We assume the same priors that we used in the parameter estimation.To account for the loudness of the signal, we estimate the comoving volume corresponding to the "horizon distance" of a posterior sample and reassign each sample a weight proportional to the inverse of the comoving volume.The horizon distance is defined as the maximum distance an optimally oriented source can be detected.

BBH Population and Merger Rate
We use the posterior samples from the individualsource parameter estimation study to infer the distributions of source-frame masses and redshifts for BBHs.The intrinsic distribution of BBH mergers can be parametrized in terms of the 'hyper-parameters' Θ pop of a population model (Thrane & Talbot 2020;Vitale et al. 2020;Abbott et al. 2020h, 2021f).The hyper parameters describe the shape of the distribution in ob-served parameters θ such as masses, spins, redshift, etc.To compare the inferred BBH distribution for masses and redshift with other catalogs, we demonstrate using models which were also employed in GWTC-3 studies (Abbott et al. 2021f).The sensitivity of our search is estimated by analyzing a simulated population that is injected throughout the O1-O3 data.Our fiducial population includes sources that extend up to a total mass of 300 M and maximum mass ratio q of 30 with constraints on individual component masses to be between 2M − 150M .In the GWTC-3 analysis, the search sensitivity is estimated semi-analytically for

Parametric Models
The observed mass distribution corrected for firstorder selection effects (as shown in figure 5) clearly shows a peak in the primary mass between 30-40 M ; this peak was previously reported in Abbott et al. (2021f).This observation motivates a parametric model where the primary mass is modelled with a powerlaw component along with added Gaussian components to account for peaks in the mass distribution.The powerlaw component of the primary mass is modelled as p(m 1 ) ∝ m α 1 .A peak in the observed distribution can be modelled with a Gaussian distribution following Ab-bott et al. (2020h, 2021f): where λ represents fraction of events in Gaussian component, C pl is the normalization for power-law component, N (x : µ M , σ M ) is the probability density function for a normal distribution for mean µ M and standard deviation σ M .An additional smoothing function is applied at the low end of the mass distribution in equation 1 i.e. m min <= m 1 < m min + δ M , where δ M is the length of the smoothing window.The other population parameters for this model are the minimum m min and maximum m max mass.The conditional probability distribution for the secondary mass m 2 is parameterized in terms of mass ratio q pl = m 2 /m 1 such that p(q pl |m 1 ) ∝ q β pl with same smoothing function at lower end of m 2 .We distinguish this definition of mass ratio q pl with a subscript 'pl' (power-law) from the other definition q = m 1 /m 2 used elsewhere in the paper.Any additional peak in the primary mass distribution can be added as additional Gaussian component in 1.
The redshift evolution of the BBH merger rate in local comoving time and space coordinates is also given as a power law (Fishbach et al. 2018) where R 0 defines the local merger rate density at z = 0 and κ represents the power law index.

Selection Effects
A search's sensitivity to BBH mergers depends on various factors such as intrinsic source parameters (such as component masses, spins), extrinsic parameters (such as orientation and sky location of the binary), the detector sensitivity, and behavior of the search to non-Gaussian transient noise.To infer the intrinsic population distribution, one needs to calculate the detection efficiency of the search i.e. the fraction of the events that can be detected with a population described by the parameters Θ pop .To estimate the detection efficiency, we add a known fiducial population to the data across the observing time of the detector network and use our search to identify them.We use an up-to-date waveform model which includes higher modes and precession, IMRPhe-nomXPHM (Pratten et al. 2020;LALSuite 2020), to simulate the fiducial population.The reference population follows a power-law mass model on m 1 and mass ratio q pl where p(m 1 ) ∝ m −α 1 and p(q pl |m 1 ) ∝ (q pl ) β with α = 2.35, β = 0. We use an isotropic spin distribution and the redshift evolution model p(z) ∝ 1 1+z dVc dz (1 + z) κ with κ = 0.The selection function is estimated using a Monte Carlo integral using the sources identified by the search as described in (Tiwari 2018;Farr 2019;Mandel et al. 2019).
where N inj are the total number of simulated sources drawn from the fiducial population, N f ound inj are the number of sources found by the search, p draw (θ k ) is probability to draw a source from the reference population with parameters θ k while π(θ k |Θ pop ) represents the conditional probability distribution of the parameter θ k given the population hyper parameter Θ pop .

Parameter Estimation for Population Parameters
Given the total number of detections N det , data d, population model described by hyper-parameters Θ pop , the population likelihood function can be described by inhomogenous Poisson processes given as (Thrane & Talbot 2020;Mandel et al. 2019;Loredo 2004;Abbott et al. 2021f), where θ is the set of individual source parameters.L(d i |θ) is the likelihood function of i-th observation.N exp is the expected number of detections over the full observation period.The total number of mergers N during the same period is given by N = N exp /ξ(Θ pop ).
We select candidates which pass the following selection criteria: IFAR > 100 years or P astro > 0.9.A total of 69 BBH events in 4-OGC catalog satisfy this criteria.We further identify GW190814 211039 as an outlier from the observed population and exclude it from further analysis due to its very high mass ratio and very low secondary component mass.We use flat priors on all population hyper-parameters except for a log uniform prior on rate parameter R 0 .We allow the maximum primary mass to vary up to 150M .We use three models for primary mass: i) pure powerlaw component ii) powerlaw + peak , and iii) powerlaw + two peaks.We use a uniform prior for the location of the first peak between 20M − 45M .For the additional peak in the powerlaw+peak model, we use a uniform prior between 45M − 80M .
Figure 6 shows the estimated differential merger rate as a function of the primary mass.The powerlaw+peak model shows a clear peak between 30M − 40M .The model with two peaks is consistent and does not indicate the presence of a second peak.The results also broadly agree with the GWTC-3 analysis (Abbott et al. 2021f).The small differences can be attributed to the selected observations, differences in the selection function of both searches, and the explored prior ranges.For example, we allow the maximum mass parameter m max to vary up to 150 M compared to 100 M in GWTC-3 analysis (Abbott et al. 2021f).This explain the differences in the differential merger rate plot at high masses (> 100 M ).
Figure 7 shows the redshift evolution of the BBH merger rate for the powerlaw+peak model.Since the majority of mergers are found near redshift z ∼ 0.2, as Table 2. BBH merger rates (90% credible interval) for various parametric models at redshift z=0 and 0.2 are shown.We use the parametrization for evolution of merger rates as redshift to be R(z) = R0(1 + z) κ .The numbers are quoted for z = 0.2 along with z = 0 as the constraints on R(z) are more stronger at z=0. expected, the merger rate is best constrained at this redshift.Table 2 shows the merger rates estimated for each model at redshifts z = 0 and at z = 0.2.

Neutron Star Binaries and Neutron Star Black Hole Binaries
GW170817 124104 (Abbott et al. 2017a) and GW190425 081805 (Abbott et al. 2020a) remain the only confidently observed mergers with BNS-compatible masses.GW170817 124104 is observed in both the LIGO Hanford and Livingston data with a joint SNR ∼ 33.This remains the only BNS observation with unambiguous electromagnetic counterparts (Abbott et al. 2017b,a).
In addition to the γ-ray burst (GRB 170817A) (Goldstein et al. 2017;Savchenko et al. 2017;Abbott et al. 2017d), the successful electromagnetic follow-up campaign that identified the associated electromagnetic transient near the galaxy NGC 4993 (Abbott et al. 2017b), supports the interpretation of this event as a neutron star merger.GW190425 081805 is observed only in the LIGO Livingston data with SNR ∼ 11.9.However, the long duration of the signal increases the power of signal consistency tests (Usman et al. 2016); given the distinct separation from the measured background distribution, we consider this a firm detection.It's masses 1.7 +0.1 −0.1 M and 1.6 +0.1 −0.1 M (assuming spin magnitude < 0.05) are consistent with interpreting both components as neutron stars.The primary mass is less than the heaviest observed neutron star ∼ 2.2M (Cromartie et al. 2019), however, galactic BNS observations typically have lower mass (zel & Freire 2016).If light black holes in these masses are produced in abundance and form binaries, then we cannot necessarily rule out this explanation from gravitational-wave observation alone.
We observe two events GW200105 162426 and GW200115 042309 with masses consistent with an NSBH system.GW200105 162426 is observed only in LIGO Livingston data with an SNR ∼ 13.3.Whereas GW200115 042309 is observed in both Livingston and Hanford data with IFAR > 100 yr.We are unable to constrain the tidal deformability of the secondary component, which is consistent with previous results (Abbott et al. 2021b;Zhu et al. 2021b).

Merger Rate
For neutron star binaries, we use a simplified approach to merger rate estimation.We assume that mergers are Poisson distributed and use a Jeffrey's prior for the rate parameter.We make the simplifying assumption that the observed sources are representative of the mass distribution and do not attempt to fit a mass distribution to the two BNS and two NSBH observations.The observed volume-time for each source class is estimated by simulating the recovery of sources with the same observed source-frame masses and observed low effective spins.
We find that if we consider both GW170817 and GW190425 to be representative members of a common BNS population, we infer the merger rate to be 200 +309 −148 Gpc −3 yr −1 using 90% credible intervals.However, given the unexpectedly large mass of GW190425, it is useful to consider the limits from GW170817 alone under the assumption that these sources either arise from different mechanisms or GW190425 is not a BNS merger.In this scenario, we find that the BNS merger rate is instead 139 +319 −118 Gpc −3 yr −1 .Assuming both GW200105 and GW200115 are NSBH mergers, we estimate the merger rate to be 19 +30 −14 Gpc −3 yr −1 .The merger rate estimation is dominated by the large uncertainties due to the few observed events and the choice of mass distribution.Alternate choices for the assumed mass distribution can lead to different results; a larger population of observed sources is required to disentangle current rate estimates from assumed mass distributions.

DATA RELEASE
The compilation of supplementary materials is available at https://github.com/gwastro/4-ogc(Nitz & Kumar 2021).To aid in follow-up analyses, the data release contains both the gravitational wave events included in this paper and also the detailed set of subthreshold BNS, NSBH, and BBH candidates.Auxiliary information such as the parameters of an associated waveform template, ranking statistics, false alarm rate, estimated P astro , results of signal consistency tests, etc. are included where possible.The data release also contains the posterior samples from our parameter estimate for each significant candidate.The configuration files needed to reproduce both the overall search and individual source parameter estimates are provided.We also provide the data products needed to reproduce the BBH population inference with alternate models; this includes the sensitivity of our search to a reference population.Finally, we make available the posterior samples of the parametric population models used in this study.

CONCLUSIONS
4-OGC is a comprehensive catalog of gravitationalwave observations from BNS, NSBH, and BBH mergers covering the 2015-2020 time period.The 94 mergers included in the 4-OGC catalog represent the largest collection to date.These include 90 BBHs, 2 BNS, and 2 NSBH mergers.7 new BBHs during O3b have been reported here that pass our significance threshold.We estimate the merger rate of BNS, NSBH, and BBH sources at z=0 to be 200 +309 −148 , 19 +30 −14 , and 16 +10 −6 yr −1 , respectively.The release of our catalog can help enable deeper investigations of gravitational wave candidates, for instance, to look for deviations from general relativity, study formation channels of compact-binary mergers, and correlate results with other archival observations.It is expected that future catalogs may contain hundreds of sources after the ongoing upgrades are completed.The next observing run is scheduled to begin operation at the end of 2022 (Abbott et al. 2020b).The importance of gravitational-wave catalogs will only grow as gravitational-wave astronomy matures and focus is able to shift to understanding increasingly detailed population features.
Table 3. Shown are the gravitational-wave observations from the full search on data from O1-O3 with Pastro > 0.5 or IFAR > 100 years sorted by observation time.We list the GPS time for each event alongside information on the observational status of the three observatories LIGO Hanford (H), LIGO Livingston (L), and Virgo (V).We also list the detectors for which our search triggered and the corresponding recovered SNRs (ρ).For some events, where multiple detectors are operational, the search does not trigger for all observatories.This is caused by requiring consistency of the triggers between detectors and triggers in individual detectors needing to exceed the SNR threshold.For multi-detector events we provide the IFAR at the associated ranking statistic for the event.The probability of astrophysical origin is estimated for all events detected by the focused BBH search and from the single-detector analyses.Events reported here for the first time are marked by bold font.3. (Continued) Shown are the gravitational-wave observations from the full search on data from O1-O3 with Pastro > 0.5 or IFAR > 100 years sorted by observation time.We list the GPS time for each event alongside information on the observational status of the three observatories LIGO Hanford (H), LIGO Livingston (L), and Virgo (V).We also list the detectors for which our search triggered and the corresponding recovered SNRs (ρ).For some events, where multiple detectors are operational, the search does not trigger for all observatories.This is caused by requiring consistency of the triggers between detectors and triggers in individual detectors needing to exceed the SNR threshold.For multi-detector events we provide the IFAR at the associated ranking statistic for the event.The probability of astrophysical origin is estimated for all events detected by the focused BBH search and from the single-detector analyses.Events reported here for the first time are marked by bold font.Table 3. (Continued) Shown are the gravitational-wave observations from the full search on data from O1-O3 with Pastro > 0.5 or IFAR > 100 years sorted by observation time.We list the GPS time for each event alongside information on the observational status of the three observatories LIGO Hanford (H), LIGO Livingston (L), and Virgo (V).We also list the detectors for which our search triggered and the corresponding recovered SNRs (ρ).For some events, where multiple detectors are operational, the search does not trigger for all observatories.This is caused by requiring consistency of the triggers between detectors and triggers in individual detectors needing to exceed the SNR threshold.For multi-detector events we provide the IFAR at the associated ranking statistic for the event.The probability of astrophysical origin is estimated for all events detected by the focused BBH search and from the single-detector analyses.Events reported here for the first time are marked by bold font.  5. We report the Bayesian parameter estimates for all 94 detections from O1-O3 scientific runs and mark the ones reported here for the first time in bold.The various columns report median and 90% credible intervals for the source-frame component mass m1 and m2, chirp mass M, mass ratio q, effective spin χ eff , luminosity distance DL, redshift z, and remnant mass and spin M f and χ f , respectively.The SNR is computed from the maximum likelihood for BBH events and the maximum likelihood analytically marginalized over the phase for NSBH or BNS events.The quoted results for BBH, BNS and NSBH events are computed using IMRPhenomXPHM, IMRPhenomD NRTidal, and IMRPhenomNSBH waveform models, respectively.

Event
Event m1/M m2/M M/M q χ eff DL/Mpc z M f /M χ f SNR GW150914 095045 34.9 +4.7 −2.9 29.7 +2.8 Table 5. (Continued) We report the Bayesian parameter estimates for all 94 detections from O1-O3 scientific runs and mark the ones reported here for the first time in bold.The various columns report median and 90% credible intervals for the sourceframe component mass m1 and m2, chirp mass M, mass ratio q, effective spin χ eff , luminosity distance DL, redshift z, and remnant mass and spin M f and χ f , respectively.The SNR is computed from the maximum likelihood for BBH events and the maximum likelihood analytically marginalized over the phase for NSBH or BNS events.The quoted results for BBH, BNS and NSBH events are computed using IMRPhenomXPHM, IMRPhenomD NRTidal, and IMRPhenomNSBH waveform models, respectively.

Figure 1 .
Figure 1.The cumulative number of binary merger observations (top) and the fiducial 1.4-1.4M BNS merger range (bottom) of the LIGO Hanford (yellow), LIGO Livingston (blue), and Virgo (green) observatories at a SNR of 8 as a function of days since the start of the advanced detector era.The distance is averaged over sky location and orientation angles.The O1 (left), O2 (middle) and O3 (right) observing periods are shown.The break in O3 demarks the boundary between O3a and O3b.

Figure 2 .
Figure 2. The bank of templates used to identify compactbinary mergers in our search as a function of their detectorframe masses.The binary neutron star (blue), neutron starblack hole (orange), binary black hole (green), and focused binary black hole (purple) regions are shown.The templates associated with an observed merger are shown with stars.Signals often have multiple templates which will produce a candidate.Here, we only show the template which produced a candidate with the lowest false alarm rate; the parameters of the selected template can only be considered as crude point estimate of the true parameters, so they may differ significantly from the posterior estimates in Table5.

Figure 3 .
Figure 3.The stacked distributions of single-detector triggered candidates observed during O3 when a single LIGO observatory was operating (orange), our selected background (blue), and for comparison the distribution of gravitational-wave mergers observed by the multi-detector analysis (green) as a function of the ranking statistic.Results from the BNS (left), NSBH (right), and BBH (bottom) analyses are shown.The method of Nitz et al. (2020) is used to extrapolate the background distribution and estimate the probability of astrophysical origin.The BBH analysis uses the statistic λs(Nitz et al. 2020)  while the others use a re-weighted SNR statistic(Babak et al. 2013;Nitz et al. 2019c).GW190425 appears in both the BNS and NSBH analysis as candidates were produced from both regions.

Figure 4 .
Figure4.The marginalized distributions for component masses m1, m2, the effective spin χ eff and the luminosity distance DL for all events that pass our detection criteria in 4-OGC.The 5th and 95th quantile values are marked with a bar, respectively.Different colors are used to aid associating each event with its posterior estimates.

Figure 5 .
Figure 5. Distribution of the source-frame component masses of the BBH population obtained from combining posteriors from all the confident BBH events in the 4-OGC analysis.We show the observed component mass distribution (left), component mass distribution corrected for the basic selection effect (middle), and the one-dimensional marginals of the total mass distribution: observed as well as corrected for selection effects (right).To correct for selection effects, we assume a constant detection threshold and correct the distribution for signal loudness varying with component masses.tion (Clesse & Garcia-Bellido 2020; Olejak et al. 2020; Lu et al. 2020; Carr et al. 2021; Liu & Lai 2021).The observed distribution of source parameters (e.g.masses and spins) and redshift distribution from a large catalog of compact-binary mergers can be used to constrain the population models that predict intrinsic source distribution as well as the rate of mergers(Abbott et al. 2019bRoulet et al. 2021;Zhu et al. 2021a).In figure5, we show the distribution of component masses obtained from stacking all the posteriors marginalized over all other parameters, with and without taking into account the zeroth order selection effect that arises due to the loudness of the signal as a function of source parameters.To infer the source population, we use the same method as in the 3-OGC analysisNitz & Capano (2021).We combine the posterior samples from each event to obtain large mass samples.We assume the same priors that we used in the parameter estimation.To account for the loudness of the signal, we estimate the comoving volume corresponding to the "horizon distance" of a posterior sample and reassign each sample a weight proportional to the inverse of the comoving volume.The horizon distance is defined as the maximum distance an optimally oriented source can be detected.
O1+O2 and with direct search performance for O3(Abbott et al. 2021f; Collaboration et al. 2021a,b) using a fiducial BBH population set with component masses between 2M −100M .We make available the data products needed to perform more detailed follow-up analyses(Nitz & Kumar 2021); astrophysically-motivated models may provide further insights(Fragione & Banerjee 2021; Stevenson et al. 2019).

Figure 6 .Figure 7 .
Figure6.Differential merger rate as a function of source-frame primary mass for the parametric mass models considered in the study: power law (orange), powerlaw + peak (blue), and powerlaw + two peaks (purple).We marginalize over all other population parameters in each model.Solid lines show the median value for the analysis done with 4-OGC events.Shaded regions show the 90% credible intervals except; the powerlaw + two peaks model is omitted to improve visual clarity.For comparison, the powerlaw + peak result from GWTC-3(Abbott et al. 2021f) is shown with a black dashed line and grey shaded region.

Table 4 .
The selection of sub-threshold candidates with Pastro > 0.2 or IFAR > 0.5 from the full search of O1-O3 data.Candidates are sorted by the observation time.The complete set of sub-threshold candidate is available in the data release and includes a selection of full parameter estimates.Here we show the detector-frame (redshifted) parameters of the template which triggered on the candidate, along with the reported SNRs (ρ) from each detector.

Table 4 .
(Continued)The selection of sub-threshold candidates with Pastro > 0.2 or IFAR > 0.5 from the full search of O1-O3 data.Candidates are sorted by the observation time.The complete set of sub-threshold candidate is available in the data release and includes a selection of full parameter estimates.Here we show the detector-frame (redshifted) parameters of the template which triggered on the candidate, along with the reported SNRs (ρ) from each detector.Table