Non-Stationary Astrophysical Stochastic Gravitational-Wave Background: A New Probe to the High Redshift Population of Binary Black Holes

The astrophysical Stochastic Gravitational Wave Background (SGWB) originates from the mergers of compact binary objects that are otherwise undetected as individual events, along with other sources such as supernovae, magnetars, etc. The individual GW signal is time-varying over a time scale that depends on the chirp mass of the coalescing binaries. Another timescale that plays a role is the timescale at which the sources repeat, which depends on the merger rate. The combined effect of these two leads to a breakdown of the time-translation symmetry of the observed SGWB and a correlation between different frequency modes in the signal covariance matrix of the SGWB. Using an ensemble of SGWB due to binary black hole coalescence, calculated using simulations of different black hole mass distributions and merger rates, we show how the structure of the signal covariance matrix varies. This structure in the signal covariance matrix brings additional information about the sources on top of the power spectrum. We show that there is a significant improvement in the Figure of Merit by using this additional information in comparison to only power spectrum estimation for the LIGO-Virgo-KAGRA (LVK) network of detectors with the design sensitivity noise with two years of observation. The inclusion of the off-diagonal correlation in the covariance of the SGWB in the data analysis pipelines will be beneficial in the quest for the SGWB signal in LVK frequency bands as well as in lower frequencies and in getting an insight into its origin.


INTRODUCTION
Gravitational waves (GW) are a novel cosmic messenger that can help us answer some of the most important questions in astrophysics, cosmology, and fundamental physics (Thorne 1995;Sathyaprakash & Schutz 2009;Bailes et al. 2021;Perkins et al. 2021;Berti et al. 2022;Abdalla et al. 2022;Auclair et al. 2022;Mastrogiovanni et al. 2022;Adhikari et al. 2022).The first direct detection of gravitational waves (GW150914) was made by the LIGO-Virgo Collaboration (Abbott et al. 2016) in September 2015 which came from the merger of two black holes (BHs).So far, the LVK has detected close to 90 compact binary mergers, which include binary black holes (BBHs), binary neutron stars (BNSs), and neutron star-black holes (NSBHs) mergers (Abbott et al. 2021c).Most of these sources are located below the redshift of z = 1 for fiducial Planck cosmology (Abbott et al. 2021c).The present detectors can only detect sources individually that are located below the redshift of one with a matched filtering SNR above 8 (Abbott et al. 2019).All the unresolved events appear as the stochastic gravitational wave background (SGWB) signal (Apreda et al. 2001;Zhu et al. 2011;Romano & Cornish 2017).The SGWB can be detected by ground-based detectors using a technique called cross-correlation which is the conventional approach to detect the SGWB signal (Allen & Romano 1999;Thrane & Romano 2013;Thrane et al. 2009).Other search techniques include The Bayesian Search (TBS) (Smith & Thrane 2018) and Cross-Correlation Intermediate (CCI) search (Coyne et al. 2016).The Pulsar Timing Array (PTA) search (Burke-Spolaor et al. 2019;Verbiest et al. 2022;Manchester 2013), operating in the nano-Hertz range, relies on the correlated signatures in the pulse arrival times of a set of pulsars.Recently, the PTA collaborations have reported the detection of the nano-Hertz gravitational waves (Zic et al. 2023;Agazie et al. 2023;Antoniadis et al. 2023;Lee 2023).However, in the LVK, we have only been able to place an upper limit on the SGWB (Abbott et al. 2021a,b).Nevertheless, we expect to detect the SGWB signal with future upgraded detectors (Thrane & Romano 2013;Christensen 2018;Renzini et al. 2022;Mentasti et al. 2023;Suresh et al. 2021).The astrophysical SGWB is a combination of various types of sources such as compact binary mergers, supernovae, magnetars, etc (Buonanno et al. 2005;Chowdhury & Khlopov 2021).However, in this paper, we concentrate solely on examining the SGWB resulting from the coalescence of BBH detected by ground-based detectors.This work can also be extended to other types of sources as well as detectors.
The SGWB is a probe to the high redshift universe using which we can study the properties of compact objects and explore their formation channels (Mukherjee & Silk 2021b;Bavera et al. 2022;Lehoucq et al. 2023;Babak et al. 2023).One advantage of SGWB over individual source detection is that SGWB can explore deep into high redshift for sources over a large range of compact object masses ranging from sub-solar mass to super-solar mass, whose origin can be primordial or astrophysical.The merger rate of compact objects of astrophysical origin, as well as their mass and spin distributions, depends on the stellar properties of galaxies (star formation rate, stellar mass, stellar metallicity, etc.) and also the formation channels (Belczynski et al. 2002;Renzo et al. 2020;Bethe & Brown 1998;Mukherjee 2022;Spera et al. 2019;Dorozsmai & Toonen 2022;Kruckow et al. 2018;Srinivasan et al. 2023).As a result, the SGWB signal and its connection with the stellar properties of the galaxies can provide rich information about the high redshift Universe in a way complementary to the electromagnetic probes.SGWB can also help us distinguish different formation channels of BHs.It is a powerful probe to distinguish astrophysical black holes (ABHs) from primordial black holes (PBHs) (Mukherjee & Silk 2021b).The merger rate of ABHs is expected to follow the star formation history.In contrast, PBHs, being formed in the very early universe, are expected to exhibit significantly different merger rates.The merger rate of PBHs is going to dominate over ABHs at high redshifts.This will appear as a distinguishable signature on the SGWB spectrum (Mukherjee & Silk 2021b;Mukherjee et al. 2022a;Atal et al. 2022).
The number of mergers contributing to the SGWB in a given time interval is expected to follow the Poisson distribution (Mandel & O'Shaughnessy 2010;Bulik et al. 2011;Dvorkin et al. 2018).This leads to fluctuations in the number of events with time (Mukherjee & Silk 2020, 2021a;Braglia et al. 2023;Ginat et al. 2023Ginat et al. , 2020)).The amount of fluctuation will depend on the merger rate and observation time.Along with that, the mass distribution of the compact objects contributing to the background can also contribute to the temporal fluctuations in the SGWB.Recently, a new technique has been proposed to search for this signal from the background using spectrogram analysis (Dey et al. 2023).Furthermore, there are other proposed techniques aimed at distinguishing these signals based on their statistical properties (Smith & Thrane 2018;Buscicchio et al. 2023;Lawrence et al. 2023).
The measurement of the SGWB signal is typically performed through the cross-correlation of the GW strain between two detectors.Such an estimator is optimal for types of signal which is continuous and Gaussian.However, the SGWB signal of astrophysical origin is expected to exhibit temporal dependence and break stationarity due to two effects: (i) the time dependence of the individual signal, which relies on the masses of the GW sources, and (ii) the time scale over which the signal repeats.Despite the temporal variations, the SGWB signal will appear statistically timetranslation symmetric over a large timescale, as the astrophysical source population properties contributing to the signal will not change over the time scale over which observations are made.
The presence of a non-stationary SGWB signal leads to a correlation between two different frequency modes of the signal, and as a result, the signal covariance matrix no longer remains diagonal but contains non-zero off-diagonal terms.In Fig. 1 we illustrate how the non-stationary SGWB signal can give rise to correlations between different frequencies in the covariance matrix and how this spectral covariance can enhance our estimates of the parameters.The Poissonian nature of the BH merger and the chirp signal breaks the time translation symmetry of the signal.However, the noise remains stationary and therefore retains time translational symmetry1 .As a result, the spectral covariance matrix of the signal contains non-zero diagonal elements, while the noise spectral covariance matrix remains diagonal.Therefore, the off-diagonal term in the spectral covariance matrix remains uncontaminated from noise.This allows the spectral covariance to provide additional information to put a better constraint on the GW source parameters.
Using a simulation-based approach to estimate SGWB for different populations, developed in this work, we estimate the SGWB signal and show that there is an intrinsic correlation between different frequencies, which can vary depending on the ABH and PBH population and their merger rates.Our results demonstrate that the statistical distribution of the SGWB signal can serve as a robust probe of high-redshift stellar properties, which can be explored within the detectable frequency band of the LIGO-Virgo-KAGRA (LVK).The paper is organized as follows: in Sec. 2, we discuss the simulation technique developed for different models of the ABHs and PBHs; in Sec. 3, we discuss the summary statistics of the SGWB and show the existence of the non-zero off-diagonal terms in the covariance matrix of the SGWB; in Sec. 4, we carry out Fisher analysis to obtain the expected constraints on the population parameters from the additional information that can be gained using the offdiagonal terms in the covariance matrix of the SGWB signal; finally, in Sec. 5, we discuss the conclusion and the future prospects.

Simulation of the SGWB
We have developed a method to simulate the SGWB signal for different black hole population models.In this section, we briefly describe the method to simulate the SGWB signal.SGWB due to BBH coalescences is a result of the superposition of individual merger events of all BH properties and formation channels, originating from low to very high redshift.The SGWB can be simulated by simulating individual events of all source masses up to very high redshift, following the given merger rate and mass distribution.In order to simulate the SGWB, we calculate the merger rate as well as mass distribution as a function of redshift.We then divide the redshifts into multiple bins, extending up to a very high redshift.For each bin, we employ Poisson sampling to determine the number of events occurring within a given observation time.For each such event, we sample two component masses from our mass distribution.Finally, we sum the background density due to all events in the given observation time.The simulation is moderately computationally expensive, requiring approximately 90 core hours to generate 1 year of SGWB data with a short-time bin of 200 seconds.

Population models for astrophysical black holes
The ABHs are expected to follow the star formation history.Due to the non-zero time delay between the formation of stars and the merger of the BHs, the merger rate is expected to shift towards lower redshifts compared to the star formation rate (SFR).Here, we model the SFR using the Madau-Dickinson star formation rate (Madau & Dickinson 2014).

Merger Rate of Astrophysical Black Holes
The merger rate of ABHs of masses m1 and m2 at redshift z can be written as (Mukherjee & Silk (2021b); Karathanasis et al. (2022a)) where PA(m, z) is the mass distribution of the ABHs merging at a redshift of z and RABH(z) is the source frame merger rate of ABHs per unit comoving volume (Karathanasis et al. 2022b) (2) where R0 is the merger rate at z=0, RSFR(z) is the Madau-Dickinson star formation rate (Madau & Dickinson 2014), 1 + ( (1+z) 2.9 )  Cao et al. 2018;Elbert et al. 2018), which we model as For a uniform log-space distribution of the initial separation between the binaries, κ = 1 (Beniamini & Piran 2019;Vitale et al. 2019;Cao et al. 2018).The current bound on the value of κ from GWTC-3 is κ > 0.7 (Karathanasis et al. 2022b).
In Fig. 2(a), we show the merger rate of ABHs with different minimum time-delay (tmin) for local merger rate of R0= 30 Gpc −3 yr −1 .For large tmin, the peak of the curve shifts towards the lower redshift.This is because the star formed at some redshift will take a longer time to merge for larger tmin.

Mass Distribution of Astrophysical Black Hole
The upper limit on the mass of the BH is set by the pairinstability supernova (PISNe) and is known as the PISN limit (around 40 M⊙ to 50 M⊙) (Farmer et al. 2019;Heger & Woosley 2002;Rakavy & Shaviv 1967;Fraley 1968;Kasen et al. 2011).Pair-instability supernova is a type of supernova that occurs in stars with initial masses 140 M⊙ ≲ M ≲ 260 M⊙ (final helium core masses of 60 M⊙ ≲ M ≲ 140 M⊙) due to pair production (electron-positron pair) in the collision between nuclei and gamma rays produced in the core (Farmer et al. 2019;Kasen et al. 2011;Fryer et al. 2001).PISNe can completely disrupt the star, leaving no stellar remnant behind (Gilmer et al. 2017;Farmer et al. 2019;Heger et al. 2003).
The PISN mass limit of BHs has been shown, through simulation, to vary with the stellar metallicity (Farmer et al. 2019).This is due to the mass lost via winds before the star starts pulsating.Stars with high metallicity are expected to lose more mass via wind as compared to low metallicity stars, thus leaving behind a lighter remnant (Mokiem et al. 2007;van Loon 2005;Vink et al. 2001).This means that the PISN mass limit is larger for stars with low metallicity.Since the metallicity is known to decrease with redshift, the PISN limit (MPISN) is expected to increase with redshift, making mass distribution a redshift-dependent quantity.Therefore, careful modeling of the dependence of MPISN on metallicity can help us understand the metallicity evolution of the universe along with the PISN mass scale itself.We model the metallicity dependence of MPISN as linear in the log of Metallicity (Z) based on the simulations by Farmer et al. (2019).
where α is a parameter that defines the metallicity dependence of MPISN and Z0 is the metallicity at z = 0. Similarly, we model redshift (z) evolution of metallicity (Z) as (Mukherjee 2022;Madau & Dickinson 2014;Karathanasis et al. 2022b,a) log where γ is a free parameter.Combining this with Eq. ( 5) gives Due to the time delay between the formation and merger of BHs and redshift-dependent mass distribution, there is going to be a mixing of BHs from different redshifts.The origin of the component BHs can be significantly differ- ent which makes the observed mass distribution of BHs very different from the mass distribution of BHs formed at that redshift.We model the mass distribution of the secondary mass of BBHs merging at redshift z by a power law, and the primary mass by a power law with a bump near MPISN, which is in agreement with the LVK GWTC-3 observation (Abbott et al. (2021c)).The bump can arise as a result of the accumulation of the BHs formed from the star undergoing a Pulsation pair-instability supernova (PPISNe) (Woosley 2017;Farmer et al. 2019).The observed mass distribution of the BHs is given by (Mukherjee 2022; Karathanasis et al. 2022c) where W(z) is the window function that takes into account the time-delay distribution.
where Ws(z) is the window function of BH mass at redshift z, and Ps(z, m) is given by (Mukherjee 2022;Karathanasis et al. 2022c) where N (m|MPISN, σm) is a Gaussian distribution with a mean value MPISN and a standard deviation σm.Pow(m|a) is a power-law distribution with an index a.The parameter λ controls the height of the bump.In Fig. 2(b), we show the mass distribution of BHs merging at three different redshifts for tmin = 300 Myr.The mass distribution shifts towards higher mass at higher redshift.

Population models of primordial black holes
PBHs are hypothesized to have formed in the very early universe soon after the Big Bang (Carr et al. 2021;Carr 1975;Niemeyer & Jedamzik 1999;Raidal et al. 2017).The merger rate of PBHs is expected to increase with redshift (Mukherjee et al. 2022b;Mukherjee & Silk 2021b;Ng et al. 2022).We model the merger rate of the PBHs as a power law in redshift (Mukherjee & Silk 2021b;Mukherjee et al. 2022a).The merger rate along with the mass distribution of the PBHs can shed light on the properties and characteristics of dark matter (Bird et al. 2022;Carr et al. 2016;Belotsky et al. 2014).It is believed that the PBHs may constitute a notable portion of the dark matter, and their population and mass distribution can help us constrain the properties of dark matter.we characterize the mass distribution of the PBHs as log-normal distribution (Dolgov & Silk 1993;Carr et al. 2017).
where Mc is the characteristic mass scale and σp is the standard deviation of the log(m/Mc).In Table .1, we list the fiducial values of the population parameters considered in the paper.

SUMMARY STATISTICS OF THE SGWB FROM SIMULATIONS
We use simulation-based techniques to estimate the power spectrum of the SGWB, as well as the non-stationary behaviour of the SGWB.We show below the results from our simulations of the (i) SGWB power spectrum and (ii) the signal covariance matrix of the SGWB, and the effect of non-stationary signals on the correlation between different frequency modes.

SGWB power spectrum
The SGWB is defined as the energy density per logarithmic frequency interval divided by the critical energy density (ρc where RGW(z, m1, m2) is the source frame merger rate of BHs per unit comoving volume between masses m1 and m2 at redshift z, dL is the luminosity distance of the source, and fr = f (1+z) is the source frame frequency.We take Mmin= 5 M⊙, and the value of zmin is selected as the maximum redshift at which the A+ detector is sensitive to the most probable mass of the considered mass distribution with the most probable value of the orientation parameter 'Θ' (Finn & Chernoff 1993).
dEgw dfr is the energy emitted by the source per unit source frame frequency (fr), where Mc is the chirp mass of the source, and Π(fr) is given by (Ajith et al. 2008) ) 2
In Fig 3, we show how the SGWB power spectrum is impacted by the different mass scales of the BH population of the astrophysical and primordial origin in solid lines and dashed lines respectively.It can be observed that the population with higher masses exhibits a larger power density at low frequencies.This is because the presence of a Gaussian peak in the mass distribution of the BHs at MPISN predominantly leads to more contribution at the lower frequencies of the SGWB.
The presence of PBHs introduces a distinctive signature in SGWB.As the merger rate of the PBHs is expected to increase with redshift, the dominant contribution to SGWB due to PBHs will come from higher redshifts.The gravitational waves from sources at higher redshifts are going to be redshifted towards lower frequencies.This means the presence of PBHs will modify the shape of the SGWB spectrum such that we have enhanced power at low frequencies.
In Fig. 4, we compare the analytically obtained Ωgw(f ) with the Ωgw(f ) obtained using the simulation developed in this work.There is less than 1% discrepancy between the two results.

Non-stationary SGWB: off-diagonal terms in the covariance matrix
The SGWB will not be uniform over time; rather, it is expected to exhibit temporal fluctuation as discussed in Sec. 1.This fluctuation arises due to the limited number of events and the diverse properties of the gravitational wave sources.This makes SGWB a non-stationary quantity.However, the SGWB signal will exhibit time translation symmetry over a large timescale, as the astrophysical source population will remain constant over the observation time scale.So, it is important to understand at what timescales the signal homogenizes and becomes statistically stationary (time-translation symmetric).The presence of a non-stationary signal will lead to the off-diagonal terms in the covariance matrix in the frequency domain.We will discuss this below in detail for different population models.
The probability distribution of the number N of events in an interval ∆T is given by where R is the mean merger rate.The standard deviation in the number of events is given by (R ∆T ) 1/2 (Dvorkin et al. 2018;Bulik et al. 2011;Kalogera et al. 2007).Therefore the relative fluctuation will scale as (R ∆T ) −1/2 .This means there is going to be significant fluctuation in the SGWB if the merger rate is below a certain value.The distribution of SGWB and associated statistical quantities can be shown to vary with BH population parameters leading to novel non-Gaussian signatures.The fluctuation in SGWB density ∆Ωgw(f) can be defined as The fluctuation ∆Ωgw(f ) is going to be small for large merger rates.We can define three-time scales (Mukherjee & Silk 2020, 2021a): τ , which is the duration of the signal for an event, ∆tevent, which is the mean duration between two consecutive events, and ∆T , which is the observation time bin.The τ depends on the chirp mass of the source ).Typically, for BBHs τ can vary from a fraction of a second to a few seconds.The time scale ∆tevent depends on the merger rate.For a large merger rate, the interval between two events is going to be small.Hence, we are going to have overlapping events for ∆tevent < τ .The fluctuation of Ωgw(f ) under this condition is going to be very low.Therefore, we are going to encounter a fluctuating SGWB only under condition τ < ∆tevent.For a local merger rate of 30 Gpc −3 yr −1 and tmin = 100 Myr, ∆tevent ∼ 1000 seconds (for the case with only ABHs).

Model dependence of SGWB distribution
The fluctuation and distribution of the SGWB depend on the population parameters of BHs.The time dependence of the background signal can serve as additional information on the top of the power spectrum, to infer various population parameters of BHs, like merger rate and mass distribution, as well as to distinguish different formation channels such as primordial and astrophysical (Mukherjee & Silk 2020).We show in Fig. 5(a) the distribution of fluctuation in Ωgw(20 Hz) for ABHs for exposure of 10 4 seconds, for MPISN= 40 M⊙ and MPISN= 80 M⊙.It can be seen that higher mass makes the distribution of SGWB relatively more dispersed and skewed at f = 20 Hz.This is because a higher mass source generates larger strain at the low frequencies.In general, the skewness of the distribution of the background signal at different frequencies depends on the underlying population of the masses contributing to the background.A large skewness also implies there are a relatively larger number of realizations with no events.In Fig. 5 In Fig. 6, we show the distribution of the fluctuation for different observation time bins over which the signal is averaged.The distribution becomes more and more skewed as we decrease the bin size.For bins smaller than 10 4 seconds, we have a significant number of realizations where we do not have any events.The distribution of Ωgw(f ), therefore, stacks at zero, and the distribution of fluctuation in Ωgw(f ) shifts towards negative value.The key signature here is that a higher merger rate and a larger time bin will make the distribution Gaussian.Also, the change in the mass distribution will change the skewness of the distribution, as shown in Fig. 5. Therefore, the fluctuation in SGWB is going to be very critical in understanding the merger rate, mass distribution, and other population parameters including the formation channels of binaries.

Model dependence of SGWB signal covariance matrix
Apart from the fluctuation and skewness of the SGWB distribution, the SGWB is expected to exhibit correlations between signals at different frequencies.These correlations arise due to the non-stationary nature of the GW signal contributing to the background.In the case of a low merger rate, there will always be a non-zero degree of correlation between signals at different frequencies within the time scale over which the signal homogenizes.This correlation is expected to depend on population parameters, particularly the mass distribution and the merger rate.The presence of the correlation between signals at different frequencies will lead to a non-zero off-diagonal covariance matrix.Here, we provide a summary of how the signal covariance matrix is influenced by the properties of the BH population namely, (i) mass distribution and (ii) merger rate.We have restricted our analysis up to f = 100 Hz, as the signal beyond this is not well measured due to small values of the overlap reduction function (Flanagan 1993;Christensen 1992).
Impact of mass distribution on covariance matrix: If a mass distribution has lower mass BHs, we expect to see the correlation between the signals at different frequencies up to larger frequency separation as low-mass BHs emit across a broader frequency range.Likewise, the mass distribution with a small variance is also expected to show more correlation as compared to the case with a large variance in mass distribution because the sources with different masses have distinct power spectra.Therefore, the mass distribution significantly influences the overall structure of the covariance matrix.In Fig. 7 and Fig. 8, we show the covariance matrix of Ωgw(f ) for a short-time bin width of 200 seconds (over which the signal is averaged) for ABHs and PBHs respectively.For a short-time bin width of 200 seconds, only a fraction of realizations have an event.In such a case, the signals at different frequencies are highly correlated.Since the maximum frequency that sources can emit is inversely proportional to its mass, the signal covariance descends more rapidly as we move away from the diagonal term in the covariance matrix for the case with a higher mass population.Similarly, the peak in the covariance matrix shifts towards the lower frequencies for distribution with higher mass BHs because the higher mass BHs emit mostly at lower frequencies.
Impact of merger rate on covariance matrix: The merger rate has a direct impact on the magnitude of the covariance matrix.A higher merger rate corresponds to a shorter time interval, ∆tevent, between two events.The decrease in ∆tevent leads to a smaller fluctuation in SGWB, resulting in a more Gaussian distribution.However, the merger rate does not alter the shape of the covariance matrix; instead, it scales the overall magnitude of the covariance matrix.As a result, it becomes possible to explore both the mass distribution and merger rate of the sources contributing to the SGWB signal using the power spectrum and the covariance matrix.

FISHER FORECAST
In this section, we perform a Fisher analysis to show the advantage of the additional information from the covariance matrix in constraining the population parameters.The Fisher matrix analysis is a powerful statistical tool used to estimate how well a set of model parameters can be constrained based on a given set of data (Fisher 1935;Tegmark et al. 1997).The Fisher information matrix is given by where L = − ln L, where L is the likelihood function, and θi and θj are model parameters.According to the Cramér-Rao inequality (Rao 1945;Cramér 1946), the minimum uncertainty in the measurement of a model parameter is given by ∆θi ⩾ 1/ √ Fii (Tegmark et al. (1997)).For Gaussian likelihood, where ⃗ Θ ∈ {θ1, θ2, . . ., θn} represents population parameters, Ωk gw (f ) is the measured SGWB density signal from k th short-time bin, Ω m gw (f ) is the mean value of SGWB signal for parameters ⃗ Θ, and C is the covariance matrix for signal averaged over given short-time bin.The covariance matrix where CN(f, f ′ ) is the noise covariance matrix (Thrane et al. 2009;Christensen 2018), and CS(f, f ′ ) is the covariance due to intrinsic fluctuation in Ωgw(f ) that can be written as where Ω m gw (f ) is the signal averaged over given short-time bin ∆T .The intrinsic covariance of Ωgw(f ) can be useful as it provides additional information that can help us better constrain the population parameters.All the analyses of the SGWB power spectrum have disregarded the correlation between frequency bins.We demonstrate the advantage of the full covariance matrix and the impact of its off-diagonal terms on estimating source properties.
The noise covariance matrix, in this analysis, is assumed to be diagonal.The diagonal assumption is a good approximation for stationary noise (Abbott et al. (2020)).However, in the non-stationary noise limit, we can expect to see the spectral correlation in the noise power spectrum (Abbott et al. 2020;Nuttall 2018;Mozzon et al. 2022).The time scale over which the noise is stationary is much larger than the time scale of the BBH signals.As a result, the offdiagonal terms in the signal covariance matrix remain uncontaminated from the noise.Also in the presence of nonstationary noise (such as glitches), it will be uncorrelated between the different detectors (except for the Schumann resonance (Schumann 1952;P. Nguyen et al. 2021;Thrane et al. 2013Thrane et al. , 2014))).The spectrum of the Schumann resonance is different from the astrophysical SGWB signal.In this analysis, we assume the noise to be Gaussian and stationary.In future work with the LVK data, we will explore the contribution of the non-stationary noise.
The Fisher information matrix in Eq. ( 18) can be expanded as (Fisher (1935); Tegmark et al. (1997)) where T obs is the total observation time, ∆T is the shorttime bin, C,i denotes the derivative of the covariance matrix with respect to the parameters θi, and Aij= [(Ω ].The covariance matrix, C in Eq. ( 20), has contributions from two sources: (i) instrument noise and (ii) intrinsic fluctuations in Ωgw(f ).We perform the Fisher analysis on the Ωgw(f ) signal to obtain the expected constraints on the parameters.We consider two cases: (A) where we assume the signal to be stationary with no intrinsic fluctuation (power spectrum-only case), and (B) where we also include the covariance due to intrinsic fluctuation (power spectrum + covariance).We apply this analysis to the A+ sensitivity of the LIGO (Hanford & Livingston) and Virgo detectors (Abbott et al. 2018;Aasi et al. 2015;Acernese et al. 2015;Barsotti et al. 2018Barsotti et al. , 2020) ) for a total observation period of two years.In case (B), we have considered a short-time bin size of 200 seconds.In Fig. 9 and Fig. 10, we show the Gaussian Fisher posterior distribution for the above two cases.We find that for ABHs, the inclusion of the covariance matrix can improve constraints on both parameters, MPISN and tmin, by up to about 90%.For PBHs, the improvement is more significant for Mc= 30 M⊙, with an enhancement in the measurement by approximately 55%, compared to Mc = 60 M⊙, where it is only around 17%.
The Figure of Merit (FoM) can be defined as the square root of the determinant of the Fisher matrix.It is a measure of the information content and effectiveness of the parameter estimation.In Fig. 11, we show the ratio of the FoM for the power spectrum + covariance case to that for the power spectrum-only case as a function of Mc.The ratio peaks near Mc = 30 M ⊙ and decreases for both higher and lower Mc values.This is because the signal length (τ ) of low-mass sources is longer than that of high-mass sources.The longer signals make the background signal more stationary, thereby reducing the magnitude of the covariance matrix.On the other hand, high-mass sources can be individually detected up to higher redshifts, and the contribution to SGWB from higher-mass BHs comes from higher redshifts.Since the signals at higher redshifts are fainter, the strength of the covariance matrix decreases for the highermass sources.Consequently, the FoM is lower for both toosmall and too-large Mc values.In Fig. 12, we demonstrate the same ratio as a function of short-time bin, ∆T.The ratio exhibits a decreasing trend with increasing short-time bin.As the size of the short-time bin width increases, the magnitude of the spectral covariance decreases, and therefore, the spectral covariance becomes increasingly insignificant when compared to the noise power spectrum.Moreover, beyond a short-time bin width size of 10 3 seconds, the ratio converges and approaches 1.This convergence signifies that the effect of spectral covariance becomes essentially negligible beyond this threshold.The change in the bounds on the parameters with the bin size (∆T) is illustrated in Fig. 13, where we demonstrate how constraints on Mc and β change with ∆T.The constraint improves as ∆T decreases, reaching a saturation point beyond a bin size of around 200 seconds.The signal becomes statistically stationary when averaged over a large time.However, by measuring the non-stationary signal over a time window of ∆T seconds, and combining the signal from a large observational period, T obs , we can obtain a more accurate estimation of the parameters.This is equivalent to the spatial averaging of the fluctuations in the cosmic microwave temperature background.Cosmological fluctuations average to zero on average over a large spatial region leading to only a monopole signal.However, the fluctuations at small scales are prominent which can be measured from the variance of the fluctuation (Tegmark et al. 1997).

CONCLUSION
In this paper, we developed a simulation suite of the SGWB signal and demonstrated how it can result in a nonstationary SGWB signal, depending on the properties of high-redshift BHs.This can lead to a non-zero correlation between different frequencies in the signal covariance matrix.We have shown that the SGWB will have both Gaussian and non-Gaussian signatures which are impacted differently by the merger rate and mass distribution.The study of non-Gaussian and non-stationary SGWB represents a new and promising approach for probing the high-redshift properties of BBH systems.
The inclusion of the spectral correlation in SGWB in our likelihood function can significantly improve our parameter estimation.This is because GW sources with different masses and merger rates lead to different structures in the covariance matrix of the SGWB signal.A larger merger rate will make the distribution more Gaussian-like, while the mass distribution will affect both the variance and skewness of the distribution.This is a key signature in distinguishing different mass distribution and merger rate models.We found that the technique can improve the measurement by up to 90% with the A+ sensitivity of the detectors.In our analysis, we assumed that the noise at two different frequencies is uncorrelated, which is a good approximation for noise with a sufficiently long time scale over which it is stationary.It is worth noting that this analysis is moderately computa-  tionally expensive.The simulation of the Ωgw(f ) signal due to ABHs for one year of observation time with a short-time bin of ∆T = 200 seconds requires approximately 90 core hours.Making the simulation computationally affordable for practical application to data within a Bayesian framework is going to be challenging.
Using the simulation-based approach to infer the signal covariance matrix makes it accessible for the first time to utilize the additional information present in the SGWB.Although computationally expensive, incorporating this approach into the estimation of the SGWB signal from data will be highly advantageous.We have demonstrated this by applying the technique to a physical model to assess its effectiveness in understanding the underlying population.In the future, we plan to apply this technique to the LVK data from the upcoming observations to explore the stellar properties of compact objects at high redshift and also investigate its performance in the presence of glitches Our method is a versatile technique that can be applied to a variety of different situations.Observations with future advanced gravitational-wave detectors will greatly enhance our ability to study non-stationary and non-Gaussian backgrounds and extract new insights into the BH and stellar properties at high redshifts.We will implement this theoretical framework in the future to explore the feasibility of understanding the BH population and its evolution with redshift for supermassive and intermediate-mass BHs from Pulsar Timing Array (Foster & Backer 1990;McLaughlin 2013;Joshi et al. 2018;Manchester et al. 2013), space-based GW detector LISA (Amaro-Seoane et al. 2017;Hughes 2006), space-based/moon-based GW detector in Deci-Hz frequency range (Sato et al. 2017;Grimm & Harms 2020) and lowmass BHs from Einstein Telescope (Punturo et al. 2010;Sathyaprakash et al. 2012) and Cosmic Explorer (Hall 2022;Hall et al. 2021).For some GW sources, LISA, Einstein Telescope, and Cosmic Explorer will achieve exceptional multiband synergy and will be able to see deep up to a redshift of approximately z ≈ 50, and over a wide range of masses.The exploration of the time-dependent aspect of SGWB in the multi-band signal can bring a new window to understand the properties of the compact objects from high redshift which are vastly unexplored from electromagnetic observations.This technique can also be useful for distinguishing between astrophysical and cosmological SGWB signals based on the structure of the covariance matrix.

ACKNOWLEDGMENTS
The authors are thankful to Shivaraj Kandhasamy for reviewing the manuscript and providing useful comments.This work is a part of the ⟨data|theory⟩ Universe-Lab which is supported by the TIFR and the Department of Atomic Energy, Government of India.The authors would like to thank the LIGO/Virgo scientific collaboration for providing the noise curves and the computer center HPC facility at TIFR for providing computing resources.LIGO is funded by the U.S. National Science Foundation.Virgo is funded by the French Centre National de Recherche Scientifique (CNRS), the Italian Istituto Nazionale della Fisica Nucleare (INFN), and the Dutch Nikhef, with contributions by Polish and Hungarian institutes.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 would also like to acknowledge the use of the following Python packages in this work: numpy (Van Der Walt et al. 2011), scipy (Jones et al. 01 ), matplotlib (Hunter 2007), astropy (Robitaille et al. 2013;Price-Whelan et al. 2018), pygtc (Bocquet & Carter 2019), and ray (Moritz et al. 2018).

Figure 1 .
Figure 1.A schematic diagram explaining the concept behind this work.The chirping behaviour of the GW sources and its merger rate breaks the time-translation symmetry in the data.This leads to a correlation between different frequency modes in the covariance matrix in the Fourier space.The inclusion of this spectral correlation between different frequencies can improve our constraints on the parameters of the sources contributing to the background (The axes are in arbitrary units).

Figure 3 .
Figure 3. SGWB density due to ABHs for different values of the PISN mass scale M PISN (solid lines) and due to PBHs for different values of the characteristic mass scale Mc (dashed lines).The local merger rate of ABHs and PBHs is taken as 30 Gpc −3 yr −1 and 20 Gpc −3 yr −1 respectively.

Figure 5 .
Figure 5. Distribution of fluctuation in Ωgw(f ) due to (a) ABHs for different M PISN (b) PBHs for different Mc, for a short-time bin of 10 4 seconds.
(b), we also show the distribution of fluctuation in Ωgw(20 Hz) for PBHs for Mc= 30 M⊙ and Mc= 60 M⊙.

Figure 9 .
Figure 9. Corner plot showing the feasibility of measuring the M PISN and t min for (a) M PISN = 40 M ⊙ , (b) M PISN = 80 M ⊙ , and t min = 100 Myr for the ABH population.The purple color represents the case where we have assumed a stationary background and the blue color represents the case where we have included the intrinsic covariance of the SGWB.

Figure 10 .
Figure 10.Corner plot showing the feasibility of measuring Mc and β for (a) Mc= 30 M ⊙ , (b) Mc= 60 M ⊙ , and β= 1.5 for the PBH population.The orange color represents the case where we have assumed a stationary background and the blue color represents the case where we have included the intrinsic covariance of the SGWB.

Figure 11 .
Figure 11.The ratio of Figure of Merit (FoM) for the power spectrum + covariance case to that for the power spectrum only case of PBHs, as a function of Mc for a short-time bin of ∆T = 200 seconds.

Figure 12 .
Figure 12.The ratio of Figure of Merit (FoM) for the power spectrum + covariance case to that for the power spectrum-only case, as a function of the short-time bin (∆T) for PBHs with Mc= 30 M ⊙ .

Figure 13 .
Figure 13.Corner plot showing how the measurement changes with the short-time bin (∆T ) for PBHs with Mc= 30 M ⊙ and β= 1.5.
Table showing the fiducial values of the population parameters.