General and consistent statistics for cosmological observations

TL;DR

This paper develops a systematic approach for reduced angular N-point spectra in cosmology and clarifies the role of observer-position terms, enabling more accurate high-precision statistical analyses of cosmological data.

Contribution

It introduces a constructive method for arbitrary N-point spectra and addresses the justification for including observer-position terms beyond linear order.

Findings

Systematic description of reduced angular N-point spectra for any N.

Observer-position terms contribute significantly beyond linear order.

No inconsistency in ensemble-averaging fields at the observer position.

Abstract

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.

Equations

$$g_o(e_a,e_b)\equiv {\eta }_{ab},$$

$$1+z:=\frac{g(u_s,k_s)}{g(e_0,k_o)},n_i:=-\frac{g(e_i,k_o)}{\sqrt{g(e_j,k_o)g(e_j,k_o)}},$$

$$n\to R^{-1}n.$$

$$\tilde{\mathcal{O}}({\tilde{x}}_o;z,n)\equiv \mathcal{O}(x_o;z,n),$$

$$(k_a)=-\omega \left(1,n_i\right),$$

$$n\equiv \sin \vartheta \cos \phi e_1+\sin \vartheta \sin \phi e_2+\cos \vartheta e_3,$$

$$V_A={\nabla }_{A}v+{\epsilon }_A^B{\nabla }_B\tilde{v},$$

$$T_{AB}=S_{AB}T+{\epsilon }_{AB}\tilde{T}+\left[{\nabla }_A{\nabla }_B-\frac{1}{2}{S}_{AB}{\nabla }^2\right]t+{\epsilon }_{(A}^C{\nabla }_{B)}{\nabla }_C\tilde{t},$$

$$\omega :=-g(e_0,k_o),$$

$$\mathcal{O}(x_o;z,n)\equiv \mathcal{O}(x_o;z,n)[\Phi ].$$

$$\Phi =\bar{\Phi }+\varphi ,$$

$${\varphi }_{\mathfrak{a}}(t,\vec{x})\equiv {\varphi }_{\mathfrak{a}}(t,\vec{x})[{\varphi }_{\mathfrak{b}}(t_0,\vec{y})].$$

$$⟨{\varphi }_{{\mathfrak{a}}_1}(t_0,{\vec{x}}_1)\dots {\varphi }_{{\mathfrak{a}}_n}(t_0,{\vec{x}}_n)⟩:=\int D\varphi (t_0)P[\varphi (t_0)]{\varphi }_{{\mathfrak{a}}_1}(t_0,{\vec{x}}_1)\dots {\varphi }_{{\mathfrak{a}}_n}(t_0,{\vec{x}}_n),$$

$$F_{{\mathfrak{a}}_1\dots {\mathfrak{a}}_n}(x_1,\dots ,x_n)$$

$$⟨{\varphi }_{\mathfrak{a}}(x)⟩\equiv 0,$$

$$⟨\Phi (x)⟩\equiv \bar{\Phi }(x).$$

$$G^{\mathrm{th}}\left(x_o;\{z_k,n_k{\}}_{k=1}^N\right):=⟨\prod _{k=1}^N\mathcal{O}(x_o;z_k,n_k)[\Phi ]⟩,$$

$$G^{\mathrm{th}}\left(x_o;\{z_k,n_k{\}}_{k=1}^N\right)\equiv G^{\mathrm{th}}\left(t_o;\{z_k,n_k{\}}_{k=1}^N\right),$$

$$G^{\mathrm{th}}\left(t_o;\{z_k,R{n}_k{\}}_{k=1}^N\right)\equiv G^{\mathrm{th}}\left(t_o;\{z_k,n_k{\}}_{k=1}^N\right).$$

$$⟨{\mathcal{O}}_{l_1{m}_1}{\mathcal{O}}_{l_2{m}_2}⟩={\delta }_{l_1{l}_2}{\delta }_{m_1,-m_2}(-1{)}^{m_2}{C}_{l_1},$$

$$\left\langle{\cal O}_{l_{1}m_{1}}{\cal O}_{l_{2}m_{2}}{\cal O}_{l_{3}m_{3}}\right\rangle=\left(\begin{array}[]{ccc}l_{1}&l_{2}&l_{3}\\ m_{1}&m_{2}&m_{3}\end{array}\right)B_{l_{1}l_{2}l_{3}}\,,$$

$$\hat{C}_{l}:=\frac{1}{2l+1}\,\sum_{m}{\cal O}_{lm}{\cal O}^{*}_{lm}\,,\hskip 28.45274pt\hat{B}_{l_{1}l_{2}l_{3}}=\sum_{m_{1},m_{2},m_{3}}\left(\begin{array}[]{ccc}l_{1}&l_{2}&l_{3}\\ m_{1}&m_{2}&m_{3}\end{array}\right){\cal O}_{l_{1}m_{1}}{\cal O}_{l_{2}m_{2}}{\cal O}_{l_{3}m_{3}}\,,$$

$$⟨{\hat{C}}_l⟩=C_l,⟨{\hat{B}}_{l_1{l}_2{l}_3}⟩=B_{l_1{l}_2{l}_3}.$$

$$⟨{\mathcal{O}}_{l_1{m}_1}\dots {\mathcal{O}}_{l_{\mathcal{N}}{m}_{\mathcal{N}}}⟩=\sum _{m_k^{\prime }}⟨{\mathcal{O}}_{l_1{m}_1^{\prime }}\dots {\mathcal{O}}_{l_{\mathcal{N}}{m}_{\mathcal{N}}^{\prime }}⟩\times D_{l_1,m_1^{\prime }{m}_1}\dots D_{l_N,m_N^{\prime }{m}_N},$$

$$\frac{\int \mathrm{d}R{\prod }_{k=1}^N\mathcal{O}(x_o;z_k,R^{-1}{n}_k)}{\int \mathrm{d}R},$$

$$R(\alpha ,\beta ,\gamma ):=R_3(\alpha )R_2(\beta )R_3(\gamma ),$$

$$\alpha \in [0,2\pi [,\beta \in [0,\pi ],\gamma \in [0,2\pi [.$$

$$R^{-1}(\alpha ,\beta ,\gamma )\equiv R(-\gamma ,-\beta ,-\alpha ).$$

$$\mathrm{d}R(\alpha ,\beta ,\gamma )\equiv \sin \beta \mathrm{d}\alpha \mathrm{d}\beta \mathrm{d}\gamma ,$$

$${⟨f(n_1,\dots ,n_N)⟩}_{\mathrm{SO}(3)}$$

Full text

1 Introduction

Upcoming surveys such as [1, 2, 3, 4, 5] will significantly enhance the quantity and quality of the data that we have at our disposal for understanding the universe. This development requires a commensurate effort on the theoretical side for the new observational input to be interpreted correctly. In particular, one needs to consider analytical expressions for cosmological observables beyond the linear order in perturbation theory, as is already clearly reflected by the extensive amount of work on the subject in the case of CMB lensing [6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17], galaxy number counts [18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32] and cosmological distances and weak lensing [33, 34, 35, 36, 37, 38, 39, 40, 41]. The inclusion of non-linear corrections is not only required for increasing the precision of the basic quantities of interest, such as the power spectrum of observables, but also for obtaining the leading contributions to higher-order statistics (e.g. bispectrum, trispectrum, etc.). The latter contain crucial information about the early universe and the formation of structure, which will become more accessible thanks to the expected advances in observations. In this paper we wish to elaborate on two important aspects of the statistics of cosmological observables for the next stages of precision cosmology which have not yet been considered in the literature.

The first issue is the lack of generality in our description of angular higher-order statistics. Indeed, in the literature, the rotationally-invariant (or “reduced”) angular $N$-point spectra of a given observable on the sky are defined only up to the trispectrum ($N=4$) case [42] and built in a rather opaque way that does not generalize easily to arbitrary $N$. The $N=3,4$ cases are well studied and measured [43, 44, 45, 46, 47] in CMB observations, whittling down the parameter spaces of the inflationary models. It is known, however, that the hierarchy between $N$-point spectra is not necessarily trivial, e.g. there exist models with a large trispectrum, but negligible bispectrum [48, 49, 50, 51]. Therefore, the study of $N$-point spectra beyond $N=4$ could be relevant for models with both negligible bispectrum and trispectrum. In any case, it is always desirable to have a systematic theoretical framework for approaching a given problem, which is what we will present here. This is a transparent construction of the reduced angular $N$-point spectra, for arbitrary $N$, their covariance matrices for arbitrary $N$ and $N^{\prime}$, as well as a flurry of useful equations for manipulating them. The central ingredient in this framework will be the generalization of the “triangular” Wigner $3-j$ symbol to a “multilateral” one, leading in particular to a neat diagrammatic representation. Our equations will hold for the idealized case of full sky measurements, but they should not be difficult to generalize to the case of partial covering. Note, finally, that the number of components of the angular spectra increases rapidly with $N$, meaning that the new higher-order statistics we can now work with are practically impossible to store, let alone manipulate, in their entirety. Nevertheless, our construction is the indispensable prerequisite for handling these objects, as it provides the starting point for research on privileged subsets, limits, shapes and efficient estimators thereof.

The second issue we want to discuss concerns the consistency and accuracy of the theoretical predictions, especially at the non-linear orders in perturbation theory which now need to be taken into account. The standard approach to analytical computations of cosmological observables neglects certain terms, namely, those containing fields evaluated at the observer point. As we will discuss, for linear order perturbation theory these observer-dependent contributions have no effect on multipoles $l\geq 2$, but we will argue that they do affect all multipoles at the non-linear level. This is an important point, because none of the works on non-linear effects so far [6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41] considers these observer-dependent terms. Moreover, at the linear level it has already been established that the consideration of the full expressions, i.e. including the observer-dependent terms, has some important advantages, as it leads to results that are gauge-invariant and respect the equivalence principle in $k$-space (absence of spurious infrared divergences) [52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63].

We identify two possible reasons behind the reticence towards considering these observer-dependent terms in the literature. The first one is rooted in a common misconception that observer-dependent terms are synonymous with “monopole”, because they only affect the later in the power spectrum of the linear theory (or the $s$ multipole for spin $s$ observables). Since low multipoles are dominated by cosmic variance, one is then naturally led to ignore observer-dependent terms. We will show that this is an artifact of linear perturbation theory, because observer-dependent terms actually do arise in non-linear corrections to all multipoles, and there is no reason to believe that these contributions are negligible. The second possible reason for ignoring these terms might be a discomfort of conceptual nature. Simply put, it is not clear whether one can consistently ensemble average them to compute correlation functions and spectra, because, contrary to the case of source points, one does not probe several “observer points” observationally for the ergodic hypothesis to apply. A first shaky aspect of this argument is that the ensemble averages that are considered in the standard approach are statistically homogeneous, so this already wipes out the observationally privileged status of the observer point, independently of whether one chooses to include the observer-dependent terms or not. More precisely, we will see that the above objection to the statistical use of observer-dependent terms arises from a misunderstanding of the way the ergodic hypothesis enters the relation between the ensemble averages of the theorist and the geometrical averages of the observer. We will perform a careful derivation of that relation under the standard assumption of statistical homogeneity and isotropy for the stochastic fluctuations, but without any assumptions about which points are included in the ensemble averages. We will then see that the error of this theoretical prediction with respect to the data averages is nothing but cosmic variance. Thus, by carelessly ensemble averaging field products at arbitrary points on the light-cone, independently of how well we can probe them observationally, the theorist is not deviating in any new way from observation other than the fundamental limitation of cosmic variance. Let us also stress that our result is completely general, as we will work with the freshly defined $N$-point spectra for arbitrary $N$. The assumption of full sky measurement will not be problematic here, as the point we wish to make is of qualitative nature. Also, we will only consider the spectra of a single observable, but the formalism generalizes straightforwardly in the case of cross-correlations of different observables.

The paper is organized as follows. In section 2 we discuss some technical aspects regarding cosmological observables and their theoretical statistical treatment. The reader who is familiar with cosmological observables can skip this part, although some of its definitions and equations will be used later on, so it might be useful to take a quick look anyway. In section 3 we discuss the present situation regarding reduced angular $N$-point spectra in the literature and then derive our formulas for defining and handling these objects for arbitrary $N$. Section 4 focuses on the issue of observer-dependent terms. We first discuss in more detail the inconsistencies, and potential inaccuracies, that result from ignoring them in the spectra of observables. Then, using our definition of the angular correlation functions and spectra given in the previous section, we derive the relation between the statistical ones of the theorist and the geometrical ones of the observer. This demonstrates that one is allowed to (and one should) take into account the observer-dependent terms inside ensemble averages. In 5 we derive the expression for the covariance matrix of two arbitrary reduced angular spectra and determine the asymptotic behavior of cosmic variance for large $l$ values. In 6 we discuss a few subtleties one should be aware of when working with the connected parts of the statistics. Finally, in section 7 we conclude. The appendix contains technical derivations to avoid clogging the main text, as well as the set of Wigner $3-j$ symbol identities we will use in this paper.

2 Preliminary considerations

2.1 Angle, redshift, and observables

A cosmological observable associated with some localized source of light (galaxy, supernovae, etc.) is a function of two space-time points ${\cal O}(x_{o};x_{s})$, the “observer” point $x_{o}$ and the “source” point $x_{s}$, that are constrained to lie on a common light-like geodesic, with $x_{s}$ in the past of $x_{o}$. Thus, for a given observer at $x_{o}$, the set of possible light source points $x_{s}$ detected through light forms the past light-cone of $x_{o}$. The numbers $x_{s}^{\mu}$ are not observables, as they are ambiguous due to the freedom of performing coordinate transformations. One can of course consider correlation functions of the form $\langle{\cal O}(x_{o};x_{1})\dots\cal O(x_{o};x_{N})\rangle$, or the associated spectra in $k$-space, but these cannot be related to observable information, especially in the presence of non-negligible inhomogeneity effects. One must therefore work directly with the physically unambiguous relations, i.e. the relations between observables, so one must parametrize the past-light cone at $x_{o}$ in terms of observables. One of them is the incoming photon direction in the sky, $n$, leading to an angular parametrization of the sky. Thus, instead of working with the “$k$-space” spectra, we will work with the “$l$-space” ones, which are directly related to observable quantities. Finally, for the radial parametrization of the light-cone there are several choices, such as the observed redshift $z$ or the luminosity/angular distances $D_{L,A}$. Here we will consider the former, which is also the most widely used and model independent. Thus, the observer at $x_{o}$ parametrizes her light-cone in terms of the observables $z$ and $n$, and the physical information lies in the relation between ${\cal O}$ and $(z,n)$ that is the function ${\cal O}(x_{o};z,n)$.

The quantities $z$ and $n$ are defined with respect to the observer rest-frame, i.e. a tetrad $e_{a}$ at $x_{o}$

[TABLE]

whose time-component $e_{0}$ is the observer’s 4-velocity, and the source 4-velocity $u_{s}$ satisfying $g_{s}(u_{s},u_{s})\equiv-1$. More precisely, if $k$ denotes the momentum 4-vector of the photon, we have

[TABLE]

where $i\in\{1,2,3\}$ is the spatial part of the 4-dimensional index $a\in\{0,1,2,3\}$. Note that $n_{i}$ is the observed angular direction, not the propagation direction, hence the minus sign. These quantities are invariant under space-time coordinate transformations, as any observable should, since a measurement cannot depend on how we parametrize that manifold. In the absence of a tetrad, the only available numbers are the components of $k$ in a coordinate system $k^{\mu}$ and the corresponding angles in the spatial part are the ones an observer uses only if the basis $\partial_{\mu}$ is orthonormal, i.e. only if $g_{\mu\nu}(x_{o})=\eta_{\mu\nu}$. This is not the case in most of the coordinate systems used in cosmology, which is why one requires a tetrad to obtain the correct angles $n$. Since this field is defined at the observer position, the corresponding correction of the angles it induces will appear as observer terms in the expressions of observables.

The observer tetrad (2.1) is defined only up to a Lorentz transformation $e_{a}\to\Lambda_{a}^{\,\,\,b}e_{b}$, but the boosts alter the observer 4-velocity, so the only ambiguity is the orientation of the spatial frame $e_{i}\to R_{i}^{\,\,\,j}e_{j}$, leading to an SO(3) ambiguity for $n$

[TABLE]

As for the observable ${\cal O}$, the theoretical expression must also be a scalar under coordinate transformations [58]

[TABLE]

but it can be a tensor with respect to the Lorentz index $a$ of the observer tetrad. For instance, $\omega$ and $n$ are components of the Lorentz vector $k_{a}:=g(e_{a},k_{o})$ at $x_{o}$

[TABLE]

which is why they transform non-trivially under Lorentz transformations of $e_{a}$. Since the boost part here is fixed by the definite observer 4-velocity, we only have to deal with the SO(3) ambiguity (2.3). To that end, we express $n$ in terms of observed angles $(\vartheta,\varphi)$

[TABLE]

so that ${\cal O}\equiv{\cal O}(x_{o};z,\vartheta,\varphi)$ becomes a function on the unit sphere $S_{2}$. The latter is a 2-dimensional manifold with coordinates $\vartheta^{A}\in\{\vartheta,\varphi\}$ and with the admissible coordinate transformations being the ones induced by the rotation (2.3).111Note that this manifold is a subset of the tangent space at $x_{o}$ and with a parametrization that is induced by the tetrad basis $e_{a}$ of that space, so these angles are not space-time coordinates. The generic observable will therefore be a tensor field on that manifold ${\cal O}_{A_{1}\dots A_{n}}$. However, in two dimensions any tensor can be locally reduced to scalars and pseudo-scalars. For instance, a vector can be decomposed into a scalar $v$ and a pseudo-scalar $\tilde{v}$ via

[TABLE]

while a tensor can be decomposed into two scalars $T,t$ and two pseudo-scalars $\tilde{T},\tilde{t}$

[TABLE]

where $S_{AB}$, $\varepsilon_{AB}$ and $\nabla_{A}$ are the metric, totally anti-symmetric tensor in 2d and covariant derivative on the 2-sphere $S_{2}$, respectively. Further decomposing these (pseudo-)scalars into spherical harmonics leads to the decomposition of ${\cal O}_{A_{1}\dots A_{n}}$ into spin-weighted spherical harmonics. To avoid the introduction of the latter, which complicates unnecessarily the formalism, here we will assume that our observables are the (pseudo-) scalars of the above decomposition (up to possible Laplacians). For instance, in the case of the CMB polarization tensor $P_{AB}$ we would directly work with its “electric” and “magnetic” components $\nabla_{A}\nabla_{B}P^{AB}$ and $\varepsilon_{A}^{\,\,\,C}\nabla_{B}\nabla_{C}P^{AB}$, respectively.

Finally, another important observable is the incoming photon frequency

[TABLE]

so the most general parametrization of light-based observables is a spectral distribution ${\cal O}(x_{o};z,n,\omega)$. In the case of observables associated with diffuse sources (e.g. the CMB), i.e. that do not have a particular emission moment and therefore no associated redshift, then we simply have ${\cal O}(x_{o};n,\omega)$. For the sake of simplicity, here we will assume that $\omega$ is either fixed, or that it is integrated over with some given spectral distribution, as in the case of the CMB temperature for instance. We will therefore work with observable functions of the form ${\cal O}(x_{o};z,n)$, but the inclusion of $\omega$ is straightforward and basically enters our equations exactly as the redshift dependence. A more detailed account of this subsection’s content, and in particular of the underlying geometrical constructions, can be found in [64].

2.2 Ensemble averaged $N$-point correlation functions

The theorist works with tensor fields on the full space-time manifold which we collectively denote by $\Phi(x)$. The observable ${\cal O}(x_{o};z,n)$ is an integro-differential functional of these fields

[TABLE]

Within cosmological perturbation theory, one splits $\Phi$ into a homogeneous and isotropic “background” solution $\bar{\Phi}$ and a fluctuation $\phi$

[TABLE]

where the latter is subject to a gauge ambiguity. This is the manifestation of infinitesimal coordinate transformations in this framework, since the background is kept fixed under these transformation. The observable ${\cal O}(x_{o};z,n)$ should, by definition, be independent of the choice of gauge, so one can fix the latter. Moreover, one can use any constraint equations (e.g. the Poisson equation) to reduce the number of fields down to those carrying degrees of freedom, which we denote by $\phi_{\mathfrak{a}}(x)$, i.e. the “$\mathfrak{a}$” index includes both space-time and internal indices. Thus, the $\phi_{\mathfrak{a}}(t_{0},\vec{x})$ data at a given time $t_{0}$ uniquely determine the ones at any other time $t$

[TABLE]

Here it is assumed that the $\phi_{\mathfrak{a}}$ also contain the momenta/velocities, for fields obeying second-order equations in time. Next, in order to compare with observations, one needs to consider a statistical ensemble of solutions. The $\phi_{\mathfrak{a}}$ are therefore promoted to stochastic fields with a corresponding probability distribution functional (pdf) associating a probability density to each field solution $\phi_{\mathfrak{a}}(x)$. Since the latter are completely determined by their configuration at some reference time $\phi_{\mathfrak{a}}(t_{0},\vec{x})$, it suffices to define that pdf on these field configurations $P\equiv P[\phi(t_{0})]$ (dropping the index $\mathfrak{a}$ and the position $\vec{x}$ for notational simplicity). One can then define the moments, i.e. the statistical averages of field products

[TABLE]

which completely determine the functional $P$, and with these one can define the field correlation functions (FCF)

[TABLE]

and use the linearity of the averaging operation to express this as a functional of the (2.13). In particular, $P$ is chosen such that

[TABLE]

which can be alternatively stated as

[TABLE]

With the above definitions one can now perform statistical averages of arbitrary functionals of the $\phi_{\mathfrak{a}}(x)$. In particular, the theoretical correlation functions (TCF) of the observables are defined by

[TABLE]

which are therefore ultimately a functional of $\bar{\Phi}$ and the FCFs (2.14). Now note that, although strict homogeneity and isotropy are lost as soon as $\phi\neq 0$, these notions can be reintroduced at the statistical level by imposing the corresponding symmetries on $P[\phi(t_{0})]$. Thus, we require that, if the two configurations $\phi_{\mathfrak{a}}(t_{0})$ and $\phi^{\prime}_{\mathfrak{a}}(t_{0})$ are related by an isometry of the background geometry, then $P[\phi(t_{0})]=P[\phi^{\prime}(t_{0})]$. As a result, the FCFs and TCFs are invariant under isometries, meaning in particular that $G^{\rm th}$ is independent of the observer position $\vec{x}_{o}$

[TABLE]

and invariant under a common rotation of its $N$ directions

[TABLE]

3 The general theory of reduced angular $N$-point spectra

3.1 The previous implicit approach

Consider some observable on the sky ${\cal O}(n)$, where for the purposes of the present discussion we can omit the redshift dependence. The quantity that is usually considered in theory is the angular $N$-point spectrum through ensemble averaging $\left\langle{\cal O}_{l_{1}m_{1}}\dots\cal O_{l_{N}m_{N}}\right\rangle$, where ${\cal O}_{lm}$ are the spherical harmonic components. In the standard approach where statistical isotropy is assumed, that $N$-point spectrum is invariant under rotations, so that not all of its components are independent. Its information can therefore be stored in a rotationally-invariant quantity with fewer indices called the “reduced” angular spectrum. Assuming for simplicity here $\left\langle{\cal O}_{lm}\right\rangle=0$, the lowest-order cases are the angular power spectrum

[TABLE]

where $C_{l}$ is the “reduced angular power spectrum”, and the angular bispectrum

[TABLE]

where $B_{l_{1}l_{2}l_{3}}$ is the “reduced angular bispectrum” and the first factor is the Wigner $3-j$ symbol. As one could expect from isotropy, $C_{l}$ and $B_{l_{1}l_{2}l_{3}}$ have no $m_{k}$ dependence. The relations (3.1) and (3.2) can be inverted and, in particular, allow one to build unbiased estimators for these quantities

[TABLE]

i.e.

[TABLE]

One can therefore work exclusively with the invariant and more compact reduced spectra, both at the theoretical and data level, instead of the full ones. The theory of the reduced angular $N$-point spectra is well developed in the literature up to the trispectrum case [65, 43, 44, 45, 46, 66, 67, 68, 47], based on the pioneering work in [42] which considered the following method. First, one looks for the relation between the full spectrum and its reduced counterpart (e.g. (3.1) and (3.2)) by imposing on the former the invariance under rotations

[TABLE]

where $D_{l,mm^{\prime}}$ is the Wigner matrix, all of them here depending on the same rotation angles. Using the identities of the Wigner matrices and $3-j$ symbols, one must then find the most general solution to this equation, which is (3.1) and (3.2) for the $N=2,3$ cases, respectively. Finally, one must also invert this relation to obtain its estimator (e.g. (3.3)) and, through that quantity, its explicit definition in terms of ensemble averages of the observables (e.g. (3.4)). Although this procedure generalizes straightforwardly to $N=4$ [42] and beyond222Following this procedure, in [69] the authors have found the analogue of Eqs. (3.1) and (3.2) for the quadrispectrum case $N=5$. In [70] one can find a formal expression for arbitrary $N$, and also the real space $N$-point correlation function in terms of the reduced angular spectrum, but without rigorous demonstration. In both of these references, however, the inversion to find the estimators and ensemble-averaged spectra in terms of the observables, i.e. the analogues of (3.3) and (3.4), is not given., it becomes quickly very difficult in practice, as one has to “guess” the relation between the full and reduced spectrum by inferring it from an increasingly complicated equation. Moreover, even after finding a solution, one must still prove that it is the most general one and then invert it, which is also harder to do for larger values of $N$. As a result, the generalization of (3.3) for arbitrary $N$ does not exist in the literature. Finally, with this procedure the reduced spectrum is determined only up to an overall normalization, which can depend on the $l_{i}$ numbers, since only the $m_{i}$ numbers are summed over in the defining equation (3.5).333See for instance [70, 66], where alternative normalizations are considered and are essentially related to the $l_{i}$-dependent factors between the Gaunt and $3-j$ coefficients.

3.2 A new constructive approach

Here we propose a computationally straightforward approach to the problem. Since the defining property behind the notion of reduced spectra is rotational invariance, we will first focus on this geometrical aspect, i.e. before any statistical considerations. We will therefore be constructing the generalization of the full estimators (3.3), which exclusively depend on observed data, with their statistical counterparts being obtainable by simply taking the ensemble average as in (3.4). In fact, we will refer to these estimators as the “observational” reduced spectra, as opposed to the “theoretical” ones (2.17) which are obtained though statistical ensembles that are observationally unavailable. This is to highlight the fact that the observational reduced spectra have a significance of their own, as they amount to all the rotationally-invariant information one can extract out of the observed data. This is the physical information, as it is independent of the artificial SO(3) ambiguity of the observer’s spatial reference frame $e_{i}$ in (2.6).

We start with the fact that the observer only has access to a single light-cone $x_{o}$ and therefore to the observables ${\cal O}(x_{o};z,n)$. The building blocks for rotationally-invariant functionals of ${\cal O}(x_{o};z,n)$ are the average of the products $\prod_{k=1}^{N}{\cal O}(x_{o};z_{k},n_{k})$ over all possible common rotations of the directions $n_{k}$, i.e. the average over the SO(3) group

[TABLE]

where ${\rm d}R$ is the Haar measure on SO(3) and is invariant under group multiplication. Thanks to this, if one rotates the $n_{k}$ in (3.6) with the same rotation $R^{\prime}$, then this can be reabsorbed in the dummy variable $R$ by the redefinition $R\to R^{\prime}R$, thus leaving the average invariant. To turn (3.6) into a well-defined integral, we consider an arbitrary orthonormal reference frame $e_{i}$ and parametrize $R$ as

[TABLE]

where $R_{i}(\theta)$ denotes the matrix corresponding to a rotation around $e_{i}$ with angle $\theta$, so that $\alpha,\beta,\gamma$ are the Euler angles

[TABLE]

In particular, the inverse matrix is simply

[TABLE]

The Haar measure ${\rm d}R$ on SO(3) now reads

[TABLE]

so the average over SO(3) of some function $f(n_{1},\dots,n_{N})$ is given by

[TABLE]

and thus, the $N$-point observational correlation functions (OCF) are defined by

[TABLE]

In particular, the $N=1$ case reduces to the average over the sphere

[TABLE]

which is shown by picking $e_{3}=n$. For $N>1$, three out of the $2N$ angles in $G^{\rm ob}$ are redundant, since we are free to rotate at will, or equivalently, to choose the reference frame $\{e_{i}\}_{i=1}^{3}$ arbitrarily. For instance, one can pick (assuming that $n_{1}$ and $n_{2}$ are not parallel)

[TABLE]

thus leaving us with a dependence on the $2N-3$ angles $\vartheta_{2},\dots\vartheta_{N}$ and $\varphi_{3},\dots,\varphi_{N}$ that parametrize $\{n_{k}\}_{k=2}^{N}$ in the $e_{i}$ basis (2.6).

Now note that the $N$-point OCF is a redshift-dependent function on

[TABLE]

that is symmetric under a common rotation of these $N$ spheres, just as its theoretical counterpart (2.19). In particular, expressing ${\cal O}$ inside $G^{\rm ob}$ in terms of the stochastic fields ${\cal O}\equiv{\cal O}[\Phi]$ and using the linearity of averaging, we obtain the identity

[TABLE]

showing that the $G^{\rm ob}$ are full unbiased estimators of the $G^{\rm th}$. The idea now is to decompose such functions in a basis of SO(3)-invariant functions on $S_{2}^{N}$. For $N=2$ the only SO(3)-invariant combination of the two vectors is $n_{1}\cdot n_{2}\in[-1,1]$ and a basis of functions on that interval are the Legendre polynomials $P_{l}$, so

[TABLE]

where

[TABLE]

is the observational 2-point spectrum and ${\cal O}_{lm}$ are the harmonic components of ${\cal O}$

[TABLE]

Note that we work with the unit-average normalization of the spherical harmonics (A.5), which is the natural one in this context, and that we keep the summation over $l,m$ indices implicit for notational simplicity. In the presence of both dummy and free $l,m$ indices, their nature will be inferable unambiguously by looking at both sides of the equation. The $m$ indices will always be clearly associated to some $l$ value and therefore run from $-l$ to $l$, while the $l$ indices run from $s$ to $\infty$, where $s$ is the spin of the observable under consideration.444For instance, we have $s=0$ for CMB temperature maps, while $s=2$ for the maps of the electric and magnetic components of the polarization field. For $N>2$, the decomposition is performed in detail in appendix A for generic functions $f(n_{1},\dots,n_{N})$ and the result is the following

[TABLE]

where

[TABLE]

and

[TABLE]

and where we have defined the following coefficients

[TABLE]

Although the latter appear as a sum over the $M_{k}$ indices, they are actually a single product because a Wigner $3-j$ symbol vanishes if the sum of its $m_{k}$ entries is not zero. Also for that reason, on can directly check that these coefficients are non-zero only if $\sum_{k=1}^{N}m_{k}=0$.

Equations (3.20) to (3.30) are the main result of this section. Note that they can actually also encompass the $N=1,2$ cases if we define the redundant notation

[TABLE]

[TABLE]

and

[TABLE]

The “harmonic” components of the $N$-point OCF (3.21) are the observational $N$-point spectra which generalize the structures (3.3) to arbitrary $N$. In particular, for the $N=3,4$ cases we recover the known results [42], up to a different normalization convention factor $(-1)^{L}\sqrt{2L+1}$ for $N=4$

[TABLE]

As for the functions (3.22), they form indeed a basis for SO(3)-invariant functions on $S_{2}^{N}$, as shown in appendix A. They are orthonormal (A.25) and therefore we can invert (3.20)

[TABLE]

Observe how the quantities (3.21) and (3.22) are explicitly SO(3)-invariant, since they only depend on total angular momentum numbers $l_{k}$ and $L_{k}$, making a total of $2N-3$ indices. As noticed earlier, this is indeed the number of independent angles present in the corresponding correlation functions. Note also that (3.21) and (3.22) are the very same combinations of the observable components and spherical harmonics, respectively. As we will see later on, all of the non-trivial contractions of indices will always be controlled by the coefficients (3.30) with no extra factors, which comes from the fact that they are orthonormal in the sense of (A.23). Thus, contrary to the implicit method described in the previous subsection 3.1, here there exists a normalization of the reduced spectra that appears as naturally privileged for computational convenience.

For the sake of completeness, let us also present the generalization of (3.4) by reminding that the TCFs (2.17) are SO(3)-invariant functions on $S_{2}^{N}$ (see (2.19)), so they can also be decomposed in the basis (3.22). Using (3.16) and the linearity of averaging we find

[TABLE]

Finally, to make contact with the implicit method described in subsection 3.1, let us also invert equation (3.49) to obtain the generalization of (3.1) and (3.2)

[TABLE]

This is the general solution to the equation of the implicit method (3.5) for arbitrary $N$, up to a $l_{k},L_{k}$-dependent normalization. This result highlights again the naturalness of the normalization used here for the reduced angular spectra, as all the complicated factors entering any of the above equations always combine to form the multilateral Wigner symbols exactly. Note also that, with our construction, the consideration of the full spectra (3.50) becomes superfluous, as one no longer requires them to derive the reduced ones as in the implicit method of the previous subsection.

3.3 Multilateral diagrammatic representation

For $N=3$ the coefficients defined in (3.30) simply reduce to the Wigner $3-j$ symbol

[TABLE]

This quantity contains a $\{l_{1}\,\,\,l_{2}\,\,\,l_{3}\}$ factor, for which an alternative definition from (A.22) is

[TABLE]

hence the common name “triangle delta”. Consequently, the general coefficient (3.30) contains a $\{l_{1}\dots l_{N}|L_{1}\dots L_{N-3}\}$ factor, defined in (A.24) as a product of certain $\{j_{i},j_{k},j_{m}\}$ factors. Generalizing the picture laid out in [42, 69] for the $N=4,5$ cases,555See also section 2.4.1 of [70] for the same picture, and for arbitrary $N$, but in the case of the flat sky limit where the sphere $S_{2}$ is replaced by $\mathbb{R}^{2}$. this quantity (A.24) can be interpreted as a “multilateral delta” in the following way. Being a product of $N-2$ triangle deltas, it is non-zero only if all of the involved integers $l_{1},\dots,l_{k}$ and $L_{1},\dots,L_{N-3}$ form their respective triangles. Moreover, with respect to the triangle ordering in (A.24), one of the edges of two neighboring triangles must have the same length $L_{k}$, so we can picture the triangles as being connected by common edges. The resulting shape is therefore a multilateral with $N$ edges of lengths $l_{1},\dots,l_{N}$, while the $L_{1},\dots,L_{N-3}$ integers correspond to the lengths of the $N-3$ diagonals connected to the vertex where the edges $l_{1}$ and $l_{N}$ join (see figure 1). Thus, the multilateral delta is not non-zero for any set of $l_{1},\dots,l_{N}$ that can form a multilateral, but only for those whose diagonals also have integer length. This suggests naming the symbols $W_{\dots}^{\dots}$ defined in (3.30) the “multilateral Wigner symbols”, with the $3-j$ symbols therefore corresponding to the “triangular” case.

From this multilateral picture one can also directly infer that the decomposition (3.20) is not unique for $N>3$, as was already recognized in [42] for the case $N=4$. Indeed, there is an ambiguity in choosing the vertex with respect to which the diagonals are drawn. This ambiguity, modulo symmetries of the $3-j$ symbols, corresponds to alternative orderings of the $(l_{k},m_{k})$ pairs and therefore leads us to consider permutations of the multilateral Wigner symbol indices. As shown in detail [42] for the $N=4$ case, these are achieved by taking linear combinations with Wigner $6-j$ symbols, thus allowing one to relate all possible orderings.

4 Observer terms and ensemble averaging

4.1 The issue

An important limitation in cosmology is that, observationally, we have access only to one single realization of the statistical process under consideration, the universe, and only from one single vantage point. One therefore requires some assumptions in order to relate the theoretical predictions, which are of statistical nature, to observations. In the standard approach, we consider stochastic small initial fluctuations in a homogeneous and isotropic universe, which evolve under gravity. We also assume that their probability distribution functional is homogeneous and isotropic and that averaging over positions and directions is equivalent to ensemble average, i.e. an ergodic hypothesis. With these assumptions one can effectively treat different parts of the universe as different realizations and relate the theoretical ensemble averages to observational light-cone averages.

The requirement for the ergodic hypothesis to be applicable is the observational availability of a large enough number of source points. Indeed, the common motto in the literature is that the theorist is allowed to ensemble average a given product of fields evaluated at source positions, because the observer probes several of these positions. This viewpoint is perfectly sufficient as long as one is only interested in the lowest-order approximation where a cosmological observable is entirely determined by the value of fields at the source position. In the era of precision cosmology, however, the above approximation is no longer valid, because observations are able to capture several sub-leading effects which depend on fields evaluated on the whole line of sight from the source to the observer (e.g. weak lensing). As one approaches the observer along its past light-cone, the number of observable points decreases and therefore, one would naively infer, so does the degree of applicability of the ergodic hypothesis. The extreme case is the observer point itself, for which we can only have a single measurement. Are we then allowed to perform ensemble averages of fields at, and in the vicinity of, the observer point, as we implicitly do in (2.17)? The aim of this section is to answer this question and its conclusion is in the affirmative.

Before demonstrating this statement, however, let us explain in more detail why this question is relevant for precision cosmology. To that end, note first that the typical expression for a cosmological observable ${\cal O}$ to linear order in perturbation theory is of the form

[TABLE]

where $X_{\rm o}$ and $X_{\rm s}$ denote field fluctuations evaluated at the observer and source positions, respectively, while $X_{\rm los}$ denotes an integral over fields evaluated along the background line of sight. The quantity that is relevant for comparing with observations is then, e.g., the corresponding 2-point correlation function, whose second-order contribution is

[TABLE]

Now note that the split (4.1) is not unique, as one can always transfer quantities between the three kinds of contributions by performing integrations by parts. Moreover, since we parametrize ${\cal O}$ with the observed redshift $z$ and position in the sky $n$, the function ${\cal O}^{(1)}(z,n)$ is gauge-invariant, because the full observable ${\cal O}(z,n)$ has no dependence on coordinates whatsoever [58]. The split (4.1) can then be chosen such that each of the three terms are separately gauge-invariant. In practice, however, the expression (4.1) often presents itself with a gauge-dependent observer term $X_{\rm o}$, which can therefore be set to zero with an appropriate choice of gauge, e.g. the synchronous gauge. Nevertheless, the expressions found in the literature are often in longitudinal gauge, in which case $X_{\rm o}\neq 0$.

As already mentioned, one can consistently avoid the presence of observer terms (gauge-invariant or not) by an appropriate integration by parts, but one then needs to take into account the resulting extra terms in $X_{\rm los}$ and $X_{\rm s}$ that this manipulation introduces. Instead, the standard approach in the literature is to work in the longitudinal gauge and simply discard $X_{\rm o}$ as soon as it appears in the computation. The issue of whether one can ensemble average field values at the observer position therefore never arises in the standard practice. As for the line-of-sight terms $X_{\rm los}$, they are ensemble averaged just like the source terms $X_{\rm s}$, which implicitly assumes that ergodicity applies to all points on the light-cone, independently of how close they are to the observer point.

It is now well understood that discarding the observer terms $X_{\rm o}$ as described above leads to several qualitative problems:

• •

In most cases the resulting ${\cal O}^{(1)}$ is not gauge-invariant under gauge transformations at the observer because one discards a gauge-dependent term $X_{\rm o}$. However, the gauge transformation of the observable only produces observer terms, by construction, which is therefore “ok” in this approach, since one is precisely neglecting such terms to begin with.

• •

Following the construction of section 2.2, Fourier decomposing the field configurations $\phi(t_{0})$ and using statistical homogeneity to eliminate one of the ${\rm d}^{3}k$ integrals, one obtains

[TABLE]

where the integrand is a combination of field $k$-spectra at $t_{0}$, multiplied by growth functions and $e^{i\vec{k}\cdot\vec{x}}$ factors evaluated on the light-cone. This integrand is generically divergent in the infrared, a fact which has led to confusion in the literature [71]. Moreover, the dependence on the infrared cut-off that has to be introduced to control the divergence in (4.3) breaks the equivalence principle, as the observable correlation function ends up depending on the uniform gravity mode. Taking into account $X_{\rm o}$ precisely cancels the divergent terms, thus showing that they cannot be associated with physical effects and that the equivalence principle holds as one would expect [56, 60, 63].666For a detailed demonstration and discussion of this cancellation in the 2-point case, we refer the reader to [56] for the luminosity distance and to [63] for the galaxy number density observables. As discussed at the beginning of section 2, the integrand in (4.3) is not itself an observable quantity, but it is interesting to see that the consideration of observer terms preserves the physical intuition in $k$-space.

• •

The sky parametrization $n$ does not generically correspond to the one of an actual observer.777See [64] for a detailed explanation. In particular, note that while the introduction of the Sachs basis allows one to locally decompose tensors on the sky in the basis the observer uses, it does not provide the observer parametrization of the sky itself that is required in order to compute correlation functions and spectra as the observer does. Indeed, the terms $X_{\rm o}$ contain corrections that implement the diffeomorphism relating the observer parametrization $n$ to the one induced by the local coordinate system around the observer and therefore to the gauge that is chosen. In longitudinal gauge the coordinate-induced frame in the observer’s tangent space is not orthonormal in the presence of vector and tensor modes, as the one of an observer should be, so neglecting $X_{\rm o}$ leads to spurious contributions to the “lensing” of observables in the sky.

Nevertheless, it turns out that there is no significant quantitative problem today in neglecting $X_{\rm o}$ to linear order in perturbation theory. This is because the actual observer approximates the ensemble average in (4.2) with the SO(3) averages of the data in the sky, denoted by $\langle\dots\rangle_{\rm SO(3)}$ in what follows. $X_{\rm o}$ is special with respect to this operation, because it has no angular dependence (if its spin vanishes) and can therefore be factored out. Consequently, the contribution $\langle X_{\rm o}^{2}\rangle_{\rm SO(3)}\equiv X_{\rm o}^{2}$ is a monopole term (for spin $s=0$) in the corresponding power spectrum, while the cross-terms are simply zero, e.g.

[TABLE]

If ${\cal O}$ is an observable of non-zero spin $s$ (e.g. for the observer velocity, where $s=1$, or the CMB polarization, where $s=2$), then the $\langle X_{\rm o}^{2}\rangle_{\rm SO(3)}$ will contribute to the first non-trivial multipole which is the $s$ multipole. But the overwhelming part of the relevant information actually lies in the rest of the spectrum, so ignoring the observer terms of cosmological observables is indeed irrelevant quantitatively.888Note that for some observables, such as galaxy number density and luminosity distance, the standard analysis is performed with correlation functions instead of spectra, in which case all multipoles contribute to the result and the lower-order ones are not subtracted. In that case, it has been shown that observer terms can be of the same order as other relativistic corrections [57, 60].

The caveat of the previous argumentation is that this is only true to linear order in perturbation theory. Indeed, to second order the solution of the observable will have cross-terms of the schematic form999For a concrete example, see equations (4.2) to (4.5) of [40] for the second-order part of the redshift perturbation.

[TABLE]

implying terms of the form

[TABLE]

for the next order of the 2-point correlation function. Importantly, these contributions will affect all multipoles of the third-order power spectrum with an amplitude that is a priori comparable to the standard terms such as $\sim\langle X^{2}_{\rm s}X_{\rm los}\rangle$, etc. Moreover, similar terms will also appear in higher-order correlation functions. Thus, in the era of precision cosmology one can no longer ignore the observer terms $X_{\rm o}$, in the usual gauge choices that are considered, meaning that it is relevant to discuss their inclusion inside ensemble averages.

4.2 Relating theory and observation

We now wish to derive the relation between the observational and theoretical $N$-point functions of observables (3.12) and (2.17). We already have equation (3.16), but this is a theoretical relation, as it involves ensemble averaging, and simply states that $G^{\rm ob}$ is a full unbiased estimator of $G^{\rm th}$. Rather, we are interested in relating the two $N$-point functions as one does in practice, i.e. we need $G^{\rm ob}$ to be evaluated on a definite field configuration, and this is achieved using the ergodic hypothesis. In particular, we want to see how the latter enters the derivation and how it affects the issue of the single observer point.

We start by recalling subsection 2.2 and consider the set of 3-dimensional field configurations $\phi_{\mathfrak{a}}(t_{0})$ over which we sum when performing the ensemble average in (2.13). Given the invariance of $P[\phi(t_{0})]$ under the isometry group, it is convenient to partition this ensemble into equivalence classes, where two field configurations are deemed equivalent if they can be related by an isometry. For the sake of simplicity, let us consider here the case of flat background space, so that the isometries form the Euclidean group $\mathbb{E}:={\rm SO(3)}\ltimes\mathbb{R}^{3}$, and let us also use the trivial Cartesian coordinates. Any element of a given equivalence class can be described as an isometry of some fixed representative $\hat{\phi}_{\mathfrak{a}}(t_{0})$

[TABLE]

where $R$ is a rotation matrix, $\vec{c}$ a translation vector and $M_{\mathfrak{a}}^{\,\,\mathfrak{b}}$ is the matrix that rotates tensor indices in $\mathfrak{a}$. We can therefore split the functional integration in (2.13) into an integral over the elements of a given class followed by an integral over all possible classes. By the latter we mean an integral over suitably chosen representatives $\hat{\phi}_{\mathfrak{a}}(t_{0})$ such that the corresponding functional integral is well-defined.101010The existence of such a splitting of the integration is a non-trivial mathematical assertion, whose proof, if possible, would go beyond the scope of this paper. Since the pdf is constant over all representatives of a given class, the integral in (2.13) becomes

[TABLE]

where we now integrate only over the set of representatives, and

[TABLE]

is the average over the Euclidean group action over the field configurations, $V$ is the total volume and ${\rm d}R$ is the SO(3) Haar measure. As one may expect, the ensemble average therefore contains a purely geometric average over the symmetry group of the pdf.

Apart from subsets of measure zero, the configurations $\hat{\phi}(t_{0})$ appearing in the integral (4.8) have a rich spatial dependence. In particular, they are non-periodic functions which therefore probe a large variety of local field profiles for large enough $V$. In the $V\to\infty$ limit, which we consider here, the ergodic hypothesis states that this probing is thorough enough to make the Euclidean average in (4.8) independent of the configuration $\hat{\phi}(t_{0})$. As a result, that average factorizes out of the integral, thus yielding the relation

[TABLE]

Now on the right-hand side it is a definite, although generic field configuration that is considered. We can then generalize this manipulation straightforwardly to the case of the FCFs (2.14), since the Euclidean group action is independent of the time variable $t$, and therefore also to the case of the TCF (2.17)

[TABLE]

where now we act with $\mathbb{E}$ directly on the 4-dimensional fields $\phi_{\mathfrak{a}}(x)$ instead of the $\phi_{\mathfrak{a}}(t_{0},\vec{x})$. Note also that we have not specified the $\bar{\Phi}$ dependence in (4.2) for simplicity. The above expression allows us to make contact with the OCFs (3.12). Let us first redefine the dummy variable $\vec{c}\to\vec{c}-R\vec{x}_{o}$ in (4.2) and let us also define the notation for shifted fields

[TABLE]

to get

[TABLE]

We next observe that, for a given value of $\vec{c}$, the fields $\phi_{\mathfrak{a},\vec{c}}$ are rotated by $R^{-1}$ around $\vec{x}_{o}$ in (4.13). But rotating the fields around $\vec{x}_{o}$ is tantamount to rotating the observed angles $n$ in the opposite direction, so the SO(3) average on the fields translates into an OCF

[TABLE]

Similarly, translating the fields by $\vec{x}_{o}-\vec{c}$ can be equivalently expressed as shifting the observer by $\vec{c}-\vec{x}_{o}$, so we finally obtain (after renaming $\vec{c}\to\vec{x}_{o}$)

[TABLE]

i.e. the average over the action of the Euclidean group on the fields $\phi_{\mathfrak{a}}$ is equivalent to averaging over all observer reference frames, as in $G^{\rm ob}$, but also over all observer positions. This is simply because we have imposed both statistical isotropy and homogeneity. For instance, the CMB maps one would obtain from another vantage point of the universe $\vec{x}^{\prime}_{o}$ would be different from the ones we observe on earth today, while the theoretical power spectrum which we calculate is the same for all vantage points. Since the relation between the correlation function and their associated spectra is linear, we have that (4.15) simply translates into

[TABLE]

Equations (4.15) and (4.16) are the central result of this section. The important thing to understand from them is that the application of the ergodic hypothesis singles out a definite field configuration on the right-hand side, but it does not single out the position of the actual observer, since all possible such positions are taken into account equally inside an average. Therefore, the applicability of the ergodic hypothesis depends on the representative field configuration which we choose on the right-hand side of (4.15), completely independently of how this configuration is probed observationally. In particular, as stated in the step in which ergodicity is actually used (4.10), what we need is that field configuration to have a rich profile, which is a decent assumption to make about the configuration of the actual universe.

What we also see, however, is that the use of ergodicity is necessary, but not sufficient to relate theory and observations yet, because the information from several $\vec{x}_{o}$ that is required in (4.15) and (4.16) is not observationally available. Rather, the step which actually singles out the observer point is the approximation of the $\vec{x}_{o}$ averages in the right hand-sides of (4.15) and (4.16) by the value at the actual observer position

[TABLE]

and

[TABLE]

respectively. Approximating an average over some set by a single value within that set is of course the crudest possible estimator of that average. One expects this estimation to be accurate enough only if the values of the considered set exhibit a relatively small dispersion around the average value. But this dispersion in $\vec{x}_{o}$ space cannot be measured for the very same reason that brought us to this approximation in the first place, i.e. the observational unavailability of data at other observation points. Nevertheless, as we will see in the next subsection, the variance of $G^{\rm ob}$ with respect to the $\vec{x}_{o}$ dependence is equal to its theoretical analogue, the statistical or “cosmic” variance, through the ergodic hypothesis. Thus, the error one makes in (4.17) or (4.18) is precisely cosmic variance, which is indeed negligible on small enough scales (see section 5.2).

Equations (4.17) and (4.18) allow to relate theory and observation, as both sides can be computed in their respective domains. Finally, another way of interpreting the result (4.15) is by combining it with (3.16) to obtain

[TABLE]

i.e. the ergodic hypothesis and statistical homogeneity equate the ensemble averaging with the translational averaging of a generic configuration. Thus, under these assumptions, while estimating $G^{\rm th}$ with $G^{\rm ob}$ first appears as a double approximation, i.e. the measurement of a single universe realization and from a single viewpoint, these two turn out to be the very same approximation.

4.3 Unambiguous covariance matrix and cosmic variance

Equation (4.19) shows that, assuming ergodicity and statistical homogeneity, ensemble averaging $G^{\rm ob}$ is tantamount to averaging it over $\vec{x}_{o}$ with a generic configuration of the fields. The last open question is therefore whether these two types of averaging are also equal when it comes to the variance associated with the estimation (4.17). The statistical covariance matrix is

[TABLE]

where $\alpha_{N}$ collectively denotes the set $\{z_{k},n_{k}\}_{k=1}^{N}$ for notational simplicity, while the spatial one is

[TABLE]

The former leads to the notion of cosmic variance, so, if both of them are the same, we would have shown that our indiscriminate use of ensemble averaging leads to no new uncertainties between theory and observation. This is simply achieved by noting that a product of two OCFs can be expressed as a partial SO(3) average of a single OCF

[TABLE]

Inserting this in (4.20) and (4.21) and using (4.15) and (3.16) we find that both covariances are equal

[TABLE]

and that

[TABLE]

Thus, the absolute “1-sigma” error of the estimator (4.17)

[TABLE]

is unambiguously cosmic variance. We see that our derivation provides a refined understanding of that notion, which is not usually expressed in cosmology textbooks. Given our statistical assumptions, cosmic variance is the error due to the fact that we observe a single realization of the universe and from a single vantage point $\vec{x}_{o}$. In accordance with the discussion of the previous subsection, if either of these two conditions were dropped, then there would be no cosmic variance. On one hand, if we had simultaneous access to the data of a single realization from all possible $\vec{x}_{o}$, then the ergodic hypothesis (4.15) would allow us to match the theoretical predictions exactly.111111Here we neglect the fact that this information would also require a time $\sim V^{1/3}$ to be collected by a main observer and therefore analyzed. On the other hand, if we could observe all possible universe realizations, even from a single viewpoint $\vec{x}_{o}$, then we would be able to compute directly the theoretical $N$-point functions, which are independent of $\vec{x}_{o}$. In the first case we are technically setting $\Sigma_{\rm spat.}\to 0$, whereas in the second one it is rather $\Sigma_{\rm stat.}\to 0$.

In the next section we will see that, although cosmic variance is always non-zero in practice, the associated relative uncertainty on the spectra must tend to zero when $l_{k}\to\infty$ as a consequence of statistical isotropy and the central limit theorem. Technically, this corresponds to $\Sigma_{\rm stat.}$ becoming negligible compared to $G^{\rm th}$ for large enough $l_{k}$ values, with statistical homogeneity and the ergodic hypothesis then implying that so does $\Sigma_{\rm spat.}$, as we just saw in (4.23). This is the reason why, in practice, one can still obtain statistical cosmological information even from a single observational point.

5 Covariance in $l$-space

5.1 The covariance matrix of angular $N$-point spectra

From now on we drop the dependence on $t_{o}$ for notational simplicity. Here we wish to derive the expression for the covariance matrix of the reduced $N$-point and $N^{\prime}$-point spectra, i.e. the harmonic analogue of (4.24). To that end, we first note that ${\rm Cov}\left(t_{o};\alpha_{N};\alpha^{\prime}_{N^{\prime}}\right)$ is a function on $S_{2}^{N}\times S_{2}^{N^{\prime}}$ that is invariant under independent rotations on each of these two components, so we can decompose each of the corresponding arguments in the basis (3.22). The harmonic components are then found by projecting twice

[TABLE]

To express this in terms of the reduced theoretical spectra, we use (4.24) and decompose the TCFs in the basis (3.22). By proceeding exactly as in appendix A to eliminate the SO(3) average in (4.24), and using the orthonormality relations (A.5) and (A.23), we find

[TABLE]

Now the second product of Wigner symbols can be simplified using (B.5) iteratively, then (B.6) and finally (B.1) to find

[TABLE]

We therefore obtain

[TABLE]

where we have defined the following contraction of multilateral Wigner symbols of equal order

[TABLE]

which is rotationally-invariant, as it does not depend on $m_{k}$ indices. From (5.4) we see that the covariance matrix of the $N$-point and $N^{\prime}$-point spectra is controlled by the squeezed $(N+N^{\prime})$-point spectrum, as depicted in figure 2, which arises from gluing together an $N$-point and an $N^{\prime}$-point diagram. Therefore, the squeezed $(N+N^{\prime})$-point spectrum here is already proportional to $\{l^{\prime}_{1}\dots l^{\prime}_{N^{\prime}}|L^{\prime}_{1}\dots L^{\prime}_{N^{\prime}-3}\}$, which is why we were able to discard that factor coming from (5.1).

The coefficients defined in (5.5) form a symmetric square matrix in the space of reduced $N$-point spectra and in (5.4) we see that this matrix acts on the $N$-point component of the squeezed $(N+N^{\prime})$-point spectrum. Comparing with the analogous expression for the correlation functions (4.24), we thus see that the $A$ matrix (5.5) implements the operation of partial SO(3) average on invariant functions in harmonic space. Finally, note that one can extend the expression (5.4) to the cases $N,N^{\prime}=1,2$ by using a redundant notation analogous to Eqs. (3.31), (3.32) and (3.38).

5.2 Cosmic variance in the large $l_{k}$ limit

Let us consider the case $N=N^{\prime}$, where the covariance matrix quantifies the typical error of the $G^{\rm th}\approx G^{\rm ob}(\vec{x}_{o})$ approximation that is cosmic variance, so that the $l$-space analogue of (4.25) reads

[TABLE]

We now note that this quantity has a simple asymptotic behavior at large $l_{k}$, which is completely determined by the central limit theorem [72]. Indeed, thanks to the assumption of statistical isotropy, the ${\cal O}_{lm}(\vec{x}_{o})$ for a given value of $l$ are $2l+1$ independent and equally distributed random variables. Thus, the product $\prod_{k=1}^{N}{\cal O}_{l_{k}m_{k}}(\vec{x}_{o})$ for given $l_{k}$ values consists of ${\cal O}\left(\prod_{k=1}^{N}l_{k}\right)$ independent and equally distributed random variables. The observational spectrum component $G^{\rm ob}_{l_{1}\cdots l_{N}|L_{1}\cdots L_{N-3}}(\vec{x}_{o})$ then corresponds to a weighted sum of these products over their $m_{k}$ indices (3.21). However, since the multilateral Wigner symbol (3.30) is zero unless $\sum_{k=1}^{N}m_{k}=0$, the sum actually only contains ${\cal O}\left(\prod_{k=1}^{N-1}l_{k}\right)$ independent and equally distributed random variables. Finally, noting that $G^{\rm ob}(\vec{x}_{o})$ is technically a sample average and $G^{\rm th}$ is its ensemble average, the central limit theorem [72] states that

[TABLE]

Consequently, the relative cosmic variance tends to zero as $l_{k}\to\infty$. In particular, this means that $G^{\rm ob}_{l_{1}\cdots l_{N}|L_{1}\cdots L_{N-3}}(\vec{x}_{o})$ tends to the $\vec{x}_{o}$-independent result $G^{\rm th}_{l_{1}\cdots l_{N}|L_{1}\cdots L_{N-3}}$, without the need to perform the average over $\vec{x}_{o}$ as in (4.16), thus justifying the approximation (4.17) for small enough scales. This is therefore how one can obtain accurate statistical cosmological information from a single vantage point for sufficiently large $l_{k}$. Following the discussion at the end of section 4, this loss of the $\vec{x}_{o}$ information in the $l_{k}\to\infty$ limit is a consequence of statistical homogeneity and the ergodic hypothesis.

6 Connected angular $N$-point spectra

6.1 Definitions, modified ergodic relation and covariance matrices

Here we focus on the part that is usually of most interest: the connected component of a given correlation function or associated spectrum. This is a non-linear combination of the full statistics $G$, so one should expect the linear ergodic relations (4.15) and (4.16) to be modified. Here we will illustrate this modification by considering the case of the connected 2,3, and 4-point functions. We start by defining the observable fluctuations, in the observational and theoretical cases, as the deviation of the observable from the respective 1-point functions

[TABLE]

or, in terms of the harmonic components,

[TABLE]

We note in particular that, by construction, the observational monopole is identically zero, while the theoretical one captures precisely the difference between the two 1-point functions

[TABLE]

The connected 2,3,4-point functions can then be defined by

[TABLE]

where here the star is a placeholder for “ob” and “th”, the corresponding averages are respectively $\langle\dots\rangle_{\rm SO(3)}$ and $\langle\dots\rangle$ and we have omitted the $\vec{x}_{o}$ dependencies in the “ob” case for notational simplicity. Clearly, the these functions are invariant under a common rotation of their arguments, so we can compute their harmonic components by projecting on the basis (3.22). Using the connected analogue of (3.17) and also (A.15), (B.5), (B.6) and the $3-j$ symbol symmetries, we recover the known expressions for the reduced power spectrum, bispecrum and trispectrum estimators [42]121212Note that (6.19) agrees with Eqs. (26) and (20) of [42], as the normalization convention is different by an overall factor of $(-1)^{L}\sqrt{2L+1}$.

[TABLE]

The theoretical counterparts are then obtained by replacing “ob” $\to$ “th” and ensemble averaging, again in accordance with the $N=2,3,4$ results of [42].

Let us now see how the linear ergodic relation (4.16) looks like in terms of the connected parts. First, because of (6.5), the observational spectra vanish identically whenever at least one of their $l_{k}$ entries is zero, i.e. they have vanishing “monopoles” by construction. In contrast, this is not the case for the theoretical spectra (6.6). As a concrete example, let us consider the relation between the connected 2-point correlation functions, by expressing them in terms of the full functions and using (4.15)

[TABLE]

where ${\rm Cov}\left(z_{1};z_{2}\right)$ is the covariance matrix of the 1-point function with itself (see (4.21)). Alternatively, in harmonic space

[TABLE]

The difference with (4.16) in the connected case is a monopole term, which is also exactly the multipole for which the observational spectrum vanishes identically $C^{\rm ob}_{0}(z_{1},z_{2})\equiv 0$. Therefore, a better display of (6.21) could be

[TABLE]

Thus, the ergodic relation (4.16) still holds for the components containing physical information $C^{\rm ob}_{l}\neq 0$, while the left-over equation relates the covariance matrix of the 1-point function to the monopole of the theoretical power spectrum. As a result, the absolute 1-sigma error (4.25) associated with the 1-point function approximation

[TABLE]

is simply

[TABLE]

This picture generalizes to the case of higher $N$, as one can check by expressing $\Delta^{\rm th}{\cal O}_{lm}$ in terms of $\Delta^{\rm ob}{\cal O}_{lm}$ in the theoretical spectra $B^{\rm th}$ and $T^{\rm th}$. One finds that the resulting equation splits into two sets of equations: the non-monopole ones, i.e. the ones where none of the $l_{k}$ is zero and for which (4.16) still holds, and the ones relating the theoretical monopoles to the covariance matrix and higher order analogues (skewness, kurtosis, etc.) of lower-$N$ spectra. For this reason, from now on all of our equations will hold up to monopole terms, so that we do not have to deal with this subtlety.

The next non-trivial case is the trispectrum, where one now subtracts products of 2-point functions. For the purposes of our point it suffices to describe it schematically as

[TABLE]

where the ellipses here denote the other two possible orderings of the 2-point function arguments. In complete analogy with the derivation in the $N=2$ case (6.20), using (4.15) we find

[TABLE]

so the difference is again a covariance matrix, but now the one of the 2-point function with itself. In harmonic space we then find

[TABLE]

where we recognize the same structure as in (6.19) by construction and we have the covariance matrix of the 2-point spectra

[TABLE]

A first important difference with the 2-point case (6.21) is that the extra terms here affect multipoles of arbitrary magnitude, although only a small subset, the one of pairwise equal $l_{k}$’s. A second important difference is that the components of $T^{\rm ob}$ with pairwise equal $l_{k}$’s are not identically zero in general, i.e. the equations no longer split as in the monopole case. One must therefore take (6.1) as it is and infer that the analogue of the single-observer approximation (4.18) in this case is actually

[TABLE]

Finally, just as the 2-point case (6.21) with $l=0$ provides the cosmic variance of the 1-point function for free, the 4-point case (6.1) with $l_{1}=l_{2}=:l$, $l_{3}=l_{4}=:l^{\prime}$ and $L=0$ provides the cosmic variance for the power spectrum. However, here $T^{\rm ob}_{lll^{\prime}l^{\prime}|0}\neq 0$, but rather (use (6.19) and properties of the $3-j$ symbols)

[TABLE]

so, inserting this in (6.1), using (6.1) and isolating the covariance matrix, we find

[TABLE]

This generalizes the result of [42], worked out for the CMB where there is no redshift dependence, to the case of generic redshifts along the light-cone. In particular, the cosmic variance of the power spectrum, i.e. the absolute 1-sigma error (4.25) in the approximation

[TABLE]

is

[TABLE]

Finally, for Gaussian statistics and equal redshifts, we recover the well-known result131313See e.g. [73] for the case of the CMB power spectrum.

[TABLE]

Given the discussion of subsection 5.2 and comparing (6.1) with (5.7), we find that here the central limit theorem essentially states that the trispectrum decays faster with growing $l$ than the $C_{l}^{2}$ terms.

6.2 Mind the monopole when using relative fluctuations

Until now we have considered the “absolute” fluctuations $\Delta^{\rm ob}{\cal O}$ and $\Delta^{\rm th}{\cal O}$, but in practice the most convenient ones to work with are the relative ones (if $G^{\rm th}(z)\neq 0$)

[TABLE]

This introduces a non-linear difference between the observational and theoretical definitions and will affect the covariance matrices of the corresponding spectra. Indeed, consider for instance the case of the relative power spectra

[TABLE]

which are denoted using tilded letters to distinguish them from the absolute ones $C_{l}^{\rm ob}$ and $C_{l}^{\rm th}$. The corresponding covariance matrix is given by

[TABLE]

and one would naively expect this to be simply

[TABLE]

where here ${\rm Cov}_{l;l^{\prime}}$ is the covariance matrix of the absolute 2-point spectrum (6.31). In particular, one would then infer the analogue of (6.34) for the corresponding cosmic variance in the equal redshift Gaussian case

[TABLE]

However, as we will now show, Eqs. (6.39) and (6.40) are actually only approximate, because in order to obtain them one must wrongly assume that $G^{\rm ob}(\vec{x}_{o};z)=G^{\rm th}(z)$ or, according to (6.5), that the monopole of the observable is zero at $\vec{x}_{o}$ for all $z$

[TABLE]

Indeed,

[TABLE]

where the remainder term $\tilde{R}_{ll^{\prime}}$ depends on the monopole $\delta^{\rm th}{\cal O}_{00}(\vec{x}_{o};z_{k})$ that comes from converting the $G^{\rm ob}(\vec{x}_{o},z_{k})$ into $G^{\rm ob}(z_{k})$ and which cannot go through the integral over $\vec{x}_{o}$, because of its dependence on that variable. The remainder term is therefore controlled by a combination of monopoles of the theoretical power spectrum, bispectrum and trispectrum which, unlike the extra “mon.” terms described in the previous section, contribute here to all $l$ values, not just to $l=0$. Thus, it is important to know whether one compares absolute quantities or relative ones, because the corresponding covariance matrices are different and therefore so will be the results of the corresponding Fisher forecasts.

7 Conclusion

In this paper we have addressed two potentially important issues for the next stages of precision cosmology. First, we have developed the general theory of the rotationally-invariant (“reduced”) angular $N$-point spectra. The method employed so far [42] is not suited for obtaining expressions for arbitrary $N$, because it requires solving some increasingly complicated algebraic equation as $N$ grows. As a result, one must work out each $N$ case separately and the present literature contains a detailed description only up to the $N=4$ case [42]. We have presented an alternative construction which provides all the relevant definitions and relations associated with reduced angular $N$-point spectra straightforwardly and for all $N$. This includes the covariance matrix of these spectra for arbitrary $N$ and $N^{\prime}$, which is controlled by a squeezed $(N+N^{\prime})$-point spectrum and leads to a cosmic variance of the $N$-point spectrum that decays as $l^{(N-1)/2}$ for large $l_{k}$ numbers. Our construction is based on the introduction of the “multilateral” Wigner symbols, which generalize the Wigner $3-j$ symbols of triangles to polygons. With these we have built an orthonormal harmonic basis for $N$-point correlation functions on the sphere, such that the reduced spectra appear as the coefficients in this basis. We have also discussed a corresponding diagrammatic representation, generalizing the one of [42, 69] for the cases $N=4,5$, which was also considered in [70] in the flat sky limit. Even though the determination of a generic $N$-point spectrum might be prohibitive numerically, our general framework allows easily for specialization, e.g. to squeezed spectra with only one $l$ varying while the others are fixed to some low value.

The second important part of the paper consists in motivating and justifying the consideration of observer terms of cosmological observables inside the ensemble averages of correlation functions and spectra. Through a careful derivation of the relation between theoretical predictions (“$G^{\rm th}$”) and observations (“$G^{\rm ob}$”), under the assumptions of statistical homogeneity and isotropy and ergodicity, we have shown that no special treatment of the observer point is required whatsoever on the theoretical side. Ensemble averaging field products at this point does not introduce any inconsistencies, nor does it imply any new assumptions or uncertainties. This motivates us to consider the full analytical expressions of cosmological observables in calculations, especially at non-linear order in perturbations where observer terms are intertwined with source and line-of-sight contributions in a non-trivial way, as we argued in subsection 4.1. Thus, with this conceptual ambiguity being lifted, we can now safely state that a rigorous treatment of cosmological observables should include the observer terms, or at least a check of the magnitude of their contribution in the final result if they are neglected.

We have also taken advantage of the present framework to discuss some subtle aspects of working with the connected part of spectra, as one usually does, when relating theory to observation. When working directly with fluctuations of observables around their average (1-point function) value, as is often the case in cosmology, the difference between full and connected spectra is simply a sum of products of spectra with $N>1$. However, here we pointed out that the notion of fluctuation is different in the theoretical and observational cases, since the 1-point functions are different (ensemble versus sky average). As a result, the relation between theoretical and observational connected spectra involves extra “monopole” terms. Moreover, if one works with relative, instead of absolute, fluctuations of observables, then one must take into account additional corrections which now affect all multipoles, not just the monopoles. A direct consequence of this fact is that the covariance matrices computed with the relative and absolute fluctuations are not proportional to each other, but have an additional remainder term. Thus, this can lead to different results in the corresponding Fisher analyses, depending on whether one uses absolute or relative observational fluctuations.

Finally, and remarkably, the fact that the observational and theoretical 1-point functions are not equal in general might also be relevant in the context of parameter estimation involving CMB data. Indeed, the standard treatment of the CMB temperature is to simply neglect cosmic variance and therefore consider these two quantities as equal $\bar{T}:=\left\langle T\right\rangle=\left\langle T\right\rangle_{\rm SO(3)}(\vec{x}_{o})$. This is justified by the high precision of the average CMB temperature measurement [74] and, most importantly, by the smallness of the associated relative cosmic variance $C_{0}\sim 10^{-5}$, which sets the fundamental uncertainty floor. In principle, however, one should consider $\bar{T}$ (or $\omega_{\gamma}$) as an extra cosmological parameter to be varied in the likelihood analysis. Through a careful Fisher analysis that incorporates the above-mentioned subtleties [75], the impact of an independent $\bar{T}$ on the other cosmological parameters and their errors is estimated in [76].

Acknowledgments

We are grateful to Fabien Lacasa and Thiago Pereira for pointing out references [70] and [69], respectively, and for useful comments. We acknowledge support by the Swiss National Science Foundation. J.Y. and E.M. are further supported by a Consolidator Grant of the European Research Council (ERC-2015-CoG grant 680886).

Appendix A Harmonic decomposition of ${\rm SO(3)}$-invariant functions on $S_{2}^{N}$

Consider a function $f$ on $S_{2}^{N}$ that satisfies

[TABLE]

for arbitrary rotations $R$. We can then average over the SO(3) group on the right-hand side (see (3.2))

[TABLE]

We then decompose each one of the $n_{k}$ dependencies into spherical harmonics

[TABLE]

where

[TABLE]

The summation over $l,m$ indices will be kept implicit for notational simplicity. In what follows, we will encounter both dummy and free $l,m$ indices, so their nature will be inferable by looking at both sides of the equation. The $m$ indices will always be clearly associated to some $l$ value and therefore run from $-l$ to $l$, while the $l$ indices run from [math] to $\infty$. Also, note that we use the less conventional normalization of spherical harmonics

[TABLE]

Equation (A.2) now reads

[TABLE]

where the Haar measure ${\rm d}R$ is given in (3.10). We can then extract the $R$-dependence out of the spherical harmonics by using their transformation property under rotations

[TABLE]

where the $D_{l}$ are the Wigner matrices forming the $(2l+1)$-dimensional irreducible representation of SO(3)

[TABLE]

We thus have

[TABLE]

where

[TABLE]

In the $N=2$ case, we can use the identity

[TABLE]

and the orthonormality relation

[TABLE]

to obtain

[TABLE]

and thus

[TABLE]

where the $P_{l}$ are the Legendre polynomials, arising from the identity

[TABLE]

and

[TABLE]

In the case of observable products, we obtain the well-known results (3.17) and (3.18). For the case of most interest $N>2$, we need to use iteratively the “Clebsch-Gordan” composition rule

[TABLE]

in order to reduce the product in (A.10) down to a single pair, in which case we can use again (A.11) and (A.12). The result is

[TABLE]

where the coefficients on the right-hand side are given in (3.30). Going back to (A.9), this finally implies the decomposition

[TABLE]

where

[TABLE]

and

[TABLE]

The latter form a basis for SO(3)-invariant functions on $S_{2}^{N}$. Indeed, we have just shown that they generate that space, so we must still show that they are orthogonal under the natural scalar product. Using the identity (B.5), where

[TABLE]

one obtains

[TABLE]

where we have defined

[TABLE]

and therefore the orthonormality relation

[TABLE]

Finally, one can use this to invert (A.19)

[TABLE]

Appendix B Some identities for Wigner $3-j$ symbols

Here we display all the identities satisfied by the Wigner $3-j$ symbols that are used in this paper.

• •

Special cases

[TABLE]

• •

Column exchange symmetry:

even permutations

[TABLE]

odd permutations

[TABLE]

• •

$m$-sign flip symmetry

[TABLE]

• •

Sum identities

[TABLE]

[TABLE]