© F. Lacasa 2020

1. Introduction

With the increase of galaxy density in current and forthcoming cosmic surveys, our statistical analysis of the large-scale structure of the Universe needs to be pushed towards new degrees of precision. Accurate covariance matrices are an important part of this effort. Indeed, using an incorrect covariance basically amounts to analysing a biased data set (e.g. Sellentin & Starck 2019). This effect has, indeed, been seen in current weak lensing surveys, along with changes in the covariance shifting cosmological constraints on S₈ = σ₈(Ω_m/0.3)^0.5 (Hildebrandt et al. 2017; Troxel et al. 2018), which is of particular importance in the current context of possible tensions between low- and high-redshift measurements of σ₈.

In the past, covariance matrices for the large scale structure were often estimated using jackknife or bootstrap techniques, however, this has been shown to be biased at a level inadequate for cosmological analyses (Norberg et al. 2009; Lacasa & Kunz 2017). Other analyses have used matter covariances coming from ensemble of simulations (Harnois-Déraps & Pen 2013), which correctly capture in-survey covariance, projected in 2D using the flat sky approximation (Sato & Nishimichi 2013). These matter covariances cannot be directly applied, however, to galaxy clustering, as Abramo et al. (2015) highlighted in stating that the amount of non-linearity is notably dependent on the galaxy selection; the scaling with bias of a 1-halo polyspectrum differs greatly from the perturbative polyspectrum. Most current analyses use covariances that either come from analytical computations or from dedicated simulations. In particular, analytical modeling using the halo model has risen to a state of the art level for several current galaxy surveys (Krause & Eifler 2017; Hildebrandt et al. 2017; Krause et al. 2017). This is the approach followed in this article, applying it to a galaxy clustering analysis using the angular power spectrum. I emphasise that the analysis and conclusions can be transferred to a real-space analysis using the two-point correlation function since there is a linear mapping between real space and harmonic space, where computations are much simpler to carry out.

The point of halo model covariances is to move beyond the vanilla analytical Gaussian formula. In the case of angular auto-spectra for galaxy clustering in disjoint redshift bins labeled i_z, j_z, this formula gives

$$\begin{array}{c}\hfill {\mathcal{C}}_{\ell ,{\ell }^{\prime }}^G=\frac{2C_{\ell }^{\mathrm{gal}}{(i_z)}^2}{2\ell +1}{\delta }_{\ell ,{\ell }^{\prime }}{\delta }_{i_z,j_z},\end{array}$$(1)

where throughout the article I use the short notation

$$\begin{array}{c}\hfill {\mathcal{C}}_{\ell ,{\ell }^{\prime }}\equiv \mathrm{Cov}(C_{\ell }^{\mathrm{gal}}(i_z),C_{{\ell }^{\prime }}^{\mathrm{gal}}(j_z)).\end{array}$$(2)

Using the halo model for non-Gaussian covariance terms allows not only for an adequate reproduction of super-sample covariance (SSC, Takada & Hu 2013) for power spectra, but it further allows for the inclusion of one-point statistics, such as cluster counts (Lacasa & Rosenfeld 2016), and 3-point statistics, such as the weak lensing bispectrum (Rizzato et al. 2019), as well.

Here I build upon Lacasa (2018) and its exhaustive analytical derivation of the covariance of the galaxy angular power spectrum with all non-Gaussian terms. Specifically, I implement the terms derived and argued there as potentially being of importance and I gauge their impact on the information content of galaxy clustering. To this end, I also use the halo model at tree-level for the prediction of the observable $C_{\ell }^{\text{gal}}$ and of the covariance 𝒞_{ℓ, ℓ′}.

In detail, I first carried out an analytical study of non-Gaussian covariance terms of the galaxy angular power spectrum (Sect. 2), recalling the analytical expressions from Lacasa (2018) (Sect. 2.1). Then I presented an approximation for braiding covariance, making it numerically tractable (Sect. 2.2). I then presented numerical results that first demonstrate the importance of braiding covariance (Sect. 3.1) and then show analytically that accounting for braiding covariance is necessary for the inclusion of other in-survey covariance terms, such as 2h1+3, which have important off-diagonal contributions (Sect. 3.2). I present a signal-to-noise analysis that shows that braiding and in-survey covariance have a substantial impact compared to a Gaussian covariance, although the impact is milder once super-sample covariance is also included (Sect. 3.3). Afterwards, I move to a Fisher analysis to show the impact of non-Gaussianity on parameter constraints, both for cosmology (Sect. 4.2) and for halo occupation distribution (HOD, Sect. 4.3). Finally, I discuss the results in Sect. 5 and, in particular, I consider how parameter degeneracies are generally eased up by the inclusion of non-Gaussianity. The data and a Python notebook that allows to reproduce all plots and results of the article are available online, along with a bit more information¹.

2. Analytical covariance

In this section, I first set out the equations for the non-Gaussian covariance terms, then I present a numerical approximation for the specific case of braiding covariance. For this purpose, a few definitions and notations are needed.

First, the (unobservable) angular power spectrum of matter between two redshifts (z_a, z_b) is

$$\begin{array}{c}\hfill C_{\ell }^m(z_a,z_b)=\frac{2}{\pi }{\displaystyle \int k^2\mathrm{d}kP(k|z_a,z_b)j_{\ell }(k{r}_a)j_{\ell }(k{r}_b).}\end{array}$$(3)

Second, halo model equations can be greatly simplified by introducing the integral

$$\begin{array}{cc}\hfill I_{\mu }^{\beta }(k_1,\cdots ,k_{\mu }|z)& ={\displaystyle \int \mathrm{d}M\frac{\mathrm{d}{n}_h}{\mathrm{d}M}\langle N_{\mathrm{gal}}^{(\mu )}\rangle b_{\beta }(M,z)}\hfill \\ \hfill & \times u(k_1|M,z)\cdots u(k_{\mu }|M,z)\hfill \end{array}$$(4)

where $\frac{\text{d}{n}_h}{\text{d}M}$ is the halo mass function, u(k|M, z) is the normalised halo profile, b_β(M, z) is the halo bias or order β², and $\langle N_{\mathrm{gal}}^{(n)}\rangle \equiv \langle N_{\mathrm{gal}}(N_{\mathrm{gal}}-1)\cdots (N_{\mathrm{gal}}-(n-1))\rangle$ is the number of n-tuples of galaxies, implicitly depending on halo mass.

Finally, further simplifications can be achieved by grouping integrals together:

$$\begin{array}{c}\hfill I_{\mu }^{{\mathrm{\Sigma }}_2}=\frac{17}{21}{I}_{\mu }^1+\frac{1}{2!}{I}_{\mu }^2\end{array}$$(5)

is the sum of second order contributions from perturbation theory and local bias, and

$$\begin{array}{c}\hfill I_{\mu }^{{\mathrm{\Sigma }}_3}\equiv \frac{1023}{1701}{I}_{\mu }^1+\frac{1}{3!}{I}_{\mu }^3\end{array}$$(6)

is the sum of third order contributions (Lacasa 2018).

2.1. Non-Gaussian terms

I recapitulate the equations for all the non-Gaussian covariance terms so that this article may be self-contained. The equations all stem from Lacasa (2018), with the slight modification that they are for the power spectrum of the usual galaxy density contrast, that is, $C_{\ell }^{\text{gal}}\equiv C_{\ell }({\delta }_{\text{gal}})$, instead of the absolute power spectrum C_ℓ(n_gal) used in Lacasa (2018). This is done to maintain maximal familiarity for most readers. In practice, this just changes an overall factor for power spectra and covariances and it does not change parameter constraints that are presented later in Sect. 4 nor any of the conclusion on the importance of the various terms.

The first non-Gaussian covariance term is by far the most studied (e.g. Takada & Hu 2013; Li et al. 2014a,b, 2018; Lacasa & Rosenfeld 2016; Lacasa & Kunz 2017; Lacasa et al. 2018; Akitsu & Takada 2018; Barreira et al. 2018a) and the one whose impact is already well recognised even for some current surveys (e.g. Hildebrandt et al. 2017): super-sample covariance (SSC). It takes the form of

$$\begin{array}{c}\hfill {\mathcal{C}}_{\ell ,{\ell }^{\prime }}^{\mathrm{SSC}}={\displaystyle \int \mathrm{d}{V}_a\mathrm{d}{V}_b{\mathrm{\Psi }}_{\ell }^{\mathrm{sqz}}(z_a){\mathrm{\Psi }}_{{\ell }^{\prime }}^{\mathrm{sqz}}(z_b){\sigma }^2(z_a,z_b)/N_{\mathrm{gal}}{(i_z)}^2{N}_{\mathrm{gal}}{(j_z)}^2,}\end{array}$$

where z_a ∈ i_z, z_b ∈ j_z, $\text{d}{V}_a=r^2(z_a)\frac{\text{d}r}{\text{d}z}(z_a)$,

and

$$\begin{array}{c}\hfill {\sigma }^2(z_a,z_b)=\frac{C_0^m(z_a,z_b)}{4\pi }\end{array}$$(7)

is the SSC kernel, and with the angle-independent trispectrum terms from the halo model, Lacasa (2018) find

$$\begin{array}{c}\hfill {\mathrm{\Psi }}_{\ell }^{\mathrm{sqz}}(z)=4I_1^{{\mathrm{\Sigma }}_2}(k_{\ell }|z)I_1^1(k_{\ell }|z)P(k_{\ell }|z)+I_2^1(k_{\ell },k_{\ell }|z),\end{array}$$(8)

which can be related to the more usual power spectrum response $\frac{\partial P_{\mathrm{gal}}}{\partial {\delta }_b}$ via ${\mathrm{\Psi }}_{\ell }^{\mathrm{sqz}}(z)=\frac{\partial P_{\mathrm{gal}}(k_{\ell })}{\partial {\delta }_b}\times I_1^0{(z)}^2$.

Although a fast approximation to SSC was recently presented by Lacasa & Grain (2019), I prefer to maintain an exact computation here. I checked that the quick approximation gives results within 5% to that of the full computation of Eq. (7) for all numerical results presented throughout the article.

Next we have non-Gaussian terms, coming from the diagonal-independent part of the trispectrum. The first and simplest is the 1-halo term where all galaxies of the 4-point function reside in the same halo,

$$\begin{array}{c}\hfill {\mathcal{C}}_{\ell ,{\ell }^{\prime }}^{1\mathrm{h}}=\frac{{\delta }_{i_z,j_z}}{4\pi }{\displaystyle \int \mathrm{d}V{I}_4^0(k_{\ell },k_{\ell },k_{{\ell }^{\prime }},k_{{\ell }^{\prime }}|z)/N_{\mathrm{gal}}{(i_z)}^4.}\end{array}$$(9)

Then come higher halo terms which should not be included independently, as I show in Sect. 3.2. We have the 2-halo 1+3 term, where one galaxy sits in a halo and the three others sit in another halo,

$$\begin{array}{cc}\hfill {\mathcal{C}}_{\ell ,{\ell }^{\prime }}^{2\mathrm{h}1+3}=\frac{2{\delta }_{i_z,j_z}}{4\pi }& {\displaystyle \int \mathrm{d}V{I}_1^1(k_{\ell }|z)I_3^1(k_{\ell },k_{{\ell }^{\prime }},k_{{\ell }^{\prime }}|z)P(k_{\ell }|z)/N_{\mathrm{gal}}{(i_z)}^4}\hfill \\ \hfill & +(\ell \leftrightarrow {\ell }^{\prime }),\hfill \end{array}$$(10)

the 3-halo base term,

$$\begin{array}{cc}\hfill {\mathcal{C}}_{\ell ,{\ell }^{\prime }}^{3\mathrm{h}-\mathrm{base}0}& =\frac{{\delta }_{i_z,j_z}}{4\pi }{\displaystyle \int \mathrm{d}V2{\left(I_1^1(k_{\ell }|z)P(k_{\ell }|z)\right)}^2{I}_2^{{\mathrm{\Sigma }}_2}(k_{{\ell }^{\prime }},k_{{\ell }^{\prime }}|z)/N_{\mathrm{gal}}{(i_z)}^4}\hfill \\ \hfill & +(\ell \leftrightarrow {\ell }^{\prime })\hfill \\ \hfill & +\frac{4{\delta }_{i_z,j_z}}{4\pi }{\displaystyle \int \mathrm{d}V2I_1^1(k_{\ell }|z)}{I}_1^1(k_{{\ell }^{\prime }}|z)I_2^{{\mathrm{\Sigma }}_2}(k_{\ell },k_{{\ell }^{\prime }}|z)\hfill \\ \hfill & \times P(k_{\ell }|z)P(k_{{\ell }^{\prime }}|z)/N_{\mathrm{gal}}{(i_z)}^4,\hfill \end{array}$$(11)

and the 4-halo term from third order contributions,

$$\begin{array}{cc}\hfill {\mathcal{C}}_{\ell ,{\ell }^{\prime }}^{4\mathrm{h}-3}=& \frac{2{\delta }_{i_z,j_z}}{4\pi }{\displaystyle \int \mathrm{d}V3!{\left(I_1^1(k_{\ell },z)\right)}^2}{I}_1^1(k_{{\ell }^{\prime }},z)I_1^{{\mathrm{\Sigma }}_3}(k_{{\ell }^{\prime }},z)\hfill \\ \hfill & \times P(k_{\ell }|z)P(k_{\ell }|z)P(k_{{\ell }^{\prime }}|z)/N_{\mathrm{gal}}{(i_z)}^4+(\ell \leftrightarrow {\ell }^{\prime }).\hfill \end{array}$$(12)

Finally, the most complicated case is braiding covariance, whose projection in spherical harmonics is found in Lacasa (2018). It has some similarities with SSC in that it is also a class of terms grouped together and it also takes the form of a double redshift integral with the non-linear physics encapsulated in separable elements:

$$\begin{array}{cc}\hfill {\mathcal{C}}_{\ell ,{\ell }^{\prime }}^{\mathrm{braid}}=2& {\displaystyle \int \mathrm{d}{V}_{\mathit{ab}}{\mathrm{\Psi }}_{\ell ,{\ell }^{\prime }}^{\mathrm{alt}}(z_a){\mathrm{\Psi }}_{\ell ,{\ell }^{\prime }}^{\mathrm{alt}}(z_b)}\hfill \\ \hfill & \times {\mathcal{B}}_{\ell ,{\ell }^{\prime }}(z_a,z_b)/N_{\mathrm{gal}}{(i_z)}^2{N}_{\mathrm{gal}}{(j_z)}^2,\hfill \end{array}$$(13)

where

$$\begin{array}{c}\hfill {\mathcal{B}}_{\ell ,{\ell }^{\prime }}(z_a,z_b)={\displaystyle \sum _{{\ell }_a}\frac{2{\ell }_a+1}{4\pi }{\left(\begin{array}{ccc}\ell & {\ell }^{\prime }& {\ell }_a\\ 0& 0& 0\end{array}\right)}^2{C}_{{\ell }_a}^m(z_a,z_b)}\end{array}$$(14)

is the braiding kernel and

$$\begin{array}{c}\hfill {\mathrm{\Psi }}_{\ell ,{\ell }^{\prime }}^{\mathrm{alt}}(z)=[2I_1^{{\mathrm{\Sigma }}_2}(k_{{\ell }^{\prime }}|z)I_1^1(k_{\ell }|z)P(k_{\ell }|z)+(\ell \leftrightarrow {\ell }^{\prime })]+I_2^1(k_{\ell },k_{{\ell }^{\prime }}|z)\end{array}$$(15)

encapsulates the non-linear physics.

2.2. An approximation to braiding covariance

Directly implementing Eq. (13) for braiding covariance is numerically challenging. Indeed, it would need the computation of ℬ_{ℓ, ℓ′}(z_a, z_b) for all pairs of multipoles and all pairs of redshifts. ℬ_{ℓ, ℓ′}(z_a, z_b), itself a sum over 𝒪(ℓ_max) multipoles, quickly makes it a burden for next-gen galaxy surveys where we target ℓ_max = 𝒪(10³).

To overcome this, I devised an approximation with an approach similar to that followed by Lacasa & Grain (2019) for super-sample covariance: we can approximate that ${\mathrm{\Psi }}_{\ell ,{\ell }^{\prime }}^{\mathrm{alt},\mathrm{clust}}/$ ${\overline{n}}_{}\text{gal}{(z)}^2$ varies slowly with redshift compared to ℬ_{ℓ, ℓ′}. Then

$$\begin{array}{c}\hfill {\mathcal{C}}_{\ell ,{\ell }^{\prime }}^{\mathrm{braid}}\approx 2{\mathrm{\Psi }}_{\ell ,{\ell }^{\prime }}^{\mathrm{alt},\mathrm{int}}(i_z){\mathrm{\Psi }}_{\ell ,{\ell }^{\prime }}^{\mathrm{alt},\mathrm{int}}(j_z)B_{\ell ,{\ell }^{\prime }}(i_z,j_z),\end{array}$$(16)

where

$$\begin{array}{c}\hfill {\mathrm{\Psi }}_{\ell ,{\ell }^{\prime }}^{\mathrm{alt},\mathrm{int}}(i_z)={\displaystyle \int \mathrm{d}V{\mathrm{\Psi }}_{\ell ,{\ell }^{\prime }}^{\mathrm{alt}}(z)}\end{array}$$(17)

and

$$\begin{array}{cc}\hfill B_{\ell ,{\ell }^{\prime }}(i_z,j_z)=& {\displaystyle \int \mathrm{d}{V}_{\mathit{ab}}{\overline{n}}_{\mathrm{gal}}{(z_a)}^2{\overline{n}}_{\mathrm{gal}}{(z_b)}^2{\mathcal{B}}_{\ell ,{\ell }^{\prime }}(z_a,z_b)/\left(I^{n_g^2}(i_z)I^{n_g^2}(j_z)\right)}\hfill \\ \hfill =& {\displaystyle \sum _{{\ell }_a}\frac{2{\ell }_a+1}{4\pi }{\left(\begin{array}{ccc}\ell & {\ell }^{\prime }& {\ell }_a\\ 0& 0& 0\end{array}\right)}^2{C}_{{\ell }_a}^{n_g^2}(i_z,j_z)}\hfill \end{array}$$(18)

with

$$\begin{array}{c}\hfill I^{n_g^2}(i_z)={\displaystyle {\int }_{z\in i_z}\mathrm{d}V{\overline{n}}_{\mathrm{gal}}{(z)}^2}\end{array}$$(19)

and

$$\begin{array}{c}\hfill C_{\ell }^{n_g^2}(i_z,j_z)={\displaystyle \int \mathrm{d}{V}_{\mathit{ab}}{\overline{n}}_{\mathrm{gal}}{(z_a)}^2{\overline{n}}_{\mathrm{gal}}{(z_b)}^2{C}_{\ell }^m(i_z,j_z)/\left(I^{n_g^2}(i_z)I^{n_g^2}(j_z)\right).}\end{array}$$(20)

I call this the “Bij approximation” for Braiding covariance, similarly to the name “Sij approximation” for super-sample covariance. The fact that the Sij approximation works very well (see Lacasa & Grain 2019) proves that the Bij should work equally well, if not better. Indeed, the similarity between the separable elements Ψ^sqz and Ψ^alt3 and the fact that ℬ_0, 0(z, z′) = σ²(z, z′) shows that ℬ_{ℓ, ℓ′} varies quickly enough with redshift for the Bij approximation to work at ℓ = ℓ′ = 0. And at higher multipoles, ℬ_{ℓ, ℓ′} only varies more quickly, making the approximation increasingly more precise. Indeed, from Eq. (14), at high (ℓ,ℓ′) ℬ_{ℓ, ℓ′} gets contributions from $C_{{\ell }_a}^m(z,z^{\prime })$ at high ℓ_a, which gets increasingly close to a Dirac δ(z, z′) due to Limber approximation. These analytical arguments ensure that the Bij approximation for Braiding covariance works at least as well as the Sij approximation for SSC.

3. Covariance results and the importance of braiding for positive definiteness

In this section, I first present the physical and technical assumptions I used for the computation of the galaxy angular power spectrum and its covariance terms, along with the numerical results for the covariances. Then I show why these results prove the importance of including some of the non-Gaussian terms presented in Sect. 2: braiding and 2h1+3. Finally, I present the impact of NG terms on the measurement signal to noise ratio of the galaxy angular power spectrum.

3.1. Setup and covariances

For the numerical results presented in this and later sections, I used a flat ΛCDM cosmological model with Planck 2018 (Planck Collaboration VI 2018) cosmological parameters (Ω_bh², Ω_ch², H₀, n_S, σ₈) = (0.022, 0.12, 67, 0.96, 0.81). For the halo model, I adopted the Tinker et al. (2008) halo mass function with the corresponding halo bias from Tinker et al. (2010). For the HOD, I used one similar to Zehavi et al. (2011): N_gal = N_cen + N_sat, with a Bernoulli distribution for the central galaxy with probability,

$$\begin{array}{c}\hfill P(N_{\mathrm{cen}}=1)=\frac{1}{2}(1+\mathrm{Erf}\left(\frac{{log}_{10}M-{log}_{10}{M}_{\min }}{{\sigma }_{\log M}}\right))\end{array}$$(21)

and a Poisson distribution for the satellite galaxies, conditioned to the presence of the central, with mean,

$$\begin{array}{c}\hfill \mathbb{E}[N_{\mathrm{sat}}|N_{\mathrm{cen}}=1]={\left(\frac{M}{M_{\mathrm{sat}}}\right)}^{{\alpha }_{\mathrm{sat}}}·\end{array}$$(22)

In this section, I consider a single redshift bin for the galaxies: 0.9 <  z <  1.019. For the HOD parameters, I used log₁₀M_min = 11.3, σ _{log M} = 0.5, M_sat = 10 × M_sat and α_sat = 1. These parameters predict a galaxy density at these redshifts equal to the predicted one for the Euclid photometric sample, that is, 3 galaxies/arcmin² (see Appendix A) which corresponds to a total of ∼450 M galaxies as I assume a full sky setup.

With these parameters, I computed the galaxy angular power spectrum $C_{\ell }^{\text{gal}}$ and the different non-Gaussian covariance contributions listed in Sect. 2 for nine individual multipoles distributed logarithmically in [30,3000]. The variance per multipole created by each term is shown in Fig. 1 plotted as a function of multipole ℓ.

Fig. 1. Different non-Gaussian contributions to the variance of the angular power spectrum per multipole in the redshift bin 0.9 <  z <  1.019.

We first see that the 3h-base0 and 4h-3 terms are negligible compared to all other terms. This means that the perturbative contributions to variances are excellently encapsulated inside super-sample covariance and braiding covariance. We can then focus on the other covariance terms considered in this article: braiding and 2h1+3. We see that braiding is actually the dominant NG contribution to the variance on large scales and remains non-negligible on most of the multipole range. The 2h1+3 term is subdominant everywhere, but it still is not negligible. I emphasise that these results are not enough to draw conclusions on the importance of the terms as they only show the diagonal, rather than the whole structure of the covariance matrices.

To examine the covariance matrices and be more representative of a survey analysis, I needed to consider not only a few multipoles but the full multipole range. Computing the covariance matrices for all single multipoles in this range is not desirable, however, because (i) it is very intensive numerically and (ii) it would not be representative of actual data analysis that bins multipoles together. Hence, I performed a binning of multipoles, which consisted of interpolating and binning from the nine original multipoles to 29 bins distributed logarithmically Δℓ/ℓ=cst in the range ℓ ∈ [32, 2290]. Hereafter, binned quantities are plotted with the indication of the central multipole of the bin, defined as the geometrical average of the bin stakes.

With these specifications, I show in Fig. 2 the correlation matrices: $C_{i,j}/\sqrt{C_{i,i}{C}_{j,j}}$ for each of the non-Gaussian covariance terms. Each term is normalised by its own diagonal to reveal its specific structure. I note that this is different from the more customary normalisation by the total diagonal, which lets us appreciate the relevance of the terms; however, this relevance will be addressed later, in Sects. 3.3 and 4.

Fig. 2. Correlation matrices for the different non-Gaussian covariance terms, normalised by its own diagonal. Top: SSC, 1h, Braiding. Bottom: 2h1+3, 3h-base0, 4h-3; the color bar is clipped at 7.

In the top row we see well-behaved terms which yield matrices with all eigenvalues ≥0 : SSC, 1-halo and braiding. The correlation coefficients are all in [−1,1]. In the bottom row we see the 2h1+3, 3h-base0 and 4h-3 terms for which the correlation coefficients can be > 1 (up to 7.4 for 2h1+3, 39 for 3h-base0 and 7.5 for 4h-3; the color bar is clipped to 7 in the plots for readability), indicating that these matrices have negative eigenvalues.

3.2. Importance of braiding for positive definiteness

In this section I examine the problem of the NG terms with negative eigenvalues: 2h1+3, 3h-base0, and 4h-3. I first give an analytical explanation why they yield, alone, correlation coefficients > 1, then I give a physical explanation why they cannot be included alone and argue why Braiding covariance is necessary to regulate them to obtain a well-behaved total covariance matrix, that is, positive definite.

First, let us become convinced in an analytical sense that the correlation coefficients > 1 seen in the bottom row of Fig. 2 are physical and not a bug in my computation. For this, I focus on the case of the 2h1+3 term. Both for simplicity, so as not to repeat similar computations thrice, and because it dominates the 3h-base0 and 4h-3 terms as seen in Fig. 1.

Let us evaluate the correlation coefficient

$$\begin{array}{c}\hfill r_{\ell ,{\ell }^{\prime }}=\frac{{\mathcal{C}}_{\ell ,{\ell }^{\prime }}}{\sqrt{{\mathcal{C}}_{\ell ,\ell }\times {\mathcal{C}}_{{\ell }^{\prime },{\ell }^{\prime }}}}\end{array}$$(23)

for the 2h1+3 term Eq. (10) in the following case: infinitesimally small redshift bins and k_ℓ, k_ℓ′ ≪ 1/R, where R is the typical radius of a halo, so that u(k)→1. These conditions mean that the redshift integrals can be replaced by a multiplication with Δz (which vanishes in the ratio) and that all halo model integrals $I_{\mu }^{\beta }$ are independent of ℓ, ℓ′. Then we get

$$\begin{array}{cc}\hfill {\left(r_{\ell ,{\ell }^{\prime }}^{2\mathrm{h}1+3}\right)}^2\approx & \frac{{(I_1^1{I}_3^{1}P(k_{\ell })+(\ell \leftrightarrow {\ell }^{\prime }))}^2}{2I_1^1{I}_3^{1}P(k_{\ell })\times 2I_1^1{I}_3^{1}P(k_{{\ell }^{\prime }})}=\frac{1}{4}\frac{{(P(k_{\ell })+P(k_{{\ell }^{\prime }}))}^2}{P(k_{\ell })P(k_{{\ell }^{\prime }})}·\hfill \end{array}$$(24)

Now I further take the condition k_eq <  k_ℓ ≪ k_ℓ′, where k_eq is the position of the maximum of the matter power spectrum P(k) (corresponding to matter-radiation equality) so that both wave vectors are in the decreasing part of P(k). In that case P(k_ℓ)≫P(k_ℓ′) and we get

$$\begin{array}{cc}\hfill r_{\ell ,{\ell }^{\prime }}^{2\mathrm{h}1+3}\approx & \frac{1}{2}\sqrt{\frac{P(k_{\ell })}{P(k_{{\ell }^{\prime }})}}>1.\hfill \end{array}$$(25)

So the result is physical: alone these covariance terms give correlation coefficients which can be > 1. This means that these terms yield incorrect covariance matrices if left alone: two measurements can be more than 100% correlated, or in other term the matrix restricted to these two points has a negative eigenvalue.

This result can also be understood more visually by using the diagrammatic formalism built by Lacasa et al. (2014). As shown by Lacasa (2018), the 4h-3 is part of the terms of the left diagram of Fig. 3, which quantifies how the 2-halo part of the spectrum is correlated with itself due to halos being clustered in a (non-Gaussian) matter field. The 3h-base0 is part of the terms of the central diagram, which quantifies how the 2-halo part of the spectrum is correlated with the 1-halo part due to halos being clustered in a (non-Gaussian) matter field. And the 2h1+3 is the entirety of the terms of the right diagram, which quantifies how the 2-halo part of the spectrum is correlated with the 1-halo part due to halo coincidence. From these diagrams it becomes clear that the 2h1+3 term is going to be maximal when ℓ is in the large-scale 2-halo dominated regime while ℓ′ is in the small-scale 1-halo dominated regime. So this term is going to yield high covariance when ℓ′≫ℓ and minimal covariance when ℓ = ℓ′, i.e. exactly the off-diagonal behaviour we see in Fig. 2.

Fig. 3. Diagrams for some of the trispectrum terms involved in the covariance of the galaxy angular power spectrum $\text{Cov}(C_{\ell }^{\mathrm{gal}},C_{{\ell }^{\prime }}^{\mathrm{gal}})$. From left to right: 4-halo, 3-halo, and 2-halo 1+3 term. Galaxies 1 and 2 are the source of the first power spectrum $C_{\ell }^{\text{gal}}$, while galaxies 3 and 4 are the source of the second power spectrum $C_{{\ell }^{\prime }}^{\mathrm{gal}}$.

Now this behaviour has to be regulated by another covariance term which makes the total covariance matrix well-behaved. Mathematically, the regulator cannot be the Gaussian part of the covariance, nor SSC, nor the 1h trispectrum term alone. First, it cannot be the Gaussian part of the covariance. Indeed, going to arbitrarily high redshifts, we can have arbitrarily high multipoles ℓ′ that fulfill the conditions k_ℓ′ ∼ ℓ′/r(z)≪1/R. At these multipoles, the Gaussian variance becomes negligible since it decreases as 1/(2ℓ′ + 1). Second, this cannot either be the super-sample covariance. Indeed, SSC gives a near degenerate covariance matrix with a single positive eigenvalue, the other being zero, as seen from Fig. 2 where the correlation matrix is 100% everywhere. So SSC cannot regulate a multitude of negative eigenvalues. Finally, for the same reason, the regulator cannot either be the 1h trispectrum term, which is constant on large scales.

We can find the regulator via the diagram discussion. Since the 2h1+3 term quantifies how the 2-halo part of the spectrum is correlated with the 1-halo part due to halo coincidence, it has to be regulated by a first term which quantifies how the 2-halo part of the spectrum is correlated with itself due to halo coincidence, and a second term which quantifies how the 1-halo part of the spectrum is correlated with itself due to halo coincidence. The first wanted term is part of braiding covariance: it is the 2-halo part of Braiding, which corresponds to the left diagram of Fig. 4. The second wanted term is the 1-halo trispectrum term, which corresponds to the right diagram of Fig. 4.

Fig. 4. Diagrams for some of the trispectrum terms involved in the covariance of the galaxy angular power spectrum. Left: 2-halo part of Braiding, right: 1-halo term.

With similar considerations, we can see that the regulator of the 3h-base0 and 4h-3 terms is Braiding covariance. So its is the sum of the 1h, Braiding, 2h1+3, 3h-base0, and 4h-3 terms that yield a well-behaved covariance. In the following I call this sum “other non-Gaussianity” (ONG) by contrast with the non-Gaussian covariance that has been the most studied to date: super-sample covariance. For comparison, in the previous literature non-SSC NG terms have also been called “in-survey” (e.g. Rizzato et al. 2019), “connected non-Gaussian (cNG)” (e.g. Barreira et al. 2018a), “trispectrum” (e.g. Li et al. 2014a) or “T0” (e.g. Wadekar & Scoccimarro 2019). Figure 5 shows the correlation matrix for the ONG group.

Fig. 5. Correlation matrix for ONG, the sum of non-Gaussian covariance terms other than SSC: 1h, 2h1+3, 3h-base0, 4h-3, and Braiding.

We see that the ONG indeed has all correlation coefficient ≤100%. Furthermore, numerical investigation shows that all eigenvalues are > 0. Thus, the addition of braiding covariance has correctly regulated the off-diagonal components of the 2h1+3, 3h-base0, and 4h-3 terms. I conclude that the inclusion of Braiding is necessary to go beyond the current state of the art for non-Gaussian covariances.

3.3. Impact on the signal-to-noise ratio

Alhough Braiding is necessary to include ONG covariance, the question remains of whether ONG has a significant impact on the information content of the galaxy angular power spectrum. In this section, I use the signal-to-noise ratio (S/N)

$$\begin{array}{c}\hfill {\left(\frac{S}{N}\right)}^2={\displaystyle \sum _{\ell ,{\ell }^{\prime }={\ell }_{\min }}^{{\ell }_{\max }}{C}_{\ell }^{\mathrm{gal}}{\mathcal{C}}_{\ell ,{\ell }^{\prime }}^{-1}{C}_{{\ell }^{\prime }}^{\mathrm{gal}}}\end{array}$$(26)

as a first metric to quantify this information content, as already used in the literature (e.g. Rizzato et al. 2019).

To this end, I also use the halo model at tree-level to predict the power spectrum $C_{{\ell }^{\prime }}^{\mathrm{gal}}$. This modelling allows for ∼10% precision; for future surveys, this is sufficient for the prediction of the covariance, but not for the prediction of the power spectrum. This is, however, not an issue for this analysis as my goal here is to gauge the relative impact of covariance terms. Figure 6 shows S/N plotted as a function of ℓ_max for different degree of sophistication in the computation of the covariance.

Fig. 6. Cumulative signal to noise ratio for the measurement of $C_{\ell }^{\text{gal}}$ in the bin 0.9 <  z <  1.019 as a function of maximum multipole of analysis. Left, from top to bottom: Gaussian covariance only, Gaussian + “other NG”, Gaussian + SSC, Gaussian + SSC + 1h, total covariance. Right: zoom on the three lowest curves: Gaussian + SSC, Gaussian + SSC + 1h, total covariance, all normalised by the value of the signal to noise using the full multipole range and the total covariance.

If the analysis is carried out on the full range of multipoles, as is scheduled, for instance, for Euclid, then non-Gaussian covariance terms have a large impact on the information content. Compared to the Gaussian case, ONG alone decreases S/N by a factor 1.7. This is clearly a large impact, and one must go beyond Gaussian covariances. Now the current state of the art includes super-sample covariance, and that term has a larger impact: SSC alone decreases S/N by a factor 3.1. Finally, when accounting for the total covariance: Gaussian+SSC+ONG, S/N decreases by a factor 3.4 compared to the Gaussian case. So ONG has a 9.4% impact on top of SSC. The 1h covariance has a negligible impact on top of SSC, so the bulk of the 9.4% impact comes from the Braiding and 2h1+3 terms.

Thus, including ONG seems fairly important (if SSC is already accounted for) given, for example, that Euclid has a requirement of 10% precision on error bars. First, I argue that ONG should still be accounted for because it makes the information systematically lower and, thus, error bars become systematically larger. Second, this section used the S/N in a single redshift bin as a metric and the question remains open of the impact on parameter constraints when summed over the entire redshift range. This is the subject of the next section.

4. Impact on parameter constraints

4.1. Setup

I use survey specifications representative of the Euclid photometric galaxy sample (Euclid Collaboration 2019): sky coverage f_SKY = 0.36, a galaxy redshift distribution of

$$\begin{array}{c}\hfill n(z)\propto {\left(\frac{z}{z_0}\right)}^{2}exp[-{\left(\frac{z}{z_0}\right)}^{3/2}],\end{array}$$(27)

where $z_0=z_{\mathrm{m}}/\sqrt{2,}$ with z_m = 0.9 the median redshift (Laureijs et al. 2011). The total density is 30 gals⋅arcmin⁻² in the redshift range [0,2.5]. Following Euclid Collaboration (2019), the sample is divided into 10 equi-populated redshift bins, whose bin stakes are z = 0.001,0.418,0.56,0.678,0.789,0.9,1.019,1.155,1.324,1.576,2.5⁴. To reproduce this redshift distribution with the halo model, I use the Halo Occupation Distribution described in Sect. 3.1, further including a redshift dependence of M_min in the form:

$$\begin{array}{c}\hfill M_{\min }(z)=M_{\min }^a+M_{\min }^{b}z+M_{\min }^c{z}^2+M_{\min }^d{z}^3.\end{array}$$(28)

As shown in Appendix A, this parametrisation allows to reproduce the Euclid-expected galaxy counts to 2.5% precision, and predicts a galaxy bias consistent with simulations.

In this section, I quantify the impact of covariances on parameter constraints using the methodology of Fisher forecasts. To this end, I use both Fisher matrices in a given redshift bin:

$$\begin{array}{c}\hfill F_{\alpha ,\beta }(i_z)={\displaystyle \sum _{\ell ,{\ell }^{\prime }={\ell }_{\min }}^{{\ell }_{\max }}{\partial }_{\alpha }{C}_{\ell }^{\mathrm{gal}}(i_z){\mathcal{C}}_{\ell ,{\ell }^{\prime }}^{-1}(i_z,i_z){\partial }_{\beta }{C}_{{\ell }^{\prime }}^{\mathrm{gal}}(i_z)}\end{array}$$(29)

and summed over all bins:

$$\begin{array}{c}\hfill F_{\alpha ,\beta }={\displaystyle \sum _{i_z,j_z}\sum _{\ell ,{\ell }^{\prime }={\ell }_{\min }}^{{\ell }_{\max }}{\partial }_{\alpha }{C}_{\ell }^{\mathrm{gal}}(i_z){\mathcal{C}}_{\ell ,{\ell }^{\prime }}^{-1}(i_z,j_z){\partial }_{\beta }{C}_{{\ell }^{\prime }}^{\mathrm{gal}}(j_z),}\end{array}$$(30)

where α, β are model parameters, that is, cosmological and HOD parameters in the following : (Ω_bh², Ω_ch², H₀, n_S, σ₈, w₀) and $({\alpha }_{\text{sat}},{\sigma }_{\log M},M_{\text{ratio}},M_{\text{min}}^a,M_{\text{min}}^b,M_{\text{min}}^c,M_{\text{min}}^d)$; ∂_α is the derivative of the observable w.r.t. parameter α.

4.2. Impact on cosmological parameters

We can first look at the Fisher matrix elements in a given redshift bin. For the purposes of illustration, I chose the bin 0.9 <  z <  1.019, which is the same bin as in Sect. 3, containing the median redshift of the galaxy sample and whose results I found representative of the whole sample. Figure 7 shows, as a function of the maximum multipole of analysis ℓ_max, the square root of the Fisher elements for each cosmological parameter of the wCDM model. This quantity is the inverse of the error bar on the considered parameter if all other (cosmological and HOD) parameters were perfectly known.

Fig. 7. (Square root of the) cumulative Fisher elements for the cosmological parameters in the considered redshift bin, as a function of the maximum multipole of the analysis.

If the analysis is carried out on the full range of multipoles then non-Gaussian covariance terms would have a mild impact on the information content for the three first parameters: Ω_bh², Ω_ch², and h, with ONG being more significant than SSC. By contrast, non-Gaussian terms have a large impact on the three last parameters: σ₈, n_S, and w₀. These three latter parameters are arguably the most interesting to constrain with surveys of the large scale structure. The measurement of σ₈ is interesting in the context of the current tension between local measurements and the CMB. The parameter n_S helps to constrain inflation and can be seen as representative of parameters in a more extended model that would change the shape of the power spectrum, for example, a running of the spectral index or massive neutrinos. Finally, the equation of state of dark energy is one of the main science drivers of current and future galaxy surveys.

Compared to the Gaussian Fisher matrix, ONG alone decreases the Fisher content on dark energy $F_{\mathit{w},\mathit{w}}^{1/2}$ by a factor 1.8; for other parameters, the factor ranges between 1.08 (for h) and 1.8 (for σ₈). Super-sample covariance decreases the information on dark energy by a factor 2.9; for other parameters the factor ranges between 1.01 (for Ω_bh²) and 2.6 (for σ₈). The total NG decreases $F_{\mathit{w},\mathit{w}}^{1/2}$ by a factor 3.3; for other parameters the factor ranges between 1.08 (for h) and 3.1 (for σ₈). When compared to Gaussian+SSC, ONG has a 14% impact on $F_{\mathit{w},\mathit{w}}^{1/2}$ ; for other parameters, the impact ranges between 5.6% (for h) and 16% (for σ₈). As in the case of Sect. 3.3, the 1h covariance has a negligible impact on top of SSC so the bulk of the ONG impact comes from the braiding and 2h1+3 terms.

In a second step, I compute the Fisher matrix summed over all redshift bins. This represents the full constraining power of the mock survey; it allows for the breaking of parameter degeneracies, in particular, between parameters for the redshift dependence of the HOD which are nearly completely degenerate in a single bin. In Fig. 8 I plot the marginalised error bars ${\sigma }_i=\sqrt{{(F^{-1})}_{\mathit{ii}}}$ for each cosmological parameter as a function of the maximum multipole of analysis ℓ_max.

Fig. 8. Marginalised error bars on each cosmological parameters using all redshift bins as a function of the maximum multipole of analysis.

When using the full multipole range, non-Gaussian covariance terms have a large impact on the information content for all cosmological parameters. Compared to the Gaussian case, ONG alone increases the error bar on w by 50%; for other parameters, the impact ranges between 14% (for n_S) and 41% (for h). SSC increases σ_w by 88%; for other parameters, the impact ranges between 1.6% (for n_S) and 65% (for h). The total NG increases σ_w by 117%; for other parameters, the impact ranges between 15% (for n_S) and 79% (for h). When compared to Gaussian+SSC, ONG has a 15% impact on σ_w ; for other parameters the impact ranges between 5.7% (for Ω_bh²) and 13% (for n_S). The ONG impact exceeds the threshold of 10% (Euclid precision requirement) for two parameters: n_S and w (σ₈ being affected at 9.6%).

It is interesting to note that ONG has a larger impact than SSC on n_S. This happens because at first order, SSC erases information on the amplitude of the power spectrum (and the redshift dependence of this amplitude) as SSC is 100% correlated. Once we have marginalised over σ₈, this amplitude erasing does not affect n_S, hence, the small (1.6%) impact of SSC on n_S. By contrast, the ONG correlation matrix has a more complex structure and contains terms that couple large and small scale measurements. This affects the lever arm necessary to constrain n_S more heavily. Thus, we can anticipate that other parameters which affect the shape of the matter power spectrum, such as a running of the spectral index or massive neutrinos, would also be more affected by ONG than by SSC.

Finally, Fig. 9 shows the Fisher plot with parameter probability distribution functions (PDFs) and 2σ ellipses that allow for parameter degeneracies to be seen for cosmological constraints using the full multipole range and marginalised over HOD parameters with flat priors. For readability, I did not include the case of Gaussian+SSC+1h, which gives curves nearly identical to the Gaussian+SSC case.

Fig. 9. Fisher ellipses on cosmological parameters, using all redshift bins and the full multipole range. The color coding is identical to the other figures: blue=Gaussian, orange=Gaussian+ONG, green=Gaussian+SSC, violet=total covariance.

We see that PDFs are progressively widened by non-Gaussianities. Furthermore, parameter degeneracies can be affected, sometimes in non-trivial way. For instance the direction of the degeneracy between w and Ω_ch² reverses, though the degeneracy is weak. Additionally, for the strength of the degeneracy, as evidenced by the ellipticity of the Fisher ellipses, it decreases slightly when including NG between n_S and w, but it increases significantly between Ω_bh², Ω_ch² and h. This latter effect dominates the total amount of degeneracy as measured by the condition number of the Fisher matrix, which increases from 4.8 × 10⁷ in the Gaussian case to 1.4 × 10⁸ in the full non-Gaussian case. This is discussed in more details in Sect. 5.

4.3. Impact on halo occupation distribution parameters

We first look at the Fisher matrix elements in the redshift bin 0.9 <  z <  1.019. Figure 10 shows, as a function of the maximum multipole of analysis ℓ_max, the square root of the Fisher elements for each HOD parameter. This quantity is the inverse of the error bar on the considered parameter if all other (cosmological and HOD) parameters were perfectly known.

Fig. 10. (Square root of the) cumulative Fisher elements for the HOD parameters in the considered redshift bin as a function of the maximum multipole of analysis.

If the analysis is carried out on the full range of multipoles, then non-Gaussian covariance terms have a large impact on the information content for all parameters. Compared to the Gaussian case, ONG alone decreases the Fisher content on α_sat, $F_{{\alpha }_{\text{sat}},{\alpha }_{\text{sat}}}^{1/2}$, by a factor 1.8; for other parameters, this factor is the same to the first decimal, ranging between 1.76 and 1.79. SSC decreases the information on α_sat by a factor 2.6; for other parameters, the factor ranges between 2.7 (for M_ratio) and 3.1 (all parameters for the redshift dependence of M_min). The total NG decreases $F_{{\alpha }_{\text{sat}},{\alpha }_{\text{sat}}}^{1/2}$ by a factor 3; for other parameters, the factor ranges between 3.1 (for M_ratio) and 3.4 (all parameters for the redshift dependence of M_min). When compared to Gaussian+SSC, ONG has a 17% impact on $F_{{\alpha }_{\text{sat}},{\alpha }_{\text{sat}}}^{1/2}$; for other parameters the impact ranges between 10% (all parameters for the redshift dependence of M_min) and 16% (for M_ratio). As in the case of cosmological parameters (Sect. 4.2) and the S/N (Sect. 3.3), the 1h covariance has a negligible impact on top of SSC, so the bulk of the ONG impact comes from the braiding and 2h1+3 terms.

I now move to the Fisher matrix summed over all redshift bins. In Fig. 11, I plot the marginalised error bars ${\sigma }_i=\sqrt{{(F^{-1})}_{\mathit{ii}}}$ for each HOD parameter, as a function of the maximum multipole of analysis ℓ_max.

Fig. 11. Marginalised error bars on each HOD parameters using all redshift bins as a function of the maximum multipole of analysis.

When using the full multipole range, non-Gaussian covariance terms have a large impact on the information content for all HOD parameters. Compared to the Gaussian case, ONG alone increases the error bar on α_sat by 19% ; for other parameters, the impact ranges between 12% (for M_ratio) and 39% (for $M_{\text{min}}^d$). SSC increases σ_αsat by 9%; for other parameters, the impact ranges between 7% (for σ _{log M}) and 64% (for $M_{\text{min}}^d$). The total NG increases σ_αsat by 24%; for other parameters, the impact ranges between 17% (for M_ratio) and 85% (for $M_{\text{min}}^d$). When compared to Gaussian+SSC, ONG has a 13% impact on σ_αsat ; for other parameters. the impact ranges between 7.4% (for M_ratio) and 13% (for $M_{\text{min}}^d$). The ONG impact is generally stronger than for cosmological parameters, exceeding the threshold of 10% (Euclid precision requirement) for four parameters: α_sat, $M_{\text{min}}^b$, $M_{\text{min}}^c$ and $M_{\text{min}}^d$.

Interestingly, the impact of ONG is greater than that of SSC for four parameters: α_sat, σ _{log M}, M_ratio, and $M_{\text{min}}^a$. Furthermore. we can note that, for $M_{\text{min}}^b$, $M_{\text{min}}^c$, and $M_{\text{min}}^d$, the inclusion of the 1h covariance makes a visible difference on top of SSC for once, although the rest of in-survey covariance and braiding are necessary to reproduce the full error bar.

Finally, Fig. 12 shows the Fisher plot with parameter PDFs and 2σ ellipses for HOD constraints using the full multipole range and marginalised over cosmological parameters with flat priors.

Fig. 12. Fisher ellipses on HOD parameters, using all redshift bins and the full multipole range. The color coding is identical to the other figures: blue=Gaussian, orange=Gaussian+ONG, green=Gaussian+SSC, violet=total covariance.

Again, PDFs are progressively widened by non-Gaussianities. Furthermore, the strength of parameter degeneracies is generally eased by NG. This is evidenced by the condition number of the Fisher matrix, which decreases from 5.6 × 10⁷ in the Gaussian case to 3.9 × 10⁷ in the full non-Gaussian case.

5. Discussion

By way of a summary of previous results, I have developed an implementation of non-Gaussian covariance terms for galaxy clustering that were initially derived in Lacasa (2018). I developed a numerically tractable approximation for braiding covariance and shown that this class of terms is necessary to include other in-survey covariance terms. Grouping braiding and in-survey under the term ONG covariance, I then studied its impact on S/N analysis and Fisher forecast on the wCDM model with the angular power spectrum with Euclid-like galaxy specifications.

ONG by itself has a large impact on all astrophysical and cosmological parameters, ranging between 12% and 50%. This impact is lowered to some extent by the other NG contender: SSC, which is already included in some current analyses. Compared to this Gaussian+SSC state of the art, ONG still has a significant impact on the covariance, a result in agreement with Barreira et al. (2018b) for weak lensing ; it can even dominate SSC in some configurations, a result that is in agreement with Wadekar & Scoccimarro (2019), which appeared after the first version of this article came out. For parameter constraints, the impact on marginalised error bars ranges between 6% and 15% ; it exceeds 10% –Euclid precision requirement– for the majority of HOD parameters and a couple of cosmological parameters of the wCDM model.

A parameter of particular interest is n_S, whose constraints are significantly affected by ONG. As SSC mostly impacts information on the power spectrum amplitude in opposition to its shape, I expect that ONG should also affect other extensions of the standard cosmological model that change the shape of the matter power spectrum, such as massive neutrinos, warm dark matter and a running of the spectral index.

Interestingly, the increase of error bars due to NG is stronger when the other parameters are fixed, and less strong after marginalisation⁵, as evidenced by Table 1. This happens because the Gaussian Fisher matrix generally has stronger parameter degeneracies compared to the non-Gaussian covariance. The inclusion of NG often increases the minor axis of the Fisher ellipses more than the major axis, leading to a decrease of ellipticity. This is evidenced by the condition number of the whole Fisher matrix (HOD+cosmological parameters) which decreases from 1.0 × 10⁹ in the Gaussian case to 6.5 × 10⁸ in the full non-Gaussian case⁶. Physically, what happens is that with a Gaussian covariance, we erroneously attribute very small error bars to the small scales; so the constraining power is located in a small number of small-scale measurements, leading to parameter degeneracies. By contrast, when NG is accounted for, error bars are increased on small scales so the constraining power is distributed more evenly among scales.

The only exceptions to this argument are Ω_bh², Ω_ch², and h, where the strength of degeneracies is increased by NG. First, I checked that this degeneracy is not an effect of the choice of parameters; it is still present if I use (Ω_b, Ω_c, h) instead of (Ω_bh², Ω_ch², h). Second, this increase of degeneracy happens because these parameters are mostly constrained by the redshift dependence of the clustering amplitude. This information is heavily affected by SSC. In terms of the likelihood approached to SSC developed in Lacasa & Grain (2019), these parameters become degenerate with the redshift evolution of the background change δ_b(z). Indeed, we see from Fig. 9 that the largest increase of the degeneracy comes from SSC.

In looking at the condition numbers, I find that it is worsened by NG for cosmological parameters (4.8 × 10⁷ → 1.4 × 10⁸) and slightly improved by NG for HOD parameters (5.6 × 10⁷ → 3.9 × 10⁷). This means the bulk of the improvement for the whole cosmo+HOD matrix comes from the change in the off-diagonal block, meaning the degeneracies between cosmological and HOD parameters. Visually inspecting the full Fisher matrices, I indeed found that several degeneracies are improved by NG, in particular, those between w and HOD. This means that NG eases up the sensitivity of Dark Energy constraints on HOD parameters and possible modelling uncertainties. This comes from the structure of the covariance and cannot be mimicked, for example, by rescaling the Gaussian covariance by an arbitrary factor which would leave degeneracies untouched.

In conclusion, including braiding and in-survey covariances is a necessity for future high-density galaxy clustering analyses. This is both because it impacts error bars at a level above the precision requirements and also because it renders cosmological constraints more robust for astrophysical uncertainties.

Acknowledgments

I thank Isaac Tutusaus for private communication on the galaxy bias of the Euclid expected photometric sample. Part of this work was supported by funds of the Département de Physique Théorique, Université de Genève. Part of this work was supported by a postdoctoral grant from Centre National d’Études Spatiales (CNES).

References

Appendix A: Redshift dependent halo occupation distribution

The specification Eq. (27) for the galaxy redshift distribution n(z) corresponds to a magnitude-limited sample, and not a volume-limited sample as is normally required for a HOD analysis. To overcome this, I fit the HOD parameters at each redshift. As n(z) is the only “data” available, I can fit a single HOD parameter and I have to assume relations for the others. Specifically, I fit the M_min parameter, assume that the ratio M_ratio = M_sat/M_min = 10 is constant and that σ _{log M} = 0.5 and α_sat = 1 are constant. I then find that I can reproduce the fitted M_min(z) with a fourth order polynomial:

$$\begin{array}{c}\hfill M_{\min }(z)=M_{\min }^a+M_{\min }^{b}z+M_{\min }^c{z}^2+M_{\min }^d{z}^3,\end{array}$$(A.1)

with parameter values (rounded to the third decimal) $M_{\text{min}}^a=11.020$, $M_{\text{min}}^b$ = −0.143, $M_{\text{min}}^c$ = 0.549 and $M_{\text{min}}^d$ = −0.105.

Using this polynomial redshift-dependent HOD (hereafter, polynomial HOD), Fig. A.1 shows the predicted galaxy redshift distribution n(z) compared to the original specification.

Fig. A.1. Left: specification for the galaxy redshift distribution from Laureijs et al. (2011) (blue solid) and prediction from the polynomial HOD (red dashed). Right: ratio of the two distributions.

We see that the redshift distribution is reproduced to better than 2.5% accuracy over the whole redshift range. Now this means that given the galaxy bias at any order, the non-linear power spectrum (etc.) can all be predicted without any additional free parameter. For instance, Fig. A.2 shows the predicted first order galaxy bias.

Fig. A.2. Galaxy bias predicted from the polynomial HOD.

I find agreement between this galaxy bias and preliminary results from Euclid internal simulations (Isaac Tutusaus, priv. comm.), which shows that this simple parametrisation is indeed capable of capturing the redshift evolution of the expected Euclid photometric sample to a satisfactory extent.

All Tables

All Figures

	Fig. 1. Different non-Gaussian contributions to the variance of the angular power spectrum per multipole in the redshift bin 0.9 <  z <  1.019.
In the text

	Fig. 2. Correlation matrices for the different non-Gaussian covariance terms, normalised by its own diagonal. Top: SSC, 1h, Braiding. Bottom: 2h1+3, 3h-base0, 4h-3; the color bar is clipped at 7.
In the text

Fig. 3. Diagrams for some of the trispectrum terms involved in the covariance of the galaxy angular power spectrum $\text{Cov}(C_{\ell }^{\mathrm{gal}},C_{{\ell }^{\prime }}^{\mathrm{gal}})$. From left to right: 4-halo, 3-halo, and 2-halo 1+3 term. Galaxies 1 and 2 are the source of the first power spectrum $C_{\ell }^{\text{gal}}$, while galaxies 3 and 4 are the source of the second power spectrum $C_{{\ell }^{\prime }}^{\mathrm{gal}}$.

In the text

	Fig. 4. Diagrams for some of the trispectrum terms involved in the covariance of the galaxy angular power spectrum. Left: 2-halo part of Braiding, right: 1-halo term.
In the text

	Fig. 5. Correlation matrix for ONG, the sum of non-Gaussian covariance terms other than SSC: 1h, 2h1+3, 3h-base0, 4h-3, and Braiding.
In the text

Fig. 6. Cumulative signal to noise ratio for the measurement of $C_{\ell }^{\text{gal}}$ in the bin 0.9 <  z <  1.019 as a function of maximum multipole of analysis. Left, from top to bottom: Gaussian covariance only, Gaussian + “other NG”, Gaussian + SSC, Gaussian + SSC + 1h, total covariance. Right: zoom on the three lowest curves: Gaussian + SSC, Gaussian + SSC + 1h, total covariance, all normalised by the value of the signal to noise using the full multipole range and the total covariance.

In the text

	Fig. 7. (Square root of the) cumulative Fisher elements for the cosmological parameters in the considered redshift bin, as a function of the maximum multipole of the analysis.
In the text

	Fig. 8. Marginalised error bars on each cosmological parameters using all redshift bins as a function of the maximum multipole of analysis.
In the text

	Fig. 9. Fisher ellipses on cosmological parameters, using all redshift bins and the full multipole range. The color coding is identical to the other figures: blue=Gaussian, orange=Gaussian+ONG, green=Gaussian+SSC, violet=total covariance.
In the text

	Fig. 10. (Square root of the) cumulative Fisher elements for the HOD parameters in the considered redshift bin as a function of the maximum multipole of analysis.
In the text

	Fig. 11. Marginalised error bars on each HOD parameters using all redshift bins as a function of the maximum multipole of analysis.
In the text

	Fig. 12. Fisher ellipses on HOD parameters, using all redshift bins and the full multipole range. The color coding is identical to the other figures: blue=Gaussian, orange=Gaussian+ONG, green=Gaussian+SSC, violet=total covariance.
In the text

	Fig. A.1. Left: specification for the galaxy redshift distribution from Laureijs et al. (2011) (blue solid) and prediction from the polynomial HOD (red dashed). Right: ratio of the two distributions.
In the text

	Fig. A.2. Galaxy bias predicted from the polynomial HOD.
In the text