The impact of braiding covariance and in-survey covariance on next-generation galaxy surveys

As galaxy surveys improve their precision thanks to lower levels of noise and the push toward small, non-linear scales, the need for accurate covariances beyond the classical Gaussian formula becomes more acute. Here I investigate the analytical implementation and impact of non-Gaussian covariance terms that I had previously derived for the galaxy angular power spectrum. Braiding covariance is such an interesting class of such terms and it gets contributions both from in-survey and super-survey modes, the latter proving difficult to calibrate through simulations. I present an approximation for braiding covariance which speeds up the process of numerical computation. I show that including braiding covariance is a necessary condition for including other non-Gaussian terms, namely the in-survey 2-, 3-, and 4-halo covariance. Indeed these terms yield incorrect covariance matrices with negative eigenvalues if considered on their own. I then move to quantify the impact on parameter constraints, with forecasts for a survey with Euclid-like galaxy density and angular scales. Compared with the Gaussian case, braiding and in-survey covariances significantly increase the error bars on cosmological parameters, in particular by 50% for the dark energy equation of state w. The error bars on the halo occupation distribution (HOD) parameters are also affected between 12% and 39%. Accounting for super-sample covariance (SSC) also increases parameter errors, by 90% for w and between 7% and 64% for HOD. In total, non-Gaussianity increases the error bar on w by 120% (between 15% and 80% for other cosmological parameters) and the error bars on HOD parameters between 17% and 85%. Accounting for the 1-halo trispectrum term on top of SSC, as has been done in some current analyses, is not sufficient for capturing the full nonGaussian impact: braiding and the rest of in-survey covariance have to be accounted for. Finally, I discuss why the inclusion of nonGaussianity generally eases up parameter degeneracies, making cosmological constraints more robust for astrophysical uncertainties. I released publicly the data and a Python notebook reproducing the results and plots of the article.


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 8 = σ 8 (Ω 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 σ 8 .
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 The data and the Python notebook are available at https:// github.com/fabienlacasa/BraidingArticle 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 Hildebrandt 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 autospectra for galaxy clustering in disjoint redshift bins labeled i z , j z , this formula gives where throughout the article I use the short notation 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 treelevel for the prediction of the observable C gal and of the covariance C , .
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 supersample 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 1 .

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 Second, halo model equations can be greatly simplified by introducing the integral where dn h dM is the halo mass function, u(k|M, z) is the normalised halo profile, b β (M, z) is the halo bias or order β 2 , and N (n) gal ≡ N gal (N gal − 1) · · · (N gal − (n − 1)) is the number of n-tuples of galaxies, implicitly depending on halo mass. 1 https://github.com/fabienlacasa/BraidingArticle 2 The terms considered here only involve local bias up to the third order, so β = 0, 1, 2, 3 with b 0 = 1. More generally, we could have non-local bias such as coming from the tidal field at second order: b s 2 .
Finally, further simplifications can be achieved by grouping integrals together: is the sum of second order contributions from perturbation theory and local bias, and is the sum of third order contributions (Lacasa 2018).

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 gal ≡ C (δ 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. 2014aLi et al. ,b, 2018Lacasa & 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 where z a ∈ i z , z b ∈ j z , dV a = r 2 (z a ) dr dz (z a ), and is the SSC kernel, and with the angle-independent trispectrum terms from the halo model, Lacasa (2018) find which can be related to the more usual power spectrum response 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, 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, the 3-halo base term, and the 4-halo term from third order contributions, 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: where is the braiding kernel and encapsulates the non-linear physics.

An approximation to braiding covariance
Directly implementing Eq. (13) for braiding covariance is numerically challenging. Indeed, it would need the computation of B , (z a , z b ) for all pairs of multipoles and all pairs of redshifts. B , (z a , z b ), itself a sum over O( max ) multipoles, quickly makes it a burden for next-gen galaxy surveys where we target max = O(10 3 ).
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 Ψ alt,clust , n gal (z) 2 varies slowly with redshift compared to B , . Then where and with and C n 2 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 Ψ alt 3 and the fact that B 0,0 (z, z ) = σ 2 (z, z ) shows that B , varies quickly enough with redshift for the Bij approximation to work at = = 0. And at higher multipoles, B , only varies more quickly, making the approximation increasingly more precise. Indeed, from Eq. (14), at high ( , ) B , gets contributions from C m a (z, z ) 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.

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. and a Poisson distribution for the satellite galaxies, conditioned to the presence of the central, with mean,

Setup and covariances
In this section, I consider a single redshift bin for the galaxies: 0.9 < z < 1.019. For the HOD parameters, I used log 10 M min = 11.3, σ logM = 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 2 (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 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 .
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 / 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.
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.

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 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 β µ are independent of , . Then we get r 2h1+3 , 2 ≈ I 1 1 I 1 3 P(k ) + ( ↔ ) 2 2I 1 1 I 1 3 P(k ) × 2I 1 1 I 1 3 P(k ) = 1 4 (P(k ) + P(k )) 2 P(k ) P(k ) · (24) Fig. 3. Diagrams for some of the trispectrum terms involved in the covariance of the galaxy angular power spectrum Cov C gal , C 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 gal , while galaxies 3 and 4 are the source of the second power spectrum C gal . 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 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. 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.
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.
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.

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) 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 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.
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.

Impact on parameter constraints
where z 0 = z m / √ 2, with z m = 0.9 the median redshift (Laureijs et al. 2011). The total density is 30 gals·arcmin −2 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 4 . To reproduce this redshift distribution with the halo model, I use the Halo Occupation Fig. 6. Cumulative signal to noise ratio for the measurement of C 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. Distribution described in Sect. 3.1, further including a redshift dependence of M min in the form: 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: and summed over all bins: where α, β are model parameters, that is, cosmological and HOD parameters in the following : (Ω b h 2 , Ω c h 2 , H 0 , n S , σ 8 , w 0 ) and (α sat , σ log M , M ratio , M a min , M b min , M c min , M d min ); ∂ α is the derivative of the observable w.r.t. parameter α.

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.
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: Ω b h 2 , Ω c h 2 , and h, with ONG being more significant than SSC. By contrast, non-Gaussian terms have a large impact on the three last parameters: σ 8 , n S , and w 0 . These three latter parameters are arguably the most interesting to constrain with surveys of the large scale structure. The measurement of σ 8 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 1/2 w,w by a factor 1.8; for other parameters, the factor ranges between 1.08 (for h) and 1.8 (for σ 8 ). Super-sample covariance decreases the information on dark energy by a factor 2.9; for other parameters the factor ranges between 1.01 (for Ω b h 2 ) and 2.6 (for σ 8 ). The total NG decreases F 1/2 w,w by a factor 3.3; for other parameters the factor ranges between 1.08 (for h) and 3.1 (for σ 8 ). When compared to Gaussian+SSC, ONG has a 14% impact on F 1/2 w,w ; for other parameters, the impact ranges between 5.6% (for h) and 16% (for σ 8 ). 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 σ i = (F −1 ) ii for each cosmological parameter as a function of the maximum multipole of analysis max .
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 Gaus-sian+SSC, ONG has a 15% impact on σ w ; for other parameters the impact ranges between 5.7% (for Ω b h 2 ) and 13% (for n S ). The ONG impact exceeds the threshold of 10% (Euclid precision requirement) for two parameters: n S and w (σ 8 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 σ 8 , 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.
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 Ω c h 2 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 Ω b h 2 , Ω c h 2 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 7 in the Gaussian case to 1.4×10 8 in the full non-Gaussian case. This is discussed in more details in Sect. 5.

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.
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 1/2 α sat ,α sat , 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 1/2 α sat ,α sat 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 Gaus-sian+SSC, ONG has a 17% impact on F 1/2 α sat ,α sat ; 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 σ i = (F −1 ) ii for each HOD parameter, as a function of the maximum multipole of analysis max .
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 d min ). SSC increases σ α sat by 9%; for other parameters, the impact ranges between 7% (for σ log M ) and 64% (for M d min ). The total NG increases σ α sat by 24%; for other parameters, the impact ranges between 17% (for M ratio ) and 85% (for M d min ). 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 d min ). The ONG impact is generally stronger than for cosmological parameters, exceeding the threshold of 10%  (Euclid precision requirement) for four parameters: α sat , M b min , M c min and M d min . Interestingly, the impact of ONG is greater than that of SSC for four parameters: α sat , σ log M , M ratio , and M a min . Furthermore. we can note that, for M b min , M c min ,and M d min , 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.
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 7 in the Gaussian case to 3.9 × 10 7 in the full non-Gaussian case.

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 Table 1. For a few cosmological parameters, increase of the error bars compared to the Gaussian case when using the full multipole range and all redshift bins.

SSC+1h
Total NG SSC+1h Total NG 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 5 , 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 9 in the Gaussian case to 6.5 × 10 8 in the full non-Gaussian case 6 . 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 Ω b h 2 , Ω c h 2 , 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 (Ω b h 2 , Ω c h 2 , 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 7 → 1.4 × 10 8 ) and slightly improved by NG for HOD parameters (5.6 × 10 7 → 3.9 × 10 7 ). This means the bulk of the improvement for the whole cosmo+HOD matrix comes from the change in the offdiagonal 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.