Dark Energy Survey Year 3 Results: Covariance Modelling and its Impact on Parameter Estimation and Quality of Fit

We describe and test the fiducial covariance matrix model for the combined 2-point function analysis of the Dark Energy Survey Year 3 (DES-Y3) dataset. Using a variety of new ansatzes for covariance modelling and testing we validate the assumptions and approximations of this model. These include the assumption of a Gaussian likelihood, the trispectrum contribution to the covariance, the impact of evaluating the model at a wrong set of parameters, the impact of masking and survey geometry, deviations from Poissonian shot-noise, galaxy weighting schemes and other, sub-dominant effects. We find that our covariance model is robust and that its approximations have little impact on goodness-of-fit and parameter estimation. The largest impact on best-fit figure-of-merit arises from the so-called $f_{\mathrm{sky}}$ approximation for dealing with finite survey area, which on average increases the $\chi^2$ between maximum posterior model and measurement by $3.7\%$ ($\Delta \chi^2 \approx 18.9$). Standard methods to go beyond this approximation fail for DES-Y3, but we derive an approximate scheme to deal with these features. For parameter estimation, our ignorance of the exact parameters at which to evaluate our covariance model causes the dominant effect. We find that it increases the scatter of maximum posterior values for $\Omega_m$ and $\sigma_8$ by about $3\%$ and for the dark energy equation of state parameter by about $5\%$.


INTRODUCTION
Our understanding of the Universe has become much more accurate in the past decades due to a massive amount of observational data collected through different probes, such as the cosmic microwave background (CMB; see e.g. Planck Collaboration 2020), Big Bang Nucleosynthesis (BBN; see e.g. Fields et al. 2020), type IA supernovae (see e.g. Riess 2017;Smith et al. 2020), number counts of clusters of galaxies (see e.g. Mantz et al. 2014;Costanzi et al. 2019;Abbott et al. 2020), the correlation of galaxy positions, and that of their measured shape (see e.g. Abbott et al. 2018;Heymans et al. 2020). From the study of that data a standard cosmological model has emerged characterized by a small number of parameters (see e.g. Frieman et al. 2008;Peebles 2012;Blandford et al. 2020). Current spectroscopic and photometric surveys of galaxies such as the Extended Baryon Oscillation Spectroscopic Survey (eBOSS 1 ) and earlier phases of the Sloan Digital Sky Survey (SDSS), the Hyper Suprime-Cam Subaru Strategic Program (HSC-SSP 2 ), the Kilo-Degree Survey (KiDS 3 ) and the Dark Energy Survey (DES 4 ) have become instrumental in testing this standard model at a new front: the growth of density perturbations in the late-time Universe. And future surveys, such as the Dark Energy Spectroscopic Instrument (DESI 5 ), the Vera Rubin Observatory Legacy Survey of Space and Time (LSST 6 ), Euclid 7 and the Nancy Grace Roman Space Telescope 8 will push this test to a precision exceeding that provided by other cosmological probes.
An important part of this program is the Dark Energy Survey, a state-of-the-art galaxy survey that completed its six-year observational campaign in January 2019 (Diehl et al. 2019) collecting data on position, color and shape for more than 300 million galaxies. This makes DES the most sensitive and comprehensive photometric galaxy survey ever performed. The main cosmological analyses of the first year (Y1) of DES data have been concluded (Abbott et al. 2018;Abbott et al. 2019b) and analyses of the first three years of data (Y3) are under way. The study of the large-scale structure (LSS) of the Universe based on the DES-Y3 data set has the potential to become the most stringent test of our understanding of cosmological physics to date.
To achieve this goal the DES team is comparing different theoretical models characterized by a range of cosmological parameters to the measured statistics of the LSS in order to determine the model and range of parameters that are in best agreement with the data. The statistics of the LSS considered in the main DES-Y3 analysis are 2-point correlation functions of the galaxy density field (galaxy clustering), the weak gravitational lensing field (cosmic shear) and the cross-correlation functions between these fields (galaxygalaxy lensing) in real space and measured in different redshift bins. These three types of 2-point correlation functions are combined into one data vector -the so-called 3x2pt data vector.
A key ingredient in analyzing these statistics is a model for the likelihood of a cosmological model given the measured correlation functions. Under the assumption of Gaussian statistical uncertainties (which is to be validated) this likelihood is completely characterized by the covariance matrix that describes how correlated the uncertainties of different data points in the 3x2pt data vector are. Validating the quality of the covariance model for the DES-Y3 2-point analyses is the main focus of this paper.
For the DES-Y3 3x2pt analysis we adopt a theoretical covariance model as our fiducial covariance matrix. This fiducial covariance model is based on a halo model and includes a dominant Gaussian component, a non-Gaussian component (trispectrum and super-sample covariance), redshift space distortions, curved sky formalism, finite angular bin width, non-Limber computation for the clustering part, Gaussian shape noise, Poissonian shot noise and f sky approximation to treat the finite DES-Y3 survey footprint (although taking into account the exact survey geometry when computing sampling noise contributions to the covariance). In order to assess the accuracy of that model, we study the impact of several approximations and assumptions that go into it (and into 2-point function covariance models in general): • the Gaussian likelihood assumption, i.e. whether knowledge of the covariance is sufficient to calculate the likelihood; • robustness with respect to the modelling of the non-Gaussian covariance contributions, i.e. contributions from the trispectrum and super sample covariance; • treatment of the fact that 2-point functions are measured in finite angular bins; • cosmology dependence of the covariance model; • random point shot-noise; • the assumption of Poissonian shot-noise; • survey geometry and the f sky approximation; • other covariance modelling details such as flat sky vs.
curved sky calculations, Limber approximation and redshift space distortions.
We generate different types of mock data and/or analytical estimates to determine how each of these effects impacts the quality of the fit between measurements of the 3x2pt data vector and maximum posterior models (quantified by the distribution of χ 2 between the two). We also show how they impact cosmological parameter constraints derived from measurements of the 3x2pt data vector. For most of these tests we employ a linearized Gaussian like-lihood framework which allows us to analytically quantify the impact of covariance errors on the χ 2 distribution and parameter constraints. This is complemented by a set of lognormal simulations and importance sampling techniques to quickly assess large numbers of mock (non-linear) likelihood analyses. This paper is part of a larger release of scientific results from year-3 data of the Dark Energy Survey and our analysis is informed by the (in some cases preliminary) analysis choices of the other DES-Y3 studies. In addition to carving out the most stringent constraints on cosmological parameters from late-time 2-point statistics of galaxy density and cosmic shear yet, the year-3 analysis of the DES collaboration is introducing and testing numerous methodological innovations that pave the way for future experiments. Details of the DESY3 galaxy catalogs and the photometric estimation of their redshift distribution are presented by Sevilla-Noarbe et al. (2020); Hartley et al. (2020); Everett et al. (2020); Myles et al. (2020); Gatti et al. (2020a); Cawthon et al. (2020); Buchs et al. (2019); Cordero et al. (2020). The measurements of galaxy shapes and the calibration of these measurements for the purpose of cosmic weak gravitational lensing analyses are detailed by Gatti et al. (2020b); Jarvis et al. (2020); MacCrann et al. (2020). Krause et al. (2020) develop and test the theoretical modelling pipeline of the DES-Y3 3x2pt analysis, Pandey et al. (2020) outline how galaxy bias is incorporated in this pipeline, DeRose et al. (2020) validate this pipeline with the help of simulated data and Muir et al. (2020) describe how we have blinded our analysis to focus our efforts on model independent validation criteria and reduce the chance for confirmation bias. The DESY3 methodology to sample high-dimensional likelihoods and to characterize external and internal tensions is outlined by Lemos et al. (2020); Doux et al. (2020). Measurements of cosmic shear 2-point correlation functions and analyses thereof are presented by Amon et al. (2020); Secco et al. (2020), the measurement and analysis of galaxy clustering 2-point statistics is carried out by Rodríguez-Monroy et al. (2020) and 2-point cross-correlations between galaxy density and cosmic shear (galaxy-galaxy lensing) are measured and analysed by , with additional analyses of lensing magnification and shear ratios carried out by Elvin-Poole et al. (2020); Sánchez et al. (2020) and results for an alternative lens galaxy sample presented by Porredon et al. (2020, prep). Finally, in DES Collaboration et al. (2020) we present our cosmological analysis of the full 3x2pt data vector.
Our paper is structured as follows. We start by presenting a discussion of our validation strategy in Section 2, where we also summarize our main findings before plunging into the details in the remaining of the paper. In Section 3 we review the modelling and structure of the 3x2pt data vector. Section 4 describes our fiducial covariance model as well as two alternatives to it that are used to validate several modelling assumptions. In Section 5 we describe our linearized likelihood formalism and derive analytically how different covariance matrices impact parameter constraints and maximum posterior χ 2 within that formalism (including the presence of nuisance parameters and allowing for Gaussian priors on these parameters). In Section 6 we present the details of each step in our validation strategy followed by a short Section 7 presenting a simple test to corroborate some of the results from the linearized framework. We conclude with a discussion of our results in Section 8. Seven appendices describe in more detail some results used in this work.

COVARIANCE VALIDATION STRATEGY AND SUMMARY OF THE RESULTS
How should one validate the quality of a covariance model (and the associated likelihood model) for the purpose of constraining cosmological model parameters from a measured statistic? A straightforward answer seems to be that one should run a large number of accurate cosmological simulations, then measure and analyse the statistic at hand in each of the simulated data sets and test whether the true parameters of the simulations are located within the, say, 68.3% quantile of the inferred parameter constraints in 68.3% of the simulations. There are however at least 2 problems with such an approach. The first one is a conceptual problem. Consider a Bayesian analysis of a measured statisticξ with a model ξ [π] that is parametrised by model parameters π. For each value of π the statistical uncertainties in the measurement ξ will have some distribution which is also called the likelihood of the parameters given the data. If this function is known, then a Bayesian analysis will assign a posterior probability distribution to the parameters as p(π|ξ) = 1 N L(π|ξ) pr(π) . (2) Here pr(π) is a prior probability distribution that parametrises prior knowledge from other experiments (or theoretical constraints) and the normalisation constant N is fixed by demanding that p(π|ξ) be a probability distribution. The 68.3% confidence region for the parameters π would then e.g. be stated as a volume V 68.3% in parameter space that contains 68.3% of the probability. To unambiguously define that volume one can e.g. impose the additional condition that min π∈V 68.3% or, more frequently, one would directly define one dimensional intervals that satisfy the above conditions for the marginalised posterior distributions on the individual parameter axes. Unfortunately, if one performs such an analysis many times one is not guaranteed that the true parameters (e.g. of a simulation) are located within V 68.3% in 68.3% of the times. This has recently been referred to as prior volume effect (this issue is discussed in, e.g. Raveri & Hu (2019) and Abbott et al. (2019a)). One may argue that a Bayesian posterior should not be interpreted in terms of frequencies but that doesn't help for the task of validating this posterior on the basis of a large number of simulated data sets. 9 9 In order to deal with the prior volume effect Joachimi et al. (2020) proposed to report parameter constraints through what they call projected joint highest posterior density. This topic will be addressed in a separate DES paper (Raveri et al. -in prep.).
Another, more practical problem is the fact that it is not (yet) feasible to generate enough sufficiently accurate mock data sets to validate covariance matrices of large data vectors with high precision. We recall that for the DES-Y1 analysis, a total of 18 realistic simulated data sets were available to validate the inference pipeline (MacCrann et al. 2018). At the same time, the main reason why N-body simulations would be required to test the accuracy of covariance (and likelihood) models is to capture contributions to the covariance coming from the trispectrum (connected 4point function) of the cosmic density field. But for DES-like analyses it has been shown that this contribution is negligible (see e.g. Krause et al. (2017); Barreira et al. (2018)). The reason for this is twofold: first, very small scales (where the trispectrum contribution to the covariance would matter most) are often cut off from analyses because on these scales already the modelling of the data vector, ξ[π], is inaccurate. And secondly, on small scales the covariance matrix is often dominated by effects coming from sparse sampling such as shot noise and shape noise. These covariance contributions are typically easy to model (although one has to be careful when estimating effective number densities and shape-noise dispersions or when estimating the number of galaxy pairs in the presence of complex survey footprints, see Troxel et al. 2018a;Troxel et al. 2018b).
As a result of the considerations above we base our covariance validation strategy mostly on the use of a linearized likelihood (where the model ξ[π] is linear in the parameters π). In this framework the Bayesian likelihood allows for an interpretation in terms of frequencies -both for total and marginalised constraints. Also, this allows us to perform large numbers of simulated likelihood analyses very efficiently, without the need to run computationally expensive Markov Chain Monte Carlo (MCMC) codes. In addition, any leading order deviation from a linearized likelihood will be next-to-leading order for the purpose of studying the impact of covariance errors (i.e. errors on errors) on our analysis.
Within the linearized likelihood formalism we confirm the findings of Krause et al. (2017); Barreira et al. (2018) for the DES-Y3 setup: both super-sample covariance and trispectrum have a negligible impact on our analysis. This allows us to estimate the impact of other assumptions in our covariance and likelihood model either analytically or by the means of simplified mock data such as lognormal simulations (as opposed to full N-body simulations, cf. Section 4.3).
We summarize our main findings in Figure 1 and Table  1 for the busy reader. For the combined data vector of the DES-Y3 two-point function analysis (the 3x2pt data vector, see details in Section 3) the left panel of Figure 1 shows the impact of different assumptions in our likelihood model on the mean and scatter of χ 2 between maximum posterior model and measurements. To obtain the maximum posterior model we are fitting for all the 28 parameters listed in Table  3 within the linearized likelihood framework described in Section 5.1. Since we assume Gaussian priors on 13 nuisance parameters, the effective number of parameters in that fit will be between 28 and 15. Within the linearized likelihood approach we find that with a perfect covariance model the average χ 2 is expected to be about 507.6, i.e. the effective number of degrees of freedom in the fit is N param,eff ≈ 23.4. The right panel of Figure 1 shows the equivalent results when cosmic shear correlation functions are excluded from the data vector (the 2x2pt data vector). The green points in both panels denote effects that have been already accounted for in the previous year-1 analysis of DES.
What stands out in our analysis is the large effect of finite angular bin sizes on the cosmic variance and mixed terms of our covariance model (cf. Section 4 for this terminology, where we also show that it is unavoidable to take into account finite bin width in the pure shot noise and shape noise terms of the covariance). In DES-Y1 this has been dealt with in an approximate manner, by computing the covariance model for a very fine angular binning and than re-summing the matrix to obtain a coarser binning . This time we incorporate the exact treatment of finite angular bin size for all the three twopoint functions into our fiducial covariance model (cf. Section 4). The blue points in Figure 1 denote improvements that have been made in the year-3 analysis compared to the year-1 covariance model. And the red points are estimates of effects that are not taken into account in the fiducial DES-Y3 likelihood -either because they are negligible, or because an exact treatment is unfeasible (cf. Section 6 for details). Adding these effects in quadrature, our results suggest that the maximum posterior χ 2 of the DES-Y3 3x2pt analysis should be on average ≈ 4% (∆χ 2 ≈ 20.3) higher than expected if the exact covariance matrix of our data vector was known. Table 1 summarizes the offsets in χ 2 displayed in the left panel of Figure 1 and also shows how parameter constraints based on the 3x2pt data vector are impacted by assumptions of our covariance and likelihood model. We distinguish two effects here: the scatter of a maximum posterior parameter π (which we denote by σ[π]) and the width of posterior constraints inferred from our likelihood model (which we denote by σπ). For our tests of likelihood non-Gaussianity we state the changes in the difference between the fiducial parameter values and the upper (high) and lower (low) boundaries of the 68.3% quantile with respect to the standard deviation of the Gaussian likelihood. For our tests of the impact of covariance cosmology, we show the mean of all σπ obtained from our 100 different covariances and also indicate the scatter of these σπ values.
The effect that has the dominant impact on parameter constraints is that of evaluating the covariance model at a 17 set of parameters that do not represent the exact cosmology of the Universe. When computing the covariance at 100 different cosmologies that were randomly drawn from a Monte Carlo Markov chain (run around a fiducial model data vector, see Section 6.8 for details) we find that the differences between these covariances introduce an additional scatter in maximum posterior parameter values. This scatter increases by about 3% for Ωm and σ8 and by about 5% for the dark energy equation of state parameter w. This increased scatter is in fact the dominant effect, since the width of the derived parameter constraints hardly changes between the different covariance matrices. Note especially that re-running the analysis with a covariance updated to the best-fit parameters does not mitigate this effect.
In Figure 2 we take the two effects that had the largest impacts on χ 2 and show the resulting mismatch between scatter of maximum posterior values and width of the inferred contours for a wider range of parameters. All of our results take into account marginalisation over nuisance parameters (and all other parameters).
Our reason for exclusively investigating the impact of covariance errors on χ 2 and parameter constraints is that those are the two measures by which our final (on-shot) data analysis will be interpreted and judged 10 . In the remainder of this paper we detail how the above results were obtained.

THE 3X2-POINT DATA VECTOR
The combined 3x2pt data vector of the DES-Y3 analysis consists of measurements of the following 2-point correlations: • the angular 2-point correlation function w(θ) of galaxy density contrast measured for luminous red galaxies in 5 different redshift bins (see e.g. Rodríguez-Monroy et al. 2020;Cawthon et al. 2020, as well as other relevant references given in Section 1), • the auto-and cross-correlation functions ξ+(θ) and ξ−(θ) between the galaxy shapes of 4 redshift bins of source galaxies (see e.g. Amon et al. 2020;Secco et al. 2020;Myles et al. 2020;Gatti et al. 2020a) , • the tangential shear γ(θ) imprinted on source galaxy shapes around positions of foreground redMaGiC galaxies (see e.g. ).
At the time of writing this paper the exact choices for redshift intervals and angular bins considered for each 2-point function are still being determined by a careful study of their impact on the robustness of DES-Y3 parameter constraints DES Collaboration et al. 2020). For the purposes of testing the modelling of the covariance matrix we will use the most recent but possibly not final DES-Y3 analysis choices. We do not expect that our tests and conclusions will change in a significant manner with further updated analysis choices. We assume that each of the correlation functions are measured in 20 logarithmically spaced angular bins between θmin = 2.5 and θmax = 250 . Some of these bins in some of the measured 2-point functions are being cut from the analysis to ensure unbiased cosmological results, resulting in a total of 531 data points when using the preliminary DES-Y3 scale cuts.
Our starting point of modelling the different 2-point functions in the 3x2pt data vector is the 3D nonlinear matter power spectrum P (k, z) at a given wavenumber k and redshift z. We obtain it by using either of the Boltzmann solvers CLASS 11 or CAMB 12 to calculate the linear power spectrum and the HALOFIT fitting formula (Smith et al. 2003) in its updated version (Takahashi et al. 2012) to turn this into the late time nonlinear power spectrum. From this 3D power spectrum the angular power spectra required for our three 2-point functions (cosmic shear (κκ), galaxy-galaxy lensing (δgκ) and galaxy-galaxy clustering (δgδg)) in the Limber approximation are given by (e.g. Krause et al. 2017; 10 Alternatively, one could investigate the distribution of p-values (or probability to exceed, cf. Hall & Taylor 2019) as opposed to the distribution of χ 2 . 11 www.class-code.net 12 camb.info Limber 1953): where χ is the comoving radial distance, i and j denote different combinations of pairs of redshift bins and the lensing efficiency q i κ and the radial weight function for clustering q i δ are given by Here H0 is the Hubble parameter today, Ωm the ratio of today's matter density to today's critical density of the Universe, a(χ) is the Universe's scale factor at comoving distance χ and b i (k, z) is a scale and redshift dependent galaxy bias. Furthermore, n i κ,g (z) denote the redshift distributions of the different DES-Y3 redshift bins of source and lens galaxies respectively, normalised such that dz n i κ,g (z) = 1 .
Note that on large angular scales the DES-Y3 analysis does not make use of the Limber approximation for galaxy clustering but instead employs the method derived in Fang et al. (2020a). The above angular power spectra are now related to the real space correlation functions w(θ), γt(θ) and ξ±(θ) as Here P are the Legendre polynomials of order , P m are the associated Legendre polynomials, x = cos θ and the functions G +,− ,2 (x) are given in appendix A (see also Stebbins 1996). Note that we only consider the auto-correlations w i (θ) for each tomographic bin since in the Y1 analysis it was shown that the cross correlations do not carry significant information (Elvin-Poole et al. 2018).
The above relations between angular power spectra and real space correlation functions can all be written in the form  This is particularly useful when deriving covariance expressions and when performing averages over finite bins in the angular scale θ. To achieve the latter, one can simply derive analytic averages of the functions F AB (θ). Both of these points will be considered in the next sections.
For most of our tests we consider the 3x2pt data vector and its covariance matrix at the fiducial cosmology described in section 5, where we also show the Gaussian priors assumed on some of these parameters when assessing the impact of covariance modelling on parameter constraints and maximum posterior χ 2 .

COVARIANCE MATRICES FOR THE 3×2PT DATA VECTOR
The covariance matrix of measurements of cosmological 2point statistics typically contains three contributions (cf. Krause & Eifler 2017;Krause et al. 2017), Here, CG is the contribution to the covariance that would be present if the cosmic matter density and cosmic shear fields where pure Gaussian random fields (see also Schneider et al. 2002;Crocce et al. 2011), CnG are contributions involving the connected 4-point function of these fields (the trispectrum) and CSSC is the so-called super-sample covariance contribution resulting from the fact that any survey only observes a finite volume of the Universe and that the mean density in that volume is subject to fluctuations due to long wavelenght modes (Takada & Hu 2013;Schaan et al. 2014).
In the fiducial DES-Y3 analysis we model all of these covariance contributions analytically. This fiducial model is described in Section 4.1. In Section 4.2 we describe an alternative model for the non-Gaussian covariance contributions that is used to test the robustness of our analysis with respect to the modelling of the trispectrum contribution. Finally, Section 4.3 describes a set of log-normal simulations (Xavier et al. 2016) and the covariance matrix of the 3x2pt data vector estimated from them. These simulations also allow us to test the accuracy of our Gaussian likelihood assumption and the treatment of masking and finite survey area in our fiducial covariance model.

Fiducial DES-Y3 Covariance
In our fiducial covariance matrix, we model the non-Gaussian covariance contributions CnG and CSSC using a halo model combined with leading-order perturbation the-ory to approximate the trispectrum of the cosmic density field and to compute the mode coupling between scales larger than the considered survey volume with scales inside that volume. These calculations are carried out using the CosmoCov code package (Fang et al. 2020a) based on the CosmoLike framework (Krause & Eifler 2017). Our modelling of these contributions has not changed with respect to the year-1 analysis of DES and we refer the reader to Krause et al. (2017) as well as to the CosmoLike papers for details. However, the modelling of the Gaussian contribution has changed as described in the following.

Gaussian covariance
Our modelling of the Gaussian covariance part has changed with respect to the year-1 analysis in the following ways:  Figure 2. Impact of covariance errors on the ratio of the standard deviation of maximum posterior parameters to the width of the posterior derived from the erroneous covariance. Green triangles shot the effect caused by non-Poissonian shot-noise and orange circles show the effect caused by the f sky approximation (cf. appendix C for our beyond-f sky treatment). These ratios have been calculated purely on the base of different analytic covariance models and within the linearized likelihood framework discussed in Section 5.1. We also show the ratio of maximum posterior parameter scatter observed from the 197 FLASK simulations to the statistical uncertainties expected from a log-normal covariance matrix matching the FLASK configuration. Within the statistical uncertainties, these ratios are consistent with 1.
• we use (and present for the first time 13 ) analytic expression for the angular bin averaging of the functions F AB (θ) (cf. Equation 10) for all 4 types of two point functions present in our data vector (see Section 6.3, this is especially relevant for the sampling-noise contribution to the covariance, cf. Troxel et al. 2018b); • we account for redshift space distortions (RSD) and also use a non-Limber calculation to obtain the galaxy-galaxy clustering power spectrum C δg δg ( ) (see Section 6.5); • we do not make use of the flat-sky approximation anymore (see Section 6.4).
To derive expressions for the Gaussian covariance part, let us first consider an all-sky survey. If a 2-point function measurementξ AB (θ) could be obtained from data on the entire sky, then for most types of 2-point correlations it would be related to power spectrum measurements C AB from a spherical harmonics decomposition of the same all-sky data through Equation 10, i.e.
A notable exception to this are the cosmic shear 2-point functionsξ± which obtain contributions from both the socalled E-mode and B-mode power spectra (Schneider et al. 2002). For these functions equation (12) in the curved sky 13 We have shared our results with Fang et al. (2020a) who have used them for their covariance calculations.
Since this is a linear equation in C( )'s, the covariance of two different 2-point function measurementsξ AB andξ CD at two different angular scales θ1 and θ2 would be given in terms of the covariance of the corresponding power spectrum measurements by Again, forξ AB (θ) =ξ±(θ) one would have to use For the auto-power spectrum of galaxy density contrast in one of our redshift bins the harmonic space covariance would be ) Here ng is the number density of the galaxies and δ 1 2 is the Kronecker symbol. To account for partial-sky surveys (such as DES) we simply divide this expression (and similar ones for the other 2-point functions) by the observed sky 0.8 z-bin 5 reduced " 2 " = 2.34 p-value = 0.7% Figure 3. Ratio of the diagonal elements of the different covariance matrices introduced in this section with respect to each other. The left panel compares the variances of measurements of ξ + (θ) while the right panel compares the variances of measurements of w(θ). To give a sense of the goodness of fit between the covariance estimated from FLASK and our fiducial analytic matrix, we treat the diagonal elements of the FLASK covariance as a multivariate Gaussian whose covariance can be inferred from the properties of the Wishart distribution (Taylor et al. 2013). The low p-value for the highest redshift bin of w(θ) most likely results from our incomplete treatment of the survey mask (cf. discussion in Section 6 and appendix C).
At this point let us introduce the following nomenclature: we will denote the terms that contain two power spectra as cosmic variance contribution to the covariance, the terms that contain no power spectrum at all as the sampling noise contributions (or shape noise and shot noise contributions) and the terms that contain contribution from one power spectrum and a sampling noise as the mixed terms. We test the accuracy of the f sky -approximation that results in Equations (16-23) in Section 6.6 by comparing it to more accurate expressions.

Analytic lognormal covariance model
To test the robustness of the CosmoLike covariance we also employ an alternative model for the connected 4-point function part of the covariance -the lognormal model. Hilbert et al. (2011) originally derived this as a model for the covariance of cosmic shear correlation function, assuming that that the lensing convergence κ can be written in terms of a Gaussian random field n as (see also Xavier et al. 2016) where it is assumed that n = 0. For given values λ > 0 and µ the power spectrum of n can be chosen such as to reproduce a desired 2-point correlation function ξκ (see Xavier et al. (2016) for caveats). Furthermore, for any given value λ > 0 one can choose µ such that κ = 0. This makes λ the only free parameter of the lognormal covariance model. Hilbert et al. (2011) show that this model leads to a number of correction terms to the Gaussian covariance model, and identify the most dominant of these terms to be Here AS is the area of the considered survey footprint and VarS(κ) is the variance of κ when averaged over the footprint. We generalise this to the covariance of 2-point correlationsξAB andξCD between arbitrary scalar fields δA, δB, δC , δD as Here, CovS(δA, δC ) is the covariance of δA and δC after the two fields have been averaged over the entire survey footprint (and likewise for the other terms appearing above). Following Hilbert et al. (2011) we use this expression even when considering non-scalar fields (i.e. the shear field) by replacing ξXY (θ) by the appropriate 2-point functions ξ+(θ), ξ−(θ), γt(θ) (or w(θ), for the scalar galaxy density contrast).
To choose the parameters λX (also called the lognormal shift parameters, cf. Xavier et al. 2016) we follow a procedure similar to the one outlined in Friedrich et al. (2018). There it is shown how the value of λX can be adjusted in order to match the re-scaled cumulant n(z) source galaxies lens galaxies Figure 4. Redshift distributions of lens galaxies (shaded regions) and source galaxies (solid lines) in our fiducial test configuration.
of the random field δX smoothed with a top-hat filter of angular radius ϑ to the value of S3 predicted by leadingorder perturbation theory for that same smoothing scale.
Since the focus in our paper is the covariance matrix of 2point statistics (hence a 4-point function), we modify their method to match instead the value of reduced fourth order cumulant The field δX here will be either projections of the 3D matter density contrast along the line-of-sight distribution of our lens galaxies or the lensing convergence fields corresponding to our 4 source redshift bins. The smoothing scale ϑ at which we use the λX to match S4 to its perturbation theory value is chosen such that it corresponds to about 10Mpc/h at the mean redshift of the line-of-sight projection kernels corresponding to the different δX . This is approximately the scale at which Friedrich et al. (2018) found the shifted lognormal model to be a good approximation of the overall PDF of density fluctuations in N-body simulations (cf. their figure 5).
Our results are shown in Table 2 , where we present the number density, galaxy bias (relevant for lenses only), shape-noise dispersion (per shear component; relevant for sources only) and the lognormal shift parameters obtained from the procedure described above. Note that for the source galaxy samples, the relevant line-of-sight projection kernel used to derive the shift parameter is the lensing kernel (and not the redshift distribution of the source galaxies). For the lens galaxies, all shift parameters come out to be > 1. As a consequence there will be pixels with negative density in our lognormal simulations. However, the fraction of such pixels is < 0.01 for all runs and all bins and setting δg = −1 in these bins has an unnoticeable effect on the statistics measured in these maps (e.g. for bin 4, which is affected most, the standard deviation of δg changes by 0.053%). Note further that at the time of completing the simulation runs presented in Section 4.3, the DES Y3 shear catalog and redshift distribution were not finalized. As a consequence, the shape noise dispersion values used for simulations differ from the values in this table.

Lognormal covariance from simulations
We also produce a test DES-Y3 covariance matrix from a set of simulations. We use the publicly available code FLASK (Full sky Lognormal Astro fields Simulation Kit) (Xavier et al. 2016) 14 to generate 800 DES-Y3 footprint sky maps of density, convergence and shear healpix maps (Górski et al. 2005) with NSIDE=8192, as well as galaxy positions catalogs, used to reproduce the DES-Y3 properties. FLASK is able to quickly produce tomographic correlated simulations of clustering and weak lensing lognormal fields based on the DES-Y3 lens and sources samples. The lognormal distribution of cosmological fields has been shown to be a good approximation (Coles & Jones 1991;Wild et al. 2005;Clerkin et al. 2017) but much less computationally expensive to generate than full N-body simulations.
As input for the simulations, we used a set of auto and cross correlated power spectrum and the lognormal field shift parameters. The theoretical input power spectrum was generated using CosmoLike, and the lognormal 14 http://www.astro.iag.usp.br/ flask/ shifts are the ones listed in Table 2. In order to reproduce the properties of shear fields, we added the shapenoise term by sampling each pixel of the simulated maps to match the correspondent shape-noise dispersion σ and number density ng of the tomographic bin. At the time of completing the simulation runs, the DES Y3 shear catalog and redshift distribution were not finalized. For this reason, the values used in the simulations are slightly diffeent from the values in Table 2. For the simulations, we set the number density for the five tomographic lens bins as 0.0227, 0.0392, 0.0583, 0.0451, 0.0278 (arcmin −2 ). The shape-noise dispersion values for the four tomographic bins of sources were set to 0.27049, 0.33212, 0.32537, 0.35037. The cosmology adopted for the theoretical power spectra is set as Ωm = 0.3, σ8 = 0.82355, ns = 0.97, Ω b = 0.048, h0 = 0.69, and Ων h 2 0 = 0.00083. We use the publicly available code TreeCorr 15 (Jarvis et al. 2004) to measure the 3x2 point correlation measurements for 200 DES-Y3 realizations. For all measurements, we used 20 log-spaced angular separation bins on scales between 2.5 and 250 arcmin. We set the bin_slop TreeCorr parameter to zero, essentially setting all estimators to bruteforce computation. In Figure 5 we show the validation of the measurements comparing with the theoretical input.
We will use the FLASK covariance mainly to estimate the impact of the survey geometry.

Comparisons among covariances
Here we present some comparisons between the different covariance matrices. In Figure 3 we show the ratio of the diag-onal elements of the different covariance matrices introduced in this section displaying both the variances of the measurements of ξ+(θ) of w(θ).
In Figure 6 we compare the covariance matrices obtained from the FLASK simulations and the analytical halo model covariance.

IMPACT OF COVARIANCE ERRORS ON A LINEARIZED GAUSSIAN LIKELIHOOD
As discussed above a full assessment of the impact of using different covariance matrices to parameter estimation becomes unfeasible due to the computational demand of running a large number of MCMC chains. Since the covariance matrices studied in this work differ by subdominant effects we do not expect large modifications in the results of the estimation of the parameters. Therefore we will bypass this difficulty by using a linearized approximation of the model data vector as a function of the parameters. The measured data is assumed to be a Gaussian multivariate variable characterized by a covariance matrix and a given prior matrix. This approach is called the Gaussian linear model (Seehars et al. 2014(Seehars et al. , 2016Raveri & Hu 2019). Within this approach we study the following impacts of different covariances: • error in the parameter estimation, characterized by the width of the contours; • the scatter of the best fit (maximum posteriors) parameters; • change in the maximum posterior χ 2 value; • error in the maximum posterior χ 2 value.
In the remainder of this section we detail this method.

Linearized likelihoods
To speed up our simulated likelihood analyses, we employ a linearized model of the data vector ξ (e.g. the DES-Y3 3x2-point function data vector). This can be considered a linear Taylor expansion of our full model around a fiducial set of parameters π 0 which is summarized in Table 3. In this approximation our model data vector becomes where the sum is over all components πα of the parameter vector π (we will use latin indices for the components of the data vector and greek indices for the components of the parameter vector). Given a 2-point function measurementξ and abbreviating Table 3. Fiducial cosmology and standard deviation of Gaussian parameter priors used in our mock likelihood analyses. A IA,i is the intrinsic alignment amplitude in the ith source redshift bin, m i is the multiplicative shear bias and ∆z s,i parametrizes systematic shifts in the photometric redshift distribution of that bin. ∆z l,i parametrizes systematic shifts in the photometric redshift distribution of the ith lens redshift bin. The Gaussian priors we choose for the parameters follow the analysis choices of Abbott et al. (2018) and we assume infinite flat priors for all other parameters.
Parameter Fiducial value σ prior our figure of merit χ 2 as a function of the parameters becomes in the linearized approximation Here we have allowed for a Gaussian prior with covariance matrix P −1 and central value π prior . To find the deviation δπ MP = π MP − π 0 from our fiducial parameters that minimizes this function (the maximum posterior value of the parameters is denoted by π MP ) we have to solve Defining a vector x such that x β = δξ T C −1 ∂ β ξ as well as the Fisher matrix F αβ = ∂ β ξ T C −1 ∂αξ this becomes We now want to consider the situation when a model covariance matrix C mod is used to calculate the likelihood in equation (30 which is different from the true covariance matrix Ctrue of the statistical uncertainties in the data vector ξ. In that case our linearized likelihood will be a Gaussian centered around π MP and with parameter covariance matrix where F mod,αβ = ∂ β ξ T C −1 mod ∂αξ is the Fisher matrix calculated from the model covariance. The actual covariance matrix of π MP includes two sources of noise. First, statistical uncertainties in the measurementξ which are described by the covariance matrix Ctrue and are represented by the first term in equation (32) that is proportional to x. And secondly, statistical uncertainties in our choice of the prior center which are described by the prior covariance matrix P −1 and are represented by the second term in equation (32) that is proportional to π prior . The latter term has the covariance matrix (F mod + P) −1 P(F mod + P) −1 (because the covariance matrix of π prior is P −1 ). Hence, the total covariance matrix of π MP can be written as For C mod = Ctrue it is easy to see that this parameter covariance Cπ,MP equals the covariance C π,like that describes the shape of our likelihood (as it should).

Impact on the width of the likelihood and scatter of best fit parameters
We can use the above findings to study the impact of different effects in covariance modelling on parameter constraints. If a covariance matrix C1 contains a noise contribution that is missing in another covariance matrix C2, then we quantify the difference between these matrices by considering two effects: • Width of likelihood contours: Denoting the Fisher matrices obtained from C1 or C2 as F1 and F2 respectively, the width of likelihood contours drawn from the different covariances are given by Hence, if the difference C1 − C2 = E represents noise contributions missing from (or miss-estimated in C2), then a comparison of C π,like, 1 and C π,like, 2 quantifies the impact of this on the width of parameter contours.
• Scatter in the center of likelihood contours: If the data vectorξ had C1 as its true covariance matrix but C2 would be used to derive the maximum posterior parameters π MP from it, then the maximum posterior parameter covariance would be given by If the difference C1 − C2 = E represents noise contributions missing from (or miss-estimated in C2), then a comparison of Cπ,MP, 2 and Cπ,MP, 1 ≡ C π,like, 1 quantifies the impact of this on the scatter in the location of parameter contours.
An inaccurate covariance model will in general have a different impact on the width and the location of parameter contours. Hence, in order to quantify the importance of different effects in covariance modelling for parameter estimation, we compare both the pair C π,like, 1 / C π,like, 2 and the pair Cπ,MP, 1 / Cπ,MP, 2.

Distribution of χ 2 when fitting for parameters
Within the linearized likelihood model developed in the previous section we now investigate how errors in the covariance model impact the distribution of χ 2 MP between measured data vectorξ and a maximum posterior model ξ MP = ξ(π MP ), We start with the case that (i) the true covariance C ofξ is known (ii) no parameter priors are used when determining the best fitting model ξ MP (iii) the true expectation valueξ ≡ ξ lies within our parameter space. I.e. there are parameters π true such that ξ(π true ) =ξ .
We will show that, as expected, in this caseχ 2 MP should follow a χ 2 -distribution with N data − Nparam degrees of freedom.
Using equations (29) and (32) (and setting again δξ ≡ ξ − ξ 0 ) one can see that the maximum posterior data vector is given by Here, the second line follows from the fact thatξ = ξ = ξ MP and we have defined the matrix It can be shown that P is an idempotent matrix (P 2 = P) and furthermore that The residual between the measurementξ and the best fitting model ξ MP can be written in terms of P aŝ Hence, the covariance matrix ofξ − ξ MP is given by This makes it straightforward to find the expectation value Similarly, the variance of χ 2 MP can be shown to be So far, we have only re-derived textbook results (Anderson 2003). Now how do χ 2 MP and Var(χ 2 MP ) change if the covariance model C mod we use to find the best fitting model ξ MP and to compute χ 2 MP is different from the true covariance matrix C ofξ?
Following similar steps as Eqs. (38) and (39) one can show that where and where the Fisher matrix F mod is computed from the model covariance C mod . Equation 45 especially shows that ξ MP is still an unbiased estimator ofξ even when C mod = C.
When deriving the moments of χ 2 MP we will still come across expectation values like (cf. Equation 43) Hence the expectation value and variance of χ 2 MP are given by where Now we are left to investigate how Equations 48 and 49 change when a Gaussian parameter prior P is included in the likelihood function (cf. Equation 30). A complication in this case is, that now ξ MP is not necessarily an unbiased estimate ofξ anymore. This is because in Equation 30 we have centered our prior around the model parameters π prior which may be different from the true parameters π true . Inserting the full expression for the maximum posterior parameters (Equation 32) into our linearized model we now get The residual betweenξ and ξ MP hence becomeŝ Treating the prior center π prior again as a random vector centered around π true , ζ also becomes a random vector with covariance Hence, along lines similar to the case without a prior, we can write the moments of χ 2 MP for a given model covariance as Notice that in the absence of priors C ζ = 0 and for the true covariance C we recover equations (43) and (44) as expected.
Equations (55) and (56) are used to produce our main result shown in Figure 1 for different covariance matrices.

EXPLORING DIFFERENT EFFECTS IN THE COVARIANCE MODELLING
Our main goal is to study the impact of including different effects in the covariance modelling on the estimation of parameters. Several covariance matrices were generated and tested under different assumptions and approximations. The main results were already shown in Section 2. We now present the details of each step in the validation strategy that was outlined in Section 5.

Gaussian likelihood assumption
A basic assumption of our framework of testing different covariance matrices is that the likelihood function of the data is Gaussian. One simple reason of why the sampling distribution of the correlation functions can not be an exact multivariate Gaussian is that this violates the positivity constraint of the power spectrum ). and when assuming that the skewness of the data points is 5 times that of the f sky approximation (green histogram). Bottom panel: Distribution of maximum posterior σ 8 when fitting the linearized model to Section 5.1 Gaussian realisations of our fiducial data vector , to lognormal realisations of our fiducial data vector (blue histogram) and to lognormal realisations with 5 times the skewness of the f sky approximation employed on Section 6.1 (orange histogram).
There are also other reasons described below. The purpose of this Subsection is to assess the impact of non-Gaussianity of the likelihood of 2-point functions in the parameter estimation. In this sense checking this basic assumption is a test of the whole framework and is different from the robustness tests for the covariance matrix modelling described in the remaining Subsections of this Section. The impact of a non-Gaussian likelihood in parameter estimation of weak lensing correlation functions has been recently studied in Lin et al. (2020) where no significant biases were found in one-dimensional posteriors of Ωm and σ8 between the multivariate Gaussian likelihood model and more complex non-Gaussian likelihood models. Also in Sellentin et al. (2018) the skewed distributions of weak lensing shear correlation functions are used to derive an analytical expression for a non-Gaussian likelihood.
We first consider a full-sky survey such that each of our 2-point function estimatorsξ AB (θ) is a harmonic transform of a harmonic space estimatorĈ AB (cf. Equation 12), i.e.
EachĈ AB is given in terms of the spherical harmonics coefficients a m , b m of two Gaussian random fields aŝ The product of two Gaussian random variables does not follow a Gaussian distribution. Therefore, in principle one would not expectĈ AB (and consequentlyξ AB (θ)) to have a Gaussian likelihood. However, at small scales, i.e. at high multipoles , the sum of the random variables a m b * m in Equation 58 will approach a Gaussian distribution by means of the central limit theorem, since there is a large number (2 + 1) of independent modes. It should be pointed out that at these small scales the galaxy density and shear fields characterized by a m and b m are themselves non-Gaussian due to the non-linear evolution of gravity.
It is hence our working hypothesis that non-Gaussianity ofĈ AB only matters at the largest scales (small 's) where both a m and b m can be considered Gaussian random variables but not their product. In the full-sky case it can then be shown that the second and third central moments ofĈ AB are given by If only a fraction f sky of the sky is being observed, these moments get divided by f sky and f 2 sky respectively. Assuming different multipoles to be uncorrelated, the corresponding moments ofξ AB (θ) can be computed as Equation 61 is of course nothing but the diagonal of the covariance matrix (cf. Equation 14).
The dominant effect of the non-Gaussianity of the C 's is a positive skewness in the distribution of our data vectors (Sellentin et al. 2018). To estimate its impact on our parameter constraints, we approximate the entire distribution of our 3x2pt data vector by a multivariate lognormal distribution. The covariance of our data vector and the skewness of each data point as given by Equation 62 are sufficient to fix the parameters of a shifted log-normal distribution. We have already discussed this in Sections 4.2 and 4.3, though with a conceptual difference: in that section we describe how to configure log-normal simulations of the cosmic density field, while here we assume measurements of the 3x2-point functions to have a multivariate log-normal distribution. To be explicit, we fix the shift parameters λ(θ) that enter the lognormal PDF of the measurementsξ AB (θ) in the different angular bins (cf. Equation 24 for the definition of λ) via the equation which relates the 2nd and 3rd central moments of log-normal random variables (Hilbert et al. 2011).
In the top panel of Figure 7 we show the impact of this non-Gaussianity on the distribution of maximum posterior χ 2 . For that figure we generated 300, 000 random realisations of our fiducial data vector from a multi-variate Gaussian distribution, 300, 000 random realisations of that data vector from a multi-variate lognormal distribution and 300, 000 random realisations from another lognormal distribution, whose skewness in each data point was increased by a factor of 5. For each of these random realisations we analytically determined the maximum posterior model within the linearized likelihood formalism of Section 5.1 and then computed the χ 2 between that model and the random realisation. The blue histogram in the top panel of Figure  7 shows the distribution of these χ 2 values for the Gaussian random realisations and the red histogram corresponds to the log-normal random realisations. The two histograms are almost identical. Hence, within the f sky -approximation employed above non-Gaussianity in the likelihood does not seem to affect our analysis. And even in the extreme scenario of enhancing the skewness of the data vector by a factor of 5 (green histogram) the increase in the scatter of χ 2 remains smaller than about 3% of the average χ 2 -which still wouldn't dominate over the other effects discussed in subsequent sections (cf. Figure 1).
The impact of non-Gaussianity on the likelihood becomes even more negligible when directly considering the distribution of maximum posterior parameters. We demonstrate this in the bottom panel of Figure 7 for the best-fit values of σ8 but find similar results for our other key cosmological parameters Ωm and w0. Therefore, we conclude that it is safe to assume a Gaussian distribution for the statistical uncertainties of the DES-Y3 2-point function measurements.

Modelling of connected 4-point function in covariance
The connected 4-point contribution to the covariance is the part that is most challenging to model analytically (Schneider et al. 2002;Hilbert et al. 2011;Sato et al. 2011;Takada & Hu 2013). This contribution is most relevant at small scales and turns out to be a small one for current large-scale structure analyses ; Barreira et al. 2018). This is for two reasons: 1) such analyses typically cut away their smallest scales because of uncertainties in the modelling of their data vectors and 2) at small scales the covariance matrix is often dominated by shape noise and shot noise which are believed to be well understood. We test whether the non-Gaussian covariance parts (by which we mean both the connected 4-point function and super-sample covariance) are a relevant contribution to our error budget by either • replacing the non-Gaussian contributions from the fiducial halo model with the lognormal covariance described in Section 4.2 • or setting it to zero, i.e. using a Gaussian covariance matrix only. Figure 1 and Table 1 show that neither of these changes has a significant impact on the distribution of χ 2 and our parameter constraints. Assuming that our halo model and lognormal recipes do not underestimate the non-Gaussian covariance parts by orders of magnitude (see e.g. Sato et al. (2009);Hilbert et al. (2011) for justifications of this assumption) this demonstrates that we are insensitive to the exact modelling of these contributions. At the same time, we want to stress that this finding holds for the specific scale cuts, redshift distributions and tracer densities of the DESY3 3x2pt analysis and cannot necessarily be generalized to other analysis setups.

Exact angular bin averaging
Equation 14 holds when measuring the 2-point correlation functions in infinitesimally small bins around the angular scales θ1 and θ2. This is unfeasible in practice and in fact also leads to divergent covariance matrices. This can for example be seen for the galaxy clustering correlation functions, where the constant term proportional to 1/n 2 g in the harmonic space covariance gives a contribution to the real space covariance of The reason for this divergence is simply the fact that the number of galaxy pairs found in an infinitesimal bin vanishes, leading to infinite shot-noise. This problem disappears when considering finite angular bins.
To analytically average over a finite angular bin [θmin, θmax], we assume that the number of galaxy pairs with angular separation θ is proportional to sin θ (corresponding to a uniform distribution of galaxies on the sky). We then replace the functions F AB (θ) in Equations 9 and 10 by For the galaxy clustering correlation function w(θ) this leads to (2 + 1)(cos θmin − cos θmax) .
The corresponding expressions for the galaxy-galaxy lensing correlation function γt(θ) and for the cosmic shear correlation functions ξ± are presented (together with derivations of all the bin averaged expressions) in appendix B.
We show below how the bin averaging solves the problem of diverging diagonal values of the covariance for w(θ). This can be seen from where Asurvey = 4πf sky is the total survey area, A bin = 2π (cos θmin − cos θmax) is the bin area and Npair the total number of galaxy pairs used to estimateŵ. The above expression is the usual formula for the shot-noise part of the real space covariance. The impact of the exact angular bin averaging for the noise and mixed terms in the Gaussian part of the covariance matrix is included for all 4 types of two point functions present in the DES-Y3 data vector and the DES-Y3 fiducial covariance.

Flat vs. Curved sky
For the Y1 analysis it was shown that the flat-sky approximation was valid for the galaxy-galaxy shear and shear-shear 2-point correlation function . In Y3 the fiducial covariance computes the full sky correlations, see equations (12) and (13). We show in Fig. (1) that the effect of including curved sky results has negligible impact on the χ 2 distribution. Table 1 shows that this is also true for parameter constraints.

RSD and Limber approximation and redshift space distortion effects
The modelling of the angular power spectrum of two tracers involves a projection from the three dimensional power spectrum that requires integrals with integrands containing the product of two spherical Bessel functions, which are highly oscillatory. The inclusion of redshift space distortion (RSD) effects in a simple linear modelling (Kaiser 1987) involves the computation of those integrals with derivatives of the Bessel functions. These integrals are notoriously difficult to perform numerically and it is usual to apply the so-called Limber approximation (Limber 1953;LoVerde & Afshordi 2008). An efficient computation of these integrals without resorting to the Limber approximation was recently implemented in the case of the angular power spectrum for galaxy clustering in Fang et al. (2020b). We use their approach to study the impact of taking into account both non-Limber computations and RSD effects in the covariance matrix. Figure 1 and Table 1 show that not taking these effects into account leads to an increase in average χ 2 of about 0.5% and an underestimation of uncertainties in key cosmological parameters by 0.6% to 1.4%.

Effect of the mask geometry
The analytical covariance models described in Section 4 make use of the so called f sky approximation, i.e. they take the covariance of an all-sky survey and divide this by the sky fraction of DES-Y3 to approximate the covariance of our partial sky data. In appendix C we show how to go beyond this approximation. First, we note there that the covariance of the 2-point function measurements between pairs of scalar random fields (δa, δ b ) and (δc, δ d ) within angular bins [θ ab − , θ ab + ] and [θ cd − , θ cd + ] is given by Here P are again the Legendre polynomials and the angular bin averaging was already evaluated. The factor N ab is the average number of pairs of random points that uniformly sample the footprint with densities na, n b (resp. nc, n d ) and whose separation falls into the angular bin [θ ab − , θ ab + ] (resp. [θ cd − , θ cd + ]). Hence, these factors describe how the exact survey geometry suppresses the number of pairs of positions in the bins [θ ab − , θ ab + ] and [θ cd − , θ cd + ] with respect to the f sky approximation. And finally, Cov{C ab 1 ,C cd 2 } is the covariance of pseudo-C estimates of the power spectra between the fields (δa, δ b ) and (δc, δ d ) (see appendix C for more details). Note that the full survey footprint modifies the covariance with respect to the f sky approximation used in Section 3 both through the factors N ab pair [θ ab − , θ ab + ]/nan b , N cd pair [θ cd − , θ cd + ]/ncn d and by changing Cov{Ĉ ab 1 ,Ĉ cd 2 } compared to Equations 16-23. One can determine the factors N ab pair [θ ab − , θ ab + ]/nan b and N cd pair [θ cd − , θ cd + ]/ncn d either by counting pairs in a set of random points that trace the survey footprint homogeneously or they can be calculated analytically from the power spectrum of the survey mask (see our appendix C as well as Troxel et al. 2018b). This will generally lead to an enhancement of statistical uncertainties (i.e. of the covariance matrix) with respect to the f sky approximation. To calculate Cov{C ab 1 ,C cd 2 } one could e.g. follow approximations made by Efstathiou (2004). We slightly modify their arguments in appendix C to arrive at (2 1 + 1)(2 2 + 1) + + C ad 1 C bc 2 + C ad 2 C bc 1 + C ad 1 C bc 1 + C ad 2 C bc 2 (2 1 + 1)(2 2 + 1) Here the mode coupling matrix M 1 2 again depends on the power spectrum of the survey mask and is also detailed in appendix C. Note that in order to keep our notation brief, we have assumed that the power spectra C ac 1 etc. in the above equation include both the underlying cosmological power spectra and contributions to the power spectra from sampling noise, such as shape-noise and shot-noise. In practice, Equation 68 and the approximations proposed by Efstathiou (2004) yield very similar results and they are both valid on scales 1, 2 which are much smaller than the typical scales of the mask W . Unfortunately, the DES-Y3 analysis mask has features and holes over a large range of scales. Hence, the angular scales of interest in the 3x2pt analysis are never strictly smaller than the scales of our mask. Hence, Equation 68 is not sufficiently accurate in our case and in fact significantly overestimates our statistical uncertainties. In Figure 8 we explain a simple scheme that can be used to correct for this. To motivate this procedure, consider how one would compute the Gaussian covariance model directly from the real space 2-point correlation functions, i.e. without taking the detour to Fourier space that was used in Section 3.
Here Ω a . . . Ω d are fours locations inside the survey mask such that the distance between Ω a and Ω b lies inside the angular bin [θ ab − , θ ab + ] and the distance between Ω c and Ω d lies inside the angular bin [θ cd − , θ cd + ]. Now the approximation of Efstathiou (2004) assumes that the 2-point functions ξ ac (θ), ξ bd (θ) are negligible on scales θ comparable to the smalles features in the mask. Schematically this amounts to making the approximation dΩ a . . . W (Ω a )ξ ac (θ ac ) ≈ξ ac dΩ a . . . W (Ω a )δ 2 Dirac (Ω a − Ω c ) whereξ ac is a suitable average of the 2-point function over different scales. Our understanding of the approximation of Efstathiou (2004) via Equation 70 is based on findings that we present in appendix D. This approximation fails when the mask contains features (e.g. holes) on scales where the 2-point function has not yet decayed. Assuming that such small scales holes cover a fraction of f mask of a more coarse version of the footprint, then this can roughly be corrected for with a multiplicative factor, i.e. by instead using the approximation dΩ a . . . W (Ω a )ξ ac (θ ac ) The right panel of Figure 8 visualizes this for the mixed terms in the covariance, where one of the correlation functions ξ ac or ξ bd is due to sampling noise such as shape-noise or shot-noise and is hence exactly proportional to a Dirac delta function. In that case, the integration is over pairs that share one end point. Now the approximation made e.g. in Efstathiou (2004) or by our Equation 68 assumes that also the correlation function between the other two end points effectively acts as a delta function with respect to the smallest scale features in the survey mask (cf. Equation 70). We find that this is not the case for the DES-Y3 mask and that it contains features on all scales relevant to our analysis. But as indicated in Equation 71, one can approximately correct for this by multiplying the mixed terms in the covariance by the fraction f mask of the coarser survey geometry that is covered by small scale holes in the mask. This can be considered a next-to-leading order correction to our Equation 68.
By applying Equation 71 twice one can see that the cosmic variance terms (terms where neither of the 2-point functions ξ ac or ξ bd are exactly proportional to delta functions) can be corrected by multiplication with f 2 mask . To implement this correction in practice we draw circles within the DES-Y3 survey footprint with radii ranging from 5arcmin to 20arcmin and measure the masking fraction in these circles. We find that this fraction is ≈ 90% across the considered scales. Multiplying the mixed terms in the covariance by that fraction and the cosmic variance terms by the square of that fraction (together with using Equation 68) we indeed find significant improvement of the maximum posterior χ 2 obtained for the FLASK simulations -as is shown in the left panel of Figure 8 (as well as in Figure 1).
In Figure 9 we use our FLASK measurements together with the technique of precision matrix expansion (PME, from inverse covariance = precision matrix Friedrich & Eifler 2018) and perform a consistency of the modelling ansatz described above by investigating the impact of masking on individual covariance terms. We find both with the PME methods and with our analytic ansatz that masking effects are most impactful in the covariance terms that depend on shape-noise of the weak lensing source galaxies (i.e. in what we called mixed terms in Section 3). This also agrees with the findings of Joachimi et al. (2020) and it further motivates the modelling of masking effects that we have described here. Nevertheless, we do not elevate this modelling ansatz to our fiducial covariance model because its motivation remains rather heuristic. But we consider it a realistic estimate for the error made by the f sky approximation and can hence use it to estimate the impact of that approximation on parameter constraints. In Figure 2 we have already shown that this impact below the 1% level, i.e. we underestimate the scatter of maximum posterior parameters by less than 1% when making the f sky approximation in our fiducial covariance model.
Note that Kilbinger & Schneider (2004); Sato et al. . Given the large area of DES-Y3 and and its numerous combinations of redshift bins, we did not find this to be feasible.

Non-Poissonian shot noise
In the Poissonian limit and in a complete region of the sky the power spectrum of shot-noise is scale-independent and Small scale holes in survey mask Figure 8. The impact of masking on the DES-Y3 covariance. The blue histogram in the right panel shows the distribution of χ 2 obtained from our FLASK data vectors when using the f sky approximation. We restrict this figure to the 2x2pt function part of the data vector since it is this part that suffers the most from masking effects (cf. Figure 1). Ansatzes in the CMB literature (e.g. Efstathiou 2004) are not sufficient to correct for this, because the DES footprint has features down to very small scales. In the main text we have motivated a possible way to correct for these small scale masking features and the orange histogram in the left panel shows that this ansatz indeed significantly improves the χ 2 obtained from our FLASK measurements. The sketch in the right panel visualises how small scale features in the mask lead to an overestimation in the covariance when using common ways to treat the impact of survey geometry on the 2-point function covariance (see main text for explanation).
given by wheren is the galaxy density per steradian. As noted in the previous subsection, in galaxy surveys not every region of the sky is fully accessible, i.e., the presence of bright stars, satellite trails, etc. lead to artificial changes in the measured galaxy density. These density changes can potentially modify the observed galaxy power spectrum and bias any cosmological analyses derived from them, and thus, they are avoided by removing certain regions of the sky where artifacts may be found. These regions are usually smaller than the resolution of the (pixelated) survey mask used to determine whether a region of the sky is within the footprint or not, since it is computationally expensive to increase the resolution. This this can be described through a fractional mask Wi = 1/fi, where i is a given pixel of the mask fi is the fractional area of the pixel unaffected by the presence of artifacts. If we compute the galaxy overdensity as wherew is the mean of the weights wi across the footprint, and Ωpix is the area of the pixels from the mask in steradians. In the case where all the pixels in the footprint are fully complete we recover Equation (72) sincen =N /Ωpix, andN = i Ni/ i 1. However, in the case that any of the pixels of the mask are not fully complete we obtain an increased shot-noise contribution by a factorW 1 (since 0 fi 1).
In previous studies the DES-Y1 lens galaxies were shown to prefer a super-Poissonian variance Gruen et al. 2018) which might be a consequence of their complex selection criteria or due to the nature of their formation and evolution (see e.g. Baldauf et al. 2013;Dvornik et al. 2018). This super-Poissonian variance leads to an enhance in shot-noise. In order to test for this effect, we proceeded to estimate the angular power spectrum, C ≈ C ,galaxies + N + δC , of DES-Y1 redmagic galaxies selected in Elvin-Poole et al. (2018) using NaMaster (Alonso et al. 2019), where N is the shot noise contribution from equation (73), and δC is the excess power which can be due to a number of factors (variations in completeness not captured by the mask, super-Poissonian shot noise, observational systematics, etc.). We also compute the power spectrum, C ,rnd of a random field with the same number of objects as the galaxy sample considered, and the probability of populating a pixel i is proportional to its weight, Wi.
We find that C ,rnd is statistically consistent with N . We then compute the ratio: In Figure 10 we show r compared to the theoretical expectation for C ,galaxies /N = C ,th /N , where C ,th is the theoretical power spectrum computed using the best-fit parameters found in Elvin-Poole et al. (2018). We also allow for a 20% variation in the linear galaxy bias, which is much larger than the uncertainty found in Elvin-Poole et al. (2018). We find that there is an excess power at 3000 that cannot be explained by an excess galaxy clustering (i.e., a larger than measured linear bias or a non-linear bias component). We identify this excess (between 2% and 6%) with δC N in equation (74) number of data points 2 of FLASK sims using different covariances (no fitting!) estimate using 1st order PME analytic estimate Figure 9. The method of precision matrix expansion (Friedrich & Eifler 2018) allows us to estimate the impact of individual covariance terms on χ 2 even when only few simulated measurements are available. The orange squares show the average χ 2 between our FLASK measurements and their mean (rescaled by a factor of N FLASK /(N FLASK − 1) to account for the correlation of individual measurements and mean) when using either no PME at all or when using PME estimates from shape-noise free sims or from the full sims. The blue dots show the corresponding χ 2 values when using the heuristically motivated analytical treatment of masking and survey geometry presented in the main text. The grey dashed line represents the number of data points and should be the average χ 2 if we had a perfect covariance model (note that for this comparison we have not performed any parameter fitting).
contribution to the covariance matrix (Philcox et al. 2020). The way we include this is by fitting a correction to the number density αn such that the excess power is compatible with zero. In order to do so we minimize the following χ 2 : where ∆C ,th is varied in the range 1.2 2 C ,th − 0.8 2 C ,th (so we are allowing for 20% uncertainty in the bias for the fit). The resulting values for αn can be found in Table 4. In Figure 1, Figure 2 and Table 1 one can see that depleting the lens galaxy densities in our fiducial covariance model by these factors has a negligible effect on both maximum posterior χ 2 and parameters constraints.  2018) allowing for a 20% uncertainty in the galaxy bias (shaded regions) for 2 redshift bins (bin 4 in blue and bin 5 in orange). Horizontal dashed lines are just to guide the eye. If the shot-noise were to be completely Poissonian, the measured and predicted ratios would agree, however we find an excess between 2% and 6%.

Cosmology dependence of the covariance model
In order to evaluate our covariance model, we choose a particular set of cosmological parameters. We do not vary these parameters when sampling our parameter posterior and this may impact the width of our parameter constraints (Hamimeche & Lewis 2008;Eifler et al. 2009;White & Padmanabhan 2015;Kalus et al. 2016). Our main reason for not sampling the covariance model along with the data model is that computing a covariance matrix is computationally too costly for this to be feasible. Recently, Carron (2013) have also indicated that it may indeed be incorrect to vary the covariance cosmology when running MCMC chains. It is only after running the MCMC chains that we can recompute the covariance at our best-fit parameters and rederive our parameter constraints -repeating this process until our constraints have converged (cf. Abbott et al. 2018, for the application of this procedure in DES Y1 data). Therefore, the cosmology at which we compute our covariance is expected to be off from the best-fit cosmology. In this subsection, we investigate how χ 2 , as well as cosmological parameter constraints, shift when computing the covariance at cosmologies that are randomly drawn from the DESY3-like posterior.
We test the robustness of our constraints against the choice of cosmological parameters at which we evaluate the covariance model by taking a set of 100 different cosmologies drawn randomly from the simulated DES Y3 3x2pt posterior and generating 100 lognormal covariance matrices. Using each of these covariances, we estimate posteriors for a given realization of simulated DES Y3 3x2pt data with noise drawn from a fiducial lognormal covariance.
Since it is prohibitively expensive to perform simulated analyses running MCMC chains for each covariance matrix, we use the technique of importance sampling. That allows us to quickly evaluate how these different likelihood modeling choices impact the derived parameter constraints without repeatedly running expensive sampling algorithms. In our importance sampling pipeline, we take a fiducial analysis as a proposal distribution, re-evaluate the likelihoods using the alternative covariance matrix, and compute importance weights as: where C fid is the fiducial covariance in the analysis and C alt the alternative one. If the changes induced by the new covariance matrix in the posterior are not too large, the reweighted samples represent the target distribution (i.e., the posterior for the alternative covariance matrix). So we have: for a function f (Xi) of the posterior samples. Here, pi is the probability of Xi under the target distribution p and qi is the probability of of Xi under the proposal distribution q (see e.g. Owen (2013) and MacKay (2002)).
To diagnose the performance of our importance sampling estimates, we use the Effective Sample Size (ESS): where N samples is the total number of posterior samples used in the estimation. The ESS as defined above quantifies the statistical power of the sample set after the re-weighting process (assuming uncorrelated samples). It is equal to the original sample size re-scaled by the ratio of the variances under each of the distributions (Martino et al. 2017), such that the error of the mean of a quantity x with standard deviation σx under the target distribution can be estimated as σx/ √ ESS. Additionally, since our proposal distribution is itself a weighted sample set, we incorporate both the original and the importance weights in our ESS estimate.
Using the fiducial lognormal covariance matrix we run the nested sampling algorithm MultiNest (Feroz & Hobson 2008;Feroz et al. 2009Feroz et al. , 2019, and perform the importance sampling procedure to estimate parameters using each of the 100 covariance matrices randomly sampled in parameter space. The (S8, Ωm) contours can be seen in Fig. (11). The effective sample sizes for the importance sampled estimates range from 16446 to 18329 (implying a standard error of the mean within 0.78% of the standard deviation for all cases), and the contours show good statistics. As the impact of covariance cosmology is barely noticeable for this range of tested parameters, we repeat the analysis for a few more extreme (and unlikely) cosmologies in appendix F.
These results all confirm that we can safely neglect the impact of the choice of covariance cosmology in DES Y3 3x2pt analysis. One caveat of this conclusion is that we have indeed only varied cosmological parameters (including galaxy bias parameters) but not nuisance parameters (multiplicative shear bias, photometric redshift uncertainties) or parameters that describe intrinsic alignment. However, the DES-Y3 shear and photo-z calibration yield tight Gaussian priors on the corresponding nuisance parameters. And intrinsic alignment is relevant only on small angular scales where the covariance matrix is dominated by sampling noise contributions. Hence, we do not expect the results of this section to change significantly had all parameters been varied.

Random point shot-noise
We also consider the effect of additional shot-noise in the measurements of galaxy clustering resulting from the use of finite numbers of random points. The Landy-Szalay estimator (Landy & Szalay 1993) is estimating the galaxy clustering correlation function inside an angular bin [θ1, θ2] aŝ where DD[θ1, θ2] is the number of galaxy pairs found within the angular bin, RR[θ1, θ2] is the (normalised) number of pairs of random points that uniformly samples the survey footprint and DR[θ1, θ2] the (normalised) number of galaxyrandom-point pairs within the angular bin. If the number density of random points nr is much larger than the number density of the galaxies ng (as is recommended for reduce sampling noise) then both RR and DR must be rescaled by factors of (ng/nr) 2 and (ng/nr) respectively. We stress that the Landy-Szalay estimator was devised at a time of very limited computational resources, where it was prohibitively costly to measure galaxy pair in a large number of random points. Hence, it was vital to minimize random point shot-noise. Nowadays, footprint geometries of photometric surveys are typically characterised by high resolution healpix maps. The most straightforward way to calculate galaxy clustering correlation function is to simply assign a value of galaxy density contrast to each of these pixels and then measure the scalar auto-correlation function of the unmasked pixels. This way, one is avoiding random point shot-noise completely.
Nevertheless, it is still very common to measure w(θ) by means of Equation 79. So we also tested what impact a finite number of random points would have on our analysis. To do so we extended expressions of Cabré & Gaztañaga (2009, see their appendix A) to the case where the same random points are used to estimate w(θ) in each of our redshift bins and also to subtract shear around random points from our galaxy-galaxy lensing correlation functions. Note that this causes a noise contribution to the 2-point function measurements that is correlated among different redshift bins. We assumed a random point density of 1.36/arcmin 2 , which is more that 20 times larger than the density of our most dense lens galaxy sample. From Table 1 it can be seen that not accounting for the random point shot-noise in the covariance leads to an increase in average χ 2 of 1% and to an underestimation of parameter uncertainties by ≈ 0.5%. Hence, this effect can be ignored for our analysis.

Effective densities and effective shape-noise
We are closing this section by spelling out an aspect of covariance modelling that may seem straightforward but which has repeatedly came up in covariace discussions.
If the tracer galaxies used to estimate 2-point correlation functions are weighted according to some weighting scheme, then this may change the effective number densities and the effective shape noise that should be used when evaluating the covariance expressions in Section 4. In the following we will derive how this can be done for each of the 2-point functions in the DESY3 3x2pt data vector.

Galaxy clustering
We start with the galaxy clustering correlation function w(θ). We assume a weighting scheme that is aimed at correcting for non-cosmological density fluctuations resulting from spatially varying observing conditions (as e.g. in Elvin-Poole et al. 2018). This means that the weights assigned to each galaxy in fact sample a weight map that spans the entire footprint.
Instead of measuring w(θ) from the weighted galaxies by means of, say, the Landy-Szalay estimator (Landy & Szalay 1993) it will be more convenient to think of the galaxy density contrast as a pixelized field on the sky. Further more, we will assume that the weight map has been normalised such that w = 1 (which can always be done without changing the outcome of the weighted measurement). Consider pixel i with galaxy count Ng,i and weight wi. If ng is the average galaxy density of the unweighted sample, then by taking expectation values with respect to many Poissonian shot-noise realisations (and hence ignoring fluctuations of the underly-ing matter density field) we get where Apix is the area of each pixel and the second to last line serves as definition of δg,i and needs the fact that we demanded w = 1. Note that these equation are only valid for an ensemble of observations that shares the same weight maps and differs only in their shot-noise realisations.
From the set of all pixels we can now estimate w(θ) within a finite angular bin [θ1, θ2] aŝ where the symbol ∆ ij [θ 1 ,θ 2 ] in the double sum over all pixels is 1 when the distance of the pixel pair i, j is within [θ1, θ2] and 0 otherwise. Note that we assume an enumeration of the pixels and that we demand i > j in the sum in order to not count any pair of pixels twice.
If shot-noise is the only source of noise, then it is straight forward to calculate the variance of this measurement as In the last line, Npair,g[θ1, θ2] is the number of unweighted galaxy pairs within the angular bin [θ1, θ2] . Note that in the presence of clustering, this should be calculated from a set of random points instead of from the actual galaxy catalog. The first factor on the right side of Equation 85 is what the shot-noise variance ofŵ should be in the absence of a weighting scheme. The second term is a 2-point function of the weight map itself. If the weight map has a whitenoise power spectrum, then this factor will be close to 1 in any angular bin that doesn't include angular distances of 0. This means that at large enough scales the last line of Equation 85 looks like the covariance for plain Poissonian shot-noise without any notion of an effective number density. This maybe surprising, but it stems from the fact that the weighting scheme we assumed does not simply multiply the galaxy density contrast field. Instead it reverses an already existing depletion of galaxy density from non-cosmological density fluctuations.
In conclusion, the effective number density that should be used to compute the covariance ofŵ[θ1, θ2] is (86)

Galaxy-galaxy lensing
We move on to consider the galaxy-galaxy lensing correlation function γt [θ1, θ2]. We assume that the lens galaxy sample comes with weights derived from a weight map w l as in the previous subsection while each source galaxy j has a weight w s j which does not come from an entire weight map but is instead the result of the individual quality of shapemeasurement for this galaxy. A measurement of γt can be constructed aŝ Here, δ l,i is the galaxy density contrast of the lenses defined in analogy to the previous subsection, t,j→i is the tangential component of the shear of source j with respect to lens galaxy i and wj is the weight of source galaxy j. Note that due to our definition of the lens galaxy density contrast this estimator already includes subtraction of shear around random points. If shot-noise and shape-noise are the only sources of noise, then it can be readily shown that (only with w s j = 1 = w l j !). Here N pair,ls [θ1, θ2] is the number of unweighted lens-source pairs in the angular bin [θ1, θ2], Ns is the total number of source galaxies and j = 1,j + i 2,j is the complex intrinsic ellipticity of source galaxy j.
Note that the final expression in Equation 88 explicitly allows for the possibility that the source weights w s j are correlated with the intrinsic ellipticities j of the source galaxies. One can interpret Equation 88 as with the effective dispersion of intrinsic ellipticity per shear component given by One subtlety here is that the above derivation requires w s j = 1. The above expressions mus be modified is this is not the case or when taking into account responses Rj of a shape catalog generated with metacalibation (Sheldon & Huff 2017). We detail what to do in the latter case in appendix G.

cosmic shear
For cosmic shear we follow Schneider et al. (2002) and construct a measurement of ξ+ from a set of sources aŝ If shape-noise is the only source of noise and if the intrinsic ellipticities of galaxies are not correlated with their weights, then the variance ofξ+ is given by Here Npair[θ1, θ2] is the number of source galaxy pairs in the bin [θ1, θ2] and we have replaced each expectation value 2 1/2,i w 2 i by σ 2 ,eff from Equation 90. Note that we again assumed wj = 1 and that this may require re-scaling of both weights and σ when using shape measurements from metacalibration.

Testing validity of effective shape noise
To test the validity of our expression for effective shape-noise in Equation 90 we run a sub-sample covariance estimator on our data (see e.g. Friedrich et al. 2016). In particular, we divide all of our source and lens galaxy samples into 200 randomly chosen sub-samples and measure the galaxygalaxy lensing correlation function of each source-lens bin combination. As a result we obtain 200 measurements ofγt in each source-lens bin combination. Since we employ completely random sub-sampling, i.e. without any regard for e.g. a division of our footprint into sub-regions, the sample covariance of these 200 measurements will almost exclusively be dominated by shape-noise and shot-noise. This is even more so, because the lens and source densities of the subsamples are very low.
In Figure 12 we show the ratio of the variances of the 200 galaxy-galaxy lensing measurementsγt in the different lens-source bin combinations to Equation 89. Assuming that the sub-sample covariances follow a Wishart distribution we find that these ratios are perfectly consistent with 1. This indicates that Equation 90 indeed yields an accurate effective shape-noise dispersion, and that one should indeed use the plain density of lens galaxies (as opposed to any notion of effective density) when evaluating covariance expressions.

A SIMPLE χ 2 TEST
In this short section we present a simple χ 2 test that does not rely on the linearized framework. However, it has the disadvantage of not addressing the impact on the estimation of parameters. Here we generate a large number of "contaminated" data vectors (we use 1, 000) by a Gaussian sampling of a given covariance matrix that includes different effects and to compute a χ 2 distribution from these data vectors using a fiducial covariance matrix. The resulting shifts in the for the shape-noise contribution to the covariance (again using Equation 90 to calculate the effective shape-noise dispersion σ ,eff ). Each row displays the variances measured for a different source redshift bin and vertical dashed lines separate points belonging to different lens redshift bins (1-5 from left to right). Assuming that the covariances estimates have a Wishart distribution we calculate the covariance matrix of these ratios (cf. Taylor et al. 2013) and find that they are consistent with 1 (both for the cosmic shear and galaxy-galaxy lensing variances).
mean value of χ 2 and their standard deviations give another benchmark for the importance of the different effects considered here. We show the results of this test in Figure 13. Note that the relative increases in χ 2 follow closely what we obtained within the linearized likelihood framework in Figure  1. This indicates that the dominant way in which covariance errors cause χ 2 offsets is not through the altered scatter of maximum posterior parameter locations but simply through using an erroneous inverse covariance when computing χ 2 . That also justifies our usage of the linearized likelihood framework since any impact of non-linear parameter dependencies on parameter fitting can be expected to be even less relevant then linear fitting in the first place.

DISCUSSIONS AND CONCLUSIONS
In this paper we have presented the fiducial covariance model of the DES-Y3 joint analysis of cosmic shear, galaxygalaxy lensing and galaxy clustering correlation functions (the 3x2pt analysis). We then investigated how the assumptions and approximations of that model (including the assumption of Gaussian statistical uncertainties) impact the distribution of maximum posterior χ 2 and maximum posterior estimates of cosmological parameters. The fiducial covariance matrix of the DES-Y3 3x2pt analysis uses the formalism of Krause & Eifler (2017) to model super-sample covariance as well as the trispectrum contribution to the covariance. The model for the Gaussian covariance part (i.e. the contributions from the disconnected 4-point function) correctly takes into account sky curvature and includes analytical averaging over the finite angular bins in which the 2-point functions are measured. Furthermore, the galaxy clustering power spectra that enter our calculation of the Gaussian covariance part are computed using the non-Limber formalism of Fang et al. (2020b) and also include modelling of redshift space distortions. The finite survey area of DES-Y3 is incorporated in the covariance model via the f sky approximation (except in the pure shape-noise and shot-noise terms where we follow Troxel et al. 2018b).
In order to perform our validation tests for the DES-Y3 covariance matrix we developed a plethora of new modelling ansatzes and testing strategies which are applicable in general. These new techniques are: • We have motivated and devised a way of drawing realisations of the 3x2pt data vector from a non-Gaussian distribution in order to test the accuracy of our Gaussian likelihood assumption. • We have derived analytic expressions for angular bin averaging of all four types of 2-point correlation functions included in the 3x2pt vector (ξ+, ξ−, γt, w). These expressions correctly account for sky curvature. To the best of our knowledge, an analytic treatment of bin averaging for cosmic shear ically derived how covariance model inaccuracies influence the distribution of maximum-posterior χ 2 and of maximumposterior parameter estimates. The results we presented also allow for the possibility of including the Gaussian priors on certain model parameters and can be used to analytically estimate the impact of covariance errors on cosmological likelihood analyses. • By fitting an effective number density to the high-plateau of galaxy clustering C measurements we have estimated how much the assumption of Poissonian shot-noise influences our likelihood analysis. This is similar in spirit to the RASCALC technique presented by Philcox et al. (2020), and we agree with those authors that non-Poissonian shot noise can be viewed as an effective description of how short-scale non-linearities in galaxy clustering influence the covariance. • We calculated covariance matrices for 100 different sets of cosmological and nuisance parameters randomly drawn from a simulated likelihood chain. This allowed us to investigate whether calculating our covariance model at a reasonable, but wrong point in parameter space significantly impacts our analysis. This was done both within our linearised likelihood framework and by using importance sampling to quickly evaluate the 100 non-linear likelihoods. • We have derived how the 2-point correlation function of weight maps influence the covariance of galaxy clustering 2point function measurements. In that context we have also found that traditional ways of deriving an effective number density for a given set of galaxy weights are erroneous when those weights are aimed at undoing a suspected depletion of galaxy density (e.g. due to variations in observing conditions).
• We have derived an expression for the effective dispersion of intrinsic source shapes for the situation when source galaxy weights are correlated with galaxy ellipticity. We have also shown how metacalibation responses (Sheldon & Huff 2017) enter that expression for effective shape-noise. • We have described a clean sub-sample covariance estimation scheme that directly measures the sampling noise contributions to the covariance from a given data set. We then used the resulting covariance estimates to test the validity of our assumed effective shape noise values. • We have employed the hybrid covariance estimation technique PME (Friedrich & Eifler 2018) to efficiently evaluate the importance of individual contributions to the covariance from only a limited set of simulated data (in our case: 200 realisations of the 3x2pt data vector including shape-noise and 100 realisations without shape noise). • We have devised a treatment of survey geometry in covariance modelling that improves upon existing approximations (of e.g. Efstathiou 2004) and we have demonstrated how to carry those approximations from harmonic space to real space.
Using these results we perform several tests for the fiducial DES-Y3 3x2pt covariance matrix and likelihood model, with the following conclusions: • The assumption of Gaussian statistical uncertainties is sufficiently accurate (cf. Section 6.1). Hence, knowledge of the covariance of the 3x2pt data vector is sufficient to model our statistical uncertainties. The main assumption made to arrive at this conclusion is that non-Gaussian error bars are primarily a large-scale problem and that at small scales the number of modes present within the DES-Y3 survey volume converges to a Gaussian distribution via the central limit theorem. • The non-Gaussian part of the covariance has a negligible impact on both maximum posterior χ 2 and parameter constraints (cf. Section 6.2). This statement is not a general one but only holds for the specific DES-Y3 3x2pt analysis setup. The main assumption made to arrive at that conclusion is that the CosmoLike model (Krause et al. 2016) or the log-normal model (Hilbert et al. 2011) for the non-Gaussian covariance do not vastly underestimate the true covariance. Given the results of Sato et al. (2009);Hilbert et al. (2011) we find this a safe assumption. • Of all covariance modelling assumptions investigated in this paper the f sky approximation (made in the mixed term and cosmic variance term of our covariance model) has the largest effect on maximum posterior χ 2 . On average it increases χ 2 between measurement and maximum posterior model by about 3.7% (∆χ 2 ≈ 18.9) for the 3x2pt data vectors and by about 5.7% (∆χ 2 ≈ 16.0) for the 2x2pt data vector (cf. Table 1). • However, neither f sky approximation nor any other covariance modelling detail tested in this paper (cf. Table 1 ; with the exception of finite bin width, see next point) has any significant impact on the location and width of constraints on the parameters Ωm, σ8, w. • The only exception to this statement is finite angular bin width which -if not taken into account in the mixed term and cosmic variance term of the covariance model -significantly increases the scatter of maximum posterior parameters (without increasing the inferred constraints accordingly). However, finite bin width has been taken into account in the past in an approximate manner -see e.g. Krause et al. (2017). • The fact that we do not know the true cosmological parameters of the Universe forces us to evaluate our covariance model at a wrong set of parameters. Even when iteratively adjusting those parameters to the maximum posterior parameters of the analysis, the parameters of the covariance model will scatter around the 'true' cosmological / nuisance parameters. We consider this an irreducible covariance error and find that it increases the maximum posterior scatter of Ωm and σ8 by about 3% and that of the dark energy equation of state parameter w by about 5% (cf. Section 6.8 and Table 1). At the same time, we find this effect to have a negligible impact on maximum posterior χ 2 .
In summary, we have shown that our fiducial covariance and likelihood model underestimates the scatter of maximum posterior parameters by about 3-5%, which is mostly caused by uncertainty in the set of cosmological and nuisance parameters at which we evaluate that model. On average, the χ 2 between maximum posterior model and measurement of the 3x2pt data vector will be ∼ 4% higher than expected with perfect knowledge of the covariance matrix. This is mainly caused by our use of the f sky approximation. We have devised an improved treatment of the full survey geometry (cf. Section 6.6) but for the reason mentioned above this was only used to test the impact of the f sky approximation on parameter constraints.
Given the small impact that we estimated from the unaccounted effects in the covariance modelling, we conclude that the fiducial covariance model is adequate to be used in the 3x2pt DES-Y3 analysis. While our validation of this covariance model has been carried out with a preliminary set of scale cuts and redshift distributions, we don't expect qualitative changes for the final DES-Y3 analysis setup.
While the DESY3 specific outcomes of our study can not straightforwardly be transferred to other surveys and analyses, our methodological innovations will be useful tools in the covariance and likelihood validation of future experiments. This work made use of the software packages GetDist (Lewis 2019), ChainConsumer (Hinton 2016), matplotlib (Hunter 2007), and numpy (Harris et al. 2020).

ACKNOWLEDGEMENTS
We would like to thank the anonymous journal referee for their helpful comments.

DATA AVAILABILITY
The DES-Y3 3x2pt covariance matrix and likelihoods will be made public upon publication of our final data analysis. C++ and python tools to configure FLASK as described in Section 4.3 are available at https://github.com/OliverFHD/ CosMomentum . Tools to compute Gaussian and halomodel covariance is available at https://github.com/CosmoLike/ CosmoCov .

APPENDIX A: CURVED SKY FORMALISM
The angular clustering correlation function of galaxies w(θ) is given in terms of the galaxy clustering power spectrum C gg as where n is the galaxy density per steradian. The term proportional to 1 n is usually ommited (cf. Ross et al. (2011)) since it sums up to 16 1 2πn 2 + 1 2 P (cos θ) = 1 2πn δD(cos θ − 1) which has to be interpreted as a 2-dimensional Dirac delta function on the sphere. According to de Putter & Takada (2010) (see also Stebbins 1996) the galaxy-galaxy lensing correlation function γt(θ) is given in terms of the galaxy-convergence cross-power spectrum C gκ as γt(θ) = 2 + 1 4π where P m are the associated Legendre Polynomials. Finally, he cosmic shear correlation functions ξ±(θ) are given by where x = cos θ, C E is the E-mode power spectrum of shear and we assume B-modes to vanish. The functions G ± ,2 (x) 16 from 2 +1 2 P (x)P (y) = δ D (x − y) -see N. Bronstein & A. Semendjajew (1979) for this and other properties of Legendre polynomials.
are defined in eq. 4.18 17 of Stebbins (1996). Eq. A4 can be expressed in terms of associated Legendre polynomials by using eq. 4.19 of Stebbins (1996), which gives In appendix B we show how to obtain A4 from the notation in Stebbins (1996). It can be seen from above that each of the considered 2-point correlation functions can be written in terms of the corresponding power spectra as
This way one directly arives at our expression for the shear correlation functions, equation A4. To account for finite bin width in equation A4 and in the covariance ofξ± one has to perform the area-weighted bin average of the functions G + ,2 (cos θ) ± G − ,2 (cos θ) . To do so one can insert the relations x 1 x 2 into equation A5. In summary this means one has to exchange the functions G + ,2 (cosθ) ± G − ,2 (cosθ) as follows: These expressions can be very efficiently calculated and pretabulated e.g. with the help of the gnu scientific library (Galassi et al. 2009).

APPENDIX C: MASKING IN REAL SPACE COVARIANCES
DES observations don't cover the entire sky, but are located within a survey mask, described by a function W (n) which is = 1 if we have observed the sky at locationn and zero otherwise 18 . In this appendix we derive how this masking changes the covariance of any measured 2-point statistics.
We start by computing, how many galaxy pairs we expect to find within our mask (assuming that galaxies do not cluster). If n1, n2 are the number densities of 2 different tracer samples, then the expected number of pairs dNpair(θ) with angular separation within [θ, θ + dθ] is given by (cf. 18 A map of the mask will come with a finally resolution, in which case the fractional values 0 < W < 1 of the mask will describe the completeness of observations within the map resolution. Troxel et al. 2018b) Using the fact that this becomes where in the last steps we have defined the angular power spectrum of the mask through (2 + 1)C W = m |W m | 2 as well as the angular 2-point function of the mask, ξ W (θ). The total number of galaxy pairs in a finite angular bin [θmin, θmax] is then The 2-point correlation function of 2 scalar random fields, From the last line it can be seen thatξ ab -apart from the normalisation by Npair[θmin, θmax]/nan b -is exactly the angular space counter part of the Pseudo-Cells estimator in the corresponding harmonic space (cf. Efstathiou 2004). Why this is the case can be understood most easily in the limit of an infinitesimal angular bin. In this limit where the second line shows that the convolution of mask and signal power spectrum in harmonic space (Efstathiou 2004) becomes a simple multiplication in angular space. Especially, normalisation by Npair(θ) is the angular analog of multiplication with the inverse mode-coupling matrix that appears in harmonic space. LetĈ ab be the pseudo-C estimator we identified in Equation C6, i.e. we writê Then the covariance of the 2-point function measurements between the fields δa&δ b and δc&δ d within angular bins [θ ab − , θ ab + ] and [θ cd − , θ cd + ] is given by Assuming that δa, δ b , δc, δ d are Gaussian random fields and defining the symbols it is straight forward to show that Cov Ĉ ab 1 ,Ĉ cd 2 is given by (cf. Efstathiou 2004) Cov Ĉ ab 1 ,Ĉ cd 2 = 1 (2 1 + 1)(2 2 + 1) At this point e.g. Efstathiou (2004); Varshalovich et al. (1988) follow with the approximation We however find that for fixed values of 1 and 2 the expression and that the area around each of these maxima contributes a similar amount to the sums over 3 and 4. Hence, we opt instead for the approximation Cov Ĉ ab 1 ,Ĉ cd 2 ≈ 1 4 C ac 1 C bd 2 + C ac 2 C bd 1 + C ac 1 C bd 1 + C ac 2 C bd 2 (2 1 + 1)(2 2 + 1) + In practice, both approximations yield very similar results and they are valid on scales 1, 2 which are much smaller than the typical scales of the mask W (Efstathiou 2004;Varshalovich et al. 1988). Unfortunately, the DES-Y3 analysis mask has features and holes over a large range of scales. Hence, the angular scales of interest in the 3x2pt analysis are never strictly smaller than the scales of our mask. Hence, Equation C15 is not sufficiently accurate in our case and infact significantly overestimates our covariance matrix. In Figure 8 we explain a simple scheme that can be used to correct for this: Calculating the covariance ofξ ab [θ ab − and θ ab + ],ξ cd [θ cd − , θ cd + ] within the Gaussian covariance model (see also Section 3) requires integration over all pairs of locations within our survey mask that fall into the angular bins  Figure 8 visualizes this for the mixed terms in the covariance, where one of the correlation functions ξ ac or ξ bd is due to sampling noise such as shape-noise of shot-noise and hence is proportional to a Dirac delta function. In that case, the integration is over pairs that share one end point. Now the approximation made e.g. in Efstathiou (2004) or by our Equation C15 assumes that also the correlation function between the other two end points effectively acts as a delta function -at least with respect to the smallest scale features in the survey mask. We find that this is not the case for the DES-Y3 mask and that it contains features on all scales relevant to our analysis. But Figure 8 also indicates a simple way to fix this: approximating the integrand over the distance of the two remaining endpoints by a delta function roughly overestimates the integral by a factor equal to one over the fraction of small scale hole in the survey footprint compared to the average scale at which the 2-point function between the 2 end points decays. And multiplying the the mixed terms in the covariance by this fraction can serve as a next-to-leading order correction to our Equation C15. By similar arguments one can deduce that the cosmic variance terms (terms where none of the end points must be joined) can be corrected by performing this multiplication twice.

APPENDIX D: MOTIVATION FOR OUR RE-SCALING ANSATZ FOR MASKING EFFECTS
In this appendix we make some of the arguments presented in Section 6.6 more precise.
We have argued in Section 6.6 that the approximation of Efstathiou (2004) for evaluating the impact of masking on the 2-point function covariance can roughly be understood as making the replacements dΩ a . . . W (Ω a )ξ ac (θ ac ) ≈ξ ac dΩ a . . . W (Ω a )δ 2 Dirac (Ω a − Ω c ) and whereξ ac is the integral of ξ ac (θ ac ) over the 2-dimensional angular distance vector θ ac andξ bd is the integral of ξ bd (θ bd ) over θ bd . Within these approximations the right side of Equation with proportionality coefficients that do not depend on the survey mask. So if our above understanding of the approximation proposed by Efstathiou (2004) is (at least approximately) correct, then its ratio with respect to the f sky approximation (i.e. the approximation where one computes the covariance of a full-sky survey and then re-scales it with the sky fraction f sky of the survey footprint) should be given by the inverse ratio of the exact value of N ab pair [θ−, θ+] to the f sky calculation N ab pair,f sky [θ−, θ+] = 4π 2 (θ 2 + − θ 2 − )f sky nan b .
In Figure D1 we show the ratio of Npair to N pair,f sky (green dashed lines) and compare it to the ratios of different covariance terms when computed with either the f sky approximation or the approximation of Efstathiou (2004) (cosmic variance term: blue dotted lines; mixed term: solid orange lines). The upper panel of the figure computes these ratios for the galaxy clustering 2-point function w(θ) in the first lens bin of our fiducial configuration while the lower panel considers the cosmic shear correlation function ξ+(θ) for our first source bins (other redshift bins and 2-point functions behave similarly). For w(θ) these different ratios indeed closely agree with each other. The agreement between the ratio of pair counts and the ratio of the different approximations for the mixed terms is especially striking. It is most likely caused by the fact that for the mixed covariance terms one of the two replacements in Equations D3 and D4 is actually exact. Even for ξ+(θ) the ratio of the different approximations for the mixed term is well described by the ratio of pair counts on most scales. For the cosmic variance of ξ+(θ) one can on the other hand observe a strong deviation. This does not necessarily signify a breakdown of our arguments but may be caused by the fact that cosmic shear is a spin-2 field. The calculations of Efstathiou (2004)  in fact only hold for scalar fields and our extension of their formulae to shear correlation functions is only approximate (see e.g. Challinor & Chon 2005, for more general calculations). Since we have identified the mixed terms to carry the strongest impact of masking on the total covariance, we do not address this any further. Instead, we consider the agreement between pair count ratios and mixed term ratios observed in Figure D1 as sufficient justification for the rescaling ansatz of the different covariance terms presented in Section 6.6.

APPENDIX E: PRECISION MATRIX EXPANSION TO INVESTIGATE THE IMPACT OF MASKING ON INDIVIDUAL COVARINCE TERMS
In this appendix we briefly summarize the PME method that went into Figure 9. The covariance of the 2x2pt (i.e. non-cosmic-shear) part of our data vector has contributions from shape-noise because of the presence of the mixed term described in Section 4. To pinpoint further which parts of our analytic covariance contribute to the elevation in χ 2 (and to further motivate our heuristic modelling ansatz for masking effect in the covariance presented in Section 6.6), we re-run 100 of the FLASK simulations with shape-noise turned off. We then use the covariances estimated from the different FLASK runs to derive corrections to our covariance model. This can be done -even with only a limited number of simulationswith the method of precision matrix expansion (PME) that was described by Friedrich & Eifler (2018). At the 1st order their expansion estimates the precision matrix Ψ (i.e. the inverse covariance matrix) aŝ Here, the matrix B model can be either the full covariance model, in which caseB is the full covariance estimated from FLASK or it could be the shape-noise free part of the covariance, in which caseB will be the covariance estimated from the shape-noise free FLASK simulations. Friedrich & Eifler (2018) have also derived a 2nd order correction to Equation E1, but given the small magnitude of our observed χ 2 elevation we restrict ourselves to the 1st order expansions which should also reduce the noise of the PME. Note that Equation E1 does not contain the inverse of any noisy matrix. This is why PME works well even in the presence of only few numerical simulations (a benefit that is even further boosted because the matrix B can be chosen to represent only sub-parts of the covariance). For each FLASK measurement of the 2x2pt data vector we estimate the 1st order PME from the remaining 196 FLASK data vectors (respectively from the ∼ 100 shape-noise free data vectors). The average resulting χ 2 values between each data vector and the mean of all data vectors are displayed in Figure 9 and compared to the χ 2 values obtained when applying the analytic masking corrections presented in Section 6.6 to either the shape-noise free covariance terms or the full covariance. The average χ 2 when using the analytic, best-guess covariance matrix is ≈ 318.8 for a total of 302 data points in the 2x2pt data vector. This corresponds to a bias in χ 2 of about 5.5%. The PME estimate of the inverse covariance manages to push this down to ≈ 307.9 (≈ 304.8 with our analytic ansatz) hence decreasing the bias in χ 2 to about 1.9% (< 1% for the analytic anasatz). If the PME correction term is computed with the shape-noise free FLASK covariance, then the bias is only slightly reduces to χ 2 ≈ 314.3 (≈ 316.6 with our analytic ansatz). Hence, the shape-noise dependent mixed terms in the covariance indeed seem to be the main cause of our remaining χ 2 offset. This was also found by Joachimi et al. (2020) for the latest analysis of the Kilo Degree Survey. These mixed terms do not depent on the connected 4-point function of the density field (cf. Section 3) and the only approximation we make in their calculation is the treatment of our survey footprint Figure F1. (S 8 , Ωm) constraints for a given noisy realization of the DES Y3 3x2pt data vector analyzed using: a fiducial covariance matrix (blue); a covariance matrix evaluated at σ 8 = σ fiducial 8 − 2σσ 8 (green); a covariance matrix evaluated at σ 8 = σ fiducial 8 + 2σσ 8 (red). The green and red posteriors were obtained by importance sampling the fiducial samples.
through the f sky approximation. Hence, we follow Joachimi et al. (2020) in our conclusion that this approximation is the main driver of the residual errors in our covariance model.

APPENDIX F: IMPACT OF EXTREME COSMOLOGIES ON PARAMETER CONSTRAINTS
To demonstrate that the importance sampling technique employed in Section 6.8 indeed manages to capture even strong changes in the likelihood, we repeat the tests presented there with covariance matrices that drastically differ from our fiducial covariance model. In particular we shift the value of σ8 for which the covariance model is evaluated by ±2σ of the marginalised σ8 constraints expected from DES-Y3. Note that is a radical change because it ignores parameter degeneracies, i.e. such a shift of σ8 without changes in other parameters would be detected at high significance. Figure  F1 shows the likelihood contours in the S8-Ωm plane obtained from both our fiducial covariance and from importance sampling with the altered covariance matrices. One can now clearly see a change in contour width. But as we have show in Section 6.8, this effect is far less significant for realistic parameter uncertainties in the covariance model.