The First Upper Bound on the Nano-Hertz Gravitational Waves and Galaxy Cross-Correlation signal using 15-year NANOGrav Data and DESI Galaxy Survey

Abstract

The recent detection of a common-spectrum stochastic signal by multiple pulsar timing array (PTA) collaborations has provided tentative evidence for a nanohertz (nHz) stochastic gravitational-wave background (SGWB). This signal can be widely interpreted as originating from a cosmic population of inspiraling supermassive black hole binaries (SMBHBs). Current PTA analyses primarily constrain the SGWB power spectrum and its auto-angular power spectrum. However, the supermassive black holes will produce an underlying correlation with the large-scale structure of the Universe, which can help in understanding the formation and evolution of the binaries. In this work, we develop a new analysis pipeline PyGxGW-PTA for studying the cross-correlation of nHz GW signal with galaxy surveys ($C^{\rm g\, GW}_\ell$) and obtain the first constraint on the SGWB and galaxy distribution cross-correlation using the NANOGrav 15-year dataset in combination with the DESI galaxy catalog. We find no statistically significant correlation between the SGWB and the large-scale distribution of DESI galaxies and using an optimal estimator we put an upper bound on $C^{\rm g\, GW}_{\ell=8} < 0.0083$ at $95\%$ C.I. This yields the first observational upper limit on the spatial correlation between the nHz SGWB and the large-scale structure of the Universe, establishing the observational groundwork for future multi-tracer analyses that will combine PTA data with next-generation galaxy surveys to unveil the SMBHB-galaxy correlation.

Figures (7)

• Figure 1: Schematic overview of the PyGxGW-PTA pipeline used to compute the optimal galaxy--GW cross-correlation estimator. The timing-residual cross-correlations between pulsar pairs are combined to construct the SGWB map, which is then masked along with the galaxy map to compute the angular cross-correlation spectrum, $C_{\ell}^{\mathrm{gGW}}$. The "Pseudo" power spectrum is not equal to the full-sky power spectrum. The mode-coupling matrix is inverted to obtain the "mask-corrected" power spectrum. The multipole moments are optimally combined using their covariance to obtain the final estimator used for setting upper limits on the galaxy--GW correlation.
• Figure 2: Validation of the cross-correlation pipeline using mock data. The dashed black curve shows the injected angular cross-power spectrum, $C_{\ell}^{\rm gGW}\,(\mathrm{Mock})$, constructed by forcing the SGWB anisotropy to follow the DESI galaxy density fluctuations. The solid green curve shows the recovered spectrum obtained after passing the noisy mock timing-residual cross-correlations through the full reconstruction pipeline, with error bars estimated from $\sim 1000$ independent realizations. We simulate 200 isotropically distributed pulsars with the timing-residual noise level set to be 100 times smaller than the median uncertainty of the NANOGrav pulsars to illustrate the recovery of the method. The close agreement between the injected and recovered spectra demonstrates that the estimator reliably reconstructs the underlying galaxy--GW correlation.
• Figure 3: Ensemble-averaged galaxy–GW cross-power spectra obtained from random, statistically independent Gaussian realizations of both the SGWB anisotropy map and the galaxy density map. All maps are masked with the DESI sky coverage, and the reconstruction settings match those used in Fig. \ref{['fig:CgGW_valid']}. The curves show the mean recovered spectrum averaged over 100, 500, and 2000 realizations. As the number of realizations increases, the mean $C_{\ell}^{\rm gGW}$ approaches zero across all multipoles, indicating that the pipeline does not introduce spurious correlations from masking, noise in the map reconstruction, or pulsar geometry. This confirms that the estimator is statistically unbiased in the absence of a true galaxy–GW correlation signal.
• Figure 4: Angular power spectra of the galaxy–GW cross-correlation obtained using the DESI galaxy catalog and the NANOGrav 15-year dataset. The square markers with error bars show the masking-corrected estimates of $\hat{C}_\ell^{\mathrm{C,gGW}}$, after averaging over multipole bins of width $\Delta\ell = 2$, $3$, and $4$. The error bars represent the uncertainties in the measured cross-correlation, derived from Monte Carlo realizations of $C_\ell^{\mathrm{C,gGW}}$ under the null hypothesis of no correlation.
• Figure 5: Logarithm of the mode–coupling matrix, $\log_{10}(M_{\ell_1\ell_2})$, induced by the sky mask. The matrix is computed from the window function of the DESI galaxy catalog, and the elements quantify how the true sky power at multipole $\ell_1$ leaks into the multipole $\ell_2$. The presence of off–diagonal structure encodes leakage introduced by incomplete sky coverage. The color scale represents the logarithmic amplitude of the coupling coefficients.