Towards an accurate multi-tracer view of the universe

The upcoming generation of cosmological surveys will map the Universe through multiple observational windows. Galaxies, Gravitational Waves, Neutral Hydrogen, each observable probes the same underlying cosmic structure through different physical processes. Extracting maximal cosmological information from these complementary datasets requires frameworks that respect observational geometry, control systematic effects, and preserve information on the underlying processes that generated them. This thesis develops a theoretical and methodological infrastructure for robust multi-messenger cosmology at the precision frontier. The first part focuses on full-sky galaxy clustering analyses. We present the first complete spherical Fourier-Bessel (SFB) power spectrum implementation with full linear general relativistic corrections without recurring to commonly (mis)used approximations. We demonstrate that neglecting these effects biases primordial non-Gaussianity measurements for Stage IV surveys. Importantly, general relativistic terms produce correlations across widely separated redshifts that provide a clean imprint of General Relativity. We extend the SFB framework to the bispectrum analysis, enabling searches for non-Gaussian signatures without flat-sky approximations. Additionally, we introduce contrastive learning as a framework to align summary statistics with cosmological parameters, discovering interpretable marked statistics that enhance constraints through optimal environment weightings without any dependence on the chosen fiducial cosmology. The second part explores gravitational waves as tracers of large-scale structure across different observational regimes. We explore how cross-correlations between SKA observations and future gravitational wave detections provide constraints on the origin and evolution of stellar-mass black holes, demonstrating that multi-messenger clustering distinguishes binary formation channels and time-delay distributions invisible to single-tracer analyses. We also demonstrate that cross-correlating the nanoHertz gravitational wave background with galaxy surveys provides a clear smoking gun for the astrophysical origin of the gravitational wave signal, likely from supermassive black hole binaries at the center of massive galaxies. We forecast that detecting this imprint is achievable with near-future pulsar timing arrays, and we show how loud individual sources can suppress the signal. However, for the first time we discover and fully characterize a systematic effect that leads to strong biasing of the reconstructed gravitational wave anisotropy maps from pulsar timing arrays. This contamination from unresolved small-scale power inflates angular power spectra by more than an order of magnitude, presenting a critical systematic for future anisotropy studies. The third part delves into novel summary statistics and field-level methods that preserve information structure. We introduce video transformer architectures with spacetime factorized attention, exploring how to embed the spacetime geometry of lightcone observations directly into neural network architectures. This approach respects causality by construction, enabling 21cm density field reconstruction without compression to power spectra and preserving phase information that allows pixel-level topology measurements inaccessible to traditional analyses. Together, these frameworks provide an infrastructure required for extracting cosmological information from heterogeneous multi-messenger datasets. This thesis investigated several improvements to match the accuracy of next-generation experiments while controlling systematics that would otherwise limit their transformative science return.