Estimating the angular power spectrum of the gravitational-wave background in the presence of shot noise

TL;DR

This paper introduces a new method to accurately estimate the angular power spectrum of the gravitational-wave background's anisotropies, effectively removing shot noise bias by combining independent data segments, thus enabling better insights into large-scale structure.

Contribution

A novel estimator for the gravitational-wave background anisotropy spectrum that is unbiased and efficient in the presence of shot noise, using combined independent data segments.

Findings

The estimator removes shot noise bias from the spectrum.

It is unbiased and achieves minimum variance with many data segments.

The method enhances the extraction of astrophysical information from GW data.

Abstract

There has been much recent interest in studying anisotropies in the astrophysical gravitational-wave (GW) background, as these could provide us with interesting new information about galaxy clustering and large-scale structure. However, this information is obscured by shot noise, caused by the finite number of GW sources that contribute to the background at any given time. We develop a new method for estimating the angular spectrum of anisotropies, based on the principle of combining statistically-independent data segments. We show that this gives an unbiased estimate of the true, astrophysical spectrum, removing the offset due to shot noise power, and that in the limit of many data segments, it is the most efficient (i.e. lowest-variance) estimator possible.

Equations

$$\Omega \text{quantity}(f,\hat{r})\equiv \frac{1}{{\rho }_{\mathrm{c}}}\frac{\text{differential}[3]{\rho }_{\mathrm{gw}}}{\text{differential}\text{quantity}(\ln f)\text{differential}[2]\hat{r}}=\frac{8\text{uppi}Gf}{3H_0^2}\frac{\text{differential}[3]{\rho }_{\mathrm{gw}}}{\text{differential}f\text{differential}[2]\hat{r}},$$

$$\begin{array}{rl}\bar{\Omega }\text{quantity}(\hat{r})& \equiv \text{expectationvalue}{\Omega \text{quantity}(\hat{r})}_{\Omega },\\ C\text{quantity}(\hat{r},{\hat{r}}^{\prime })& \equiv \text{expectationvalue}{\Omega \text{quantity}(\hat{r})\Omega \text{quantity}({\hat{r}}^{\prime })}_{\Omega }-\bar{\Omega }\text{quantity}(\hat{r})\bar{\Omega }\text{quantity}({\hat{r}}^{\prime }),\end{array}$$

$${\Omega }_{\ell m}\equiv {\int }_{S^2}\text{differential}[2]\hat{r}{Y}_{\ell m}^{\ast }\text{quantity}(\hat{r})\Omega \text{quantity}(\hat{r}).$$

$$\begin{array}{rl}\bar{\Omega }& \equiv \text{expectationvalue}{\Omega \text{quantity}(\hat{r})}_{\Omega }=\mathrm{constant},\\ C\text{quantity}(\hat{r}\text{dotproduct}{\hat{r}}^{\prime })& \equiv \text{expectationvalue}{\Omega \text{quantity}(\hat{r})\Omega \text{quantity}({\hat{r}}^{\prime })}_{\Omega }-{\bar{\Omega }}^2.\end{array}$$

$$\text{expectationvalue}{{\Omega }_{\ell m}}_{\Omega }={\delta }_{\ell 0}{\delta }_{m0}\sqrt{4\text{uppi}}\bar{\Omega }.$$

$$C_{\ell }\equiv {\int }_{S^2}\text{differential}[2]\hat{r}{P}_{\ell }\text{quantity}(\hat{r}\text{dotproduct}{\hat{r}}^{\prime })C\text{quantity}(\hat{r}\text{dotproduct}{\hat{r}}^{\prime }),$$

$$\text{derivative}\text{quantity}(\ln \ell )\mathrm{Var}\text{quantity}[\Omega {]}_{\Omega }\approx \frac{\ell \text{quantity}(\ell +1)}{2\text{uppi}}{C}_{\ell }.$$

$$\mathrm{Cov}\text{quantity}[{\Omega }_{\ell m},{\Omega }_{{\ell }^{\prime }{m}^{\prime }}{]}_{\Omega }={\delta }_{\ell {\ell }^{\prime }}{\delta }_{m{m}^{\prime }}{C}_{\ell },$$

$$\mathrm{Var}\text{quantity}[{\Omega }_{\ell m}{]}_{\Omega }=C_{\ell },$$

$$\text{expectationvalue}{\text{quantity}|{\Omega }_{\ell m}{|}^2}_{\Omega }=\text{quantity}|\text{expectationvalue}{{\Omega }_{\ell m}}_{\Omega }{|}^2+C_{\ell }={\delta }_{\ell 0}{\delta }_{m0}4\text{uppi}{\bar{\Omega }}^2+C_{\ell }.$$

$${\hat{C}}_{\ell }^m\equiv \text{quantity}|{\Omega }_{\ell m}{|}^2.$$

$$\mathrm{Var}\text{quantity}[{\hat{C}}_{\ell }^m{]}_{\Omega }=\text{quantity}(1+{\delta }_{m0})C_{\ell }^2.$$

$${\hat{C}}_{\ell }^{(\mathrm{std})}\equiv \frac{1}{2\ell +1}\sum _{m=-\ell }^{+\ell }\text{quantity}|{\Omega }_{\ell m}{|}^2.$$

$$\mathrm{Var}\text{quantity}[{\hat{C}}_{\ell }^{(\mathrm{std})}{]}_{\Omega }=\frac{2}{2\ell +1}{C}_{\ell }^2.$$

$$C_{\ell }\to C_{\ell }+\mathcal{W},$$

$$\text{expectationvalue}{\cdots }_{\Omega }\equiv \text{cosmological average},$$

$$\text{expectationvalue}{\cdots }_S\equiv \text{shot noise average}.$$

$$\begin{array}{rl}\text{expectationvalue}{{\Omega }_{\ell m}^{\mu }}_S& ={\Omega }_{\ell m},\\ \mathrm{Cov}\text{quantity}[{\Omega }_{\ell m}^{\mu },{\Omega }_{{\ell }^{\prime }{m}^{\prime }}^{\nu }{]}_S& ={\delta }_{\ell {\ell }^{\prime }}{\delta }_{m{m}^{\prime }}{\delta }_{\mu \nu }{\mathcal{W}}_{\tau }.\end{array}$$

$$\begin{array}{rl}\text{expectationvalue}{{\Omega }_{\ell m}^{\mu }}_{S,\Omega }& =0,\\ \mathrm{Cov}\text{quantity}[{\Omega }_{\ell m}^{\mu },{\Omega }_{{\ell }^{\prime }{m}^{\prime }}^{\nu }{]}_{S,\Omega }& ={\delta }_{\ell {\ell }^{\prime }}{\delta }_{m{m}^{\prime }}\text{quantity}(C_{\ell }+{\delta }_{\mu \nu }{\mathcal{W}}_{\tau }),\end{array}$$

$${\hat{C}}_{\ell }^{\mu \nu }\equiv \frac{1}{2\ell +1}\sum _{m=-\ell }^{+\ell }{\Omega }_{\ell m}^{\mu }{\Omega }_{\ell m}^{\nu \ast }.$$

$$\text{expectationvalue}\ast {{\hat{C}}_{\ell }^{\mu \nu }}_{S,\Omega }=C_{\ell }+{\delta }_{\mu \nu }{\mathcal{W}}_{\tau }.$$

$${\hat{C}}_{\ell }\equiv \frac{2}{N_{\tau }\text{quantity}(N_{\tau }-1)}\sum _{\mu =1}^{N_{\tau }}\sum _{\nu =\mu +1}^{N_{\tau }}{\hat{C}}_{\ell }^{\mu \nu }.$$

$$\mathrm{Var}\big{[}\hat{C}_{\ell}\big{]}_{S,\Omega}=\quantity[\frac{2}{N_{\tau}\quantity(N_{\tau}-1)\quantity(2\ell+1)}]^{2}\sum_{m=-\ell}^{+\ell}\sum_{m^{\prime}=-\ell}^{+\ell}\sum_{\mu=1}^{N_{\tau}}\sum_{\mu^{\prime}=1}^{N_{\tau}}\sum_{\nu=\mu+1}^{N_{\tau}}\sum_{\nu^{\prime}=\mu^{\prime}+1}^{N_{\tau}}\mathrm{Cov}\quantity[\Omega_{\ell m}^{\mu}\Omega_{\ell m}^{\nu*},\Omega_{\ell m^{\prime}}^{\mu^{\prime}}\Omega_{\ell m^{\prime}}^{\nu^{\prime}*}]_{S,\Omega}.$$

$$\mathrm{Var}\big{[}\hat{C}_{\ell}\big{]}_{S,\Omega}=\frac{2}{2\ell+1}\quantity[C_{\ell}^{2}+\frac{2\mathcal{W}_{\tau}C_{\ell}}{N_{\tau}}+\frac{\mathcal{W}_{\tau}^{2}}{N_{\tau}\quantity(N_{\tau}-1)}].$$

$$\mathrm{Var}\big{[}\hat{C}_{\ell}\big{]}_{S,\Omega}\simeq\frac{2}{2\ell+1}\quantity(C_{\ell}+\mathcal{W}_{T})^{2}.$$

$$\mathrm{Var}\big{[}\hat{C}_{\ell}\big{]}_{S,\Omega}\approx\frac{4}{2\ell+1}\mathcal{W}_{T}^{2}.$$

$$\text{expectationvalue}{{\Omega }_{\ell m}^{\mu }{\Omega }_{\ell m}^{\nu \ast }{\Omega }_{\ell m^{\prime }}^{{\mu }^{\prime }\ast }{\Omega }_{\ell m^{\prime }}^{{\nu }^{\prime }}}_{S,\Omega },$$

$$\Lambda \equiv \sum _{i=1}^N{\lambda }_i,$$

$$\text{expectationvalue}{\Lambda }_S=\text{expectationvalue}{\lambda }_S\text{expectationvalue}{N}_S=R{\tau }_{\mathrm{s}}\bar{n}\text{updelta}V\equiv \bar{\Lambda }.$$

$${\kappa }_n\equiv {\text{derivative}[n]K\text{quantity}(x)x|}_{x=0},$$

Full text

††preprint: KCL-PH-TH/2019-59

I Introduction

The first two observing runs of Advanced LIGO Aasi et al. (2015) and Advanced Virgo Acernese et al. (2015) have firmly established gravitational-wave (GW) astronomy as a new means of exploring the Universe. With an ever-growing number of compact binary coalescences (CBCs) being catalogued Abbott et al. (2019a); gra , attention is now increasingly turning toward other, as-yet-undetected, GW observables, such as the stochastic GW background (SGWB). The SGWB is a persistent, pseudorandom GW signal, formed from the incoherent superposition of many GW sources throughout cosmic history Allen (1996); Christensen (2019), which can be searched for by cross-correlating data between multiple GW detectors Romano and Cornish (2017); Smith and Thrane (2018).

Unlike analogous backgrounds of electromagnetic radiation [e.g., the cosmic microwave background (CMB)], the SGWB is not attenuated by intervening matter, and may give us observational access to a wealth of exotic early-Universe sources, such as cosmic strings, phase transitions, and inflationary tensor perturbations Caprini and Figueroa (2018). However, one expects the dominant contribution to the SGWB to be the combined emission from CBCs at much lower redshift. This astrophysical gravitational-wave background (AGWB) Regimbau (2011); Abbott et al. (2016, 2018) is one of the main observational targets of LIGO and Virgo as they approach design sensitivity. Since the components of CBCs are the results of stellar evolution, they should reside in galaxies,111Here we are neglecting the possibility that the black holes may be primordial; in this case, there would also be a SGWB from primordial black hole binaries Mandic et al. (2016); Clesse and García-Bellido (2017); Wang et al. (2018). and the AGWB should trace the distribution of galaxies throughout the local Universe. There has therefore been significant recent interest in AGWB anisotropies Contaldi (2017); Cusin et al. (2017); Jenkins and Sakellariadou (2018); Cusin et al. (2018); Jenkins et al. (2018, 2019); Jenkins and Sakellariadou (2019); Namikawa et al. (2016a, b) and associated observational searches Abbott et al. (2017, 2019b); Renzini and Contaldi (2019) and data-analysis methods Thrane et al. (2009); Gair et al. (2014); Romano et al. (2015); Ain et al. (2018); Renzini and Contaldi (2018); Conneely et al. (2018); Hotinli et al. (2019), as these might provide an entirely new probe of galaxy clustering and large-scale structure (LSS).

The AGWB has two immediate advantages over traditional LSS studies using galaxy surveys: first, it automatically has full sky coverage, thanks to the global beam pattern of the GW interferometers; second, it has (in principle, at least) unlimited depth, since galaxies are not selected down to some brightness threshold, but instead contribute to the AGWB at any redshift.222Note, however, that the nearest galaxies dominate the signal, as the GW flux decays as the inverse distance squared.

There exists, however, a serious obstacle to using AGWB measurements in this way. Due to the finite number of CBCs throughout the Universe, and the very short time they each spend in the LIGO/Virgo frequency band, the angular power spectrum of the AGWB is dominated by shot noise. As was first shown in Jenkins and Sakellariadou (2019), this shot noise is orders of magnitude larger than the “true” angular power spectrum, completely obscuring the much subtler correlations caused by LSS.

In this article, we present a simple but powerful method for estimating the true AGWB angular power spectrum in the presence of shot noise, exploiting the statistical independence of different shot noise realisations at different times. We begin by reviewing the definition and statistical properties of the AGWB angular power spectrum in Sec. II. In Sec. III we describe how shot noise affects the angular power spectrum, and describe the hierarchical relationship between the ensemble of shot noise realisations and the ensemble of LSS realisations. In Sec. IV we construct an unbiased estimator of the true, astrophysical power spectrum, and derive a simple expression for its variance (some details of this calculation are given in Appendix A). In an appropriate limit, this estimator is the minimum-variance unbiased estimator (MVUE) for the $C_{\ell}$’s in the presence of shot noise, which we show in Appendix B.

Our discussion focuses on the AGWB composed of CBCs, as this will likely be the first stochastic GW signal to be detected. However, all of our main results hold for any stochastic source that is comprised of a finite number of transient events.

II The SGWB angular power spectrum

The SGWB is described by the density parameter

[TABLE]

i.e., the GW energy density $\rho_{\mathrm{gw}}$ arriving from direction $\hat{{\bf\it r}}$ on the sky, with observed frequency in a logarithmic bin around $f$, measured in units of the cosmological critical energy density $\rho_{\mathrm{c}}\equiv 3H_{0}^{2}/\quantity(8\uppi G)$. (Here and throughout, we use units in which $c=1$.) We treat $\Omega\quantity(f,\hat{{\bf\it r}})$ as a random field on the sphere for each frequency $f$. Our analysis holds for each frequency bin independently, so we drop the explicit $f$ dependence below.

Approximating the SGWB as a Gaussian random field (GRF), it is fully specified by its first two moments,

[TABLE]

which we call the mean and the two-point correlation function (2PCF), respectively. Here we have introduced a cosmological ensemble average $\expectationvalue{\cdots}_{\Omega}$, which we discuss further in Sec. III. It is also convenient to define the spherical harmonic components (SHCs),

[TABLE]

These form a set of complex random variables that fully specify any given random realisation of the field, and allow us to analyse the statistics of the SGWB on different angular scales.

For extragalactic sources, one expects the SGWB to be statistically isotropic, so that the mean is the same everywhere on the sky, and the 2PCF depends only on the angle between the two points,

[TABLE]

An immediate consequence of the first line is that the SHCs all have zero mean, except the monopole,333It is important to stress the difference between the mean $\bar{\Omega}$ and the monopole $\Omega_{00}$. Both describe the isotropic component of the SGWB, but the former is a fixed deterministic quantity that can be predicted from theory, while the latter is a random variable that must be measured. Different cosmological observers will measure different values of the monopole, drawn from a distribution that has $\bar{\Omega}$ as its mean. This distinction is usually not made in the case of the CMB, where fluctuations in the temperature monopole over different realisations are negligible. However, the SGWB monopole can vary significantly between different random realisations Meacher et al. (2014), particularly due to shot noise.

[TABLE]

(The factor of $\sqrt{4\uppi}$ here is just the integral of $Y_{00}$ over the sphere.)

The statistical properties of the anisotropies are conveniently described by a multipole expansion of the 2PCF,

[TABLE]

which we call the angular power spectrum. Together with $\bar{\Omega}$, this completely characterises a statistically isotropic GRF, with each $\ell$ corresponding to anisotropies of angular size $\approx\uppi/\ell$. More precisely, the contribution to the total angular variance of the SGWB per logarithmic angular scale can be approximated as

[TABLE]

The $C_{\ell}$’s also have a very natural interpretation in terms of the SHCs, (3). It is possible to show that for a statistically isotropic background,444Note that here, and throughout, we define the covariance of two random variables with a complex conjugate on the second argument, i.e. $\mathrm{Cov}\quantity[X,Y]\equiv\expectationvalue{XY^{*}}-\expectationvalue{X}\expectationvalue{Y^{*}}=\mathrm{Cov}\quantity[Y,X]^{*}$. This ensures that the variance is always real, since $\mathrm{Var}\quantity[X]\equiv\mathrm{Cov}\quantity[X,X]=\expectationvalue*{\quantity|X|^{2}}-\quantity|\expectationvalue{X}|^{2}$.

[TABLE]

and in particular,

[TABLE]

i.e., each $\Omega_{\ell m}$ is uncorrelated with the others, and has variance $C_{\ell}$.

II.1 The standard $C_{\ell}$ estimator

We want to measure the SGWB angular power spectrum using a network of GW detectors. The standard method for this exploits (9), which can be rewritten as

[TABLE]

We see immediately that for $\ell>0$, we can use the absolute-square of each SHC as an estimator of the corresponding multipole,

[TABLE]

This $m$-dependent “naive” estimator is unbiased, in the sense that $\expectationvalue*{\hat{C}^{m}_{\ell}}_{\Omega}=C_{\ell}$. We can calculate its variance explicitly in the case where the SGWB is a GRF by using Isserlis’ theorem to give

[TABLE]

The $m$ index runs over $-\ell\leq m\leq+\ell$, so in fact we have $2\ell+1$ different naive estimators, each with complementary statistical information. We can therefore construct a better estimator by averaging over all of these. This defines the standard $C_{\ell}$ estimator,

[TABLE]

The variance in this case is significantly lower, particularly at large $\ell$,

[TABLE]

Note that this $1/\quantity(2\ell+1)$ scaling is exactly what one should expect for $2\ell+1$ independent measurements of the same quantity. The variance in (14) is in fact the minimum possible variance for any estimator of the angular power spectrum in the absence of shot noise, as it saturates the Cramér-Rao bound (see Appendix B).

For the case $\ell=0$, we see from (10) that the standard estimator is biased, $\expectationvalue*{\hat{C}_{0}^{(\mathrm{std})}}_{\Omega}=4\uppi\bar{\Omega}^{2}+C_{0}$. The problem here is that one wants to estimate the cosmological mean ($\sqrt{4\uppi}\bar{\Omega}$) and variance ($C_{0}$) of the observed monopole, $\Omega_{00}$, but with only a single cosmological realisation. This is impossible, as the sample variance is undefined for a single sample; the best one can do is estimate $\bar{\Omega}$. For this reason, the rest of this article will focus on higher multipoles, $\ell>0$.

Of course, in practice, the SHCs themselves must first be estimated from the strain data in the interferometers; this in itself is a nontrivial data analysis challenge Romano and Cornish (2017). However, the details of this process are unimportant for the purposes of this article.

III Shot noise and hierarchical averaging

When modelling the AGWB, one usually averages over the merger times of CBCs and the spatial locations of their host galaxies to give a smooth merger rate $R$ and galaxy number density $n$. However, any observed realisation of the AGWB will inevitably consist of a finite number of CBCs, emitted from a finite number of galaxies. (This was first investigated in the context of the AGWB monopole in Meacher et al. (2014).)

It was recently shown Jenkins and Sakellariadou (2019) that for a SGWB composed of a finite number of sources, the measured angular power spectrum becomes

[TABLE]

where $C_{\ell}$ is the angular power of the intrinsic, astrophysical anisotropy, and $\mathcal{W}$ represents the additional shot-noise power. The latter is independent of $\ell$—i.e., the shot noise is spectrally white in harmonic space. If the SGWB is composed of $N$ persistent localised sources, then schematically the shot noise scales as $\mathcal{W}\sim\bar{\Omega}^{2}/N$. If the sources are transient (such as CBCs) then this becomes $\mathcal{W}_{\tau}\sim\bar{\Omega}^{2}/\quantity(R\tau)$, with $R$ the event rate and $\tau$ the time interval used to measure the stochastic GW energy density. In the latter case, we write the shot noise power as $\mathcal{W}_{\tau}$ to emphasise that it depends on the observer’s choice of $\tau$ (subject to the condition that $\tau\gg 1/f$, so that the appropriate Fourier component can be reliably measured). Figure 1 gives a visualisation of how the shot noise power affects the GW intensity distribution on the sky.

It was shown in Jenkins and Sakellariadou (2019) that, for the AGWB, the temporal shot noise associated with finite CBCs per observation time is significantly larger than the true $C_{\ell}$ spectrum, while the spatial shot noise associated with a finite galaxy number density is significantly smaller; we therefore focus exclusively on the former in this analysis. This allows us to exploit the fact that we have observational access to multiple realisations of the temporal shot noise: for each successive observation interval $\tau$, we can obtain a set of SHCs with shot noise power $\mathcal{W}_{\tau}$. The shot-noise fluctuations in each set of SHCs is associated with a different set of CBCs in a different set of galaxies, so it is immediately clear that each successive shot noise realisation is statistically independent. This is in contrast with the spatial shot noise, for which we can only observe a single realisation (as it is tied to the random positions of galaxies).

There are two logically distinct random processes that govern the observed SHCs: the distribution of matter on large scales, and the emission of a finite number of GW signals from this matter distribution in a given observation period. We model these processes, and their corresponding ensemble averages, in a hierarchical manner:

The true $\Omega_{\ell m}$ are drawn from Gaussian distributions with variance $C_{\ell}$. This process is associated with the cosmological averaging operation from Sec. II,

[TABLE]

which can be thought of as an average over an “ensemble of Universes,” with each Universe having a distinct random realisation of LSS. Of course, we only have access to a single such realisation. 2. 2.

The true $\Omega_{\ell m}$ are modulated by shot noise, so that a set of “noisy” $\Omega^{\mu}_{\ell m}$ is drawn from a distribution with the true components $\Omega_{\ell m}$ as its mean. This draw is independent for each time interval, with different intervals being labelled by the index $\mu$. We write the associated average over shot noise realisations as

[TABLE]

Physically, this reflects the fact that a given cosmological distribution of galaxies can correspond to many different SGWB realisations, as the number and times-of-arrival of transient GW signals from each galaxy are essentially random.

Figure 2 is a schematic representation of the two averaging procedures.

More quantitatively, the first two moments of the noisy SHCs under the shot noise average for a fixed realisation of LSS are

[TABLE]

The first line is true by definition, while the second states that each SHC in each time interval has equal shot noise power, and is uncorrelated with all the others.555The lack of correlation between time intervals is due to the acausal relationship between distant GW sources, while the lack of correlation between different $\ell m$ is due to statistical isotropy. Averaging also over realisations of LSS, we find

[TABLE]

where we have used (5) and (8), and have introduced the shorthand $\expectationvalue{\cdots}_{S,\Omega}\equiv\expectationvalue{\expectationvalue{\cdots}_{S}}_{\Omega}$.666Note that $\mathrm{Cov}\quantity[X,Y]_{S,\Omega}\equiv\expectationvalue{XY^{*}}_{S,\Omega}-\expectationvalue{X}_{S,\Omega}\expectationvalue{Y^{*}}_{S,\Omega}$, which is not equal to $\expectationvalue{\mathrm{Cov}\quantity[X,Y]_{S}}_{\Omega}$. (Recall that $\ell>0$, so the SHCs all have zero mean.)

IV Mitigating the shot noise

In this section, we use (19) to define a function of the noisy SHCs $\Omega_{\ell m}^{\mu}$ that is an unbiased estimator of the true angular power spectrum $C_{\ell}$ in the presence of shot noise. We then show that (in the appropriate limit) this is the minimum-variance unbiased estimator (MVUE).

We start by modifying the standard autocorrelation estimator (13), forming a set of cross-correlations between different time intervals,

[TABLE]

These are unbiased if and only if $\mu\neq\nu$,

[TABLE]

We can combine these estimators in much the same way that we combine multiple “naive” estimators (11) to form the standard estimator (13). Suppose that our total observing time is $T$, so that there are $N_{\tau}\equiv T/\tau$ time segments. Then we have $N_{\tau}\quantity(N_{\tau}-1)/2$ pairs $\mu,\nu$ for which $\mu\neq\nu$. Summing over these, we define the combined estimator

[TABLE]

This is unbiased, $\expectationvalue*{\hat{C}_{\ell}}_{S,\Omega}=C_{\ell}$, with variance given by

[TABLE]

Evaluating (23) requires us to evaluate the fourth moment of the noisy SHCs. This would be trivial if the SHCs were all Gaussian, but we must account for the Poisson-like nature of the shot noise. In Appendix A we calculate the fourth moment using the same statistical model for the CBC rate density as in Jenkins and Sakellariadou (2019); this results in

[TABLE]

This expression is tied to the fact that we have excluded the on-diagonal terms $\mu=\nu$ when constructing (22); otherwise, there would be additional contributions to the variance (the estimator would also no longer be unbiased). Note that since $\mathcal{W}_{\tau}\propto 1/\tau$, we have $\mathcal{W}_{\tau}/N_{\tau}\propto 1/T$. This means that we can’t “win” by decreasing the length of the data segments $\tau$, only by increasing the total observing time $T$. In fact, writing $\mathcal{W}_{T}=\mathcal{W}_{\tau}\quantity(\tau/T)=\mathcal{W}_{\tau}/N_{\tau}$, we see that in the limit where $N_{\tau}\gg 1$, (24) becomes

[TABLE]

This is exactly the standard cosmic variance expression from (14), but with $C_{\ell}$ replaced by $C_{\ell}+\mathcal{W}_{T}$. In Appendix B, we show that this is in fact the minimum possible variance of any unbiased estimator for the $C_{\ell}$’s in the presence of shot noise, saturating the Cramér-Rao bound Kay (1993). The estimator (22) is therefore the MVUE in the limit $N_{\tau}\gg 1$.

At the opposite extreme, for the minimum number of segments, $N_{\tau}=2$, the variance is nearly twice as large (taking $\mathcal{W}_{T}\gg C_{\ell}$),

[TABLE]

An illustrative example for $N_{\tau}=10$ is shown in Fig. 3.

The term $\mathcal{W}_{T}$ in (25) is the same as that appearing in the mean of the standard estimator, $\expectationvalue*{C_{\ell}^{(\mathrm{std})}}_{S,\Omega}=C_{\ell}+\mathcal{W}_{T}$, so this new optimal estimator is still affected by the presence of shot noise; the crucial improvement is that the shot noise only adds to the variance of the estimator, and does not bias the spectrum as in the standard case.

V Conclusion

We have developed a new method for estimating the true, astrophysical angular power spectrum of a statistically-isotropic AGWB, in a way that is not biased by the presence of temporal shot noise. The estimator (22) is based on the simple principle of cross-correlating between statistically-independent time intervals. It is straightforward to compute from LIGO/Virgo strain data using existing pipelines, and in the limit of a large number of time intervals (regardless of their length), it achieves the lowest possible variance of any estimator on the same data. This is a very significant improvement upon the standard analysis, given that the shot noise power is expected to exceed the astrophysical angular power by a factor of $\approx 10^{4}$ Jenkins and Sakellariadou (2019). We expect that these results will be highly valuable for future GW observing campaigns, and may give us access to interesting and novel cosmological information.

Acknowledgements.

We thank Andrew Matas for reading the manuscript and providing us with useful feedback. A.C.J. thanks Eugene Lim for valuable discussions. This article has been assigned document number LIGO-P1900209. Some of the results in this article have been derived using the HEALPix package Gorski et al. (2005). A.C.J. is supported by King’s College London through a Graduate Teaching Scholarship. J.D.R. acknowledges support from National Science Foundation Award No. PHY-1505861 and start-up funds from Texas Tech University. M.S. is supported in part by the Science and Technology Facility Council (STFC), United Kingdom, under the research grant No. ST/P000258/1.

Appendix A Fourth moment of the spherical harmonic components with shot noise

Here we compute the fourth moment of the noisy SHCs,

[TABLE]

using the same statistical model for the AGWB as in Jenkins and Sakellariadou (2019).

In a given volume element $\updelta V$, the number of galaxies is modelled as a Poisson random variable $N\sim\mathrm{Pois}\quantity[\bar{n}\updelta V]$, where $\bar{n}$ is the mean galaxy number density,777A more sophisticated approach would account for the statistical properties of halos Cooray and Sheth (2002). However, this simple assumption is sufficient for the calculation here. and with the number of CBCs in each galaxy forming a set of independent and identically distributed Poisson random variables $\lambda_{i}\sim\mathrm{Pois}\quantity[R\tau_{\mathrm{s}}]$, where $R$ is the mean rate per galaxy, and $\tau_{\mathrm{s}}=\tau/\quantity(1+z)$ is the source-frame time interval of the observation. The total CBC count in the volume element is then

[TABLE]

which follows a compound Poisson distribution, with mean

[TABLE]

Here we characterise this distribution by its cumulants,

[TABLE]

where $K\quantity(x)$ is the cumulant-generating function,

[TABLE]

The first two cumulants, $\kappa_{1},\kappa_{2}$, are just the mean and the variance. Higher cumulants $\kappa_{n}$ represent the “connected components” of the $n$th moments of the distribution, and vanish if the distribution is Gaussian. The fourth moment of $\Lambda$ is given in terms of cumulants as

[TABLE]

Consider now a set of different volume elements observed at different times $\updelta V_{\mu}$, whose total CBC counts $\Lambda_{\mu}$ are independent and identically distributed. Their fourth moment is then

[TABLE]

where $\delta_{\mu\nu\rho}$ is non-zero only if $\mu=\nu=\rho$, etc. The first four cumulants are

[TABLE]

The rate per galaxy is typically on the order of $\mathrm{Myr}^{-1}$, while the source-frame observing time is at most a few years, so we we have $R\tau_{\mathrm{s}}\sim 10^{-6}$. So to a very good approximation, we can take

[TABLE]

In order to compute the fourth moment of the SHCs $\Omega_{\ell m}^{\mu}$, we translate from the CBC number count to the comoving CBC rate density,

[TABLE]

and replace $\bar{\Lambda}$ with $\bar{\mathcal{R}}\equiv\bar{n}R$. We introduce ${\bf\it\zeta}$ to represent the parameters of the CBC (masses, spins, …) and of the galaxy (star formation rate, metallicity, …). Making this replacement, we obtain

[TABLE]

where here the delta symbols are shorthand for

[TABLE]

i.e., two GW sources must coincide in space, in time, and in parameter space, in order to contribute to the shot noise.

Now, using the equation for the AGWB in terms of the CBC rate density from Jenkins and Sakellariadou (2019), we write

[TABLE]

where $Y_{\ell m}^{\mu}$ is shorthand for $Y_{\ell m}\quantity(\hat{{\bf\it r}}_{\mu})$, etc., and $\mathcal{S}\quantity(\nu_{\mathrm{s}},{\bf\it\zeta})$ encodes the GW emission of each CBC as a function of source-frame frequency $\nu_{\mathrm{s}}$, and has dimensions of $\quantity[\mathrm{time}]^{4}$ in units where $c=1$. Using (37) and the properties of the spherical harmonics, this becomes

[TABLE]

where the shorthand (38) is now no longer used. The coefficients on the 3- and 4-point terms are defined as

[TABLE]

and are just combinations of Wigner $3j$ symbols. In addition to the 2-point term $\mathcal{W}_{\tau}$ we defined before, we now have 3-point and 4-point terms that appear,

[TABLE]

Taking the cosmological average of (LABEL:eq:Omega-4pt) and subtracting the 2nd moments, we therefore find

[TABLE]

Note that the term proportional to $\mathcal{X}_{\tau}$ has vanished, due to the SHCs having zero mean for $\ell>0$ (recall that we do not consider $\ell=0$ here). The term proportional to $\mathcal{Y}_{\tau}$ is associated with the fourth cumulant $\kappa_{4}$, and is therefore a sign of the non-Gaussian nature of the shot noise fluctuations (as $\kappa_{n}=0$ for $n>2$ in the Gaussian case). However, this term does not end up contributing to the variance of the estimator $\hat{C}_{\ell}$, as this includes only off-diagonal pairs $\mu\neq\nu$ and $\mu^{\prime}\neq\nu^{\prime}$ by design.

Appendix B Minimum-variance estimation of the angular power spectrum

Here we show that the estimator in (22) is the MVUE for the true angular power spectrum $C_{\ell}$ in the presence of shot noise. To do so, we use the fact that the variance of any estimator $\hat{\theta}\quantity(x)$ for a parameter $\theta$ of a probability distribution $p\quantity(x|\theta)$ obeys the Cramér-Rao bound Kay (1993),

[TABLE]

where $\mathcal{L}\quantity(x|\theta)\equiv\ln p\quantity(x|\theta)$ is the log-likelihood.

We start by considering the case with no shot noise. Assuming Gaussianity, the joint log-likelihood for the SHCs is then

[TABLE]

and we find that

[TABLE]

This shows that the standard cosmic variance expression (14) is the best one can do in the zero-shot-noise case, as it saturates the Cramér-Rao bound.

Now we include shot noise, and consider the noisy SHCs, $\Omega_{\ell m}^{\mu}$. Though we know these are not Gaussian (see Appendix A), the Gaussian case is by far the most tractable, so we consider it first. The joint Gaussian log-likelihood is fully specified by (19),

[TABLE]

where ${\bf\it\Omega}_{\ell m}=(\Omega^{1}_{\ell m},\ldots,\Omega^{N_{\tau}}_{\ell m})$ is a vector of the noisy SHCs for a given $\ell m$, and $\mathsf{C}_{\ell}$ is the corresponding $N_{\tau}\times N_{\tau}$ covariance matrix,

[TABLE]

One can show that this has determinant

[TABLE]

and inverse

[TABLE]

using the matrix determinant lemma and the Sherman-Morrison formula, respectively. Taking the second derivative with respect to $C_{\ell}$, we find

[TABLE]

so the Cramér-Rao bound is the same as before, but with $C_{\ell}\to C_{\ell}+\mathcal{W}_{T}$. This is exactly the variance we derived for our estimator in (25).

Equation (51) was derived using the Gaussian log-likelihood (47), and one may worry about whether it holds in the case we are interested in, given that the noisy SHCs do not follow a Gaussian distribution. However, we have seen that our estimator saturates this bound in the limit $N_{\tau}\gg 1$. As shown in Jaynes (2003), the only probability distribution that saturates the Cramér-Rao bound under a given set of constraints [e.g., the constraints on the first two moments in (19)] is that which maximises the entropy under those constraints. The maximum-entropy distribution with fixed variance is a Gaussian Jaynes (2003), so this shows that our estimator must be Gaussian in the limit $N_{\tau}\gg 1$, and that the calculations above are valid in that limit. (The approach to Gaussianity for $N_{\tau}\gg 1$ can also be shown using the central limit theorem.)