Stochastic Gravitational-Wave Backgrounds: Current Detection Efforts and Future Prospects

Abstract

The collection of individually resolvable gravitational wave (GW) events makes up a tiny fraction of all GW signals that reach our detectors, while most lie below the confusion limit and are undetected. Similarly to voices in a crowded room, the collection of unresolved signals gives rise to a background that is well-described via stochastic variables and, hence, referred to as the stochastic GW background (SGWB). In this review, we provide an overview of stochastic GW signals and characterise them based on features of interest such as generation processes and observational properties. We then review the current detection strategies for stochastic backgrounds, offering a ready-to-use manual for stochastic GW searches in real data. In the process, we distinguish between interferometric measurements of GWs, either by ground-based or space-based laser interferometers, and timing-residuals analyses with pulsar timing arrays (PTAs). These detection methods have been applied to real data both by large GW collaborations and smaller research groups, and the most recent and instructive results are reported here. We close this review with an outlook on future observations with third generation detectors, space-based interferometers, and potential noninterferometric detection methods proposed in the literature.

1. Introduction

Gravitational waves (GWs) are perturbations of the spacetime metric caused by extremely energetic events throughout the Universe. Up until now, direct detections of GWs have been coherent measurements of resolved waveforms in detector datastreams which may be traced back to single, point-like sources. These make up a tiny fraction of the gravitational-wave sky: The vast collection of unresolved signals corresponding to multiple point sources or extended sources adds up incoherently, giving rise to gravitational-wave backgrounds (GWBs). A variety of different backgrounds is expected given the range of GW sources in the Universe; however, regardless of their origin, most of these are treated as stochastic, as they may be described by a non-deterministic strain signal and are, hence, referred to as stochastic gravitational-wave backgrounds (SGWBs).

Some backgrounds are stochastic by generation processes (e.g., inflationary tensor modes), whereas others are stochastic due to the characteristics or limitations of the specific detector used to observe them (e.g., the cumulative signal from many binary black hole coalescences). An SGWB of this latter nature is by definition at the threshold of detection, making it effectively a detector-dependent observable. Determining whether a signal is a background and also whether it inherently behaves as a stochastic field is then an iterative process, where first all signals present in the detector timestream undergo signal-to-noise ratio (SNR) estimation, which determines their resolvability, and then the background is estimated as the cumulative sum of only the unresolved and/or sub-threshold signals. In practice, each GW detector will be able to access qualitatively different backgrounds, depending on the noise levels and frequency ranges probed.

With this review, we intend to provide a clear and complete yet concise reference for the broader community interested in stochastic search efforts while surveying our field of research. In laying out stochastic background detection techniques, we focus on three categories of detectors: ground-based detector networks, pulsar timing arrays (PTAs) monitored by radio telescopes, and spaceborne detectors comprised of sets of satellites. We do not discuss in detail the impressive work on data acquisition and processing necessary to put the data in the form required for the implementation of the searches we report (see Section 4 and Section 5); we refer the reader to several other papers that delineate these efforts (see, e.g., [1,2,3,4,5] for recent discussions of ground-based laser interferometer detector characterization and calibration and [6,7,8] for reviews of PTA experiments).

We consider the present second generation (2G) ground-based detector network to be made up of three detectors that are fully operational: the NSF-funded Laser Interferometer Gravitational-wave Observatory (LIGO) pair of detectors in the United States [9] and the Virgo detector in Italy [10]. To these, we can add the GEO600 [11] detector in Germany, which, however, does not have comparable sensitivity to the others given its smaller size, and the Kamioka Gravitational-wave Detector (KAGRA) in Japan [12], which has not reached initial design sensitivity yet. Future detectors will include the Indian Initiative in Gravitational-wave Observations [13]. Beyond these, we discuss the potential of the upcoming third generation (3G) detector network, which will include the Einstein Telscope (ET) [14,15] and Cosmic Explorer (CE) [16,17,18]. All these monitor frequencies around ${10}^2$ Hz, spanning roughly two orders of magnitude in total, ∼10–10${}^3$ Hz, depending on specific characteristics.

Regarding nanohertz-frequency stochastic backgrounds, $f\sim {10}^{-9}$ Hz, we outline the methodology of pulsar timing array (PTA) experiments and present results from the North American Nanohertz Observatory for Gravitational Waves (NANOGrav [19]), the Parkes Pulsar Timing Array (PPTA [20]), the European Pulsar Timing Array ([EPTA [21,22]), and the consortium of these collaborations known as the International Pulsar Timing Array (IPTA [23,24]), of which the Indian PTA (InPTA [25]) is also now a member. We comment on the common-spectrum noise process reported by the above collaborations that may be related to stochastic GWs or similar noise spectra in timing array pulsars [26]. We describe standard search methods for an isotropic GWB, and we also briefly discuss methods for resolving anisotropies in the GWB. We further comment on the role of future observation efforts such as the Square Kilometer Array (SKA [27]) and the five-hundered-meter aperture spherical radio telescope (FAST [28]).

As for spaceborne interferometers, we primarily discuss the potential detection capabilities of the Laser Interferometer Space Antenna (LISA), which is an ESA-led mission in collaboration with NASA planned to launch in the mid 2030s. LISA will scan a broad range of frequencies, $f\sim {10}^{-5}$–${10}^{-1}$ Hz, with a maximum sensitivity around a millihertz, allowing us to tune into a new, vast range of the GW sky at unprecedented depth and volume [29].

We will start by introducing the theoretical underpinnings of GWB science in Section 2. In Section 3, we review expected sources of GWBs typically considered in the literature. We outline the different properties of the backgrounds in relation to the detectors that measure GWs today and those that will probe them in the foreseeable future. We offer a detailed explanation of stochastic search methods in Section 4, differentiating between different observing strategies and detection regimes. In Section 5, we summarise current detection results, reporting the most stringent upper limits yet on SGWBs while highlighting how close we are to detection. Finally, we close this review in Section 6 with an overview of future prospects with the upcoming GW detectors.

2. Theory of Stochastic Backgrounds

In this section, we provide the theoretical foundations of GW science essential to the description of generation and propagation of GWs in the Universe. First, we review the mathematical description of stochastic backgrounds in a common framework, regardless of their astrophysical or cosmological origin, which lies at the basis of (almost) all stochastic search methods reviewed in Section 4. Then, we introduce the fractional energy density of GWs, which is the main observable targeted by stochastic searches.

2.1. Gravitational-Wave Strain and Stokes Parameters

We work in the linearised regime of general relativity, such that the spacetime metric is close to that of flat Minkowski space, with small perturbations that encode GWs, $g_{\mu \nu }={\eta }_{\mu \nu }+h_{\mu \nu }$. Neglecting terms of order $h^2$, the vacuum Einstein field equation can then be written as a wave equation in the De Donder gauge [30],

(1)$\square (h_{\mu \nu }-\frac{1}{2}{\eta }_{\mu \nu }{h}_{\alpha }^{\alpha })=0,$

where the box operator is the D’Alembert operator, defined as $\square ={\partial }_{\mu }{\partial }^{\mu }$. Solutions to this equation may be constructed as linear superpositions of plane waves,

(2)$h_{\mu \nu }=A{e}_{\mu \nu }{e}^{i\left(k_{\alpha }{x}^{\alpha }\right)},$

where $e_{\mu \nu }$ is the normalised polarisation tensor, and $k^{\alpha }$ is the wave four-vector. By inserting this solution in the equation above, one recovers the dispersion relation $k_{\alpha }{k}^{\alpha }=0$. This implies that the waves travel at the speed of light in the linearised regime. We may then pick the transverse-traceless (TT) gauge and reduce $h_{\mu \nu }$ to a purely spatial perturbation, $h_{ij}$, which carries two degrees of freedom or equivalently two independent, spin-2, polarisation states.

Generalising Equation (2) as in [31], for example, the GW strain at time t and position vector $x$ given by an infinite superposition of plane waves incoming from all directions on the sky $\widehat{\mathit{n}}$,

(3)$\widehat{\mathit{n}}=(sin\theta cos\varphi ,sin\theta sin\varphi ,cos\theta ),$

in standard angular coordinates $(\theta ,\varphi )$ on the 2-sphere may be written as

(4)$h_{ij}(t,\mathit{x})={\int }_{-\infty }^{+\infty }df{\int }_{S^2}d\widehat{\mathit{n}}\sum _{P=+,\times }{h}_P(f,\widehat{\mathit{n}}){ϵ}_{ij}^P\left(\widehat{\mathit{n}}\right)e^{i2\pi f(\widehat{\mathit{n}}·\mathit{x}-t)}.$

The spatial wave vector is written explicitly as $\mathit{k}=2\pi f\widehat{\mathit{n}}$ (we will use $c=1$ throughout unless otherwise noted), P labels the polarisation states, and $h_P$ are the respective amplitudes assigned to each state. Choosing the linear polarisation basis $\{+,\times \}$, the orthogonal polarisation base tensors ${ϵ}^P$ may be written as

(5)$\begin{array}{cc}\hfill {\mathbf{ϵ}}^{+}& ={\mathit{e}}_{\theta }\otimes {\mathit{e}}_{\theta }-{\mathit{e}}_{\varphi }\otimes {\mathit{e}}_{\varphi },\hfill \end{array}$

(6)$\begin{array}{cc}\hfill {\mathbf{ϵ}}^{\times }& ={\mathit{e}}_{\theta }\otimes {\mathit{e}}_{\varphi }+{\mathit{e}}_{\varphi }\otimes {\mathit{e}}_{\theta },\hfill \end{array}$

with

(7)$\begin{array}{cc}\hfill {\mathit{e}}_{\theta }& =(cos\theta cos\varphi ,cos\theta sin\varphi ,-sin\theta ),\hfill \end{array}$

(8)$\begin{array}{cc}\hfill {\mathit{e}}_{\varphi }& =(-sin\varphi ,cos\varphi ,0),\hfill \end{array}$

as defined in [32]. For a stochastic signal, $h_P$ represents a random complex number. Most often, these are assumed to be drawn from a Gaussian probability distribution. This is likely a good approximation for a cosmological GWB, as we will see in Section 3. For a stochastic GWB of astrophysical sources, this assumption may break down but the central limit theorem will guarantee that the statistics approach with respect to a Gaussian random field if the signal is sourced by a sufficiently large number of independent events and that any high signal-to-noise outliers have been subtracted from the detector’s timestreams. Under the Gaussian assumption, the statistical properties of the amplitudes are then characterised solely by the second moments of the distribution $⟨h_P^{}(f,\widehat{\mathit{n}})h_{P^{\prime }}^{\star }(f^{\prime },{\widehat{\mathit{n}}}^{\prime })⟩$, which, assuming statistical homogeneity, correspond to ensemble averages as [33]

(9)$\begin{array}{cc}\hfill \left(\begin{array}{cc}⟨h_{+}^{}(f,\widehat{\mathit{n}})h_{+}^{\star }(f^{\prime },{\widehat{\mathit{n}}}^{\prime })⟩& ⟨h_{+}^{}(f,\widehat{\mathit{n}})h_{\times }^{\star }(f^{\prime },{\widehat{\mathit{n}}}^{\prime })⟩\\ ⟨h_{\times }^{}(f,\widehat{\mathit{n}})h_{+}^{\star }(f^{\prime },{\widehat{\mathit{n}}}^{\prime })⟩& ⟨h_{\times }^{}(f,\widehat{\mathit{n}})h_{\times }^{\star }(f^{\prime },{\widehat{\mathit{n}}}^{\prime })⟩\end{array}\right)={\delta }^{\left(2\right)}(\mathit{n},{\mathit{n}}^{\prime })\delta (f-f^{\prime })\times & \\ \hfill \left(\begin{array}{cc}I(f,\widehat{\mathit{n}})+Q(f,\widehat{\mathit{n}})& U(f,\widehat{\mathit{n}})-iV(f,\widehat{\mathit{n}})\\ U(f,\widehat{\mathit{n}})+iV(f,\widehat{\mathit{n}})& I(f,\widehat{\mathit{n}})-Q(f,\widehat{\mathit{n}})\end{array}\right)& ,\hfill \end{array}$

where the Stokes parameters I, the intensity, Q and U, giving the linear polarisation, and V, the circular polarisation, have been introduced. The four Stokes parameters completely describe the polarisation of the observed signal in analogy with electromagnetic Stokes parameters for the photon [34]. The difference here is that whilst the electromagnetic Q and U Stokes parameters transform as spin-2 quantities with respect to rotations, their GW strain counterparts transform as spin-4. In both cases, intensity I behaves as a scalar under rotations, while V transforms as a pseudo-scalar. More details on these and their spin-weighted behaviour may be found in [35].

2.2. The Energy Density of Gravitational Waves

The evolution of our Universe is described in terms of its expansion rate, also referred to as the Hubble rate,

(10)$H(t)\equiv \frac{\mathrm{d}}{\mathrm{d}t}lna,$

where $a\left(t\right)$ is the cosmological scale factor. Using general relativity, we can write this rate in terms of the different forms of energy density present in the Universe (this is most conveniently performed as a function of redshift z),

(11)$H\left(z\right)=H_0{({\Omega }_{\mathrm{R}}{(1+z)}^4+{\Omega }_{\mathrm{M}}{(1+z)}^3+{\Omega }_k{(1+z)}^2+{\Omega }_{\Lambda })}^{1/2},$

where $H_0\approx 70$ km s${}^{-1}$ Mpc${}^{-1}$ is the Hubble rate today. In Equation (11), ${\Omega }_i={\rho }_i/{\rho }_{\mathrm{c}}$ are the density parameters of each energy component of the Universe: radiation $\mathrm{R}$, which includes photons and relativistic neutrinos; matter $\mathrm{M}$, which includes baryons and cold dark matter; curvature k; and the cosmological constant $\Lambda$. The critical energy density ${\rho }_{\mathrm{c}}$ is a normalisation that can be interpreted as the energy density required to close our (flat) Universe today, which is explicitly

(12)${\rho }_{\mathrm{c}}=\frac{3H_0^2}{8\pi G},$

where $G=6.67\times {10}^{-11}$ m${}^3$ kg${}^{-1}$ s${}^{-1}$ is the universal gravitational constant.

At this stage, one may ask the following: where do GWs fit in? In fact, GWs feature in Equation (11), as a part of the radiation energy density, as they are considered to be carried by massless particles propagating at the speed of light.¹ This allows us to write

(13)${\Omega }_{\mathrm{R}}={\Omega }_{\gamma }+{\Omega }_{\nu }+{\Omega }_{\mathrm{GW}}+\cdots ,$

with the terms in the sum representing the energy density in photons, neutrinos, and GWs, respectively. While the present day value of the radiation density is very small (${\Omega }_{\mathrm{R}}\sim {10}^{-4}$), we see from Equation (11) that at very high redshift, $z\gtrsim {10}^3$, it dominated the Universe’s expansion. As a result, any additional contributions to the radiation density, including those from primordial GWs, will leave an imprint on cosmological observables generated at these early epochs: principally, the cosmic microwave background (CMB) and the light element abundances predicted by Big Bang nucleosynthesis (BBN) [37]. These imprints are identical for all forms of additional radiation density; thus, by convention they are parameterised in terms of the effective number of neutrino species, $N_{\mathrm{eff}}$. Constraints on $N_{\mathrm{eff}}$ can straightforwardly be converted into constraints on GWB energy density. The most up-to-date analysis combines CMB and BBN data with external datasets such as measurements of the baryon acoustic oscillation scale, producing [38]

(14)${\Omega }_{\mathrm{GW}}<1.2\times {10}^{-6}{\left(\frac{H_0}{100\mathrm{km}{\mathrm{s}}^{-1}{\mathrm{Mpc}}^{-1}}\right)}^{-2}.$

Given these strong constraints on ${\Omega }_{\mathrm{GW}}$, it is reasonable to consider GWs as perturbations in the Universe, which do not influence the bigger picture but obey the propagation rules, the structure, and the energy content of the Universe that have been set out.

The density parameter appearing in Equation (14) encapsulates the total energy density of GWs at all frequencies,² as this is what determines the Hubble expansion rate. However, all other searches for gravitational waves focus on particular frequency bands; thus, it is convenient to define the GW energy density frequency spectrum [39],

(15)${\Omega }_{\mathrm{GW}}\left(f\right)=\frac{1}{{\rho }_{\mathrm{c}}}\frac{d{\rho }_{\mathrm{GW}}\left(f\right)}{dlnf};$

note that the frequency f here is the measured frequency, assuming an observer and the intervention of an instrument which mediates the true signal and the measurement. Furthermore, beyond the isotropic value of GWB over the entire sky it is possible to include any anisotropy by adding a directional dependence ${\Omega }_{\mathrm{GW}}\left(f\right)\to {\Omega }_{\mathrm{GW}}(f,\widehat{\mathit{n}})$ where $\widehat{\mathit{n}}$ is the unit vector of a line of sight.

It is now possible to recover the cosmological information contained in the ${\Omega }_{\mathrm{GW}}$ parameter, which is effectively the cumulative GW energy history weighted per redshifted frequency bin, as presented in [40]. This may be observed by noting that homogeneity and isotropy throughout the universe imply

(16)${\rho }_{\mathrm{GW}}={\int }_0^{\infty }\frac{df}{f}{\int }_0^{\infty }dz\frac{N\left(z\right)}{1+z}{f}_{\mathrm{r}}\frac{d{E}_{\mathrm{GW}}}{d{f}_{\mathrm{r}}},$

where z is redshift, $N\left(z\right)$ is the number of GW emitters as a function of redshift, and $f_{\mathrm{r}}$ is the emission frequency of GWs in the rest frame of the emitter, $f_{\mathrm{r}}=f(1+z)$. $d{E}_{\mathrm{GW}}/d{f}_{\mathrm{r}}$ is then the spectral energy density of the source population. Equation (15) may be rewritten as

(17)${\rho }_{\mathrm{GW}}={\int }_0^{\infty }\frac{df}{f}{\rho }_{\mathrm{c}}{\Omega }_{\mathrm{GW}}\left(f\right),$

hence, equating Equations (15) and (16) frequency by frequency yields

(18)${\Omega }_{\mathrm{GW}}\left(f\right)=\frac{1}{{\rho }_{\mathrm{c}}}{\int }_0^{\infty }dz\frac{N\left(z\right)}{1+z}{\left[f_{\mathrm{r}}\frac{d{E}_{\mathrm{GW}}}{d{f}_{\mathrm{r}}}\right]}_{f_{\mathrm{r}}=f(1+z)}.$

Equation (18) implies that the fractional energy density of GWs per measured logarithmic frequency intervals at redshift zero, ${\Omega }_{\mathrm{GW}}$, is proportional to the integral over the cosmic history of the comoving number density of the sources multiplied by the emitted energy fraction of each source, in the appropriately redshifted frequency bin. This directly relates ${\Omega }_{\mathrm{GW}}$ to the source event rate as a function of redshift, $N\left(z\right)$. In the case of different types of sources, this may be extended to a sum over all contributors,

(19)$N\left(z\right){\left[f_{\mathrm{r}}\frac{d{E}_{\mathrm{GW}}}{d{f}_{\mathrm{r}}}\right]}_{f_{\mathrm{r}}=f(1+z)}⟶\sum _i{N}_i\left(z\right){\left[f_{\mathrm{r}}\frac{d{E}_{\mathrm{GW},i}}{d{f}_{\mathrm{r}}}\right]}_{f_{\mathrm{r}}=f(1+z)}.$

${\Omega }_{\mathrm{GW}}\left(f\right)$ may also be re-written in terms of the intensity parameter I defined in Equation (9) [32].

(20)${\Omega }_{\mathrm{GW}}\left(f\right)=\frac{4{\pi }^2{f}^3}{{\rho }_{\mathrm{c}}G}I\left(f\right).$

We can derive this relationship using the Isaacson formula for the energy density of the GWB [41,42],

(21)${\rho }_{\mathrm{GW}}=\frac{1}{32\pi G}⟨{\dot{h}}_{ij}{\dot{h}}_{ij}⟩,$

where the angle brackets indicate an average over many wavelengths. Inserting the plane-wave decomposition in Equation (4) and rewriting the two-point statistics of each plane wave in terms of the Stokes parameters in Equation (9), we obtain

(22)${\rho }_{\mathrm{GW}}=\frac{\pi }{2G}{\int }_0^{\infty }\mathrm{d}f{\int }_{S^2}\mathrm{d}\hat{\mathit{n}}{f}^{2}I(f,\hat{\mathit{n}}).$

Note that we have assumed ergodicity here such that the spatial average in Equation (21) is equivalent to the ensemble average in Equation (9). If we treat intensity as isotropic, $I(f,\hat{\mathit{n}})\to I\left(f\right)$, then we immediately recover Equation (20). This fundamental relation allows us to connect GWB observations to cosmological implications of the background and may easily be extended to a directional statement as

(23)${\Omega }_{\mathrm{GW}}(f,\widehat{\mathit{n}})=\frac{4{\pi }^2{f}^3}{{\rho }_{\mathrm{c}}G}I(f,\widehat{\mathit{n}}).$

Note that this implies

(24)${\Omega }_{\mathrm{GW}}\left(f\right)=\frac{1}{4\pi }{\int }_{S^2}\mathrm{d}\widehat{\mathit{n}}{\Omega }_{\mathrm{GW}}(f,\widehat{\mathit{n}}),$

however, one must take care of factors of $4\pi$ as other conventions are sometimes used. Note that we are using the conventions as in [32], where all power spectra are one-sided.

As for the other Stokes’ parameters, these do not contribute to ${\Omega }_{\mathrm{GW}}$ and are equal to zero in the case of an unpolarised background. This may be easily gleaned by their definitions, in Equation (9), as the second moments of the $h_{+}$ and $h_{\times }$ fields should be equal and the expectation value of any cross correlations between the two should vanish. ${\Omega }_{\mathrm{GW}}$, thus, inherits the spectral properties from the GWB intensity as per Equation (20). As will emerge in the following sections, the spectral dependence is often reduced to a power law, in which case it is common practice to write³

(25)${\Omega }_{\mathrm{GW}}\left(f\right)={\Omega }_{\mathrm{GW}}\left(f_{\mathrm{ref}}\right){\left(\frac{f}{f_{\mathrm{ref}}}\right)}^{\alpha },$

where $\alpha$ is the spectral index of the signal, and $f_{\mathrm{ref}}$ is a reference frequency. This assumption is often extended to be valid in each direction independently, such that the spectral and directional dependence of the GW signal factor out as

(26)$I(f,\widehat{\mathit{n}})=I(f_{\mathrm{ref}},\widehat{\mathit{n}}){\left(\frac{f}{f_{\mathrm{ref}}}\right)}^{\alpha -3}.$

This is particularly useful when performing multiparameter estimations such as map-making as it allows one to integrate out the spectral dependence of the reconstructed signal and solve for an intensity sky map. Note that in the ensemble, averaging is required to recover the Stokes’ parameters, and all temporal phase information in the GWs is discarded; as argued in [43], no physical information may be recovered from the coherent phase component.

3. Sources of Stochastic Backgrounds

GWBs fall broadly into two distinct categories based on the underlying generation mechanism. The first is a superposition of signals of astrophysical origin, galactic and/or extra-galactic, that are not individually detected or resolved. This type of background will, thus, be unresolved and incoherent in that the temporal phase of the signal is expected to be random and is expected to be stochastic in the limit where the number of sources is very large. The second category is cosmological or primordial backgrounds. These include GWBs generated by the evolution of vacuum fluctuations during an inflationary epoch that end up as superhorizon tensor modes at the end of inflation [44,45]. It also includes GWBs generated by nonlinear phenomena at early times such as phase transitions or via emission by topological defects [46,47,48]. For the most part, primordial backgrounds are formed via stochastic processes, resulting in stochastic signals; on top of this, propagation effects contribute to isotropising a primordial signal, in most cases erasing any initial anisotropies [49,50]. In particular, stochasticity here implies no initial information is retained by the temporal phase of the background [51]. Cosmological GWBs can also be stochastic due to many overlapping signals, for example, in the case of cosmic strings. Assuming that the orientation of individual astrophysical sources is random on the sky, these are in general expected to be unpolarised, as an equal distribution in the polarisation modes will cancel out any local asymmetry. Primordial signals are also mostly assumed to be unpolarised, as the underlying stochastic process is not expected to select a preferential frame, with some exceptions as we explain below. What then characterises backgrounds of different origin is the spectral dependence of the source energy density in Equation (18): the power spectrum of the signal depends on the specific GW emission mechanism, which is imprinted in the measured strain. Note that this measurement includes non-negligible selection effects, as qualitatively different backgrounds contribute from different redshift shells and from different directions.

In this section, we review both astrophysical and cosmological GWBs, providing the necessary background for the targeted searches discussed in Section 5. We also comment on the observational properties of the signal which are essential to understand when building an optimal search method. The various sources are also summarised in Figure 1, which includes the sensitivity of several GW detection efforts for reference.

3.1. Astrophysical Backgrounds

Astrophysical GWBs are the collection of all GWs generated by astrophysical processes which are individually unresolved by your GW detector. These can be either individual subthreshold signals, or they can be so numerous that they add up incoherently and form a continuous signal in the timestream.

Perhaps the most studied signal in the literature is a background sourced by a collection of inspiralling and merging compact binary systems. These include black hole binaries, neutron star binaries, white dwarf binaries, and systems counting a mixed pair of these objects. Black hole binaries in particular are a vast category of sources, as the mass of each black hole in the binary ranges between 5 and 10${}^{10}{M}_{\odot }$. The lower mass of $5M_{\odot }$ here refers to the lack of black holes between the Tolman–Oppenheimer–Volkoff limit (which sets the maximum mass of a neutron star to $2.9M_{\odot }$ [57]) and $5M_{\odot }$, as seen in low-mass X-ray binary observations [58,59]⁴. Specifically, a BH is considered massive ($M_{\mathrm{BH}}\sim {10}^5$–${10}^{10}{M}_{\odot }$), intermediate ($M_{\mathrm{BH}}\sim {10}^2$–${10}^5{M}_{\odot }$), and stellar mass ($M_{\mathrm{BH}}\sim 5$–${10}^2{M}_{\odot }$) [61]. There are binaries that have a mass ratio close to one, where the two black holes in the binary fall within the same mass category, and there are also so-called extreme mass ratio (EMR) binaries where the two objects belong to different families, including extreme mixed binaries such as massive black hole–neutron star pairs. There is a direct relation between binary masses and GW emission frequencies; however, overall, the GW power spectra have a similar shape. In particular, in the early emission stage of the binary, i.e., the inspiral phase when the energy emission increases adiabatically as the binary’s orbit shrinks very gradually, energy emission can be described by very few parameters and is well-approximated by a power law as explained below. The masses also impact the characteristic amplitude of a GW signal, which in turn may be linked to the detection threshold via the signal-to-noise ratio (SNR) of each signal. The latter sets the maximum expected redshift depth at which each type of signal may be observed by the detectors probing that frequency range. A good example of this is the background generated by double white dwarfs, which is expected to be observable by the LISA detector only at very low redshift [62] such that LISA will only be sensitive to white dwarfs from the Milky Way and, perhaps, to those in the local group [63].

Other astrophysical sources of GWs include asymmetric supernovae explosions and rotating neutron stars which are not spherically symmetric and may produce triaxial emission [64]. These are considered minor sources of GWs compared to compact binaries as their characteristic signal is weaker. Furthermore, they are also considerably harder to model; hence, the discussion that follows is mainly focused on the former; nevertheless, many equations presented are general and may be applied to these secondary sources as well. In principle, the fact that these are harder to model also implies that they will necessarily be targets of stochastic searches, as the latter may be considered as generic, “un-modelled” searches of GW power.

It may seem odd to use the fractional parameter ${\Omega }_{\mathrm{GW}}$ to describe GW energy from astrophysical events, as historically this is considered a cosmological quantity; however, several authors and collaborations have deemed it a useful convention allowing scientists from different fields to easily compare results. With this in mind, it may be considered as a compact measure of the GW history of the low-redshift Universe, after galaxy and star formation have taken place. To see this in more detail, following [40], let N be the number of events of a particular signal type in a comoving volume. This may be decomposed in event rate $\dot{N}$ multiplied by cosmic time slice $d{t}_r/dz$,

(27)$N\left(z\right)=\dot{N}\frac{d{t}_r}{dz}\equiv \frac{\dot{N}}{(1+z)H\left(z\right)}.$

Then, as in [65], the fractional energy density may be rewritten as

(28)${\Omega }_{\mathrm{GW}}=\frac{f}{{\rho }_{\mathrm{c}}}{\int }_0^{z_{\max }}dz\frac{\dot{N}\left(z\right)}{(1+z)H\left(z\right)}⟨\frac{d{E}_{\mathrm{GW}}}{d{f}_{\mathrm{r}}}{|}_{f_{\mathrm{r}}=f(1+z)}⟩,$

where the angle brackets indicate an averaging over particular source population parameters $\theta$, such that

(29)$⟨\frac{d{E}_{\mathrm{GW}}}{d{f}_{\mathrm{r}}}⟩=\int d\theta p\left(\theta \right)\frac{d{E}_{\mathrm{GW}}(\theta ;f_{\mathrm{r}})}{d{f}_{\mathrm{r}}},$

and $p\left(\theta \right)$ is the probability distribution of source parameters. In the case of a stochastic background, one would expect the parameter space to be fully sampled by the population; however, note that selection effects imposed by the detector may break this limit. One would expect the event rate $\dot{N}$ to locally increase with redshift, as the volume of each redshift shell increases. In fact, the event rate for stellar mass compact binary coalescences has been found to increase with redshift when analysing the direct detections reported by LIGO [66]. However, this rate cannot increase monotonically out to infinity and must have a turnaround point at some peak redshift $\widehat{z}$ and then decay to zero as star formation ceases in the early Universe. In general, this may be modelled phenomenologically with [67]

(30)$\dot{N}\left(z\right)=\mathcal{C}({\alpha }_{\mathrm{r}},{\beta }_{\mathrm{r}},\widehat{z})\frac{{\dot{N}}_0{(1+z)}^{{\alpha }_{\mathrm{r}}}}{1+{\left(\frac{1+z}{1+\widehat{z}}\right)}^{{\alpha }_{\mathrm{r}}+{\beta }_{\mathrm{r}}}},$

such that at low redshift, the rate increases as $\dot{N}\propto {(1+z)}^{{\alpha }_{\mathrm{r}}}$, then reaches a maximum at $z\sim \widehat{z}$, and finally decreases at high redshift as $\dot{N}\propto {(1+z)}^{-{\beta }_{\mathrm{r}}}$. $\mathcal{C}$ is the normalisation constant

(31)$\mathcal{C}({\alpha }_{\mathrm{r}},{\beta }_{\mathrm{r}},\widehat{z})=1+{(1+\widehat{z})}^{-({\alpha }_{\mathrm{r}}+{\beta }_{\mathrm{r}})}$

which ensures $\dot{N}(z=0)={\dot{N}}_0$. Constraining $\widehat{z},{\alpha }_{\mathrm{r}}$ and ${\beta }_{\mathrm{r}}$ provides insight on binary formation history and, consequently, star formation history as well. In practice, however, stochastic searches probe the redshift-integrated observable ${\Omega }_{GW}\left(f\right)$ only. The redshift history may be reconstructed with careful modelling and cross-correlation with direct measurements of individual events for which redshift is estimated; see the approach presented in [65]. Note that selection effects due to the mass distribution of sources need to be taken into account in this sort of analysis. Alternatively to this phenomenological approach, the redshift-dependent event rate may also be probed by conducting a careful study of the specific astrophysical conditions which favour compact binary inspirals and mergers; for example, the evolution of galaxy star formation and chemical enrichment can be combined with prescriptions for merging of compact binaries obtained via stellar evolution simulations [68,69].

To study the spectral dependence of different backgrounds, we can plug in different source models in Equation (18) and calculate the spectral shape of ${\Omega }_{\mathrm{GW}}$. This is particularly useful in the case of an astrophysical background as N and $E_{\mathrm{GW}}$ are directly related to the star formation history of the Universe, and astrophysical GWs from different source populations will contribute to ${\Omega }_{\mathrm{GW}}$ in different amounts at different times. This may result in a method that differentiates between astrophysical backgrounds and also test different stellar and galaxy evolution models [70].

An example of astrophysical background with a simple expected spectral dependence is that of compact binaries in the adiabatic inspiral phase [71]. In this case, it is sufficient to model GW radiation at leading quadrupole order [72],

(32)$\frac{d{E}_{\mathrm{GW}}}{d\ln f_{\mathrm{r}}}=\frac{{\pi }^{2/3}}{3}{G}^{2/3}{\mathcal{M}}^{5/3}{f}_{\mathrm{r}}^{2/3},$

where $\mathcal{M}$ is the chirp mass, given by

(33)$\mathcal{M}=\frac{{\left(m_1{m}_2\right)}^{3/5}}{{(m_1+m_2)}^{1/5}},$

where $m_1$ and $m_2$ are the masses of the two compact objects in the binary. The co-moving number density will depend on chirp masses in the binary distribution and may be rewritten as an integral over $\mathcal{M}$,

(34)$N={\int }_{{\mathcal{M}}_{\min }}^{{\mathcal{M}}_{\max }}\frac{dN}{d\mathcal{M}}d\mathcal{M},$

where the integration extremes will depend on source population characteristics. Rearranging the terms in Equation (28), one obtains

(35)${\Omega }_{\mathrm{GW}}^{\mathrm{CB}}\propto f^{2/3}{\int }_0^{z_{\max }}dz\frac{1}{{(1+z)}^{1/3}}{\int }_{{\mathcal{M}}_{\min }}^{{\mathcal{M}}_{\max }}d\mathcal{M}\frac{dN}{d\mathcal{M}}{\mathcal{M}}^{5/3},$

where the average over the population parameters is simply reduced to the integral over all possible chirp masses. Here, the superscript CB labels compact binaries generically. This provides a simple spectral dependence for ${\Omega }_{\mathrm{GW}}^{\mathrm{CB}}$, which may be condensed as⁵

(36)${\Omega }_{\mathrm{GW}}^{\mathrm{CB}}={\Omega }_{\mathrm{GW}}^{\mathrm{CB}}\left(f_{\mathrm{ref}}\right){\left(\frac{f}{f_{\mathrm{ref}}}\right)}^{2/3},$

where ${\Omega }_{\mathrm{GW}}^{\mathrm{CB}}\left(f_{\mathrm{ref}}\right)$ is the appropriately normalised energy density calculated at a reference frequency $f_{\mathrm{ref}}$.

This simple model may be extended to all general inspiral-dominated GWBs as the same principles apply and has been reprised all over the literature and used in almost every GWB search with data available. The same principle, for example, has been used in the estimation of the expected background from stellar mass black hole and neutron star binaries in the LIGO detector frequency band [73], with some tweaks to account for the chirp signal from merging black holes which has a substantially different energy output per frequency. This result is shown in Figure 2.

3.2. Primordial Backgrounds

There are a range of potential physical processes occurring in the primordial Universe which would result in a rich variety of GW signals. We report the most noteworthy here; however, this is a nonexhaustive list of all sources which have been studied in the literature. We focus especially on sources which have corresponding search pipelines and for which constraints have been derived with existing GW data. For a more in-depth review of the subject, see [76].

First and foremost, all inflationary models include an irreducible GWB as a byproduct of inflation [76,77]. Any inflationary process results in the stretching of quantum fluctuations which are parametrically amplified into classical density perturbations, producing particularly tensor perturbations which we identify as GWs. The amplification process causes perturbation wavelengths to expand until they reach the size of the horizon, when they exit causal contact and their evolution remains frozen until the Universe’s horizon returns to their scale in the post-inflationary era and they may re-enter. At time of re-entry, the GWB is formed by standing waves, and their phases are correlated between modes with opposite momenta [78]. In most models, the perturbations are adiabatic, Gaussian, and approximately scale-invariant, resulting in a similarly defined GWB. In this case, the fractional energy density will be approximately frequency independent at frequencies $f\gg H_0\sim {10}^{-18}\mathrm{Hz}$,

(37)${\Omega }_{\mathrm{GW}}^{\mathrm{Infl}.}\left(f\right)\approx \mathrm{const}.,$

and will be primarily characterised by its amplitude. The perturbations are also expected to imprint a pattern of B-modes in the polarization of the CMB [79], which has not been detected yet and is the aim of numerous CMB experiments [80,81,82,83]. This nondetection imposes an upper bound on the inflationary GWB amplitude, setting it well below the reach of present or planned GW detectors. The prevalent opinion in the community is that precision measurements of the CMB are still the best chance to garner decisive evidence of inflation and discriminate between models.

Beyond the irreducible background, however, the inflationary era may also source a second GWB which in certain cases dominates over the irreducible component, often presenting strongly scale-dependent features or characteristic peaks in the power spectrum [76]. For example, the presence of additional fields during the inflationary epoch other than the inflaton could have an effect on GW emission, as these may have interactions resulting in strong particle production or may act as spectator fields and exhibit a sub-luminal propagation speed [84,85]. The presence of additional symmetries in the inflationary sector would also have an impact, as these would result in breaking of space reparametrizations, allowing the graviton to acquire a mass. Alternative theories of gravity (other than general relativity) governing the inflationary era would have an impact on all aspects of primordial GWs. The scalar density perturbations generated during inflation may also source GWs; this may occur either at second order in perturbation theory [86] or by acting as seeds for primordial black holes [87,88]. Finally, GWs can also be produced during preheating [89,90,91,92] due to the violent transfer of energy from the inflationary sector to Standard Model plasma. All these possible inflationary scenarios are effectively extremely hard to probe with electromagnetic observations; hence, GWs present a unique window into this epoch. See for example [93] for a discussion of inflationary backgrounds in the LISA band and how a positive detection would revolutionise the current science of inflation. For more details and references to the specific models that may give rise to the effects discussed above, refer to the review [76].

After inflation, primordial GWs may have been generated as a consequence of phase transitions which may have occurred as the primordial Universe cooled down in its post-inflationary adiabatic expansion. First-order phase transitions in particular would proceed via the nucleation of true-vacuum bubbles, sourcing GWs through the collision of these bubbles [44,46], as well as through the resulting acoustic and/or turbulent motion of the primordial plasma [55]. The characteristic power spectrum of a background formed by bubble collisions presents a sharp peak at the frequency related to the energy scale at which the transition occurred, $T_{*}$,

(38)$f_{\mathrm{peak}}\sim {10}^{-6}\mathrm{Hz}\times \frac{\beta }{H_{*}}\times \frac{T_{*}}{200\mathrm{GeV}},$

where $\beta /H_{*}$ is the rate of the transition in units of the Hubble rate at the transition epoch. For typical values $\beta /H_{*}\sim {10}^3$, we see that this peak lies in the megahertz (mHz) band probed by LISA for transitions at the electroweak scale $T_{*}\sim 200\mathrm{GeV}$, allowing us to test a myriad of particle physics models which predict such signals.

Phase transitions may also give rise to stable topological defects, the most prominent example of which are cosmic strings. These line-like topological defects are generated through spontaneous symmetry breaking in the early Universe [47,48] and are a generic prediction of many theories beyond the Standard Model of particle physics [94]. On macroscopic scales, cosmic strings are effectively described by a single parameter, the dimensionless string tension $G\mu$, which is determined by the energy scale at which they were formed. Once formed, dynamical evolution results in a dense network of cosmic string loops which oscillate at relativistic speeds, producing copious GWs which combine to produce a strong SGWB signal for which its amplitude and spectral shape are set by $G\mu$.

As mentioned above, a potential GW signal of cosmological origin is the SGWB from primordial black holes (PBHs); i.e., black holes formed in the early Universe rather than through stellar evolution. Since black holes are massive and non-baryonic and interact only through gravity, PBHs are very natural and well-motivated cold dark matter (CDM) candidates. The cosmic mass density of PBHs is, therefore, usually expressed as a fraction of the total CDM mass density, $f_{\mathrm{PBH}}\equiv {\rho }_{\mathrm{PBH}}/{\rho }_{\mathrm{CDM}}$. Similarly to stellar-origin black holes, PBHs can form binaries (either by forming in close proximity to each other or through dynamical encounters) which then inspiral and emit GWs or they can emit GWs through close hyperbolic encounters [95]. The cumulative signal from many such binaries produces a SGWB spectrum which is determined by the PBH mass spectrum and DM mass fraction $f_{\mathrm{PBH}}$.

In general, GWBs emitted in the primordial Universe are considered stochastic under the ergodic assumption, which equates the ensemble average of an observable with its time and/or spatial average. Essentially, we assume that by observing large enough regions of the Universe today, or a given region for long enough time, we have access to many independent realisations of the system. This is realistic given two fundamental premises of cosmology, namely that the Universe is statistically homogeneous and isotropic, such that the initial conditions of the GW-generating process may be considered the same at every point in space and that the GW source fulfills causality, operating at a time when the typical size of a region of causal contact in the Universe was smaller than today. It can be shown [76] that the correlation scale today of a GW signal from the early Universe is tiny comparable to the present Hubble scale, verifying this requirement.

Stochasticity would also imply the lack of a preferred polarisation mode in the signal, as all mode configurations should be fully sampled in the observation set. While this is intuitively true for astrophysical backgrounds, the argument for primordial GWs is more delicate, and while a background may be stochastic for the reasons described above, there may be specific formation channels which will violate parity and induce an asymmetry in the polarisation modes, effectively generating circularly polarised waves. This effect may be found in alternative theories of gravity, such as Chern–Simons gravity [96] or certain quantum gravity models [97]. The study of parity-violating theories of gravity is inspired by grand unification, as the electro-weak sector of the standard model is chiral and, thus, maximally violates parity [98,99].

3.3. Anisotropies in Stochastic Backgrounds

We can distinguish two root causes of anisotropy in stochastic signals. Firstly, there may be inhomogeneities in the GW source mechanisms, such as a particular distribution of the sources on the sky, which necessarily produces an anisotropic signal. Secondly, as GWs propagate, they accumulate line-of-sight effects [100], crossing different matter density fields which are (if mildly) inhomoheneously distributed in our Universe. For astrophysical backgrounds, the spatial distribution of GW events is a biased tracer of the underlying light matter distribution, while primordial backgrounds will present anisotropies linked to the particular physical processes that occurred in the early Universe. In order to investigate these effects, it is useful to define the fractional GW energy density per solid angle $d\widehat{\mathit{n}}$,

(39)${\Omega }_{\mathrm{GW}}(f,\widehat{\mathit{n}})=\frac{1}{{\rho }_c}\frac{d{\rho }_{\mathrm{GW}}}{dlnfd\widehat{\mathit{n}}},$

such that

(40)${\overline{\Omega }}_{\mathrm{GW}}\left(f\right)=\frac{1}{4\pi }{\int }_{S^2}d\widehat{\mathit{n}}{\Omega }_{\mathrm{GW}}(f,\widehat{\mathit{n}}),$

as in [101,102]. ${\overline{\Omega }}_{\mathrm{GW}}\left(f\right)$ is the GWB monopole, where the bar is introduced to highlight the fact that it is an isotropic all-sky average over the directional fluctuations. The latter may be described by the GW density contrast,

(41)${\delta }_{\mathrm{GW}}(f,\widehat{\mathit{n}})=\frac{{\Omega }_{\mathrm{GW}}(f,\widehat{\mathit{n}})-{\overline{\Omega }}_{\mathrm{GW}}\left(f\right)}{{\overline{\Omega }}_{\mathrm{GW}}\left(f\right)},$

in analogy with the formalism used to quantify the level of anisotropy in CMB temperature [103] and in galaxy clustering [104]. This dimensionless parameter provides an unambiguous metric to compare different results.

At first order in the perturbations, the density contrast may be decomposed into three distinct contributions related to the origin of the anisotropies,

(42)${\delta }_{\mathrm{GW}}\simeq {\delta }_{\mathrm{GW}}^{\mathrm{s}}+{\delta }_{\mathrm{GW}}^{\mathrm{los}}+\mathcal{D}\widehat{\mathit{n}}·{\widehat{\mathit{v}}}_{\mathrm{obs}}.$

Here, ${\delta }_{\mathrm{GW}}^{\mathrm{s}}$ is the source term resulting from anisotropies at the GW source, for example, due to the specific distribution of events for astrophysical backgrounds, while ${\delta }_{\mathrm{GW}}^{\mathrm{los}}$ is the propagation term [100], which encapsulates accumulated line-of-sight effects. The final term is the dipole of magnitude $\mathcal{D}$ induced by the observer’s peculiar velocity ${\mathit{v}}_{\mathrm{obs}}$. Note that the propagation term includes contributions similar to those encountered in CMB anisotropy calculations: the Sachs–Wolfe effect, due to the local curvature at emission, and the integrated Sachs–Wolfe effect, due to the curvature perturbations encountered along the line-of-sight [105]; a Doppler term due to the source’s peculiar velocity; and higher order effects such as gravitational lensing and redshift-space distortions [106]. All these terms have been the subject of intense study in recent years and both astrophysical and primordial contributions have been modelled extensively (see for example [107,108,109] for independent calculations of astrophysical background anisotropies and [50,102,110] for examples of primordial ones). The propagation term is often assumed to be negligible in the case of astrophysical backgrounds or it is directly included as part of the source term; in [111,112], it is expressly estimated and shown to account for at most ∼10% of the total anisotropy. In the case of a truly primordial background dating back as far as or before the CMB, however, this term will be considerably larger and will be comparable to the intrinsic anisotropies at large scales as shown in [100].

We analyse here the general scenario where the density contrast at the source behaves as a Gaussian random field, as this is mostly assumed in the literature as well. In this case, ${\delta }_{\mathrm{GW}}^{\mathrm{s}}$ is characterised by its two-point correlation function,

(43)$C_{\mathrm{GW}}(f,cos\theta )=⟨{\delta }_{\mathrm{GW}}^{\mathrm{s}}(f,\widehat{\mathit{n}}){\delta }_{\mathrm{GW}}^{\mathrm{s}}(f,{\widehat{\mathit{n}}}^{\prime })⟩,$

in analogy with CMB anisotropies, as the first moment $⟨{\delta }_{\mathrm{GW}}^{\mathrm{s}}⟩$ is zero by definition, and all higher order moments will either vanish or may be expressed in terms of $C_{\mathrm{GW}}$ [102]. Here, $cos\theta =\widehat{\mathit{n}}·{\widehat{\mathit{n}}}^{\prime }$ and the angular brackets represent an averaging over all $\widehat{\mathit{n}},{\widehat{\mathit{n}}}^{\prime }$ direction pairs of aperture $\theta$. $C_{\mathrm{GW}}(f,cos\theta )$ will then be different from zero where there are consistent spatial correlations on the sky associated to $\theta$, at frequency f. To observe this formally, we expand the two-point function in spherical harmonic space as

(44)$C_{\mathrm{GW}}(f,cos\theta )=\sum _{\ell =0}^{\infty }\frac{2\ell +1}{4\pi }{C}_{\ell }^{\delta }\left(f\right)P_{\ell }(cos\theta ),$

where $P_{\ell }$ is the ℓth Legendre polynomial. By inverting this equation, we obtain the angular power spectrum $C_{\ell }^{\delta }\left(f\right)$ for the GW energy density contrast,

(45)$C_{\ell }^{\delta }\left(f\right)=2\pi {\int }_{-1}^{+1}dcos\theta C_{\mathrm{GW}}(f,cos\theta )P_{\ell }(cos\theta ),$

which may be estimated from observations. Note that quantity $\ell (\ell +1)C_{\ell }/2\pi$ is roughly the contribution per logarithmic ℓ bin to the variance of ${\delta }^{\mathrm{s}}$; hence, this will often be the quantity shown in plots quantifying the $C_{\ell }^{\delta }$ power spectrum from different signal models and in GW data sets.

Alternatively, we can also calculate the angular power spectrum of the ${\Omega }_{\mathrm{GW}}$ parameter itself, first by decomposing it into spherical harmonic components,

(46)${\Omega }_{\mathrm{GW}}(f,\widehat{\mathit{n}})=\sum _{\ell =0}^{\infty }\sum _{m=-\ell }^{+\ell }{\Omega }_{\ell m}^{\mathrm{GW}}\left(f\right)Y_{\ell m}\left(\widehat{\mathit{n}}\right),$

and then by defining

(47)$C_{\ell }^{\Omega }=\frac{1}{2\ell +1}\sum _{m=-\ell }^{+\ell }⟨{\Omega }_{\ell m}^{\mathrm{GW}}\left(f\right){\Omega }_{\ell m}^{\mathrm{GW}}\left(f\right)⟩,$

where angle brackets indicate an ensemble average over random realizations of ${\Omega }_{\mathrm{GW}}$. The information encapsulated in $C_{\ell }^{\Omega }$ and $C_{\ell }^{\delta }$ is the same, and they are both dimensionless quantities; however, note that the numerical scale of $C_{\ell }^{\Omega }$ is roughly that of ${\Omega }_{\mathrm{GW}}^2$, whereas $C_{\ell }^{\delta }$ shows the fractional contribution of each mode to the power spectrum with respect to the monopole ${\overline{\Omega }}_{\mathrm{GW}}$.

To quote the expected levels of anisotropy for astrophysical backgrounds from calculations in the references mentioned above, one may compare Figure 2 in [107] and Figure 4 in [108]. It has been pointed out in the literature [102,113,114,115] that there is a tension between these two predictions; the difference in the estimates highlighted here stems probably from different treatments of the source effects at very low redshift, where the clustering at small scales induces non-linear effects in the field which the overall anisotropy level is very sensitive to. Cosmological backgrounds are expected to present CMB-like anisotropies of ${\delta }_{\mathrm{GW}}^{\mathrm{s}}\sim {10}^{-5}$ [100]. A comparison between these estimates and the expected sensitivity-per-multipole of various GW detectors may be found in [116].

Finally, note that to properly assess the sourced GW anisotropies from direct measurements, it is necessary to subtract out the kinematic dipole term induced by the observer’s peculiar velocity as seen in Equation (42). This may be performed similarly to CMB analysis, potentially using the independent estimate of ${\mathit{v}}_{\mathrm{obs}}$ from CMB data themselves.

3.4. Observational Properties of Stochastic Backgrounds

The observational properties of the SGWB such as time-domain characteristics and selection effects imposed by the observation method merit special attention as these are crucial for developing effective search methods. For the sake of simplicity, this discussion is focused on astrophysical backgrounds where individual events are clearly defined and are not extremely model dependent. However, most of the following considerations generally apply to primordial backgrounds as well and have intuitive extension to backgrounds from topological defects and primordial black holes for example.

Three main regimes are identified in the literature based on the frequency of occurrence of the events in time. These are referred to as Poisson noise regime, where GWs appear as isolated bursts in the timestream, Gaussian regime, where the timestream is packed with numerous overlapping signals and the central limit theorem applies, and popcorn or intermittent regime, which lies somewhere in between. A visualisation of these is provided in Figure 3, generated using mock-data generation open source package presented in [117,118]. The discriminating factor between regimes is encapsulated in the duty cycle, which is the mean ratio between event duration and the time interval between consecutive events. This will also be equal to the average number of events present at the detector at a given observation time. While it is intuitively clear that the Gaussian background in Figure 3 is stochastic, this is not obvious of the other regimes. Operationally, a signal may be considered stochastic if it is more effective, in the Bayesian sense of the word, to search for it in the data using a stochastic model for the waveform as opposed to a deterministic one [119]. Furthermore, a signal is defined as resolvable if it can be decomposed into non-overlapping and individually detectable signals, in either the time or frequency domain. With this in mind, the only way to know whether a signal is stochastic, is to search for it and see what method yields the best result or more stringent upper limit. In the case of astrophysical GW signals, there is a continuous distribution of events where the highest SNRs are expected to be resolved, and beneath a certain detector-imposed cutoff, the events accumulate to form the background signal. The boundary between the resolved and unresolved population is not obvious, and a few bright signals stand out above the confusion background. When these are resolvable, deterministic signals, they should be subtracted from data, leaving a residual background that may be analysed with stochastic tools. In [119], the authors present simulations designed to establish whether certain signals are best searched for with stochastic or deterministic methods. They find, similarly to [120] and somewhat unsurprisingly, that deterministic signal models are favoured whenever the background is made up of a small number of sources, and stochastic models are preferred for larger source numbers. The boundary between regimes however is not well defined, and hybrid models comprised of a deterministic signal model for individual bright sources and a Gaussian–stochastic signal model for the confusion background are generally most effective to extract the population parameters.

These analyses highlight the distinction between event rate and detection rate. The event rate is the actual rate at which binaries are merging in the visible Universe, while the detection rate is the rate with which we can detect them, which is plagued by selection effects depending on both the binary’s intrinsic and extrinsic parameters and the detector’s characteristics.

There is a further degree of complication caused by these unresolved, bright signals in the timestream which behave as shot noise and may hinder the correct interpretation of a detection. These signals represent random, Poissonian fluctuations in the finite population of sources and will induce a white noise component in the angular power spectrum. The observational effect that arises is induced by the fact that nearby, bright sources approximately randomly distributed around the observer will mask the underlying structure of the Universe. This is analysed and quantified in [121,122,123] for the monopole and anisotropies of an astrophysical background, respectively. Here, Poisson white noise is modelled with a free parameter representing the cutoff distance from Earth beyond which the GW signal may be considered a stochastic background by the definition given above. The cutoff naturally depends on the sensitivity of the detector array under consideration. This is analysed in [122,123] for the anisotropic background, and it is found that the white noise component dominates current LIGO and Virgo detectors; however, the planned upgrades of these detectors and the future 3G detectors may be sufficient to curb this effect. As discussed in [124,125], another strategy to mitigate the shot noise is to cross-correlate the GW signal with a galaxy catalogue, as this cross-correlation spectrum has a much lower level of shot noise. Furthermore, as GW sources are finite in time, then Poisson statistics dictate that this noise decays as the inverse of the observation time; hence, long observing runs will progressively mitigate this effect. Quantifying this distinction is an evergrowing effort. The current event rate estimates based on astrophysical models informed by the LIGO detections show that there is approximately one binary neutron star merger every 13 s and one binary black hole merger every 223 s [73]. Clearly, the vast majority of these lie well beneath the threshold for detection. This is quantified for example in [126] where the detection rate per year is encoded as a function of the merging binary masses. It is shown that there are several astrophysical models for compact binaries which fit the detection rate of ∼10 events per year from the aLIGO O1 and O2 data runs. It is clear that the rate of “intermediate” subthreshold events from an astrophysical population following these models is low and remains far from the stochastic limit as defined above.

It is worth pointing out that these observational features of astrophysical backgrounds discussed here may themselves provide insight in compact binary population parameters such as binary coalescence rates for different binary species, as detailed in [127]. Furthermore, the time dependence of astrophysical backgrounds in general may be a powerful discriminating factor between these and backgrounds of primordial origin.

4. Detection Approaches and Methodologies

Given its stochastic nature, the main challenge in measuring an SGWB is confidently distinguishing signal from the (stochastic) noise in GW detectors. The interplay between signal and noise bears a central role in stochastic search methods, and in particular the validity of a given method may depend on the expected signal-to-noise regime or ratio (SNR): whether the noise dominates the measurement, i.e., the detector operates in the low signal regime, or whether the signal is comparable to or stronger than the detector noise, i.e., high signal regime. Here, we classify search methods based on assumptions made primarily on the signal, such as Gaussianity and isotropy, starting from minimal assumptions on SNR, and we highlight how the latter then impacts the formulation of the optimal estimators.

Cross-correlation searches have cast a wide net, targeting a variety of stochastic backgrounds within the same style of search and also refining methods to extract specific signals from the data—see, e.g., the dark photon search (e.g., [128]) or the cosmic string search (e.g., [129]). We divide searches into three broad categories. We present the methodology to search for an isotropic stochastic background, which relies on the isotropic assumption to simplify the estimator considerably. Then, we show what happens when this assumption is relaxed in searches for anisotropic backgrounds. Finally, we review special searches which target elusive signatures which may imply deviation from general relativity and/or the standard model. However, first, let us review the common denominators in stochastic searches, namely the cross-correlation statistic and the correlated detector response.

4.1. Interdetector and Spatial Correlations

The first step in most stochastic searches is to consider the auto- and cross-correlation of datastreams from one or more detectors, as this operation discards all phase information and targets the signal power directly. To then deconvolve the response of the detectors to the signal, many searches rely on an overlap function, which determines the cross-power response of a pair of detectors to a stochastic background. The most common approach is to estimate the significance of excess correlated power in multiple detectors, integrated over all available observing time, often assuming noise in different detectors not to be correlated⁶. Throughout this section, these methodologies will be reviewed in detail; first, however, let us formally introduce the cross-correlation statistic and derive the overlap functions for different detector setups as these will be the building blocks for the rest of this review.

To start, consider a Gaussian, stationary GWB observed through two distinct detectors, 1 and 2. Each detector timestream may be decomposed into signal and noise components as

(48)$d_1\left(t\right)=R_1\left(t\right)\ast h\left(t\right)+n_1\left(t\right),d_2\left(t\right)=R_2\left(t\right)\ast h\left(t\right)+n_2\left(t\right),$

where $R_1$ and $R_2$ represent the individual detector response functions to GW signals, and “*” represents a convolution. These are characterised by functional dependencies on intrinsic and extrinsic parameters of the GWs such as polarisation and incoming direction; however, we omit these now as they may be easily included at a later stage.

We consider both h and n fields in Equation (48) as mean zero Gaussian variates; hence, these are fully described by their second-order moments. Let us work in the Fourier domain, where the assumption of stationarity implies that the noise covariance is diagonal, i.e., the noise in each frequency bin is uncorrelated⁷,

(49)$⟨n_i\left(f\right)n_i\left(f^{\prime }\right)⟩={\delta }_{f{f}^{\prime }}{P}_{f,i},$

with $P_i$ as the noise power spectrum in detector i. We define the discrete Fourier transforms of the data streams starting at timestamp $\tau$ as ${\tilde{d}}_1^{\tau }$, ${\tilde{d}}_2^{\tau }$, and drop the time label immediately for the sake of conciseness. We then take the cross-correlation of these quantities as

(50)$C_{12}\left(f\right)={\tilde{d}}_1^{}\left(f\right){\tilde{d}}_2^{\star }\left(f\right),$

which, assuming uncorrelated noise between detectors, has expectation value of

(51)$⟨C_{12}\left(f\right)⟩=R_1^{}\left(f\right)R_2^{\star }\left(f\right)⟨{\tilde{h}}_1^{}\left(f\right){\tilde{h}}_2^{\star }\left(f\right)⟩=T_{\mathrm{obs}}{\Gamma }_{12}\left(f\right)I_{\mathrm{GW}}\left(f\right).$

${\Gamma }_{12}$ is the overlap reduction function, where the reduction comes from the fact that it is a sky-integrated quantity (more details to follow), and $I_{\mathrm{GW}}\left(f\right)$ is the sky-integrated power spectrum of the background. Note the normalisation $T_{\mathrm{obs}}$ as we assume the observation lasts a finite amount of time. The covariance of the correlated measurement is as

(52)$\mathrm{Cov}\left[C_{12f,f^{\prime }}^{}\right]=⟨C_{12}\left(f\right)C_{12}^{\star }\left(f^{\prime }\right)⟩-⟨C_{12}^{}\left(f\right)⟩⟨C_{12}^{\star }\left(f^{\prime }\right)⟩\approx T_{\mathrm{obs}}{\delta }_{f{f}^{\prime }}{P}_{1,f}{P}_{2,f},$

taking the weak signal limit in the final equality. In this case, the cross-correlation statistic has been shown to be a sufficient statistic [131], making it a good choice for stochastic searches. This is because there is little information to be drawn out from the auto-correlated measurements, as these are effectively measurements of the detector noise power; conversely, to use auto-correlations to estimate the signal, one would need an independent and accurate measurement of the noise power to model and subtract out the noise. This approach would require an iterative estimator where both signal and noise are included in the fit [132]. Note that using the cross-correlation component only of the data and, thus, ignoring the auto-correlations allows for considerable simplification of the optimal estimators and computations necessary in data analysis.

The functional form of the overlap function of a detector network can be derived from Equation (51) by appropriately combining the individual response functions of each detector. In the case of laser interferometers, the basic response function is the arm response function, and arms are then combined to form the detector. This may be performed physically, by using laser interferometry as is performed in LIGO/Virgo [133], or in post-processing, with time-delayed interferometry (TDI) as in the case of LISA [134,135]. Considering arm ${\mathit{u}}_{PS}$ spanning from point P to point S, the one-way arm response $\overrightarrow{u}:P\to S$ in the frequency domain is described as [33]

(53)$R_{\overrightarrow{u}}^A(f,\widehat{\mathit{n}})=\frac{1}{2}\mathit{u}\otimes \mathit{u}:{\mathbf{ϵ}}^A{\mathcal{T}}_{\overrightarrow{u}}(f,\widehat{\mathit{n}}·\mathit{u})e^{i\frac{2\pi f}{c}\widehat{\mathit{n}}·{\mathit{x}}_S},$

with the one-way timing transfer function relative to a photon travelling along ${\mathit{u}}_{PS}$, ${\mathcal{T}}_{\mathit{u}}$, as

(54)${\mathcal{T}}_{\overrightarrow{u}}(f,\widehat{\mathit{n}}·\mathit{u})=\frac{L}{c}{e}^{-\frac{\pi fL}{c}(1+\widehat{\mathit{n}}·\widehat{\mathit{u}})}\sin \mathrm{c}\left(\frac{\pi fL}{c}(1+\widehat{\mathit{n}}·\widehat{\mathit{u}})\right).$

In the case of LIGO-like detectors, which are ${90}^{\circ }$-arm, Michelson–Morley style interferometers, a laser photon performs a round trip along one arm before interfering with the beam in the other arm; hence, the relevant response function is the two-way arm response function along $\overleftrightarrow{u}:P\leftrightarrow S$, which we can write in terms of the one-way timing transfer function as

(55)$R_{\overleftrightarrow{u}}^A(f,\widehat{\mathit{n}})=\frac{1}{2}\mathit{u}\otimes \mathit{u}:{\mathbf{ϵ}}^A({\mathcal{T}}_{\overrightarrow{u}}(f,\widehat{\mathit{n}}·\mathit{u})e^{i\frac{2\pi f}{c}(\widehat{\mathit{n}}·{\mathit{x}}_S-L)}+{\mathcal{T}}_{\overleftarrow{u}}(f,\widehat{\mathit{n}}·\mathit{u})e^{i\frac{2\pi f}{c}\widehat{\mathit{n}}·{\mathit{x}}_P});$

note the extra $-L$ term which accounts for the extra phaseshift incurred when travelling over the same distance twice. Then, assuming that the detector is not moving during the measurement, we can substitute in ${\mathit{x}}_S={\mathit{x}}_P+\mathit{u}$, and rewrite the above as

(56)$R_{\overleftrightarrow{u}}^A(f,\widehat{\mathit{n}})=\frac{1}{2}\mathit{u}\otimes \mathit{u}:{\mathbf{ϵ}}^A{\mathcal{T}}_{\overleftrightarrow{u}}(f,\widehat{\mathit{n}}·\mathit{u})e^{i\frac{2\pi f}{c}\widehat{\mathit{n}}·{\mathit{x}}_P},$

where the two-way timing transfer function along the same arm is derived as

(57)$\begin{array}{cc}\hfill {\mathcal{T}}_{\overleftrightarrow{u}}(f,\widehat{\mathit{n}}·\mathit{u})& =\frac{L}{c}{e}^{-i\frac{2\pi fL}{c}}[e^{-i\frac{\pi fL}{c}(1-\widehat{\mathit{n}}·\widehat{\mathit{u}})}\mathrm{sinc}\left(\frac{\pi fL}{c}(1+\widehat{\mathit{n}}·\widehat{\mathit{u}})\right)\hfill \\ \hfill & +e^{i\frac{2\pi fL}{c}(1+\widehat{\mathit{n}}·\widehat{\mathit{u}})}\mathrm{sinc}\left(\frac{\pi fL}{c}(1-\widehat{\mathit{n}}·\widehat{\mathit{u}})\right)].\hfill \end{array}$

Note that in the small antenna limit, i.e., when the scale of the detector is much smaller than the wavelength probed, $L\ll c/f_{\mathrm{GW}}$, we have that $|{\mathcal{T}}_{\overleftrightarrow{u}}|⟶1$, implying that all the frequency dependence of the detector response is encoded in the phase term of Equation (56). Consider now detector 1 which is made up of arms $\mathit{u}$ and $\mathit{v}$, joined at point P; its polarisation response in frequency space, in the small antenna limit, will simply be

(58)$F_1^A(f,\widehat{\mathit{n}})=R_{\overleftrightarrow{u}}^A(f,\widehat{\mathit{n}})-R_{\overleftrightarrow{v}}^A(f,\widehat{\mathit{n}})=\frac{1}{2}({\mathit{u}}_1\otimes {\mathit{u}}_1-{\mathit{v}}_1\otimes {\mathit{v}}_1):{\mathbf{ϵ}}^A{e}^{i\frac{2\pi f}{c}\widehat{\mathit{n}}·{\mathit{x}}_1},$

where we have substituted in ${\mathit{x}}_1$ for ${\mathit{x}}_P$ as we consider it to be the location of the detector itself. Finally, the unpolarised overlap function⁸ of a detector pair is given by combining the response functions of the two detectors [33],

(59)${\gamma }_{12}(f,\widehat{\mathit{n}})=F_1^{+}(f,\widehat{\mathit{n}})F_2^{+}(f,\widehat{\mathit{n}})+F_1^{\times }(f,\widehat{\mathit{n}})F_2^{\times }(f,\widehat{\mathit{n}}),$

while the overlap reduction function $\Gamma$ is obtained by sky-integrating the above,

(60)${\Gamma }_{12}\left(f\right)={\int }_{S^2}d\widehat{\mathit{n}}{\gamma }_{12}(f,\widehat{\mathit{n}}).$

Note that polarisation F functions here defined specifically for ground-based interferometers correspond to the individual detector response functions seen in Equation (51). We have chosen to incorporate the frequency-dependent phase terms into the response functions directly; however, in the literature, these are sometimes left out of F functions to underline the fact that, in the small antenna limit, there is no interaction between GW frequency and detector geometry [33]. In the overlap function, the phase terms add, giving rise to a time delay term representing the frequency-dependent phase shift between the two detectors; explicitly,

(61)${\gamma }_{12}(f,\widehat{\mathit{n}})={\overline{\gamma }}_{12}\left(\widehat{\mathit{n}}\right)e^{i\frac{2\pi f}{c}\widehat{\mathit{n}}·\mathit{b}},$

where $\overline{\gamma }$ is the geometric component of the overlap function, and $\mathit{b}$ is the baseline vector that points from detector 1 to detector 2. The numerical values of $\overline{\gamma }$ as a function of direction then depend specifically on the coordinate basis chosen to represent the function; often, the choice is to represent it in an Earth-centered celestial frame, where the function will be time-dependent, as the arms sweep round and round on the sky following the Earth’s rotation on its axis. This behaviour gives rise to the specific scan strategy of the detector array, which is essential for directional stochastic searches. On the other hand, the frequency-dependent phase term highlights the role of the baseline length in a stochastic search: the ratio between $\left|\mathit{b}\right|$ and GW wavelength will determine what frequencies will be suppressed in the search; the longer the baseline, the higher the frequencies accessible will be. The overlap reduction functions for the LIGO Hanford and Livingston detector pair (L–H), LIGO Hanford and Virgo (H–V), and LIGO Hanford and Kagra (L–K) are shown in Figure 4a, while the respective (instantaneous) overlap functions are plotted on the sky in Figure 5 at a particular common time.

As mentioned above, long arm length, spaceborne detectors such as the future LISA three-satellite constellation, oriented in roughly an equilateral triangle, will measure GWs by combining the one-way arm responses (53) of several arms at different times composing TDI configurations or channels⁹. The most basic configuration can be thought of as a TDI–Michelson–Morley interferometer, where one constructs a response such as that in Equation (58). However, the great advantage of TDI is that it is possible to construct channels that suppress detector noise and enhance the signal [138]; this is a very active area of study [139,140,141,142,143,144]. For the purposes of this review, let us cite the TDI $\{X,Y,Z\}$ set of channels, which are based on the Michelson–Morley style measurement in Equation (58), where nominally X is centered around spacecraft 1, referred to hereafter as $S_1$; Y around spacecraft 2, $S_2$; and Z around spacecraft 3, $S_3$. However, to minimise noise, each light path is “reflected” back and forth between spacecrafts multiple times such that the noise interferes destructively. The total number of reflections is then linked to a version number for the TDI channel; e.g., “TDI 1 X” refers to a channel constructed by “interfering” the following light paths:

(62)$(S_1\to S_2\to S_1\to S_3\to S_1)-(S_1\to S_3\to S_1\to S_2\to S_1).$

For more details on this calculation, see also official LISA documentation [145]. Another useful set of channels is the orthonormal set which approximately diagonalise the $X,Y,Z$ noise-correlation matrix, $A,E,T$ [146]. This set is often used for simplicity as one may consider the auto-correlations of A and E only, treating them as separate measurements, and use T as a Sagnac channel which is mostly sensitive to detector noise and blind to the signal [147,148]. However, in practice diagonalizing $X,Y,Z$ is not a simple task as the detectors continuously move with respect to the measurement. The auto- and cross-correlated overlap functions for LISA are constructed as in Equation (59), where now the polarisation responses are referred to TDI channels. Note that LISA does not operate in the small antenna limit, hence the timing transfer function remains frequency dependent, inducing oscillations in TDI responses at high frequency, as observed in Figure 4b.

In the case of pulsar timing arrays, we can use the one-way transfer function defined in Equation (54), as in this case the arm vector $\widehat{\mathit{u}}$, points from the solar system barycenter to the pulsar, and the arm length L is the distance from the solar system barycenter to the pulsar. One can expand the sinc function using Euler’s identity and isolate the exponentials with arguments that depend on $fL/c$. Those exponential terms vary rapidly with sky direction for $fL/c\gg 1$ averaging down to zero in Equation (60) for Earth–pulsar baselines greater than 100 parsec in the nanohertz band (e.g., Figure 1 in [149]). From here, it is straightforward to calculate the (now frequency-independent) overlap reduction function by considering the cross-correlation of data from two pulsars, which serve as one-way timing beams (see, e.g., [33,149,150] for a more complete discussion). The overlap reduction function for a pair of pulsars a and b now depends on the angle between Earth–pulsar baselines, ${\zeta }_{ab}$ [150]:

(63)${\Gamma }_{\mathrm{PTA}}\left({\zeta }_{ab}\right)=\frac{1}{2}-\frac{1}{4}\left(\frac{1-cos{\zeta }_{ab}}{2}\right)+\frac{3}{2}\left(\frac{1-cos{\zeta }_{ab}}{2}\right)ln\left(\frac{1-cos{\zeta }_{ab}}{2}\right).$

where ${\Gamma }_{\mathrm{PTA}}$ is commonly referred to as the Hellings–Downs curve and may be observed in Figure 6.

4.2. Isotropic Background Search Methods

In this section, we review the case of isotropic backgrounds searches. We start by presenting the approach for ground-based interferometers, which may be considered the simplest case as these operate in the low-signal regime. We then lay out the detection strategy of PTAs. The application of these methods on real data and corresponding search results are reviewed in Section 5.

4.2.1. Ground-Based Detectors

The flagship LVK stochastic search targets a Gaussian, isotropic, stationary, and unpolarised background [56]. This set of assumptions results in a simple, ready-to-use estimator which results in a relatively rapid analysis of LIGO-Virgo data. All other collaboration searches use, as a starting point, the breakdown of one of these assumptions, and most extend the estimator to a particular case.

As detailed in Section 3, a stochastic signal which satisfies all assumptions listed above is fully described by its spectral energy density, which may be expressed in terms of the fractional energy density ${\Omega }_{\mathrm{GW}}\left(f\right)$. Under the assumption the latter presents a power-law spectrum as in Equation (25), we can use the cross-correlation statistic presented in Equations (51) and (52) to write down an optimal filter to apply to the data and draw out an estimate for the amplitude ${\Omega }_{\mathrm{GW}}\left(f_{\mathrm{ref}}\right)$ at reference frequency $f_{\mathrm{ref}}$. To start, we remind the reader of Equation (20) relating the measured GWB intensity to the energy density, which we may rewritten as

(64)${\Omega }_{\mathrm{GW}}\left(f_{\mathrm{ref}}\right)=\mathcal{H}(f,f_{\mathrm{ref}})I\left(f\right),$

with

(65)$\mathcal{H}(f,f_{\mathrm{ref}})=\frac{4{\pi }^2}{{\rho }_{c}G}{f}_{\mathrm{ref}}^3{\left(\frac{f}{f_{\mathrm{ref}}}\right)}^{3-\alpha }.$

By inspection, we already oberve that $\mathcal{H}(f,f_{\mathrm{ref}})$ is an unbiased filter to use on the data to derive an estimate of the background, in the case where the GWB obeys Equations (20) and (51). In general, there are several methods to show that this filter is optimal and write an estimator for ${\Omega }_{\mathrm{GW}}$, and different assumptions on the signal and noise components result in distinct solutions. We report the derivation shown in [151]. The approach here is to find the function which maximises the SNR of the estimated GWB signal, defined as

(66)$\mathrm{SNR}=\frac{{\widehat{\Omega }}_{\mathrm{GW}}{|}_{f_{\mathrm{ref}}}}{{\sigma }_{\Omega }^2},$

where

(67)${\sigma }_{\Omega }^2=\mathrm{Cov}\left[{\widehat{\Omega }}_{\mathrm{GW}}{|}_{f_{\mathrm{ref}}}\right].$

Re-writing the expectation value for the data correlation $C\left(f\right)$ in terms of ${\Omega }_{\mathrm{GW}}\left(f_{\mathrm{ref}}\right)$ using Equations (20) and (51),

(68)$⟨C\left(f\right)⟩=T_{\mathrm{obs}}{\Gamma }_{12}\left(f\right){\mathcal{H}}^{-1}(f,f_{\mathrm{ref}}){\Omega }_{\mathrm{GW}}\left(f_{\mathrm{ref}}\right),$

and assuming the low-signal limit following Allen and Romano [151], one obtains

(69)${\widehat{\Omega }}_{\mathrm{GW}}{|}_{f_{\mathrm{ref}}}=F_{f_{\mathrm{ref}}}^{-1}{\int }_0^{\infty }df\frac{{\Gamma }_{12}\left(f\right)\mathcal{H}(f,f_{\mathrm{ref}})}{P_1\left(f\right)P_2\left(f\right)}C\left(f\right),$

where $F_{f_{\mathrm{ref}}}$ may be thought of as a normalisation which depends on the chosen reference frequency,

(70)$F_{f_{\mathrm{ref}}}=T_{\mathrm{obs}}{\int }_0^{\infty }df\frac{{\Gamma }_{12}^2\left(f\right){\mathcal{H}}^2(f,f_{\mathrm{ref}})}{P_1\left(f\right)P_2\left(f\right)}.$

In fact, in deriving the variance ${\sigma }_{\Omega }^2$ as formulated in Equation (67), it turns out that $F_{f_{\mathrm{ref}}}$ is equal to the inverse variance of the estimator, up to a factor of 2,

(71)${\sigma }_{\Omega }^2=\frac{1}{2}{\left[{\int }_0^{\infty }df\frac{{\Gamma }_{12}^2\left(f\right){\mathcal{H}}^2(f,f_{\mathrm{ref}})}{P_1\left(f\right)P_2\left(f\right)}\right]}^{-1}\equiv \frac{F_{f_{\mathrm{ref}}}}{2}.$

where $P_1$ and $P_2$ are the noise power spectra in detectors 1 and 2, respectively, as defined in Equation (49). In principle, these are not known and need to be estimated carefully so as not to bias the search. When operating in the weak signal, it is reasonable to calculate the noise power spectra from the data directly; however, one should be cautious not to use these estimates on the very same data stretch they have been estimated from. In practice, for stochastic searches using ground-based interferometers, the data are segmented into short time segments and the full estimator takes the average of ${\Omega }_{\mathrm{GW}}$ over the entire dataset. To estimate ${\Omega }_{\mathrm{GW}}$ in a single data segment $\tau$, one may use noise power spectra estimated from other segments within a time window over which detector noise is stationary; for example, in the LVK stochastic searches, the average of $P_i$ in the adjacent segments is taken. Otherwise, it is also possible to estimate the noise power spectrum using a parametric fit, as in [152,153,154], or by using a wavelet expansion using methods as shown for example in [155].

Another approach to construct an estimator for the SGWB is to construct a likelihood function for the signal given the data, and maximise it to obtain a maximum likelihood detection statistic. Let us switch to a compact vector notation, where the data vector $\mathit{d}=(d_1,d_2,...,d_N)$ spans the detector space comprising N detectors, and data covariance is the generalised $N\times N$ correlation matrix $\mathit{C}$, where each entry is $C_{ij}=d_i{d}_j$. Starting with a zero-mean Gaussian assumption for both the signal h and noise $\mathit{n}$ terms in the data, we can write

(72)$\mathcal{L}\left(\mathit{d}\right)\propto \prod _{f,\tau }\frac{1}{{\left|\mathit{C}\right|}^{1/2}}{e}^{\frac{1}{2}{\mathit{d}}^{†}{\mathit{C}}^{-1}\mathit{d}},$

where all observed time segments $\tau$ and frequency modes f are considered independent. A method based on this likelihood employs both auto-correlations and cross-correlations of the datastreams to reconstruct SGWB and is discussed in the context of ground-based detectors in [156]. Using auto-correlations is not recommended in scenarios where the detector noise terms are not independently measured, as these will be inextricable from the auto-power spectra. This likelihood is in fact proposed for future LISA data analysis in [147,148], making full use of the independent estimate of the noise from GW-insensitive channels (see details on Sagnac channels here [157]). A more simple, alternative route can be taken by noting that, in the low-signal approximation, the Gaussian assumption can actually be applied to the average of the two-detector residual $C_{ij}-⟨C_{ij}⟩$, leading to a different formulation,

(73)$\mathcal{L}\left(C_{ij}\right)\propto \prod _{f,\tau }\frac{1}{{\sigma }^2}{e}^{\frac{1}{2}(C_{ij}-⟨C_{ij}⟩){\sigma }^{-2}{(C_{ij}-⟨C_{ij}⟩)}^{\star }}.$

This holds in the limit where data are averaged over many segments such that the distribution of the residuals approximates a Gaussian; otherwise, a single residual will not follow this distribution [158]. In [131], Matas and Romano show that likelihood (72) is well approximated by (73) in the low-signal limit. In this case, maximising either likelihood with respect to ${\Omega }_{\mathrm{GW}}\left(f_{\mathrm{ref}}\right)$ yields the same solution as above. The extension to include multiple baselines is straightforward, assuming each $C_{ij}$ may be considered as an independent measurement. The full likelihood $\mathcal{L}\left(C\right)$ becomes the product over the individual baseline likelihoods as

(74)$\mathcal{L}\left(\mathit{C}\right)\propto \prod _b\mathcal{L}\left({\mathit{C}}_b\right),$

where the baseline index b cycles over the detector pairs $ij$ in the set, without double counting baselines.

The approach can be taken one step further to formulate a hybrid frequentist–Bayesian analysis as proposed in [159] by re-expressing the likelihood in terms of $\Omega$ and the model parameters we would like to constrain,

(75)$\mathcal{L}\left({\widehat{\Omega }}_{\mathrm{GW}}\left(f\right)\right|\Theta )\propto \exp \left[-\frac{1}{2}\sum _{b,f}{\left(\frac{{\widehat{\Omega }}_{\mathrm{GW}}^b\left(f\right)-{\Omega }_{\mathrm{M}}\left(f\right|\Theta )}{{\sigma }_{\Omega ,b}^2\left(f\right)}\right)}^2\right],$

where ${\Omega }_M$ is the fractional energy density of the model M which is being tested, $\Theta$ are the parameters on which the model depends, and ${\widehat{\Omega }}_{\mathrm{GW}}^b\left(f\right)$ is an estimator for the strength of the SGWB in a single frequency bin (see, e.g., [160]).

Equation (75) is used as the basis for parameter estimation of models for the GWB, and Bayesian model selection [159]. In the case of parameter estimation, one estimates the posterior distribution of the parameters of the model, $\Theta$,

(76)$\begin{array}{c}\hfill p(\Theta |{\widehat{\Omega }}_{GW}\left(f\right),M)=\frac{\mathcal{L}\left({\widehat{\Omega }}_{\mathrm{GW}}\left(f\right)\right|\Theta )p(\Theta )}{p\left[{\widehat{\Omega }}_{\mathrm{GW}}\left(f\right)\right|M]},\end{array}$

where now we have specified model M explicitly, $p(\Theta )$ is the prior on the model parameters, and

(77)$\begin{array}{c}\hfill p\left({\widehat{\Omega }}_{\mathrm{GW}}\left(f\right)\right)=\int \mathcal{L}\left({\widehat{\Omega }}_{\mathrm{GW}}\left(f\right)\right|\Theta )p(\Theta )d\Theta \end{array}$

is the model evidence or “marginal likelihood”. The posterior is typically estimated either by brute force estimation or, in high dimensional cases, by MCMC methods.

In addition to estimating the posterior in Equation (76), one can perform Bayesian model selection by comparing the marginal likelihood for two separate models. This is used to construct an odds ratio between two models $\mathcal{A}$ and $\mathcal{B}$,

(78)$\begin{array}{c}\hfill {\mathcal{O}}_{\mathcal{B}}^{\mathcal{A}}=\frac{\pi \left(\mathcal{A}\right)}{\pi \left(\mathcal{B}\right)}\frac{p\left[{\widehat{\Omega }}_{\mathrm{GW}}\left(f\right)\right|\mathcal{A}]}{p\left[{\widehat{\Omega }}_{\mathrm{GW}}\left(f\right)\right|\mathcal{B}]}.\end{array}$

The first term on the right-hand side is the prior odds, where one specifies any prior information about the preference for one model over the other. The second term is simply the ratio of evidences. Odds ratios, or Bayes factors (simply the second term on the right-hand side, i.e., equal prior odds), can be used to distinguish between models for the GWB. A general heuristic guide for their use is discussed in Chapter 3 of [33]. Bayes factors have been used recently as the detection statistic in searches that seek to measure alternative polarizations of GWs [56,160,161], in searching for a background from superradiant instabilities around black holes [162,163], as well as searching for models of cosmological backgrounds [164,165,166]. This hybrid-Bayesian methodology has also been proposed in attempting to distinguish between a GWB and globally correlated magnetic noise [167], and simultaneously search for multiple GWB models contributing at once [168]. We discuss the problem of globally correlated magnetic noise for current and future detectors in more detail in Section 5. In addition, a Fisher matrix formalism has also been proposed to simultaneously search for multiple GWB contributions [169].

The detection approaches outlined above have been validated both through mock data challenges [170] and hardware injection campaigns [171]. The latter consisted in generating an isotropic Gaussian stochastic signal by manipulating the test masses at both LIGO sites and successfully recovering the injection via a cross-correlation stochastic search. However, the mock data challenge in [170] considers the case of a semirealistic astrophysical background of binary black holes and neutron stars. Here, the authors prove that the cross-correlation statistical approach assuming a Gaussian background remains valid even in the case of an intermittent astrophysical signal, although it is suboptimal. The fully optimal search methods for this sort of signal are described in Section 4.4.

4.2.2. Pulsar Timing Arrays

In describing the typical Bayesian analysis that is currently performed by most PTAs to search for an isotropic GWB, we follow the methods laid out in, e.g., [172,173,174,175,176], and refer to more specific studies throughout our discussion. Unlike for ground-based interferometers where the detector strain time series is uniformly sampled with equivalent time stamps at each detector, the data collected by PTAs are the arrival times of pulses from an ensemble of pulsars on the sky. Hence, it is preferable to work in the time domain directly, and the typical starting point is the timing residuals of all pulsars, $\delta \mathit{t}$, which are left over after subtracting off the best-fit timing model constructed using deterministic parameters of the pulsars (e.g., rotation frequency, spin-down, binary parameters, sky position, etc.). The likelihood of obtaining n timing residuals given model parameters $\mathbf{\theta }$ for a given pulsar is a multivariate Gaussian as

(79)$\mathcal{L}\left(\delta \mathit{t}\right|\mathbf{\theta })\propto \frac{1}{{\left|\mathit{C}\right|}^{1/2}}{e}^{-\frac{1}{2}{(\delta \mathit{t}-\mathit{s})}^{\mathrm{T}}{\mathit{C}}^{-1}(\delta \mathit{t}-\mathit{s})},$

where $\mathit{s}$ includes contributions from deterministic signals with explicitly modelled time dependence, including the evolution of arrival times as described by the timing model. The covariance matrix $\mathit{C}$ represents contributions from stochastic processes. Note that while we have used “∝”, the specific normalization does matter in this analysis in general, as the spectra of some of the stochastic processes that contribute to $\mathit{C}$ are themselves described by a set of free parameters.

The diagonal elements of $\mathit{C}$ yield “white” noise that is uncorrelated in time whereas off-diagonal elements of $\mathit{C}$ manifest as time-correlated “red” processes, including the GWB. Red processes are modeled in the frequency domain with a power-law for the power spectral density of timing residuals

(80)$P\left(f\right|A,\gamma )=\frac{A^2}{12{\pi }^2}{\mathrm{yr}}^3{\left(f\mathrm{yr}\right)}^{-\gamma },$

where f is the GW or noise frequency, and A and $-\gamma$ are the power-law amplitude and spectral index, respectively. We discuss sources of red noise that affect pulsars individually in more detail in Section 5.2.2. For the GW case, amplitude A corresponds to the GW strain amplitude at the reference frequency $f_{\mathrm{ref}}={\mathrm{yr}}^{-1}$ and $\gamma$ corresponds to $\alpha$ in Equation (25) and ${\alpha }^{\prime }$ explained in footnote ³: $\gamma =3-2{\alpha }^{\prime }=5-\alpha$ = 13/3. The lowest frequency and the size of the frequency bin $\Delta f$ are usually selected to be equal to the inverse of the total observation time $T_{\mathrm{obs}}$, which is typically more than a decade. The data are collected every few weeks. The observed power spectral density of timing residuals $P\left(f\right)$ [s${}^3$] is derived from the GW strain power spectral density $S\left(f\right)$ (s) in Section A4 in [177]. Neglecting frequencies f other than multiples of $T_{\mathrm{obs}}^{-1}$, the densities and strain are related as

(81)$P\left(f\right)=\frac{S\left(f\right)}{12{\pi }^2{f}^2}=\frac{h{\left(f\right)}^2}{12{\pi }^2{f}^3}.$

With each pulsar modelled by tens of timing model parameters and the necessity to invert the covariance matrix that models both temporal and spatial correlations for every likelihood evaluation, in practice it is necessary to marginalise the likelihood over nuisance parameters, which include timing model parameters and Gaussian coefficients that yield a specific realization of red noise. Below, we describe the marginalization procedure [178,179] that is employed by recent searches (e.g., [26,54,180]). First, the likelihood is rewritten as

(82)$\begin{array}{cc}\hfill \mathcal{L}\left(\delta \mathit{t}\right|\mathit{ϵ},\mathit{A},\mathit{\eta })& \propto \frac{1}{{\left|\mathit{N}\right|}^{1/2}}{e}^{-\frac{1}{2}{(\delta \mathit{t}-\mathit{s})}^T{\mathit{N}}^{-1}(\delta \mathit{t}-\mathit{s})},\hfill \end{array}$

(83)$\begin{array}{cc}\hfill \mathrm{where}\mathit{s}& =\mathit{M}\mathit{ϵ}+\mathit{F}\mathit{a}+{\mathit{s}}^{\prime }.\hfill \end{array}$

The covariance matrix $\mathit{N}$ is now diagonal and, thus, represents the white noise associated with each observation; $\mathit{ϵ}$ represents deterministic timing model parameters and the “design matrix” $\mathit{M}$ maps them to the timing residuals; $\mathit{F}\mathit{a}$ represents low frequency “red” noise. The latter is modelled by a Fourier series with a Fourier amplitudes vector $\mathit{a}$ and a matrix $\mathit{F}$ of unit-norm frequency Fourier modes. Typically, the frequencies of the Fourier series are multiples of $1/T_{\mathrm{obs}}$ [174]. $\mathit{\eta }$ represents hyperparameters that govern the shape of the red noise spectrum defined by $\mathit{F}$ and $\mathit{a}$, and ${\mathit{s}}^{\prime }$ denotes deterministic terms that are not marginalised over. To be precise, the matrix $\mathit{N}$ is block-diagonal, consisting of the nominal error in pulse arrival time, but parameterised by “EFAC,” which multiplies arrival time errors; “EQUAD”, which is added to the arrival time errors in quadrature; and “ECORR”, which accounts for correlated measurement uncertainties in contemporaneous, multi-band observations and contributions to the off-diagonal terms [54,179,180].

Further compactifying the equations, we express

(84)$\begin{array}{cc}\hfill \mathit{T}& =\left(\begin{array}{cc}\mathit{M}& \mathit{F}\end{array}\right),\hfill \end{array}$

(85)$\begin{array}{cc}\hfill \mathit{b}& =\left(\begin{array}{c}ϵ\\ \mathit{A}\end{array}\right).\hfill \end{array}$

Gaussian priors are then placed on $\mathit{b}$,

(86)$\begin{array}{c}\hfill \pi \left(\mathit{b}\right|\mathit{\eta })=\frac{exp\left(-\frac{1}{2}{\mathit{b}}^{\mathrm{T}}{\mathit{B}}^{-1}\mathit{b}\right)}{\sqrt{2\pi det\left\{\mathit{B}\right\}}},\end{array}$

with covariance $\mathit{B}$, for which its parameters have previously been referred to as $\mathit{\eta }$,

(87)$\begin{array}{c}\hfill \mathit{B}=\left(\begin{array}{cc}\infty & 0\\ 0& \mathit{\phi }\end{array}\right),\end{array}$

where we have left the prior on the timing model parameters unconstrained (the upper left corner indicates a block with diagonal entries of infinity). The posterior is then proportional to the original multivariate likelihood described by Equation (79) times the prior described by Equation (86),

(88)$\begin{array}{c}\hfill \mathcal{P}(\mathit{b},\mathit{\eta }|\delta \mathit{t})\propto \mathcal{L}(\delta \mathit{t}\left|\mathit{b}\right)\times \pi \left(\mathit{b}\right|\mathit{\eta })\times \pi (\mathit{\eta }).\end{array}$

The initial timing model parameters and the design matrix are obtained with least-squares fitting. The matrix $\mathit{\phi }$ contains information about the spectrum of low frequency noise, including intrinsic red noise for each pulsar and the GWB. If we label pulsars as a and b and frequency bins (multiples of $1/T_{\mathrm{obs}}$) as i and j, then the elements of this matrix are

(89)$\begin{array}{c}\hfill {\left[\mathit{\phi }\right]}_{\left(ai\right),\left(bj\right)}={\Gamma }_{ab}{\rho }_i{\delta }_{ij}+{\eta }_{ai}{\delta }_{ab}{\delta }_{ij},\end{array}$

where ${\rho }_i$ is the spectrum of the GWB (and is, therefore, present in all cross-correlations and auto-correlations), and ${\eta }_{ai}$ is the spectrum of intrinsic red noise of pulsar a in frequency bin i (and is, therefore, only present on the diagonals). The vectors of maximum-likelihood values of $\mathit{b}$, $\widehat{\mathit{b}}$, are obtained by solving

(90)$\begin{array}{cc}\hfill \widehat{\mathit{b}}& ={\Sigma }^{-1}\mathit{d},\hfill \end{array}$

(91)$\begin{array}{cc}\hfill \Sigma & ={\mathit{T}}^{\mathrm{T}}{\mathit{N}}^{-1}\mathit{T}+{\mathit{\phi }}^{-1},\hfill \end{array}$

(92)$\begin{array}{cc}\hfill \mathit{d}& ={\mathit{T}}^{\mathrm{T}}{\mathit{N}}^{-1}\delta \mathit{t}.\hfill \end{array}$

Alternatively, one can explicitly marginalise over $\mathit{b}$ as

(93)$\begin{array}{cc}\hfill \mathcal{P}\left(\mathit{\eta }\right|\delta \mathit{t})& =\int d\mathit{b}\mathcal{P}(\mathit{b},\mathit{\eta }|\delta \mathit{t})\hfill \\ \hfill & =\frac{exp(-\frac{1}{2}\delta {\mathit{t}}^T{\mathit{C}}^{-1}\delta \mathit{t})}{\sqrt{2\pi detC}},\hfill \end{array}$

where $\mathit{C}=\mathit{N}+\mathit{T}\mathit{B}{\mathit{T}}^T$. $\mathit{C}$ can be efficiently inverted using the Woodbury formula to increase computational efficiency.

In principle, $\mathit{\phi }$ is block diagonal with independent frequency bins but correlated noise between pulsars when a GWB is present. However, a GWB signal is expected to reveal itself as a strong common red noise process before cross-correlations are evident. Recent analyses have, therefore, first evaluated the term in Equation (89) with ${\Gamma }_{ab}={\delta }_{ab}$ [178] as the null hypothesis and ${\Gamma }_{ab}$ given by Equation (63) as the signal hypothesis. As we discuss in Section 5.2, there is strong evidence for a common spectrum process (${\Gamma }_{ab}={\delta }_{ab}$) compared to the intrinsic red-noise only hypothesis (${\Gamma }_{ab}=0$), but little compelling evidence for cross-correlations compared to a common spectrum process.

The ${\eta }_{ai}$ and ${\rho }_i$ spectra are modeled as a power law,

(94)$\begin{array}{cc}\hfill {\eta }_{a,i}& =P_{\mathrm{red}}\left(f_i\right|A_{a,\mathrm{red}},{\gamma }_{a,\mathrm{red}})\Delta f,\hfill \end{array}$

(95)$\begin{array}{cc}\hfill {\rho }_i& =P_{\mathrm{GW}}\left(f_i\right|A_{\mathrm{GW}},{\gamma }_{\mathrm{GW}})\Delta f,\hfill \end{array}$

and so instead of having $N_{\mathrm{freq}}$ parameters per pulsar (plus $N_{\mathrm{freq}}$ more parameters for the GWB), one has $2(N_{\mathrm{pulsar}}+1)$ free hyper-parameters that specify the shape of noise processes. White noise parameters, such as the typical EFAC, EQUAD, and ECORR parameters are measured on a per-pulsar (and per-observatory or per-observing back-end) basis. Typical pulsar datasets are modelled by tens of such parameters; thus, white noise parameters are often fixed to their maximum likelihood values that were calculated by analyzing pulsars individually.

Bayesian inference with the marginalised likelihood typically yields posterior samples for the $2(N_{\mathrm{pulsar}}+1)$ hyperparameters, and the Bayesian evidence, which is the integral of the likelihood over the prior. The evidence is used to select a model that best describes the data. In some recent analyses, explicit Bayes factors (the ratio of evidences between models) are calculated using a “product-space” approach that foregoes individual evidence calculations in favor of sampling from multiple models simultaneously [54,180].

It is also common to consider an optimal statistic which is built similar to the one presented for LVK searches [149,181,182]. We follow the conventions and methods in [182], which were implemented in [54]. We start from the auto-covariance and cross-covariance matrices calculated from the timing residuals,

(96)$\begin{array}{cc}\hfill {\mathit{C}}_a& =⟨\delta {\mathit{t}}_a\otimes \delta {\mathit{t}}_a⟩,\hfill \end{array}$

(97)$\begin{array}{cc}\hfill {\mathit{S}}_{ab}& =⟨\delta {\mathit{t}}_a\otimes \delta {\mathit{t}}_b⟩{|}_{a\ne b}.\hfill \end{array}$

For a GWB with PSD given by $P_{\mathrm{GW}}\left(f\right)$, the cross-covariance matrix is given by

(98)$\begin{array}{cc}\hfill {\mathit{S}}_{ab}& =P_{\mathrm{GW}}\left(f\right)\left(F_a{\Gamma }_{ab}{F}_b^T\right).\hfill \end{array}$

Specifically, for a background dominated by supermassive black hole binaries that evolve solely due to GW emission, $\gamma =13/3$. The optimal statistic, ${\widehat{A}}^2$, is then given by [149,181,182]

(99)${\widehat{A}}^2=\frac{{\displaystyle \sum _{ab}}\delta {\mathit{t}}_a^T{C}_a^{-1}{\tilde{\mathit{S}}}_{ab}{C}_b^{-1}\delta {\mathit{t}}_b}{{\displaystyle \sum _{ab}}{\mathit{C}}_a^{-1}{\tilde{\mathit{S}}}_{ab}{\mathit{C}}_b^{-1}{\tilde{\mathit{S}}}_{ba}},$

with

(100)$A_{\mathrm{GW}}^2{\tilde{\mathit{S}}}_{ab}={\mathit{S}}_{ab}.$

In the second line, we define the amplitude-independent cross-correlation matrix, which makes $⟨\widehat{A^2}⟩=A_{\mathrm{GW}}^2$. If $A_{\mathrm{GW}}=0$, then the variance of the optimal statistic is given by

(101)$\begin{array}{c}\hfill {\sigma }_0={\left[\sum _{ab}{\mathit{C}}_a^{-1}{\tilde{\mathit{S}}}_{ab}{\mathit{C}}_b^{-1}{\tilde{\mathit{S}}}_{ba}\right]}^{-1/2},\end{array}$

which allows us to construct an SNR $\rho ={\widehat{A}}^2/{\sigma }_0$. The SNR significance can be found by estimating its noise-only distribution using simulations of the data with pulsar sky positions randomly assigned or random phase shifts applied to the data to remove a signal [179,183].

The implementation of the optimal statistic requires evaluating ${\mathit{C}}_a$ and ${\mathit{S}}_{ab}$ for a given choice of hyperparameters that define $\mathit{a}$, which are amplitudes of sinusoids on which the red noise is decomposed. One typically uses hyperparameters that correspond to $\widehat{\mathit{b}}$ in Equation (92). One can also use the maximum a posteriori estimates of the red noise parameters from the Bayesian analysis, $A_{\mathrm{red},a}^{\max }$, ${\gamma }_{\mathrm{red},a}^{\max }$, $A_{\mathrm{GW}}^{\max }$, and ${\gamma }_{\mathrm{GW}}^{\max }$, but this can result in a large bias in the estimate of ${\widehat{A}}^2$ because the individual red noise parameters and the common red noise parameters are highly correlated. In practice, the red noise parameters are drawn from the posterior samples from the Bayesian analysis and for each draw the optimal statistic and its SNR are calculated. One can estimate a value of ${\widehat{A}}^2$ from this distribution, which is significantly less biased than if one calculates the optimal statistic once from the maximum a posteriori estimates of the noise parameters [182].

While our focus in this section has been on methods for detecting an isotropic GWB with PTAs, note that there is also a recent paper that aims to measure both continuous sources from individual supermassive black hole binaries as well as an isotropic confusion background [184].

4.3. Anisotropic Background Detection Methods

Mapping an SGWB is a distinct problem from that of tracing back an individual source position. Individual source mapping leverages the coherence of radiation arriving from the source, using the phase information of the signal to estimate both angular and frequency-dependent phase along with other intrinsic and extrinsic source parameters. In either case, however, a basic requirement for directional analyses is multiple observations of the same signal. These multiple observations can be made in different ways. Multiple detectors at different locations will impose a known phase shift on a coherent signal, which can be used to constrain the position on the sky. In the case of a persistent, coherent source, unequal-time auto-correlation of a single detector signal can yield angular phase information, and in the case of some sources of GWs, matched filtering can be used to pinpoint the location on the sky with a significantly larger effective baseline [185]. For the case of a persistent and incoherent source, such as an SGWB, one may use equal-time correlations of data from different detectors and the motion of individual detectors with respect to the frame in which the signal is stationary to reconstruct angular information. We focus on this approach here.

Let us start by discussing the weak signal limit, which applies to the 2G ground-based detector network (e.g., LIGO and Virgo). As we did for the isotropic case, we can assume that the residuals for the data cross-correlation spectra and our model are Gaussian, similarly to Equation (73). The log-likelihood may be written as a function of an anisotropic signal I explicitly as

(102)$ln\mathcal{L}\left(C^{\tau }\left(f\right)\right|I(f,\widehat{\mathit{n}}))\propto \frac{1}{2}\left[(C^{\tau }\left(f\right)-⟨C^{\tau }\left(f\right)⟩){\mathbf{\sigma }}^{-2}{(C^{\tau }\left(f\right)-⟨C^{\tau }\left(f\right)⟩)}^{\star }+\mathrm{Tr}ln{\mathbf{\sigma }}^{\mathbf{2}}\right],$

where we generically model the correlated data as

(103)$⟨C^{\tau }\left(f\right)⟩={\int }_{S^2}d\widehat{\mathit{n}}{A}^{\tau }(f,\widehat{\mathit{n}})I(f,\widehat{\mathit{n}}).$

The operator A projects the sky signal I into the datastream, and contains all the necessary correlated response information of the detectors. For ground- and space-based detectors, A is given by Equation (59) up to appropriate normalisation factors. We have reintroduced the timestamp label $\tau$ here as it serves an important purpose: tracking the time-dependence of the detector response on the sky; the sky signal I is considered absolutely stationary¹⁰ and, thus, $\tau$-independent. For ground-based interferometric detectors, operator A depends only upon the sidereal time, which means the data stream can be folded on the sidereal day before we invert Equation (103) to solve for I [186,187,188]. To probe the effects of the projection operator A on the GW observations, it is useful to simulate a signal as a single point source on the sky and “project” it analytically (or numerically) into cross-correlated mock data, which produces the point-spread function of the detector pair. This has often been studied to glean the effective angular resolution of the detectors to a stochastic source—we refer the reader to [189] for a complete presentation of the subject.

We can now derive a maximum-likelihood mapping solution by maximising the likelihood with respect to the sky map. To best illustrate the structure of the solution, let us use a compact matrix notation where operators in bold generally span multiple dimensions, i.e., time, frequency, direction, and baseline, such as $A_b^{\tau }(f,\widehat{\mathit{n}})\to \mathit{A}$. These are contracted appropriately to achieve the desired result,

(104)$\tilde{\mathit{I}}={\left({\mathit{A}}^{†}{\Sigma }^{-1}\mathit{A}\right)}^{-1}{\mathit{A}}^{†}{\Sigma }^{-1}\mathit{C}.$

In the simple case of a single frequency and single time segment of data, we can think of $\mathit{C}$ as a vector containing cross-correlations between different baselines, $\Sigma$ is the noise covariance of the cross-correlation of the detector pairs, and $\mathit{A}$ is now a matrix that contains geometric factors that project different sky directions onto correlations between different detectors. Explicitly, $\Sigma$ entries are of the form

(105)${\Sigma }_b^{\tau }(f,f^{\prime })=\mathrm{Cov}{\left[P_i^f{P}_j^{f^{\prime }}\right]}_{\tau },$

where $i,j$ label the two detectors in the baseline b. More details on the calculation may be found in [35,189,190].

In practice, this solution may need to be implemented iteratively as the noise covariance $\Sigma$ may not be known a priori and, hence, may need to be estimated from the data themselves. In the low signal limit however, as is the case with current ground-based networks, the noise covariance is estimated independently of the signal and Equation (104) is, thus, a closed-form mapping solution. The first term on the right-hand side is the Fisher information matrix,

(106)$\mathit{F}={\mathit{A}}^{†}{\Sigma }^{-1}\mathit{A},$

which is inverted and applied to the projection map,

(107)$\mathit{z}={\mathit{A}}^{†}{\Sigma }^{-1}\mathit{C},$

to extract the maximum-likelihood solution $\tilde{\mathit{I}}$. z is often referred to as the dirty map in the literature, as it effectively represents a first, naïve representation of the cross-spectral density on the sky, while $\tilde{\mathit{I}}$ is considered the clean map as it is the result of a deconvolution. The Fisher matrix contains the pixel–pixel covariance information necessary to deconvolve the detector response from the dirty map. In principle, this method can also be employed for anisotropic searches in a strong signal regime; in this case, if the signal is Gaussian, the map point-estimate remains accurate; however, covariance $\Sigma$ requires extra terms as the signal contribution may not be ignored.

However, in the presence of a strong signal term, it is preferable to include auto-correlation terms in the likelihood, as these contain precious information. This is the case for, e.g., LISA and PTAs. Hence, we consider again the Gaussian likelihood (72), and as in [191,192,193,194], we extend it to include directional information,

(108)$\ln \mathcal{L}\left({\mathit{d}}^{\tau }\left(f\right)\right|I(f,\widehat{\mathit{n}}))\propto \frac{1}{2}\left[{\mathit{d}}^{\tau }\left(f\right){\mathit{C}}^{-1}{\mathit{d}}^{\tau \star }\left(f\right)+\mathrm{Tr}ln\mathit{C}\right],$

where the correlation matrix ${\mathit{C}}^{\tau }\left(f\right)$ (re-introducing the time and frequency dependence explicitly here for clarity) may be estimated as

(109)$\mathit{C}\left(f\right)={\int }_{S^2}d\widehat{\mathit{n}}\mathit{A}(f,\widehat{\mathit{n}})I(f,\widehat{\mathit{n}})+\mathit{N},$

with $\mathit{N}$ as the noise covariance. Note that $\mathit{A}$ is an $N\times N$ matrix representing the correlated response of the full network.

The above notation implies that we are discussing LISA because we have specified the frequency domain and that $\mathit{A}$ depends on $\tau$. In principle, PTA analysis works analogously—note the similarity of Equation (109) to Equation (89), with the added complication of a dependence on direction. However, it is worth briefly reiterating specific differences. First, the “detector response” is time-independent for PTAs, and in general we assume it is frequency independent. In addition, PTAs leave data in the time domain but expand the GW signal and intrinsic red noise in the frequency domain using a set of Fourier basis functions with an appropriate prior to impose a power-law spectrum and correlations between pulsars (see Equations (87)–(92)). It is the correlations between pulsars that are directly affected, with the anisotropy of the background causing a deviation from the Hellings–Downs curve that depends upon the specific realization of the background. This means that Equation (89) changes:

(110)$\begin{array}{c}\hfill {\Gamma }_{ab}{\rho }_i\to {\int }_{S^2}d\widehat{\mathit{n}}{\Gamma }_{ab}\left(\widehat{\mathit{n}}\right){\rho }_i\left(\widehat{\mathit{n}}\right),\end{array}$

with i once again specifying the frequency index, a and b specifying pulsars we are correlating, and ${\rho }_i\left(\widehat{n}\right)$ specifying the GWB amplitude at frequency index i in direction $\widehat{\mathit{n}}$. In this case, ${\Gamma }_{ab}\left(\widehat{\mathit{n}}\right)$ can be considered a direct analogue to $\mathit{A}(f,\widehat{\mathit{n}})$.

In [193], the authors show that extracting the maximum-likelihood map solution from Equation (108) yields

(111)$\begin{array}{cc}\hfill \tilde{\mathit{I}}& ={\mathit{F}}^{-1}·\frac{1}{2}\mathrm{Tr}\left[{\mathit{C}}^{-1}\frac{\partial \mathit{C}}{\partial \tilde{\mathit{I}}}{\mathit{C}}^{-1}(\mathit{d}\otimes {\mathit{d}}^{\star }-\mathit{N})\right],\hfill \end{array}$

(112)$\begin{array}{cc}\hfill \mathit{F}& =\frac{1}{2}\mathrm{Tr}\left[{\mathit{C}}^{-1}\frac{\partial \mathit{C}}{\partial \tilde{\mathit{I}}}{\mathit{C}}^{-1}\frac{\partial \mathit{C}}{\partial \tilde{\mathit{I}}}\right],\hfill \end{array}$

where $\mathit{F}$ is the Fisher information matrix. This solution is in fact structurally identical to the Bond–Jaffe–Knox mapping procedure in CMB [195,196] and requires an iterative implementation: first, a “guess” map is injected in $\mathit{C}$ (Equation (109)), then a first solution $\tilde{\mathit{I}}$ is calculated using Equation (111), which updates the estimate of $\mathit{C}$ and so on until convergence. Beyond SNR considerations, this iterative method is also useful in the case of correlated noise terms as these may be included in the noise model $\mathit{N}$, whereas solution (104) rests upon the assumption that the noise is completely uncorrelated between different detectors/measurements.

Alternatively, in the strong or intermediate signal regime one can employ a fully Bayesian approach to the problem as presented in, e.g., [192,194,197]. In [192], the authors develop a fully Bayesian search pipeline to use PTAs to map the GW sky using the formalism detailed in [191]. In [194], the authors present a Bayesian search for anisotropies in the GWB using LISA, where they decompose the sky onto the first three spherical harmonics, parameterize the noise covariance curve, and use MCMC sampling to map the diffuse sky and estimate the noise properties of the detector simultaneously.

Several different approaches can be taken when applying these map-making methods to real data [157,186,189,190,191,192,197,198,199,200,201,202,203,204] based on specific necessities. The spectral dependence of the map is often factored out (although not always, see, e.g., [154,169,202,205]), under the assumption in Equation (26) and integrated over to obtain broadband maps of the GW sky. Then, the direction-dependent component of the signal may be expanded and mapped in spherical harmonics,

(113)$\begin{array}{c}\hfill I\left(\widehat{\mathit{n}}\right)=\sum _{\ell =0}^{\infty }\sum _{m=-\ell }^{\ell }{I}_{\ell m}{Y}_{\ell m}\left(\widehat{\mathit{n}}\right),\end{array}$

or discretized directly on the sky in pixel space,

(114)$\begin{array}{c}\hfill I\left(\widehat{\mathit{n}}\right)=\sum _{p^{\prime }=1}^{N_{\mathrm{pix}}}{I}_{p^{\prime }}{\delta }_{p{p}^{\prime }},\end{array}$

where $Y_{\ell m}$ are the standard spherical harmonics, and we have used p and $p^{\prime }$ to enumerate pixels on the sky that are specified using some choice of pixelization (e.g., HealPix [206])¹¹. Other works have argued for the use of curl and gradient spherical harmonics, as in CMB analyses [199,200,201], or bipolar spherical harmonics [207], although to our knowledge these have not yet been employed on real data.

In both the pixel and spherical harmonic bases, the Fisher matrix $\mathit{F}$ is usually very poorly conditioned for a single baseline, as the sky coverage is incomplete; this is the case even after time-integration, which takes into account the sky-scanning strategy of the baseline. Regularisation techniques are, therefore, needed to estimate maximum likelihood solutions such as those presented in Equations (112) and (106). Examples of regularization techniques include cutting off the sum in (113) a priori at some maximum ℓ [190], singular value decomposition (SVD) (see, e.g., [208]), and the drastic approach of limiting the calculation to the diagonal of $\mathit{F}$, ignoring all off-diagonal correlations [198]. In some cases, combinations of these approaches are used.

The spherical harmonic basis has mostly been employed when searching for a diffuse background, which dominates the low-ℓ modes, while pixel space has been employed when searching for point-like stochastic sources, although it has been shown that there is a one-to-one mapping between the two [209] (see, e.g., Figures 1 and 2 of [192] for an illustrative example of how specific anisotropic maps correspond to different spherical harmonic decompositions). There may however be subtle differences between the conditioning in spherical harmonic and pixel space; in interferometry in general, the choice has often been to obtain a spherical harmonic domain solution since the measurements directly constrain this domain (see e.g., [210] for an application to CMB interferometry). The reconstruction of the actual “image” is then left as a separate and straightforward conversion problem.

However, as discussed in [152], in the case of GW detectors, there is a non-trivial coupling of spherical modes in the sky response due to the non-compactness of the beam, indicating that adopting this approach may have significant drawbacks. In particular, there is no isotropic kernel for a primary beam describing the coupling of individual pairs of directions, hence limiting the number of spherical modes in the solution in a piece-wise fashion could make the problem more ill-conditioned than it actually is. Finally, in the case of PTAs it has been shown that, when working in the spherical harmonics basis, the non-uniformity of pulsars on the sky makes it impossible to decouple the monopole from other multipole moments, making the search less sensitive to isotropic backgrounds [203,204]. In those studies, the authors propose using a basis formed with principal component analysis of the direction-dependent sensitivity of the specific PTA.

In fully Bayesian examples, using too many spherical harmonics terms or pixels on the sky can make the problem intractable or, in a detection case, could widen the parameter space beyond what is possible to resolve. Techniques that can be used to make the full Bayesian problem more tractable might be to only search over the first few spherical harmonics, as in, e.g., [192,194]. Other techniques have also been proposed, one of which constructs a set of basis “sky vectors” from only the columns of the SVD matrices that are associated with non-zero singular values of the Fisher matrix. One then samples over their coefficients to reconstruct the map using a lower-dimensional representation of the sky [197,201].

Although notably more challenging when attempting map-making with a singular $\mathit{F}$, the spectral dependence may be solved for as well, resulting in narrowband directional searches. These have been carried out by the LVK collaboration targeting point-like sources which may have an interesting spectral emission (e.g., asymmetric neutron stars, supernova remnants, the galactic bulge) and have been performed in pixel space [187,188,211]. As point sources are not correlated with other points on the sky, this type of search ignores the off-diagonal correlations in the Fisher matrix to avoid the inversion problem. Another approach employed on LIGO data in [154] is that of frequency-banding, where one solves for maps in separate frequency bands spanning several Hz, making minimal assumptions about the spectrum of the signal in each band. This allows both inverting the full Fisher matrix and obtaining full sky maps in each band and attempt a model-free reconstruction of the frequency spectrum. In the case of pulsar timing arrays, the signal is dominated by a small ($\lesssim 10$) number of frequencies, and so in the fully Bayesian analysis in [202], they estimate different spherical harmonic coefficients at different frequencies to allow for maps that are dominated by point sources in different directions at different frequencies.

4.4. The Approach towards Non-Gaussian Backgrounds

The cross-correlation searches reviewed above have been shown to be optimal in the case of a continuous background, which, for example, presents a timestream as in the third row of Figure 3; however, these will be suboptimal in the case of a background which is intermittent, such as the first row of Figure 3. As argued in [212], an intermittent background will not satisfy the central limit theorem; thus, searching for such a signal under the Gaussian assumption will reduce the sensitivity of the search, resulting in delays until detection (although it will not bias the point estimate). Given the observations of GWs accumulated up until now [213] and the lack of a detection of a Gaussian stochastic background, it is clear that a signal of sub-threshold GW bursts from compact binary coalescences is present in LIGO and Virgo data and may be the first stochastic signal that can be detected with ground-based detectors.

A few optimised search methods for non-Gaussian backgrounds have been proposed: A maximum likelihood approach [212] where each isolated GW making up the background is modelled as a “stochastic burst” of GW strain, correlated between detectors; a Bayesian search method which leverages deterministic parameter estimation [214,215], fitting each burst with a compact binary waveform; and, soon, a fully Bayesian stochastic search [216] which extends the formalism of [212]. This method has also been extended to estimate the anisotropic properties of the subthreshold CBC population [217].

While substantially different in formulation and implementation, these approaches share the premise of re-parametrisation of the background signal based on a Gaussianity parameter, $\xi$, which quantifies the degree of non-Gaussianity. Effectively, this is the same as the duty cycle introduced in Section 3.4 and is linked to the duration of each individual GW burst and the frequency with which they appear in the detector. If normalised, $0\le \xi \le 1$ is the probability that a burst is present in the detector at any given time, and if $\xi =1$, the background is Gaussian. The detector noise, on the other hand, is usually still assumed to be Gaussian¹². Then, that assuming we segment data into short segments labelled by i, the likelihood function for the data vector ${\mathit{d}}_i$ (for a set of detectors) given the data covariance $\mathit{C}$ and signal model h is

(115)$\mathcal{L}\left({\mathit{d}}_i\right)\propto \prod _{f,\tau }\frac{1}{{\left|\mathit{C}\right|}^{1/2}}{e}^{\frac{1}{2}({\mathit{d}}_i-{\mathit{h}}_i){\mathit{C}}^{-1}{({\mathit{d}}_i-{\mathit{h}}_i)}^{\star }}.$

For an intermittent background, model ${\mathit{h}}_i$ in the likelihood will be either a GW burst or a CBC waveform, if the data segment ${\mathit{d}}_i$ has a signal, or simply equal to zero if it does not [215]. The full likelihood ${\mathcal{L}}_{\mathrm{full}}$ of the data with an intermittent signal may then be written as the mixture of a signal-and-noise likelihood ${\mathcal{L}}_s$ and a noise-only likelihood ${\mathcal{L}}_n$ as

(116)${\mathcal{L}}_{\mathrm{full}}\left({\mathit{d}}_i\right)=\xi {\mathcal{L}}_s\left({\mathit{d}}_i\right|{\mathit{h}}_i)+(1-\xi ){\mathcal{L}}_n\left({\mathit{d}}_i\right|0),$

featuring the duty cycle $\xi$ as a parameter. Maximising this likelihood (with any of these approaches) results in posterior distributions for both the signal parameters that construct $h_i$, and $\xi$; furthermore, it is also possible to include noise parameters here such that these also are estimated.

5. Current Detection Efforts of SGWBs

A variety of searches based on the methods described in Section 4 have been implemented in interferometric and timing data, in the case of ground-based detectors and pulsar timing arrays (PTA), respectively, to either detect or constrain a wide range of stochastic backgrounds. Beyond these, many investigations have been carried out to assess the detection potential of future detectors on the ground and in space and their ability to distinguish different background sources.

In this Section, we start by discussing the most recent analysis results from the ground-based detector LIGO–Virgo–KAGRA (LVK) collaboration, obtained with the third observing run (O3) of the 2G detectors: the Advanced LIGO and the Advanced Virgo. These results build on the previous data runs O1 and O2 of the 2G detectors. The overall sensitivity of the 2G detector network may be observed in Figure 7. We then discuss future detection prospects for the 2G network, which should reach design sensitivity with the next few upgrades of the instruments [218]. When representing the sensitivity to a GWB, we use “power-law integrated” sensitivity curves, which show the sensitivity to a GWB at each frequency (at a chosen SNR level) with a power law spectral shape that is tangent to the curve at that frequency [219].

For PTAs, we provide details about existing and proposed collaborations across the world and report on their searches for nanohertz GWBs. We review the upper limits on the amplitude of the isotropic GWB from circular supermassive black hole binaries. Next, we briefly outline main sources of noise that limit the searches and discuss the future of timing array experiments. We close by briefly discussing other interesting GW searches, such as resonant bar experiments and the Moon-based proposals, and highlight their bearing on stochastic searches.

5.1. Searches with Ground-Based Laser Interferometers

The ground-based detector communities have put in a continued, concerted effort towards detecting a stochastic GW background between ∼10–${10}^3$ Hz. The principal candidate for a stochastic background in the ground-based detector band is one given by inspiralling and merging binary black holes and neutron stars. This is clear given the sensitivity of the LIGO-Virgo network, as shown for example in Figure 7, and the fact that these detectors have already detected a large set of CBCs [74]. This background has been exclusively searched for via a cross-correlation search, tailored to the astrophysical background in particular by reweighting the data in frequency according to the expected power-law spectral dependence as derived in Equation (36). First efforts of isotropic stochastic searches date back to 2003 and the first science run [228], followed by several other science runs of the LIGO detectors [16,229,230,231], and the Advanced detector era since 2016 [161,232], and, finally, the inclusion of Virgo [56]. Beyond these, several other searches have been carried out by targeting specific GWBs; notably, searches which allow for anisotropy in the signal have regularly been carried out [211,233,234,235,236], as well as searches for cosmic string networks [129,237,238]. Other targeted searches include searches for non-GR polarisation modes [56,161,239]. Furthermore, several stochastic search efforts have been carried out by small teams outside the LVK collaboration; let us cite here a set of directional searches complementary to the LVK ones [153,154,240], a search for correlations between the anisotropic GWB and galaxy catalogues [241], searches for ultralight vector bosons [162,163], a search for a primordial inflationary background [242], and a search for parity-violating stochastic signals [166].

We focus on current search results in this section, detailing the detector characterisation issues that the Advanced detectors have faced up to now. We discuss future challenges for ground-based detectors in Section 6.1, where we explore SGWB detection strategies with third generation interferometers.

5.1.1. Search Results for an Isotropic Background by LVK

Applying the cross-correlation recipe described in Section 4 to the real LIGO-Virgo datasets requires firstly identifying the valid cross-correlation times at which different detectors are simultaneously online and fully operational and, subsequently, Fourier transforming the measured timestreams to the frequency domain. In practice, it is convenient to divide these into smaller time segments and fast-Fourier transform (FFT) each segment, which is treated as an independent measurement. This reduces the cost of Fourier transforming, excludes frequencies that are far into the $1/f$ tail, and ensures that noise properties are stationary within each segment.

In the most recent isotropic stochastic analysis [56], raw data are segmented into 192 s segments which are then Hann-windowed and FFTed. The data are downsampled from the native sampling rate of 16 kHz to 4096 Hz, and the full frequency band included in the analyses spans 20–1726 Hz; this choice is made so as to exclude the highest and lowest frequencies available where noise becomes prohibitively large. Data quality cuts are applied based on a measure of noise stationarity between segments or when there are known detector artifacts in the data, and these cuts reduce viable data by ∼20% in O3. In addition, frequencies where narrow spectral artifacts (“lines”) occur and are caused by known instrumental disturbances or show coherence between detectors are removed [243] and typically cut ∼15–20% of the frequency band. More details on the methodology, exact values of data excised for each baseline in O3, and implementation of the quality checks can be found in [56,243]. Note that this is the first advanced-detector era stochastic analysis to include multiple baselines¹³; we highlight in this section the important implications of adding two long baselines to the network. Specifically, O3 analysis spans Hanford–Livingston (H–L), Hanford–Virgo (H–V), and Livingston–Virgo (L–V).

Several new analysis features were implemented in [56], as the LVK collaboration gears up for a possible SGWB detection. The results are presented by using a log-uniform prior, assigning equal weight to each order of magnitude in ${\Omega }_{\mathrm{GW}}$, in line with our current state of knowledge. Furthermore, a correlated noise analysis targeting magnetic noise resonances was carried out, using the framework set up by the likelihood (75). The analysis confirmed the lack of correlated magnetic noise across the interferometer network and laid out a clear methodology, which will prove essential for similar analyses with future, more sensitive detectors (e.g., 3G). For the first time, short gates (∼1 s) were employed to remove frequent, loud glitches which populated the LIGO detector O3 data [244,245]. These were not as frequent nor problematic in previous data runs, and investigations are still under way to determine the sources of these noise artifacts. Without their removal, loud glitches would bias the PSD of a 192-second segment, which would then not pass the stationarity cut, reducing the effective live time of the LIGO detectors to less than 50%, whereas gated datasets are reduced by about 1% of the original total length.

The H–L baseline contributes the most to overall network sensitivity, while Virgo baselines remain important for higher frequencies and, hence, especially for searches employing higher values of the spectral index $\alpha$. The combined H–L–V spectrum for all three data runs shown in Figure 8 is consistent with Gaussian noise, with symmetric fluctuations around zero. Note that the Virgo baselines contribute information around 64 Hz, where the overlap reduction function of H–L is null. In general, the three baseline combinations are more sensitive given that each detector has different consistent monochromatic noise lines (e.g., the harmonics of the suspension mechanism) which are excluded from analyses, inducing gaps in the frequency spectra.

As there is no clear evidence of a GWB, as LVK places upper limits on the three fixed power law models: $\alpha =0$, typically associated with a scale-invariant, cosmological model [76]; $\alpha =2/3$, which describes a population of inspiralling compact binaries [64]; and $\alpha =3$, which corresponds to a flat strain power [151]. These are astrophysically motivated further in Section 3. With 95% confidence, it was found that ${\Omega }_{\mathrm{GW}}{|}_{25\mathrm{Hz}}<0.6\pm 1.7\times {10}^{-8}(\alpha =0)$, $0.3\pm 1.2\times {10}^{-8}(\alpha =2/3)$, and $0.4\pm 1.3\times {10}^{-9}(\alpha =3)$. In recent analyses, $\alpha$ is also allowed to vary as a parameter of the GWB model in likelihood (75). In this case, a mean-zero Gaussian prior is assumed, with standard deviation $3.5$. Marginalising over $\alpha$ with this method yields a 95% confidence upper limit of ${\Omega }_{\mathrm{GW}}{|}_{25\mathrm{Hz}}<0.7\pm 2.7\times {10}^{-8}$. When comparing the hypotheses of signal and noise to noise-only, the ${\log }_{10}$ Bayes Factor is found to be $-0.3$. The posteriors for the joint fit of $\alpha$ and ${\Omega }_{\mathrm{GW}}{|}_{25\mathrm{Hz}}$ are shown here in Figure 9.

The results discussed so far assume that the stochastic GW signal is entirely composed of the two transverse-traceless (TT) polarisation modes predicted by GR. However, modified gravity theories generally contain additional polarisations, up to a maximum of six: the two TT modes of GR, one scalar “breathing” mode, and three longitudinal modes—one scalar and two vector [247]. A detection of non-TT polarisation content would provide compelling evidence for the breakdown of GR. Using the method first developed in Ref. [160], LVKs have conducted searches for non-GR polarisations in GWB, using the distinct overlap reduction functions associated with each different set of modes (scalar, vector, and tensor) and marginalising over the power-law slope $\alpha$ as described above [56,239]. The resulting upper limits on the GWB intensity are broadly similar in each case, giving ${\Omega }_{\mathrm{T}}<6.4\times {10}^{-9}$ for tensors, ${\Omega }_{\mathrm{V}}<7.9\times {10}^{-9}$ for vectors, and ${\Omega }_{\mathrm{S}}<2.1\times {10}^{-8}$ for scalars with the O3 dataset [56].

LVK assessed future detection prospects by projecting the expected astrophysical background signal, given the current GW detections, onto the expected sensitivities of upcoming detector networks. This comparison is present in Figure 2 in Section 3.1, where we discuss the different contributions to the stochastic background. The expected signal includes contributions from BBHs, BNSs, and BHNSs, all of which have representative direct detections in the advanced LIGO-Virgo runs. We discuss details of the astrophysical background model in Section 3.1; in essence, the total GWB at the selected reference frequency is expected to be ${\Omega }_{\mathrm{GW}}^{\mathrm{CBC}}{|}_{25\mathrm{Hz}}<1.9\times {10}^{-9}$, which is about an order of magnitude away from the upper limit set with $\alpha =2/3$. Taking ${\Omega }_{\mathrm{GW}}^{\mathrm{CBC}}$ as an upper limit, it is compared with the 2$\sigma$ integrated sensitivity curve for O3,2-year projections for the Advanced LIGO-Virgo network operating at design sensitivity [218], and a 2-year, 50% duty cycle A+ network at design sensitivity [248]. Figure 2 shown in Section 3 suggests that by the time design sensitivity is reached for Advanced LIGO-Virgo detectors, these will be marginally sensitive to the background. By the A+ phase, the vast majority of the background signal should be within reach of detection.

Both the astrophysical GWB spectrum and the set of individual BBH signals resolved by LIGO/Virgo are determined by the cosmic BBH merger rate history $\dot{N}\left(z\right)$. The non-observation of the GWB from this category of sources sets an upper bound on this function, while the set of resolved signals provides a lower bound at low redshift. By combining both pieces of information, we can, therefore, infer something about this merger history [56,65], potentially revealing interesting astrophysical information about BBH abundances and formation channels, such as whether BBHs are predominantly formed through isolated binary evolution or through dynamical assembly in dense stellar environments. As was recently pointed out in Refs. [249,250,251], these inferences about the BBH merger rate at high redshift can also be used to set an upper bound on the fraction of individual detections that are expected to undergo strong gravitational lensing.At present (i.e., post-O3), the uncertainty on $\dot{N}\left(z\right)$ is still large, spanning several orders of magnitude at redshifts $z\gtrsim 1$ [56], but we can look forward to significant improvements with future LVK observing runs. The inferred redshift shape of $\dot{N}\left(z\right)$ is currently consistent with the BBH merger rate directly tracking the cosmic star formation rate (see Figure 10), but there are some tentative hints of a shallower growth at low redshift, as one would expect if there is a non-zero delay time between star formation and the eventual merger event. Note that the star formation rate reported here is arbitrarily normalised.

LVK stochastic searches have also been used to search for various cosmological sources of GWs that have been proposed in the literature, which we reviewed in Section 3. Among the most promising of these sources are cosmic strings, and LVK stochastic searches can be used in particular to place upper limits on the cosmic string tension $G\mu$. The precise constraint depends on the model adopted for the string loop network, but the strongest constraint from the first three observing runs is $G\mu <4\times {10}^{-15}$ [129]. This is the strongest existing constraint on a cosmic string model from any experiment; for comparison, CMB constraints are typically at the level of $G\mu \lesssim {10}^{-7}$ [253]. Another potential signal of cosmological origin is SGWB from PBHs. The non-observation of this signal by LIGO/Virgo implies that PBHs in the mass range $2M_{\odot }–200M_{\odot }$ cannot make up the entirety of the cosmic CDM budget [254,255], $f_{\mathrm{PBH}}\lesssim 1$; by the time LIGO/Virgo reach design sensitivity, this constraint is expected to improve to $f_{\mathrm{PBH}}\lesssim {10}^{-3}$. Searches have also been conducted for the “scalar-induced” SGWB generated by the collapsing overdensities during the process of PBH formation [165]. For GWs in the LVK frequency band, this corresponds to the formation of very light PBHs, with masses $\lesssim {10}^{13}\mathrm{kg}$. With future third-generation interferometers, these searches are expected to dominate over existing constraints coming from electromagnetic signatures of PBH evaporation.

We have aimed here to provide only a taster of some of the numerous cosmological GW sources that have been probed with LVK stochastic searches. Many other examples exist, including first-order phase transitions [164], superradiant bosonic condensates around black holes [162,163], and even sources which are not associated with GWs, but which couple directly to the LIGO/Virgo test masses to mimic a SGWB signal, such as dark photons [128].

5.1.2. Search Results for an Anisotropic Background by LVK

The anisotropic search methods outlined in Section 4.3 were most recently applied to LIGO and Virgo detector data in [211]. For the first time, data from Virgo was included in the search, leading to a three-baseline network, similarly to the isotropic search presented in Section 4.2.1. The data were pre-processed as in the latter, starting from the gated LIGO data-sets. To increase the efficiency of the search, the data were folded to one sidereal day before analysis [186]; this effectively reduces the cost of the search by a factor equal to the number of days in the dataset, with extremely little loss of information.

Three different, complementary approaches were employed in the search: A broadband, pixel-based analysis targeting point sources, referenced as the broad-band radiometer (BBR) analysis; a broadband, spherical harmonic decomposition (SHD) analysis, sensitive to extended stochastic sources [190]; and, finally, a narrowband analysis targetting three specific sky locations where a GWB signal may be expected, dubbed the narrowband radiometer (NBR) search. Notably, the so-called “radiometer” searches are both carried out in pixel space and ignore off-diagonal correlations present in the Fisher information matrix (Equation (106)) under the assumption that each direction is fully uncorrelated. The SHD search on the other hand is aimed at diffuse sources which span multiple directions and, thus, requires tackling the Fisher matrix inversion problem. The spherical harmonic domain is chosen here to reduce the effective number of degrees of freedom to a small number of modes on the sky, assuming the measurement is diffraction-limited, i.e., the maximum resolvable scale on the sky is directly related to frequency f and aperture D of an interferometric measurement by the following (see, e.g., [235] for a more detailed discussion),

(117)${\ell }_{\max }\sim \frac{2\pi Df}{c}.$

In a cross-correlation GW analysis with pairs of detectors, D is the length of the baseline for each detector pair. To obtain an estimate for a sensible ${\ell }_{\max }$ to use in the search, frequency f is assumed to be the frequency at which the detectors are most sensitive to the stochastic signal; in practice, three different spectral indices were employed in the SHD search in [211], which approach LIGO/Virgo sensitivity curves at different “peak” frequencies. These suggest a band-limiting of ${\ell }_{\max }=3,$ 4, and 16 for $\alpha =0,$ 2/3, and 3, respectively.

The BBR analysis yields sky maps consistent with mean zero Gaussian noise and, thus, provides directional upper limits, as seen in the signal-to-noise ratio (SNR) maps presented in Figure 11. Note that these are obtained by assuming each direction is fully independent and, thus, should not be read as maps at all, but rather each pixel value should be interpreted as the upper limit on stochastic GWs in that particular direction in the case where that direction is the only one contributing to a stochastic signal. Similarly, SHD analysis produces spherical harmonic coefficients for GWB, band-limited according to the prescription above, which are then converted to the pixel domain yielding maps consistent with noise. The angular power spectrum $C_{\ell }^{\mathrm{GWB}}$ as defined by Equation (47) can be directly estimated from the coefficients; however, as the measurement is noise-dominated, the unbiased estimator is [154,190]

(118)${\widehat{C^{\prime }}}_{\ell }^{\mathrm{GWB}}={\widehat{C}}_{\ell }^{\mathrm{GWB}}-{\mathcal{N}}_{\ell },{\mathcal{N}}_{\ell }=\frac{1}{2\ell +1}\sum _m{F}_{\ell m,\ell m}^{-1}.$

Here, ${\widehat{C}}_{\ell }^{\mathrm{GWB}}$ is the naive angular power spectrum estimator calculated as in Equation (47) using the spherical harmonic coefficients obtained from the data. ${\mathcal{N}}_{\ell }$ is an estimate of the noise covariance of the GW measurement, containing the appropriately normalised inverse Fisher matrix F calculated in spherical harmonic coordinates and is taken to be a measure of the uncertainty. Employing the unbiased estimator results in an angular power spectrum estimate fully consistent with zero; thus 95%, upper limits are placed on the spectrum as shown in Figure 12.

The NBR analysis described in [211] sets upper limits on frequency spectra from three astrophysically relevant sky locations: the direction of the X-ray emitting binary system Scorpius X-1 [257]; the direction of the remnant of Supernova 1987A [258]; and the center of our Galaxy, which may contain relevant GW sources [259]. No overdensity of GWs is found in any of these directions; hence, narrowband upper limits are placed in each frequency bin. Recently, LVK extended this type of search to an “All–Sky, All–Frequency” (ASAF) search [260], which produces sky maps for each frequency bin, in the pixel domain and at high pixelisation, assuming each pixel is uncorrelated. This type of search is best suited to quickly identify candidates for continuous, quasi-monochromatic GW wave sources coming from point sources on the sky, e.g., rotating neutron stars. The results of this search found a small set of outliers which may be investigated with continuous waves searches, for example, with methods proposed in [261].

5.2. Stochastic Searches with Pulsar Timing Arrays

Unlike ground-based interferometers, the first astrophysical signal observed with PTAs is expected to be a stochastic signal [120] from supermassive black hole binaries. For an insightful review of the astrophysics that can be performed with PTAs, see [262], and for a broader overview of the field, see [263]. The concept of such experiments was proposed by Sazhin [264] and further discussed by Detweiler [265]. To reach nanohertz frequencies, PTAs perform arrival time measurements from an ensemble of millisecond pulsars for a minimum of several years. The first three PTAs are the EPTA [21], NANOGrav [19], and PPTA [20]. The dependence of the Hellings–Downs curve in Equation (63) on the angular coordinates of pulsars implies that timing arrays benefit from observing pulsar pairs across the entire sky. Specifically, it was shown by Siemens et al. [266] that it is important to increase the number of observed pulsar pairs to resolve Earth-term spatial correlations of the background from the spatially uncorrelated pulsar-term component of the background, which manifests as time-correlated noise.

While EPTA and NANOGrav mostly observe pulsars from the northern hemisphere, PPTA has access to pulsars from the southern hemisphere. Recently, two more efforts have emerged: the Indian Pulsar Timing Array (InPTA [25], which is the newest member of IPTA); and the MeerTime project [267], which is based on the MeerKAT [268] facility in South Africa, a precursor of SKA [27]. A recent proposal sees the new telescope FAST [28], the largest “dish” radio telescope in the world, laying the foundations for a PTA collaboration in China. The IPTA [23] is the global effort to accelerate the detection of the nanohertz GWB by combining data sets across the previously mentioned member collaborations.

5.2.1. Search Results for an Isotropic Nanohertz Background

Searches for isotropic GWB constrain the power-law strain amplitude A of the background at the GW frequency $f=1{\mathrm{yr}}^{-1}$. Following Equation (36) and the relation between spectral indices explained in footnote ³, the strain spectrum of the background has the form $h\left(f\right)=A{(f/1{\mathrm{yr}}^{-1})}^{-2/3}$. The limits on A, denoting either confidence or credibility, are reportedly between $1.1\times {10}^{-14}$ and $1.0\times {10}^{-15}$ at $95\%$ [173,180,269,270,271,272,273,274].

As discussed in Section 4.2.2, it is expected that, prior to the detection of Hellings–Downs spatial correlations, temporal correlations with the same spectra will be observed in pulsar data sets. This phenomenon is now known as the “common-spectrum process”, and it was first reported in the search for the background with the NANOGrav 12.5-year data set [54], where the authors for the first time did not place limits on the background amplitude but instead provided the measurement of the amplitude of the common-spectrum process.

After that, the PPTA collaboration reported [26] statistical evidence for the common-spectrum process in a similar fashion, while also suggesting that current analyses do not consider a model where spectra are similar but not common. Furthermore, the authors simulated data set showed evidence for the common-spectrum process, whereas amplitudes of the injected noise power spectra were different by several orders of magnitude. Thus, it was shown that the evidence for the common-spectrum process might arise from pulsar noise described by parameters in a similar range or by model mis-specification.

On the other hand, compared with NANOGrav results [54], PPTA [26] measured the spectral slope $-\gamma$ of the common-spectrum process to be even more consistent with the theoretical prediction for the GWB from supermassive black hole binaries, $\gamma =13/3$. The relation between this value and the strain spectral index ${\alpha }^{\prime }=-2/3$ is explained in Section 4.2.2. Furthermore, under certain analysis conditions, constraints on spatial correlations as a function of pulsar-Earth baselines obtained by the PPTA seem to be even in a better agreement with the Hellings–Down curve. The PPTA measurement is presented here in Figure 6 and the NANOGrav measurement may be seen in Figure 7 of [54]. The PPTA measurement is affected by the pulsar with strong unmodelled excess noise [26]. The EPTA collaboration has also reported evidence for the common-spectrum process with a 24-year-long data set from an observation of six pulsars [275]. Measurements of power-law parameters of the common-spectrum process by NANOGrav, PPTA, and EPTA are presented in Figure 13.

The results of searches for the isotropic background are summarised in Table 1. The observation time span $T_{\mathrm{obs}}$ does not directly increase with publication year due to the inclusion or noninclusion of older data where the limited radio frequency band is not sufficient to model “chromatic” red noise arising from pulse propagation effects in the interstellar medium. One also might notice that a few upper limits on the GW strain amplitude are below the amplitude of the common-spectrum process. This is partly caused by the adoption of the older Solar System ephemerides in earlier publications, which were known to contain systematic errors. The origin of the common-spectrum process is to be clarified in future work through better constraints on spatial and temporal correlations. For a more thorough discussion about the noise with references please refer to the following subsections.

5.2.2. Challenges in GWB Searches with PTAs

PTAs, as galactic-scale gravitational wave detectors, are affected by a wide variety of noise sources. These include instrumental noise artifacts from the instruments and the hardware that processes raw telescope data; effects induced by the Solar System planets and/or the interstellar medium; and jittering of the neutron stars themselves—sources of radio waves—and their gravitationally bound companions such as other neutron stars, white dwarfs, or black holes. Both instrumental and environmental noise sources play an important role. As discussed in Section 4.2.2, stochastic processes in PTAs can be uncorrelated in time, i.e., white, and time-correlated, i.e., red. Some noise processes introduce spatial correlations between pulse arrival times in pulsars across the sky. Below we briefly discuss prospects for future PTAs given the current knowledge of the noise.

Environmental white noise associated with pulse shape changes is referred to as jitter (e.g., [277,278]). Without jitter noise, measurement errors are limited by the radio pulse signal-to-noise ratio. New scientific instruments such as FAST and SKA will significantly improve upon previous radiometer (white instrumental) noise levels, thus expanding the number of pulsars in the array and improving noise power spectra of the observed pulsars. Therefore, the timing precision of the brightest pulsars will be limited by environmental properties. In particular, MPTA reported [267] an extraordinary timing precision for some of the millisecond pulsars previously observed by the PPTA. PSR J1909–3744 showed the root-mean-square of timing residuals (rms¹⁴) of 66 ns with 11 months of observations, whereas in the second data release of PPTA [279], this pulsar showed rms residuals of 240 ns. Pulse integration times may be increased in the future to further lower white noise levels. Following Equation (81), a pulsar with an achievable timing precision of 10 ns observed biweekly, in 10 years, would be able to probe a GW strain of $5\times {10}^{-17}$ at 3 nHz. In reality, however, the detection of a nanohertz stochastic background requires timing an ensemble of pulsars, and there are sources of red noise that limit timing arrays at low frequencies.

Red noise is resolved in most pulsars after they are observed for a couple of years. A comprehensive study of this topic for the first data release of the IPTA can be found in [280]. Normally, the strongest red noise in recycled pulsars appears at low radio frequencies $\nu$ from variations in the dispersion measure, electron column density towards the line of sight. Due to dependence $P\left(f\right)\propto {\nu }^{-2}$, this noise is easy to model and it is straightforward to separate its contribution from the one of the GWB. Some red noise was found to be affecting only certain radio frequency bands or back-end observing systems. Independent of $\nu$, “achromatic” red noise is also referred to as spin noise as it is associated with pulsar rotational irregularities. Spin noise with a power spectral density which continuously increases towards low frequencies with a power-law might eventually turn over, but no such spectra were found in millisecond pulsars [281,282]. One may find the term “timing noise” in the literature, which is used to denote spin noise and glitches, i.e., sudden changes in pulsar spin periods followed by an exponential relaxation (e.g., [283]). GWBs are separated from spin noise through spatial correlations, although incorrect red noise models may result in incorrect inferences of the GWB [284]. For the EPTA second data release, it was shown [282] that tailored pulsar red noise models with band- and system-related noise marginally increase evidence for Hellings–Downs correlations, whereas PPTA showed that the inclusion of the brightest PSR J0437–4715 (which shows significant red noise and residual excess noise [285]) results in an increase in the inferred correlation coefficients $\Gamma \left({\zeta }_{ab}\right)$ [26]. Before spatial correlations are resolved, both the GWB and spin noise can yield evidence for the common-spectrum process [26]. Inferences of the GWB based on the common-spectrum process at this point in time may be treated as speculative, especially given that spin noise in pulsars may yield a power-law $P\left(f\right)$ with $\gamma =4$ [286], which is close to $\gamma =13/3=4.3\left(3\right)$ for the GWB. Some low-frequency noise can also introduce spatial correlations which can still be distinguished from the Hellings–Downs prediction [287]. Time-correlated errors in terrestrial time standards, also known as clock errors, can yield monopolar spatial correlations between pulsars a and b, $\Gamma \left({\zeta }_{ab}\right)=1$. Pulse arrival times are referenced to the position of the Solar System barycenter; thus, systematic errors in the barycenter position can yield red noise with dipolar spatial correlations $\Gamma \left({\zeta }_{ab}\right)=cos\left({\zeta }_{ab}\right)$. This noise can be mitigated byt improving Solar System ephemerides, and when improvements are not available, it can be modelled as a purely deterministic process [288,289]. In an analysis with the NANOGrav 11 year data set [180], it was demonstrated that the results are affected by the choice of Solar System ephemeris, whereas deterministic modeling of systematic errors in ephemerides renders consistent results. The good news is that the contributions of ephemeris noise are found to be rather weak for current data sets and ephemerides. In particular, in the search for a stochastic background in the second data release of the PPTA [26], the authors computed Bayesian evidence for errors in every parameter of the deterministic ephemeris model: masses and orbital Keplerian parameters of planets from Mars to Saturn, with orbital frequencies in the PTA band. They found marginal evidence for systematic errors in older ephemerides and no evidence for errors in present-day ephemerides used for GW analyses between 2020 and 2021. Red noise is mostly environmental; thus, it might be worth allocating more time to fainter pulsars that do not show evidence of such noise.

With the above in mind, we provide sensitivity projections for a few scenarios of future PTAs in Figure 14. When evaluating the performance of future PTAs, several simplifications are usually made. For example, often, red noise and effects of fitting for the timing model are neglected, as in Figure 1 and Figure 7. We addressed these points in Figure 14 by adopting the formalism from [290] and the relevant code [291]. The sensitivity estimates can be further refined based on more realistic noise properties and observing scenarios—a subject of future work. For a more detailed review of measurement errors in pulsar timing arrays please refer to [292]. Challenges for GW detection with pulsar timing arrays are also discussed in [293].

5.2.3. Search Results for an Anisotropic Nanohertz Background

To date, one search for an anisotropic GWB with PTAs has been performed [202]. In the latter, the authors use six pulsars from the EPTA with a timing baseline of $T_{\mathrm{obs}}=17.7\mathrm{years}$. While the results indicate that the data are unable to improve upon a physically motivated prior, the search setup, motivation, and results are useful to review.

The authors use a spherical harmonic decomposition with ${\ell }_{\max }=4$ based on analytical arguments of the PTA resolution established in [294]. They perform three parameterizations and one independent analysis. First, they use a single set of anisotropy coefficients across the entire frequency band (lowest five harmonics of $1/T_{\mathrm{obs}}$) and chose a fixed spectral index consistent with a population of supermassive black hole binaries. Second, they allow different anisotropy coefficients at each of the lowest five harmonics of $1/T_{\mathrm{obs}}$. Third, they allow different spherical harmonic coefficients for each of the lowest four harmonics of $1/T_{\mathrm{obs}}$ and use a single set of spherical harmonics coefficients to describe up to the fiftieth harmonic. In all cases, they impose that the power on the sky in each pixel is positive, as negative power would be unphysical. Fourth, they independently estimate the cross-correlations between individual pulsars and then map these onto a chosen basis.

The results yielded no detection of a GWB, and limits were set on the angular power spectrum $C_{\ell }$ for each of the individual parameterizations (and for each individual frequency, where appropriate). The limits on strain amplitude in multipoles with $\ell >0$ are $\lesssim 40$% of the monopole. However, these limits do not improve upon the physically motivated prior. Current astrophysical estimates suggest that the strain amplitude in multipoles with $\ell >0$ is around 20% of the monopole [295]. In the case of the fourth analysis, the physical prior is not applied, and the authors observed that sensitivity to the dipole is reduced due to the distribution of the six pulsars on the sky, highlighting the desire for using a large number of pulsars that are uniformly distributed on the sky.

5.2.4. Search Results for a Nanohertz Background Not Related to Supermassive Black Hole Binaries

While quite a few sources of a GWB that can be measured by PTAs have been proposed, outside of the generic isotropic background, which is generally assumed to be dominated by unresolved supermassive blackhole binaries, only a few explicit searches have been carried out.

Recently, both the PPTA and NANOGrav have placed constraints on first-order phase transitions in the early Universe [296,297]. In [296], the authors use 45 pulsars to search for a signal from first-order phase transitions using three different models. The common spectrum process first reported in NANOGrav data in [54] is also consistent with a first-order phase transition at temperatures below the electroweak scale, with $T^{\star }\lesssim 10\mathrm{MeV}$ at the $1\sigma$. They calculate posterior distributions on the temperature and bubble nucleation rate when phase transition occurs, as well as the strength of the phase transition and the friction between bubble walls and the plasma. However, it is important to note that these posteriors assume only a GWB from phase transitions is present, while the data cannot distinguish between this model and one described by a GWB from supermassive binary black holes. In addition, no spatial correlations between pulsars are found; hence, no detection of a GWB is claimed. The results of [297] are similar—no evidence of spatial correlations between pulsars is reported, and constraints on the temperature of a potential phase transition are found to be $1\mathrm{MeV}\lesssim T^{\star }\lesssim 100\mathrm{MeV}$. For a recent discussion on separating a background from supermassive black hole binaries and first order phase transitions (or other potential sources of a nHz GWB) using a multi-messenger approach, see [298].

Several groups have recently searched NANOGrav, PPTA, and IPTA data sets for evidence of alternative polarizations of GWs. In [299,300], the authors claim a Bayes factor of ∼100 in preference of background from spatially correlated scalar-transverse gravitational waves compared to a simple common spectrum red noise process using both NANOGrav and IPTA data. They claim similar significance when comparing a scalar transverse background to a background from tensor transverse gravitational waves. However, they found no such significant evidence when searching PPTA data [301]. In [302], the authors find a preference for a scalar transverse background, but with less significant evidence. However, this significance goes away when the authors account for uncertainties in the solar system ephemeris using [289] or when removing one pulsar, PSR J0030+0451.

5.3. Other Stochastic Background Searches

There are many GW experiments that we have omitted from the discussion above. In this section, we very briefly highlight a few that are of particular interest, some of which are shown in Figure 7. Our focus here is on the nanohertz–kilohertz frequency band; however, see the review article [303] and references therein for a summary of detection methods at MHz frequencies and above.

The oldest method of GW detection, pioneered by Weber [304,305] in the 1960s with his “resonant bar” experiments, is to monitor the resonant frequencies of some macroscopic test object. When acted upon by a GW with a frequency matching one of the object’s normal modes, the stretching and squeezing action they induce excites vibrations in these modes which can be amplified and detected. Numerous such resonant-mass experiments have been operated since the 1960s [306], but none has had sufficient sensitivity to place meaningful bounds on the GWB below the $\Omega \left(f\right)<1$ level required to prevent GWs from over-closing the Universe. However, the same idea has been applied with success on a much larger scale, using the Earth itself as a resonant mass. By monitoring the Earth’s normal mode frequencies using gravimeter data, Coughlin and Harms [220] placed an upper limit on the GWB in the megahertz frequency band. The corresponding PI curve has a minimum of $\Omega \left(f\right)\approx 1.3\times {10}^{-2}$ at $f\approx 1.6\mathrm{mHz}$. Recently, a trio of similar experiments have been proposed in which gravimeters installed on the Lunar surface could be used to monitor the Moon’s normal modes, providing a cleaner probe of GWs in the megahertz frequency band [307,308].

Another GW experiment which has successfully constrained the GWB below the $\Omega \left(f\right)<1$ level is the work of Armstrong et al. [221], who carried out Doppler-tracking observations of the Cassini spacecraft. The principle here is very similar to IFOs and PTAs in that it tracks the trajectories of photons in the presence of GWs; however, rather than measuring changes in the light-travel time, one attempts to measure changes in the photon frequency due to GW-induced Doppler shifts between the two test masses (this is essentially the time derivative of the perturbation to the light-travel time). This analysis probed GWs in the $1–100\mathsf{\mu }\mathrm{Hz}$ band, with the corresponding PI curve reaching a minimum of $\Omega \left(f\right)\approx 9.9\times {10}^{-3}$ at $f\approx 3.9\mathsf{\mu }\mathrm{Hz}$. This is a particularly interesting frequency range, as it straddles a significant gap between PTA constraints at nanohertz frequencies and future constraints from LISA and other space-based interferometers in the megahertz band; however, the Cassini constraint in this “$\mathsf{\mu }$Hz gap” is not strong enough to be of astrophysical or cosmological interest. Recently, it has been pointed out that much more competitive GWB constraints could be obtained in this frequency band using high-precision observations of binary systems (in particular, binary pulsars and the Earth–Moon system), by searching for the imprints of stochastic GWs on their orbits [225,226]. Other proposals for accessing these frequencies include an extremely long-baseline space-based interferometer such as the $\mathsf{\mu }$Ares concept [309].

In addition to changing the frequency and/or time-of-arrival of photons, GWs can also cause angular deflections in photon trajectories, such that distant astronomical objects appear displaced from their true positions on the sky. These GW-induced deflections have a distinctive geometric pattern across the sky, allowing one to search for a GWB signal with high-precision astrometry from observatories such as GAIA and THEIA [310,311,312]. The resulting PI curve lies in roughly in the same nanohertz frequency band probed by PTAs, albeit with improved high-frequency sensitivity. Forecasts in Ref. [312] indicate that next-generation astrometry with an observatory such as THEIA may be sensitive to $\Omega \left(f\right)\sim {10}^{-16}$ at nanohertz frequencies, surpassing even the reach of pulsar-timing observations with SKA.

Finally, a GW detection technique which is subject to growing interest in the community is atom interferometry, in which, rather than using interference between photons to measure the phase difference between two spatial paths, one exploits the wave-like behaviour of matter on the quantum scale to measure interference between atoms. While first-generation atom interferometry experiments such as AION [224] and MAGIS [313,314] are in their early stages and are not expected to have significant sensitivity to GWs, proposals for future kilometer-scale experiments are forecast to provide extremely useful GWB constraints in the $0.1–1\mathrm{Hz}$ band. In particular, the AION-km proposal has a forecast sensitivity of $\Omega \left(f\right)\approx 2.7\times {10}^{-12}$ at $f\approx 0.12\mathrm{Hz}$.

6. Stochastic Background Detection Prospects with Future Gravitational-Wave Detectors

We close this review with an outlook on the detection capabilities of upcoming GW detectors. We first discuss the potential of 3G detectors, which will be considerably more sensitive to GWs in the same frequency bands as LVK instruments, as may be observed in Figure 7, and specifically prospects of performing component separation with ET and CE. We then review the upcoming LISA mission and lay out different stochastic search methods which may be employed with such a detector.

6.1. Stochastic Searches with Third Generation Interferometers

By the mid-late 2030s, it is expected that ground-based gravitational wave detectors will evolve to the third generation of laser interferometers comprising ET and CE. These observatories are expected to outperform current 2G network sensitivity by several orders of magnitude in strain, extending the BBH detection horizon out to the earliest epochs of star formation at $z\approx 20$ and beyond [52,315]. This would allow 3G detectors to individually resolve more than $99.9\%$ of all stellar origin BBH signals in the Universe [316]. This eliminates the motivation to search for the background of stellar origin BBHs with standard cross-correlation methods described in Section 4 because only very few of such signals will remain unresolved, whereas searches for the non-Gaussian background outlined in Section 4.4 might be more appropriate. In between 2G and 3G detectors, there is also the proposed “2.5G” Australian detector NEMO [317], which will use very high laser power to target high frequency signals from tidal effects in the late stages of binary neutron star inspirals.

However, cross-correlation based searches might still be relevant for BNS and BHNS backgrounds because individual events will only be resolved up to redshifts $⪆1$ [318]. The approach will be then to subtract out all the resolved sources from the timestream and run a stochastic search on the remaining data. This may be performed iteratively and/or assuming different contributions relative to the signal to perform component separation [319].

Furthermore, cross-correlation based stochastic searches are suitable for searches for cosmological gravitational wave backgrounds, following a subtraction of individual compact binary coalescence signals [316,320]. The noise formed by residuals from the subtraction of best-fit compact binary coalescence waveforms can be reduced to negligible levels following the demonstration by [321].

Third-generation detectors also present a unique detection opportunity for primordial black holes. The merger rate of primordial black hole binaries per unit volume is predicted to grow continuously towards redshifts of 50 and beyond [322], producing a high redshift source of stochastic GWs which could be probed with 3G cross-correlation searches. Note that this remains subject to the theoretically uncertain rate of such events that the nature provided us with; let us cite several investigations on the matter [254,323,324,325]. Attribution of a merged black hole binary to primordial and not stellar origin based on redshift alone will not be statistically strong even with 40-kilometer Cosmic Explorer, 20-kilometer Cosmic Explorer South, and Einstein Telescope [326]. This means that even in the scenario where primordial black holes are abundant, studies of these events will require a combination of fitting for their ensemble properties as well as cross-correlation-based searches for the stochastic background that are sensitive to subthreshold signals. In the case of competing stochastic backgrounds, studies have been performed to carry out simultaneous Bayesian fitting of both astrophysical and primordial components [168,327], and methods have been proposed to distinguish primordial black holes from astrophysical ones [328,329].

Correlating anisotropies in the GWB which electromagnetic observations has also been proposed as a target for future detectors. In some cases, this means using galaxy catalogues to guide the model used for anisotropic GWB searches [241], while in other cases, the proposal is to correlate sub-threshold GW signals with EM transients to probe fundamental physics [127,330].

There are quite a few technological challenges in the development of 3G detectors that are directly relevant for GWB searches, for which its sensitivity to ${\Omega }_{\mathrm{gw}}$ scales with $f^{-3}$ and, therefore, benefit from improved low frequency sensitivity. Below, we discuss two important and probable sources of low frequency noise that are important for GWB searches.

First, globally correlated magnetic noise from, e.g., Schumann resonances [331,332] will, if not subtracted or dealt with in some manner, provide a lower limit on our ability to measure a GWB in the frequency band from ∼8–40 Hz [167,333,334,335,336,337,338,339,340,341], and could even be an issue at higher frequencies [341]. The assumption that a GWB will look like correlated noise between multiple detectors means that correlated magnetic noise could masquerade as a GWB if we are not careful. There are several avenues we can take to avoid this issue. One of the most promising is Wiener filtering using low noise magnetometers near the detectors [334,335,336,338]. Another option is spectral separation using the hybrid-Bayesian method discussed at the end of Section 4.2, where one models the contribution of the correlated magnetic noise to the cross-correlation statistic and includes it with a model for the GWB [167]. Finally, reducing the coupling of magnetic fields into the strain channel of the detectors will also be vital [341]. In the end, some combination of all three of these ideas will likely be needed.

Second, low frequency noise due to gravity gradients caused by temperature fluctuations in the atmosphere (infrasound) or seismic waves, known as “Newtonian noise,” is likely a limiting noise source at frequencies up to ∼30 Hz [14,18,342,343,344,345]. For a complete review, see [344]. One cannot shield against Newtonian noise and, therefore, need to either subtract [346,347,348,349,350] it or cleverly design the local environment of the detectors to reduce surface seismic waves and atmospheric perturbations [351]. Estimating the specific contribution of Newtonian noise to the detectors from, e.g., seismic waves, is calculable based on a realisation of the seismic wavefield [344], but is very difficult to achieve in practice. Proof of principle tests of using an array of seismometers or geophones to create Wiener filters or use machine learning to clean a “target” seismometer channel are promising, and work on optimal design of local sensor arrays is mature [346,347,348,349,350,352]. However, cleaning a seismometer is quite different from cleaning Newtonian noise from a gravitational wave detector. In the end, it is likely that for 3G detectors, an array of geophones, broadband seismometers, and tiltmeters (which can be used to separate seismic tilt from seismic waves) will be needed to subtract Newtonian noise due to seismic waves [344,353]. Meanwhile, an array of infrasound microphones or LIDAR detectors, along with potentially using wind fences around facilities, will be needed to reduce and cancel Newtonian noise due to atmospheric fluctuations [354,355].

6.2. Stochastic Searches with the Laser Interferometer Space Antenna

Stochastic searches with LISA present distinct challenges from those discussed for PTAs and ground-based interferometer networks, primarily due to the fact that LISA is considered to be a single detector; hence, it is not possible to use a cross-correlation search with a second detector which is not co-located with LISA. It is, however, possible to auto- and cross-correlate TDI channels, as we explain here.

The antenna itself will be comprises three spacecrafts which will be positioned on three heliocentric orbits and will specify a plane at a 60${}^{\circ }$ angle with the ecliptic, in an equilateral triangle configuration which will be maintained throughout the entire duration of the mission (4+ years) [29]. Each spacecraft will send a laser beam to each other spacecraft such that, in total, six laser measurements will be made, comprising the links of the detector. The links are then combined to form TDI channels, as detailed in Section 4.1. Basic LISA studies start from considering sets of three channels, where each channel is centered around one spacecraft; however, more sophisticated studies are already under way [138,140,356]. These channels produce a set of observations, which we can treat as timestreams for the purposes of this discussion (although these are effectively combined in post-processing and are not the read-out of the detector), $X\left(t\right),Y\left(t\right),Z\left(t\right)$, as introduced in Section 4.1. Taking the same approach as for real interferometer timestreams, we can decompose each of these in their signal and noise components,

(119)$X\left(t\right)=R_X\left(t\right)\ast h\left(t\right)+n_X\left(t\right),$

where $R_x$ is the specific response of the X channel, and $n_X$ is its noise term. The noise terms for the three channels are not independent, as each shares (possibly multiple) links; hence, a residuals analysis such as that based on Likelihood (73) is out of the question. Hence, LISA stochastic analyses are based on Gaussian likelihoods for the data as Likelihood (72), where the data are now a set of TDI channels,

(120)$\mathcal{L}\left(\mathit{d}\right)\propto \prod _{f,\tau }\frac{1}{{\left|\mathit{C}\right|}^{1/2}}{e}^{\frac{1}{2}\mathit{T}{\mathit{C}}^{-1}{\mathit{T}}^{\star }},$

with $\mathit{T}={(X,Y,Z)}^T$, and $\mathit{C}$ is the full covariance matrix of the data, $\mathit{C}=⟨\mathit{T}\otimes {\mathit{T}}^{\star }⟩$. The key to detecting a stochastic background with LISA is then understanding how instrument noise and GW signals manifest in different TDI channels and finding an accurate modelling approach to plug into Likelihood (120). The first step consists in picking the specific TDI channels and assumptions to work in—there are several options provided in the literature [134,135,146,357,358,359].

Both maximum-likelihood methods and Bayesian approaches which compute the full posteriors have been proposed to extract and characterise an SGWB from LISA data. In the first case, it is possible to produce a maximum-likelihood solution by fixing a noise model and iterating over the data until the signal converges to the maximum-likelihood estimate; this is essentially the same approach as that laid out in Section 4.3 for anisotropic searches, based on [193]. Otherwise, one may take a full Bayesian approach, parametrizing signal and noise components in the data and assuming appropriate priors for the two sets of parameters, then sampling the Likelihood to obtain informative posteriors. This is performed for example in [147] for an isotropic stochastic background, where the data were given by the LISA Data Challenge (LDC). Here, the authors show that, using the noise orthogonal $A,E,T$ set of TDI channels, they are able to recover the stochastic signal and independently measure the relevant LISA noise parameters for each arm of the interferometer. Note, however, that this is performed in the absence of competing signals.

In fact, component separation will be a major challenge for LISA data analysis, as a variety of signals will contribute to the same frequency band and overlap in time [360]. In the case of stochastic backgrounds, as detailed in Section 3, LISA will be sensitive to an anisotropic galactic white dwarf binary background which traces the shape of the Milky Way, a mostly isotropic background of either primordial or stellar-origin black hole and neutron star binaries [361], and also potentially other primordial backgrounds, such as those arising from inflation [93], first-order phase transitions [55,362], and cosmic strings [53,363], all of which may be either isotropic or anisotropic. Galactic binaries will by far dominate the measurement, and in fact several approaches have been proposed towards component separation, relying on characteristic, observable differences of each component such as their spectral shape [148,364], or their distribution on the sky, as seen in [193,194]. In any case, it will first be necessary to disentangle the resolvable white dwarf binaries from the unresolvable white dwarf binary background in the data streams, as detailed in [132].

Anisotropic stochastic searches with LISA leverage the motion of the LISA constellation around the sun to map the signal in the solar system barycenter frame. The method proposed in [193] fixes the noise model, simulates mock data assuming the LISA response functions in [365] and uses the Bond–Jaffe–Knox maximum-likelihood solution in Equation (111) to test the mapping capabilities of the detector. It was found that the angular resolution of LISA for stochastic signals highly depends on their SNR as a function of frequency, as the detector response varies greatly as a function of the latter. The stochastic signals compete with the detector noise, which at low SNR dominates the modes of the Fisher information matrix (112). Beyond this, the angular resolution at which LISA “sees” the signal is diffraction-limited; hence, the higher the frequency, the better the angular resolution will be. In the best-case scenario of a high SNR signal dominating at frequencies ∼0.1 Hz, the angular resolution will be of the order ${\ell }_{\max }\sim 16$—see Figure 6 of [193]. Another approach proposed in [194] tackles the mapping problem by setting priors on the spherical harmonic modes of the signal and solving for their amplitudes in a Bayesian framework, using a single TDI channel. This is especially useful when the sky distribution of the signal is well known, as is the case for the galactic white dwarf binary background. The technique proposed here is based on the Clebsch–Gordan expansion, which allows one to constrain the spherical harmonics with a non-negative distribution. This method has proven to be effective for the recovery of a mock galactic background and yields an informative set of posteriors on the coefficients of the signal—see Figure 3 of [194]. However, it has only been tested at the low resolution of ${\ell }_{\max }=4$ as the solver scales significantly in computational cost with the number of modes one tries to recover.

Author Contributions

Conceptualization: A.I.R., B.G., A.C.J. and P.M.M. Data curation: A.I.R., B.G., A.C.J. and P.M.M. Software: A.I.R., B.G., A.C.J. and P.M.M. Making of figures and tables: A.I.R., B.G. and A.C.J. Writing—original draft: A.I.R., B.G., A.C.J. and P.M.M. Writing—review & editing: A.I.R., B.G., A.C.J. and P.M.M. All authors have read and agreed to the published version of the manuscript.

Funding

AIR acknowledges the support of the National Science Foundation and the LIGO Laboratory. BG is supported by the Italian Ministry of Education, University and Research within the PRIN 2017 Research Program Framework, n. 2017SYRTCN. PMM was supported by the NANOGrav Physics Frontiers Center, National Science Foundation (NSF), award number 2020265. Parts of this research were conducted by the Australian Research Council Centre of Excellence for Gravitational Wave Discovery (OzGrav), through project number CE170100004.

Data Availability Statement

This manuscript shows figures created using the LVK open data [246,256], and publicly available PTA results presented in [26,54,275,285].

Acknowledgments

The authors thank Katerina Chatziioannou and Joseph Romano for precious comments on our original draft.

Conflicts of Interest

The authors declare no conflict of interest.

Footnotes

Notes

1. There are several theories that postulate modifications to this statement, but these modifications are strongly constrained by multimessenger observations of the binary neutron star merger GW170817 [36].

2. Strictly speaking, GWs with wavelengths larger than the size of the cosmological horizon are frozen out by Hubble friction and, thus, do not contribute to the effective energy density. As a result, the density parameter appearing on the left-hand side of Equation (14) should be interpreted as an integral over ${\Omega }_{\mathrm{GW}}\left(f\right)$ from a minimum frequency of $f\sim {10}^{-15}\mathrm{Hz}$, which corresponds to the size of the horizon at the epoch when the CMB photons were emitted.

3. Often in the literature, the letter $\alpha$ is used to denote a different spectral index, hereby denoted ${\alpha }^{\prime }$, of the characteristic strain spectrum of the stochastic gravitational wave background $h\left(f\right)\equiv \sqrt{fS\left(f\right)}\propto f^{{\alpha }^{\prime }}$, where $S\left(f\right)$ is the power spectral density of the background. This leads to the power-law index of $2{\alpha }^{\prime }+2$ for ${\Omega }_{\mathrm{GW}}$ in Equation (25) instead of $\alpha$.

4. There is still much uncertainty around this black hole “lower mass gap” and it remains a main focus of research in the field; see, e.g., the discussions in [58,59,60].

5. As mentioned in the footnote 3, the energy density spectral index $\alpha =2/3$ is related to the strain amplitude spectral index ${\alpha }^{\prime }$ such that $\alpha =2{\alpha }^{\prime }+2$. Therefore ${\alpha }^{\prime }=-2/3$, which is the quantity often referred to as $\alpha$ in publications. Due to the parametric choices made in the different search pipelines, $\alpha$ is most frequently used in LVK literature, while ${\alpha }^{\prime }$ is most frequently used in the PTA literature.

6. There are, in fact, sources of spatially correlated noise. They can be distinguished from SGWBs by either a deterministic component or by an overlap reduction function different to those of SGWBs. We provide more details in Section 5.

7. Time segments with non-stationary noise are removed from analyses of real data [56], the stationarity is usually determined empirically. Note that when applying a non-trivial windowing function to time-domain data, e.g., for computing a discrete Fourier transform, we introduce correlations between frequency bins. This effect and the methods to mitigate it are described in [130].

8. Note that this is result is derived from the detector response to a fully unpolarised stochastic background; in the case of a polarised signal, extra terms need to be considered. More details on this can be found, e.g., in [136,137].

9. Note that directly interfering laser beams in the case of LISA is impossible due to energy dispersion along the arm [29].

10. Note that this is no longer the case in the presence of temporal shot noise, e.g., for the astrophysical GWB from compact binary coalescences [122]. In this case, each time segment will have random GWB intensity fluctuations due to the finite number of sources. The statistical independence of these fluctuations at different times can be leveraged to mitigate the impact of shot noise on measurements of the angular power spectrum [123].

11. Note that in these two cases $I_p$ and $I_{\ell m}$ have different units because in (114) the basis function carry units ${\mathrm{sr}}^{-1}$.

12. In fact, in [215] the authors discuss the extension of their optimal search method in the presence of non-Gaussian noise as well. The implementation of the method becomes substantially more involved, however practical considerations about the nature and behaviour of the non-Gaussian noise can simplify it considerably.

13. Initial-detector era searches on LIGO Science Run 5 (S5)/Virgo Science Run 1 (VSR1) [16] and LIGO Science Run 6, Virgo Science Runs 2 and 3 (VSR23) [231] used the full Hanford, Livingston, Virgo detector network.

14. Often used in the literature to quantify noise in a given pulsar. For rms ${\sigma }^2$ of Gaussian white noise, the power spectral density of residuals $P=2\pi {\sigma }^2\Delta t$, where $\Delta t$ is the time between observations, which is typically a couple of weeks. Note, rms residuals may increase with longer data spans due to more sources of noise becoming prominent at longer time scales. Sometimes rms residuals are provided for a particular source of noise.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Table

Figure 1. An overview of potential GWB signals across the frequency spectrum. The light blue curve shows the prediction for single-field slow-roll inflation with a canonical kinetic term, with tensor-to-scalar ratio [Forumla omitted. See PDF.] [52]. The pink curve shows a GWB from Nambu–Goto cosmic strings, using “model 2” of the loop network, with a dimensionless string tension of [Forumla omitted. See PDF.] [53]. The brown curve shows a GWB from inspiralling supermassive BBHs, with the amplitude and shaded region shown here corresponding to the common noise process in the NANOGrav 12.5-year data set [54]. The two grey curves show GWBs generated by first-order phase transitions at the electroweak scale (∼200 GeV) and the QCD scale (∼200 MeV), respectively [55]. The yellow curve shows a GWB generated by stellar-mass compact binaries, based on the mass distributions and local merger rates inferred by LVK detections [56]. The dashed curves show various observational constraints, as described further in Section 5 (this includes the PPTA constraint, which intersects the possible NANOGrav SMBBH signal); the dotted curve shows the integrated constraint from measurements of [Forumla omitted. See PDF.], which cannot be directly compared with the frequency-dependent constraint curves but is shown here for indicative purposes.

Figure 2. Projections for the CBC stochastic background using the third GW transient catalogue (GWTC-3), first presented in [74]. On the right, the total expected background is compared to the integrated sensitivity of current and future configurations of the 2G detector network. This plot was obtained using open data published in [75].

Figure 3. Examples of interferometer timestreams populated by GWBs with different time-domain properties, inspired by [64]: Poisson noise (top panel), popcorn (middle panel), and Gaussian (bottom panel). The strain on the y-axis is normalised to show the fractional contribution of each burst/wavelet to the signal, while the [Forumla omitted. See PDF.] parameter is the duty cycle, i.e., the mean ratio between event duration and the time interval between consecutive events.

Figure 4. Overlap reduction functions for different GW interferometer detectors: ground-based (left) and space-based (right). Note that, in panel (b), the positive and negative parts of the LISA TDI overlap function [Forumla omitted. See PDF.] are plotted separately, as the y-axis is log-scaled.

Figure 5. Geometric components of the overlap functions for different ground-based detector combinations mapped on the sky in geocentric celestial coordinates at a common time.

Figure 6. The dashed line is the Hellings–Downs curve described by Equation (63). [Forumla omitted. See PDF.] is the angle between Earth–pulsars baselines, and [Forumla omitted. See PDF.] measures the spatial correlation. At [Forumla omitted. See PDF.], the function assumes the pulsars located at different distances. When distances are the same, [Forumla omitted. See PDF.]. Shaded regions correspond to constraints on spatial correlations obtained with the second data release of the Parkes Pulsar Timing Array except PSR J0437–4715 [26].

Figure 7. A survey of constraints (all at 95% confidence) on the GWB across the frequency spectrum. Solid curves indicate existing results from LIGO/Virgo’s first three observing runs [56], monitoring of the Earth’s normal modes [220], Doppler tracking of the Cassini satellite [221], pulsar timing observations by the PPTA [222], and CMB temperature and polarisation spectra measured by Planck [222,223]. Dashed curves are forecast constraints for LIGO at A+ sensitivity, Einstein Telescope [14], AION-km [224], LISA [29], binary resonance searches [225,226], and pulsar timing with the Square Kilometre Array [227]. The dotted curve indicates the level of the integrated constraint from measurements of [Forumla omitted. See PDF.] [38]; note that this is a constraint on the total GW energy density over a broad frequency range and cannot be directly compared to the other constraints. Note also that both the Planck and [Forumla omitted. See PDF.] constraints apply only to primordial GWs emitted before the epoch of BBN. See Figure 1 for various GWB signal predictions in relation to these constraints.

Figure 8. Combined cross-correlation spectra for the three O3 baselines. The grey lines here mark [Forumla omitted. See PDF.] contours. The data, by eye, appears to follow the Gaussian distribution. Note that adding the Virgo baselines has filled in the gap in the H–L data present around the 60 Hz power line (this broad line is present in both USA-based detectors, but not in Italy). The original version of this plot is presented in [56]; the one shown here was obtained using open data published in [246].

Figure 9. Posterior distributions for the GWB energy density [Forumla omitted. See PDF.] at reference frequency [Forumla omitted. See PDF.] Hz and spectral index [Forumla omitted. See PDF.]. The dashed grey lines show the priors used in parameter estimation. Note that a log-uniform prior was chosen here on [Forumla omitted. See PDF.], and a mean–zero Gaussian was chosen for [Forumla omitted. See PDF.] as the detector noise PSD is approximately flat. The original version of this plot is presented in [56]; the one shown here was obtained using the open data published in [246].

Figure 10. Collection of all posterior draws for the merger rate as a function of redshift [Forumla omitted. See PDF.] inferred from the combination of the resolved BBH detections in the GWTC-2 catalog [252] and the upper limits on the stochastic background obtained with LIGO/Virgo O3 data. The original version of this plot is presented in [56]; the one shown here was obtained using open data published in [246].

Figure 11. SNR maps obtained with the “broad–band radiometer” search method applied to O1, O2, and O3 data, for different fixed values of the spectral index [Forumla omitted. See PDF.], as explained in the text. It is clear there is no significant high SNR signal in the maps. The original version of this plot is presented in [211]; the one shown here was obtained using open data published in [256].

Figure 12. Angular power spectrum upper limits obtained with the “spherical harmonic decomposition” search method applied to O1, O2, and O3 data, for different fixed values of the spectral index [Forumla omitted. See PDF.], as explained in the text. The original version of this plot is presented in [211]; the one shown here was obtained by using open data published in [256].

Figure 13. Power-law amplitude A and the negative of the spectral index, [Forumla omitted. See PDF.] (see Equation (80)), of the common component of power spectral density of timing residuals [Forumla omitted. See PDF.] as reported NANOGrav (orange [54]), PPTA (blue [26]), and EPTA (green [275]).

Figure 14. Power-law-integrated sensitivity [290,291] at the level of SNR = 3 of simulated PTAs to the SGWB strain [Forumla omitted. See PDF.]. We demonstrate the effects of timing precision, red noise, and a number of array pulsars. Pulsars are distributed isotropically across the sky and observed biweekly. Orange lines represent variations of the 15-year-long PPTA DR2, with real pulsar [Forumla omitted. See PDF.], white noise levels based on the rms residuals ranging from 240 ns to 24,050 ns [279] and best-fit red noise parameters from [285]. The dashed line represents PPTA DR2 extended by 5 more years of observation, whereas the dotted line represents PPTA DR2 with a doubled number of pulsars with the same noise properties. Red line is an example of the future IPTA with 20 years of observations of 100 pulsars with white noise rms drawn uniformly between 20 and 500 ns. Red noise [Forumla omitted. See PDF.] is drawn uniformly from [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] drawn uniformly from [Forumla omitted. See PDF.], a fiducial choice based on [285]. Green lines represent SKA observing 50 pulsars with white noise rms of 30 ns for 15 years, with (solid line) and without (dashed line) red noise assumed for the IPTA.

Table 1

Results of searches for the nanohertz-frequency isotropic GWB with the amplitude A and the energy-density spectral index $\alpha =2/3$ in chronological order. The first three columns list details of corresponding publications. Publications that adopted alphabetical author list are denoted by asterisks. The fourth column contains 95% upper limits on A for the GWB and measurements of A of the common-spectrum process (CP) with credible levels of $1\sigma$ in [26] and 5–95% in [54,275,276]. The values in the last two columns are the characteristics of the data set: the total observation time $T_{\mathrm{obs}}$ and the number of pulsars in the array $N_{\mathrm{psr}}$.

Publication	Collaboration	Year	$\mathit{A}(\times {10}^{-15})$	${\mathit{T}}_{\mathbf{obs}}\left(\mathbf{yr}\right)$	${\mathit{N}}_{\mathbf{psr}}$
Jenet et al. [269]	PPTA	2006	<11	20	7
van Haasteren et al. [173]	EPTA	2011	<6	11	5
Demorest et al. [270]	NANOGrav	2013	<7	5	17
Shannon et al. [271]	PPTA	2013	<2.4	25	6
Lentati et al. [272]	EPTA	2015	<3.0	18	6
* Arzoumanian et al. [178]	NANOGrav	2015	<1.5	9	37
Shannon et al. [274]	PPTA	2015	<1.0	11	4
Verbiest et al. [24]	IPTA	2016	<1.7	21	4
* Arzoumanian et al. [180]	NANOGrav	2018	<1.45	11.4	45
* Arzoumanian et al. [54]	NANOGrav	2020	CP: $1.9_{-0.6}^{+0.8}$	12.5	45
Goncharov et al. [26]	PPTA	2021	CP: $2.2_{-0.3}^{+0.4}$	15.0	26
Chen et al. [275]	EPTA	2021	CP: $3.0_{-0.7}^{+0.9}$	24.0	6
* Antoniadis et al. [276]	IPTA	2022	CP: $2.8_{-0.8}^{+1.2}$	30.2	53

References

1. Davis, D.; Walker, M. Detector Characterization and Mitigation of Noise in Ground-Based Gravitational-Wave Interferometers. Galaxies; 2022; 10, 12. [DOI: https://dx.doi.org/10.3390/galaxies10010012]

2. Sun, L.; Goetz, E.; Kissel, J.S.; Betzwieser, J.; Karki, S.; Viets, A.; Wade, M.; Bhattacharjee, D.; Bossilkov, V.; Covas, P.B. et al. Characterization of systematic error in Advanced LIGO calibration. Class. Quant. Grav.; 2020; 37, 225008. [DOI: https://dx.doi.org/10.1088/1361-6382/abb14e]

3. Sun, L.; Goetz, E.; Kissel, J.S.; Betzwieser, J.; Karki, S.; Bhattacharjee, D.; Covas, P.B.; Datrier, L.E.; Kandhasamy, S.; Lecoeuche, Y.K. et al. Characterization of systematic error in Advanced LIGO calibration in the second half of O3. arXiv; 2021; arXiv: 2107.00129

4. Acernese, F.; Agathos, M.; Ain, A.; Albanesi, S.; Allocca, A.; Amato, A.; Andrade, T.; Andres, N.; Andrić, T.; Ansoldi, S. et al. Calibration of Advanced Virgo and reconstruction of detector strain h(t) during the Observing Run O3. Class. Quant. Grav.; 2022; 39, 045006. [DOI: https://dx.doi.org/10.1088/1361-6382/ac3c8e]

5. Akutsu, T.; Ando, M.; Arai, K.; Arai, Y.; Araki, S.; Araya, A.; Aritomi, N.; Asada, H.; Aso, Y.; Bae, S. et al. Overview of KAGRA: Calibration, detector characterization, physical environmental monitors, and the geophysics interferometer. PTEP; 2021; 2021, 05A102. [DOI: https://dx.doi.org/10.1093/ptep/ptab018]

6. Verbiest, J.P.W.; Osłowski, S.; Burke-Spolaor, S. Pulsar Timing Array Experiments. Handbook of Gravitational Wave Astronomy; Bambi, C.; Katsanevas, S.; Kokkotas, K.D. Springer: Singapore, 2020; pp. 1-42. [DOI: https://dx.doi.org/10.1007/978-981-15-4702-7_4-1]

7. Tiburzi, C. Pulsars Probe the Low-Frequency Gravitational Sky: Pulsar Timing Arrays Basics and Recent Results. PASA; 2018; 35, e013. [DOI: https://dx.doi.org/10.1017/pasa.2018.7]

8. Hobbs, G.; Dai, S. Gravitational wave research using pulsar timing arrays. Natl. Sci. Rev.; 2017; 4, pp. 707-717. [DOI: https://dx.doi.org/10.1093/nsr/nwx126]

9. Abbott, B.P.; Abbott, R.; Abbott, T.D.; Abernathy, M.R.; Acernese, F.; Ackley, K.; Adams, C.; Adams, T.; Addesso, P.; Adhikari, R.X. et al. GW150914: The Advanced LIGO Detectors in the Era of First Discoveries. Phys. Rev. Lett.; 2016; 116, 131103. [DOI: https://dx.doi.org/10.1103/PhysRevLett.116.131103]

10. Accadia, T.; Acernese, F.; Antonucci, F.; Astone, P.; Ballardin, G.; Barone, F.; Barsuglia, M.; Basti, A.; Bauer, T.S.; Bebronne, M. et al. Status of the Virgo project. Class. Quantum Gravity; 2011; 28, 114002. [DOI: https://dx.doi.org/10.1088/0264-9381/28/11/114002]

11. Affeldt, C.; Danzmann, K.; Dooley, K.L.; Grote, H.; Hewitson, M.; Hild, S.; Hough, J.; Leong, J.; Lück, H.; Prijatelj, M. et al. Advanced techniques in GEO 600. Class. Quantum Gravity; 2014; 31, 224002. [DOI: https://dx.doi.org/10.1088/0264-9381/31/22/224002]

12. Castelvecchi, D. Japan’s pioneering detector set to join hunt for gravitational waves. Nature; 2019; 562, pp. 9-10. [DOI: https://dx.doi.org/10.1038/d41586-018-07867-z] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30602755]

13. LIGO-India, Proposal of the Consortium for Indian Initiative in Gravitational-Wave Observations (IndIGO). Available online: https://dcc.ligo.org/LIGO-M1100296/public (accessed on 1 November 2021).

14. Punturo, M.; Abernathy, M.; Acernese, F.; Allen, B.; Andersson, N.; Arun, K.; Barone, F.; Barr, B.; Barsuglia, M.; Beker, M. et al. The Einstein Telescope: A third-generation gravitational wave observatory. Class. Quant. Grav.; 2010; 27, 194002. [DOI: https://dx.doi.org/10.1088/0264-9381/27/19/194002]

15. Sathyaprakash, B.; Abernathy, M.; Acernese, F.; Amaro-Seoane, P.; Andersson, N.; Arun, K.; Barone, F.; Barr, B.; Barsuglia, M.; Beker, M. et al. Scientific Potential of Einstein Telescope. Proceedings of the 46th Rencontres de Moriond on Gravitational Waves and Experimental Gravity; La Thuile, Italy, 20–27 March 2011; pp. 127-136.

16. Abadie, J.; Abbott, B.P.; Abbott, R.; Abbott, T.D.; Abernathy, M.; Accadia, T.; Acernese, F.; Adams, C.; Adhikari, R.; Affeldt, C. et al. Upper limits on a stochastic gravitational-wave background using LIGO and Virgo interferometers at 600–1000 Hz. Phys. Rev. D; 2012; 85, 122001. [DOI: https://dx.doi.org/10.1103/PhysRevD.85.122001]

17. Reitze, D.; Adhikari, R.X.; Ballmer, S.; Barish, B.; Barsotti, L.; Billingsley, G.; Brown, D.A.; Chen, Y.; Coyne, D.; Eisenstein, R. et al. Cosmic Explorer: The U.S. Contribution to Gravitational-Wave Astronomy beyond LIGO. Bull. Am. Astron. Soc.; 2019; 51, 035.

18. Hall, E.D.; Kuns, K.; Smith, J.R.; Bai, Y.; Wipf, C.; Biscans, S.; Adhikari, R.X.; Arai, K.; Ballmer, S.; Barsotti, L. et al. Gravitational-wave physics with Cosmic Explorer: Limits to low-frequency sensitivity. Phys. Rev. D; 2021; 103, 122004. [DOI: https://dx.doi.org/10.1103/PhysRevD.103.122004]

19. McLaughlin, M.A. The North American Nanohertz Observatory for Gravitational Waves. Class. Quantum Gravity; 2013; 30, 224008. [DOI: https://dx.doi.org/10.1088/0264-9381/30/22/224008]

20. Manchester, R.N.; Hobbs, G.; Bailes, M.; Coles, W.A.; van Straten, W.; Keith, M.J.; Shannon, R.M.; Bhat, N.D.R.; Brown, A.; Burke-Spolaor, S.G. et al. The Parkes Pulsar Timing Array Project. PASA; 2013; 30, e017. [DOI: https://dx.doi.org/10.1017/pasa.2012.017]

21. Kramer, M.; Champion, D.J. The European Pulsar Timing Array and the Large European Array for Pulsars. Class. Quantum Gravity; 2013; 30, 224009. [DOI: https://dx.doi.org/10.1088/0264-9381/30/22/224009]

22. Desvignes, G.; Caballero, R.N.; Lentati, L.; Verbiest, J.P.W.; Champion, D.J.; Stappers, B.W.; Janssen, G.H.; Lazarus, P.; Osłowski, S.; Babak, S. et al. High-precision timing of 42 millisecond pulsars with the European Pulsar Timing Array. Mon. Not. R. Astron. Soc.; 2016; 458, pp. 3341-3380. [DOI: https://dx.doi.org/10.1093/mnras/stw483]

23. Hobbs, G.; Archibald, A.; Arzoumanian, Z.; Backer, D.; Bailes, M.; Bhat, N.D.R.; Burgay, M.; Burke-Spolaor, S.; Champion, D.; Cognard, I. et al. The International Pulsar Timing Array project: Using pulsars as a gravitational wave detector. Class. Quantum Gravity; 2010; 27, 084013. [DOI: https://dx.doi.org/10.1088/0264-9381/27/8/084013]

24. Verbiest, J.P.W.; Lentati, L.; Hobbs, G.; van Haasteren, R.; Demorest, P.B.; Janssen, G.H.; Wang, J.B.; Desvignes, G.; Caballero, R.N.; Keith, M.J. et al. The International Pulsar Timing Array: First data release. MNRAS; 2016; 458, pp. 1267-1288. [DOI: https://dx.doi.org/10.1093/mnras/stw347]

25. Joshi, B.C.; Arumugasamy, P.; Bagchi, M.; Bandyopadhyay, D.; Basu, A.; Dhanda Batra, N.; Bethapudi, S.; Choudhary, A.; De, K.; Dey, L. et al. Precision pulsar timing with the ORT and the GMRT and its applications in pulsar astrophysics. J. Astrophys. Astron.; 2018; 39, 51. [DOI: https://dx.doi.org/10.1007/s12036-018-9549-y]

26. Goncharov, B.; Shannon, R.M.; Reardon, D.J.; Hobbs, G.; Zic, A.; Bailes, M.; Curyło, M.; Dai, S.; Kerr, M.; Lower, M.E. et al. On the Evidence for a Common-spectrum Process in the Search for the Nanohertz Gravitational-wave Background with the Parkes Pulsar Timing Array. ApJ; 2021; 917, L19. [DOI: https://dx.doi.org/10.3847/2041-8213/ac17f4]

27. Dewdney, P.E.; Hall, P.J.; Schilizzi, R.T.; Lazio, T.J.L.W. The Square Kilometre Array. IEEE Proc.; 2009; 97, pp. 1482-1496. [DOI: https://dx.doi.org/10.1109/JPROC.2009.2021005]

28. Nan, R.; Li, D.; Jin, C.; Wang, Q.; Zhu, L.; Zhu, W.; Zhang, H.; Yue, Y.; Qian, L. The Five-Hundred Aperture Spherical Radio Telescope (fast) Project. Int. J. Mod. Phys. D; 2011; 20, pp. 989-1024. [DOI: https://dx.doi.org/10.1142/S0218271811019335]

29. Amaro-Seoane, P.; Audley, H.; Babak, S.; Baker, J.; Barausse, E.; Bender, P.; Berti, E.; Binetruy, P.; Born, M.; Bortoluzzi, D. et al. Laser Interferometer Space Antenna. arXiv; 2017; arXiv: 1702.00786

30. Wald, R.M. General Relativity; University of Chicago Press: Chicago, IL, USA, 1984.

31. Cornish, N.J. Mapping the gravitational-wave background. Class. Quantum Gravity; 2001; 18, pp. 4277-4291. [DOI: https://dx.doi.org/10.1088/0264-9381/18/20/307]

32. Allen, B.; Ottewill, A.C. Detection of Anisotropies in the Gravitational-Wave Stochastic Background. Phys. Rev. D; 1996; 56, pp. 545-563. [DOI: https://dx.doi.org/10.1103/PhysRevD.56.545]

33. Romano, J.D.; Cornish, N.J. Detection methods for stochastic gravitational-wave backgrounds: A unified treatment. Living Rev. Relativ.; 2017; 20, pp. 1-223. [DOI: https://dx.doi.org/10.1007/s41114-017-0004-1]

34. Jackson, J.D. Classical Electrodynamics; 3rd ed. Wiley: New York, NY, USA, 1999.

35. Renzini, A.I. Mapping the Gravitational-Wave Background. Ph.D. Thesis; Imperial College London: London, UK, 2020; [DOI: https://dx.doi.org/10.25560/83803]

36. Abbott, B.P.; Abbott, R.; Abbott, T.D.; Acernese, F.; Ackley, K.; Adams, C.; Adams, T.; Addesso, P.; Adhikari, R.X.; Adya, V.B. et al. Gravitational Waves and Gamma-rays from a Binary Neutron Star Merger: GW170817 and GRB 170817A. Astrophys. J. Lett.; 2017; 848, L13. [DOI: https://dx.doi.org/10.3847/2041-8213/aa920c]

37. Kolb, E.W.; Turner, M.S. The Early Universe; Addison-Wesley: Boston, MA, USA, 1990; Volume 69.

38. Pagano, L.; Salvati, L.; Melchiorri, A. New constraints on primordial gravitational waves from Planck 2015. Phys. Lett. B; 2016; 760, pp. 823-825. [DOI: https://dx.doi.org/10.1016/j.physletb.2016.07.078]

39. Christensen, N. Stochastic Gravitational Wave Backgrounds. Rep. Prog. Phys.; 2019; 82, 016903. [DOI: https://dx.doi.org/10.1088/1361-6633/aae6b5] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30462612]

40. Phinney, E.S. A Practical Theorem on Gravitational Wave Backgrounds. arXiv; 2001; arXiv: 0108028

41. Isaacson, R.A. Gravitational Radiation in the Limit of High Frequency. I. The Linear Approximation and Geometrical Optics. Phys. Rev.; 1968; 166, pp. 1263-1271. [DOI: https://dx.doi.org/10.1103/PhysRev.166.1263]

42. Isaacson, R.A. Gravitational Radiation in the Limit of High Frequency. II. Nonlinear Terms and the Ef fective Stress Tensor. Phys. Rev.; 1968; 166, pp. 1272-1279. [DOI: https://dx.doi.org/10.1103/PhysRev.166.1272]

43. Renzini, A.I.; Romano, J.D.; Contaldi, C.R.; Cornish, N.J. A comparison of maximum likelihood mapping methods for gravitational-wave backgrounds. Phys. Rev. D; 2022; 105, 023519. [DOI: https://dx.doi.org/10.1103/PhysRevD.105.023519]

44. Grishchuk, L. Amplification of gravitational waves in an isotropic universe. Sov. J. Exp. Theor. Phys.; 1975; 40, 409.

45. Maggiore, M. Gravitational wave experiments and early universe cosmology. Phys. Rep.; 2000; 331, pp. 283-367. [DOI: https://dx.doi.org/10.1016/S0370-1573(99)00102-7]

46. Hogan, C.J. Gravitational radiation from cosmological phase transitions. Mon. Not. R. Astron. Soc.; 1986; 218, pp. 629-636. [DOI: https://dx.doi.org/10.1093/mnras/218.4.629]

47. Battye, R.A.; Caldwell, R.R.; Shellard, E.P.S. Gravitational waves from cosmic strings. Topological Defects in Cosmology; World Scientific Pub. Co. Inc.: Singapore, 1997; pp. 11-31.

48. Vilenkin, A.; Shellard, E.S. Cosmic Strings and Other Topological Defects; Cambridge University Press: Cambridge, UK, 2000.

49. Bartolo, N.; De Luca, V.; Franciolini, G.; Lewis, A.; Peloso, M.; Riotto, A. Primordial Black Hole Dark Matter: LISA Serendipity. Phys. Rev. Lett.; 2019; 122, 211301. [DOI: https://dx.doi.org/10.1103/PhysRevLett.122.211301]

50. Bartolo, N.; Bertacca, D.; De Luca, V.; Franciolini, G.; Matarrese, S.; Peloso, M.; Ricciardone, A.; Riotto, A.; Tasinato, G. Gravitational wave anisotropies from primordial black holes. J. Cosmol. Astropart. Phys.; 2020; 02, 028. [DOI: https://dx.doi.org/10.1088/1475-7516/2020/02/028]

51. Margalit, A.; Contaldi, C.R.; Pieroni, M. Phase Decoherence of Gravitational Wave Backgrounds. Phys. Rev. D; 2020; 102, 083506. [DOI: https://dx.doi.org/10.1103/PhysRevD.102.083506]

52. Maggiore, M. Gravitational Waves. Vol. 2: Astrophysics and Cosmology; Oxford University Press: Oxford, UK, 2018.

53. Auclair, P.; Blanco-Pillado, J.J.; Figueroa, D.G.; Jenkins, A.C.; Lewicki, M.; Sakellariadou, M.; Sanidas, S.; Sousa, L.; Steer, D.A.; Wachter, J.M. et al. Probing the gravitational wave background from cosmic strings with LISA. J. Cosmol. Astropart. Phys.; 2020; 04, 034. [DOI: https://dx.doi.org/10.1088/1475-7516/2020/04/034]

54. Arzoumanian, Z.; Baker, P.T.; Blumer, H.; Bécsy, B.; Brazier, A.; Brook, P.R.; Burke-Spolaor, S.; Chatterjee, S.; Chen, S.; Cordes, J.M. et al. The NANOGrav 12.5 yr Data Set: Search for an Isotropic Stochastic Gravitational-wave Background. ApJ; 2020; 905, L34. [DOI: https://dx.doi.org/10.3847/2041-8213/abd401]

55. Caprini, C.; Hindmarsh, M.; Huber, S.; Konstandin, T.; Kozaczuk, J.; Nardini, G.; No, J.M.; Petiteau, A.; Schwaller, P.; Servant, G. et al. Science with the space-based interferometer eLISA. II: Gravitational waves from cosmological phase transitions. JCAP; 2016; 04, 001. [DOI: https://dx.doi.org/10.1088/1475-7516/2016/04/001]

56. Abbott, R.; Abbott, T.D.; Abraham, S.; Acernese, F.; Ackley, K.; Adams, A.; Adams, C.; Adhikari, R.X.; Adya, V.B.; Affeldt, C. et al. Upper limits on the isotropic gravitational-wave background from Advanced LIGO and Advanced Virgo’s third observing run. Phys. Rev. D; 2021; 104, 022004. [DOI: https://dx.doi.org/10.1103/PhysRevD.104.022004]

57. Kalogera, V.; Baym, G. The maximum mass of a neutron star. Astrophys. J. Lett.; 1996; 470, pp. L61-L64. [DOI: https://dx.doi.org/10.1086/310296]

58. Bailyn, C.D.; Jain, R.K.; Coppi, P.; Orosz, J.A. The Mass Distribution of Stellar Black Holes. ApJ; 1998; 499, pp. 367-374. [DOI: https://dx.doi.org/10.1086/305614]

59. Özel, F.; Psaltis, D.; Narayan, R.; McClintock, J.E. The Black Hole Mass Distribution in the Galaxy. ApJ; 2010; 725, pp. 1918-1927. [DOI: https://dx.doi.org/10.1088/0004-637X/725/2/1918]

60. Gupta, A.; Gerosa, D.; Arun, K.G.; Berti, E.; Farr, W.M.; Sathyaprakash, B.S. Black holes in the low mass gap: Implications for gravitational wave observations. Phys. Rev. D; 2020; 101, 103036. [DOI: https://dx.doi.org/10.1103/PhysRevD.101.103036]

61. Celotti, A.; Miller, J.C.; Sciama, D.W. Astrophysical evidence for the existence of black holes. Class. Quantum Gravity; 2000; 16, A3. [DOI: https://dx.doi.org/10.1088/0264-9381/16/12A/301]

62. Cooray, A.R. Gravitational wave background of neutron star-white dwarf binaries. Mon. Not. R. Astron. Soc.; 2004; 354, pp. 25-30. [DOI: https://dx.doi.org/10.1111/j.1365-2966.2004.08152.x]

63. Korol, V.; Belokurov, V.; Moore, C.J.; Toonen, S. Weighing Milky Way Satellites with LISA. Mon. Not. R. Astron. Soc.; 2021; 502, pp. L55-L60. [DOI: https://dx.doi.org/10.1093/mnrasl/slab003]

64. Regimbau, T. The astrophysical gravitational wave stochastic background. Res. Astron. Astrophys.; 2011; 11, 369. [DOI: https://dx.doi.org/10.1088/1674-4527/11/4/001]

65. Callister, T.; Fishbach, M.; Holz, D.; Farr, W. Shouts and Murmurs: Combining Individual Gravitational-Wave Sources with the Stochastic Background to Measure the History of Binary Black Hole Mergers. Astrophys. J. Lett.; 2020; 896, L32. [DOI: https://dx.doi.org/10.3847/2041-8213/ab9743]

66. Fishbach, M.; Holz, D.E.; Farr, W.M. Does the Black Hole Merger Rate Evolve with Redshift?. Astrophys. J. Lett.; 2018; 863, L41. [DOI: https://dx.doi.org/10.3847/2041-8213/aad800]

67. Madau, P.; Dickinson, M. Cosmic Star Formation History. Annu. Rev. Astron. Astrophys.; 2014; 52, pp. 415-486. [DOI: https://dx.doi.org/10.1146/annurev-astro-081811-125615]

68. Boco, L.; Lapi, A.; Goswami, S.; Perrotta, F.; Baccigalupi, C.; Danese, L. Merging Rates of Compact Binaries in Galaxies: Perspectives for Gravitational Wave Detections. Astrophys. J.; 2019; 881, 157. [DOI: https://dx.doi.org/10.3847/1538-4357/ab328e]

69. Boco, L.; Lapi, A.; Chruslinska, M.; Donevski, D.; Sicilia, A.; Danese, L. Evolution of Galaxy Star Formation and Metallicity: Impact on Double Compact Objects Mergers. Astrophys. J.; 2021; 907, 110. [DOI: https://dx.doi.org/10.3847/1538-4357/abd3a0]

70. Callister, T.; Sammut, L.; Qiu, S.; Mandel, I.; Thrane, E. The limits of astrophysics with gravitational-wave backgrounds. Phys. Rev. X; 2016; 6, 031018. [DOI: https://dx.doi.org/10.1103/PhysRevX.6.031018]

71. Sesana, A.; Vecchio, A.; Colacino, C.N. The stochastic gravitational-wave background from massive black hole binary systems: Implications for observations with Pulsar Timing Arrays. Mon. Not. R. Astron. Soc.; 2008; 390, pp. 192-209. [DOI: https://dx.doi.org/10.1111/j.1365-2966.2008.13682.x]

72. Thorne, K.; Hawking, S.; Israel, W. Three Hundred Years of Gravitation; Cambridge University Press: Cambridge, UK, 1987; 330.

73. Abbott, B.P.; Abbott, R.; Abbott, T.D.; Acernese, F.; Ackley, K.; Adams, C.; Adams, T.; Addesso, P.; Adhikari, R.X.; Adya, V.B. et al. GW170817: Implications for the Stochastic Gravitational-Wave Background from Compact Binary Coalescences. Phys. Rev. Lett.; 2018; 120, 091101. [DOI: https://dx.doi.org/10.1103/PhysRevLett.120.091101] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/29547330]

74. LIGO Scientific CollaborationVirgo CollaborationKAGRA Scientific Collaboration. The population of merging compact binaries inferred using gravitational waves through GWTC-3. arXiv; 2021; arXiv: 2111.03634

75. LIGO Scientific CollaborationVirgo CollaborationKAGRA Collaboration. The Population of Merging Compact Binaries Inferred Using Gravitational Waves through GWTC-3—Data Release. 2021; Available online: https://zenodo.org/record/5655785 (accessed on 1 November 2021).

76. Caprini, C.; Figueroa, D.G. Cosmological backgrounds of gravitational waves. Class. Quantum Gravity; 2018; 35, 163001. [DOI: https://dx.doi.org/10.1088/1361-6382/aac608]

77. Guzzetti, M.C.; Bartolo, N.; Liguori, M.; Matarrese, S. Gravitational waves from inflation. Riv. Nuovo Cim.; 2016; 39, pp. 399-495. [DOI: https://dx.doi.org/10.1393/ncr/i2016-10127-1]

78. Contaldi, C.R.; Magueijo, J. Unsqueezing of standing waves due to inflationary domain structure. Phys. Rev. D; 2018; 98, 043523. [DOI: https://dx.doi.org/10.1103/PhysRevD.98.043523]

79. Kamionkowski, M.; Kosowsky, A.; Stebbins, A. A Probe of primordial gravity waves and vorticity. Phys. Rev. Lett.; 1997; 78, pp. 2058-2061. [DOI: https://dx.doi.org/10.1103/PhysRevLett.78.2058]

80. COrE Collaboration Bouchet, F.R. COrE (Cosmic Origins Explorer) A White Paper. arXiv; 2011; arXiv: 1102.2181

81. Matsumura, T.; Akiba, Y.; Borrill, J.; Chinone, Y.; Dobbs, M.; Fuke, H.; Ghribi, A.; Hasegawa, M.; Hattori, K.; Hattori, M. et al. Mission design of LiteBIRD. J. Low Temp. Phys.; 2014; 176, 733. [DOI: https://dx.doi.org/10.1007/s10909-013-0996-1]

82. Abazajian, K.N.; Adshead, P.; Ahmed, Z.; Allen, S.W.; Alonso, D.; Arnold, K.S.; Baccigalupi, C.; Bartlett, J.G.; Battaglia, N.; Benson, B.A. et al. CMB-S4 Science Book, First Edition. arXiv; 2016; arXiv: 1610.02743

83. Ade, P.; Aguirre, J.; Ahmed, Z.; Aiola, S.; Ali, A.; Alonso, D.; Alvarez, M.A.; Arnold, K.; Ashton, P.; Austermann, J. et al. The Simons Observatory: Science goals and forecasts. JCAP; 2019; 02, 056. [DOI: https://dx.doi.org/10.1088/1475-7516/2019/02/056]

84. Kuroyanagi, S.; Chiba, T.; Sugiyama, N. Precision calculations of the gravitational wave background spectrum from inflation. Phys. Rev. D; 2009; 79, 103501. [DOI: https://dx.doi.org/10.1103/PhysRevD.79.103501]

85. Chung, D.J.H.; Kolb, E.W.; Riotto, A.; Tkachev, I.I. Probing Planckian physics: Resonant production of particles during inflation and features in the primordial power spectrum. Phys. Rev. D; 2000; 62, 043508. [DOI: https://dx.doi.org/10.1103/PhysRevD.62.043508]

86. Ananda, K.N.; Clarkson, C.; Wands, D. The Cosmological gravitational wave background from primordial density perturbations. Phys. Rev. D; 2007; 75, 123518. [DOI: https://dx.doi.org/10.1103/PhysRevD.75.123518]

87. Clesse, S.; García-Bellido, J.; Orani, S. Detecting the Stochastic Gravitational Wave Background from Primordial Black Hole Formation. arXiv; 2018; arXiv: 1812.11011

88. Garcia-Bellido, J.; Peloso, M.; Unal, C. Gravitational Wave signatures of inflationary models from Primordial Black Hole Dark Matter. JCAP; 2017; 09, 013. [DOI: https://dx.doi.org/10.1088/1475-7516/2017/09/013]

89. Khlebnikov, S.Y.; Tkachev, I.I. Relic gravitational waves produced after preheating. Phys. Rev. D; 1997; 56, pp. 653-660. [DOI: https://dx.doi.org/10.1103/PhysRevD.56.653]

90. Easther, R.; Lim, E.A. Stochastic gravitational wave production after inflation. JCAP; 2006; 04, 010. [DOI: https://dx.doi.org/10.1088/1475-7516/2006/04/010]

91. Garcia-Bellido, J.; Figueroa, D.G. A stochastic background of gravitational waves from hybrid preheating. Phys. Rev. Lett.; 2007; 98, 061302. [DOI: https://dx.doi.org/10.1103/PhysRevLett.98.061302]

92. Easther, R.; Giblin, J.T.; Lim, E.A. Gravitational Wave Production at the End of Inflation. Phys. Rev. Lett.; 2007; 99, 221301. [DOI: https://dx.doi.org/10.1103/PhysRevLett.99.221301]

93. Bartolo, N.; Caprini, C.; Domcke, V.; Figueroa, D.G.; Garcia-Bellido, J.; Guzzetti, M.C.; Liguori, M.; Matarrese, S.; Peloso, M.; Petiteau, A. et al. Science with the space-based interferometer LISA. IV: Probing inflation with gravitational waves. J. Cosmol. Astropart. Phys.; 2016; 12, 026. [DOI: https://dx.doi.org/10.1088/1475-7516/2016/12/026]

94. Jeannerot, R.; Rocher, J.; Sakellariadou, M. How generic is cosmic string formation in SUSY GUTs. Phys. Rev. D; 2003; 68, 103514. [DOI: https://dx.doi.org/10.1103/PhysRevD.68.103514]

95. García-Bellido, J.; Jaraba, S.; Kuroyanagi, S. The stochastic gravitational wave background from close hyperbolic encounters of primordial black holes in dense clusters. arXiv; 2021; arXiv: 2109.11376

96. Alexander, S.; Marciano, A.; Smolin, L. Gravitational origin of the weak interaction’s chirality. Phys. Rev. D; 2014; 89, 065017. [DOI: https://dx.doi.org/10.1103/PhysRevD.89.065017]

97. Freidel, L.; Minic, D.; Takeuchi, T. Quantum gravity, torsion, parity violation and all that. Phys. Rev. D; 2005; 72, 104002. [DOI: https://dx.doi.org/10.1103/PhysRevD.72.104002]

98. Alexander, S.H. Isogravity: Toward an Electroweak and Gravitational Unification. arXiv; 2007; arXiv: 0706.4481

99. Smolin, L. The Plebanski action extended to a unification of gravity and Yang-Mills theory. Phys. Rev. D; 2009; 80, 124017. [DOI: https://dx.doi.org/10.1103/PhysRevD.80.124017]

100. Contaldi, C.R. Anisotropies of Gravitational Wave Backgrounds: A Line Of Sight Approach. Phys. Rev. Lett.; 2017; B771, pp. 9-12. [DOI: https://dx.doi.org/10.1016/j.physletb.2017.05.020]

101. Jenkins, A.C.; Sakellariadou, M. Anisotropies in the stochastic gravitational-wave background: Formalism and the cosmic string case. Phys. Rev. D; 2018; 98, 063509. [DOI: https://dx.doi.org/10.1103/PhysRevD.98.063509]

102. Jenkins, A.C.; O’Shaughnessy, R.; Sakellariadou, M.; Wysocki, D. Anisotropies in the astrophysical gravitational-wave background: The impact of black hole distributions. Phys. Rev. Lett.; 2019; 122, 111101. [DOI: https://dx.doi.org/10.1103/PhysRevLett.122.111101]

103. Durrer, R. The Cosmic Microwave Background; Cambridge University Press: Cambridge, UK, 2008; [DOI: https://dx.doi.org/10.1017/CBO9780511817205]

104. Peebles, P.J.E. The Large-Scale Structure of the Universe; Princeton University Press: Princeton, NJ, USA, 2020.

105. Sachs, R.; Wolfe, A. Perturbations of a cosmological model and angular variations of the microwave background. Astrophys. J.; 1967; 147, pp. 73-90. [DOI: https://dx.doi.org/10.1086/148982]

106. Bartolo, N.; Bertacca, D.; Matarrese, S.; Peloso, M.; Ricciardone, A.; Riotto, A.; Tasinato, G. Anisotropies and non-Gaussianity of the Cosmological Gravitational Wave Background. Phys. Rev. D; 2019; 100, 121501. [DOI: https://dx.doi.org/10.1103/PhysRevD.100.121501]

107. Cusin, G.; Dvorkin, I.; Pitrou, C.; Uzan, J.P. First predictions of the angular power spectrum of the astrophysical gravitational wave background. Phys. Rev. Lett.; 2018; 120, 231101. [DOI: https://dx.doi.org/10.1103/PhysRevLett.120.231101] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/29932710]

108. Jenkins, A.C.; Sakellariadou, M.; Regimbau, T.; Slezak, E. Anisotropies in the astrophysical gravitational-wave background: Predictions for the detection of compact binaries by LIGO and Virgo. Phys. Rev. D; 2018; 98, 063501. [DOI: https://dx.doi.org/10.1103/PhysRevD.98.063501]

109. Capurri, G.; Lapi, A.; Baccigalupi, C.; Boco, L.; Scelfo, G.; Ronconi, T. Intensity and anisotropies of the stochastic gravitational wave background from merging compact binaries in galaxies. JCAP; 2021; 11, 032. [DOI: https://dx.doi.org/10.1088/1475-7516/2021/11/032]

110. Geller, M.; Hook, A.; Sundrum, R.; Tsai, Y. Primordial Anisotropies in the Gravitational Wave Background from Cosmological Phase Transitions. Phys. Rev. Lett.; 2018; 121, 201303. [DOI: https://dx.doi.org/10.1103/PhysRevLett.121.201303]

111. Bertacca, D.; Ricciardone, A.; Bellomo, N.; Jenkins, A.C.; Matarrese, S.; Raccanelli, A.; Regimbau, T.; Sakellariadou, M. Projection effects on the observed angular spectrum of the astrophysical stochastic gravitational wave background. Phys. Rev. D; 2020; 101, 103513. [DOI: https://dx.doi.org/10.1103/PhysRevD.101.103513]

112. Bellomo, N.; Bertacca, D.; Jenkins, A.C.; Matarrese, S.; Raccanelli, A.; Regimbau, T.; Ricciardone, A.; Sakellariadou, M. CLASS_GWB: Robust modeling of the astrophysical gravitational wave background anisotropies. arXiv; 2021; arXiv: 2110.15059

113. Cusin, G.; Dvorkin, I.; Pitrou, C.; Uzan, J.P. Comment on the article “Anisotropies in the astrophysical gravitational-wave background: The impact of black hole distributions” by A.C. Jenkins et al. [arXiv:1810.13435]. arXiv; 2018; arXiv: 1811.03582

114. Jenkins, A.C.; Sakellariadou, M.; Regimbau, T.; Slezak, E.; O’Shaughnessy, R.; Wysocki, D. Response to Cusin et al’s comment on arXiv:1810.13435. arXiv; 2019; arXiv: 1901.01078

115. Pitrou, C.; Cusin, G.; Uzan, J.P. Unified view of anisotropies in the astrophysical gravitational-wave background. Phys. Rev. D; 2020; 101, 081301. [DOI: https://dx.doi.org/10.1103/PhysRevD.101.081301]

116. Alonso, D.; Contaldi, C.R.; Cusin, G.; Ferreira, P.G.; Renzini, A.I. The N_ℓ of gravitational wave background experiments. arXiv; 2020; arXiv: 2005.03001

117. Regimbau, T.; Dent, T.; Del Pozzo, W.; Giampanis, S.; Li, T.G.F.; Robinson, C.; Van Den Broeck, C.; Meacher, D.; Rodriguez, C.; Sathyaprakash, B.S. et al. Mock data challenge for the Einstein Gravitational-Wave Telescope. Phys. Rev. D; 2012; 86, 122001. [DOI: https://dx.doi.org/10.1103/PhysRevD.86.122001]

118. Regimbau, T.; Meacher, D.; Coughlin, M. Second Einstein Telescope mock science challenge: Detection of the gravitational-wave stochastic background from compact binary coalescences. Phys. Rev. D; 2014; 89, 084046. [DOI: https://dx.doi.org/10.1103/PhysRevD.89.084046]

119. Cornish, N.J.; Romano, J.D. When is a gravitational-wave signal stochastic?. Phys. Rev. D; 2015; 92, 042001. [DOI: https://dx.doi.org/10.1103/PhysRevD.92.042001]

120. Rosado, P.A.; Sesana, A.; Gair, J. Expected properties of the first gravitational wave signal detected with pulsar timing arrays. Mon. Not. R. Astron. Soc.; 2015; 451, pp. 2417-2433. [DOI: https://dx.doi.org/10.1093/mnras/stv1098]

121. Meacher, D.; Thrane, E.; Regimbau, T. Statistical properties of astrophysical gravitational-wave backgrounds. Phys. Rev. D; 2014; 89, 084063. [DOI: https://dx.doi.org/10.1103/PhysRevD.89.084063]

122. Jenkins, A.C.; Sakellariadou, M. Shot noise in the astrophysical gravitational-wave background. Phys. Rev. D; 2019; 100, 063508. [DOI: https://dx.doi.org/10.1103/PhysRevD.100.063508]

123. Jenkins, A.C.; Romano, J.D.; Sakellariadou, M. Estimating the angular power spectrum of the gravitational-wave background in the presence of shot noise. Phys. Rev. D; 2019; 100, 083501. [DOI: https://dx.doi.org/10.1103/PhysRevD.100.083501]

124. Alonso, D.; Cusin, G.; Ferreira, P.G.; Pitrou, C. Detecting the anisotropic astrophysical gravitational wave background in the presence of shot noise through cross-correlations. Phys. Rev. D; 2020; 102, 023002. [DOI: https://dx.doi.org/10.1103/PhysRevD.102.023002]

125. Canas-Herrera, G.; Contigiani, O.; Vardanyan, V. Cross-correlation of the astrophysical gravitational-wave background with galaxy clustering. Phys. Rev. D; 2020; 102, 043513. [DOI: https://dx.doi.org/10.1103/PhysRevD.102.043513]

126. Dvorkin, I.; Uzan, J.P.; Vangioni, E.; Silk, J. Exploring stellar evolution with gravitational-wave observations. Mon. Not. R. Astron. Soc.; 2018; 479, pp. 121-129. [DOI: https://dx.doi.org/10.1093/mnras/sty1414]

127. Mukherjee, S.; Silk, J. Time-dependence of the astrophysical stochastic gravitational wave background. Mon. Not. R. Astron. Soc.; 2020; 491, pp. 4690-4701. [DOI: https://dx.doi.org/10.1093/mnras/stz3226]

128. Abbott, R.; Abbott, T.D.; Acernese, F.; Ackley, K.; Adams, C.; Adhikari, N.; Adhikari, R.X.; Adya, V.B.; Affeldt, C.; Agarwal, D. et al. Constraints on dark photon dark matter using data from LIGO’s and Virgo’s third observing run. arXiv; 2021; arXiv: 2105.13085

129. Abbott, R.; Abbott, T.D.; Abraham, S.; Acernese, F.; Ackley, K.; Adams, A.; Adams, C.; Adhikari, R.X.; Adya, V.B.; Affeldt, C. et al. Constraints on Cosmic Strings Using Data from the Third Advanced LIGO–Virgo Observing Run. Phys. Rev. Lett.; 2021; 126, 241102. [DOI: https://dx.doi.org/10.1103/PhysRevLett.126.241102] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/34213926]

130. Talbot, C.; Thrane, E.; Biscoveanu, S.; Smith, R. Inference with finite time series: Observing the gravitational Universe through windows. Phys. Rev. Res.; 2021; 3, 043049. [DOI: https://dx.doi.org/10.1103/PhysRevResearch.3.043049]

131. Matas, A.; Romano, J.D. Frequentist versus Bayesian analyses: Cross-correlation as an approximate sufficient statistic for LIGO-Virgo stochastic background searches. Phys. Rev. D; 2021; 103, 062003. [DOI: https://dx.doi.org/10.1103/PhysRevD.103.062003]

132. Littenberg, T.; Cornish, N.; Lackeos, K.; Robson, T. Global Analysis of the Gravitational Wave Signal from Galactic Binaries. Phys. Rev. D; 2020; 101, 123021. [DOI: https://dx.doi.org/10.1103/PhysRevD.101.123021]

133. Aasi, J.; Abbott, B.P.; Abbott, R.; Abbott, T.; Abernathy, M.R.; Ackley, K.; Adams, C.; Adams, T.; Addesso, P.; Adhikari, R.X. et al. Advanced LIGO. Class. Quantum Gravity; 2015; 32, 074001. [DOI: https://dx.doi.org/10.1088/0264-9381/32/7/074001]

134. Armstrong, J.W.; Estabrook, F.B.; Tinto, M. Time-Delay Interferometry for Space-based Gravitational Wave Searches. Astrophys. J.; 1999; 527, pp. 814-826. [DOI: https://dx.doi.org/10.1086/308110]

135. Tinto, M.; Armstrong, J.W. Cancellation of laser noise in an unequal-arm interferometer detector of gravitational radiation. Phys. Rev. D; 1999; 59, 102003. [DOI: https://dx.doi.org/10.1103/PhysRevD.59.102003]

136. Seto, N.; Taruya, A. Measuring a parity-violation signature in the early universe via ground-based laser interferometers. Phys. Rev. Lett.; 2007; 99, 121101. [DOI: https://dx.doi.org/10.1103/PhysRevLett.99.121101] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/17930491]

137. Seto, N.; Taruya, A. Polarization analysis of gravitational-wave backgrounds from the correlation signals of ground-based interferometers: Measuring a circular-polarization mode. Phys. Rev. D Part. Fields, Gravit. Cosmol.; 2008; 77, pp. 1-27. [DOI: https://dx.doi.org/10.1103/PhysRevD.77.103001]

138. Bayle, J.B.; Hartwig, O.; Staab, M. Adapting time-delay interferometry for LISA data in frequency. Phys. Rev. D; 2021; 104, 023006. [DOI: https://dx.doi.org/10.1103/PhysRevD.104.023006]

139. Bayle, J.B.; Lilley, M.; Petiteau, A.; Halloin, H. Effect of filters on the time-delay interferometry residual laser noise for LISA. PRD; 2019; 99, 084023. [DOI: https://dx.doi.org/10.1103/PhysRevD.99.084023]

140. Vallisneri, M.; Bayle, J.B.; Babak, S.; Petiteau, A. Time-delay interferometry without delays. Phys. Rev. D; 2021; 103, 082001. [DOI: https://dx.doi.org/10.1103/PhysRevD.103.082001]

141. Bayle, J.B.; Vallisneri, M.; Babak, S.; Petiteau, A. On the matrix formulation of time-delay interferometry. arXiv; 2021; arXiv: 2106.03976

142. Tinto, M.; Dhurandhar, S.; Joshi, P. Matrix representation of time-delay interferometry. PRD; 2021; 104, 044033. [DOI: https://dx.doi.org/10.1103/PhysRevD.104.044033]

143. Page, J.; Littenberg, T.B. Bayesian time delay interferometry. PRD; 2021; 104, 084037. [DOI: https://dx.doi.org/10.1103/PhysRevD.104.084037]

144. Baghi, Q.; Thorpe, J.I.; Slutsky, J.; Baker, J. Statistical inference approach to time-delay interferometry for gravitational-wave detection. PRD; 2021; 103, 042006. [DOI: https://dx.doi.org/10.1103/PhysRevD.103.042006]

145. LISA Data Challenge Manual. 2018; Available online: https://lisa-ldc.lal.in2p3.fr/static/data/pdf/LDC-manual-001.pdf (accessed on 1 November 2021).

146. Tinto, M.; Dhurandhar, S.V. TIME DELAY. Living Rev. Relativ.; 2005; 8, 4. [DOI: https://dx.doi.org/10.12942/lrr-2005-4] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/28163648]

147. Adams, M.R.; Cornish, N.J. Discriminating between a stochastic gravitational wave background and instrument noise. Phys. Rev. D; 2010; 82, 022002. [DOI: https://dx.doi.org/10.1103/PhysRevD.82.022002]

148. Adams, M.R.; Cornish, N.J. Detecting a stochastic gravitational wave background in the presence of a galactic foreground and instrument noise. Phys. Rev. D; 2014; 89, 022001. [DOI: https://dx.doi.org/10.1103/PhysRevD.89.022001]

149. Anholm, M.; Ballmer, S.; Creighton, J.D.E.; Price, L.R.; Siemens, X. Optimal strategies for gravitational wave stochastic background searches in pulsar timing data. PRD; 2009; 79, 084030. [DOI: https://dx.doi.org/10.1103/PhysRevD.79.084030]

150. Hellings, R.W.; Downs, G.S. Upper limits on the isotropic gravitational radiation background from pulsar timing analysis. ApJ; 1983; 265, pp. L39-L42. [DOI: https://dx.doi.org/10.1086/183954]

151. Allen, B.; Romano, J.D. Detecting a stochastic background of gravitational radiation: Signal processing strategies and sensitivities. Phys. Rev. D; 1999; 59, 102001. [DOI: https://dx.doi.org/10.1103/PhysRevD.59.102001]

152. Renzini, A.I.; Contaldi, C.R. Mapping Incoherent Gravitational Wave Backgrounds. Mon. Not. R. Astron. Soc.; 2018; 481, pp. 4650-4661. [DOI: https://dx.doi.org/10.1093/mnras/sty2546]

153. Renzini, A.I.; Contaldi, C.R. Gravitational Wave Background Sky Maps from Advanced LIGO O1 Data. Phys. Rev. Lett.; 2019; 122, 081102. [DOI: https://dx.doi.org/10.1103/PhysRevLett.122.081102]

154. Renzini, A.; Contaldi, C. Improved limits on a stochastic gravitational-wave background and its anisotropies from Advanced LIGO O1 and O2 runs. Phys. Rev.; 2019; D100, 063527. [DOI: https://dx.doi.org/10.1103/PhysRevD.100.063527]

155. Chatziioannou, K.; Haster, C.J.; Littenberg, T.B.; Farr, W.M.; Ghonge, S.; Millhouse, M.; Clark, J.A.; Cornish, N. Noise spectral estimation methods and their impact on gravitational wave measurement of compact binary mergers. Phys. Rev. D; 2019; 100, 104004. [DOI: https://dx.doi.org/10.1103/PhysRevD.100.104004]

156. Cornish, N.J.; Romano, J.D. Towards a unified treatment of gravitational-wave data analysis. Phys. Rev. D; 2013; 87, 122003. [DOI: https://dx.doi.org/10.1103/PhysRevD.87.122003]

157. Cornish, N.J. Detecting a stochastic gravitational wave background with the Laser Interferometer Space Antenna. Phys. Rev. D; 2002; 65, 022004. [DOI: https://dx.doi.org/10.1103/PhysRevD.65.022004]

158. Coughlin, M.; Christensen, N.; Gair, J.; Kandhasamy, S.; Thrane, E. Method for estimation of gravitational-wave transient model parameters in frequency–time maps. Class. Quantum Gravity; 2014; 31, 165012. [DOI: https://dx.doi.org/10.1088/0264-9381/31/16/165012]

159. Mandic, V.; Thrane, E.; Giampanis, S.; Regimbau, T. Parameter Estimation in Searches for the Stochastic Gravitational-Wave Background. Phys. Rev. Lett.; 2012; 109, 171102. [DOI: https://dx.doi.org/10.1103/PhysRevLett.109.171102]

160. Callister, T.; Biscoveanu, A.S.; Christensen, N.; Isi, M.; Matas, A.; Minazzoli, O.; Regimbau, T.; Sakellariadou, M.; Tasson, J.; Thrane, E. Polarization-based Tests of Gravity with the Stochastic Gravitational-Wave Background. Phys. Rev. X; 2017; 7, 041058. [DOI: https://dx.doi.org/10.1103/PhysRevX.7.041058]

161. Abbott, B.P.; Abbott, R.; Abbott, T.D.; Abraham, S.; Acernese, F.; Ackley, K.; Adams, C.; Adya, V.B.; Affeldt, C.; Agathos, M. et al. Search for the isotropic stochastic background using data from Advanced LIGO’s second observing run. Phys. Rev. D; 2019; 100, 061101. [DOI: https://dx.doi.org/10.1103/PhysRevD.100.061101]

162. Tsukada, L.; Callister, T.; Matas, A.; Meyers, P. First search for a stochastic gravitational-wave background from ultralight bosons. Phys. Rev. D; 2019; 99, 103015. [DOI: https://dx.doi.org/10.1103/PhysRevD.99.103015]

163. Tsukada, L.; Brito, R.; East, W.E.; Siemonsen, N. Modeling and searching for a stochastic gravitational-wave background from ultralight vector bosons. Phys. Rev. D; 2021; 103, 083005. [DOI: https://dx.doi.org/10.1103/PhysRevD.103.083005]

164. Romero, A.; Martinovic, K.; Callister, T.A.; Guo, H.K.; Martínez, M.; Sakellariadou, M.; Yang, F.W.; Zhao, Y. Implications for First-Order Cosmological Phase Transitions from the Third LIGO-Virgo Observing Run. Phys. Rev. Lett.; 2021; 126, 151301. [DOI: https://dx.doi.org/10.1103/PhysRevLett.126.151301]

165. Romero-Rodriguez, A.; Martinez, M.; Pujolàs, O.; Sakellariadou, M.; Vaskonen, V. Search for a scalar induced stochastic gravitational wave background in the third LIGO-Virgo observing run. Phys. Rev. Lett.; 2022; 128, 051301. [DOI: https://dx.doi.org/10.1103/PhysRevLett.128.051301]

166. Martinovic, K.; Badger, C.; Sakellariadou, M.; Mandic, V. Searching for parity violation with the LIGO-Virgo-KAGRA network. Phys. Rev. D; 2021; 104, L081101. [DOI: https://dx.doi.org/10.1103/PhysRevD.104.L081101]

167. Meyers, P.M.; Martinovic, K.; Christensen, N.; Sakellariadou, M. Detecting a stochastic gravitational-wave background in the presence of correlated magnetic noise. PRD; 2020; 102, 102005. [DOI: https://dx.doi.org/10.1103/PhysRevD.102.102005]

168. Martinovic, K.; Meyers, P.M.; Sakellariadou, M.; Christensen, N. Simultaneous estimation of astrophysical and cosmological stochastic gravitational-wave backgrounds with terrestrial detectors. Phys. Rev. D; 2021; 103, 043023. [DOI: https://dx.doi.org/10.1103/PhysRevD.103.043023]

169. Parida, A.; Suresh, J.; Mitra, S.; Jhingan, S. Component separation map-making for stochastic gravitational wave background. arXiv; 2019; arXiv: 1904.05056

170. Meacher, D.; Coughlin, M.; Morris, S.; Regimbau, T.; Christensen, N.; Kandhasamy, S.; Mandic, V.; Romano, J.D.; Thrane, E. Mock data and science challenge for detecting an astrophysical stochastic gravitational-wave background with Advanced LIGO and Advanced Virgo. Phys. Rev. D; 2015; 92, 063002. [DOI: https://dx.doi.org/10.1103/PhysRevD.92.063002]

171. Biwer, C.; Barker, D.; Batch, J.C.; Betzwieser, J.; Fisher, R.P.; Goetz, E.; Kandhasamy, S.; Karki, S.; Kissel, J.S.; Lundgren, A.P. et al. Validating gravitational-wave detections: The Advanced LIGO hardware injection system. Phys. Rev. D; 2017; 95, 062002. [DOI: https://dx.doi.org/10.1103/PhysRevD.95.062002]

172. van Haasteren, R.; Levin, Y.; McDonald, P.; Lu, T. On measuring the gravitational-wave background using Pulsar Timing Arrays. MNRAS; 2009; 395, pp. 1005-1014. [DOI: https://dx.doi.org/10.1111/j.1365-2966.2009.14590.x]

173. van Haasteren, R.; Levin, Y.; Janssen, G.H.; Lazaridis, K.; Kramer, M.; Stappers, B.W.; Desvignes, G.; Purver, M.B.; Lyne, A.G.; Ferdman, R.D. et al. Placing limits on the stochastic gravitational-wave background using European Pulsar Timing Array data. MNRAS; 2011; 414, pp. 3117-3128. [DOI: https://dx.doi.org/10.1111/j.1365-2966.2011.18613.x]

174. Lentati, L.; Alexander, P.; Hobson, M.P.; Taylor, S.; Gair, J.; Balan, S.T.; van Haasteren, R. Hyper-efficient model-independent Bayesian method for the analysis of pulsar timing data. PRD; 2013; 87, 104021. [DOI: https://dx.doi.org/10.1103/PhysRevD.87.104021]

175. van Haasteren, R.; Vallisneri, M. New advances in the Gaussian-process approach to pulsar-timing data analysis. PRD; 2014; 90, 104012. [DOI: https://dx.doi.org/10.1103/PhysRevD.90.104012]

176. Lentati, L.; Alexander, P.; Hobson, M.P.; Feroz, F.; van Haasteren, R.; Lee, K.J.; Shannon, R.M. TEMPONEST: A Bayesian approach to pulsar timing analysis. MNRAS; 2014; 437, pp. 3004-3023. [DOI: https://dx.doi.org/10.1093/mnras/stt2122]

177. Hobbs, G.; Jenet, F.; Lee, K.J.; Verbiest, J.P.W.; Yardley, D.; Manchester, R.; Lommen, A.; Coles, W.; Edwards, R.; Shettigara, C. TEMPO2: A new pulsar timing package—III. Gravitational wave simulation. MNRAS; 2009; 394, pp. 1945-1955. [DOI: https://dx.doi.org/10.1111/j.1365-2966.2009.14391.x]

178. Arzoumanian, Z.; Brazier, A.; Burke-Spolaor, S.; Chamberlin, S.J.; Chatterjee, S.; Christy, B.; Cordes, J.M.; Cornish, N.J.; Crowter, K.; Demorest, P.B. et al. The NANOGrav Nine-year Data Set: Limits on the Isotropic Stochastic Gravitational Wave Background. ApJ; 2016; 821, 13. [DOI: https://dx.doi.org/10.3847/0004-637X/821/1/13]

179. Taylor, S.R.; Lentati, L.; Babak, S.; Brem, P.; Gair, J.R.; Sesana, A.; Vecchio, A. All correlations must die: Assessing the significance of a stochastic gravitational-wave background in pulsar timing arrays. PRD; 2017; 95, 042002. [DOI: https://dx.doi.org/10.1103/PhysRevD.95.042002]

180. Arzoumanian, Z.; Baker, P.T.; Brazier, A.; Burke-Spolaor, S.; Chamberlin, S.J.; Chatterjee, S.; Christy, B.; Cordes, J.M.; Cornish, N.J.; Crawford, F. et al. The NANOGrav 11 Year Data Set: Pulsar-timing Constraints on the Stochastic Gravitational-wave Background. ApJ; 2018; 859, 47. [DOI: https://dx.doi.org/10.3847/1538-4357/aabd3b]

181. Chamberlin, S.J.; Creighton, J.D.E.; Siemens, X.; Demorest, P.; Ellis, J.; Price, L.R.; Romano, J.D. Time-domain implementation of the optimal cross-correlation statistic for stochastic gravitational-wave background searches in pulsar timing data. PRD; 2015; 91, 044048. [DOI: https://dx.doi.org/10.1103/PhysRevD.91.044048]

182. Vigeland, S.J.; Islo, K.; Taylor, S.R.; Ellis, J.A. Noise-marginalise optimal statistic: A robust hybrid frequentist-Bayesian statistic for the stochastic gravitational-wave background in pulsar timing arrays. PRD; 2018; 98, 044003. [DOI: https://dx.doi.org/10.1103/PhysRevD.98.044003]

183. Cornish, N.J.; Sampson, L. Towards robust gravitational wave detection with pulsar timing arrays. PRD; 2016; 93, 104047. [DOI: https://dx.doi.org/10.1103/PhysRevD.93.104047]

184. Bécsy, B.; Cornish, N.J. Joint search for isolated sources and an unresolved confusion background in pulsar timing array data. Class. Quantum Gravity; 2020; 37, 135011. [DOI: https://dx.doi.org/10.1088/1361-6382/ab8bbd]

185. Prix, R.; Itoh, Y. Global parameter-space correlations of coherent searches for continuous gravitational waves. Class. Quantum Gravity; 2005; 22, pp. S1003-S1012. [DOI: https://dx.doi.org/10.1088/0264-9381/22/18/S14]

186. Ain, A.; Dalvi, P.; Mitra, S. Fast gravitational wave radiometry using data folding. Phys. Rev. D Part. Fields, Gravit. Cosmol.; 2015; 92, 022003. [DOI: https://dx.doi.org/10.1103/PhysRevD.92.022003]

187. Ain, A.; Suresh, J.; Mitra, S. Very fast stochastic gravitational wave background map making using folded data. PRD; 2018; 98, 024001. [DOI: https://dx.doi.org/10.1103/PhysRevD.98.024001]

188. Goncharov, B.; Thrane, E. All-sky radiometer for narrowband gravitational waves using folded data. PRD; 2018; 98, 064018. [DOI: https://dx.doi.org/10.1103/PhysRevD.98.064018]

189. Mitra, S.; Dhurandhar, S.; Souradeep, T.; Lazzarini, A.; Mandic, V.; Bose, S.; Ballmer, S. Gravitational wave radiometry: Mapping a stochastic gravitational wave background. Phys. Rev. D; 2008; 77, 042002. [DOI: https://dx.doi.org/10.1103/PhysRevD.77.042002]

190. Thrane, E.; Ballmer, S.; Romano, J.D.; Mitra, S.; Talukder, D.; Bose, S.; Mandic, V. Probing the anisotropies of a stochastic gravitational-wave background using a network of ground-based laser interferometers. Phys. Rev. D Part. Fields, Gravit. Cosmol.; 2009; 80, 122002. [DOI: https://dx.doi.org/10.1103/PhysRevD.80.122002]

191. Mingarelli, C.M.F.; Sidery, T.; Mandel, I.; Vecchio, A. Characterizing gravitational wave stochastic background anisotropy with pulsar timing arrays. PRD; 2013; 88, 062005. [DOI: https://dx.doi.org/10.1103/PhysRevD.88.062005]

192. Taylor, S.R.; Gair, J.R. Searching for anisotropic gravitational-wave backgrounds using pulsar timing arrays. PRD; 2013; 88, 084001. [DOI: https://dx.doi.org/10.1103/PhysRevD.88.084001]

193. Contaldi, C.R.; Pieroni, M.; Renzini, A.I.; Cusin, G.; Karnesis, N.; Peloso, M.; Ricciardone, A.; Tasinato, G. Maximum likelihood map-making with the Laser Interferometer Space Antenna. Phys. Rev. D; 2020; 102, 043502. [DOI: https://dx.doi.org/10.1103/PhysRevD.102.043502]

194. Banagiri, S.; Criswell, A.; Kuan, T.; Mandic, V.; Romano, J.D.; Taylor, S.R. Mapping the gravitational-wave sky with LISA: A Bayesian spherical harmonic approach. MNRAS; 2021; 507, pp. 5451-5462. [DOI: https://dx.doi.org/10.1093/mnras/stab2479]

195. Bond, J.R.; Jaffe, A.H.; Knox, L. Estimating the power spectrum of the cosmic microwave background. Phys. Rev.; 1998; D57, pp. 2117-2137. [DOI: https://dx.doi.org/10.1103/PhysRevD.57.2117]

196. Rocha, G.; Contaldi, C.; Bond, J.; Gorski, K. Application of XFaster power spectrum and likelihood estimator to Planck. Mon. Not. R. Astron. Soc.; 2011; 414, pp. 823-846. [DOI: https://dx.doi.org/10.1111/j.1365-2966.2010.17980.x]

197. Cornish, N.J.; van Haasteren, R. Mapping the nano-Hertz gravitational wave sky. arXiv; 2014; arXiv: 1406.4511

198. Ballmer, S.W. A radiometer for stochastic gravitational waves. Class. Quantum Gravity; 2006; 23, pp. S179-S185. [DOI: https://dx.doi.org/10.1088/0264-9381/23/8/S23]

199. Gair, J.; Romano, J.D.; Taylor, S.; Mingarelli, C.M.F. Mapping gravitational-wave backgrounds using methods from CMB analysis: Application to pulsar timing arrays. PRD; 2014; 90, 082001. [DOI: https://dx.doi.org/10.1103/PhysRevD.90.082001]

200. Gair, J.R.; Romano, J.D.; Taylor, S.R. Mapping gravitational-wave backgrounds of arbitrary polarisation using pulsar timing arrays. PRD; 2015; 92, 102003. [DOI: https://dx.doi.org/10.1103/PhysRevD.92.102003]

201. Romano, J.D.; Taylor, S.R.; Cornish, N.J.; Gair, J.; Mingarelli, C.M.F.; van Haasteren, R. Phase-coherent mapping of gravitational-wave backgrounds using ground-based laser interferometers. PRD; 2015; 92, 042003. [DOI: https://dx.doi.org/10.1103/PhysRevD.92.042003]

202. Taylor, S.R.; Mingarelli, C.M.F.; Gair, J.R.; Sesana, A.; Theureau, G.; Babak, S.; Bassa, C.G.; Brem, P.; Burgay, M.; Caballero, R.N. et al. Limits on Anisotropy in the Nanohertz Stochastic Gravitational Wave Background. PRL; 2015; 115, 041101. [DOI: https://dx.doi.org/10.1103/PhysRevLett.115.041101]

203. Ali-Haïmoud, Y.; Smith, T.L.; Mingarelli, C.M.F. Fisher formalism for anisotropic gravitational-wave background searches with pulsar timing arrays. PRD; 2020; 102, 122005. [DOI: https://dx.doi.org/10.1103/PhysRevD.102.122005]

204. Ali-Haïmoud, Y.; Smith, T.L.; Mingarelli, C.M.F. Insights into searches for anisotropies in the nanohertz gravitational-wave background. PRD; 2021; 103, 042009. [DOI: https://dx.doi.org/10.1103/PhysRevD.103.042009]

205. Suresh, J.; Agarwal, D.; Mitra, S. Jointly setting upper limits on multiple components of an anisotropic stochastic gravitational-wave background. PRD; 2021; 104, 102003. [DOI: https://dx.doi.org/10.1103/PhysRevD.104.102003]

206. Gorski, K.M.; Hivon, E.; Banday, A.J.; Wandelt, B.D.; Hansen, F.K.; Reinecke, M.; Bartelman, M. HEALPix—A Framework for high resolution discretization, and fast analysis of data distributed on the sphere. Astrophys. J.; 2005; 622, pp. 759-771. [DOI: https://dx.doi.org/10.1086/427976]

207. Hotinli, S.C.; Kamionkowski, M.; Jaffe, A.H. The search for statistical anisotropy in the gravitational-wave background with pulsar timing arrays. Open J. Astrophys.; 2019; 2, 8. [DOI: https://dx.doi.org/10.21105/astro.1904.05348]

208. Ivezić, Ž.; Connolly, A.; VanderPlas, J.; Gray, A. Statistics, Data Mining, and Machine Learning in Astronomy: A Practical Python Guide for the Analysis of Survey Data; Princeton Series in Modern Observational Astronomy; Princeton University Press: Princeton, NJ, USA, 2014.

209. Suresh, J.; Ain, A.; Mitra, S. Unified mapmaking for an anisotropic stochastic gravitational wave background. PRD; 2021; 103, 083024. [DOI: https://dx.doi.org/10.1103/PhysRevD.103.083024]

210. Myers, S.T.; Contaldi, C.R.; Bond, J.R.; Pen, U.L.; Pogosyan, D.; Prunet, S.; Sievers, J.L.; Mason, B.S.; Pearson, T.J.; Readhead, A.C.S. et al. A fast gridded method for the estimation of the power spectrum of the CMB from interferometer data with application to the cosmic background imager. Astrophys. J.; 2003; 591, pp. 575-598. [DOI: https://dx.doi.org/10.1086/375509]

211. Abbott, R.; Abbott, T.D.; Abraham, S.; Acernese, F.; Ackley, K.; Adams, A.; Adams, C.; Adhikari, R.X.; Adya, V.B.; Affeldt, C. et al. Search for anisotropic gravitational-wave backgrounds using data from Advanced LIGO and Advanced Virgo’s first three observing runs. Phys. Rev. D; 2021; 104, 022005. [DOI: https://dx.doi.org/10.1103/PhysRevD.104.022005]

212. Drasco, S.; Flanagan, E.E. Detection methods for non-Gaussian gravitational wave stochastic backgrounds. Phys. Rev.; 2003; D67, 082003. [DOI: https://dx.doi.org/10.1103/PhysRevD.67.082003]

213. Abbott, R.; Abbott, T.D.; Acernese, F.; Ackley, K.; Adams, C.; Adhikari, N.; Adhikari, R.X.; Adya, V.B.; Affeldt, C.; Agarwal, D. et al. GWTC-3: Compact Binary Coalescences Observed by LIGO and Virgo During the Second Part of the Third Observing Run. arXiv; 2021; arXiv: 2111.03606

214. Thrane, E. Measuring the non-Gaussian stochastic gravitational-wave background: A method for realistic interferometer data. Phys. Rev. D; 2013; 87, 043009. [DOI: https://dx.doi.org/10.1103/PhysRevD.87.043009]

215. Smith, R.; Thrane, E. Optimal Search for an Astrophysical Gravitational-Wave Background. Phys. Rev. X; 2018; 8, 021019. [DOI: https://dx.doi.org/10.1103/PhysRevX.8.021019]

216. Lawrence, J.; Turbang, K.; Matas, A.; Renzini, A.; Van Remortel, N.; Romano, J. A stochastic-signal-based search for intermittent gravitational-wave backgrounds. 2022; in preparation

217. Banagiri, S.; Mandic, V.; Scarlata, C.; Yang, K.Z. Measuring angular N-point correlations of binary black hole merger gravitational-wave events with hierarchical Bayesian inference. Phys. Rev. D; 2020; 102, 063007. [DOI: https://dx.doi.org/10.1103/PhysRevD.102.063007]

218. Barsotti, L.; Fritschel, P.; Evans, M.; Gras, S. Advanced LIGO Design Sensitivity Curve. 2018; Available online: https://dcc.ligo.org/LIGO-T1800044/public (accessed on 1 November 2021).

219. Thrane, E.; Romano, J.D. Sensitivity curves for searches for gravitational-wave backgrounds. Phys. Rev.; 2013; D88, 124032. [DOI: https://dx.doi.org/10.1103/PhysRevD.88.124032]

220. Coughlin, M.; Harms, J. Constraining the gravitational wave energy density of the Universe using Earth’s ring. Phys. Rev. D; 2014; 90, 042005. [DOI: https://dx.doi.org/10.1103/PhysRevD.90.042005]

221. Armstrong, J.W.; Iess, L.; Tortora, P.; Bertotti, B. Stochastic gravitational wave background: Upper limits in the 10**-6-Hz 10**-3-Hz band. Astrophys. J.; 2003; 599, pp. 806-813. [DOI: https://dx.doi.org/10.1086/379505]

222. Lasky, P.D.; Mingarelli, C.M.; Smith, T.L.; Giblin, J.T., Jr.; Thrane, E.; Reardon, D.J.; Caldwell, R.; Bailes, M.; Bhat, N.R.; Burke-Spolaor, S. et al. Gravitational-wave cosmology across 29 decades in frequency. Phys. Rev. X; 2016; 6, 011035. [DOI: https://dx.doi.org/10.1103/PhysRevX.6.011035]

223. Akrami, Y.; Arroja, F.; Ashdown, M.; Aumont, J.; Baccigalupi, C.; Ballardini, M.; Banday, A.J.; Barreiro, R.B.; Bartolo, N.; Basak, S. et al. Planck 2018 results. X. Constraints on inflation. Astron. Astrophys.; 2020; 641, A10. [DOI: https://dx.doi.org/10.1051/0004-6361/201833887]

224. Badurina, L.; Bentine, E.; Blas, D.; Bongs, K.; Bortoletto, D.; Bowcock, T.; Bridges, K.; Bowden, W.; Buchmueller, O.; Burrage, C. et al. AION: An Atom Interferometer Observatory and Network. JCAP; 2020; 05, 011. [DOI: https://dx.doi.org/10.1088/1475-7516/2020/05/011]

225. Blas, D.; Jenkins, A.C. Detecting stochastic gravitational waves with binary resonance. arXiv; 2021; arXiv: 2107.04063

226. Blas, D.; Jenkins, A.C. Bridging the μHz gap in the gravitational-wave landscape with binary resonance. arXiv; 2021; arXiv: 2107.04601

227. Janssen, G.H.; Hobbs, G.; McLaughlin, M.; Bassa, C.G.; Deller, A.T.; Kramer, M.; Lee, K.J.; Mingarelli, C.M.F.; Rosado, P.A.; Sanidas, S. et al. Gravitational wave astronomy with the SKA. PoS; 2015; AASKA14, 037. [DOI: https://dx.doi.org/10.22323/1.215.0037]

228. Abbott, B.; Abbott, R.; Adhikari, R.; Ageev, A.; Allen, B.; Amin, R.; Anderson, S.B.; Anderson, W.G.; Araya, M.; Armandula, H. et al. Analysis of first LIGO science data for stochastic gravitational waves. Phys. Rev. D; 2004; 69, 122004. [DOI: https://dx.doi.org/10.1103/PhysRevD.69.122004]

229. Abbott, B.; Abbott, R.; Adhikari, R.; Agresti, J.; Ajith, P.; Allen, B.; Allen, J.; Amin, R.; Anderson, S.B.; Anderson, W.G. et al. Upper limits on a stochastic background of gravitational waves. Phys. Rev. Lett.; 2005; 95, 221101. [DOI: https://dx.doi.org/10.1103/PhysRevLett.95.221101] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/16384203]

230. Abbott, B.; Abbott, R.; Adhikari, R.; Agresti, J.; Ajith, P.; Allen, B.; Amin, R.; Anderson, S.B.; Anderson, W.G.; Araya, M. et al. Searching for a Stochastic Background of Gravitational Waves with LIGO. Astrophys. J.; 2007; 659, pp. 918-930. [DOI: https://dx.doi.org/10.1086/511329]

231. Aasi, J.; Abbott, B.P.; Abbott, R.; Abbott, T.; Abernathy, M.R.; Accadia, T.; Acernese, F.; Ackley, K.; Adams, C.; Adams, T. et al. Improved Upper Limits on the Stochastic Gravitational-Wave Background from 2009–2010 LIGO and Virgo Data. Phys. Rev. Lett.; 2014; 113, 231101. [DOI: https://dx.doi.org/10.1103/PhysRevLett.113.231101]

232. Abbott, B.P.; Abbott, R.; Abbott, T.D.; Abernathy, M.R.; Acernese, F.; Ackley, K.; Adams, C.; Adams, T.; Addesso, P.; Adhikari, R.X. et al. Upper Limits on the Stochastic Gravitational-Wave Background from Advanced LIGO’s First Observing Run. Phys. Rev. Lett.; 2017; 118, 121101.Erratum on 2017, 119, 029901 [DOI: https://dx.doi.org/10.1103/PhysRevLett.118.121101]

233. Abbott, B.; Abbott, R.; Adhikari, R.; Agresti, J.; Ajith, P.; Allen, B.; Amin, R.; Anderson, S.B.; Anderson, W.G.; Arain, M. et al. Upper limit map of a background of gravitational waves. Phys. Rev. D; 2007; 76, 082003. [DOI: https://dx.doi.org/10.1103/PhysRevD.76.082003]

234. Abadie, J.; Abbott, B.P.; Abbott, R.; Abernathy, M.; Accadia, T.; Acernese, F.; Adams, C.; Adhikari, R.; Ajith, P.; Allen, B. et al. Directional Limits on Persistent Gravitational Waves Using LIGO S5 Science Data. Phys. Rev. Lett.; 2011; 107, 271102. [DOI: https://dx.doi.org/10.1103/PhysRevLett.107.271102]

235. Abbott, B.P.; Abbott, R.; Abbott, T.D.; Abernathy, M.R.; Acernese, F.; Ackley, K.; Adams, C.; Adams, T.; Addesso, P.; Adhikari, R.X. et al. Directional Limits on Persistent Gravitational Waves from Advanced LIGO’s First Observing Run. Phys. Rev. Lett.; 2017; 118, 121102. [DOI: https://dx.doi.org/10.1103/PhysRevLett.118.121102]

236. Abbott, B.P.; Abbott, R.; Abbott, T.D.; Abraham, S.; Acernese, F.; Ackley, K.; Adams, C.; Adhikari, R.X.; Adya, V.B.; Affeldt, C. et al. Directional limits on persistent gravitational waves using data from Advanced LIGO’s first two observing runs. Phys. Rev. D; 2019; 100, 062001. [DOI: https://dx.doi.org/10.1103/PhysRevD.100.062001]

237. Aasi, J.; Abadie, J.; Abbott, B.P.; Abbott, R.; Abbott, T.; Abernathy, M.R.; Accadia, T.; Acernese, F.; Adams, C.; Adams, T. et al. Constraints on Cosmic Strings from the LIGO-Virgo Gravitational-Wave Detectors. Phys. Rev. Lett.; 2014; 112, 131101. [DOI: https://dx.doi.org/10.1103/PhysRevLett.112.131101]

238. Abbott, B.P.; Abbott, R.; Abbott, T.D.; Acernese, F.; Ackley, K.; Adams, C.; Adams, T.; Addesso, P.; Adhikari, R.X.; Adya, V.B. et al. Constraints on cosmic strings using data from the first Advanced LIGO observing run. Phys. Rev. D; 2018; 97, 102002. [DOI: https://dx.doi.org/10.1103/PhysRevD.97.102002]

239. Abbott, B.P.; Abbott, R.; Abbott, T.D.; Acernese, F.; Ackley, K.; Adams, C.; Adams, T.; Addesso, P.; Adhikari, R.X.; Adya, V.B. et al. Search for Tensor, Vector, and Scalar Polarizations in the Stochastic Gravitational-Wave Background. Phys. Rev. Lett.; 2018; 120, 201102. [DOI: https://dx.doi.org/10.1103/PhysRevLett.120.201102]

240. Agarwal, D.; Suresh, J.; Mitra, S.; Ain, A. Upper limits on persistent gravitational waves using folded data and the full covariance matrix from Advanced LIGO’s first two observing runs. PRD; 2021; 104, 123018. [DOI: https://dx.doi.org/10.1103/PhysRevD.104.123018]

241. Yang, K.Z.; Mandic, V.; Scarlata, C.; Banagiri, S. Searching for Cross-Correlation Between Stochastic Gravitational Wave Background and Galaxy Number Counts. Mon. Not. R. Astron. Soc.; 2020; 500, pp. 1666-1672. [DOI: https://dx.doi.org/10.1093/mnras/staa3159]

242. Kapadia, S.J.; Lal Pandey, K.; Suyama, T.; Kandhasamy, S.; Ajith, P. Search for the Stochastic Gravitational-wave Background Induced by Primordial Curvature Perturbations in LIGO’s Second Observing Run. Astrophys. J. Lett.; 2021; 910, L4. [DOI: https://dx.doi.org/10.3847/2041-8213/abe86e]

243. Covas, P.B.; Effler, A.; Goetz, E.; Meyers, P.M.; Neunzert, A.; Oliver, M.; Pearlstone, B.L.; Roma, V.J.; Schofield, R.M.; Adya, V.B. et al. Identification and mitigation of narrow spectral artifacts that degrade searches for persistent gravitational waves in the first two observing runs of Advanced LIGO. Phys. Rev. D; 2018; 97, [DOI: https://dx.doi.org/10.1103/PhysRevD.97.082002]

244. LIGO Document T2000384-v4. Available online: https://dcc.ligo.org/T2000384/public (accessed on 1 November 2021).

245. Matas, A.; Dvorkin, I.; Regimbau, T.; Romero, A. Applying Gating to Stochastic Searches in O3. 2021; Available online: https://dcc.ligo.org/P2000546/public (accessed on 1 November 2021).

246. LIGO Scientific CollaborationVirgo CollaborationKAGRA Scientific Collaboration. Data for Upper Limits on the Isotropic Gravitational-Wave Background from Advanced LIGO’s and Advanced Virgo’s Third Observing Run. 2021; Available online: https://dcc.ligo.org/LIGO-G2001287 (accessed on 1 November 2021).

247. Will, C.M. The Confrontation between General Relativity and Experiment. Living Rev. Relativ.; 2014; 17, 4. [DOI: https://dx.doi.org/10.12942/lrr-2014-4]

248. Abbott, B.P.; Abbott, R.; Abbott, T.D.; Abraham, S.; Acernese, F.; Ackley, K.; Adams, C.; Adya, V.B.; Affeldt, C.; Agathos, M. et al. Prospects for observing and localizing gravitational-wave transients with Advanced LIGO, Advanced Virgo and KAGRA. Living Rev. Relativ.; 2020; 23, 3. [DOI: https://dx.doi.org/10.1007/s41114-020-00026-9]

249. Mukherjee, S.; Broadhurst, T.; Diego, J.M.; Silk, J.; Smoot, G.F. Inferring the lensing rate of LIGO-Virgo sources from the stochastic gravitational wave background. Mon. Not. R. Astron. Soc.; 2021; 501, pp. 2451-2466. [DOI: https://dx.doi.org/10.1093/mnras/staa3813]

250. Buscicchio, R.; Moore, C.J.; Pratten, G.; Schmidt, P.; Bianconi, M.; Vecchio, A. Constraining the lensing of binary black holes from their stochastic background. Phys. Rev. Lett.; 2020; 125, 141102. [DOI: https://dx.doi.org/10.1103/PhysRevLett.125.141102]

251. Buscicchio, R.; Moore, C.J.; Pratten, G.; Schmidt, P.; Vecchio, A. Constraining the lensing of binary neutron stars from their stochastic background. Phys. Rev. D; 2020; 102, 081501. [DOI: https://dx.doi.org/10.1103/PhysRevD.102.081501]

252. Abbott, R.; Abbott, T.D.; Abraham, S.; Acernese, F.; Ackley, K.; Adams, A.; Adams, C.; Adhikari, R.X.; Adya, V.B.; Affeldt, C. et al. GWTC-2: Compact Binary Coalescences Observed by LIGO and Virgo During the First Half of the Third Observing Run. Phys. Rev. X; 2021; 11, 021053. [DOI: https://dx.doi.org/10.1103/PhysRevX.11.021053]

253. Ade, P.A.; Aghanim, N.; Armitage-Caplan, C.; Arnaud, M.; Ashdown, M.; Atrio-Barandela, F.; Aumont, J.; Baccigalupi, C.; Banday, A.J.; Barreiro, R.B. et al. Planck 2013 results. XXV. Searches for cosmic strings and other topological defects. Astron. Astrophys.; 2014; 571, A25. [DOI: https://dx.doi.org/10.1051/0004-6361/201321621]

254. Vaskonen, V.; Veermäe, H. Lower bound on the primordial black hole merger rate. Phys. Rev. D; 2020; 101, 043015. [DOI: https://dx.doi.org/10.1103/PhysRevD.101.043015]

255. Carr, B.; Kohri, K.; Sendouda, Y.; Yokoyama, J. Constraints on Primordial Black Holes. Rep. Prog. Phys.; 2021; 84, 116902. [DOI: https://dx.doi.org/10.1088/1361-6633/ac1e31]

256. LIGO Scientific CollaborationVirgo CollaborationKAGRA Scientific Collaboration. Data Products and Supplemental Information for O3 Stochastic Directional Paper. 2021; Available online: https://dcc.ligo.org/LIGO-G2002165 (accessed on 1 November 2021).

257. Abbott, B.; Abbott, R.; Abbott, T.; Acernese, F.; Ackley, K.; Adams, C.; Adams, T.; Addesso, P.; Adhikari, R.; Adya, V. et al. Search for gravitational waves from Scorpius X-1 in the first Advanced LIGO observing run with a hidden Markov model. Phys. Rev. D; 2017; 95, 122003. [DOI: https://dx.doi.org/10.1103/PhysRevD.95.122003]

258. Sun, L.; Melatos, A.; Lasky, P.D.; Chung, C.T.Y.; Darman, N.S. Cross-correlation search for continuous gravitational waves from a compact object in SNR 1987A in LIGO Science run 5. Phys. Rev. D; 2016; 94, 082004. [DOI: https://dx.doi.org/10.1103/PhysRevD.94.082004]

259. Aasi, J.; Abadie, J.; Abbott, B.P.; Abbott, R.; Abbott, T.; Abernathy, M.R.; Accadia, T.; Acernese, F.; Adams, C.; Adams, T. et al. Directed search for continuous gravitational waves from the Galactic center. Phys. Rev. D; 2013; 88, 102002. [DOI: https://dx.doi.org/10.1103/PhysRevD.88.102002]

260. Abbott, R.; Abbott, T.D.; Acernese, F.; Ackley, K.; Adams, C.; Adhikari, N.; Adhikari, R.X.; Adya, V.B.; Affeldt, C.; Agarwal, D. et al. All-sky, all-frequency directional search for persistent gravitational-waves from Advanced LIGO’s and Advanced Virgo’s first three observing runs. arXiv; 2021; arXiv: 2110.09834

261. Tenorio, R.; Keitel, D.; Sintes, A.M. Application of a hierarchical MCMC follow-up to Advanced LIGO continuous gravitational-wave candidates. Phys. Rev. D; 2021; 104, 084012. [DOI: https://dx.doi.org/10.1103/PhysRevD.104.084012]

262. Burke-Spolaor, S.; Taylor, S.R.; Charisi, M.; Dolch, T.; Hazboun, J.S.; Holgado, A.M.; Kelley, L.Z.; Lazio, T.J.W.; Madison, D.R.; McMann, N. et al. The astrophysics of nanohertz gravitational waves. A&A Rev.; 2019; 27, 5. [DOI: https://dx.doi.org/10.1007/s00159-019-0115-7]

263. Taylor, S.R. The Nanohertz Gravitational Wave Astronomer. arXiv; 2021; arXiv: 2105.13270

264. Sazhin, M.V. Opportunities for detecting ultralong gravitational waves. Sov. Astron.; 1978; 22, pp. 36-38.

265. Detweiler, S. Pulsar timing measurements and the search for gravitational waves. ApJ; 1979; 234, pp. 1100-1104. [DOI: https://dx.doi.org/10.1086/157593]

266. Siemens, X.; Ellis, J.; Jenet, F.; Romano, J.D. The stochastic background: Scaling laws and time to detection for pulsar timing arrays. Class. Quantum Gravity; 2013; 30, 224015. [DOI: https://dx.doi.org/10.1088/0264-9381/30/22/224015]

267. Bailes, M.; Jameson, A.; Abbate, F.; Barr, E.D.; Bhat, N.D.R.; Bondonneau, L.; Burgay, M.; Buchner, S.J.; Camilo, F.; Champion, D.J. et al. The MeerKAT telescope as a pulsar facility: System verification and early science results from MeerTime. PASA; 2020; 37, e028. [DOI: https://dx.doi.org/10.1017/pasa.2020.19]

268. Jonas, J.L. MeerKAT - The South African Array With Composite Dishes and Wide-Band Single Pixel Feeds. IEEE Proc.; 2009; 97, pp. 1522-1530. [DOI: https://dx.doi.org/10.1109/JPROC.2009.2020713]

269. Jenet, F.A.; Hobbs, G.B.; van Straten, W.; Manchester, R.N.; Bailes, M.; Verbiest, J.P.W.; Edwards, R.T.; Hotan, A.W.; Sarkissian, J.M.; Ord, S.M. Upper Bounds on the Low-Frequency Stochastic Gravitational Wave Background from Pulsar Timing Observations: Current Limits and Future Prospects. ApJ; 2006; 653, pp. 1571-1576. [DOI: https://dx.doi.org/10.1086/508702]

270. Demorest, P.B.; Ferdman, R.D.; Gonzalez, M.E.; Nice, D.; Ransom, S.; Stairs, I.H.; Arzoumanian, Z.; Brazier, A.; Burke-Spolaor, S.; Chamberlin, S.J. et al. Limits on the Stochastic Gravitational Wave Background from the North American Nanohertz Observatory for Gravitational Waves. ApJ; 2013; 762, 94. [DOI: https://dx.doi.org/10.1088/0004-637X/762/2/94]

271. Shannon, R.M.; Ravi, V.; Coles, W.A.; Hobbs, G.; Keith, M.J.; Manchester, R.N.; Wyithe, J.S.B.; Bailes, M.; Bhat, N.D.R.; Burke-Spolaor, S. et al. Gravitational-wave limits from pulsar timing constrain supermassive black hole evolution. Science; 2013; 342, pp. 334-337. [DOI: https://dx.doi.org/10.1126/science.1238012]

272. Lentati, L.; Taylor, S.R.; Mingarelli, C.M.F.; Sesana, A.; Sanidas, S.A.; Vecchio, A.; Caballero, R.N.; Lee, K.J.; van Haasteren, R.; Babak, S. et al. European Pulsar Timing Array limits on an isotropic stochastic gravitational-wave background. MNRAS; 2015; 453, pp. 2576-2598. [DOI: https://dx.doi.org/10.1093/mnras/stv1538]

273. NANOGrav Collaboration Arzoumanian, Z.; Brazier, A.; Burke-Spolaor, S.; Chamberlin, S.; Chatterjee, S.; Christy, B.; Cordes, J.M.; Cornish, N.; Crowter, K. et al. The NANOGrav Nine-year Data Set: Observations, Arrival Time Measurements, and Analysis of 37 Millisecond Pulsars. ApJ; 2015; 813, 65. [DOI: https://dx.doi.org/10.1088/0004-637X/813/1/65]

274. Shannon, R.M.; Ravi, V.; Lentati, L.T.; Lasky, P.D.; Hobbs, G.; Kerr, M.; Manchester, R.N.; Coles, W.A.; Levin, Y.; Bailes, M. et al. Gravitational waves from binary supermassive black holes missing in pulsar observations. Science; 2015; 349, pp. 1522-1525. [DOI: https://dx.doi.org/10.1126/science.aab1910] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/26404832]

275. Chen, S.; Caballero, R.N.; Guo, Y.J.; Chalumeau, A.; Liu, K.; Shaifullah, G.; Lee, K.J.; Babak, S.; Desvignes, G.; Parthasarathy, A. et al. Common-red-signal analysis with 24-yr high-precision timing of the European Pulsar Timing Array: Inferences in the stochastic gravitational-wave background search. MNRAS; 2021; 508, pp. 4970-4993. [DOI: https://dx.doi.org/10.1093/mnras/stab2833]

276. Antoniadis, J.; Arzoumanian, Z.; Babak, S.; Bailes, M.; Bak Nielsen, A.S.; Baker, P.T.; Bassa, C.G.; Bécsy, B.; Berthereau, A.; Bonetti, M. et al. The International Pulsar Timing Array second data release: Search for an isotropic gravitational wave background. MNRAS; 2022; 510, pp. 4873-4887. [DOI: https://dx.doi.org/10.1093/mnras/stab3418]

277. Shannon, R.M.; Osłowski, S.; Dai, S.; Bailes, M.; Hobbs, G.; Manchester, R.N.; van Straten, W.; Raithel, C.A.; Ravi, V.; Toomey, L. et al. Limitations in timing precision due to single-pulse shape variability in millisecond pulsars. MNRAS; 2014; 443, pp. 1463-1481. [DOI: https://dx.doi.org/10.1093/mnras/stu1213]

278. Lam, M.T.; McLaughlin, M.A.; Arzoumanian, Z.; Blumer, H.; Brook, P.R.; Cromartie, H.T.; Demorest, P.B.; DeCesar, M.E.; Dolch, T.; Ellis, J.A. et al. The NANOGrav 12.5 yr Data Set: The Frequency Dependence of Pulse Jitter in Precision Millisecond Pulsars. ApJ; 2019; 872, 193. [DOI: https://dx.doi.org/10.3847/1538-4357/ab01cd]

279. Kerr, M.; Reardon, D.J.; Hobbs, G.; Shannon, R.M.; Manchester, R.N.; Dai, S.; Russell, C.J.; Zhang, S.; van Straten, W.; Osłowski, S. et al. The Parkes Pulsar Timing Array project: Second data release. PASA; 2020; 37, e020. [DOI: https://dx.doi.org/10.1017/pasa.2020.11]

280. Lentati, L.; Shannon, R.M.; Coles, W.A.; Verbiest, J.P.W.; van Haasteren, R.; Ellis, J.A.; Caballero, R.N.; Manchester, R.N.; Arzoumanian, Z.; Babak, S. et al. From spin noise to systematics: Stochastic processes in the first International Pulsar Timing Array data release. MNRAS; 2016; 458, pp. 2161-2187. [DOI: https://dx.doi.org/10.1093/mnras/stw395]

281. Goncharov, B.; Zhu, X.J.; Thrane, E. Is there a spectral turnover in the spin noise of millisecond pulsars?. MNRAS; 2020; 497, pp. 3264-3272. [DOI: https://dx.doi.org/10.1093/mnras/staa2081]

282. Chalumeau, A.; Babak, S.; Petiteau, A.; Chen, S.; Samajdar, A.; Caballero, R.N.; Theureau, G.; Guillemot, L.; Desvignes, G.; Parthasarathy, A. et al. Noise analysis in the European Pulsar Timing Array data release 2 and its implications on the gravitational-wave background search. MNRAS; 2022; 509, pp. 5538-5558. [DOI: https://dx.doi.org/10.1093/mnras/stab3283]

283. Haskell, B.; Melatos, A. Models of pulsar glitches. Int. J. Mod. Phys. D; 2015; 24, 1530008. [DOI: https://dx.doi.org/10.1142/S0218271815300086]

284. Hazboun, J.S.; Simon, J.; Siemens, X.; Romano, J.D. Model Dependence of Bayesian Gravitational-wave Background Statistics for Pulsar Timing Arrays. ApJ; 2020; 905, L6. [DOI: https://dx.doi.org/10.3847/2041-8213/abca92]

285. Goncharov, B.; Reardon, D.J.; Shannon, R.M.; Zhu, X.J.; Thrane, E.; Bailes, M.; Bhat, N.D.R.; Dai, S.; Hobbs, G.; Kerr, M. et al. Identifying and mitigating noise sources in precision pulsar timing data sets. MNRAS; 2021; 502, pp. 478-493. [DOI: https://dx.doi.org/10.1093/mnras/staa3411]

286. Meyers, P.M.; O’Neill, N.J.; Melatos, A.; Evans, R.J. Rapid parameter estimation of a two-component neutron star model with spin wandering using a Kalman filter. MNRAS; 2021; 506, pp. 3349-3363. [DOI: https://dx.doi.org/10.1093/mnras/stab1952]

287. Tiburzi, C.; Hobbs, G.; Kerr, M.; Coles, W.A.; Dai, S.; Manchester, R.N.; Possenti, A.; Shannon, R.M.; You, X.P. A study of spatial correlations in pulsar timing array data. MNRAS; 2016; 455, pp. 4339-4350. [DOI: https://dx.doi.org/10.1093/mnras/stv2143]

288. Guo, Y.J.; Li, G.Y.; Lee, K.J.; Caballero, R.N. Studying the Solar system dynamics using pulsar timing arrays and the LINIMOSS dynamical model. MNRAS; 2019; 489, pp. 5573-5581. [DOI: https://dx.doi.org/10.1093/mnras/stz2515]

289. Vallisneri, M.; Taylor, S.R.; Simon, J.; Folkner, W.M.; Park, R.S.; Cutler, C.; Ellis, J.A.; Lazio, T.J.W.; Vigeland, S.J.; Aggarwal, K. et al. Modeling the Uncertainties of Solar System Ephemerides for Robust Gravitational-wave Searches with Pulsar-timing Arrays. ApJ; 2020; 893, 112. [DOI: https://dx.doi.org/10.3847/1538-4357/ab7b67]

290. Hazboun, J.S.; Romano, J.D.; Smith, T.L. Realistic sensitivity curves for pulsar timing arrays. PRD; 2019; 100, 104028. [DOI: https://dx.doi.org/10.1103/PhysRevD.100.104028]

291. Hazboun, J.; Romano, J.; Smith, T. Hasasia: A Python package for Pulsar Timing Array Sensitivity Curves. J. Open Source Softw.; 2019; 4, 1775. [DOI: https://dx.doi.org/10.21105/joss.01775]

292. Verbiest, J.P.W.; Shaifullah, G.M. Measurement uncertainty in pulsar timing array experiments. Class. Quantum Gravity; 2018; 35, 133001. [DOI: https://dx.doi.org/10.1088/1361-6382/aac412]

293. Lommen, A.N. Pulsar timing arrays: The promise of gravitational wave detection. Rep. Prog. Phys.; 2015; 78, 124901. [DOI: https://dx.doi.org/10.1088/0034-4885/78/12/124901] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/26564968]

294. Sesana, A.; Vecchio, A. Measuring the parameters of massive black hole binary systems with pulsar timing array observations of gravitational waves. PRD; 2010; 81, 104008. [DOI: https://dx.doi.org/10.1103/PhysRevD.81.104008]

295. Mingarelli, C.M.F.; Lazio, T.J.W.; Sesana, A.; Greene, J.E.; Ellis, J.A.; Ma, C.P.; Croft, S.; Burke-Spolaor, S.; Taylor, S.R. The local nanohertz gravitational-wave landscape from supermassive black hole binaries. Nat. Astron.; 2017; 1, pp. 886-892. [DOI: https://dx.doi.org/10.1038/s41550-017-0299-6]

296. Arzoumanian, Z.; Baker, P.T.; Blumer, H.; Bécsy, B.; Brazier, A.; Brook, P.R.; Burke-Spolaor, S.; Charisi, M.; Chatterjee, S.; Chen, S. et al. Searching for Gravitational Waves from Cosmological Phase Transitions with the NANOGrav 12.5-Year Dataset. PRL; 2021; 127, 251302. [DOI: https://dx.doi.org/10.1103/PhysRevLett.127.251302] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/35029450]

297. Xue, X.; Bian, L.; Shu, J.; Yuan, Q.; Zhu, X.; Bhat, N.D.R.; Dai, S.; Feng, Y.; Goncharov, B.; Hobbs, G. et al. Constraining Cosmological Phase Transitions with the Parkes Pulsar Timing Array. PRL; 2021; 127, 251303. [DOI: https://dx.doi.org/10.1103/PhysRevLett.127.251303] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/35029430]

298. Moore, C.J.; Vecchio, A. Ultra-low-frequency gravitational waves from cosmological and astrophysical processes. Nat. Astron.; 2021; 5, pp. 1268-1274. [DOI: https://dx.doi.org/10.1038/s41550-021-01489-8]

299. Chen, Z.C.; Yuan, C.; Huang, Q.G. Non-tensorial gravitational wave background in NANOGrav 12.5-year data set. Sci. China Physics, Mech. Astron.; 2021; 64, 120412. [DOI: https://dx.doi.org/10.1007/s11433-021-1797-y]

300. Chen, Z.C.; Wu, Y.M.; Huang, Q.G. Searching for Isotropic Stochastic Gravitational-Wave Background in the International Pulsar Timing Array Second Data Release. arXiv; 2021; arXiv: 2109.00296

301. Wu, Y.M.; Chen, Z.C.; Huang, Q.G. Constraining the Polarization of Gravitational Waves with the Parkes Pulsar Timing Array Second Data Release. arXiv; 2021; arXiv: 2108.10518[DOI: https://dx.doi.org/10.3847/1538-4357/ac35cc]

302. Arzoumanian, Z.; Baker, P.T.; Blumer, H.; Bécsy, B.; Brazier, A.; Brook, P.R.; Burke-Spolaor, S.; Charisi, M.; Chatterjee, S.; Chen, S. et al. The NANOGrav 12.5-year Data Set: Search for Non-Einsteinian Polarization Modes in the Gravitational-wave Background. ApJ; 2021; 923, L22. [DOI: https://dx.doi.org/10.3847/2041-8213/ac401c]

303. Aggarwal, N.; Aguiar, O.D.; Bauswein, A.; Cella, G.; Clesse, S.; Cruise, A.M.; Domcke, V.; Figueroa, D.G.; Geraci, A.; Goryachev, M. et al. Challenges and opportunities of gravitational-wave searches at MHz to GHz frequencies. Living Rev. Relativ.; 2021; 24, 4. [DOI: https://dx.doi.org/10.1007/s41114-021-00032-5]

304. Weber, J. Detection and Generation of Gravitational Waves. Phys. Rev.; 1960; 117, pp. 306-313. [DOI: https://dx.doi.org/10.1103/PhysRev.117.306]

305. Weber, J. Evidence for discovery of gravitational radiation. Phys. Rev. Lett.; 1969; 22, pp. 1320-1324. [DOI: https://dx.doi.org/10.1103/PhysRevLett.22.1320]

306. Aguiar, O.D. The Past, Present and Future of the Resonant-Mass Gravitational Wave Detectors. Res. Astron. Astrophys.; 2011; 11, pp. 1-42. [DOI: https://dx.doi.org/10.1088/1674-4527/11/1/001]

307. Harms, J.; Ambrosino, F.; Angelini, L.; Braito, V.; Branchesi, M.; Brocato, E.; Cappellaro, E.; Coccia, E.; Coughlin, M.; Della Ceca, R. et al. Lunar Gravitational-wave Antenna. Astrophys. J.; 2021; 910, 1. [DOI: https://dx.doi.org/10.3847/1538-4357/abe5a7]

308. Jani, K.; Loeb, A. Gravitational-Wave Lunar Observatory for Cosmology. J. Cosmol. Astropart. Phys.; 2021; 2021, 044. [DOI: https://dx.doi.org/10.1088/1475-7516/2021/06/044]

309. Sesana, A.; Korsakova, N.; Sedda, M.A.; Baibhav, V.; Barausse, E.; Barke, S.; Berti, E.; Bonetti, M.; Capelo, P.R.; Caprini, C. et al. Unveiling the gravitational universe at μ-Hz frequencies. Exper. Astron.; 2021; 51, pp. 1333-1383. [DOI: https://dx.doi.org/10.1007/s10686-021-09709-9]

310. Book, L.G.; Flanagan, E.E. Astrometric Effects of a Stochastic Gravitational Wave Background. Phys. Rev. D; 2011; 83, 024024. [DOI: https://dx.doi.org/10.1103/PhysRevD.83.024024]

311. Moore, C.J.; Mihaylov, D.P.; Lasenby, A.; Gilmore, G. Astrometric Search Method for Individually Resolvable Gravitational Wave Sources with Gaia. Phys. Rev. Lett.; 2017; 119, 261102. [DOI: https://dx.doi.org/10.1103/PhysRevLett.119.261102]

312. Garcia-Bellido, J.; Murayama, H.; White, G. Exploring the Early Universe with Gaia and THEIA. J. Cosmol. Astropart. Phys.; 2021; 2021, 023. [DOI: https://dx.doi.org/10.1088/1475-7516/2021/12/023]

313. Graham, P.W.; Hogan, J.M.; Kasevich, M.A.; Rajendran, S.; Romani, R.W. Mid-band gravitational wave detection with precision atomic sensors. arXiv; 2017; arXiv: 1711.02225

314. Abe, M.; Adamson, P.; Borcean, M.; Bortoletto, D.; Bridges, K.; Carman, S.P.; Chattopadhyay, S.; Coleman, J.; Curfman, N.M.; DeRose, K. et al. Matter-wave Atomic Gradiometer Interferometric Sensor (MAGIS-100). Quantum Sci. Technol.; 2021; 6, 044003. [DOI: https://dx.doi.org/10.1088/2058-9565/abf719]

315. Vitale, S.; Farr, W.M.; Ng, K.K.Y.; Rodriguez, C.L. Measuring the Star Formation Rate with Gravitational Waves from Binary Black Holes. ApJ; 2019; 886, L1. [DOI: https://dx.doi.org/10.3847/2041-8213/ab50c0]

316. Regimbau, T.; Evans, M.; Christensen, N.; Katsavounidis, E.; Sathyaprakash, B.; Vitale, S. Digging deeper: Observing primordial gravitational waves below the binary black hole produced stochastic background. Phys. Rev. Lett.; 2017; 118, 151105. [DOI: https://dx.doi.org/10.1103/PhysRevLett.118.151105] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/28452538]

317. Ackley, K.; Adya, V.B.; Agrawal, P.; Altin, P.; Ashton, G.; Bailes, M.; Baltinas, E.; Barbuio, A.; Beniwal, D.; Blair, C. et al. Neutron Star Extreme Matter Observatory: A kilohertz-band gravitational-wave detector in the global network. Publ. Astron. Soc. Austral.; 2020; 37, e047. [DOI: https://dx.doi.org/10.1017/pasa.2020.39]

318. Maggiore, M.; Van Den Broeck, C.; Bartolo, N.; Belgacem, E.; Bertacca, D.; Bizouard, M.A.; Branchesi, M.; Clesse, S.; Foffa, S.; García-Bellido, J. et al. Science Case for the Einstein Telescope. JCAP; 2020; 03, 050. [DOI: https://dx.doi.org/10.1088/1475-7516/2020/03/050]

319. Cutler, C.; Harms, J. BBO and the neutron-star-binary subtraction problem. Phys. Rev. D; 2006; 73, 042001. [DOI: https://dx.doi.org/10.1103/PhysRevD.73.042001]

320. Sachdev, S.; Regimbau, T.; Sathyaprakash, B.S. Subtracting compact binary foreground sources to reveal primordial gravitational-wave backgrounds. PRD; 2020; 102, 024051. [DOI: https://dx.doi.org/10.1103/PhysRevD.102.024051]

321. Sharma, A.; Harms, J. Searching for cosmological gravitational-wave backgrounds with third-generation detectors in the presence of an astrophysical foreground. PRD; 2020; 102, 063009. [DOI: https://dx.doi.org/10.1103/PhysRevD.102.063009]

322. Raidal, M.; Spethmann, C.; Vaskonen, V.; Veermäe, H. Formation and evolution of primordial black hole binaries in the early universe. J. Cosmol. Astropart. Phys.; 2019; 2019, 018. [DOI: https://dx.doi.org/10.1088/1475-7516/2019/02/018]

323. Mandic, V.; Bird, S.; Cholis, I. Stochastic Gravitational-Wave Background due to Primordial Binary Black Hole Mergers. Phys. Rev. Lett.; 2016; 117, 201102. [DOI: https://dx.doi.org/10.1103/PhysRevLett.117.201102] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/27886472]

324. Ali-Haïmoud, Y.; Kovetz, E.D.; Kamionkowski, M. Merger rate of primordial black-hole binaries. Phys. Rev. D; 2017; 96, 123523. [DOI: https://dx.doi.org/10.1103/PhysRevD.96.123523]

325. De Luca, V.; Franciolini, G.; Pani, P.; Riotto, A. The Minimum Testable Abundance of Primordial Black Holes at Future Gravitational-Wave Detectors. J. Cosmol. Astropart. Phys.; 2021; 2021, 039. [DOI: https://dx.doi.org/10.1088/1475-7516/2021/11/039]

326. Ng, K.K.Y.; Chen, S.; Goncharov, B.; Dupletsa, U.; Borhanian, S.; Branchesi, M.; Harms, J.; Maggiore, M.; Sathyaprakash, B.S.; Vitale, S. On the single-event-based identification of primordial black hole mergers at cosmological distances. arXiv; 2021; arXiv: 2108.07276

327. Biscoveanu, S.; Talbot, C.; Thrane, E.; Smith, R. Measuring the primordial gravitational-wave background in the presence of astrophysical foregrounds. Phys. Rev. Lett.; 2020; 125, 241101. [DOI: https://dx.doi.org/10.1103/PhysRevLett.125.241101] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33412041]

328. Mukherjee, S.; Silk, J. Can we distinguish astrophysical from primordial black holes via the stochastic gravitational wave background?. Mon. Not. R. Astron. Soc.; 2021; 506, pp. 3977-3985. [DOI: https://dx.doi.org/10.1093/mnras/stab1932]

329. Mukherjee, S.; Meinema, M.S.P.; Silk, J. Prospects of discovering sub-solar primordial black holes using the stochastic gravitational wave background from third-generation detectors. arXiv; 2021; arXiv: 2107.02181

330. Mukherjee, S.; Silk, J. Fundamental physics using the temporal gravitational wave background. Phys. Rev. D; 2021; 104, 063518. [DOI: https://dx.doi.org/10.1103/PhysRevD.104.063518]

331. Schumann, W.O. Über die strahlungslosen Eigenschwingungen einer leitenden Kugel, die von einer Luftschicht und einer Ionosphärenhülle umgeben ist. Z. Naturforschung Teil A; 1952; 7, pp. 149-154. [DOI: https://dx.doi.org/10.1515/zna-1952-0202]

332. Schumann, W.O.; König, H. Über die Beobachtung von “atmospherics” bei geringsten Frequenzen. Naturwissenschaften; 1954; 41, pp. 183-184. [DOI: https://dx.doi.org/10.1007/BF00638174]

333. Christensen, N. Measuring the stochastic gravitational radiation background with laser interferometric antennas. Phys. Rev. D; 1992; 46, pp. 5250-5266. [DOI: https://dx.doi.org/10.1103/PhysRevD.46.5250] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/10014914]

334. Thrane, E.; Christensen, N.; Schofield, R.M.S. Correlated magnetic noise in global networks of gravitational-wave detectors: Observations and implications. PRD; 2013; 87, 123009. [DOI: https://dx.doi.org/10.1103/PhysRevD.87.123009]

335. Thrane, E.; Christensen, N.; Schofield, R.M.S.; Effler, A. Correlated noise in networks of gravitational-wave detectors: Subtraction and mitigation. PRD; 2014; 90, 023013. [DOI: https://dx.doi.org/10.1103/PhysRevD.90.023013]

336. Coughlin, M.W.; Christensen, N.L.; De Rosa, R.; Fiori, I.; Gołkowski, M.; Guidry, M.; Harms, J.; Kubisz, J.; Kulak, A.; Mlynarczyk, J. et al. Subtraction of correlated noise in global networks of gravitational-wave interferometers. Class. Quantum Gravity; 2016; 33, 224003. [DOI: https://dx.doi.org/10.1088/0264-9381/33/22/224003]

337. Himemoto, Y.; Taruya, A. Impact of correlated magnetic noise on the detection of stochastic gravitational waves: Estimation based on a simple analytical model. PRD; 2017; 96, 022004. [DOI: https://dx.doi.org/10.1103/PhysRevD.96.022004]

338. Coughlin, M.W.; Cirone, A.; Meyers, P.; Atsuta, S.; Boschi, V.; Chincarini, A.; Christensen, N.L.; De Rosa, R.; Effler, A.; Fiori, I. et al. Measurement and subtraction of Schumann resonances at gravitational-wave interferometers. PRD; 2018; 97, 102007. [DOI: https://dx.doi.org/10.1103/PhysRevD.97.102007]

339. Himemoto, Y.; Taruya, A. Correlated magnetic noise from anisotropic lightning sources and the detection of stochastic gravitational waves. PRD; 2019; 100, 082001. [DOI: https://dx.doi.org/10.1103/PhysRevD.100.082001]

340. Himemoto, Y.; Nishizawa, A.; Taruya, A. Impacts of overlapping gravitational-wave signals on the parameter estimation: Toward the search for cosmological backgrounds. PRD; 2021; 104, 044010. [DOI: https://dx.doi.org/10.1103/PhysRevD.104.044010]

341. Janssens, K.; Martinovic, K.; Christensen, N.; Meyers, P.M.; Sakellariadou, M. Impact of Schumann resonances on the Einstein Telescope and projections for the magnetic coupling function. PRD; 2021; 104, 122006. [DOI: https://dx.doi.org/10.1103/PhysRevD.104.122006]

342. Saulson, P.R. Terrestrial gravitational noise on a gravitational wave antenna. PRD; 1984; 30, pp. 732-736. [DOI: https://dx.doi.org/10.1103/PhysRevD.30.732]

343. Hughes, S.A.; Thorne, K.S. Seismic gravity-gradient noise in interferometric gravitational-wave detectors. PRD; 1998; 58, 122002. [DOI: https://dx.doi.org/10.1103/PhysRevD.58.122002]

344. Harms, J. Terrestrial gravity fluctuations. Living Rev. Relativ.; 2019; 22, 6. [DOI: https://dx.doi.org/10.1007/s41114-019-0022-2]

345. Amann, F.; Bonsignorio, F.; Bulik, T.; Bulten, H.J.; Cuccuru, S.; Dassargues, A.; DeSalvo, R.; Fenyvesi, E.; Fidecaro, F.; Fiori, I. et al. Site-selection criteria for the Einstein Telescope. Rev. Sci. Instruments; 2020; 91, 094504. [DOI: https://dx.doi.org/10.1063/5.0018414] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33003778]

346. Driggers, J.C.; Harms, J.; Adhikari, R.X. Subtraction of Newtonian noise using optimized sensor arrays. PRD; 2012; 86, 102001. [DOI: https://dx.doi.org/10.1103/PhysRevD.86.102001]

347. Coughlin, M.; Harms, J.; Christensen, N.; Dergachev, V.; DeSalvo, R.; Kandhasamy, S.; Mandic, V. Wiener filtering with a seismic underground array at the Sanford Underground Research Facility. Class. Quantum Gravity; 2014; 31, 215003. [DOI: https://dx.doi.org/10.1088/0264-9381/31/21/215003]

348. Coughlin, M.; Mukund, N.; Harms, J.; Driggers, J.; Adhikari, R.; Mitra, S. Towards a first design of a Newtonian-noise cancellation system for Advanced LIGO. Class. Quantum Gravity; 2016; 33, 244001. [DOI: https://dx.doi.org/10.1088/0264-9381/33/24/244001]

349. Coughlin, M.W.; Harms, J.; Driggers, J.; McManus, D.J.; Mukund, N.; Ross, M.P.; Slagmolen, B.J.J.; Venkateswara, K. Implications of Dedicated Seismometer Measurements on Newtonian-Noise Cancellation for Advanced LIGO. PRL; 2018; 121, 221104. [DOI: https://dx.doi.org/10.1103/PhysRevLett.121.221104]

350. Badaracco, F.; Harms, J.; Bertolini, A.; Bulik, T.; Fiori, I.; Idzkowski, B.; Kutynia, A.; Nikliborc, K.; Paoletti, F.; Paoli, A. et al. Machine learning for gravitational-wave detection: Surrogate Wiener filtering for the prediction and optimized cancellation of Newtonian noise at Virgo. Class. Quantum Gravity; 2020; 37, 195016. [DOI: https://dx.doi.org/10.1088/1361-6382/abab64]

351. Harms, J.; Hild, S. Passive Newtonian noise suppression for gravitational-wave observatories based on shaping of the local topography. Class. Quantum Gravity; 2014; 31, 185011. [DOI: https://dx.doi.org/10.1088/0264-9381/31/18/185011]

352. Badaracco, F.; Harms, J. Optimization of seismometer arrays for the cancellation of Newtonian noise from seismic body waves. Class. Quantum Gravity; 2019; 36, 145006. [DOI: https://dx.doi.org/10.1088/1361-6382/ab28c1]

353. Harms, J.; Venkateswara, K. Newtonian-noise cancellation in large-scale interferometric GW detectors using seismic tiltmeters. Class. Quantum Gravity; 2016; 33, 234001. [DOI: https://dx.doi.org/10.1088/0264-9381/33/23/234001]

354. Creighton, T. Tumbleweeds and airborne gravitational noise sources for LIGO. Class. Quantum Gravity; 2008; 25, 125011. [DOI: https://dx.doi.org/10.1088/0264-9381/25/12/125011]

355. Fiorucci, D.; Harms, J.; Barsuglia, M.; Fiori, I.; Paoletti, F. Impact of infrasound atmospheric noise on gravity detectors used for astrophysical and geophysical applications. PRD; 2018; 97, 062003. [DOI: https://dx.doi.org/10.1103/PhysRevD.97.062003]

356. Muratore, M.; Vetrugno, D.; Vitale, S. Revisitation of time delay interferometry combinations that suppress laser noise in LISA. Class. Quantum Gravity; 2020; 37, 185019. [DOI: https://dx.doi.org/10.1088/1361-6382/ab9d5b]

357. Estabrook, F.B.; Tinto, M.; Armstrong, J.W. Time-delay analysis of LISA gravitational wave data: Elimination of spacecraft motion effects. Phys. Rev. D; 2000; 62, 042002. [DOI: https://dx.doi.org/10.1103/PhysRevD.62.042002]

358. Prince, T.A.; Tinto, M.; Larson, S.L.; Armstrong, J. The LISA optimal sensitivity. Phys. Rev. D; 2002; 66, 122002. [DOI: https://dx.doi.org/10.1103/PhysRevD.66.122002]

359. Sylvestre, J.; Tinto, M. Noise characterization for LISA. Phys. Rev. D; 2003; 68, 102002. [DOI: https://dx.doi.org/10.1103/PhysRevD.68.102002]

360. Crowder, J.; Cornish, N.J. LISA source confusion. Phys. Rev. D; 2004; 70, 082004. [DOI: https://dx.doi.org/10.1103/PhysRevD.70.082004]

361. Chen, Z.C.; Huang, F.; Huang, Q.G. Stochastic Gravitational-wave Background from Binary Black Holes and Binary Neutron Stars and Implications for LISA. Astrophys. J.; 2019; 871, 97. [DOI: https://dx.doi.org/10.3847/1538-4357/aaf581]

362. Caprini, C.; Chala, M.; Dorsch, G.C.; Hindmarsh, M.; Huber, S.J.; Konstandin, T.; Kozaczuk, J.; Nardini, G.; No, J.M.; Rummukainen, K. et al. Detecting gravitational waves from cosmological phase transitions with LISA: An update. JCAP; 2020; 03, 024. [DOI: https://dx.doi.org/10.1088/1475-7516/2020/03/024]

363. Boileau, G.; Jenkins, A.C.; Sakellariadou, M.; Meyer, R.; Christensen, N. Ability of LISA to detect a gravitational-wave background of cosmological origin: The cosmic string case. Phys. Rev. D; 2022; 105, 023510. [DOI: https://dx.doi.org/10.1103/PhysRevD.105.023510]

364. Crowder, J.; Cornish, N. A Solution to the Galactic Foreground Problem for LISA. Phys. Rev. D; 2007; 75, 043008. [DOI: https://dx.doi.org/10.1103/PhysRevD.75.043008]

365. LISA Data Challenge Manual. Available online: https://lisa-ldc.lal.in2p3.fr/static/data/pdf/LDC-manual-Sangria.pdf (accessed on 1 November 2021).