Subtracting compact binary foreground sources to reveal primordial gravitational-wave backgrounds

Detection of primordial gravitational-wave backgrounds generated during the early universe phase transitions is a key science goal for future ground-based detectors. The rate of compact binary mergers is so large that their cosmological population produces a confusion background that could masquerade the detection of potential primordial stochastic backgrounds. In this paper we study the ability of current and future detectors to resolve the confusion background to reveal interesting primordial backgrounds. The current detector network of LIGO and Virgo and the upcoming KAGRA and LIGO-India will not be able to resolve the cosmological compact binary source population and its sensitivity to stochastic background will be limited by the confusion background of these sources. We find that a network of three (and five) third generation (3G) detectors of Cosmic Explorer and Einstein Telescope will resolve the confusion background produced by binary black holes leaving only about 0.013\% (respectively, 0.00075\%) unresolved; in contrast, as many as 25\% (respectively, 7.7\%) of binary neutron star sources remain unresolved. Consequently, the binary black hole population will likely not limit observation of primordial backgrounds but the binary neutron star population will limit the sensitivity of 3G detectors to $\Omega_{\rm GW} \sim 10^{-11}$ at 10 Hz (respectively, $\Omega_{\rm GW} \sim 3\times 10^{-12}$).


I. INTRODUCTION
With the continued detections of gravitational waves from binary black hole mergers [1][2][3][4][5][6] and binary neutron star inspirals [7,8], the LIGO Scientific and Virgo Collaborations have kept up to their promise of taking us into an era of gravitational-wave astronomy. In addition to these loud and nearby sources that are seen as isolated transient events, there is a population of weak, unresolved sources at higher redshifts [9][10][11][12][13]. The superposition of these sources is expected to be the main contributor to the astrophysical stochastic background which may be detectable in the next few years as the Advanced LIGO [14] and Virgo detectors [15] reach their design sensitivity and accumulate more data [16,17]. Assuming the most probable rate for compact binary mergers at the time (103 +110 −63 Gpc −3 yr −1 [3] for BBH and 1540 +3200 −1220 Gpc −3 yr −1 [7] for BNS), it has been shown that the total background may be detectable with a signalto-noise-ratio of 3 after 40 months of total observation time, based on the expected timeline for Advanced LIGO and Virgo to reach their design sensitivity [17]. The astrophysical background potentially contains a wealth of information about the history and evolution of a population of point sources, but it is a confusion noise background that obscures the observation of the primordial gravitational-wave background (PGWB) produced in the very early stages of the Universe. Proposed theoretical cosmological models include the amplification of vacuum a szs1416@psu.edu b tania.regimbau@lapp.in2p3.fr c bss25@psu.edu fluctuations during inflation [18][19][20], pre Big-Bang models [21][22][23], cosmic (super) strings [24][25][26][27] , or phase transitions [28][29][30]. For a comprehensive discussion of cosmological gravitational-wave backgrounds, we refer the reader to reviews by Maggiore and [31] and Binétruy et al. [32].
Detection of the primordial gravitational-wave background would create a unique window on the very first instants of the Universe, up to the limits of the Planck era, and on the physical laws that apply at the highest energy scales. Needless to say that such a detection would have a profound impact on our understanding of the evolution of the Universe.
The current detector network of LIGO and Virgo and the upcoming KAGRA and LIGO-India will not be able to resolve the cosmological compact binary source population and its sensitivity to stochastic background will be limited by the confusion background of these sources [48]. With the increased sensitivity of the third generation gravitational-wave detectors, such as the Einstein Telescope (ET) [49] and the Cosmic Explorer (CE) [50], it will be possible to detect and resolve almost all of the binary black hole mergers, even the ones at high redshifts. In this work, we explore the possibility of probing the cosmological gravitational-wave background with the third generation detectors, after removing the astro-physical background from compact binary mergers from the data. This work is an extension to [48], where the authors have shown the level at which we can expect amplitude of background from unresolved, subthreshold signals from compact binary coalescences (CBC) using different detector networks. We extend the previous study to also provide an estimate of errors we introduce while subtracting the signals above threshold for the most optimistic network of detectors considered by [48]. The idea of subtracting foreground signals to extract stochastic backgrounds was already explored [51] in the context of the the Big Bang Observer [52], including a noise projection method that could reduce errors due to imperfect subtraction [53].
Data from gravitational-wave detectors are dominated by environmental and instrumental backgrounds. Consequently, it is not possible to identify even deterministic signals without sophisticated data processing such as matched filtering [54]. Stochastic backgrounds cannot be reliably detected in a single detector-they are found by cross-correlating the data from a pair of detectors. Indeed, the stochastic background present in one of the detectors acts as a matched filter for the data in the other detector [55][56][57]. Unfortunately, this means that any common noise in a pair of detectors could masquarade as stochastic background [58]. If detectors are geographically well separated then the risk of common noise of terrestrial origin is greatly reduced. Additionally, certain backgrounds of terrestrial origin could be measured and subtracted [59]. Even in the absence of any terrestrial background, a pair of detectors would see the same astrophysical background, which would show up as correlated 'noise' although detectors might be geographically well separated. As a result, the only possible way to improve the sensitivity of a detector network to primordial backgrounds is to subtract foreground astrophysical signals.
The rest of the paper is organized as follows. In Sec. II, we describe the basic method that we use to calculate the gravitational-wave spectrum from the error introduced by imperfect subtraction of CBC signals. In Sec. III, we describe the framework used to estimate the deviations of the estimated parameters of the CBC sources from their true values. We discuss the simulation of a population of binaries in Sec. IV, discuss the result of the imperfect subtraction of such signals in Sec. V, and we discuss our results in Sec. VI.

II. METHOD
The energy-density spectrum in gravitational waves is described by the dimensionless quantity [57], where dρ GW is the energy density in the frequency interval f to f + df , ρ c = 3H 2 0 c 2 /8πG is the closure energy density, and H 0 is the Hubble constant equal to 67.8±0.9 km/c/Mpc [60]. The gravitational-wave energy spectrum density can be written as a sum of contribution from the astrophysical and cosmological energy densities, Taking the contribution of the compact binary coalescences out of the astrophysical background, and writing it explicitly, we have, Here Ω astro, r is the remaining astrophysical background after taking out the contribution from the CBC sources. When estimating the parameters of a binary source, by using Monte Carlo methods, or nested sampling, we invariably end up with parameters that deviate from the true values because of the noise in the detector. Therefore when we subtract the recovered CBC signals from the data, we introduce an additional background due to the error in subtraction, Ω error .
where Ω cbc, rec is the background from the recovered CBC sources that we can subtract from our data, Ω error is the background because of the error introduced from such a subtraction, Ω cbc, unres is the background from the unresolved CBC sources which are not detected as foreground events. Let us assume that we have an experiment where we have detected a list of CBC sources and subtracted them from the data. Now we are left with the gravitational-wave backgrounds, Ω error , Ω cbc, unres , on top of the cosmological and astrophysical (from sources other than the CBCs) backgrounds. We want to answer the question of whether the cosmological or astrophysical backgrounds from sources other than CBCs can stand above the residual background after removal of the CBC sources. That is, In order for us to be able to detect the gravitationalwave background from cosmological sources or that from different astrophysical sources, we would need Ω residual = Ω error + Ω cbc, unres to lie below these. The gravitational-wave energy density from a population of compact binary sources is given by [48], where F (f ) is the total flux, sum of individual contributions where N is the number of sources in the Monte Carlo sample, and T −1 assures that flux has the correct dimension, T being the total time of the data sample. h +,k (f ) andh ×,k (f ) are the Fourier domain waveforms for the two polarizations, and the index k runs over all the sources. We calculate Ω error as, where, To get an estimate of Ω error , we need to estimate the quantities,h recovered

III. ESTIMATING THE DEVIATION FROM TRUE VALUE OF THE MEASURED SOURCE PARAMETERS
Ideally we want the full Bayesian posteriors to estimate the deviation from the true value of parameters. However, at present it is unfeasible to compute the full posterior probability distribution functions of all 15 binary parameters for the hundreds of thousands of sources that we simulate up to a redshit of 10 in the following section. The Fisher matrix provides a computationally inexpensive method to estimate the errors in the case when the posteriors are Gaussian, which is, unfortunately, not true in general. Nevertheless, for the purpose of building a proof-of-principle concept the Fisher matrix method is adequate and the only practical approach to obtain the magnitude of errors in the estimation of parameters. To this end, we follow the framework described in [61] and calculate the errors in estmating the parameters of the compact binary system using the Fisher matrix method.
According to the post-Newtonian expansion formalism [62], the gravitational-wave strain from a compact binary coalescence in frequency domain is given bỹ where A is the amplitude of the waveform, and Ψ(f ) is the phase given by Here t c is the time of coalescence, φ c is the coalescence phase, ν = (πM f ) 1/3 , M is the total mass (M = m 1 +m 2 ), η is the symmetric mass ratio (η = m 1 m 2 /M 2 ) of the system, and the α k terms are known as the post-Newtonian (PN) coefficients. In this work, we restrict ourselves to 0-PN approximation (or the Newtonian approximation, k = 0), which will be justified below. For the Fisher matrix study, we choose a set of independent parameters θ for describing the gravitational waveform, where f 0 is a reference frequency needed to keep the parameters for the Fisher matrix dimensionless. M is the dimensionless chirp mass, and is defined as M = η 3/5 M/M . Writing the phase of the waveform in terms of these parameters, we have, or equivalently, In going from Eq. (13) to Eq. (14), we have truncated the expansion at α 0 term, plugged in the value α 0 = 1, and we have introduced the Newton's constant G, the speed of light c, and solar mass M , explicitly to keep all quantities in the Eq. (13) dimensionless, and defined masses in solar mass units.
The Fisher matrix elements are given by, are the partial derivatives of the waveform with respect to θ i , the parameters of the waveforms, and S n (f ) is the single-sided power spectral density of the detector. The partial derivatives of the waveform can be calculated analytically: and, h θ3 (f ; θ) = Af −7/6 e i(Ψ(f ; θ)−π/2) 5 128 The variance-covariance matrix, or simply the covariance matrix, defined as the inverse of the Fisher information matrix, is given by Once we have the covariance matrix, we use a multivariate normal random number generator to generate observed values of the parameters, P O , based on the multivariate dirstribution with the mean equal to the true value of the parameters, P T and covariance matrix as Σ. The error in parameter estimation is then given by where

IV. POPULATION SYNTHESIS FOR MULTIPLE DETECTORS
We simulate a population of binary black hole and binary neutron star systems up to a redshift of 10, and then calculate an estimate of Ω cbc, rec and Ω error as outlined in Sec. II and Sec. III. The list of compact binaries (neutron star binaries or black hole binaries) is generated following a Monte Carlo procedure described in [48,[63][64][65], and using the fiducial model of [17] for the distribution of the parameters (masses, redshift, position on the sky, polarization and inclination angle of the binary). In particular, we assume a redshift distribution which is derived from the star formation rate (SFR) of [66] and accounts for a delay between the formation of the progenitors and the merger. We further consider the median rates estimated from the first LIGO observation run.
1. For BBHs, the intrinsic masses m 1 , m 2 (in the source frame) are selected from the power-law distribution (Saltpeter initial mass function [67]) considered in [3,68] of the primary (i.e., the larger mass) companion p(m 1 ) ∝ m −2.35 1 and from a uniform distribution of the secondary companion. In addition, we require that the component masses take values in the range 5-50 M . For BNSs, the intrinsic masses m 1 , m 2 (in the source frame) are both drawn from a Gaussian distribution centred around 1.33 M with a standard deviation of 0.09 M .
2. The redshift z is drawn from a probability distribution p(z) given by obtained by normalizing the merger rate of binaries in the observer frame, R z (z) per interval of redshift, over the range z ∈ [0, 10]. We choose to cut off the redshift integral at z max = 10, since redshifts larger than 5 contribute little to the background [17]. The merger rate in the observer frame is 1 where dV /dz is the comoving volume element and R m (z) is the rate per comoving volume in the source frame, given by where R f (z f ) is the binary formation rate as a function of the redshift at formation time, z f = z(t f ) is the source redshift at formation, p(t d ) is the distribution of the time delay t d between the formation and merger of the binary, z = z(t m ) is the source redshift at merger. The integration in Eq. 26 over z f is performed for all the redshifts corresponding to t f such that t m = t f + t d .
We assume that the binary formation rate R f (z f ) scales with the SFR. We follow the the cosmic star formation model of [66] which uses the Springer-Hernquist functional form [78] to fit to the GRB-based high-redshift SFR data of [79] but normalized based on the procedure described in [80,81]. This fit results in ν = 0.146M /yr/Mpc 3 , z m = 1.72, a = 2.80, and b = 2.46 [66]. The value of R m (z = 0) is chosen as the local merger rate estimate from the LIGO-Virgo observations. For the rate of BBH mergers, we use the most recent published result associated with the power-law mass distribution 56 +44 −27 Gpc −3 yr −1 [6]. For the BNS case, we set R m (z = 0) to 920 +2220 −790 Gpc −3 yr −1 also from [6]. Massive black holes are formed preferentially in low-metallicity environments [16,82]. For systems where at least one black hole has a mass larger than 30M , we re-weight the star formation rate R f (z) by the fraction of stars with metallicities less than half the solar metallicity [17]. Following [16,17], we use the mean metallicity-redshift relation of [83], and scale it upwards by a factor of 3 to account for local observations [66,84].
3. The location on the sky, the cosine of the inclination angle, the polarization, and the coalescence phase are drawn from uniform distributions.

A. Detector Network
We consider two networks of third generation detectors: one with three total detectors, out of which two have the sensitivity of CE located at LIGO Hanford and LIGO Livingston locations and one with the sensitivity of ET located at the location of Virgo; and a five-detector network with one detector with the sensitivity of ET at the location of Virgo, and detectors with CE sensitivity at locations of LIGO Hanford, LIGO Livingston, LIGO India, and KAGRA. We choose these configurations for the detector-networks because it was shown in [48], that the astrophysical "confusion" background from unresolved BBH sources is decreased by orders of magnitude, reaching Ω GW (10Hz) = 10 −14 − 10 −13 and Ω GW (10Hz) = 10 −16 − 10 −14 respectively.

V. SIMULATIONS
We simulate a population of BBH and BNS mergers according to the procedure described in Sec. IV for a year of data. There are 76,107 BBH and 1,438,835 BNS signals in our simulation. For each source, we calculate the expected network SNR assuming perfect template match, given by where index i runs over all the sources, and ρ det is the SNR for each source and detector pair (i, det), andh det i (f ) = F det +hi,+ + F det ×hi,× is the Fourier domain waveform projected on the detector.
We considered a source as resolvable and a part of the "foreground", whenever ρ net i ≥ ρ thresh = 12.0. We use the 0 order PN approximation for waveforms, since the results from that and a full inspiral-merger-ringdown model agree to a great extent below 100 Hz. It has been shown for various detector combinations that frequencies below 100 Hz account for more than 99% of the SNR for the stochastic search [64,69]. Therefore for calculating Ω error , we only consider the 0th-PN model to compute the Fisher matrix for each source in our simulation.
We calculate the Fisher matrices (and the variancecovariance matrices) for all the sources in our simulation, and recover a set of parameters in order to calculate Ω residual,BNS and Ω residual,BBH .
Our results are shown plotted in Fig. 1. For the threedetector case, we find that 49% of the BNS sources are unresolved (with a network SNR < 12), whereas only 0.013% of the BBH sources are unresolved. For the fivedetector case, we find that 25% of the BNS sources are unresolved while only 0.00075% of the BBH sources remain unresolved. We show the results for network SNR threshold of 12 in the first two rows of Fig. 1. The first row shows the results for BBH (left: for a 3 detector 3G network, right: for a 5 detector 3G network) and the second row shows the results for the BNS (left: for a 3 detector 3G network, right: for a 5 detector 3G network). We can see that the Ω residual = Ω error + Ω cbc, unres depends on the network SNR threshold. The higher the network SNR threshold, the lower the Ω error but higher the Ω cbc, unres . Thus, the network SNR threshold can be varied to minimize the Ω residual .
For the BBH case, we have not tried to optimize the Ω residual, BBH , since it lies much below the Ω residual, BNS . For the BNS case, we can see from the second row of Fig. 1, that we may be able to lower the residual background by decreasing the network SNR threshold, since the residual is dominated by the unresolved sources. We decided to lower the network SNR threshold to 8 (the threshold at which we should be able to resolve signals in case of Gaussian noise); these results are shown in the last row of Fig. 1. With a network SNR threshold of 8, the number of unresolved BNS sources for a three (and five) network of 3G detectors reduces to 25% from 49% (7.7% from 25%). We have managed to lower the BNS residual background by lowering the detector network SNR threshold. The residual background from the BNS sources still dominates over the BBH background and is the limiting factor for the primordial backgrounds we can observe. An alternative would be to follow the noise projection method described in Ref. [51], which does not require the SNR optimization procedure described here.

VI. DISCUSSION
Conclusions of our study are summarized in Fig. 2. The figure plots the energy density in gravitational waves Ω(f ) from axion inflation [85], a network of cosmic strings [24][25][26][27], a background produced during post-inflation by oscillations of a fluid with an equation-of-state stiffer than radiation [86], and from post-inflation preheating scenarios [87,88] aided by parametric resonance [23,89]. For reference, we show the strength of the stochastic background from vacuum fluctuations during standard inflation [18][19][20], although this will not be detectable by any of the foreseen ground-based detector networks; oth-      , the background that remains after imperfect subtraction of resolved sources (dashed, red lines) and the sum of the latter two (solid, deep-red lines). The left panels are for a network of three 3G detectors and the right panels are for a network of five 3G detectors. We deem a source is resolved if the signal-to-noise it produces is ≥ 12 for the top four panels, and ≥ 8 for the bottom two panels.  The figure also shows the sensitivity of a network of three (and five) 3G detectors to stochastic backgrounds assuming a one-year integration but in the absence of confusion backgrounds from compact binaries or other astrophysical populations. It is immediately apparent that the residual background, after (imperfect) subtraction of the foreground sources, from binary neutron stars will limit the strength of primordial backgrounds that could be detected by 3G detectors. With a network of three (and five) 3G detectors, the sensitivity will be limited to Ω GW ≥ 10 −11 at 10 Hz (respectively, Ω GW ≥ 3 × 10 −12 at 15 Hz). The binary black hole population, on the other hand, can be fully resolved and the residual from that population has negligible effect on the raw sensitivity to stochastic backgrounds. The rate of binary neutron stars could be larger or smaller than the median rate of R m (z = 0) = 920 +2220 −790 Gpc −3 yr −1 assumed in this paper, which would correspondingly increase or decrease the confusion background of these sources. Finally, increasing the number of 3G detectors from three to five improves the sensitivity to stochastic backgrounds by about factor of 5. This is accounted by the ability of the five-detector network to detect and subtract a greater number of sources; the volume reach for a five-detector network increases by a factor (5/3) 3 ∼ 4.6 relative to a three-detector network.
Keeping in mind that the strengths of the primordial backgrounds depend on the specific model parameters that are not known, and the residual background could vary based on the uncertainty in rate of compact binary mergers and the their mass distribution, among other things, the figure shows the most promising primordial background sources that this subtraction scheme could reveal: cosmic strings, background from fluids with stiff EOS, and axion inflation.