Disconnected pseudo-$C_\ell$ covariances for projected large-scale structure data
Carlos García-García David Alonso Emilio Bellini

  Signatures of the orbital angular momentum of an infrared light beam in the two-photon transition matrix element: A step toward attosecond chronoscopy of photoionization
    Gopal Dixit
    Received: date / Revised version: date
=========================================================================================================================================================================



  

§ INTRODUCTION

    The two-point correlation of different fields projected on the celestial
    sphere is one of the most common observables used in the analysis of large
    datasets in astrophysics, from studies of the Cosmic Microwave Background
    <cit.>
    to large-scale structure and weak lensing surveys
    <cit.>.
    Using these two-point functions, one achieves a high level of data
    compression (with respect to the size of the raw datasets – time-ordered
    data, images or catalogs). They can also be directly used to constrain
    cosmological and astrophysical parameters assuming that one can model
    their likelihood. This is usually done by assuming that the two-point
    functions are Gaussianly distributed, which is often a good approximation
    due to the central limit theorem
    <cit.>. In this case, the only
    obstacle that remains is being able to estimate the covariance matrix of a
    set of two-point correlators. Since the form of this covariance directly
    affects the posterior parameter uncertainties, a precise determination of
    it is of paramount importance. In large-scale structure experiments, this
    has often been resolved by making use of one's own data through resampling
    techniques
    <cit.>,
    or by generating a large number of mock realizations
    <cit.>. With
    the advent of the larger current and future surveys, the increasing size
    of the data vector and of the volume to be simulated has made this
    solution impractical, and fully analytical and hybrid estimators are now
    being used.
  
    The problem of producing accurate analytical estimates of the covariance matrix for large-scale structure data has seen significant progress in the last few years <cit.>. As described in <cit.>, the covariance matrix recieves three main contributions:
    
      
  * Gaussian covariance: this is the contribution to the covariance from the disconnected part of the trispectrum of the different fields involved (also called the “disconnected” covariance). In simpler terms, this is the covariance matrix one would obtain if all fields involved were Gaussianly distributed.
      
  * Connected non-Gaussian covariance: this is the contribution from the connected trispectrum (which would vanish if all fields were Gaussianly distributed).
      
  * Super-sample covariance: this is the additional coupling between different scales induced by density fluctuations on scales larger than the volume mapped. This term also vanishes for Gaussian fields.
    
  
    On most scales relevant for cosmological studies, the Gaussian contribution dominates the error budget, although the connected and super-sample terms cannot be neglected <cit.>. The Gaussian contribution is trivial to compute for fields observed over the full sky:
    
    Cov(C^ab_ℓ,C^cd_ℓ)=δ^K_ℓℓ'C^ac_ℓ C^bd_ℓ+C^ad_ℓ C^bc_ℓ/2ℓ+1,

    where C^xy_ℓ is the angular power spectrum between two maps x and y on multipole ℓ. Unfortunately, the presence of a sky mask in general induces non-trivial couplings between different ℓs, which must be accurately estimated in order to produce unbiased evaluations of the parameter likelihood <cit.>.
  
    In this paper, we will present and generalize methods developed in the context of CMB experiments to account for the impact of survey geometry on the Gaussian part of the power spectrum covariance matrix <cit.>[See also <cit.> for a similar application to the problem of 3D power spectrum covariances.], and will study in detail the performance of these methods for large-scale structure and weak lensing datasets. We have also implemented these approximations in the public code NaMaster[<https://github.com/LSSTDESC/NaMaster>] <cit.>, making the computation of accurate Gaussian covariance matrices significantly simpler for the community.

    The paper is structured as follows: Section <ref> presents the methods and approximations used to calculate accurate covariances. In Section <ref> we test the methods against Gaussian simulations and study their performance as well as their impact on the final cosmological parameter estimation. We then summarize our results and conclude in Section <ref>. Appendix <ref> presents the performance of these methods in the flat-sky approximation, and we provide technical details of the software implementation in Appendix <ref>.

  

§ ANALYTICAL GAUSSIAN COVARIANCES

    

 §.§ Preliminaries

      We will deal with spin-0 and spin-2 fields defined on a 2-dimensional space. In two dimensions, spin-s fields in general have two components a( x)=(a_1( x),a_2( x))[E.g. for CMB polarization, a spin-2 field, these components are the Stokes parameters (Q,U), while for cosmic shear the two components are usually labeled (γ_1,γ_2).]. Forming a complex number from these components, a_1+ia_2, spin-s fields transform, under a coordinate rotation with angle ψ, as a_1+ia_2→ (a_1+ia_2)exp(i sψ). Thus, spin-0 fields are invariant under rotations, and are usually expressed as real-valued fields with a single component.
      
      Given a field a( x), with 1 (spin-0) or 2 components (spin-2), defined on the coordinates x, we define its generalized Fourier coefficients as
      
    a_ k=x †kx  a( x),

      where the operator x denotes an integral or sum over all values of the coordinates x, and †kx are a set of orthogonal functions. We will also assume that the †kx are a complete set of basis functions, in which case:
      
    x †kxlx=1Δ^x( k, l),
       k kx†ky=1Δ^k( x, y),

      where k denotes an integral over all possible generalized Fourier coefficients k, and Δ^x and Δ^k are generalized delta functions, defined through their action on functions of x or k:
      
    k f( k)Δ^x( k, l)≡ f( l),
       x f( x)Δ^k( x, y)≡ f( y).

      For a spin-s quantity, lx can be written in terms of two spin-raising and spin-lowering operators, ƪ and ƪ̅, and a set of scalar orthogonal functions q( l, x) as:
      
    lx=-β_ℓ,s/2(
            [     ƪ^s+ƪ̅^s  i(ƪ^s-ƪ̅^s); -i(ƪ^s-ƪ̅^s)     ƪ^s+ƪ̅^s ]) q( l, x),

      where β_ℓ,s is a normalization factor defined in Table <ref>.
    
      Finally, we will assume that all fields are Gaussian stochastic fields that are additionally statistically isotropic. As a consequence of the latter, different generalized Fourier modes are uncorrelated:
      
    ⟨ a_ k b^†_ l⟩≡ C^ab_ℓ K Δ^x( k, l),

      where K is a volume factor (see below) and C^ab_k is the power spectrum. Defined this way, the power spectrum is a matrix, with elements
      
    ⟨ a^α_ k(b^β_ l)^*⟩≡( C^ab_ℓ)_αβ K Δ^x( k, l),

      where a^α is the α-th element of field a. It will often be useful in what follows to think of C^ab_ℓ as a 1-dimensional vector that we will denote by vec( C^ab_ℓ). To do so, we simply map the two indices (α,β) into a single number A, such that vec( C^ab_ℓ)_A = ( C^ab_ℓ)_αβ.
    
      

      All the functions and operators above can be specialized to fields defined on the sphere or the 2D plane (flat sky approximation) as described in Table <ref>.

    

 §.§ The pseudo-C_ℓ method

      This section provides a very brief introduction to the pseudo-C_ℓ power spectrum estimator. Further details can be found in e.g. <cit.>. 
      In any practical situation we do not have access to maps of a given field a over the full sky, but rather to a weighted or masked version of them
      
    ã( x)≡ w_a( x) a( x),

      where w_a is commonly called the “mask”. Due to the convolution theorem, the generalized Fourier coefficients of the masked field will be a convolution of the mask and true field coefficients:
      
    ã_ l=   k[x w^a( x) †lx kx] a_ k
    ≡   k ^a M_ l k  a_ k,

      where we have defined the mode-coupling coefficients ^a M_ l k in the second line.
    
      Correlating the generalized Fourier coefficients of two masked fields therefore yields a mode-coupled version of their true underlying power spectrum:
      
    ⟨ã_ lb̃^†_ l⟩=   kq ^a M_ l k⟨ a_ k b^†_ q⟩ ^b M^†_ l q
    
                                                                                  =   K k ^a M_ l k C^ab_k ^b M^†_ l k

    
      The pseudo-C_ℓ estimator then proceeds in two steps:
      
        
  * We first bin different l modes into sets of them called bandpowers (typically bands of similar ℓ or annuli of flat-sky Fourier modes spanning a range of radii). Let us denote a given bandpower by its index q. We must note that it is more appropriate to use bandpower-averaged spectra when the underlying spectrum does not vary much within each ℓ bin. When this is not the case, it is often useful to apply ℓ-weights (such as D_ℓ≡ℓ(ℓ+1)C_ℓ/(2π)). The large-scale structure spectra discussed here are sufficiently flat that the binning used is appropriate. The binned pseudo-power spectrum is:
        
    C̃^ab_q=∑_ l∈ q B_q^ l ã_ lb̃^†_ l,

        where the bandpower weights are normalized such that ∑_ l∈
          qB_q^ l=(KΔ^x( 0))^-1.
        
  * Then, the correlation between bandpowers induced by the mode-coupling coefficients is partially reversed by multiplying C̃^ab_q by the so-called binned “mode-coupling matrix” ℳ, giving the final estimator
        
    vec(Ĉ^ab_q) = ∑_q'(ℳ^-1)_qq'  vec(C̃^ab_q').


        The main advantage of the pseudo-C_ℓ estimator is that the mode-coupling matrix ℳ is directly related to the coupling coefficients ^a M_ l l', and can be computed analytically making use of methods that scale like ℓ_ max^3 (see e.g. <cit.>).
      
    
      For completeness, the mode-coupling matrices for flat-sky and curved-sky fields are given by <cit.>:
      
        
  * Curved sky. After averaging over the harmonic number m, the mode-coupling matrices are:
        
    ⟨1/2ℓ+1∑_m=-ℓ^ℓ vec[ã_ lb̃^†_ l]⟩=∑_ℓ' M^s_as_b_ℓℓ' vec[ C^ab_ℓ'],

        with
        
    M^00_ℓℓ'=(2ℓ'+1) Ξ^00_ℓℓ'(w_a,w_b), M^02_ℓℓ'=(2ℓ'+1) Ξ^0+_ℓℓ'(w_a,w_b) 1,
        M^22_ℓℓ'=(2ℓ'+1)
              (
              [ Ξ^++_ℓℓ'        0        0  Ξ^–_ℓℓ';        0 Ξ^++_ℓℓ' -Ξ^–_ℓℓ'        0;        0 -Ξ^–_ℓℓ' Ξ^++_ℓℓ'        0;  Ξ^–_ℓℓ'        0        0 Ξ^++_ℓℓ' ]),

        where
        
    Ξ^00_ℓℓ'(w,v)  ≡∑_ℓ”P^wv_ℓ”/4πℓℓ'ℓ”000^2
       Ξ^0+_ℓℓ'(w,v)  ≡∑_ℓ”P^wv_ℓ”/4πℓℓ'ℓ”000ℓℓ'ℓ”2-20
       Ξ^±±_ℓℓ'(w,v) ≡∑_ℓ”P^wv_ℓ”/4πℓℓ'ℓ”2-20^21±(-1)^ℓ+ℓ'+ℓ”/2.

        Here the 2-by-3 matrix-like quantities are the Wigner 3-j symbols, and
        
    P^vw_ℓ≡∑_m=-ℓ^ℓ v_ℓ mw^*_ℓ m.

        
  * Flat sky. In this case the averaging over the Fourier-space azimuth happens while binning into bandpowers, and therefore the unbinned mode-coupling matrix is defined before binning. Assuming flat bandpowers, such that B_q^ l=(2π)^2/(L_xL_yN_q), where N_q is the number of Fourier-space modes in the q-th bandpower:
        
    ⟨ vec[C̃^ab_q]⟩=∑_ l∈ q1/N_q∑_ k M^s_as_b_ l k  vec[ C^ab_k]

        with
        
    M^00_ l  k   ≡Ξ̅^00_ l k,
     M^02_ l  k   ≡(
              [  Ξ̅^0+_ l k -Ξ̅^0-_ l k;  Ξ̅^0-_ l k  Ξ̅^0+_ l k ])
     M^22_ l  k   ≡(
              [  Ξ̅^++_ l k -Ξ̅^+-_ l k -Ξ̅^+-_ l k   Ξ̅^–_ l k;  Ξ̅^+-_ l k  Ξ̅^++_ l k  -Ξ̅^–_ l k -Ξ̅^+-_ l k;  Ξ̅^+-_ l k  -Ξ̅^–_ l k  Ξ̅^++_ l k -Ξ̅^+-_ l k;   Ξ̅^–_ l k  Ξ̅^+-_ l k  Ξ̅^+-_ l k  Ξ̅^++_ l k ]),

        where
        
    Ξ̅^00_ l k=(2π/L_xL_y)^2 (w_a)_ l- k(w_b)^*_ l- k,
       Ξ̅^0+_ l k=(2π/L_xL_y)^2 (w_a)_ l- k(w_b)^*_ l- kcos2Δφ,
       Ξ̅^0-_ l k=(2π/L_xL_y)^2 (w_a)_ l- k(w_b)^*_ l- ksin2Δφ,
       Ξ̅^++_ l k=(2π/L_xL_y)^2 (w_a)_ l- k(w_b)^*_ l- kcos^22Δφ,
       Ξ̅^+-_ l k=(2π/L_xL_y)^2 (w_a)_ l- k(w_b)^*_ l- kcos2Δφ sin2Δφ,
       Ξ̅^–_ l k=(2π/L_xL_y)^2 (w_a)_ l- k(w_b)^*_ l- ksin^22Δφ,

        and Δφ is the relative angle between l and k.
      

      Before we move on to covariances, it is worth considering the case of
      unmasked field (i.e w^a( x)=1 everywhere). In this case ^a
        M_ l l'=1Δ^x( l, l'), and therefore different modes are uncorrelated (as should have been obvious). In a non-ideal case where the mask is still sufficiently well behaved (i.e. masks without too much small-scale structure), we can still expect the coupling coefficients ^a M_ l l' to be sharply peaked around l= l'.
  
    

 §.§ Covariance matrices

      So far we have not assumed anything about the statistics of the fields, other than the fact that they are isotropic (Eq. <ref>). This section presents a method to estimate the disconnected part of the power spectrum covariance for the pseudo-C_ℓ estimator.
    
      Let A and F be the vector indices corresponding to the pairs of
      field indices (α,β) and (ϕ,γ), respectively, and let us start by considering the covariance
      
    Σ^AF_ l l'≡⟨ã^α_ lb̃^β *_ lf̃^ϕ_ l'g̃^γ *_ l'⟩-⟨ã^α_ lb̃^β *_ l⟩⟨f̃^ϕ_ l'g̃^γ *_ l'⟩

      The covariance of the binned bandpowers C̃ can then be computed as
      
    Cov( vec(C̃^ab_q)_A, vec(C̃^fg_q')_F)=∑_ l∈ q B^ l_q∑_ l'∈ q' B^ l'_q' Σ^AF_ l l'.

      which can then be used to estimate the covariance of the mode-decoupled bandpowers multiplying it by the inverse mode-coupling matrix twice. I.e., schematically:
      
    Cov(Ĉ)=ℳ^-1· Cov(C̃)·(ℳ^-1)^T,

      where we have suppressed all indices for simplicity. The problem of estimating the pseudo-C_ℓ covariance therefore reduces to estimating Σ^AC_ l l'.
    
      We now make use of Wick's theorem, which states that, for Gaussian fields, ⟨ a b f g⟩=⟨ a b⟩⟨ f g⟩+⟨ a f⟩⟨ b g⟩+⟨ a g⟩⟨ b f⟩. In this case, the expression for Σ^AF_ l l' reads:
      
    Σ^AF_ l l'=   [Kk ^aM^αα'_ l k ^gM^γγ'*_ l' kC^ag,(α'γ')_k][Kq ^bM^ββ'*_ l q ^fM^ϕϕ'_ l' qC^bf,(β'ϕ')_q]^*
       +((g,γ)↔(f,ϕ)),

      where we implicitly sum over repeated indices (e.g. α'), and the second term is equivalent to the first one after swapping the roles of fields f and g. Without any further approximations, for each pair ( l, l'), we would need to perform two 2-dimensional integrals, and therefore the calculation would scale like ℓ_ max^6, quickly becoming unfeasible.
    
      Under the assumption that the coupling coefficients M_ l k are sharply peaked around l= k, we can simplify the expression above approximating the power spectra as constants within the support of the coupling coefficients <cit.>. Explicitly, we approximate 
      
    C^ag,(α'γ')_kC^bf,(β',ϕ')_q≃ C^ag,(α'γ')_(ℓC^bf,(β',ϕ')_ℓ')≡1/2(C^ag,(α'γ')_ℓ C^bf,(β',ϕ')_ℓ'+C^ag,(α'γ')_ℓ' C^bf,(β',ϕ')_ℓ).

      In this case, the expression for Σ^AF_ l l' simplifies to
      
    Σ^AF_ l l'=K^2C^ag,(α'γ')_(ℓC^bf,(β',ϕ')_ℓ') ^agW^αγ,α'γ'_ l l'( ^bfW^βϕ,β'ϕ'_ l l')^*+((g,γ)↔(f,ϕ)),

      where we have defined the covariance coupling coefficients
      
    ^abW^αβ,α'β'_ l l'=k ^aM^αα'_ l k( ^bM^ββ'_ l' k)^*.


      In order to compute these coefficients, let us start by defining the quantities
      
    ^abI^± s_a,± s_b_ l l'≡kxyw_a( x)w_b( y) q^± s_a_ l k( x) [q^± s_b_ l' k( y)]^*,

      where s_a is the spin of field a, and
      
    q^± s_ l k( x)≡β_ℓ,sβ_k,s/2([ƪ^s q( l, x)]^* ƪ^sq( k, x)±[ƪ̅^sq( l, x)]^*ƪ̅^sq( k, x)).


      Now, in what follows, we will be concerned with the auto- and
      cross-correlations of spin-0 and spin-2 fields. Thus, to simplify the
      notation, we will enumerate the different types of coupling coefficients
      that exist for a spin-0 field with a single component that we will call
      δ, in analogy to the projected galaxy overensity, and for a
      spin-2 field, γ, with E and B components, γ_E and γ_B, in analogy to the cosmic shear field. With this setup, all the possible non-zero W^αβ,α'β'_ l l' can be expressed in terms of the I^± s_a,± s_b_ l l' as follows:
      
    W^δδ,δδ_ l, l'=I^0,0_ l l'; 
       W^δγ_E,δγ_E_ l, l'=W^δγ_B,δγ_B_ l, l'=I^0,+2_ l, l';
            W^δγ_E,δγ_B_ l, l'=-W^δγ_B,δγ_E_ l, l'=-i I^0,-2_ l, l';
       W^γ_Eγ_E,γ_Eγ_E_ l, l'=W^γ_Eγ_B,γ_Eγ_B_ l, l'=W^γ_Bγ_E,γ_Bγ_E_ l, l'=W^γ_Bγ_B,γ_Bγ_B_ l, l'=I^+2,+2_ l l'
       W^γ_Eγ_E,γ_Bγ_B_ l, l'=W^γ_Bγ_B,γ_Eγ_E_ l, l'=-W^γ_Bγ_E,γ_Eγ_B_ l, l'=-W^γ_Eγ_B,γ_Bγ_E_ l, l'=I^-2,-2_ l l'
       W^γ_Eγ_E,γ_Eγ_B_ l, l'=-W^γ_Eγ_B,γ_Eγ_E_ l, l'=W^γ_Bγ_E,γ_Bγ_B_ l, l'=-W^γ_Bγ_B,γ_Bγ_E_ l, l'=-i I^+2,-2_ l l'
       W^γ_Eγ_E,γ_Bγ_E_ l, l'=W^γ_Eγ_B,γ_Bγ_B_ l, l'=-W^γ_Bγ_E,γ_Eγ_E_ l, l'=-W^γ_Bγ_B,γ_Eγ_B_ l, l'=i I^-2,+2_ l l'.

      Thus, in principle, we only need to compute 7 different types of terms
      (I^0,0, I^0,±2, I^±2,±2 and I^±2,∓2)). In
      order to simplify these expressions further, we follow <cit.> and neglect all gradients of the masks[It is worth noting that it is possible in principle to avoid this approximation, as demonstrated in <cit.>.], which allows us to relate the different I^± s_a,± s_b through the following set of identities:
      
    xw( x)(ƪ̅^2q( l, x))^* ƪ̅^2q( k, x)
               =xw( x)ƪ^2 q^*( l, x) ƪ̅^2q( k, x)
       =xq^*( l, x) ƪ^2(ƪ̅^2q( k, x)w( x))
       ≃xq^*( l, x) (ƪ^2ƪ̅^2q( k, x))w( x)
       =1/β^2_k,2xq^*( l, x) q( k, x)w( x)
       =xq^*( l, x) (ƪ̅^2ƪ^2q( k, x))w( x)
       =xƪ̅^2(w( x) q^*( l, x)) ƪ^2q( k, x)
       ≃xw( x) (ƪ^2 q( l, x))^* ƪ^2q( k, x),

      where we have made repeated use of integration by parts.
    
      Using these identities together with the completeness relation of the basis functions, it is possible to simplify the expressions for the I^± s_a,± s_b:
      
    ^abI^0,0_ l l'= ^abI^0,+2_ l l'= ^abJ^0_ l l';  ^abI^+2,+2_ l l'= ^abJ^+_ l l';
        ^abI^+2,-2_ l l'= ^abI^-2,+2_ l l'= ^abJ^-_ l l';  ^abI^0,-2_ l l'= ^abI^-2,-2_ l l'=0,

      where  we have defined
      
    ^abJ^0_ l l'=x(w_a w_b)( x) q( l, x) [q( l', x)]^*,
        ^abJ^±_ l l'=x(w_a w_b)( x) q^±2_ l l'( x).

      Thus, the only surviving non-zero coupling coefficients are
      
    W^δδ,δδ_ l, l'=W^δγ_E,δγ_E_ l, l'=W^δγ_B,δγ_B_ l, l'=J^0_ l l';
       W^γ_Eγ_E,γ_Eγ_E_ l, l'=W^γ_Eγ_B,γ_Eγ_B_ l, l'=W^γ_Bγ_E,γ_Bγ_E_ l, l'=W^γ_Bγ_B,γ_Bγ_B_ l, l'=J^+_ l l'
       W^γ_Eγ_E,γ_Eγ_B_ l, l'=-W^γ_Eγ_B,γ_Eγ_E_ l, l'=W^γ_Bγ_E,γ_Bγ_B_ l, l'=-W^γ_Bγ_B,γ_Bγ_E_ l, l'=
       =-W^γ_Eγ_E,γ_Bγ_E_ l, l'=-W^γ_Eγ_B,γ_Bγ_B_ l, l'=W^γ_Bγ_E,γ_Eγ_E_ l, l'=W^γ_Bγ_B,γ_Eγ_B_ l, l'=-i J^-_ l l'.

    
      We have reduced the problem of computing the covariance in Eq. <ref> to the problem of computing the coupling coefficients <ref> entering Eq. <ref>, and we have now shown that there are only 3 independent coefficients, given by Eq. <ref>. In order to simplify the calculation further, it is now useful to inspect these results for the specific case of fields defined on the sphere.
    
      

  §.§.§ Covariances for curved skies

        As described in Section <ref>, in the curved-sky
        estimator it is common to first average over the harmonic number m
        as part of the bandpower binning operation. Let us, therefore, define
        
    C̃^a_α b_β_ℓ=1/2ℓ+1∑_m=-ℓ^ℓ (ã_α)_ℓ m (b̃_β)_ℓ m^*.

        The covariance of these objects, under the approximation in Eq. <ref>, can be computed as:
        
    Cov(C̃^a_α b_β_ℓ,C̃^f_ϕ g_γ_ℓ') = C^a_α' f_ϕ'_(ℓC^b_β' g_γ'_ℓ') ^af_bg𝒲^αϕ,α'ϕ'_βγ,β'γ'(ℓ,ℓ')+((g,γ)↔(f,ϕ)),

        where we have defined the symbols
        
    ^af_bg𝒲^αϕ,α'ϕ'_βγ,β'γ'(ℓ,ℓ')≡∑_m=-ℓ^ℓ∑_m'=-ℓ'^ℓ' ^afW^αϕ,α'ϕ'_ l l'( ^bgW^βγ,β'γ'_ l l')^*/(2ℓ+1)(2ℓ'+1).

        These quantities involve terms of the form:
        
    ^af_bg𝒥^XY_ℓℓ'≡∑_m=-ℓ^ℓ∑_m'=-ℓ'^ℓ' ^afJ^X_ l l'( ^bgJ^Y_ l l')^*/(2ℓ+1)(2ℓ'+1)

        Making use of the completeness relation for the Wigner 3-j symbols, it is possible to show, that these can be written as:
        
    ^af_bg𝒥^00_ℓℓ'=Ξ^00_ℓℓ'(w_aw_f,w_bw_g), ^af_bg𝒥^0+_ℓℓ'=Ξ^0+_ℓℓ'(w_aw_f,w_bw_g),
        ^af_bg𝒥^++_ℓℓ'=Ξ^++_ℓℓ'(w_aw_f,w_bw_g), ^af_bg𝒥^–_ℓℓ'=Ξ^–_ℓℓ'(w_aw_f,w_bw_g),
        ^af_bg𝒥^+-_ℓℓ'=
                ^af_bg𝒥^-+_ℓℓ'=
                ^af_bg𝒥^0-_ℓℓ'=0,

        where the Ξ^XY_ℓℓ' are defined in Eq. <ref>.
        Thus, computationally speaking, the problem of computing covariance matrices reduces to that of computing the same coupling coefficients needed for the computation of the pseudo-C_ℓ power spectra themselves, except now they involve product of two masks, rather than the masks alone. Using the same notation as in Eq. <ref>, the only non-zero 𝒲 coefficients are:
        
    𝒲^δδ  ,δδ_δδ  ,δδ=
               𝒲^δδ  ,δδ_δγ_X,δγ_X=
               𝒲^δγ_X,δγ_X_δγ_Y,δγ_Y=𝒥^00
       𝒲^δδ  ,δδ_γ_Xγ_Y,γ_Xγ_Y=
               𝒲^δγ_Z,δγ_Z_γ_Xγ_Y,γ_Xγ_Y=𝒥^0+
       𝒲^γ_Xγ_Y,γ_Xγ_Y_γ_Wγ_Z,γ_Wγ_Z=𝒥^++
       𝒲^γ_Xγ_E,γ_Xγ_B_γ_Yγ_E,γ_Yγ_B=
               𝒲^γ_Xγ_E,γ_Xγ_B_γ_Bγ_Y,γ_Eγ_Y=
               𝒲^γ_Bγ_X,γ_Eγ_X_γ_Bγ_Y,γ_Eγ_Y=
       𝒲^γ_Xγ_B,γ_Xγ_E_γ_Yγ_B,γ_Yγ_E=
               𝒲^γ_Xγ_B,γ_Xγ_E_γ_Eγ_Y,γ_Bγ_Y=
               𝒲^γ_Eγ_X,γ_Bγ_X_γ_Eγ_Y,γ_Bγ_Y=-𝒥^–
       𝒲^γ_Xγ_B,γ_Xγ_E_γ_Yγ_E,γ_Yγ_B=
               𝒲^γ_Xγ_B,γ_Xγ_E_γ_Bγ_Y,γ_Eγ_Y=
               𝒲^γ_Eγ_X,γ_Bγ_X_γ_Bγ_Y,γ_Eγ_Y=
       𝒲^γ_Xγ_E,γ_Xγ_B_γ_Yγ_B,γ_Yγ_E=
               𝒲^γ_Xγ_E,γ_Xγ_B_γ_Eγ_Y,γ_Bγ_Y=
               𝒲^γ_Bγ_X,γ_Eγ_X_γ_Eγ_Y,γ_Bγ_Y=𝒥^–,

        where (X,Y,Z,W) stand for either E or B, and where we have suppressed all redundant indices (including ℓℓ'). Any pseudo-C_ℓ covariance element can then be found by replacing these results in Eq. <ref>. Some explicit examples for common terms can be found in Appendix A.3.3 of <cit.>.

      

  §.§.§ Covariances for flat skies

        Similar results hold in the case of flat skies. As mentioned in
        Section <ref>, in this case averaging over the
        Fourier-space azimuth happens while binning into bandpowers. Under the
        assumption that the underlying power spectra are roughly constant
        within each bandpower, in this case, the covariance matrix of the bandpowers defined in Eq. <ref> takes the form:
        
    Cov(C̃^a_α b_β_q,C̃^f_ϕ g_γ_q') = C^a_α' f_ϕ'_(qC^b_β' g_γ'_q')∑_ l∈ q1/N_q∑_ l' ∈ q'1/N_q'  ^af_bg𝒲̅^αϕ,α'ϕ'_βγ,β'γ'( l, l')+((g,γ)↔(f,ϕ)),

        where the coefficients 𝒲̅ are related to the
        mode-coupling coefficients Ξ̅ defined in
        Eq. <ref>, in the same way that the curved-sky coefficients 𝒲 were related to the Ξ:
        
    𝒲̅^δδ  ,δδ_δδ  ,δδ=
               𝒲̅^δδ  ,δδ_δγ_X,δγ_X=
               𝒲̅^δγ_X,δγ_X_δγ_Y,δγ_Y=𝒥̅^00
       𝒲̅^δδ  ,δδ_γ_Xγ_Y,γ_Xγ_Y=
               𝒲̅^δγ_Z,δγ_Z_γ_Xγ_Y,γ_Xγ_Y=𝒥̅^0+
       𝒲̅^γ_Xγ_Y,γ_Xγ_Y_γ_Wγ_Z,γ_Wγ_Z=𝒥̅^++
       𝒲̅^γ_Xγ_E,γ_Xγ_B_γ_Yγ_E,γ_Yγ_B=
               𝒲̅^γ_Xγ_E,γ_Xγ_B_γ_Bγ_Y,γ_Eγ_Y=
               𝒲̅^γ_Bγ_X,γ_Eγ_X_γ_Bγ_Y,γ_Eγ_Y=
       𝒲̅^γ_Xγ_B,γ_Xγ_E_γ_Yγ_B,γ_Yγ_E=
               𝒲̅^γ_Xγ_B,γ_Xγ_E_γ_Eγ_Y,γ_Bγ_Y=
               𝒲̅^γ_Eγ_X,γ_Bγ_X_γ_Eγ_Y,γ_Bγ_Y=-𝒥̅^–
       𝒲̅^γ_Xγ_B,γ_Xγ_E_γ_Yγ_E,γ_Yγ_B=
               𝒲̅^γ_Xγ_B,γ_Xγ_E_γ_Bγ_Y,γ_Eγ_Y=
               𝒲̅^γ_Eγ_X,γ_Bγ_X_γ_Bγ_Y,γ_Eγ_Y=
       𝒲̅^γ_Xγ_E,γ_Xγ_B_γ_Yγ_B,γ_Yγ_E=
               𝒲̅^γ_Xγ_E,γ_Xγ_B_γ_Eγ_Y,γ_Bγ_Y=
               𝒲̅^γ_Bγ_X,γ_Eγ_X_γ_Eγ_Y,γ_Bγ_Y=𝒥̅^–,

        and, as before,
        
    ^af_bg𝒥̅^00_ l l'=Ξ̅^00_ l l'(w_aw_f,w_bw_g), ^af_bg𝒥̅^0+_ l l'=Ξ̅^0+_ l l'(w_aw_f,w_bw_g),
        ^af_bg𝒥̅^++_ l l'=Ξ̅^++_ l l'(w_aw_f,w_bw_g), ^af_bg𝒥̅^–_ l l'=Ξ̅^–_ l l'(w_aw_f,w_bw_g),
        ^af_bg𝒥̅^+-_ l l'=
                ^af_bg𝒥̅^-+_ l l'=
                ^af_bg𝒥̅^0-_ l l'=0.

        
    

 §.§ Approximate covariances

      When presenting our results in Section <ref>, we will compare the true covariance matrix, estimated from a large number of Gaussian simulations, with the analytical covariance estimated under different approximations. In descending order of complexity, these are:
      
        
  * The narrow-kernel approximation (labeled NKA here), described in the previous sections. This approximation assumes that the support of the harmonic-space masks (represented by the mode-coupling coefficients in e.g. Eq. <ref>) is small compared to the variation of the true power spectrum with ℓ. Additionally, it neglects all derivatives of the sky mask when accounting for the spin nature of the fields involved.
        
  * The spin-0 approximation corresponds to a simplified version of the NKA in which the spin nature of all fields involved is completely ignored, and all fields, including the E- and B-mode components of a spin-2 field, are treated as spin-0 quantities.
        
  * The mode-counting approximation (labeled MC here),
          commonly known as the Knox formula <cit.>,
          which applies the result found for full-sky observations (Eq.
          <ref>) to masked fields, corrected by an overall factor that accounts to the loss of modes due to the sky mask. In this case, the covariance is simply given by:
        
    Cov(C^ab_q,C^cd_q')=δ^K_qq'C^ad_qC^bc_q+C^ac_qC^bd_q/(2ℓ_q+1) f_ sky N_q,

        where ℓ_q is the mean multipole in the q-th bandpower, N_q is the number of multipoles assigned to it, and f_ sky is the available sky fraction.
      

  

§ RESULTS

    

 §.§ Simulations

      In order to quantitatively study the performance of the analytical approximations for the power spectrum covariance matrix introduced in Section <ref>, we generate a large number of Gaussian simulations including a number of realistic observational effects.
      
      Each simulation is a set of maps corresponding to a number of spin-0 and spin-2 fields that are drawn as Gaussian random fields following a set of input power spectra that include all relevant cross-correlations between different fields. We generate simulations for two types of fields modeled after the two main large-scale structure observables of photometric redshift surveys:
      
        
  * Spin-0 fields, corresponding to maps of the overdensity of galaxies within a given redshift bin projected on the sphere. In keeping with the notation introduced in Section <ref>, we will label these fields as δ^a, where the index a denotes the redshift bin.
        
  * Spin-2 fields, corresponding to maps of the cosmic shear measured from the projected shapes of galaxies in a given redshift bin. We will label these fields γ^a=(γ_E^a,γ_B^a) where, again, the index a denotes the redshift bin.
      
      The cross-correlation between two of these fields can be written as:
      
    C^ab_ℓ = S^ab_ℓ+N^ab_ℓ,

      where S and N are the power spectra of the cosmological signal and noise respectively. We model the signal part for δ and γ_E as:
      
    S^ab_ℓ = ∫ dχW^a_ℓ(χ) W^b_ℓ(χ)/χ^2 P(k=ℓ+1/2/χ,z),

      where χ is the comoving radial distance, z≡ z(χ) is the corresponding redshift in the lightcone, P(k,z) is the matter power spectrum. The window functions are given by <cit.>[See <cit.> and references therein for details about these calculations.]:
      
    W^δ_ℓ(χ)=b(z) H(z) p_z(z),
       W^γ_ℓ(χ)=f_ℓ 3H_0^2Ω_M/2(1+z)χ ∫ dz' p_z(z') χ(z')-χ/χ(z'),

      where H(z) is the expansion rate in units where the speed of light is c=1, H_0≡ H(z=0), b(z) is the linear galaxy bias, p_z is the normalized redshift distribution of galaxies within the redshift bin, and
      
    f_ℓ≡√((ℓ+2)(ℓ+1)ℓ(ℓ-1))/(ℓ+1/2)^2.

      We assume zero signal for the shear B-modes.
      
            
      
      

      The noise power spectrum is diagonal (i.e. zero between different fields and redshift bins), and is given by:
      
    N_ℓ^δδ = 1/n_Ω,
            N_ℓ^γ_Eγ_E = N_ℓ^γ_Bγ_B = σ_γ^2/n_Ω,

      where n_Ω is the mean number density of galaxies in units of sterad^-2 (see below), and σ_γ=0.28 is the intrinsic shape scatter per ellipticity component.

      We consider the case of auto- and cross-correlations between all fields in two redshift bins, with redshift distributions modeled as Gaussians with width σ_z=0.13 centered around redshifts 0.75 and 0.95, with number densities n_Ω=7.5  arcmin^-2. The corresponding redshift distributions p_z(z) are shown in Figure <ref>. For simplicity we use a constant bias b(z)=1. For these specifications, we generate signal power spectra using the Core Cosmology Library (CCL <cit.>) for cosmological parameters (Ω_c,Ω_b,h,A_s,n_s)=(0.27,0.045,0.67,2.1×10^-9,0.96). The resulting non-zero signal and noise power spectra are shown in Figure <ref>.

      After generating a set of Gaussian maps, we mask them making use of a realistic sky mask. This mask has three main components: a cut in declination based on the expected sky coverage of LSST <cit.>, a more conservative Galactic cut using the dust reddening data from <cit.> and a set of 100 randomly positioned holes with a radius of 1 degree. To explore the case of cross-correlations between fields with different masks, we generated two masks for the two redshift bins above, consisting of two different sets of random holes. These masks are shown in Figure <ref>.

      We have also explored the impact of the presence of sky contaminants on the estimate of the covariance matrix. As described in <cit.>, the presence of a small contamination from observational systematics in the data can be accounted for through a technique known as mode deprojection <cit.>. In this method, the contamination is modeled as a linear contribution at the map level from contaminants with a known template. Mode deprojection then consists on projecting the data onto the subspace of modes that are perpendicular to those templates, effectively removing all modes from the map that “look like” any of the contaminants. This removal of modes has been shown to provide unbiased estimates of the power spectrum, however it could potentially affect the power spectrum uncertainties due to the loss of statistical power. To study this effect, we have also generated contaminant maps of two types:
      
        
  * Large-scale contaminants: Gaussian random maps with a red spectrum of the form C_ℓ∝ (ℓ+1)^β, where β is a random number chosen within the range β∈(-1,-3).
        
  * Small-scale contaminants: Gaussian random maps with a flat spectrum C_ℓ= const.
      
      In both cases, we fixed the amplitude of the contaminant power spectrum such that it would yield a 10% contamination in the data power spectrum at ℓ=400. For spin-2 fields, we assumed the same power spectrum for E and B modes, with no cross-correlation between them. The resulting contaminant power spectra are shown in Figure <ref>. When exploring the effects of mode deprojection, we generated 100 contaminant maps of both types and added them to each simulated realization. We then deprojected the full set of 100 contaminant templates from the simulated maps and computed the corresponding unbiased power spectra. Note that, in what follows, our fiducial results do not include the effects of mode deprojection. These are discussed separately.

      For our fiducial results, we generated a set of 20,000 random
      simulations. This number was chosen in order to recover the covariance
      matrix for all possible auto- and cross-correlations with sufficient
      accuracy. All simulations were generated as HEALPix[<http://healpix.sourceforge.net>] <cit.> maps with resolution N_ side=512. We then used NaMaster to compute all possible power spectra for each simulation using narrow bandpowers of width Δℓ=3 from ℓ=2 to ℓ=1023. Finally, we used the power spectra from the simulations to estimate the sample covariance matrix
      
    Cov( C) =1/N_ sim-1∑_i=1^N_ sim( C_i-C̅)·( C_i-C̅)^T,

      where C_i is the vector of all possible power spectra for the i-th realization, and C̅ is the mean of this vector over all realizations. The comparison of this sample covariance with the analytical approximations described above is presented in the next sections.

    

 §.§ Qualitative comparison

      
      
      
      As a first step, we visually compare the main properties of the sample covariance matrix estimated from the simulations and the NKA, Spin-0 and MC approximations described in the previous section. Figure <ref> shows four rows of the covariance matrix of different auto- and cross-correlations. The rows correspond to bandpowers centered on multipoles ℓ_q≃60, 90, 120 and 150. The upper panels show results for the non-zero power spectra (δ-δ, δ-γ_E and γ_E-γ_E), with the solid blue, dashed orange and dotted green lines showing results for the sample covariance matrix, the analytical covariance using the NKA approximation and their difference, respectively. For comparison, the black stars show the diagonal covariance matrix elements predicted by the MC approximation (Eq. <ref>). We find an excellent agreement between the simulated and analytical covariances, with very small deviations in the amplitude of the diagonal and first few off-diagonal elements.
      
      The bottom panel in Figure <ref> shows the same rows of the covariance matrix for power spectra involving B-modes (and therefore with zero signal expectation value). In this case we find significant differences, at the level of 30-50%, on the covariance matrix elements, with the analyticial prediction underestimating the error bars overall. This is expected and can be understood as follows: the presence of a sky mask mixes E and B modes. Although this mixing can be accounted for at the level of the power spectrum through the pseudo-C_ℓ estimator, the leaked modes contribute to the variance. This is particularly significant for power spectra involving B-modes, since the E-mode amplitude is significantly larger, especially at ℓ≲200, as can be seen in Fig. <ref>. Thus, if the effects of E-B mixing caused by the sky mask are not accurately accounted for in the estimation of the covariance matrix for power spectra involving B-mode maps, we can expect a misestimation of the contribution to the covariance from the leaked E modes that would underpredict the uncertainties. This is not a problem for power spectra involving only E modes, since the only B modes that leak into them are those associated with noise, and they have the same amplitude as the noise E modes move into the B-mode map.

      This is further illustrated by Figure <ref>. The figure shows, for the case of a single redshift bin, the difference between the correlation matrices associated with the sample covariance matrix and the NKA estimate[The correlation matrix is defined as r_ij =  Cov_ij/√( Cov_ii Cov_jj)]. While the differences between both matrices are small for all elements involving δ and γ_E, all terms involving B-modes show a significant disagreement, particularly the  γ_Eγ_E-γ_Bγ_B, γ_Eγ_B-γ_Eγ_B and γ_Bγ_B-γ_Bγ_B boxes.

      We thus conclude that while the NKA estimator is able to recover the covariance matrix for the non-zero power spectrum elements (i.e. those involving δ and γ_E) with high accuracy, a more sophisticated approach would be needed in order to obtain a precise estimate of the uncertainties for components involving B-modes. This is not a major concern, since B-mode power spectra are predominantly used as null tests, while cosmological parameter constraints are driven by the analysis of δ and γ_E. For completeness, Figure <ref> shows the difference between the correlation matrices for the sample covariance and the NKA estimator for all non-zero cross-correlations between different bins in the case of two redshift bins with different small-scale masks, where we find a similarly good agreement.

    

 §.§ Quantitative comparison

      
      In order to quantify the validity of our analytical approximations, we need to compare the NKA covariance with the sample covariance estimated from the simulations. However, comparing two matrices is not as straightforward as comparing their elements one by one. The covariance between far-away bandpowers is expected to be very close to zero, and therefore a direct comparison of those elements would easily yield large relative differences simply due to the statistical noise in the sample covariance matrix. We will therefore quantify the differences between the different covariances making use of scalar quantities formed from them. The impact of the analytical approximations on the final parameter constraints will then be described in detail in Section <ref>.
      
      As a first test to quantify the differences between covariance matrices, we compute the relative difference between their eigenvalues. This is shown in Figure <ref> for a data vector combining all auto- and cross-correlations between δ and γ_E for a single redshift bin. The eigenvalues of both matrices are roughly similar, with relative differences of about 5%. The figure also shows the eigenvalues of the covariance matrix estimated using the Spin-0 approximation, which achieves a similar level of precision (even marginally higher in some cases).

      Another scalar quantity that can be used to compare different covariances is the χ^2. For a given random data vector d with mean m and covariance matrix Cov, this is given by
      
    χ^2=( d- m)^T· Cov^-1·( d- m).

      
      
      We compute this quantity for a data vector d composed of different auto- and cross-correlations for each of the 20,000 Gaussian simulations, with m given by the mean over all simulations and different choices of covariance matrix. Figure <ref> shows the distribution of χ^2 values for the three non-zero power spectra in the case of a single bin: δ-δ, δ-γ_E and γ_E-γ_E. The histograms show the distribution for the sample covariance matrix (blue), and the analytical NKA and Spin-0 estimators (orange and green, respectively). We additionally plot the theoretical χ^2 distribution under the assumption that the underlying data vector is Gaussianly distributed (red dashed lines). In the simplest case of purely spin-0 quantities (leftmost panel), we find an excellent agreement between the different distributions. In the cases involving the spin-2 fields, we see noticeable differences between the distributions found with the sample covariance and the approximate ones. These differences are small, corresponding to less of a 2 and 4% shift in the mean χ^2 for the δ-γ_E and γ_E-γ_E cases, respectively, and a negligible variation in the width of the distributions. We therefore expect these differents to have a negligible effect on the posterior parameter distributions, as we show explicitly in Section <ref>. The fact that these differences appear only for power spectra involving spin-2 quantities indicate that the NKA and spin-0 methods are imperfect at describing the additional mode coupling caused in the presence of a mask for higher-spin fields. This is rather obvious in the case of the spin-0 approximation but, interestingly, we find that the NKA and Spin-0 predictions yield results that are almost indistinguishable from each other. We therefore conclude that the additional approximation made in the NKA method for spin-2 fields – neglecting the spatial derivatives of the mask – is effectively equivalent to ignoring the spin nature of the fields involved. Note, however, that this is not the case for B-modes, where the NKA estimator outperforms the spin-0 approximation by up to one order of magnitude, even though its accuracy is very poor (as we described in the previous section). For completeness, Figure <ref> shows the distribution of χ^2 values for a data vector composed of all possible auto- and cross-correlations of δ and γ_E for the case of two redshift bins, where similar conclusions hold.

      We have also explored the impact of contaminant deprojection on the different covariance matrix estimates. The loss of modes due to deprojection can potentially increase the variance of the power spectrum estimates, affecting the accuracy with which an analytical estimator would be able to recover the covariance. Figure <ref> shows the diagonal of the covariance matrix for the δ-δ,      δ-γ_E and γ_E-γ_E power spectra. Each panel displays the diagonal for the sample covariance matrix estimated from simulations without contaminants or contaminant deprojection (blue
      line), the sample covariance from simulations with contaminants and contaminant deprojection (dashed orange line) and for the NKA covariance   (green line). We see that the power spectrum uncertainties are almost indistinguishable with or without deprojection, and that those relative differences are much smaller than the differences between the sample covariance and the NKA estimator. We therefore conclude that, except in the case where a very large set of contaminant maps are deprojected (comparable with the number of unmasked pixels in the map), the analytical approximation to the covariance matrix should be as accurate as in the absence of contaminants (i.e. accurate enough).

      
      Figure <ref> serves also to illustrate another important
      point. In all our tests we find that the largest differences between the
      sample and NKA covariances occur on large scales. For δ-δ
      correlations, the effect is limited to the first few multipoles
      (ℓ≲10), while for γ_E-γ_E we are only able to
      recover the sample covariance errors within 5% for ℓ≳40.
      This is the main source of the small mismatch observed in
      Figure <ref>. Note, however, that most of the cosmological
      information is obtained from the higher multipoles, due to their higher
      statistical weight, and, as we will show in the next section, the effect
      on the final parameter constraints is negligible. If accurate
      covariance matrix elements are needed on these large scales for cosmic shear, they
      can be estimated alternatively making use of fast, low resolution
      simulations (e.g. HEALPix N_ side=64 maps), or computed
      exactly as in Ref. <cit.>.

    

 §.§ Impact on parameter estimation

      
      Ultimately, the most important test to judge the accuracy of the
      analytical covariance matrix estimators implemented here is to study
      their impact on the posterior distribution of cosmological parameters
      derived from power spectrum measurements. Assuming flat priors and
      a fiducial Gaussian likelihood approximation <cit.>, in which the covariance matrix is computed only once for the fiducial model, the posterior distribution for parameters θ⃗ is simply given by:
      
    -2log p(θ⃗|Ĉ)=(Ĉ- C(θ⃗))^T· Cov^-1·(Ĉ- C(θ⃗))+ const.,

      where Ĉ is a vector of power spectrum measurements, Cov is their covariance matrix, and C(θ⃗) is their theoretical prediction for parameters θ⃗. It is worth noting that this Gaussian likelihood is not accurate on large scales, where the small number of modes invalidates the application of the central limit theorem. Even in this regime, the likelihood can be approximated through the method described in <cit.>, which still requires an accurate estimate of the power spectrum covariance.

      We explore log p(θ⃗|Ĉ) for the two parameters θ⃗≡(Ω_M,σ_8), for a data vector composed of all possible auto- and cross-correlations between δ and γ_E in the case of two redshift bins described in Section <ref>. In this simple two-dimensional scenario, we simply sample the distribution in a regular grid of 100 by 100 points for each parameter. We construct a data vector Ĉ from the theoretical prediction for the experimental setup described in Section <ref>, and produce theoretical predictions for it at each grid point using CCL. Note that, when evaluating the posterior for the sample covariance matrix, one needs to correct for the finite number of simulations used to construct the covariance. In most situations this can be done simply by rescaling the inverse covariance matrix by a factor given by <cit.>
      
    Cov^-1 ⟶N_s-2-N_ data/N_s-1  Cov^-1,

      where N_s=20,000 is the number of samples used to estimate the
      covariance, and N_ data=3510 is the number of data points.
      
      The 68% and 95% confidence level contours associated with the posterior distributions for the sample covariance and the NKA estimator are shown in Figure <ref>. We find that both distributions agree with each other remarkably well, and that the 1σ errors for each parameter agree for both covariances up to 0.3%. Note that, since we have not included any statistical noise in the data vector Ĉ, the relative difference in the means of both distributions is zero by construction (since Eq. <ref> is bounded from below by zero). When adding Gaussian statistical noise compatible with the sample covariance matrix, we observe small differences (smaller than 0.3σ) in the best fit parameters found with both covariance matrices. These differences, however, are not systematic, and we have verified that averaging over several noise realizations does not yield biased best-fit parameters.
      
      We therefore conclude that the analytical approximations for the power spectrum covariance matrix explored here are able to reproduce the true posterior distribution of cosmological parameters to very high accuracy.

  

§ DISCUSSION

    Estimating accurate covariance matrices for projected two-point
    correlators is an ubiquitous problem in modern cosmology
    <cit.>, but it
    is particularly relevant for large-scale structure datasets. This problem
    is further complicated in this case, in comparison with e.g. CMB
    experiments, by two factors: 
    the fact that the fields involved are
    non-Gaussian at some level and the arguably higher complexity of the sky masks used in optical datasets.
    The main impact of survey geometry is the statistical
    coupling it induces between different Fourier/harmonic modes, which must
    be accurately characterized in order to obtain reliable estimates of the
    posterior distribution for cosmological parameters. In this paper we have
    focused on the impact of survey geometry on the dominant Gaussian (i.e.
    disconnected) part of the covariance matrix.
    
    We have described and generalized existing analytical approaches to
    estimate the covariance matrix for pseudo-C_ℓ power spectrum
    estimators
    <cit.>, and
    implemented them in the public code NaMaster
    <cit.>, making it straightforward to fully account for
    the effects of survey geometry on the data uncertainties. With these
    approximations, computationally speaking, the problem of estimating a
    covariance matrix is as complex as that of computing the power spectrum
    itself, and scales with the number of pixels in the map as N_
      pix^3/2. We leave for future work the study of the impact of
      extra complications, as position-dependent noise, on their performance.
    
    The main finding of this paper is the excellent performance of the
    analytical methods described in Section <ref>. We have
    shown that
    
    the NKA estimator (see Section <ref>) is able to
    recover the covariance matrix for all power spectra with non-zero signal
    expectation value (i.e. those involving the galaxy overdensity δ or
    the E-mode shear γ_E), as well as the posterior distribution for
    cosmological parameters, to a high degree of accuracy.

    More in detail, we have also found that the impact of contaminant
    deprojection on the covariance matrix, through the corresponding loss of
    modes, is negligible unless a very large number of contaminant templates
    are removed. This simplifies the procedure to estimate covariance matrices
    for galaxy clustering data, which are particularly sensitive to a large
    number of astrophysical and observational systematics. Additionally, we
    have found that, although the NKA estimator is accurate enough, it is not
    able to perfectly capture the additional effects of mode coupling that are
    present for spin-2 fields, and that a simpler approach treating the shear
    E modes as a spin-0 object is able to reach similar levels of accuracy.
    Due to the imperfect treatment of the E/B mixing caused by the sky mask
    in the NKA estimator, we also find that the predicted covariance matrix
    for any power spectra involving B modes differs significantly from the
    true sample covariance, and, therefore, this approach cannot be used to
    reliably estimate the uncertainties of B-mode power spectra. Likewise,
    we find that the NKA estimator yields inaccurate estimates of the power
    spectrum uncertainties on the largest scales (ℓ≲50). Although
    these modes carry a substantially smaller statistical weight, if more
    accurate covariances are needed on these large scales, they can be easily
    computed making use of fast low-resolution simulations or exactly
    following Refs. <cit.>.
    
    In spite of these shortcomings, we find that the approximations described in this paper are able to provide estimates of the power spectrum covariance matrix that are sufficiently accurate for current and future tomographic large-scale structure cosmological datasets. The main advantage of this approach is the computational cost, which is comparable to that of estimating the power spectrum in the first place, and which  scales with pixel resolution in a similar way. This is, therefore, significantly less time-consuming than generating large numbers of mock datasets (even simple Gaussian or log-normal realizations), and more reliable than the traditional jackknife resampling techniques. The method, in all its generality, is currently implemented in the public code NaMaster. Future extensions to this work will focus on improving the estimator for power spectra involving B modes, and incorporating the impact of E/B-mode purification <cit.>.
    
  We thank Thibaut Louis, Eva-Maria Mueller, Francisco Javier Sánchez and
  Anže Slosar for useful comments and discussion. CGG is supported the
  Spanish grant BES-2016-077038, partially funded by the ESF and by
  AYA2015-67854-P from the Ministry of Industry, Science and Innovation of
  Spain and the FEDER funds. He was partially supported by a Balzan Fellowship
  while in Oxford. He would like to thank New College and the Department of
  Physics at Oxford for their hospitality. DA acknowledges support from STFC
  through an Ernest Rutherford Fellowship, grant reference ST/P004474/1. EB
  acknowledges support from the Beecroft Trust.  Some of the results in this
  paper have been derived using the HEALPix <cit.> package.
  
  

§ FLAT SKY

    

    
    We have repeated the analysis described in Section <ref> on flat-sky realizations, making use of the flat-sky implementation of NaMaster. We generate Gaussian realizations of the galaxy overdensity and shear maps making use of flat-sky extensions of the methods described in Section <ref>. In this case we use a high-resolution mask constructed from the bright-object mask distributed with the first data release of the HSC collaboration <cit.> for the VVDS field, which is shown in Figure <ref>. 
    
    We find similar levels of accuracy in the NKA and spin-0 estimators compared to the curved-sky case. The higher resolution and smaller area of these simulations allow us to focus on the small-scale galaxy clustering and lensing power spectra, covering the range of multipoles ℓ∈(120,17640) in constant bandpowers of width Δℓ=240. On these small scales, the shear power spectrum is more dominated by noise (e.g. see Fig. <ref>), and therefore there is roughly the same power in E and B modes. This reduces the sensitivity of the method to an inaccurate treatment of E/B leakage, and the agreement between the sample covariance matrix and the NKA estimator improves significantly. This can be seen in Figure <ref>, which shows the χ^2 distributions for all possible power spectra (including those including B modes) in the case of a single redshift bin. In all cases we find a good agreement between the χ^2 distributions derived from all covariance matrix estimates (sample covariance, NKA and spin-0 in blue, orange and green respectively), which also accurately follow the expected χ^2 distribution for the corresponding number of degrees of freedom (red dashed line).
    
  

§ SOFTWARE IMPLEMENTATION

    In addition to the code functionality described in Section 3 of <cit.>, we have now included the capability to estimate Gaussian covariance matrices using the NKA method. This functionality is structured around a python class called NmtCovarianceWorkspace. These objects are used to compute and store the covariance mode-coupling coefficients in Eqs. <ref> and <ref>. They are initialized from two pairs of fields corresponding to the two power spectra for which the covariance is required. Once initialized, these coefficients can be reused for any other set of fields with the same combination of sky masks. NaMaster then provides routines to estimate covariance matrix elements making use of the coupling coefficients stored in a NmtCovarianceWorkspace object and best-guess power spectra for the fields involved using Eqs. <ref> and <ref>. Further details about the implementation and practical examples can be found in <https://namaster.readthedocs.io/en/latest/sample_covariance.html>.

  JHEP
  

