This paper focuses on two aspects of the statistics of cosmological observables that are important for the next stages of precision cosmology. First, we note that the theory of reduced angular $N$-point spectra has only been developed in detail up to the trispectrum case and in a fashion that makes it difficult to go beyond. To fill this gap, here we present a constructive approach that provides a systematic description of reduced angular $N$-point spectra and their covariance matrices, for arbitrary $N$. Second, we focus on the common practice in the literature on cosmological observables, which consists in simply discarding a part of the expression, namely, the terms containing fields evaluated at the observer position. We point out that this is not justified beyond linear order in perturbation theory, as these terms contribute to all the multipoles of the corresponding spectra and with a magnitude that is of the same order as the rest of the non-linear corrections. We consider the possibility that the reason for neglecting these terms is a conceptual discomfort when using ensemble averages, which originates in an apparent tension between the ergodic hypothesis and the privileged position of the observer on the light-cone. We clarify this subtle issue by performing a careful derivation of the relation between the theoretical statistical predictions and the observational estimators for all $N$. We conclude that there is no inconsistency whatsoever in ensemble-averaging fields at and near the observer position, thus clearing the way for consistent and robust high-precision calculations.

We present the galaxy power spectrum in general relativity. Using a novel approach to computing the power spectrum, we derive the galaxy power spectrum taking into account all the relativistic effects in observations. We first demonstrate that the presence of an infrared divergence in previous calculations of the galaxy power spectrum is inconsistent with the equivalence principle, and we show that this unphysical behavior on large scales vanishes if we consider all the relativistic effects consistently. As a consequence, the relativistic effects in the galaxy power spectrum yield no corrections to the measurement of the primordial non-Gaussianity. We present our numerical computation of the full galaxy power spectrum, with particular emphasis on the deviations from the standard redshift-space power spectrum on large scales. We conclude that, as relativistic effects significantly alter the galaxy power spectrum at $k\lesssim k_{eq}$, they need to be taken into account in the analysis of large-scale data.

We perform theoretical and numerical studies of the full relativistic two-point galaxy correlation function, considering the linear-order scalar and tensor perturbation contributions and the wide-angle effects. Using the gauge-invariant relativistic description of galaxy clustering and accounting for the contributions at the observer position, we demonstrate that the complete theoretical expression is devoid of any long-mode contributions from scalar or tensor perturbations and it lacks the infrared divergences in agreement with the equivalence principle. By showing that the gravitational potential contribution to the correlation function converges in the infrared, our study justifies an IR cut-off $(k_{\text{IR}} \leq H_0)$ in computing the gravitational potential contribution. Using the full gauge-invariant expression, we numerically compute the galaxy two-point correlation function and study the individual contributions in the conformal Newtonian gauge. We find that the terms at the observer position such as the coordinate lapses and the observer velocity (missing in the standard formalism) dominate over the other relativistic contributions in the conformal Newtonian gauge such as the source velocity, the gravitational potential, the integrated Sachs-Wolf effect, the Shapiro time-delay and the lensing convergence. Compared to the standard Newtonian theoretical predictions that consider only the density fluctuation and redshift-space distortions, the relativistic effects in galaxy clustering result in a few percent-level systematic errors beyond the scale of the baryonic acoustic oscillation. Our theoretical and numerical study provides a comprehensive understanding of the relativistic effects in the galaxy two-point correlation function, as it proves the validity of the theoretical prediction and accounts for effects that are often neglected in its numerical evaluation.

The remarkable properties of the geodesic light-cone (GLC) coordinates allow analytic expressions for the light-cone observables, providing a new non-perturbative way for calculating the effects of inhomogeneities in our Universe. However, the gauge-invariance of these expressions in the GLC formalism has not been shown explicitly. Here we provide this missing part of the GLC formalism by proving the gauge-invariance of the GLC expressions for the light-cone observables, such as the observed redshift, the luminosity distance, and the physical area and volume of the observed sources. Our study provides a new insight on the properties of the GLC coordinates and it complements the previous work by the GLC collaboration, leading to a comprehensive description of light propagation in the GLC representation.

Comparing the luminosity distance measurements to its theoretical predictions is one of the cornerstones in establishing the modern cosmology. However, as shown in Biern & Yoo, its theoretical predictions in literature are often plagued with infrared divergences and gauge-dependences. This trend calls into question the sanity of the methods used to derive the luminosity distance. Here we critically investigate four different methods --- the geometric approach, the Sachs approach, the Jacobi mapping approach, and the geodesic light cone (GLC) approach to modeling the luminosity distance, and we present a unified treatment of such methods, facilitating the comparison among the methods and checking their sanity. All of these four methods, if exercised properly, can be used to reproduce the correct description of the luminosity distance.