Probing the large scale structure using gravitational-wave observations of binary black holes

Third generation gravitational-wave (GW) detectors are expected to detect a large number of binary black holes (BBHs) to large redshifts, opening up an independent probe of the large scale structure using their clustering. This probe will be complementary to the probes using galaxy clustering -- GW events could be observed up to very large redshifts ($z \sim 10$) although the source localization will be much poorer at large distances ($\sim$ tens of square degrees). We explore the possibility of probing the large scale structure from the spatial distribution of the observed BBH population, using their two-point (auto)correlation function. We find that we can estimate the bias factor of population of BBH (up to $z \sim 0.7$) with a few years of observations with these detectors. Our method relies solely on the source-location posteriors obtained from the GW events and does not require any information from electromagnetic observations. This will help in identifying the type of galaxies that host the BBH population, thus shedding light on their origins.


I. INTRODUCTION
Gravitational-wave (GW) observations by LIGO and Virgo have opened a new era of astronomy [1].On the completion of the third observing run, the LIGO-Virgo-KAGRA collaboration published a combined catalog of GW transients (GWTC-3 [2][3][4]), which reported ∼ 90 significant detections of GWs from compact binaries mergers.Independent analyses of the public LIGO-Virgo data have revealed a few additional events in the same data set [5][6][7][8].Detection of GW events has become a routine now and GW sky is filling up rapidly.
The dominant sources of GW signals for the LIGO-Virgo-KAGRA detectors [9][10][11] is the merger of compact objects.Several hundreds to thousands of such observations are expected in the next few years [12].LIGO-India [13,14] is expected to come online sometime during next decade, coinciding with the upgraded Advanced LIGO (A+) detectors [15].These additional detectors will significantly improve the localization of binary mergers.There are several ongoing efforts to build the next generation of ground based detectors.Proposals for next generation detectors include that of i) LIGO Voyager, which is expected to observe binary neutron stars (BNSs) up to a horizon redshift of z ∼ 0.5 [16,17], ii) Einstein Telescope (ET), which is expected to have a BNS horizon of z ∼ 2 [16,18], and iii) Cosmic Explorer (CE) with an expected BNS horizon of z ∼ 20 [16,19].Third generation (3G) detectors like CE and ET will have sensitivity that will be an order of magnitude better than that of Advanced LIGO and will be sensitive to frequencies as low as 1Hz.
Although the network configuration and the sensitivity of the proposed 3G detectors are not finalized, studies of various configurations and their implication on source localization and parameter estimation suggest that these detectors will be able to observe binary black hole (BBH) mergers up to large redshifts (detection horizon up to z ∼ 100) [16].For the redshift range z ∈ [0, 3], a significant fraction of BBH mergers can be localized to within 1 square degree [20].These observations will, in essence, create a "survey" of well-localized GW events across the sky, akin to galaxy surveys.Since this survey would trace the underlying large-scale structure of the universe, it is natural to ask what properties of the large-scale structure GW events track, and what that tells us about the underlying astrophysics describing the binaries.In this work, we investigate the possibility of detecting one feature of the large-scale structure (LSS) tracked by GW events, the large-scale bias, using 3G detector network solely from GW observations.This paper is organized as follows: In section II, we discuss the current status of cosmology using GW observations.section III describes the methods of probing LSS features using galaxy observations and how can we extend it to GW observations with 3G detectors.In section IV, we describe the simulations performed in this study and their results.Section V summarizes the results.

II. COSMOLOGY USING GWS
During the second observing run, LIGO and Virgo detected GWs from a BNS merger, GW170817 [21], for the first time.Electromagnetic (EM) counterparts of this event were also detected by several telescopes [22][23][24], which enabled the identification of the host galaxy of the merger.This led to a precise measurement of the redshift of GW170817 and the first measurement of the Hubble constant H 0 from GW observations [25].BNS detections expected in the near future with EM counterparts should improve the precision of this measurement, potentially contributing to resolving the apparent tension between the Planck measurement of H 0 [26] and that from type Ia supernovae [27].Some studies also explore the techniques of cross-matching or cross-correlating galaxy catalogs with BBH observations to constrain H 0 [28][29][30][31][32][33][34][35][36].
With a large number of GW detections expected in the near future, we will have a population of BBH and BNS mergers distributed over a large redshift range, providing a new tracer of the large scale structure.Recent studies show that by cross correlating the GW events with galaxy catalogs, the large scale structure can be probed by estimating the linear bias [37,38] or by the lensing of GWs [39].In this work, we explore the possibility of probing the clustering of BBHs by estimating their two-point (auto)correlation function.If these mergers happen in specific types of galaxies, the clustering of the BBHs should trace that of such galaxies.If, for some reason, BBHs are predominantly distributed outside galaxies, their clustering information should reveal this.Thus, an independent estimation of the clustering of BBHs offers an interesting probe of not only the large scale structure, but also the astrophysical environment of the mergers.

III. LARGE SCALE STRUCTURES OF THE UNIVERSE
The two-point correlation function (2PCF) ξ(r) is related to the excess probability δP (r), above what is expected for a random distribution, of finding a pair of objects (e.g., galaxies or, in the context of this work, BBH mergers) separated by distance r.This can be expressed as where n is the number of objects per unit volume and dV is the volume element.For the matter overdensity field δ(x) := ρ(x)/ρ − 1, where ρ(x) is the local matter density and ρ the mean matter density of the Universe, the 2PCF is given by where angle brackets denote the ensemble average which, in turn, can be estimated by averaging over a large volume.The above equation assumes statistical homogeneity and isotropy of the universe, hence ξ is only a function of the magnitude r of the separation vector y − x between the two points x and y.In general, the 2PCF is also a function of the redshift z.However, when we restrict ourselves to a relatively narrow redshift bin ∆z, it can be assumed to be a constant within that redshift range.The distribution of the galaxies in the Universe is expected to trace the underlying matter distribution.At large scales, to a good approximation, the 2PCF of the galaxies ξ gal (r) is related to that of matter ξ m (r) through a simple relation [40] where b gal is the galaxy bias, taken to be scale-independent.Usually, the value of b gal depends on the luminosity and color type of galaxies [41].Similarly, we can also define a bias which quantifies the clustering of the observed BBH population: If we are able to measure b BBH from GW observations, this would allow us to compare it against b gal estimated from other observations (e.g., EM galaxy surveys), thus providing hints to the host environments of the BBH mergers.

A. Estimating the BBH correlation function from GW observations
The interpretation of ξ(r) as the excess probability of finding points separated by a distance r allows one to construct fast estimators of the correlation function from data.The Landy-Szalay (LS) estimator [42] is the most commonly-used estimator, and is given by, Here, DD(r) denotes the number of point-pairs in the data (galaxy catalogs or BBH mergers) separated by a distance r, RR(r) denotes the number of point-pairs in an equal-sized simulated random catalog separated by a distance r, and DR(r) denotes the number of data-random point pairs separated by a distance r.Since the size of the simulated random data set is something we have control over, we can choose to have more number of points in the random catalog.If N D , N R are the number of points in the data and random catalogs respectively, for general N D , N R , Eq.( 5) gets modified [42] to, The correlation function of the galaxies from a survey can be estimated using Eq.( 6).With next generation GW detectors like ET and CE, we expect to detect BBH mergers up to large redshifts.If the number of detections are sufficiently large, we can use their localization information to study how these GW events are clustered by estimating the correlation function ξ BBH (r). .The radial direction corresponds to comoving distance and the angular direction corresponds to RA (the dec coordinate is projected out).BBH events are distributed according to the input power spectrum and bias factor (= 1.5) in each redshift bin.The errors on localization are drawn from a probability distribution described in the text.

B. Smearing of the correlation function due to GW localization errors
The challenge in estimating ξ BBH (r) is that the precision in the GW source localization (sky location and distance) will be poor as compared to the galaxy localization (which can effectively be described as a point in the survey volume).Due to the large statistical uncertainties in the GW localization, the observed correlation function of BBHs will be modified from the actual correlation function -the poor source localization distributes weights from the points of actual location to a smeared field around those points.The "smearing" of the correlation function will depend on the distribution of the GW localization uncertainties from the population.The smeared correlation function (Figure 1) can be computed by convolving the actual correlation function with the ensemble averaged localization posteriors obtained from GW data.We describe this below.
In the absence of any measurement errors, the probability distribution P tr (µ) of the location µ of BBH mergers is given by where δ (3) is the three dimensional Dirac delta function, µ i denotes the three-dimensional location of BBH i, and N is the total number of BBHs in the survey volume V such that

∆N(z)
3. Solid curve with shaded region shows the total number of merger events as a function of redshift in shell of thickness ∼ 350 h −1 Mpc and ∼ 500 h −1 Mpc in comoving distance.These numbers are calculated by assuming the redshift distribution of BBHs from [44] and the local merger rates of BBHs estimated in [45].The dashed lines show the average number of mergers in the shell of given thickness for which the errors in sky localization are within a degree square and errors in estimating the comoving distance are ≤ 90 h −1 Mpc for a network of three 3G detectors.field is given by 1 Above, Ptr = 1/V denotes the volume-averaged probability density.The correlation function between two points µ and ν in the field of density contrast is given by where ⟨⟩ denotes ensemble averages.Using Eq.( 7), we can write (10) Now we investigate how the true correlation function ξ tr (µ, ν) gets smeared by the presence of measurement uncertainties.Assuming that the localization posteriors follow Gaussian distributions, where µ i is the true location of the i th BBH, C i is the covariance matrix for the corresponding localization posterior (assumed to be diagonal), and ∆µ i is the scatter induced by the detector noise.In the absence of systematic biases ∆µ i will be distributed according to a Gaussian distribution of mean zero and covariance matrix C i .We now marginalize P i (x − µ i , ∆µ i ) over ∆µ i : This averaging can be performed on the posterior (as opposed to the final correlation function) since the noise-induced shifts ∆µ i are uncorrelated with the BBH locations µ i .The resulting posterior P i (x − µ i ) is a Gaussian distribution with mean µ i and covariance matrix 2C i .Using the property of Dirac delta function, P i (x − µ i ) can be rewritten as and the probability distribution of the location of a population of BBH mergers is given by The correlation function between two points x and y of this probability field is given by Note that each term in the sum over i, j is equal to the joint posterior probability of the BBH mergers i and j to take the positions x and y, respectively.In a frequentist interpretation, this is equivalent to the joint probability of drawing two samples of x and y from the posteriors of the two events.The relation between this and our simulations should be apparent now.We only want to consider correlations between two different BBH mergers; thus, we restrict the sum to i ̸ = j.Now, Using Eq.( 13), this can be rewritten as Now, we make the following assumptions: 1. Assuming that the posterior distributions are uncorrelated with the actual location of mergers (uniform sky coverage assumption), we can write: 2. Since P i and P j are posterior probability distributions estimated from two independent GW events (uncorrelated noise), 3. Motivated by the homogeneity of space, we assume Using these assumptions Eq.( 16) can be rewritten as where we have used Eq.( 10) for the last step.The smeared correlation function of the probability density contrast field δ P (x) := P (x)/ P − 1 is given by Using Eqs.( 17) and ( 9), this can be rewritten as dV ν P (x−µ) P (y−ν) ξ tr (µ, ν) (19) This can be used to compute the smeared correlation function ξ(x, y) from the true correlation function ξ tr (µ, ν).Essentially we convolve the true correlation function ξ tr by a smoothing function (ensemble-averaged localization posteriors).Due to the homogeneity and isotropy of space, the true correlation function only depends on the magnitude of the difference of its arguments ξ tr (|ν − µ|).We can exploit this by transforming to new integration variables, a := x − µ and b := y − ν to get: where This shows that the smeared correlation function also depends only on the separation of points x and y, i.e., ξ(y − x).However, unlike the true correlation function, the direction also matters unless posterior functions are isotropic.In the case of BBH mergers, we expect the radial uncertainty to be much larger than the angular uncertainties and therefore we cannot demand isotropy and thus the orientation of x − y matters.To deal with this, we average Eq.( 20) over all orientations for a given r := |x − y| to find the spherically averaged correlation function ξ(r).We choose the volume of interest to be large enough to permit every possible orientation with minimal bias.
In this work, we simulate the combined probability distribution of BBHs by placing GW posteriors around the true BBH locations, after introducing a noise-induced scatter in the mean of the posteriors.The posterior distributions of right ascension (RA), declination (dec) and comoving distance (d), estimated from N simulated events are combined to create a normalized combined posterior probability field , where x = {RA, dec, d}.Assuming that the localization posteriors follow Gaussian distributions, ) where µ i is the true location of the i th BBH and C i is the covariance matrix of the corresponding localization posterior (assumed to be diagonal), while N = 1/ (2π) 3 |C i | is a normalization constant.Note that the individual posteriors will in general not be centered around the true BBH locations, because of the scatter ∆µ i introduced by the detector noise.This random scatter is drawn from a mean-zero Gaussian distribution of covariance matrix C i .Figure 2 shows the P (x) from a simulated catalog of BBH observations.

IV. SIMULATIONS AND RESULTS
To put our method to test, we use the publicly available code LOGNORMAL_GALAXIES [46], to simulate galaxy catalogs at various redshifts with input power spectrum taken as the matter power spectrum approximated by the fitting function of Eisenstein and Hu [43] and consistent with the Planck-18 cosmological parameters [26].This code enables one to generate mock galaxy catalogs assuming lognormal probability density function for the matter field and galaxies.We assume that GW events occur in any random subsample of the galaxies in the catalog, which essentially implies b BBH = b gal .We simulate three different catalogs having input linear bias b gal = [1.2,1.5, 2.0].Since it is possible to directly infer the distance (within the localization errors) to the BBH from the GW observations, peculiar velocities of galaxies will not play any role (unlike EM galaxy surveys where the distance is inferred from the redshift).Hence, while generating the catalogs, we switched off peculiar velocities in the code.
We then simulate the mock BBH catalogs using the steps outlined below and check whether we are able to recover the bias consistent with the input value.For simplicity, we have assumed the input b gal to be redshift-independent, however, our conclusions on the recovery of the bias would remain unchanged even if we use an evolving bias.These are the steps involved: 1. Choose a shell of thickness 350 h −1 Mpc around the given redshift.The value was chosen so that we have enough events in the shell, and the actual correlation function does not vary appreciably within the redshift bin.The extent ∆z of the redshift bin corresponding to this shell thickness at redshift z = 0.3 (1.0) turns out to be 0.13 (0.2)2 .
2. Randomly select N galaxies from this shell as proxy for GW events and put localization error bars on each event assuming Gaussian posteriors (see below for details regarding uncertainties).
3. Select one point from each of the N posteriors.This simulates a particular realization of galaxy locations.
Use LS estimator to estimate the correlation function.
We repeat this process 1000 times and take the average to get ξ(r).
4. To estimate the variance, we create 50 galaxy catalogs corresponding to different realization of the cosmic matter field to account for cosmic variance.For each of these catalogs, we select 20 sub-catalogs N of random galaxies each to account for fluctuations due to sampling, thus amounting to a total of 1000 sub-catalogs.One subcatalog catalog was taken as realization of our universe and ξ BBH (r) was estimated using steps described above.
Error bars on ξ BBH (r) was placed making use of the scatter estimated from other sub-catalogs.
5. Estimate the bias factor b BBH by comparing the recovered correlation function ξ BBH (r) with the smeared matter correlation function ξ m sm (r).We estimate the correlation function in the range of comoving distance r ∈ [10, 50] ∼ h −1 Mpc as this is well within the chosen shell thickness and the linear bias approximation is valid in this range.To find the fit for the bias factor we then define a χ 2 function, where is the smeared matter correlation function for the given distribution of localization errors, b is the bias factor, ξ est (r i ) is the correlation function estimated from the simulated catalog for each r-bin r i , and Σ ij is the covariance matrix between i-th and j-th bin.We estimate the covariance matrix using 1000 sub-catalogs we generated for this study.We define the likelihood function as, and use Bayesian analysis to estimate posterior distribution corresponding to the likelihood function Eq. ( 23) for a uniform prior for parameter b in range b ∈ [0, 5].
We use the open source nested sampling based sampler DYNESTY [47] for this purpose.Since we used the galaxies as proxy for GW events, we expect the recovered b BBH to be consistent with the input bias factor used for simulating the galaxy distribution.Clearly, this method will be valid only when the errors in localization do not exceed the range of comoving distances we are trying to probe.This translates into the requirement that errors in RA, dec should be within a degree and errors in the comoving distances should not exceed few tens of Mpc.To find if this requirement can be fulfilled with 3G detectors, we perform GW parameter estimation studies using a population of BBH events distributed up to redshift ∼ 1.2 with 3G detector network 2CE-ET (CE locations: one in USA and one in Australia.ET location: proposed one in Europe).In Table I, we list the location and and low frequency cutoff used for 3G detectors network.In our simulations, we use BBH population with the powerlaw plus a peak distribution of primary masses p(m 1 ) ∝ m −α 1 with α = 2.3 [48].We use the IMR-PHENOMPV2 [49] waveform available in the LALSUITE [50]  software package along with the appropriate detector PSDs [51,52] to simulate our signals, and use the PYCBCINFER-ENCE package [53] to determine distribution of localization errors.We find that a significant fraction of events up to z ≃ 1 fulfills this requirement.Figure 3 shows number of expected BBH mergers (using the BBH merger rate given in [48]) at various redshifts for one year of observations along with fraction of events that are expected to be localized well enough for this type of study.In our simulations, this selection introduces no significant biases; however possible selection effects need to be considered for the actual analysis.The distribution of the widths of the 68% credible regions of the marginalized posteriors on RA, dec and comoving distance can be approximated by truncated Gaussian distributions with mean {µ RA = 0. We neglect the correlations between the errors in RA, dec and distance.
Figure 4 shows the smeared correlation function compared to the estimated correlation function from a simulation using 5000 GW observation in a shell around the redshift z = 0.3.Figure 5 shows the bias factor recovered from different redshift bins using different observation duration (5, 7 and 10 years).The estimated b BBH , in general, are consistent with the simulated bias within error-bars.The small number of events where the actual value is outside the error bars is consistent with statistical fluctuations.Note that even with a moderate observational time of five years, we can recover the bias to within ∼ 20% for z ≲ 0.7.Due to large spread of localization volumes, the bias recovery becomes difficult for high redshift 3 .We would like to point out that if the localization volumes with 3G detectors can be improved by either better sensitivity or improvements in waveform modelling, the results presented in this study will improve.
In order to test the robustness of this method, we simulated 1000 catalogues with galaxies distributed with b gal = 1.2, 1.5, and 2.0.For reference, we also generated a catalog with no underlying correlation function i.e. galaxies are distributed randomly (corresponding to b = 0) at redshift z = 0.5.In Figure 6, we show the stacked posterior samples from all the individual runs.The width of stacked posterior samples depends on the localization volume distribution, and sampling errors.We note that the stacked posteriors are peaked around the injected values.These reference distribution, for given localization volumes, can be used to assign significance to a particular bias recovery measurement with respect to the random distribution.
There are various ways one can use the recovered b BBH (z) to understand the properties of the event hosts.For example, one can compare the clustering properties of the galaxies as measured from the optical surveys with b BBH (z) and obtain insights on the type of galaxies that host these merger events.Further, the recovered bias can also be related to the host dark matter halo mass [57].In general, if one assumes that the typical masses of the haloes hosting these GW events do not evolve with redshift, one can predict the redshift-dependence of b BBH (z) for a given cosmological model.This then can be compared with the observations to understand the formation channels of the BBHs.
We would like to emphasize that in order to convert luminosity distance samples obtained from GW localization volumes, to either comoving volume, or redshift, we assume an underlying cosmological model: ΛCDM model with Planck 2018 [26] values as implemented in ASTROPY [58,59].The results presented herein are intrinsically tied to this chosen cosmological model.In order to pursue a model independent approach, one needs to marginalize over cosmological parameters.Furthermore, when dealing with real-world data, additional factors and effects come into play.These may include the development of optimal methodologies for estimating the covariance ma-trix between distance bins, incorporation of detector response functions, consideration of weak lensing effects, and more.We intend to investigate these effects in our future work.

V. SUMMARY
In this work we explored the possibility of probing large scale structure with BBH observations using third generation GW detectors.We showed that bias factor can be estimated using clustering information of BBH events with 5-10 years of observations.This can be achieved solely from the GW observations, without requiring EM counterparts or galaxy catalogs.The bias factor b BBH estimated from various redshifts will enable us to find whether the BBH mergers track the distribution of specific types of galaxies, or dark matter halos.Although the statistical precision of the estimated bias b BBH is weaker than that of the galaxy bias obtained from EM galaxy surveys, it is important to note that the GW-based analysis probes the underlying matter distribution using a novel astrophysical tracer, thus enabling an independent probe of the large scale structure.We intend to extend this analysis to include effects such as selection bias, and method of cross correlating with galaxy catalogs to probe higher redshift, etc.

9 FIG. 1 .
FIG.1.The "smeared" correlation function (dashed lines) and the "true" correlation function (solid lines) for various redshifts.The true correlation function is simply the matter correlation function calculated using Eisenstein-Hu prescription[43] for the standard model of cosmology.The "smearing" of the correlation function due to measurement errors is calculated assuming that the distribution of errors in localization of GW population follow a Gaussian distribution with mean {µ RA = 0.5 • , µ dec = 0.5 • , µ d = 50 h −1 Mpc} and standard deviation {σ RA = 0.5 • , σ dec = 0.5 • , σ d = 20h −1 Mpc}.

FIG. 2 .
FIG.2.An example of probability field obtained from localization posteriors from a realization of simulated catalog of BBH observations in redshift range z ∈ [0.1, 1.1].The radial direction corresponds to comoving distance and the angular direction corresponds to RA (the dec coordinate is projected out).BBH events are distributed according to the input power spectrum and bias factor (= 1.5) in each redshift bin.The errors on localization are drawn from a probability distribution described in the text.

FIG. 4 .
FIG.4.Smeared correlation function for a given distribution of localization errors is plotted along with the one recovered from simulated events at redshift 0.3 and input bias factor of 1.5.Smeared correlation function is scaled with input bias for comparison.We used 5000 simulated events distributed in a shell of thickness 350 h −1 Mpc around the given redshift.

FIG. 5 .FIG. 6 .
FIG.5.The recovered bias factor b BBH from various redshifts bins (with shell thickness of ∼ 350 h −1 Mpc).The catalogs were created using the matter power spectrum of Eisenstein-Hu with different values of linear bias 1.2 (left) 1.5 (middle) and 2.0 (right).Each subplot shows the estimated bias factor, along with the corresponding error bars (68% confidence regions), using GW observations of BBHs over a period of 5, 7 and 10 years.

TABLE I .
[54][55][56]ations of 3G detector network (location, low frequency cutoff f low ) considered in this study.We use the design sensitivity noise curves for ET and CE.These detectors configuration for CE and ET are also used in previous works[54][55][56].