# It takes two to know one: computing accurate one-point PDF covariances from effective two-point PDF models

Cora Uhlemann, Oliver Friedrich, Aoife Boyle, Alex Gough, Alexandre

Barthelemy, Francis Bernardeau, Sandrine Codis

► To cite this version:

Cora Uhlemann, Oliver Friedrich, Aoife Boyle, Alex Gough, Alexandre Barthelemy, et al.. It takes two to know one: computing accurate one-point PDF covariances from effective two-point PDF models. Open Journal of Astrophysics, 2023, 6, 10.21105/astro.2210.07819 . cea-04240674

HAL Id: cea-04240674

https://cea.hal.science/cea-04240674

Submitted on 13 Oct 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
IT TAKES TWO TO KNOW ONE: COMPUTING ACCURATE ONE-POINT PDF COVARIANCES
FROM EFFECTIVE TWO-POINT PDF MODELS

CORA UHLEMANN<fnref n="1" />,*, OLIVER FRIEDRICH<fnref n="2" />, AOIFE BOYLE<fnref n="3" />, ALEX GOUGH<fnref n="1" />, ALEXANDRE BARTHELEMY<fnref n="2" />,
FRANCIS BERNARDEAU<fnref n="4" />,<fnref n="5" />, AND SANDRINE CODIS<fnref n="3" />

Version January 4, 2023

Abstract

One-point probability distribution functions (PDFs) of the cosmic matter density are powerful cosmological probes that extract non-Gaussian properties of the matter distribution and complement two-point statistics. Computing the covariance of one-point PDFs is key for building a robust galaxy survey analysis for upcoming surveys like Euclid and the Rubin Observatory LSST and requires good models for the two-point PDFs characterising spatial correlations. In this work, we obtain accurate PDF covariances using effective shifted lognormal two-point PDF models for the mildly non-Gaussian weak lensing convergence and validate our predictions against large sets of Gaussian and non-Gaussian maps. We show how the dominant effects in the covariance matrix capturing super-sample covariance arise from a large-separation expansion of the two-point PDF and discuss differences between the covariances obtained from small patches and full sky maps. Finally, we describe how our formalism can be extended to characterise the PDF covariance for 3D-dimensional spectroscopic fields using the 3D matter PDF as an example. We describe how covariances from simulated boxes with fixed overall density can be supplemented with the missing super-sample covariance effect by relying on theoretical predictions validated against separate-universe style simulations.

# 1. INTRODUCTION

The recent rise of non-Gaussian statistics for galaxy clustering and weak lensing has highlighted the need for accurate models of the likelihood function and, in particular, the covariance matrices of such statistics. Aside from the challenge to model non-Gaussian statistics and capture their response to changing cosmological parameters in the nonlinear regime, getting accurate covariance matrices can prove difficult and require thousands of numerical simulations to ensure convergence (Taylor & Joachimi 2014; Colavincenzo et al. 2019). For statistics relying on higher-order correlation functions involving *N* points, modelling the covariance generally requires assumptions for the shape of 2*N*-point correlations. Even for Gaussian fields, this task has only been achieved recently at all orders *N* (Hou et al. 2021). To go beyond the assumption of Gaussianity, log-normal models of the cosmic density field have become a tool of choice in higher-order analyses such as studies of the bispectrum (Martin et al. 2012; Halder et al. 2021), moments of cosmic random fields (Gatti et al. 2019), density split statistics (Friedrich et al. 2018; Gruen et al. 2018), the full probability distribution function (PDF) of cosmic random fields (Uhlemann et al. 2020; Friedrich et al. 2020; Boyle et al. 2021), and even map-based inference approaches (Sarma Boruah et al. 2022). The reason why (shifted) log-normal simulations (Xavier et al. 2016) have become so popular is that they obey a similar hierarchy between variance and skewness as the physical cosmic density field. In particular, they can be tuned to match both a desired power spectrum at all scales and a desired

skewness at one given scale (c.f. Section 5.2 and Figure 4 of Friedrich et al. 2020). This however also highlights a limitation of those simulations: they potentially fail to realistically capture the hierarchy of moments beyond the skewness (for an extension see Baratta et al. 2020) as well as the full shape of higher order N-point functions.

With the large data vectors that future multi-probe analyses of upcoming surveys like Euclid (Laureijs et al. 2011; Amendola et al. 2018) and the Rubin Observatory LSST (Ivezić et al. 2019) are aiming for, estimated covariances - whether from log-normal simulations or other mock data - can lead to a huge degradation of cosmological constraining power (Dodelson & Schneider 2013a; Friedrich & Eifler 2017; Percival et al. 2022). This makes it desirable to have an analytic understanding of measurement uncertainties. One-point statistics such as the PDF of the weak lensing convergence have the advantage that both their dependence on cosmological parameters and their covariance, being given in terms of the two-point PDF, can be modelled theoretically (Codis et al. 2016b; Uhlemann et al. 2017). In this work, we will focus on the weak lensing convergence smoothed on mildly nonlinear scales (of order $10^\prime$ at source redshifts of $z \simeq 1-2$), where the bulk of the PDF and its cosmology dependence can be predicted at the percent level using large deviation statistics (Bernardeau & Valageas 2000; Bernardeau & Reimberg 2016; Barthelemy et al. 2020a; Boyle et al. 2021) and its covariance is well-approximated by a shifted lognormal distribution (Boyle et al. 2021). In the opposite regime for small smoothing scales, the high-convergence tail of the weak lensing convergence PDF and its covariance have been modelled using a halo-model formalism in Thiele et al. (2020), finding qualita-

<footnote n="1">¹ School of Mathematics, Statistics and Physics, Newcastle University, Herschel Building, NE1 7RU Newcastle-upon-Tyne, U.K.</footnote>

<footnote n="2">² Universitäts-Sternwarte, Fakultät für Physik, Ludwig-Maximilians Universität München, Scheinerstr. 1, 81679 München, Germany</footnote>

<footnote n="3">³ AIM, CEA, CNRS, Université Paris-Saclay, Université Paris Diderot, Sorbonne Paris Cité, 91191 Gif-sur-Yvette, France</footnote>

<footnote n="4">⁴ Université Paris-Saclay, CNRS, CEA, Institut de physique théorique, 91191, Gif-sur-Yvette, France and</footnote>

<footnote n="5">⁵ CNRS & Sorbonne Université, UMR 7095, Institut d'Astrophysique de Paris, 75014, Paris, France</footnote>

<footnote>*cora.uhlemann@newcastle.ac.uk</footnote>
tive agreement with simulations but significant quantitative discrepancies attributed to strong sensitivity to small-scale effects.

For galaxy clustering, the error on and correlation among the factorial moments of galaxy counts in cells $F_k = \langle N(N-1)\dots(N-k) \rangle$ have been computed for hierarchical models in Szapudi & Colombi (1996) and Szapudi et al. (1999). Those studies included cosmic errors due to a finite volume and discreteness errors due to a finite number of tracers relevant for clustering but not the weak lensing convergence, which is instead affected by shape noise due to shear measurements. In Repp & Szapudi (2021a), the variance and covariance among galaxy counts in cells was obtained for the case of an averaged cell correlation, which has previously been considered in Codis et al. (2016b). As we explain here, and derive in Appendix B, the large-separation expansion in terms of bias functions adopted in Codis et al. (2016b); Uhlemann et al. (2017); Repp & Szapudi (2021a); Repp & Szapudi (2021b) can be generalised to the physically more realistic case of a separation-dependent correlation. The large-separation two-point PDF and the resulting super-sample covariance effect can also be related to the idea of the position-dependent matter PDF studied in Jamieson & Loverde (2020). Recent work of Bernardeau (2022) determined covariances of density PDFs in hierarchical models, which can be used to understand the dominant terms in the covariance matrix of cosmological density fields in the strongly non-Gaussian regime. In this work, we build on this foundation but focus on the weakly non-Gaussian weak lensing PDF for which hierarchical models are less suitable and better results can be achieved using shifted lognormal models.

## STRUCTURE

In Section 2, we discuss the general covariance modelling for the two-dimensional case of the weak lensing (or photometric galaxy clustering) PDF based on a model for the joint two-point PDF. In Section 3, we use a large-separation expansion for this two-point PDF and link it to dominant terms in an eigendecomposition of the covariance matrix. In Section 4, we validate our covariance models using several quantitative tests. In Section 5, we present the basics for generalising our findings to the case of the three-dimensional matter (or spectroscopic galaxy clustering) PDF. We conclude in Section 6 and provide an outlook for how our results can be applied to related one-point observables, include tomography and information from two-point statistics.

## EXECUTIVE SUMMARY

We develop a theoretical model for the covariance of measurements of the one-point weak lensing convergence PDF on mildly nonlinear scales in finite bins. Our model is based on an integration over the two-point PDF (6) at different separations and we rely on a joint shifted lognormal model for that two-point PDF (28) tuned to match the PDF skewness and shown to capture the density-dependence of the two-point correlation (12) in simulations, as illustrated in Figure 8. We establish expressions for the finite sampling (11) and super-sample covariance (13) contributions to the PDF covariance and their scaling with the inverse of the survey area. We

<figure><img src="image_2.png" /><figcaption>FIG. 1.— PDF covariance matrix for circular apertures of radius $\theta_s = 7.5'$ at source redshift $z_s = 2$ comparing the measurement from 108 full sky Takahashi lensing maps (upper triangle) to the shifted lognormal prediction (lower triangle).</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>FIG. 2.— Fisher forecasts validating the predicted shifted lognormal PDF covariance (purple) against measured covariances from 108 realistic full sky simulated maps (red) and 500 associated mock FLASK maps (green). The parameter errors obtained from the raw measured covariances have additionally been corrected by the Hartlap factor (72). To highlight the variation that can result from small samples, we also show contours for two subsamples of 100 of the 500 mock FLASK maps (green dotted).</figcaption></figure>

validate our predictions against large sets of simulated convergence maps both from N-body simulations (Takahashi et al. 2017) and using the FLASK tool for generating log-normal random fields (Xavier et al. 2016). To verify the quality of our covariance model we employ a range of tests including visual comparisons, an eigen-decomposition of the covariance, a band-decomposition
of the precision matrix, $\chi^2$ tests and Fisher forecasts, with Figures 1 and 2 summarising two main results. In Figure 1 we show the measured covariance for the central region of the binned $\kappa$-PDF for a mildly nonlinear smoothing radius of $\theta_s = 7.5'$ at source redshift $z_s = 2$ from all the $N_{sim} = 108$ simulated full sky $\kappa$ maps from Takahashi et al. (2017) (upper triangle) along with the shifted lognormal prediction (lower triangle, using equations (6) and (28)). For the predictions, we used the cosmological parameters of the simulation along with the shift parameter $s \simeq 0.085$ obtained from the second and third moments measured from the PDF. In Figure 8 in the main text we show that the environment-dependent clustering that enters the eigendecomposition of the covariance is well captured by the shifted lognormal model.

Aside from the good visual agreement, we quantify their consistency using two tests. First, we test the compatibility of the diagonal elements of the two covariances shown in Figure 1 with a $\chi^2$ test, taking into account their uncertainty and correlation among each other<fnref n="1" /> (covariance of the covariance). We use the covariance between the numerically estimated diagonal covariance elements $\hat{C}_{ii}$ and $\hat{C}_{jj}$ as given in equation (36) following Taylor et al. (2013), where $N_{sim} = 108$ is the number of simulations used to estimate and the covariance elements $C_{ij}$ are evaluated with our theoretical model. Using the above expression to compute a $\chi^2$ between the diagonal of the estimated and predicted covariance we find $\chi^2 = 36.8$ for a number of 26 bins. This corresponds to a marginal $\sim 1.75\sigma$ deviation between the two matrix diagonals.

Finally, we assess the quality of the off-diagonal elements of our covariance model by performing a Fisher forecast. Figure 2 compares Fisher contours on the parameters $\sigma_8$ and $\Omega_{cdm}$ derived from a number of different covariance matrices: our covariance model (purple), the covariance estimated from the 108 N-body simulations of Takahashi et al. (2017) (red) and the covariance estimated from 500 log-normal simulations (green). For the simulations we use the maximum number of full sky maps available ($N_{sim} = 108$) which is not ideal for a PDF data vector of length 26. To mitigate this limitation, we include the Hartlap correction factor (72) for the Fisher forecasts. Additionally, we estimated the covariance from 500 mock full sky maps generated from FLASK and checked compatibility of the results. The differences in the width and alignment of the constraints obtained from these matrices is negligible compared to the overall parameter uncertainties. Hence, our covariance model is sufficiently accurate for application to observational data. We demonstrate the flexibility of our covariance model by adapting it to predict super sample covariance effects for the 3-dimensional matter PDF paving the way to spectroscopic galaxy clustering.

## 2. WEAK LENSING PDF COVARIANCE

The weak lensing convergence PDF is a well-suited complementary summary statistic that can extract non-Gaussian information from cosmic shear observations, thus complementing two-point statistics. In this work,

we focus on mildly nonlinear scales, where weak lensing convergence PDF and its response to fundamental cosmological parameters can be predicted at percent-level accuracy from large-deviation statistics as shown in Barthelemy et al. (2020a); Boyle et al. (2021). In the same work it was shown how one can use shifted log-normal convergence fields generated by FLASK (Xavier et al. 2016) to match the convergence PDF and obtain precise covariance estimates in good agreement with maps from N-body simulations, of which often only a very limited number is available. Following this work, we will study the covariance of the central region of the one-point PDF of the weak lensing convergence smoothed with a top-hat filter of radius $\theta_s = 7.3'$ at a fixed source redshift $z_s = 2$ using 26 linearly spaced bins in the range $-0.024 < \kappa < 0.026$. In this Section we will use both Gaussian and shifted lognormal fields generated with FLASK to test our theoretical framework for predicting the weak lensing PDF covariance based on different models for the two-point PDF of the weak lensing convergence in cells of a given separation. While we focus on the two-dimensional case of the weak lensing PDF here, the formalism can be adapted to describe the three-dimensional clustering PDF as we lay out in Section 5.

### 2.1. Covariance prediction from joint 2-point PDF

Before calculating the covariance matrix of a PDF measurement, let us first specify how this measurement is conducted. Let us assume that we have observed the gravitational lensing convergence field on a part of the sky with the area $A_{survey}$. Now let us smooth this convergence field with a circular aperture of some angular radius $\theta_s$ (potentially cutting away parts of the edges of the survey, so that each aperture falls fully inside the survey area), and let us pixelise the resulting map. Each pixel of that map contains the value of the average convergence within a radius $\theta_s$ around the pixel center.

We would now like to estimate the probability density of a certain value $\kappa$ for the smoothed convergence field. In practice one will need to compute this PDF within a set of finite bins, say around a set of central values $\kappa_i$ and with constant bin width $\Delta\kappa$, which results in the bins $\text{bin}_i = [\kappa_i - \Delta\kappa/2, \kappa_i + \Delta\kappa/2]$. The PDF at each $\kappa_i$ can then be approximated as

$$ \mathcal{P}(\kappa_i) \approx \frac{P(\kappa \in [\kappa_i - \Delta\kappa/2, \kappa_i + \Delta\kappa/2])}{\Delta\kappa}, \quad (1) $$

where $\mathcal{P}(\kappa)$ is the PDF (obtained in the limit $\Delta\kappa \to 0$) and $P(A)$ is the probability of statement $A$. With our pixelised map we can estimate the right-hand side of this to obtain the measurement

$$ \hat{\mathcal{P}}(\kappa_i) \equiv \frac{\#\{\text{pix with } |\kappa - \kappa_i| \le \Delta\kappa/2\}}{\#\{\text{pix}\}\Delta\kappa}. \quad (2) $$

Following Bernardeau (2022), we introduce the weight function

$$ \chi_i(\kappa) = \begin{cases} 1 & \text{if } \kappa \in [\kappa_i - \Delta\kappa/2, \kappa_i + \Delta\kappa/2] \\ 0 & \text{else} \end{cases}, \quad (3) $$

such that the estimator $\hat{\mathcal{P}}(\kappa_i)$ is given by

$$ \hat{\mathcal{P}}(\kappa_i) = \frac{\sum_P \chi_i(\kappa_P)}{N_P \Delta\kappa} \approx \frac{1}{A_{\text{survey}}} \int_{\text{survey}} d^2\theta \frac{\chi_i(\kappa(\theta))}{\Delta\kappa}, \quad (4) $$

<footnote n="1">¹ Note that typically only the diagonal elements of estimated covariances are well described with multi-variate Gaussian noise.</footnote>
<figure><img src="image_4.png" /><figcaption>FIG. 3.— Comparison of the continuous distribution of relative distances in a square Cartesian survey patch of side length $\Theta_{side} = 5\text{deg} = 300'$ (black) and the discrete version for a coarse grid corresponding to non-overlapping circular cells of radius 7.3' (blue). We also show the distribution of relative distances on the full sky with maximum $\Theta_{side} = \pi = 180\text{deg}$ (red).</figcaption></figure>

where P labels the $N_P = A_{survey}/A_{pix}$ pixels of area $A_{pix}$, $\theta$ are angular coordinates on the sky, and $\kappa(\theta)$ is the smoothed, continuous convergence field at that location. The second moment of our above estimator is given by

$$ \begin{aligned} & \langle \hat{P}(\kappa_i) \hat{P}(\kappa_j) \rangle \\ & \approx \frac{1}{A_{\text{surv.}}^2} \int_{\text{surv.}} d^2\theta_1 \int_{\text{surv.}} d^2\theta_2 \frac{\langle \chi_i(\kappa(\theta_1)) \chi_j(\kappa(\theta_2)) \rangle}{(\Delta\kappa)^2} \\ & \approx \frac{1}{A_{\text{surv.}}^2} \int_{\text{surv.}} d^2\theta_1 \int_{\text{surv.}} d^2\theta_2 \mathcal{P}(\kappa_i, \kappa_j; |\boldsymbol{\theta}_1 - \boldsymbol{\theta}_2|), \end{aligned} \quad (5) $$

where $\mathcal{P}(\kappa_i, \kappa_j; \theta)$ is the joint PDF of two points of the smoothed convergence field at distance $\theta$. This equation can be simplified to

$$ \langle \mathcal{P}(\kappa_i) \mathcal{P}(\kappa_j) \rangle = \int d\theta P_d(\theta) \mathcal{P}(\kappa_i, \kappa_j; \theta), \quad (6a) $$

with $P_d(\theta)$ indicating the distribution of angular distances in a given survey area. The covariance of the PDF estimator can then be computed from this moment via

$$ \operatorname{cov}(\mathcal{P}(\kappa_i), \mathcal{P}(\kappa_j)) = \langle \hat{P}(\kappa_i) \hat{P}(\kappa_j) \rangle - \hat{P}(\kappa_i) \hat{P}(\kappa_j), \quad (6b) $$

where $\hat{P} \equiv \langle \hat{P} \rangle \approx \mathcal{P}$. We see that covariance predictions rely on two key ingredients: the distance distribution of angular separations and the two-point PDF of the lensing convergence at those separations, which will be discussed in more detail in Section 2.3.

**Distance distributions.** The relevant distance distributions for weak lensing describe a two-dimensional field. For the full sky, the uniform distribution of angular separations (measured in radians) reads

$$ P_{d, \text{full sky}}(\theta) = \frac{1}{2} \sin \theta, \quad (7) $$

with a maximum at $\theta_{max} = \pi/2$. This distribution when written in terms of distances relative to the maximal distance $\Theta_{side} = \pi$ is illustrated as red line in Figure 3.

For a square survey patch, the distribution of angular distances relative to the size of the survey patch, $\hat{\theta} = \theta/\Theta_{side}$, is

$$ P_{d, \text{square}}(\hat{\theta}) = \begin{cases} 2\hat{\theta}[\pi + (\hat{\theta} - 4)\hat{\theta}] & \hat{\theta} \in [0, 1[, \\ 2\hat{\theta}[\pi - 2 + 4\sqrt{\hat{\theta}^2 - 1} - \hat{\theta}^2 - 4\sec^{-1}(\hat{\theta})] & \hat{\theta} \in [1, \sqrt{2}]. \end{cases} \quad (8) $$

This distribution is illustrated as black line in Figure 3, with a maximum at roughly half the patch size where $\hat{\theta}_{max} \approx 0.479$ and a linear growth for small distances, $P_{d,square}(\hat{\theta}) \approx 2\pi\hat{\theta}$.

Before looking at PDF covariance matrix examples in Section 2.2 and introducing models for the two-point PDF in Section 2.3, let us discuss the main universal effects affecting the PDF covariance: i) the shot-noise effect related to the finite sampling of cells and ii) the super-sample covariance effect related to a finite survey area.

### 2.1.1. Finite sampling (shot noise) effects

Previously we discussed the covariance of the PDF in the continuum limit. In the more realistic discrete case, we estimate the PDF in a bin centred at $\kappa_i$ with width $\Delta_i$ (in our case linearly spaced bins $\Delta_i = \Delta$), so $\kappa \in [\kappa_i - \Delta_i/2, \kappa_i + \Delta_i/2]$ using a finite sample of $N_T$ cells that are placed on a regular grid in the survey area. The histogram-value of the PDF in a bin is obtained from averaging the continuous PDF $\mathcal{P}(\kappa)$ over that bin

$$ \mathcal{P}_i = \int_{-\Delta_i/2}^{\Delta_i/2} \frac{\mathcal{P}(\kappa_i + \delta\kappa_i)}{\Delta_i} \approx \mathcal{P}(\kappa_i), \quad (9) $$

and for small bin sizes close to the value of the continuous PDF at the bin center $\bar{\kappa}_i$. To compute the covariance for the discrete case, we can discretise equation (6) by averaging the continuous joint PDF over the $\kappa$-bins and replacing the integration over the distance distribution by a sum-based average over the set of $N_d = N_T^2 = \sum_{\theta \in \Theta} N_d(\theta)$ distances $\Theta$ of all grid cells with multiplicities $N_d(\theta)$ that are part of the survey patch

$$ \begin{aligned} \langle \mathcal{P}_i \mathcal{P}_j \rangle &= \sum_{\theta \in \Theta} \frac{N_d(\theta)}{N_d} \int_{-\Delta_i/2}^{\Delta_i/2} \int_{-\Delta_j/2}^{\Delta_j/2} \frac{\mathcal{P}(\kappa_i + \delta\kappa_i, \kappa_j + \delta\kappa_j; \theta)}{\Delta_i \Delta_j} \\ &= \sum_{\theta \in \Theta} \frac{N_d(\theta)}{N_d} \mathcal{P}_{ij}(\theta), \end{aligned} \quad (10) $$

where for small enough bins $\mathcal{P}_{ij}(\theta) \approx \mathcal{P}(\kappa_i, \kappa_j; \theta)$. In Figure 3 we show a comparison between the continuous distribution of distances in a square patch and the discrete case for a coarser grid of non-overlapping cells. We now split this sum over cell separations showing that the contribution from zero separation leads to a finite sampling or shot-noise term, and the remainder resembles the super-sample covariance contribution that dominates for the case of well-separated cells. Since the one-point PDF can be obtained from the two-point PDF by marginalisation, $\mathcal{P}(\kappa_i) = \int d\kappa_j \mathcal{P}(\kappa_i, \kappa_j; \theta)$, for zero separation we require that $\mathcal{P}(\kappa_i, \kappa_j; \theta = 0) = \mathcal{P}(\kappa_i)\delta_D(\kappa_i - \kappa_j)$. In the discrete case the joint two-point PDF at zero separation
is given by $\mathcal{P}_{ij}(\theta = 0) = \bar{\mathcal{P}}_i \delta_{ij} / \Delta_i$, required by the relation $\bar{\mathcal{P}}_i = \sum_j \Delta_j \mathcal{P}_{ij}(\theta = 0)$. The zero separation multiplicity equals the cell number $N_d(\theta = 0) = N_T$, such that the zero-separation term in equation (10) produces a finite sampling (or shot-noise) covariance contribution<fnref n="2" />

$$ \operatorname{cov}_{\text{FS}}(\mathcal{P}_i, \mathcal{P}_j) = \frac{N_T}{N_d} \mathcal{P}_{ij}(\theta = 0) = \frac{\bar{\mathcal{P}}_i}{\Delta_i N_T} \delta_{ij}. \quad (11) $$

Unsurprisingly, this finite sampling term is inversely proportional to the number of cells. If the cells are placed at a fixed separation, an increase of the area corresponds to a suppression of the shot-noise covariance term as $\sigma_{\text{FS}}^2(\mathcal{P}) \propto N_T^{-1} \propto \Theta_{\text{side}}^{-2} \propto A_{\text{survey}}^{-1}$, which is inversely proportional to the area.

While in principle the finite sampling term can be rendered arbitrarily small by using an increasing number of cells, this comes at the price of heavy cell overlaps at small separations. In this regime, a small-scale expansion with small $\Delta_{\xi}^{1/2}(\theta) = \sqrt{\sigma^2 - \xi(\theta)} \ll \xi(\theta) \ll 1$ can be used to capture effects of heavy cell overlaps as discussed in Bernardeau (2022). In the closed-form two-point PDF models discussed in Section 2.3 those overlap effects will be automatically included.

### 2.1.2. Super-sample covariance effects

Beyond assuming specific (log-)normal models for the bivariate PDF, one can obtain accurate estimates of the two-point PDF in the large separation regime illustrated in Figure 4. A large-scale expansion with small correlation $\xi(\theta) \ll \Delta_{\xi}(\theta) \ll 1$ is suitable for well-separated cells where the joint PDF at leading order can be approximated as Codis et al. (2016b)

$$ \mathcal{P}(\kappa_i, \kappa_j; \theta \gtrsim 2\theta_s) \approx \mathcal{P}(\kappa_i)\mathcal{P}(\kappa_j)[1 + \xi(\theta)b_1(\kappa_i)b_1(\kappa_j)], \quad (12) $$

with some bias function $b_1$, which describes how the two-point correlation at a given separation is modulated by the local weak-lensing convergences. This functional has been found to be very robust and the previous result extends what is expected in cosmological density fields (Bernardeau 1996). When cells are at sufficiently large separations and the initial conditions are Gaussian, this modulation is independent of separation. We will introduce a discretised estimator in equation (16c) and derive the functional form of those bias functions for specific forms of the two-point PDF discussed in Section 2.3 by an expansion in the correlation $\xi(\theta)$. For the related case of densities in spheres and cylinders, this bias function can be accurately predicted by spherical or cylindrical collapse, respectively. The large-deviation statistics techniques used to predict the one-point weak lensing PDF (Barthelemy et al. 2020a; Boyle et al. 2021) based on cylindrical collapse for thin redshift slices could be

extended to derive the functional form of this bias, although this is beyond the scope of this work. This large-separation expansion predicts a leading order super-sample covariance term and further sub-leading contributions that we will describe in Section 3. When all cells are at sufficiently wide separations $\theta \gtrsim 2\theta_s$, the covariance is obtained from integrating the large-separation limit of the two-point PDF (12) to obtain a super-sample covariance (SSC) covariance term

$$ \operatorname{cov}_{\text{SSC}}(\mathcal{P}(\kappa_i), \mathcal{P}(\kappa_j)) = \bar{\xi}(b_1 \mathcal{P})(\kappa_i)(b_1 \mathcal{P})(\kappa_j), \quad (13a) $$

where we defined the mean correlation

$$ \bar{\xi} = \int^{\theta_{\text{max}}} d\theta P_d(\theta)\xi(\theta), \quad (13b) $$

with $P_d(\theta)$ being the probability distribution of angular distances and $\theta_{\text{max}}$ being the maximum distance set by the size of the survey patch, respectively. This mean correlation agrees with the variance of the mean of $\kappa$ measured across different patches, so $\bar{\xi} = \sigma_{\kappa}^2$.

For a square patch of side length $\Theta_{\text{side}}$, the maximal distance in the patch is $\theta_{\text{max}} \simeq \sqrt{2}\Theta_{\text{side}}$ where the normalised distance distribution $P_d(\hat{\theta})$ is given in Equation (8), by conservation of probability we have

$$ P_d(\theta) = P_d(\hat{\theta}) \frac{d\hat{\theta}}{d\theta} = P_d\left(\hat{\theta} = \frac{\theta}{\Theta_{\text{side}}}\right) \frac{1}{\Theta_{\text{side}}} \stackrel{\hat{\theta} \ll 1}{\propto} \Theta_{\text{side}}^{-2}. \quad (14) $$

The main contribution in the mean correlation integral comes from small angles $\theta$, where the correlation function is largest and $P_d(\hat{\theta}) \propto \hat{\theta}$. Hence, the SSC covariance term from equation (13) scales as $\operatorname{cov}_{\text{SSC}}(\mathcal{P}_i, \mathcal{P}_j) \propto \Theta_{\text{side}}^{-2} \propto A_{\text{survey}}^{-1}$, which is inversely proportional to the survey area.

In the discrete case, the covariance can be expressed as a sum of the previously discussed finite sampling term (11) and a super-sample covariance like term

$$ \operatorname{cov}(\mathcal{P}_i, \mathcal{P}_j) = \langle \mathcal{P}_i \mathcal{P}_j \rangle - \bar{\mathcal{P}}_i \bar{\mathcal{P}}_j \quad (15a) $$

$$ = \operatorname{cov}_{\text{FS}}(\mathcal{P}_i, \mathcal{P}_j) \quad (15b) $$

$$ + \frac{1}{N_d} \sum_{\theta \in \Theta \setminus \{0\}} \left[ \mathcal{P}_{ij}(\theta) - \bar{\mathcal{P}}_i \bar{\mathcal{P}}_j \underbrace{\frac{N_d}{N_T(N_T - 1)}}_{\approx 1} \right]. $$

If cells at separation $\theta$ were completely independent, such that $\mathcal{P}_{ij}(\theta) = \bar{\mathcal{P}}_i \bar{\mathcal{P}}_j$, then only the shot-noise term from the finite sampling of cells would remain. The second term in equation (15) includes significant contributions from cells at wide separation $\theta \gtrsim 2\theta_s$, where we can use the large-separation PDF expansion (12). This leads to a super-sample covariance term analogous to equation (13)

$$ \operatorname{cov}_{\text{SSC}}(\mathcal{P}_i, \mathcal{P}_j) = \bar{\xi}(b_1 \mathcal{P})_i (b_1 \mathcal{P})_j, \quad (16a) $$

where $(b_1 \mathcal{P})_i$ are obtained from a bin-average of the involved functions resembling equation (9) and the mean correlation can be obtained as

$$ \bar{\xi} = \sum_{\theta \in \Theta \setminus \{0\}} \frac{N_d(\theta)}{N_d} \xi(\theta). \quad (16b) $$

Using the relationship $\langle \kappa'_j | \kappa_i \rangle_{\kappa'_j} = b_1(\kappa_i) / \bar{\xi}$, the value

<footnote n="2">² This shot-noise term was already described in (Codis et al. 2016b). Even if there was no spatial correlation in the sample, in the Poisson limit the covariance of finding $N_i$ cells in convergence bin $\operatorname{cov}_{\text{SN}}(N_i, N_j) = \bar{N}_i \delta_{ij}$, where $\bar{N}_i$ is the mean number of cells expected to be in the convergence bin centred on $\bar{\kappa}_i$. The number of cells, $N_i$ in a given convergence bin of width $\Delta_i$ centered on $\kappa_i$ is related to the histogram-value of the binned PDF from equation (9) such that $N_i = N_T \Delta_i \mathcal{P}_i$ and similarly $\bar{N}_i = N_T \Delta_i \bar{\mathcal{P}}_i$.</footnote>
<figure><img src="image_5.png" /><figcaption>FIG. 4.— Two-point correlation of circular cells of radius $\theta_s \simeq 7.3'$ (solid black) in comparison to the matter correlation (dashed black), the variance in cells (dashed blue) and the difference $\Delta_{\xi}(\theta)$ (red line) that vanishes in the limit of small separations. The blue dotted lines indicate the point at which the cell separation equals twice the smoothing scale, while the green dotted lines indicate the point where the correlation is half of the variance, corresponding to slightly smaller scales.</figcaption></figure>

of the bias function in the bin $\kappa_i$ can be estimated by comparing the mean of the convergence at separation $\theta = 2\theta_s$, $\kappa'_j$, from this given convergence $\kappa_i$ to the average correlation $\bar{\xi}$ using a sum over the 6 nearest neighbour (NN) pixels at that separation

$$ \hat{b}_1(\kappa_i) = \frac{\sum_{\mathrm{P}} \sum_{\mathrm{P}' \in \mathrm{NN}(2\theta_s)} \chi_i(\kappa_{\mathrm{P}}) \kappa_{\mathrm{P}'} }{6 \sum_{\mathrm{P}} \chi_i(\kappa_{\mathrm{P}})} . \quad (16c) $$

This definition is the analogue of the estimator for densities in spheres and cylinders (Codis et al. 2016b; Uhlemann et al. 2017, 2018).

Since the correlation function decreases with increasing angular separation and the probability distribution of cell distances peaks at roughly half the patch size, this super-sample covariance term is the most relevant for small patches while it is heavily suppressed for full sky maps. In some cases, like for simulation-based analysis or estimation of covariances, it is desirable to construct small patches which closely resemble the full sky behaviour. The main contribution from the SSC covariance term can be reduced by subtracting the mean $\kappa$ in every small patch, thus reducing the mean correlation $\bar{\xi}$ and modifying the shape of the bias function $b$ as we will discuss in Section 3.2.

## 2.2. Illustrative covariance measurements

To validate our covariance models, we generate a large set of Gaussian as well as shifted lognormal (FLASK) $\kappa$ maps (with the shift parameter tuned to fit the $\kappa$-PDF skewness). We start from 10 full-sky convergence maps (healpix $N_{\mathrm{side}} = 4096$ corresponding to 0.85' resolution) that we each divide into 1200 non-overlapping square patches with side-length of 5deg. We chose the patch size to resemble the typical size of simulated convergence maps for which differing cosmologies are available (Liu et al. 2018; Giocoli et al. 2018; Harnois-Déraps et al. 2018). We build a data vector $\vec{S}$ from the values of the

<figure><img src="image_6.png" /><figcaption>FIG. 5.— Diagonal of the covariance matrix for the non-Gaussian FLASK PDFs (solid lines) and the Gaussian version (dashed lines) measured from small patches with non-overlapping cells (red), overlapping cells (blue) and overlapping cells with mean-subtraction in every patch (purple). The Gaussian case shows a symmetric covariance contribution as expected. The corresponding correlation matrices are shown in Figure 6. Additionally, we show our theoretical predictions for the Gaussian case (dotted lines, partially obscured by the dashed lines). The dotted vertical lines indicate the range of $\kappa$ values considered for the covariance matrices and Fisher forecasts.</figcaption></figure>

PDF histogram in the central region $-0.024 < \kappa < 0.026$ measured from the top-hat smoothed $\kappa$ map by considering the regular grid of all cells (overlapping) or only cells separated by at least twice the smoothing radius (non-overlapping). The covariance matrix of the binned $\kappa$-PDF across those patches is obtained using the formula

$$ C_{ij} = \langle (S_i - \bar{S}_i)(S_j - \bar{S}_j) \rangle, \quad \bar{S}_i = \langle S_i \rangle, \quad (17) $$

where $\langle \cdot \rangle$ indicates an ensemble average over different map realisations for the full sky or patches, respectively. We rescale all covariances to mimic a total Euclid-like survey area $A_{\mathrm{survey}} = 15,000 \mathrm{ deg}^2$ using the rescaling factor $A_{\mathrm{patch}}/A_{\mathrm{survey}}$ for the patches and $A_{\mathrm{fullsky}}/A_{\mathrm{survey}}$ for the full-sky maps. We verified the validity of this rescaling using covariance measurements from full sky maps in addition to the theoretical arguments in Section 2.1.2. For visualisation purposes, it is often useful to normalise the covariance matrix, with components $C_{ij}$, by its diagonal components, which defines the correlation matrix $R$ such that

$$ R_{ij} = \frac{C_{ij}}{\sqrt{C_{ii}C_{jj}}}. \quad (18) $$

In Figures 5 and 6 we illustrate two effects that impact the overall size and structure of the PDF covariance and the derived correlation matrix, finite sampling (shot noise) and super-sample covariance. In Figure 5 we show the diagonal of the covariance matrix for FLASK patches (solid lines) and Gaussian patches (dashed lines). In Figure 6 we show the corresponding correlation matrices ordered by decreasing covariance, with the upper triangle showing the FLASK patches and the lower triangle showing the Gaussian patches.

We consider three cases for the covariance in small patches of side length 5deg:
<figure><img src="image_7.png" /><figcaption>FIG. 6.— Correlation matrix between the bins of the weak lensing convergence $\kappa$ PDF measured at a smoothing scale $\approx 7.3'$ measured from $10 \times 1200$ (5 deg)$^2$ patches of 10 full sky FLASK maps of $\kappa$ including the nonzero mean in the patches (left and middle panel) with non-overlapping (NO) cells (left) and overlapping cells (middle) and $\kappa - \bar{\kappa}$ after subtracting the mean with overlapping cells (right panel). The upper triangles labelled by FLASK show results from shifted-lognormal realisations, while the lower triangles those of Gaussian realisations. The diagonal of the covariance matrix is shown in Figure 5.</figcaption></figure>

* raw $\kappa$ measured from non-overlapping cells (red lines in Figure 5, left panel in Figure 6),

* raw $\kappa$ measured from overlapping cells (blue lines in Figure 5, middle panel in Figure 6), and

* mean-subtracted $\bar{\kappa} = \kappa - \bar{\kappa}$ from overlapping cells (purple lines in Figure 5, right panel in Figure 6).

We see that while the lognormal maps have a visible asymmetry in the covariance diagonal compared to the Gaussian maps, both lead to very similar correlation matrices. Comparing the red and blue lines shows that given a fixed total area, increasing the number of cells for the PDF measurement decreases the finite sampling (or shot noise term) that determines the covariance matrix diagonal, as we explained in Section 2.1.1. The suppression of the diagonal for the case of an increased number of cells comes at the price of making cells overlap which creates a strong positive correlation band around the diagonal, as can be seen from comparing the left and middle panel in Figure 6. Comparing the blue and purple lines in Figure 5 illustrates that subtracting the mean in small patches significantly reduces the impact of the super-sample covariance term discussed in Section 2.1.2 and causes the 4-tiled pattern in the middle panel of Figure 6 to transform to a 9-tiled pattern in the right panel. Additionally, the removal of this dominant term causes the overlap contribution in the band along the diagonal to become more apparent. As we will see later, the mean-subtracted convergence in small patches provides a proxy for the convergence on the full sky, as evident from comparing the right panel of Figure 6 to Figures 1 and 9. While most of our discussion is focused on the underlying weak lensing convergence field, we discuss the inclusion of shape noise arising from cosmic shear measurements in Section 2.3.4.

## 2.3. Specific models for the joint two-point PDF

As evident from equation (6), the PDF covariance can be computed directly if the two-point PDF $P(\kappa_i, \kappa_j; \theta)$ is known a priori for all separations. In the following we will discuss three possible cases relying on closed-form expressions for that joint PDF, a bivariate Gaussian (normal) PDF, a bivariate shifted lognormal PDF and a bivariate PDF built from a marginal one-point PDF and a correlation structure, for example specified by an expansion for large separations or heavy overlaps.

### 2.3.1. Bivariate Gaussian

A Gaussian joint PDF is parametrised by the one-point PDF variance $\sigma^2$ and the two-point correlation between cells $\xi = \xi_{12}(\theta)$ such that

$$P_G(\kappa_1, \kappa_2; \theta) = P_G(\kappa_1, \kappa_2; \sigma^2, \xi_{12}(\theta)) \quad (19)$$

$$= \frac{\exp\left[-\frac{\sigma^2}{2d}[\kappa_1^2 + \kappa_2^2] + \frac{\xi_{12}(\theta)}{d}\kappa_1\kappa_2\right]}{2\pi\sqrt{d}}, \quad (20)$$

with the determinant $d = \sigma^4 - \xi_{12}^2 = \sigma^4[1 - (\xi_{12}/\sigma^2)^2]$. One can perform a large-separation expansion for $\xi_{12}/\sigma^2 \ll 1$ valid in the regime of no overlaps as shown in Figure 4. At leading order, one can approximate $d \approx \sigma^4$ and expand the remaining $\xi_{12}(\theta)$ term in the exponential to obtain

$$\frac{P_G(\kappa_1, \kappa_2; \theta)}{P_G(\kappa_1)P_G(\kappa_2)} = 1 + \xi_{12}(\theta) \frac{\kappa_1 \kappa_2}{\sigma^2 \sigma^2} + \mathcal{O}(\xi^2), \quad (21)$$

which is a meaningful expansion in the regime $\kappa \lesssim \mathcal{O}(\sigma)$ and reproduces the linear ‘Kaiser’ bias function (Kaiser 1984; Codis et al. 2016a; Uhlemann et al. 2017)

$$b_{1,G}(\kappa) = \kappa/\sigma^2. \quad (22)$$

We will discuss a large-separation expansion up to higher orders in $\xi(\theta)/\sigma^2$ in Section 3.1.2.
**The impact of mean-subtraction.** For small patches, the mean convergence $\kappa_s$ will significantly fluctuate from patch to patch with a sizeable average correlation $\bar{\xi} = \sigma_{\kappa_s}^2$ causing the super-sample covariance. To better emulate measurements of modern surveys on a significant fraction of the full sky, subtracting this mean convergence in each patch in the spirit of the peak-background split (Bardeen et al. 1986; Mo & White 1996; Sheth & Tor men 1999) can be desirable. The two-point PDF of the mean-subtracted lensing convergence $\tilde{\kappa}_i = \kappa_i - \bar{\kappa}$ is ob tained from a trivariate Gaussian PDF for $\kappa_1, \kappa_2$ and the background value in the patch $\kappa_s = \bar{\kappa}$, which is then integrated out

$$ \mathcal{P}_G(\tilde{\kappa}_1, \tilde{\kappa}_2) = \int d\kappa_s \mathcal{P}_G(\kappa_1 = \tilde{\kappa}_1 + \kappa_s, \kappa_2 = \tilde{\kappa}_2 + \kappa_s, \kappa_s). \qquad (23a) $$

The covariance matrix is of the form

$$ \Sigma_{\kappa_1, \kappa_2, \kappa_s} = \begin{pmatrix} \sigma^2 & \xi_{12} & \xi_{1s} \\ \xi_{12} & \sigma^2 & \xi_{2s} \\ \xi_{1s} & \xi_{2s} & \sigma_s^2 \end{pmatrix}, \qquad (23b) $$

where we assume $\bar{\xi} = \sigma_s^2 = \xi_{1s} = \xi_{2s}$, which holds for large enough patches compared to the cell size. The integration yields a bivariate Gaussian PDF for the mean-subtracted convergence $\mathcal{P}(\tilde{\kappa}_1, \tilde{\kappa}_2)$ with reduced variance $\sigma_{\tilde{\kappa}}^2 = \sigma^2 - \bar{\xi}$ and reduced correlation $\xi_{\tilde{\kappa}}(\theta) = \xi(\theta) - \bar{\xi}$, where the reduction is set by the variance at the patch scale, $\bar{\xi} = \sigma_s^2$. The mean-subtracted correlation $\xi - \bar{\xi}$ now leads to a vanishing leading order super-sample covariance term from equation (24) and creates a new effective leading order term that we will discuss in Section 3.2.1. Obviously this above described formalism has a long history which goes back to early papers on constrained Gaussian random fields e.g (Bertschinger 1987) and (Hoffman & Ribak 1991)

**Covariance.** When obtaining the covariance according to equation (6) from the large-separation expansion of the Gaussian (21), the leading 1 disappears and the correlation $\xi_{12}(\theta)$ is integrated over the distance distribution to yield $\bar{\xi}$ which sets the amplitude of a super-sample covariance like term we will discuss further in Section 2.1.2

$$ \operatorname{cov}_G(\mathcal{P}(\kappa_i), \mathcal{P}(\kappa_j)) = \bar{\xi}(b_1 \mathcal{P})_G(\kappa_i)(b_1 \mathcal{P})_G(\kappa_j) + \mathcal{O}(\bar{\xi}^2). \qquad (24) $$

In Figure 7 we compare the correlation matrices as measured (upper triangle) and predicted (lower triangle) for Gaussian fields of the three cases discussed before (generated with FLASK, with the same initial power spectrum and smoothing as before), from non-overlapping cells with raw $\kappa$ (left) along with overlapping cells of raw $\kappa$ (middle) and mean-subtracted $\tilde{\kappa} = \kappa - \bar{\kappa}$ (right panel). In the non-overlapping case on the left-hand side, we can see how the finite sampling term (11) causes the diagonal to be most pronounced with a sub-dominant 4-tile pattern showing the leading order super-sample covariance like term from equation (24). When making cells overlap as shown in the middle panel, the finite sampling term becomes irrelevant and the 4-tile pattern [due to the factorisation of that covariance into a product of two independent bias functions with one zero-crossing each] as seen in equation (24)] becomes more pronounced, while the resulting cell overlap causes a significant increase of the correlation in nearby bins as can

be seen in a band along the diagonal. Considering the mean-subtracted convergence shown in the right panel will remove the leading order, super-sample covariance term from equation (24) and instead lead to a 9-tile pattern [due to the factorisation of that covariance contribution into a products of two bias functions with two zero-crossings each, see equation (56)].

### 2.3.2. Bivariate shifted lognormal

A shifted lognormal joint PDF for the lensing convergence $\kappa$ is parametrised by the variance $\sigma_{\kappa}^2$ and two-point correlation $\xi_{\kappa} = \xi_{12,\kappa}(\theta)$ along with a shift parameter $s$ related to the skewness, which we will keep fixed throughout here. We can define a Gaussian-distributed variable $g(\kappa)$ with zero-mean

$$ g(\kappa) = \ln\left(\frac{\kappa}{s} + 1\right) - \mu_G = \ln\left(\frac{\kappa}{s} + 1\right) + \frac{\sigma_G^2}{2}, \qquad (25) $$

where $\mu_G = -\sigma_G^2/2$ is fixed to ensure $\langle \kappa \rangle = 0$ and the Gaussian variance is $\sigma_G^2 = \ln(1 + \sigma_{\kappa}^2/s^2)$. This can be used to construct the shifted lognormal univariate PDF as

$$ \mathcal{P}_{\text{LN}}(\kappa; \sigma_G^2, s) = \frac{\Theta(\kappa + s)}{\sqrt{2\pi}\sigma_G(\kappa + s)} \exp\left[\frac{-g(\kappa)^2}{2\sigma_G^2}\right], \qquad (26) $$

where $\Theta$ indicates the Heaviside step function, such that the PDF is zero if $\kappa \le -s$. The shift parameter, $s$, can be set to replicate the desired skewness (see the formula for $\lambda$ in Xavier et al. 2016)

$$ \begin{gathered} s = \frac{\sigma_{\kappa}}{\tilde{\mu}_3} (1 + y(\tilde{\mu}_3)^{-1} + y(\tilde{\mu}_3)) - \langle \kappa \rangle \approx 3 \frac{\sigma_{\kappa}^4}{\langle \kappa^3 \rangle}, \qquad (27) \\ y(\tilde{\mu}_3) = \sqrt[3]{\frac{2 + \tilde{\mu}_3^2 + \tilde{\mu}_3 \sqrt{4 + \tilde{\mu}_3^2}}{2}} \approx 1, \end{gathered} $$

where $\langle \kappa \rangle$ is the desired mean, $\sigma$ the target variance and $\tilde{\mu}_3 = \langle \kappa^3 \rangle / \sigma_{\kappa}^3$ the skewness. The shift parameter is typically positive and for our fiducial weak lensing convergence case we use $s = 0.1145$ throughout the main text as done in Boyle et al. (2021). In the executive summary, we use a different shift parameter because the Takahashi et al. (2017) simulations have a different cosmology. The Gaussian limit can be reproduced by sending $s/\sigma_{\kappa} \to \infty$.

The bivariate shifted lognormal distribution also depends on the correlation $\xi_G = \ln(1 + \xi_{\kappa}/s^2)$ and the determinant $d = \sigma_G^4 - \xi_G^2$

$$ \begin{gathered} \mathcal{P}_{\text{LN}}(\kappa_1, \kappa_2; \sigma_G^2, \xi_G, s) = \frac{\Theta(\kappa_1 + s)\Theta(\kappa_2 + s)}{2\pi\sqrt{d}(\kappa_1 + s)(\kappa_2 + s)} \qquad (28) \\ \times \exp\left[-\frac{\sigma_G^2}{2d}[g(\kappa_1)^2 + g(\kappa_2)^2] + \frac{\xi_G}{d}g(\kappa_1)g(\kappa_2)\right]. \end{gathered} $$

The joint shifted lognormal PDF can be expanded at large separations where $\xi_G(\xi_{\kappa}) / \sigma_G^2 \ll 1$ to obtain up to leading order in $\xi_{\kappa}$

$$ \frac{\mathcal{P}_{\text{LN}}(\kappa_1, \kappa_2; \theta)}{\mathcal{P}_{\text{LN}}(\kappa_1)\mathcal{P}_{\text{LN}}(\kappa_2)} = 1 + \frac{\xi_{\kappa}(\theta)}{s^2} \frac{g(\kappa_1)}{\sigma_G^2} \frac{g(\kappa_2)}{\sigma_G^2} + \mathcal{O}(\xi_{\kappa}^2). \qquad (29) $$
<figure><img src="image_8.png" /><figcaption>FIG. 7.— Correlation matrix between the bins of the Gaussian weak lensing convergence $\kappa$ PDF at a smoothing scale $\approx 7.3'$ measured from 10×1200 (5 deg)$^2$ patches from 10 full sky Gaussian maps (upper triangle) and the theory prediction (lower triangle) for $\kappa$ including the nonzero mean in the patches (left and middle panel) with non-overlapping (NO) cells (left) and overlapping cells (middle) and $\kappa - \bar{\kappa}$ after subtracting the mean with overlapping cells (right panel). The predictions for the diagonal of the covariance matrices are shown in Figure 5.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>FIG. 8.— $\kappa$-dependent clustering as measured by the cell bias from one full-sky map in the Takahashi simulation in comparison to the theoretical prediction for a shifted lognormal field from equation (30).</figcaption></figure>

In particular, it is useful to define a leading order bias function resembling the Gaussian expression $b_{1,G} = \kappa/\sigma^2$

$$b_{1,\text{LN}}(\kappa) = \frac{g(\kappa)}{s\sigma_G^2} \approx \frac{sg(\kappa)}{\sigma_\kappa^2}, \quad (30)$$

where for $\kappa \ll s$ we have $sg(\kappa) \approx \kappa + \sigma_\kappa^2/(2s)$ and hence reproduce the Gaussian result modulo a small shift. For larger $\kappa$ the shifted lognormal result deviates from the Gaussian one to closely resemble measurements from simulated $\kappa$ maps (Takahashi et al. 2017) using the estimator (16c) as illustrated in Figure 8.

### The impact of mean-subtraction.

When describing small patches which will feature different background ‘sample’ convergence values $\kappa_s$, one needs to integrate over a trivariate PDF of the convergences $\kappa_{1/2}$ and the sample value $\kappa_s$, similarly as for the

mean-subtraction in the Gaussian case (23a)

$$\mathcal{P}(\tilde{\kappa}_1, \tilde{\kappa}_2) = \int d\kappa_s \mathcal{P}_{\text{LN}}(\kappa_1 = \tilde{\kappa}_1 + \kappa_s, \kappa_2 = \tilde{\kappa}_2 + \kappa_s, \kappa_s). \quad (31)$$

To evaluate this, let us assume that $\tilde{\kappa}$ is a shifted log-normal random variable (Hilbert et al. 2011), i.e. given in terms of a zero-mean Gaussian random variable $\tilde{g} = \ln(\tilde{\kappa}/\tilde{s} + 1) + \tilde{\sigma}_G^2/2$ with variance $\tilde{\sigma}_G^2 = \ln(1 + \sigma_{\tilde{\kappa}}^2/\tilde{s}^2)$ and a new lognormal shift $\tilde{s}$. The correlation of the Gaussian variables is given by $\tilde{\xi}_G = \ln(1 + \xi_{\tilde{\kappa}}/\tilde{s}^2)$. To find an effective shifted lognormal approximation for $\tilde{\kappa}$, we must determine variance, correlation and the shift parameter $\tilde{s}$. Recall that we know the variance of $\kappa \equiv \kappa_s + \tilde{\kappa}$ from any procedure for calculating the matter power spectrum. In the same way, we can also know the covariance $\langle \tilde{\kappa}_i \tilde{\kappa}_s \rangle = \xi_{is}$ and the variance $\langle \tilde{\kappa}_s^2 \rangle = \sigma_s^2$, which we will approximate as identical $\xi_{is} = \sigma_s^2 = \tilde{\xi}$. So we get the following variance and correlation for the mean-subtracted $\tilde{\kappa}$

$$\sigma_{\tilde{\kappa}}^2 = \langle \tilde{\kappa}^2 \rangle = \langle \kappa^2 \rangle - 2\langle \kappa \kappa_s \rangle + \langle \kappa_s^2 \rangle \approx \sigma_{\kappa}^2 - \tilde{\xi}, \quad (32a)$$

$$\begin{aligned} \xi_{\tilde{\kappa}} &= \langle \tilde{\kappa}_1 \tilde{\kappa}_2 \rangle = \langle \kappa_1 \kappa_2 \rangle - \langle \kappa_1 \kappa_s \rangle - \langle \kappa_2 \kappa_s \rangle + \langle \kappa_s^2 \rangle \\ &\approx \xi_{\kappa} - \tilde{\xi}, \end{aligned} \quad (32b)$$

and hence the Gaussian variance $\tilde{\sigma}_G^2 = \ln[1 + (\sigma_{\tilde{\kappa}}^2 - \tilde{\xi})/\tilde{s}^2]$ and correlation $\tilde{\xi}_G = \ln[1 + (\xi_{\tilde{\kappa}} - \tilde{\xi})/\tilde{s}^2]$. The shift parameter can be related to the skewness as described in equation B.11 of Hilbert et al. (2011)

$$\langle \tilde{\kappa}^3 \rangle = \frac{3}{\tilde{s}} \langle \tilde{\kappa}^2 \rangle^2 + \frac{1}{\tilde{s}^3} \langle \tilde{\kappa}^2 \rangle^3 \approx \frac{3}{\tilde{s}} \langle \tilde{\kappa}^2 \rangle^2 \Rightarrow \tilde{s} \approx \frac{3 \langle \tilde{\kappa}^2 \rangle^2}{\langle \tilde{\kappa}^3 \rangle}. \quad (33)$$

Following Friedrich et al. (2018) we use the shift parameter $\tilde{s}$ to match the skewness of $\tilde{\kappa}$, which is obtained from the third moment of the difference

$$\langle \tilde{\kappa}^3 \rangle = \langle (\kappa - \kappa_s)^3 \rangle = \langle \kappa^3 \rangle - 3\langle \kappa^2 \kappa_s \rangle + 3\langle \kappa \kappa_s^2 \rangle. \quad (34)$$
<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>FIG. 9.— (Left) Correlation matrix between the bins of the shifted lognormal weak lensing convergence $\kappa$ PDF at a smoothing scale $\approx 7.3'$ comparing FLASK measurements (upper triangles) to the theory predictions (lower triangles). (left panel) Measurements from 20×1200 square patches of (5deg)$^2$ size cut from 20 full sky realisations in comparison to the lognormal theory prediction in patches. While the first order bias function $b_{1,LN}$ (30) creates the 4-tile pattern, cell overlaps create the positive correlation band along the diagonal. (right panel) Measurements from 500 full sky realisations in comparison to the lognormal theory prediction on the full sky. The 9-tile pattern is set by the second order bias function $b_{2,LN}$ (58), while cell overlaps create the strong correlation along the diagonal. (Right) Diagonal of the covariance matrix for the non-Gaussian FLASK PDFs measured from the full sky (green solid) and as predicted by the lognormal theory (green dotted). We also show the result for the FLASK patches without/with mean subtraction (purple/blue solid) and the lognormal theory computed for the patches without/with mean-subtraction (purple/blue dotted). The shaded bands are indicating estimated errors)</figcaption></figure>

The third moment of $\kappa$ can be expressed in terms of the variance $\sigma_{\kappa}$ and original shift parameter $s$ using $\langle \kappa^3 \rangle \approx 3\sigma_{\kappa}^4/s$. Additionally, we can compute the relevant mixed moments of a zero-mean normal random variable $Y$ and a lognormal variable $Z = s_z(e^X - 1)$ given in terms of a normal variable $X$. The relevant relations are $\langle Y(e^X - 1) \rangle = \xi_{XY}$ and $\langle Y^2(e^X - 1) \rangle = \xi_{XY}^2 (\sigma_X^2 + 1 - \sigma_X^{-2})$ which yields $\langle Y(e^X - 1)^2 \rangle = 0$. This can be used to obtain $\langle \kappa^2 \kappa_s \rangle = 0$ and $\langle \kappa \kappa_s^2 \rangle \approx \bar{\xi}^2 (\sigma_G^2 + 1 - \sigma_G^{-2})/s$. Inserting this into equations (33) and (34) leads to a relation between the shift parameters $s$ and $\tilde{s}$, which we found to be practically identical in our case.

**Covariance.** Again, we can integrate the large-separation expansion of the shifted lognormal (57a) over the distance distribution to obtain the covariance from equation (6) and eventually get at leading order

$$ \text{cov}_{\text{LN}}(\mathcal{P}(\kappa_i), \mathcal{P}(\kappa_j)) = \bar{\xi}(b_1 \mathcal{P})_{\text{LN}}(\kappa_i)(b_1 \mathcal{P})_{\text{LN}}(\kappa_j) + \mathcal{O}(\bar{\xi}^2), \quad (35) $$

which is in form analogous to the Gaussian result (24) despite the modified ingredients. This explains why the correlation matrix from small patches (with sizable $\bar{\xi}$) in the FLASK case closely resembles the Gaussian case, as seen in Figure 6 while the diagonal is modified due to the non-Gaussian shape of the final convergence PDF as seen in Figure 5.

In Figure 9 we compare the measurements from FLASK maps to the theoretical predictions for shifted lognormal fields. In the left panel on the left side we show the correlation matrix obtained from 1200 FLASK patches of 5deg each cut from 20 full sky realisations (upper triangle) along with the theoretical prediction (lower triangle) finding excellent agreement. We can see that the correlation matrix has a 4-tile structure. This is because the mean correlation $\bar{\xi}$ in the patches is sizeable such that the first order covariance term (35) dominates the visual appearance. In the right panel on the left side we show the correlation matrix obtained from 500 full sky

FLASK maps (upper triangle) along with the theoretical prediction (lower triangle) finding excellent agreement. We can see that the correlation matrix closely resembles the case of the mean-subtracted convergence measured in small patches, shown in the lower panel of Figure 6. This is because the mean correlation $\bar{\xi}$ on a full sky is so small that the first order covariance term (35) becomes subdominant compared to the second order term that we will discuss in Section 3.1.4. The latter term also appears as leading order contribution in the mean-subtracted case addressed in Section 3.2.3. In the right panel we focus on the diagonal of the covariance matrix, comparing FLASK maps and theoretical predictions for different cases. This plot shows that while the PDF covariance of the raw $\kappa$ in small patches looks significantly different, the mean-subtracted $\bar{\kappa}$ in small patches resembles the full sky covariance. The theoretical prediction captures all three behaviours quite well, considering relative 1-$\sigma$ errors on the diagonal of the covariance matrix (shaded bands) following from the covariance between the numerically estimated diagonal covariance elements $\hat{C}_{ii}$ and $\hat{C}_{jj}$ as given in Taylor et al. (2013)

$$ \text{Cov}(\hat{C}_{ii}, \hat{C}_{jj}) = \frac{2C_{ij}^2}{N_{\text{sim}} - 1}. \quad (36) $$

We notice that while the lognormal prediction for the full sky and the raw $\kappa$ in patches agree extremely well with the predictions, there are some minor discrepancies for the mean-subtracted $\bar{\kappa}$ in small patches. As this discrepancy was absent for Gaussian fields shown in Figure 5, we suspect that it is caused by a more complicated patch-to-patch correlation induced by cutting multiple patches from one full sky, which is unaccounted for in the theory describing the mean-subtracted convergence as another shifted lognormal random field.

### 2.3.3. Bivariate Gaussian copula with general marginals

In case neither the bivariate normal nor the shifted lognormal distributions can provide a good description
of the target PDF as its univariate marginal, one can design a bivariate PDF by specifying the marginal and correlation structure independently. A joint PDF can be built from a given marginal one-point PDF $P(\kappa)$ and a correlation structure encoded in a copula density $c_{12}$ which quantifies the correlation between ranked variables obtained from the cumulative distribution function $C(\kappa) = \int^\kappa d\tilde{\kappa} P(\tilde{\kappa})$, such that

$$P(\kappa_1, \kappa_2; \theta) = P(\kappa_1)P(\kappa_2)c_{12}(C(\kappa_1), C(\kappa_2); \theta), \quad (37)$$

where $\theta$ indicates angular separation as before. If the $\kappa_i$ and $\kappa_2$ are independent, then the copula density is unity $c_{12,\text{indep}} = 1$. A Gaussian copula only depends on the distance through the pairwise correlation $\xi_{12}(\theta)$, $c_{12,G}(C(\kappa_1), C(\kappa_2); \xi_{12}(\theta))$ and is a known function which can be computed numerically although it has no analytical expression. Beyond those two cases, systematically constructing the dependence structure is possible in the large-separation regime, where the copula can be expanded in powers of the pairwise correlation $\xi_{12}(\theta)$ as we will further discuss in Section 3.

### 2.3.4. Inclusion of shape noise

Observationally, the weak lensing convergence map is obtained from cosmic shear measurements, which are subject to shape noise often constituting the dominant source of noise. Since galaxies are intrinsically elliptical, the observed shear contains a contribution from this intrinsic signal, quantified by the variance of the intrinsic ellipticity $\sigma_\epsilon$. In simulated maps, shape noise can be included by adding a white noise – modeled by a Gaussian random map $\kappa_{SN}$ – to the ‘raw’ simulated convergence map $\kappa$ before smoothing. The standard deviation of this Gaussian is typically assumed to be

$$\sigma_{SN} = \frac{\sigma_{\epsilon}}{\sqrt{n_g \cdot \Omega_{\theta_s}}} \propto \theta_s^{-1}, \quad (38)$$

with the shape-noise parameter $\sigma_\epsilon$ (here taken to be $\sigma_\epsilon = 0.30$), the galaxy density $n_g$ (here taken to be $30 \text{ arcmin}^{-2}$) and the solid angle $\Omega_{\theta_s}$ (here taken to correspond to circular apertures of radius $\theta_s = 7.3'$).

Smoothing induced shape noise correlations and covariance. At the raw $\kappa$ map level, the shape noise values in different pixels are uncorrelated and $\xi \equiv 0$. Top-hat smoothing the field at angular size $\theta_s$ will correlate shape noise values in overlapping cells (being at separation $\theta < 2\theta_s$). Their correlation is related to the area of the cell overlap (which is a symmetric lens)

$$A(\theta < 2\theta_s) = 2\theta_s^2 \arccos\left(\frac{\theta}{2\theta_s}\right) - \frac{\theta}{2}\sqrt{4\theta_s^2 - \theta^2}. \quad (39a)$$

The resulting correlation of the top-hat smoothed shot noise contribution to the weak lensing convergence is

$$\xi_{\text{SN}}(\theta < 2\theta_s) = \sigma_{\text{SN}}^2 A(\theta) / (\pi\theta_s^2), \quad \xi_{\text{SN}}(\theta \ge 2\theta_s) = 0. \quad (39b)$$

To obtain the covariance of the shot noise contribution alone, we integrate the bivariate Gaussian with variance $\sigma_{\text{SN}}^2$ and covariance (39) over distances following equa-

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /><figcaption>FIG. 10.— PDF covariance diagonal (upper panel) and correlation matrix (lower panel) for the pure shape-noise lensing convergence in circular apertures of radius $\theta_s$ comparing the prediction from equation (40) to the measurement from shape noise realisations. The slight apparent asymmetry between positive and negative $\kappa$ is caused by the noise due to the finite number of realisations.</figcaption></figure>

tion (6) to get

$$\langle \mathcal{P}_{\text{SN}}(\kappa_i) \mathcal{P}_{\text{SN}}(\kappa_j) \rangle = \int_0^{2\theta_s} d\theta P_d(\theta) \mathcal{P}_G(\kappa_i, \kappa_j; \sigma_{\text{SN}}^2, \xi_{\text{SN}}(\theta)) + \mathcal{P}_G(\kappa_i; \sigma_{\text{SN}}^2) \mathcal{P}_G(\kappa_j; \sigma_{\text{SN}}^2) \int_{2\theta_s}^{\infty} d\theta P_d(\theta),$$

where $P_d(\theta)$ indicates the distribution of angular distances in the survey area, which is approximately linear for small distances $P_d(\theta) \approx 2\pi\theta$. This gives the PDF covariance for the case of pure shape noise as

$$\text{cov}_{\text{SN}}(\mathcal{P}(\kappa_i), \mathcal{P}(\kappa_j)) = \int_0^{2\theta_s} d\theta P_d(\theta) \mathcal{P}_G(\kappa_i, \kappa_j; \sigma_{\text{SN}}^2, \xi_{\text{SN}}(\theta)) - \mathcal{P}_G(\kappa_i; \sigma_{\text{SN}}^2) \mathcal{P}_G(\kappa_j; \sigma_{\text{SN}}^2) \int_0^{2\theta_s} d\theta P_d(\theta). \quad (40)$$
In Figure 10 we compare this theoretical result to a measurement across 1000 shape noise realisations for square maps of length 5deg with a top-hat smoothing of 7.3' finding excellent agreement.

### Shape noise induced convergence PDF convolution.

Shape noise impacts the one-point weak lensing convergence PDF as if it was convolved with a zero-centred Gaussian (Clerkin et al. 2016)

$$ 
\begin{aligned}
\mathcal{P}_{\kappa+\text{SN}}(\hat{\kappa}) &= (\mathcal{P}_{\text{G}}(\sigma_{\text{SN}}) * \mathcal{P})(\hat{\kappa}) \\
&= \frac{1}{\sqrt{2\pi}\sigma_{\text{SN}}} \int d\kappa \exp\left(-\frac{(\hat{\kappa} - \kappa)^2}{2\sigma_{\text{SN}}^2}\right) \mathcal{P}(\kappa),
\end{aligned}
\quad (41) $$

where the shape noise standard deviation $\sigma_{SN}$ is given in equation (38). For the two-point PDF of the smoothed weak lensing we similarly have $\hat{\kappa}_i = \kappa_i + \kappa_{i,SN}$ and the bivariate PDF is given in terms of a 2D-convolution

$$ \mathcal{P}_{\kappa+\text{SN}}(\hat{\kappa}_1, \hat{\kappa}_2; \theta) = (\mathcal{P}_{\text{G}}(\sigma_{\text{SN}}, \xi_{\text{SN}}(\theta)) * \mathcal{P})(\hat{\kappa}_1, \hat{\kappa}_2; \theta), \quad (42) $$

where the shape noise contribution is described by a bivariate Gaussian PDF from equation (19) with standard $\sigma_{SN}$ (38) and correlation $\xi_{SN}(\theta)$ (39). If the joint PDF of raw $\kappa_i$ was Gaussian, then the noisy PDF of $\hat{\kappa}_i$ arising from the convolution would be another zero-mean Gaussian distribution with covariance matrix $\Sigma_{\kappa+\text{SN}} = \Sigma_{\kappa} + \Sigma_{\text{SN}}$.

For the case of a shifted lognormal $\kappa$ PDF, we recompute the covariance from the convolved two-point PDF (42) to find that the main impact of adding shape noise is to change the diagonal of the covariance matrix, shown in Figure 11, while leaving the correlation matrix largely similar. This result was obtained from $N_{\text{sim}} = 500$ FLASK maps of the full sky with a smoothing of 7.3' applied to the raw $\kappa$ and the shape-noise added $\hat{\kappa}$, respectively. Considering relative $1-\sigma$ errors on the diagonal of the covariance matrix (36) evaluated for $N_{\text{sim}}$ the predictions are in good agreement with the measurements.

## 3. LARGE SEPARATION COVARIANCE EXPANSIONS

In this Section we will perform an expansion of the joint two-point PDF in powers of the two-point cell correlation function $\xi_{12}(\theta)$, partially with the help of the cumulant generating function (CGF). We will pay particular attention to the leading and next-to-leading order contributions to the covariance matrix, which reveals a super-sample covariance term that is controlled by the average correlation in the survey patch. We will compute the expansion to all orders in $\xi_{12}$ for a Gaussian field and the two leading order terms for the general non-Gaussian case, illustrated by results for normal and shifted lognormal fields. This will naturally lead to an expansion of the covariance matrix that resembles an eigendecomposition, which we will discuss in Section 4.1.

### 3.1. Covariance expansion for raw $\kappa$

#### 3.1.1. Leading order determined by $\kappa$-dependent clustering

Let us first recall that at leading order in the large-separation regime, the joint two-point PDF can be approximated in terms of the one-point PDF and the

<figure><img src="image_14.png" /><figcaption>FIG. 11.— PDF covariance diagonal for the noise-free lensing convergence as measured from 500 FLASK full sky maps (green dashed, shaded band indicating estimated errors) and predicted for the shifted lognormal (green dotted) in comparison to the noisy lensing convergence as measured from FLASK (green solid, shaded band indicating estimated errors) and predicted from equation (6) with the noisy two-point PDF (42) (green dot-dashed).</figcaption></figure>

‘bias’ functions indicating the density-dependent two-point clustering

$$ \mathcal{P}(\kappa_1, \kappa_2, \theta) \approx \mathcal{P}(\kappa_1)\mathcal{P}(\kappa_2)[1 + \xi(\theta)b(\kappa_1)b(\kappa_2)]. \quad (43) $$

Hence, we have for the covariance of the binned probability

$$ \operatorname{cov}(\mathcal{P}(\kappa_i), \mathcal{P}(\kappa_j)) = \bar{\xi}(\bar{\mathcal{P}}b)(\kappa_i)\bar{\mathcal{P}}b(\kappa_j) + \delta_{ij} \frac{\bar{\mathcal{P}}(\kappa_i)}{\Delta_i N_T}, \quad (44) $$

where $N_T$ is the total number of cells and the mean correlation $\xi$ is given by equation (13b). Note that the shot-noise term only acts on the diagonal. For the diagonal, the leading-order prediction for covariance diagonal is then

$$ \operatorname{var}(\mathcal{P}(\kappa)) = \bar{\xi}(b\bar{\mathcal{P}})^2(\kappa) + \frac{\bar{\mathcal{P}}(\kappa)}{\Delta N_T}, \quad (45) $$

where $\Delta$ is the PDF bin width. For a Gaussian field, the bias function is given by linear Kaiser bias $b(\kappa) = \kappa/\sigma_{\kappa}^2$, which seems to be a decent approximation for the central region of $\kappa$-values, see Figure 8. Note that the cell bias $b(\kappa)$ has at least one zero crossing, for the case of raw $\kappa$, around $\kappa \simeq 0$, while for the mean-subtracted case this is modified as described in Section 3.2.

Let us also recall that for non-overlapping cells, typically the ‘finite sampling’ effect dominates (well-described by Poisson). For overlapping cells, the number of cells is extremely large thus effectively removing the finite sampling term, but requiring the treatment of overlaps. The impact of overlapping cells can be captured by either closed-form expressions for the two-point PDF valid at small distances (like for the Gaussian and shifted-lognormal case discussed here) or a heavy overlap expansion complementary to the large separation one as described in Bernardeau (2022).

#### 3.1.2. Beyond leading order: Gaussian case

For a Gaussian field of zero mean, we have the following simple structure of the joint two-point cumulant
generating function in terms of the two-point correlation $\xi = \xi_{12}(\theta)$ as follows

$$ \varphi_G(\lambda_1, \lambda_2) = \varphi_{0,G}(\lambda_1) + \varphi_{0,G}(\lambda_2) + \xi\varphi_{1,G}(\lambda_1)\varphi_{1,G}(\lambda_2), \quad (46) $$

where $\varphi_{0,G} = \lambda^2\sigma^2/2$, $\varphi_{1,G}(\lambda) = \lambda$ and $\varphi_{n\ge 2,G}(\lambda) = 0$. By expanding the last term $\exp[\xi\lambda_1\lambda_2]$ in a power series and converting $\lambda_i$ to derivatives w.r.t. $\kappa_i$, we obtain

$$ \frac{\mathcal{P}_G(\kappa_1, \kappa_2, \theta)}{\mathcal{P}_G(\kappa_1)\mathcal{P}_G(\kappa_2)} = 1 + \sum_{n=1}^{\infty} \frac{\xi^n(\theta)}{n!} b_{n,G}(\kappa_1)b_{n,G}(\kappa_2) \quad (47) $$

$$ b_{n,G}(\kappa) = \frac{(-1)^n}{\mathcal{P}_G(\kappa)} \frac{\partial^n \mathcal{P}_G(\kappa)}{\partial \kappa^n} = \frac{1}{\sigma^n} He_n\left(\frac{\kappa}{\sigma}\right), \quad (48) $$

which can be written in terms of the probabilist's Hermitian polynomials $He_n$

$$ He_n(x) = (-1)^n \exp\left[\frac{x^2}{2}\right] \frac{d}{dx^n} \exp\left[-\frac{x^2}{2}\right], \quad (49) $$

when using $x = \kappa/\sigma$. The first two terms read (Desjacques et al. 2018, see also)

$$ b_{1,G}(\kappa) = \frac{\kappa}{\sigma^2}, \quad b_{2,G}(\kappa) = \frac{\kappa^2 - \sigma^2}{\sigma^4} = b_{1,G}^2 - \frac{1}{\sigma^2}, \quad (50) $$

which can be also obtained from a direct expansion of the joint PDF (19). This leads to the following expansion for the covariance following equation (6)

$$ \operatorname{cov}(\mathcal{P}_G(\kappa_1), \mathcal{P}_G(\kappa_2)) = \sum_{n=1}^{\infty} \frac{\bar{\xi}^n}{n!} \frac{\partial^n \mathcal{P}_G(\kappa_1)}{\partial \kappa_1^n} \frac{\partial^n \mathcal{P}_G(\kappa_2)}{\partial \kappa_2^n} \quad (51) $$

$$ = \sum_{n=1}^{\infty} \frac{\bar{\xi}^n}{n!\sigma^{2n}} He_n\left(\frac{\kappa_1}{\sigma}\right) \mathcal{P}_G(\kappa_1) He_n\left(\frac{\kappa_2}{\sigma}\right) \mathcal{P}_G(\kappa_2), $$

$$ \bar{\xi}^n := \int d\theta P_d(\theta) \xi^n(\theta) \approx \sum_{\theta \in \Theta \setminus \{0\}} \frac{N_d(\theta)}{N_d} \xi^n(\theta), \quad (52) $$

where in the second line we have defined the average of powers of the two-point correlation function over the survey patch in analogy to the mean from equation (13b). The averages of powers of the correlation function $\bar{\xi}^n$ compared to powers of the variance $\sigma^{2n}$ control the large-separation expansion of the covariance and their hierarchy is illustrated in Figure 12. This plot also shows that the covariance contributions scale inverse proportional to the survey area $\bar{\xi}^n \propto A^{-1}$.

Due to the hierarchical ordering of averages of the correlation function we expect the two leading order terms in the regime $\kappa \lesssim \mathcal{O}(\sigma)$ to be

$$ \frac{\operatorname{cov}(\mathcal{P}_G(\kappa_1), \mathcal{P}_G(\kappa_2))}{\mathcal{P}_G(\kappa_1)\mathcal{P}_G(\kappa_2)} = \frac{\bar{\xi} \frac{\kappa_1}{\sigma^2} \frac{\kappa_2}{\sigma^2}}{\sigma^2} \quad (53a) $$

$$ + \frac{\bar{\xi}^2}{2\sigma^4} \frac{\kappa_1^2 - \sigma^2}{\sigma^2} \frac{\kappa_2^2 - \sigma^2}{\sigma^2} \quad (53b) $$

$$ + \mathcal{O}\left(\frac{\bar{\xi}^3}{\sigma^6}\right), $$

averages 2-point cell correlation vs. variance powers

<figure><img src="image_15.png" /><figcaption>FIG. 12.— Scaling of the average of powers of the two-point cell correlation function compared to the variance for the raw $\kappa$ maps in 5deg patches (red plusses) and the full sky result rescaled by a factor of $A_{\text{patch}}/A_{\text{fullsky}}$ (blue plusses). Aside from the removal of the first order term, those hierarchies closely resemble each other as well as the mean-subtracted $\bar{\kappa}$ maps in patches (red crosses).</figcaption></figure>

which predicts a super-sample covariance leading order term similar to the separate universe picture (81), being proportional to the variance of the mean $\kappa$ in different patches, $\xi = \sigma^2(\kappa_s)$, and the derivative of the PDF. This large-separation expansion of the covariance will be connected to an eigendecomposition in Section 4.1. Here, the hierarchical ordering of $\bar{\xi}^n/\sigma^{2n}$ will reflect in the eigenvalues and the bias functions multiplied by the Gaussian PDF evaluated at a given $\kappa$-value will yield the eigenvectors (after orthogonalisation and appropriate normalisation). In this context, the predicted zero-crossings of the first bias function $b_{1,G}(\kappa)$ at $\kappa = 0$ and of the second $b_{2,G}(\kappa)$ at $\kappa = \sigma^2$ will translate into the approximate zero crossings of the first two eigenvectors as we will show later in Figure 15.

### 3.1.3. Beyond leading order: general non-Gaussian case

Higher-order terms are obtained from expanding the joint two-point cumulant generating function in powers of the two-point correlation $\xi = \xi_{12}(\theta)$ as follows

$$ \varphi(\lambda_1, \lambda_2) = \varphi_0(\lambda_1) + \varphi_0(\lambda_2) + \xi\varphi_1(\lambda_1)\varphi_1(\lambda_2) + \frac{\xi^2}{2}[\varphi_1^2(\lambda_1)\varphi_2(\lambda_2) + \varphi_2(\lambda_1)\varphi_1^2(\lambda_2)] + \mathcal{O}(\xi^3) \quad (54) $$

We are then in position to compute the (connected part) of the joint PDF,

$$ \mathcal{P}(\kappa_1, \kappa_2) - \mathcal{P}(\kappa_1)\mathcal{P}(\kappa_2) = \int \frac{d\lambda_1}{2\pi i} \int \frac{d\lambda_2}{2\pi i} \exp[-\lambda_1\kappa_1 - \lambda_2\kappa_2 + \varphi(\lambda_1, -\lambda_1) + \varphi(\lambda_2, -\lambda_2)] \times \quad (55a) $$

$$ + \frac{\xi^2}{2}[\varphi_1^2(\lambda_1)\varphi_2(\lambda_2) + \varphi_2(\lambda_1)\varphi_1^2(\lambda_2)] + \mathcal{O}(\xi^3) \quad (55b) $$

We obtain the covariance according to equation (6) after integrating over $\lambda_i$ and the distance distribution and
symmetrising

$$ \frac{\mathrm{cov}(\mathcal{P}(\kappa_1), \mathcal{P}(\kappa_2))}{\mathcal{P}(\kappa_1)\mathcal{P}(\kappa_2)} = \xi b_1(\kappa_1)b_1(\kappa_2) \quad (56a) $$

$$ + \frac{\xi^2}{2} [(b_2 + q_1)(\kappa_1)(b_2 + q_1)(\kappa_2) - q_1(\kappa_1)q_1(\kappa_2)] \quad (56b) \\ + \mathcal{O}(\xi^3), $$

where we defined

$$ (b_1 \mathcal{P})(\kappa) = \int \frac{d\lambda}{2\pi i} \varphi_1(\lambda) \exp[-\lambda\kappa + \varphi(\lambda)] \quad (56c) $$

$$ (b_2 \mathcal{P})(\kappa) = \int \frac{d\lambda}{2\pi i} \varphi_1(\lambda)^2 \exp[-\lambda\kappa + \varphi(\lambda)] \quad (56d) $$

$$ (q_1 \mathcal{P})(\kappa) = \int \frac{d\lambda}{2\pi i} \varphi_2(\lambda) \exp[-\lambda\kappa + \varphi(\lambda)] \quad (56e) $$

For a saddle-point approximation of the integral, we would typically expect $b_2 \approx b_1^2$, although that does not reproduce the Gaussian or shifted lognormal result. The expression (56) consists of dyadic products that resemble the structure of an eigendecomposition, which we will look at in Section 4.1, although eigenvectors are additionally pairwise orthogonal.

### 3.1.4. Beyond leading order: shifted lognormal case

The recipe described in the previous paragraph for general non-Gaussian distributions is not valid for a lognormal PDF, which is not well-described in terms of a CGF. Here we detail an approach that can be used to obtain lognormal bias functions by applying a large-separation expansion to the shifted lognormal two-point PDF (28). This bias expansion is similar in spirit to the Gaussian case, but requires an additional expansion of the Gaussian field correlation $\xi_G$ in terms of the underlying $\kappa$ correlation $\xi_\kappa$. The joint shifted lognormal PDF can be expanded at large separations where $\xi_G(\xi_\kappa)/\sigma_G^2 \ll 1$ to obtain up to next-to-leading order in $\xi_\kappa$

$$ \begin{aligned} \frac{\mathcal{P}_{\mathrm{LN}}(\kappa_1, \kappa_2; \xi_\kappa)}{\mathcal{P}_{\mathrm{LN}}(\kappa_1)\mathcal{P}_{\mathrm{LN}}(\kappa_2)} &= 1 + \xi_G(\xi_\kappa) \frac{g(\kappa_1)}{\sigma_G^2} \frac{g(\kappa_2)}{\sigma_G^2} \quad (57a) \\ &\quad + \frac{\xi_G^2(\xi_\kappa)}{2} \frac{(g(\kappa_1)^2 - \sigma_G^2)}{\sigma_G^4} \frac{(g(\kappa_2)^2 - \sigma_G^2)}{\sigma_G^4} \\ &\quad + \mathcal{O}(\xi_G^3) \\ &= 1 + \left(\frac{\xi_\kappa}{s^2} - \frac{\xi_\kappa^2}{2s^4}\right) \frac{g(\kappa_1)}{\sigma_G^2} \frac{g(\kappa_2)}{\sigma_G^2} \quad (57b) \\ &\quad + \frac{\xi_\kappa^2}{2s^4} \frac{(g(\kappa_1)^2 - \sigma_G^2)}{\sigma_G^4} \frac{(g(\kappa_2)^2 - \sigma_G^2)}{\sigma_G^4} \\ &\quad + \mathcal{O}(\xi_\kappa^3), \end{aligned} $$

In particular, it is useful to define two additional bias functions $q_{1,LN}$ and $b_{2,LN}$ following the previous expansion (56), with their combination resembling the Gaussian expression from equation (50)

$$ (b_2 + q_1)_{\mathrm{LN}}(\kappa) = \frac{g(\kappa)^2 - \sigma_G^2}{s^2 \sigma_G^4} \approx \frac{[sg(\kappa)]^2 - \sigma_\kappa^2}{\sigma_\kappa^4} \quad (58a) $$

$$ q_{1,LN}(\kappa) = \frac{1}{s} b_{1,LN}(\kappa) = \frac{g(\kappa)}{s^2 \sigma_G^2} \approx \frac{g(\kappa)}{\sigma_\kappa^2}, \quad (58b) $$

<figure><img src="image_16.png" /><figcaption>FIG. 13.— A comparison of a rescaled version of the second order bias functions from the shifted lognormal model used here $b_{2,LN}$ (solid) along with $(b_2 + q_1)_{LN}$ (dot-dashed) from equations (58) (colours for different variances) and the limiting behaviour (59) (black solid) in comparison to the Gaussian result (black dashed). We also show lines of constant $\sigma^2$ (dotted) to indicate the zero crossings of $b_{2,LN}(\kappa)$ which appear in the second eigenvectors for the $\kappa$ PDF covariance in the lower panel of Figure 16.</figcaption></figure>

where we used the result for $b_{1,LN}$ from equation (30). We see that in the limit of $s \to \infty$ we recover $sg(\kappa) \to \kappa$ and the Gaussian result with $q_{1,G} = 0$, while for mildly non-Gaussian fields, we typically have $q_1 \ll b_2$. In a small variance expansion where $\sigma/s \ll 1$, we obtain the limiting behaviour

$$ b_{2,LN}(\kappa) \xrightarrow{\frac{\sigma}{s} \ll 1} \frac{s^2 \ln \left(1 + \frac{\kappa}{s}\right)^2 - \sigma^2}{\sigma^4} - \frac{1}{4s^2} \quad (59) $$

Using the all-order results for $b_{n,G}$ (48) together with the expansion of $\xi_G(\xi_\kappa) = \sum_{n=1} \frac{(-1)^{n-1}}{n} \left(\frac{\xi_\kappa}{s^2}\right)^n$ allows to compute all order expressions for the lognormal bias functions. Note that due to the difference between $\xi_G$ and $\xi_\kappa$, the higher-order terms receive contributions from lower-order bias functions, as predicted by the CGF-based expansion from before. In Figure 13 we show a comparison of the lognormal terms $b_2$ (solid) and $b_2 + q_1$ (dashed) to the Gaussian result for the shift parameter adopted in this work and different variances. We can see that for our case (closely corresponding to the green line), the $q_1$ contribution is negligible. The crossing between the solid $b_2$ line and the dotted line of constant $\sigma^2$ indicates the zero-crossing of $b_2(\kappa)$, which is reflected in the second eigenvector for the $\kappa$ covariance shown in the lower panel of Figure 16. For the case of a survey area being a significant portion of the sky, the mean correlation $\xi$ will be extremely small and this second order bias function $b_2(\kappa)$ will determine the covariance structure. The two zero-crossings of $b_2$ cause the 9-tiled correlation matrix structure displayed earlier, notably in the upper right panel of Figure 9.

## 3.2. Covariance expansion for mean-subtracted

$\bar{\kappa} = \kappa - \kappa_s$

To determine the PDF covariance for the case of a mean-subtracted $\kappa$ field, one needs to consider the 3 variables $\kappa_1$, $\kappa_2$ and $\kappa_s$, where $\kappa_{1/2}$ are as before and $\kappa_s$ is the mean convergence in the survey patch. Like the con-
vergence field, those variables all have zero expectation value and we are then interested in the joint PDF of the mean-subtracted $\tilde{\kappa}_i$

$$ \tilde{\kappa}_1 = \kappa_1 - \kappa_s, \quad \tilde{\kappa}_2 = \kappa_2 - \kappa_s. \qquad (60) $$

Let us note that this tilde variable is formally equivalent to the slope variable that was introduced in Bernardeau et al. (2014, 2015). The CGF for the three variables defined above is written here up to linear order in the two-point function at the survey scale, $\xi_s$, and up to quadratic order in the two-point correlation at the given cell separation $\xi_{12}$

$$ \begin{aligned} \varphi(\lambda_1, \lambda_2, \lambda_s) &= \frac{1}{2}\sigma_s^2\lambda_s^2 + \varphi_0(\lambda_1) + \varphi_0(\lambda_2) + \xi_{12}\varphi_1(\lambda_1)\varphi_1(\lambda_2) \\ &\quad + \lambda_s[\xi_{1s}\varphi_1(\lambda_1) + \xi_{2s}\varphi_1(\lambda_2)] \\ &\quad + \frac{\xi_{12}^2}{2}[\varphi_1^2(\lambda_1)\varphi_2(\lambda_2) + \varphi_2(\lambda_1)\varphi_1^2(\lambda_2)] + \mathcal{O}(\xi_{12}^3) \\ \varphi(\lambda_1, \lambda_s) &= \frac{1}{2}\xi_s\lambda_s^2 + \varphi_0(\lambda_1) + \lambda_s\xi_{1s}\varphi_1(\lambda_1), \end{aligned} \qquad (62) $$

where we assumed that the patches are big enough such that the $\kappa_s$ PDF is Gaussian. The one- and two-point PDFs are obtained from this CGF by performing the change of variable (60) and then marginalising over the background $\kappa_s$, leading to

$$ \begin{aligned} \mathcal{P}(\tilde{\kappa}_1) &= \int \frac{d\lambda_1}{2\pi i} e^{-\lambda_1 \tilde{\kappa}_1 + \varphi(\lambda_1, -\lambda_1)} \\ \mathcal{P}(\tilde{\kappa}_1, \tilde{\kappa}_2) &= \iint \frac{d\lambda_1}{2\pi i} \frac{d\lambda_2}{2\pi i} e^{-\lambda_1 \tilde{\kappa}_1 - \lambda_2 \tilde{\kappa}_2 + \varphi(\lambda_1, \lambda_2, -\lambda_1 - \lambda_2)}. \end{aligned} \qquad (64) $$

When all computations are made at linear order in $\xi$, we can set $\sigma_s^2 = \xi_{1s} = \xi_{2s} = \xi_{12} = \xi(\theta)$. To obtain the second order, one needs to start from $\bar{\xi} = \sigma_s^2 = \xi_{1s} = \xi_{2s} \neq \xi_{12} = \xi(\theta)$, which we will use to illustrate the similarities and differences between the Gaussian and non-Gaussian cases.

### 3.2.1. Gaussian case

For a Gaussian field, the two terms $\lambda_s^2$ and $\lambda_s[\varphi_1(\lambda_1) + \varphi_1(\lambda_2)]$ can be summarised and act to reduce the variance and correlation function by $\bar{\xi}$, such that joint PDF $\mathcal{P}(\tilde{\kappa}_1, \tilde{\kappa}_2)$ is a Gaussian with variance $\sigma_{\tilde{\kappa}}^2 = \sigma^2 - \bar{\xi}$ and correlation $\xi_{\tilde{\kappa}}(\theta) = \xi(\theta) - \bar{\xi}$. Hence, the covariance can be obtained in full analogy to the previous computation that led to equation (51) which is modified to read

$$ \operatorname{cov}(\mathcal{P}_G(\tilde{\kappa}_1), \mathcal{P}_G(\tilde{\kappa}_2)) = \sum_{n=2}^{\infty} \frac{(\bar{\xi} - \bar{\xi})^n}{n!} \frac{\partial^n \mathcal{P}_G(\tilde{\kappa}_1)}{\partial \tilde{\kappa}_1^n} \frac{\partial^n \mathcal{P}_G(\tilde{\kappa}_2)}{\partial \tilde{\kappa}_2^n}, \qquad (65) $$

where due to the definition of the average of powers of the correlation function (52) the $n=1$ term disappears. Note that $\mathcal{P}_G(\tilde{\kappa}_1)$ does not denote the same PDF as $\mathcal{P}_G(\kappa_1)$, for the Gaussian case they differ by the variance being $\sigma_{\tilde{\kappa}}^2 = \sigma_{\kappa}^2 - \bar{\xi}$ rather than $\sigma_{\kappa}^2$. Aside from this (small) difference and the removal of the $n=1$ term, we now have a different prefactor consisting of averages of powers of $\bar{\xi} - \bar{\xi}$ rather than $\bar{\xi}$. Those averages of powers of the correlation function $\bar{\xi}^n$ and $(\bar{\xi} - \bar{\xi})^n$ compared

to powers of the variances $\sigma^2$ and $(\sigma^2 - \bar{\xi})^n$ control the large-separation expansion of the covariance (in the Gaussian case given in equation (51) for raw $\kappa$ and in equation (65) for the mean-subtracted $\tilde{\kappa} = \kappa - \bar{\kappa}$) follow similar hierarchies as illustrated in Figure 12.

In particular, the leading order term will be

$$ \begin{aligned} \frac{\operatorname{cov}(\mathcal{P}_G(\tilde{\kappa}_1), \mathcal{P}_G(\tilde{\kappa}_2))}{\mathcal{P}_G(\tilde{\kappa}_1)\mathcal{P}_G(\tilde{\kappa}_2)} &= \frac{(\bar{\xi} - \bar{\xi})^2}{2\sigma_{\tilde{\kappa}}^4} \underbrace{\left(\frac{\tilde{\kappa}_1^2}{\sigma_{\tilde{\kappa}}^2} - 1\right)\left(\frac{\tilde{\kappa}_2^2}{\sigma_{\tilde{\kappa}}^2} - 1\right)}{\sigma_{\tilde{\kappa}}^4 b_{2,G}(\tilde{\kappa}_1)b_{2,G}(\tilde{\kappa}_2)} \\ &\quad + \mathcal{O}\left(\frac{(\bar{\xi} - \bar{\xi})^3}{\sigma_{\tilde{\kappa}}^6}\right), \end{aligned} \qquad (66) $$

which predicts a first eigenvalue/-vector similar to the second eigenvalue/-vector of the Gaussian case without mean-subtraction (53b) because $(\bar{\xi} - \bar{\xi})^2/\sigma_{\tilde{\kappa}}^4 \simeq \bar{\xi}^2/\sigma^4$ as shown in Figure 12. The zero-crossings of the first eigenvector will be at $\tilde{\kappa} \approx \sigma^2 - \bar{\xi} \approx \sigma^2$, closely following the result for full sky maps where $\bar{\xi}$ is negligibly small and eigenvectors are shown as green/olive lines in the middle panel of Figure 16, as discussed in Section 4.1. Note that as before the joint two-point PDF of the mean-subtracted $\tilde{\kappa}$ values could be written in terms of the marginal Gaussian PDFs and a series of Hermite polynomials as in equation (48) with prefactors set by averages of powers of $(\bar{\xi} - \bar{\xi})/\sigma_{\tilde{\kappa}}^2$.

### 3.2.2. Leading order: general non-Gaussian case

Similarly as for the Gaussian case, we will rewrite the joint CGF of the three variables $\lambda_1, \lambda_2, \lambda_s$ in terms of the joint CGF for just two variables $\lambda_i$ and $\lambda_s$, where $\lambda_s = -(\lambda_1 + \lambda_2)$. We will first discuss the leading order result which reads

$$ \begin{aligned} & \varphi(\lambda_1, \lambda_2, -\lambda_1 - \lambda_2) - \varphi(\lambda_1, -\lambda_1) - \varphi(\lambda_2, -\lambda_2) \qquad (67) \\ &= \xi\varphi_1(\lambda_1)\varphi_1(\lambda_2) + \bar{\xi}[\lambda_1\lambda_2 - \lambda_1\varphi_1(\lambda_2) - \lambda_2\varphi_1(\lambda_1)] + \mathcal{O}(\xi^2) \\ &= (\bar{\xi} - \bar{\xi})\varphi_1(\lambda_1)\varphi_1(\lambda_2) + \bar{\xi}[\varphi_1(\lambda_1) - \lambda_1][\varphi_1(\lambda_2) - \lambda_2] + \mathcal{O}(\xi^2), \end{aligned} $$

where we used $\xi = \xi(\theta)$ for brevity and reorganised the formula such that the for a Gaussian only the first term remains, because $\varphi_1(\lambda) = \lambda$. We are then in position to compute the connected part of the joint PDF, given by

$$ \begin{aligned} \mathcal{P}(\tilde{\kappa}_1, \tilde{\kappa}_2) - \mathcal{P}(\tilde{\kappa}_1)\mathcal{P}(\tilde{\kappa}_2) &= \iint \frac{d\lambda_1}{2\pi i} \frac{d\lambda_2}{2\pi i} e^{-\lambda_1 \tilde{\kappa}_1 - \lambda_2 \tilde{\kappa}_2 + \varphi(\lambda_1, -\lambda_1) + \varphi(\lambda_2, -\lambda_2)} \\ &\quad \times \left[ (\bar{\xi} - \bar{\xi})\varphi_1(\lambda_1)\varphi_1(\lambda_2) + \bar{\xi}[\varphi_1(\lambda_1) - \lambda_1][\varphi_1(\lambda_2) - \lambda_2] + \mathcal{O}(\bar{\xi}^2, \bar{\xi}^2, (\bar{\xi} - \bar{\xi})^2) \right]. \end{aligned} $$

When integrated over the distance distribution, the term linear in $\bar{\xi} - \bar{\xi}$ will vanish and the term involving products of $\varphi_1(\lambda_i) - \lambda_i$ will create a new bias function $\tilde{b}_{1,NG}(\tilde{\kappa})$ for the mean-subtracted convergence $\tilde{\kappa}$ only present for non-Gaussian fields. In terms of this, the PDF covariance can be written as

$$ \operatorname{cov}(\mathcal{P}(\tilde{\kappa}_1), \mathcal{P}(\tilde{\kappa}_2)) = \bar{\xi}(\tilde{b}_{1,NG}\mathcal{P})(\tilde{\kappa}_1)(\tilde{b}_{1,NG}\mathcal{P})(\tilde{\kappa}_2) \quad (68a) \\ + \mathcal{O}(\{\bar{\xi}^2, \bar{\xi}^2, (\bar{\xi} - \bar{\xi})^2\}), $$
<figure><img src="image_17.png" /><figcaption>FIG. 14.— Additional first-order bias $\tilde{b}_{1,NG}(\tilde{\kappa})$ (green) of the mean-subtracted $\tilde{\kappa} = \kappa - \tilde{\kappa}$ obtained from the sum of the ‘raw’ bias $b_1$ (green) and the response of the PDF $\partial \log \tilde{P}/\partial \kappa$ (red) as given in equation (68b). We also show predictions from the shifted lognormal case (dashed) relying on the PDF from equation (26) and the bias $b_{1,LN}$ from equation (30)</figcaption></figure>

where we introduced the leading order composite bias

$$ \tilde{b}_{1,NG}(\tilde{\kappa}) = b_1(\tilde{\kappa}) + \frac{\partial}{\partial \tilde{\kappa}} \log \mathcal{P}(\tilde{\kappa}). \quad (68b) $$

In Figure 14 we illustrate the first composite non-Gaussian bias function (green) along with its ingredients, the first order bias $b_1$ (blue) and the logarithmic derivative of the PDF (red) as measured in one FLASK realisation (data points) and predicted for a shifted log-normal (dashed lines). We can see that for a mildly non-Gaussian field, this composite non-Gaussian bias function is more than an order of magnitude smaller than the raw bias $b_1$, such that there is no guarantee that our result at leading order in $\tilde{\xi}$ constitutes the most relevant correction. This can be established analytically for a shifted lognormal distribution using equations (26) and (30), which gives

$$ \tilde{b}_{1,LN}(\tilde{\kappa}) = -\frac{1}{\tilde{\kappa} + s} + \frac{x(\tilde{\kappa})}{s\sigma_G^2(\sigma_{\tilde{\kappa}}^2)} \left(1 - \frac{1}{1 + \tilde{\kappa}/s}\right). \quad (68c) $$

The minimum value is $\tilde{b}_{1,LN}(\tilde{\kappa} \approx 0) \approx -1/s$ and in the limit of $s \to \infty$ we recover the Gaussian result with vanishing $\tilde{b}_1$.

### 3.2.3. Next-to-leading order: general non-Gaussian case

We now generalise the previous result up to next-to-leading order, where we find

$$ \begin{aligned} & \varphi(\lambda_1, \lambda_2, -\lambda_1 - \lambda_2) - \varphi(\lambda_1, -\lambda_1) - \varphi(\lambda_2, -\lambda_2) \quad (69) \\ &= (\xi - \bar{\xi})\varphi_1(\lambda_1)\varphi_1(\lambda_2) + \bar{\xi}[\varphi_1(\lambda_1) - \lambda_1][\varphi_1(\lambda_2) - \lambda_2] \\ & \quad + \frac{\xi^2}{2}[\varphi_1^2(\lambda_1)\varphi_2(\lambda_2) + \varphi_2(\lambda_1)\varphi_1^2(\lambda_2)] + \mathcal{O}(\xi^3). \end{aligned} $$

Again, we compute the connected part of the joint PDF, given by

$$ \begin{aligned} & \mathcal{P}(\tilde{\kappa}_1, \tilde{\kappa}_2) - \mathcal{P}(\tilde{\kappa}_1)\mathcal{P}(\tilde{\kappa}_2) \\ &= \iint \frac{d\lambda_1}{2\pi i} \frac{d\lambda_2}{2\pi i} e^{-\lambda_1 \tilde{\kappa}_1 - \lambda_2 \tilde{\kappa}_2 + \varphi(\lambda_1, -\lambda_1) + \varphi(\lambda_2, -\lambda_2)} \\ & \quad \times \left[ (\xi - \bar{\xi})\varphi_1(\lambda_1)\varphi_1(\lambda_2) + \bar{\xi}[\varphi_1(\lambda_1) - \lambda_1][\varphi_1(\lambda_2) - \lambda_2] \right. \\ & \quad \left. + \frac{\bar{\xi}^2}{2}[\varphi_1(\lambda_1) - \lambda_1]^2[\varphi_1(\lambda_2) - \lambda_2]^2 \right. \\ & \quad \left. + \bar{\xi}(\xi - \bar{\xi})\varphi_1(\lambda_1)\varphi_1(\lambda_2)[\varphi_1(\lambda_1) - \lambda_1][\varphi_1(\lambda_2) - \lambda_2] \right. \\ & \quad \left. + \frac{\xi^2}{2}[\varphi_1^2(\lambda_1)\varphi_2(\lambda_2) + \varphi_2(\lambda_1)\varphi_1^2(\lambda_2)] \right. \\ & \quad \left. + \frac{(\xi - \bar{\xi})^2}{2}\varphi_1^2(\lambda_1)\varphi_1^2(\lambda_2) + \mathcal{O}(\xi^3) \right]. \end{aligned} \quad (70a) $$

When integrated over the distance distribution, all terms linear in $\tilde{\xi} - \bar{\xi}$ will vanish. The remaining contributions are

a) terms in the first two lines, involving powers of $\varphi_1(\lambda) - \lambda$ creating new bias functions, $\tilde{b}_{n,NG}$, only present for non-Gaussian fields.<fnref n="3" />

b) terms in the last two lines containing products of powers of $\varphi_n$ that yield the bias functions encountered for the covariance of the raw $\kappa$ PDF, at second order $b_2$ from equation (56d) and $q_1$ from equation (56e).

This yields

$$ \begin{aligned} \operatorname{cov}(\mathcal{P}(\tilde{\kappa}_1), \mathcal{P}(\tilde{\kappa}_2)) ={}& \bar{\xi}(\tilde{b}_{1,NG}\mathcal{P})(\tilde{\kappa}_1)(\tilde{b}_{1,NG}\mathcal{P})(\tilde{\kappa}_2) \\ & + \frac{\bar{\xi}^2}{2}(\tilde{b}_{2,NG}\mathcal{P})(\tilde{\kappa}_1)(\tilde{b}_{2,NG}\mathcal{P})(\tilde{\kappa}_2) \quad (71a) \\ & + \frac{\bar{\xi}^2}{2}[(b_2\mathcal{P})(\tilde{\kappa}_1)(q_1\mathcal{P})(\tilde{\kappa}_2) + (q_1\mathcal{P})(\tilde{\kappa}_1)(b_2\mathcal{P})(\tilde{\kappa}_2)] \\ & + \frac{(\bar{\xi} - \bar{\xi})^2}{2}(b_2\mathcal{P})(\tilde{\kappa}_1)(b_2\mathcal{P})(\tilde{\kappa}_2) \quad (71b) \\ & + \mathcal{O}(\{\bar{\xi}^3, \bar{\xi}\bar{\xi}^2, \bar{\xi}^3, \bar{\xi}(\bar{\xi} - \bar{\xi})^2, (\bar{\xi} - \bar{\xi})^3\}), \end{aligned} $$

where we defined the following NLO composite bias function only present for non-Gaussian fields

$$ \tilde{b}_{2,NG}(\tilde{\kappa}) = b_2(\tilde{\kappa}) + 2 \frac{1}{\mathcal{P}(\tilde{\kappa})} \frac{\partial(b_1 \mathcal{P})}{\partial \tilde{\kappa}} + \frac{1}{\mathcal{P}(\tilde{\kappa})} \frac{\partial^2 \mathcal{P}}{\partial \tilde{\kappa}^2}. \quad (71c) $$

For a Gaussian field, the first three lines in equation (71) do not contribute because $\varphi_1(\lambda) = \lambda$ and $\varphi_{n \ge 2} = 0$, so we recover the leading order of equation (65). Note that while we ordered terms according to powers of $\tilde{\xi}$ here, there is no guarantee that the leading order term proportional to $\tilde{\xi}$ will be larger than the next-to-leading order terms proportional to $\bar{\xi}^2$ and $(\bar{\xi} - \bar{\xi})^2$, as highlighted by the Gaussian case where the leading order term vanishes and the shifted lognormal case where we found the term to be subdominant. Typically, power of averages of the mean correlation are much smaller than averages

<footnote n="3">³ At higher orders, this will generalise to a class of terms with products of positive powers of $\varphi_1(\lambda - \lambda)$ and powers of $\varphi_1(\lambda)^2$ and/or $\varphi_{n \ge 2}(\lambda)$.</footnote>
of powers, so $\bar{\xi}^2 \ll \bar{\xi}^2 \approx (\bar{\xi} - \bar{\xi})^2$ as evident from Figure 12 and for a weakly non-Gaussian field the $\tilde{b}_{2,NG}$ term will only give a small contribution. For the mildly non-Gaussian $\kappa$ PDF we are considering here, we observe that the first order term is effectively removed by the mean-subtraction, as evident from the visual inspection of the covariance matrix shown in Figure 6 and the dominant eigenvectors, which resemble the full sky case shown in green/olive in the lower panel of Figure 16. Given the similarity between the Gaussian and FLASK correlation matrices, we can deduce that the $b_2$ term given by equation (58) for the shifted lognormal case is in fact the dominant contribution.

Our results from this section can largely be carried over to the covariance modelling for the 3-dimensional matter PDF on mildly nonlinear scales (roughly 10Mpc/h at low redshifts), as we will discuss in Section 5 below. Typically, the matter PDF features stronger non-Gaussianities, such that the lognormal shift parameter will be close to unity $s_{3D} \lesssim 1$. In Appendix A we additionally discuss tree models as a useful class of non-Gaussian models for which all order covariance expansion results can be obtained.

## 4. QUANTITATIVE COVARIANCE MODEL TESTS

Beyond performing a visual comparison between the covariances obtained from simulated maps and the theoretical models, we quantitatively compare their accuracy using

1. an eigendecomposition of the covariance validating the set of eigenvalues and leading order eigenvectors,

2. a comparison of the dominant bands in the almost band-diagonal precision matrix as inverse of the covariance,

3. a $\chi^2$-test (advocated in Friedrich et al. 2021) to establish the Gaussianity of the PDF data vector and accuracy of the model covariance, and

4. a Fisher forecast for two cosmological parameters $\{\Omega_m, \sigma_8\}$.

When comparing our measurements to the predictions in tests involving the inverse covariance or precision matrix, we account for that fact that matrix inversion is a non-linear operation, so the noise in the estimate of the covariance elements will lead to a bias in the elements of the precision matrix. This bias from covariance estimation noise turns out to be just a factor multiplying the entire inverse covariance matrix, the so called Kaufman-Hartlap factor (Kaufman 1967; Hartlap et al. 2006),

$$h = (N_{\text{sim}} - N_d - 2)/(N_{\text{sim}} - 1), \quad (72)$$

where $N_{\text{sim}}$ is the number of samples used and $N_d$ is the length of the data vector. This effect closely resembles the impact of the non-Gaussian distribution of a data vector that is a consequence of a finite number of realisations (Sellentin & Heavens 2017) and of size $(N_{\text{sim}} - N_d + N_p - 1)/(N_{\text{sim}} - 1)$, where $N_p$ is the number of parameters. Another effect is due to the nonlinear relationship between the precision matrix and the resulting parameter covariance matrix (Percival et al. 2014),

which can partially counteract the previous effect, but was found to be negligible in our setting. Beyond direct effects for the width of the likelihood contours, noise also leads to an uncertainty for their location in parameter space (Dodelson & Schneider 2013b), which is only negligible if $N_{\text{sim}} - N_d \gg N_d - N_p$ as is the case for our simplified setting with a short data vector.

### 4.1. Eigendecomposition of the covariance matrix

The eigendecomposition of the covariance matrix of dimension $n \times n$ (set by the length $n$ of the data vector) can be written in terms of a complete set of eigenvalues $\lambda_i$ (for convenience ordered in descending order, so $\lambda_1 > \dots > \lambda_n$) and associated normalised eigenvectors $e_i$ arranged in an orthogonal matrix $O$

$$\begin{align} \operatorname{cov}(\mathcal{P}, \mathcal{P}) &= O \cdot \operatorname{diag}(\lambda_1, \dots, \lambda_n) \cdot O^T, \quad O = (e_1, \dots, e_n) \nonumber \\ &= \sum_{k=1}^{n} \lambda_k e_k e_k^T. \tag{73} \end{align}$$

In terms of components, we can write this as

$$\operatorname{cov}(\mathcal{P}_i, \mathcal{P}_j) = \sum_{k=1}^{n} \lambda_k (\mathbf{e}_k \mathbf{e}_k^T)_{ij} = \sum_{k=1}^{n} \lambda_k (\mathbf{e}_k)_i (\mathbf{e}_k)_j, \quad (74)$$

where $(\mathbf{e}_k)_i$ denotes the $i$-th component of the eigenvector $\mathbf{e}_k$. The eigendecomposition can be understood in the context of the analytical results from the large-separation expansion in Section 3. For the Gaussian case, we obtained an expansion consisting of dyadic products with building blocks of the form $He_n(x) \exp(-x^2/2)$ with $x = \kappa/\sigma$ corresponding to $(b_n\mathcal{P})(\kappa)$. While the Hermite polynomials are orthogonal with respect to the Gaussian weight function $\exp(-x^2/2)$ corresponding to the PDF $\mathcal{P}_G(\kappa)$, the dyadic product building blocks are not generally orthogonal as they carry a combined weight of $\exp(-x^2)$. However, the orthogonality of neighbouring bias contributions $b_n\mathcal{P}$ and $b_{n+1}\mathcal{P}$ is ensured by the alternating symmetry of the Hermite functions.

Considering the binned PDF means that there is just a finite set of eigenvectors, which receive contributions from a whole series of $b_n\mathcal{P}$ terms. We illustrate that in Figure 15 by comparing the first three eigenvectors (solid/dashed/dotted) of the PDF covariance for a Gaussian $\kappa$ field obtained from the full theory expression (cyan) with the normalised bias functions (light grey). The leading order bias contribution $b_1\mathcal{P}$ agrees with the first eigenvector almost perfectly, while $b_2\mathcal{P}$ resembles the overall shape of the second eigenvector with slightly modified zero crossings. This agreement reflects that there is a hierarchy in the contributions from different powers of the correlation function compared to the variance $\bar{\xi}^n/\sigma^{2n}$ as shown in Figure 12. The slight difference between $b_2\mathcal{P}$ and the second eigenvector is due to contributions from higher-order terms as we illustrate by showing the eigenvectors of the Hermite expansions (51) up to a finite $n_{\text{max}}$ (shaded from dark grey to cyan with increasing $n_{\text{max}}$). The large difference between $b_3\mathcal{P}$ and the third eigenvector comes from its non-orthogonality with $b_1\mathcal{P}$ as can be seen from the difference between the dotted lines in light and darker grey. For the case of the full sky (or the mean-subtracted $\kappa$ in patches),
<figure><img src="image_18.png" /><figcaption>FIG. 15.— First (solid), second (dashed) and third (dotted) eigenvectors of the Gaussian $\kappa$ PDF covariance predicted for small patches (cyan) and the normalised bias functions $(b_n P)_G(\kappa)$ (light grey). We also show the eigenvectors obtained from a Hermite expansion (51) up to order $n_{max}$ (shaded from dark grey to cyan). The vertical black dotted lines indicate where $\kappa = \pm\sigma$.</figcaption></figure>

the smallness (or vanishing) of the mean correlation $\bar{\xi}$ will practically remove the contribution of $b_1 P$ such that the first and second eigenvectors will be resembling $b_2 P$ and $b_3 P$, respectively. This qualitative statement remains true for mildly non-Gaussian fields such as our shifted lognormal (FLASK) case. We show the eigenvalues and dominant eigenvectors for the Gaussian and FLASK maps in Figure 16. The eigenvalues in the upper panel of Figure 16 show a good agreement between measurements and theory and suggest a clear ordering between contributions from different bias functions arising from a large-separation expansion, qualitatively resembling the behaviour seen in Figure 12 for the averages of powers of the cell correlation in comparison with the powers of the variance. The middle and lower panels of Figure 16 show the leading three eigenvectors for Gaussian (middle) and lognormal (bottom) $\kappa$ fields in small patches (blue/cyan) and the full sky (green/olive). The first eigenvector (solid) is determined by $b_1 P$ and hence shows one zero-crossing close to $\kappa \simeq 0$ as predicted by $b_{1,G}$ (22) and $b_{1,LN}$ (30). The second eigenvector (dashed) is driven by $b_2 P$ with two zero-crossings around $\kappa \simeq \sigma_\kappa$ as predicted by $b_{2,G}$ (50) and $b_{2,G}$ (58). The third eigenvector (dotted) involves $b_3 P$, which is predicted to have three zero-crossings for the Gaussian case from equation (48), but is also sensitive to higher-order terms as we illustrated in Figure 15. For the covariance of the $\kappa$ PDF on the full sky we show the first two eigenvectors (measurements in green, predictions in olive). We observe that the first eigenvector closely resembles the second eigenvector of the $\kappa$ PDF from small patches (measurements in blue, predictions in cyan), as predicted for Gaussian fields in equation (65) and for mildly non-Gaussian fields in equation (71). Overall, the theoretical predictions are in good agreement with the measurements from Gaussian and FLASK maps. The Gaussian and the lognormal cases are qualitatively similar with a few important differences appearing mostly as shifts towards negative $\kappa$ (reflecting the PDF peak location) and an asymmetry between positive and negative $\kappa$ reflect-

<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /><figcaption>FIG. 16.— (Upper panel) Eigenvalues of the Gaussian $\kappa$ PDF covariance for $\kappa$ from small patches (measurement in blue, prediction in cyan) and the full sky (measurement in green, prediction in olive). (Middle and lower panel) The dominant eigenvectors of the Gaussian (middle) and FLASK (lower) $\kappa$ PDF covariance, including the 1st eigenvector (solid lines), 2nd eigenvector (dashed lines) and 3rd eigenvector for raw $\kappa$ from small patches (dotted lines). The move from small patches to the full sky (solid and dashed lines) essentially removes the first term in the eigendecomposition and replaces it with the second term as evident from the agreement of the first green cross with the second blue cross, and the agreement of the solid green and the dashed blue lines. The FLASK results are qualitatively similar to the Gaussian case except for subtle shifts in the location of the zero-crossings and the asymmetry between positive and negative $\kappa$, both related to the non-Gaussianity of the underlying field.</figcaption></figure>

ing the skewness of the PDF. Due to the nature of the eigendecomposition, the contribution of the finite sampling term on the diagonal of the covariance matrix (11) to individual eigenvectors is typically small, such that the leading order results for the non-overlapping case (not shown) resemble the overlapping case aside from
additional noise caused by the small number of sampling cells.

In principle, one can reconstruct the covariance matrix from the terms in the eigendecomposition, although in the presence of an incomplete set of eigenvectors (a number less than the length of the data vector), this matrix will not be invertible and hence unsuitable for statistical analysis. For the case of non-overlapping cells, there is a large contribution to the diagonal from shot noise, so one can attempt a covariance reconstruction using only the first few eigenvectors along with a matched diagonal mostly arising from shot noise. In practice, we opt for the approach to use an invertible covariance baseline model, such as the ones obtained from the Gaussian or shifted lognormal model. If this was found to be insufficiently accurate, it could be successively improved using the precision matrix expansion (Friedrich & Eifler 2017) or iteratively corrected by adjusting the dominant eigenvectors to their desired values as predicted by a large-scale expansion discussed in Section 3.

## 4.2. Bands of the precision matrix

The accuracy of cosmological parameter estimation relies on an accurate modelling of the precision (inverse covariance) matrix which enter the $\chi^2$ tests and Fisher forecasts (Taylor et al. 2013). We show this precision matrix in the case of shifted lognormal fields from the full sky (rescaled to the survey area) in the upper panel of Figure 17 as measured in FLASK maps (upper triangle) and predicted (lower triangle). We notice that the inverse covariance is close to band-diagonal with alternating signs. In the lower panel we focus on the three dominant bands (black, blue, red) comparing the prediction (solid) to the measurements (dotted). We find similar results for the precision matrix shape obtained from FLASK patches, and for both the raw $\kappa$ and the mean-subtracted $\tilde{\kappa} = \kappa - \bar{\kappa}$ from overlapping cells. We can qualitatively understand the agreement between $\kappa$ and $\tilde{\kappa}$ in terms of the previously discussed eigendecomposition. In principle, the precision matrix can be obtained from the eigenvectors and eigenvalues of the covariance matrix as

$$ \mathrm{cov}^{-1}(\mathcal{P}, \mathcal{P}) = O \cdot \mathrm{diag}(\lambda_1^{-1}, \dots, \lambda_n^{-1}) \cdot O^T, $$

which means eigenvectors are retained but their order is reversed. From the upper panel of Figure 16 we expect a flat hierarchy of the highest eigenvalues that determine the inverse, with the lowest order eigenvectors (associated with the largest eigenvalues) playing a minor role for the shape of the precision matrix.

The band-diagonal structure of the precision matrix is mainly driven by the strong correlations of neighbouring bins induced by the overlap of cells arising at small separations. This can be inferred from a comparison with the non-overlapping case and could be formalised by evaluating the small separation contributions described in Bernardeau (2022) for the analogous case of densities. We can illustrate this effect in a simpler setting by computing the covariance between 1D self-centered cells for a Poisson (constant) power spectrum with no correlation between non-overlapping cells. Thus for two cells of half-length $x_i$ and $x_j$, their covariance is given by the integral over the (constant) power spectrum multiplied by the 1D

<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /><figcaption>FIG. 17.— (Upper panel) Impression of the approximately band-diagonal precision matrix obtained as inverse of the $\kappa$ covariance matrix for overlapping cells for shifted lognormal fields as predicted (lower triangle) and measured in 500 FLASK full sky realisations corrected by the Hartlap factor (72) (upper triangle). This matrix is close to penta-diagonal, featuring 5 central bands, one diagonal band (black) and 2 neighbouring bands above/below (blue, red). We highlight those bands with solid/dotted lines for the measured and predicted case. On the diagonal we show the mean of the FLASK measurements and the prediction. The measurement from a finite number of realisations in the upper panel exhibits additional noise off the diagonal. (Lower panel) Comparison of the bands in the precision matrix between the shifted lognormal theory prediction (dotted lines) and the measurement from FLASK realisations (solid lines).</figcaption></figure>

top-hat window function of each cell. We thus have

$$ \mathrm{cov}_{1D}(x_i, x_j) \propto \int dk \frac{\sin(kx_i) \sin(kx_j)}{kx_i kx_j} \propto [\max(x_i, x_j)]^{-1}. \qquad (75a) $$

This special covariance shape is called a hook matrix

$$ \mathrm{cov}_{1D} \propto \begin{pmatrix} 1/x_1 & 1/x_2 & 1/x_3 & \dots & 1/x_n \\ 1/x_2 & 1/x_2 & 1/x_3 & \dots & 1/x_n \\ 1/x_3 & 1/x_3 & 1/x_3 & \dots & 1/x_n \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1/x_n & \dots & 1/x_n \end{pmatrix}. \qquad (75b) $$

This matrix can be inverted by obtaining the cofac-
tor matrix $C$ which determines the inverse as $cov_{1D}^{-1} = C^T/\det cov_{1D}$. The elements of the cofactor matrix are given by Cramer's rule in terms of determinants of sub-matrices obtained from $cov_{1D}$ by deleting the $i$th row and $j$th column

$$C_{i,j} = (-1)^{i+j} \det(\mathrm{cov}_{1D}|_{i,j}). \quad (76)$$

It can then be shown that indeed the cofactor elements vanish everywhere except on the diagonal and the directly adjacent bands to the main diagonal. This is because the removal of a row $i$ and column $j$ in the covariance (75) leads to linearly dependent rows/columns and hence $C_{i,j} = 0$ for $|i-j| < 2$. As the covariance is symmetric, so is the cofactor matrix and hence the inverse of the covariance will have a band-diagonal structure. When the values of the radii $x_i$ (resembling $\kappa_i$) are an increasing sequence, the diagonal will have a positive sign and the band adjacent to the diagonal will have a negative sign.

### 4.3. $\chi^2$ tests

We will use a $\chi^2$-test (advocated in Friedrich et al. 2021) to establish the Gaussianity of the data vector and assess the accuracy of our predicted covariances. Let us consider the $\chi^2$ associated with the original $N$ samples of the data vector. For every sample, labelled by $i$, of the data vector, $D_{\mathrm{org}}[i]$, of length $n$ (encoding the number of degrees of freedom), one computes the $\chi^2$

$$\chi^2_{\mathrm{org}}[i] = (D_{\mathrm{org}}[i] - \mu)^T \cdot C^{-1} \cdot (D_{\mathrm{org}}[i] - \mu), \quad (77)$$

using the mean data vector, $\mu$ (also of length $n$) and data covariance $C$ (of dimension $n \times n$). The distribution of $\chi^2$ follows a $\chi^2$ distribution with the appropriate number of degrees of freedom $n$ (here the number of PDF bins). To test the Gaussianity of the data vector distribution, one can create another sample of data vectors, $D_G$, by drawing $N$ times from a multivariate Gaussian using the mean data vector, $\mu$ and covariance $C$. Then one recomputes the $\chi^2$ for every sample using the covariance to check whether it follows the same distribution

$$\chi^2_G[i] = (D_G[i] - \mu)^T \cdot C^{-1} \cdot (D_G[i] - \mu). \quad (78)$$

If so, one can assert that the data vector is 'close enough' to Gaussian distributed around the mean with the measured covariance. Both of those quantities will have the mean $\overline{\chi^2} = \overline{\chi^2_{\mathrm{org}}} = n$, since the covariance matrix is identical and is reproduced in the expectation value of the fluctuations of the data vector around the mean.

One can also use a $\chi^2$ test to compare the closeness of two covariances, where $C$ denotes the original 'true' (or target) covariance and $C_{\mathrm{mod}}$ the model (or approximated) covariance. By recomputing the $\chi^2$ using the model covariance

$$\chi^2_{\mathrm{mod}}[i] = (D_{\mathrm{org}}[i] - \mu)^T \cdot (C_{\mathrm{mod}})^{-1} \cdot (D_{\mathrm{org}}[i] - \mu) \quad (79)$$

and checking its distribution compared to the original one $\chi^2_{\mathrm{org}}$ mentioned above, one can assess how well the model and true covariance agree. We found that the covariance of non-overlapping cells from the FLASK patches is well captured in both the Gaussian and shifted lognormal theory, as it is dominated by the finite sampling term and thus not shown here. The covariance of

<figure><img src="image_23.png" /></figure>

<figure><img src="image_24.png" /><figcaption>FIG. 18.— $\chi^2$-tests for the validation of our predicted covariances and the Gaussianity of our data vector showing the $\chi^2$ distribution (black line) and the result for a data vector drawn from a Gaussian with the measured covariance (black histogram). (Upper panel) Covariance matrix test for the PDF of the convergence $\kappa$ (blue) or mean-subtracted convergence $\kappa - \bar{\kappa}$ (purple) measured from overlapping cells computed in small patches comparing the FLASK measurements (solid) to the predictions from a shifted lognormal (dotted). (Lower panel) Covariance matrix test for overlapping cells computed for the full sky comparing the measurements from 500 FLASK realisations (solid) to the lognormal prediction (dotted).</figcaption></figure>

heavily overlapping cells from FLASK patches was visually different from the Gaussian case, mainly driven by the asymmetry in the diagonal of the covariance matrix shown in Figure 5. In Figure 18 we show the $\chi^2$ test results for the overlapping case with FLASK covariances for patches (upper panel) and the full sky (lower panel). The agreement of the $\chi^2$ obtained from the measured covariances (coloured solid) and a draw from a Gaussian with the same variance (black) shows that the data vectors are close to Gaussian distributed. We find good agreement between the FLASK covariance measurements (coloured solid) and the shifted lognormal predictions (coloured dotted). The results for patches (upper panel) show a slight (6-9%) bias of the theory towards higher $\chi^2$, which could be caused by the patch-to-patch correlation arising from cutting multiple patches from each FLASK full sky that is unaccounted for in the theory. Interestingly, we find no evidence for a corresponding discrepancy in the Fisher forecasts for those two cases, demon-
strating the complementarity of the two tests adopted in Friedrich et al. (2021).

### 4.4. Fisher forecast validation

To forecast the errors on a set of cosmological parameters, $\vec{p}$, we use the Fisher matrix formalism. The Fisher matrix given a (set of) summary statistics in the data vector $\vec{D}$ is defined as

$$F_{ij} = \sum_{\alpha, \beta} \frac{\partial D_{\alpha}}{\partial p_i} (C^{-1})_{\alpha \beta} \frac{\partial D_{\beta}}{\partial p_j}, \quad (80)$$

where $D_\alpha$ is the $\alpha$-th element of the data vector $\vec{D}$, $C^{-1}$ denotes the precision matrix introduced in Section 4.2, which is the matrix-inverse of the data covariance $C$ given in equation (17) and is assumed that the covariance matrix is independent of cosmology. The Fisher formalism rests on the assumption that the summary data vector is Gaussian distributed which we have explicitly checked in the last subsection for the central region of the weak lensing PDF considered here. The parameter covariance matrix $C(\vec{p})$ is then obtained as inverse of the Fisher matrix, $C_{ij}^P = (F^{-1})_{ij}$. In the Fisher formalism, marginalisation over a subset of parameters is achieved by simply selecting the appropriate sub-elements of the parameter covariance.

To further validate our theoretical covariances, we perform simple Fisher forecasts for two cosmological parameters {$\Omega_{cdm}, \sigma_8$} with the measured and predicted covariance matrices, while keeping the data vector and its dependence on cosmological parameters in (80) as predicted from large-deviation statistics (see Boyle et al. 2021).<fnref n="4" /> The Fisher forecasts in Figure 19 demonstrate that measuring non-overlapping cells only (cyan) leads to much wider parameter errors than including cell overlaps (blue), such that capturing cell overlaps is imperative to avoid artificially inflating the error budgets (as already pointed out in Szapudi & Colombi 1996). Additionally we show that when measuring covariances from small patches obtained from FLASK and subtracting the mean in each small patch (purple) this reproduces the result obtained from full sky maps rescaled to the survey area (green). In Figure 20 we show that our shifted lognormal predicted covariances (dotted) for small patches with mean-subtraction (purple) and the full sky (green) agree well with the FLASK measurements (solid). We also illustrate the impact of shape noise (orange) and validate our convolution approach from equation (42) against FLASK measurements with added shape noise.

## 5. EXTENSION TO 3D CLUSTERING PDF COVARIANCE

In this Section, we describe the essential steps needed to adapt our covariance model to 3-dimensional fields using the 3D matter PDF as an example. While the line-of-sight projection involved in computing the weak

<figure><img src="image_25.png" /><figcaption>FIG. 19.— Fisher forecasts comparing different FLASK-generated covariances using non-overlapping cells (cyan), overlapping cells (blue), overlapping cells with mean-subtraction (purple) and overlapping cells on full sky maps (green). All covariances were rescaled to mimic the same survey area.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>FIG. 20.— Fisher forecasts validating the predicted covariances (dotted) against FLASK measurements (solid) for shifted lognormal fields on small patches (purple), the full sky (green) and including shape noise (orange).</figcaption></figure>

lensing convergence and photometric galaxy density on mildly nonlinear scales tends to Gaussianise the underlying field, 3D matter and spectroscopic galaxy densities typically show a greater degree of non-Gaussianity on mildly nonlinear scales. This might complicate the two-

<footnote n="4">⁴ The predicted derivatives rely on large-deviations theory for the convergence PDF (Barthelemy et al. 2020a). Using CLASS (Blas et al. 2011), linear and non-linear (Halofit, Peacock & Smith 2014) matter density variances are calculated to rescale the matter density CGF in each redshift slice. The convergence CGF is obtained by integration over the redshift slices, and translated to the final PDF by an inverse Laplace transform.</footnote>
<figure><img src="image_27.png" /><figcaption>FIG. 21.— 3D matter density PDF at redshift $z = 0$ and smooth- ing radius $R = 10$ Mpc/h as measured from the Quijote simula- tions (data points, with error bars indicating the sample variance), predicted by large-deviation theory (black line) and approximated with a lognormal distribution (black dotted) and a shifted lognor- mal distribution with tuned shift parameter (black dashed).</figcaption></figure>

point PDF modelling, but we will show that the previ- ously discussed shifted lognormal model from Section 2.3 has flexibility to match the skewness of the one-point distribution and capture the clustering properties. Fi- nally, we will focus on the theoretical modelling of the super-sample covariance, which can be used to comple- ment simulated covariances from periodic boxes missing this effect.

### 5.1. Applicability of the shifted lognormal model

The three-dimensional matter PDF is typically ex- tracted in terms of the normalised cell density $\hat{\rho} = \rho/\rho_s$ as ratio of the physical cell density $\rho$ and the mean sample density $\rho_s$, or the associated zero-mean density contrast $\delta = \hat{\rho} - 1$. For illustration purposes, we use previously obtained measurements for the matter PDF at redshift $z = 0$ in spheres of radius $R = 10$ Mpc/h (Uhlemann et al. 2020) in the Quijote simulation suite (Villaescusa-Navarro et al. 2019). The one-point PDF of the normalised density<fnref n="5" /> $\hat{\rho}$ was extracted in logarithmi- cally spaced bins designed to equally capture under- and overdensities in the strongly non-Gaussian distribution and means and covariances obtained from 15,000 realisa- tions at the fiducial cosmology.

In Figure 21 we give a visual impression of the mean PDF for the normalised matter density $\hat{\rho}$ as measured from the simulations (data points) and large-deviation theory prediction (black line, following Uhlemann et al. 2016, 2020). To describe clustering observables with the shifted lognormal model, we can simply replace the weak lensing convergence $\kappa$ with the density contrast $\delta$ in equa- tions (26) and (28). Choosing a unit shift parameter $s = 1$ corresponds to the usual lognormal approximation (dotted black line, Coles & Jones 1991). This already provides a decent phenomenological fit to first princi- ple predictions for the matter PDF from large deviation statistics (Bernardeau et al. 2014; Uhlemann et al. 2016,

2020), and can be improved by including a shift factor below unity (dashed black line).<fnref n="6" /> The good match for the fiducial PDF in Figure 21 and the leading order bias in Figure 23 suggests that the shifted lognormal model can be used to obtain accurate analytical PDF covari- ance estimates. This can be done by computing the joint PDF (82) following a similar lognormal matching proce- dure as sketched for the mean-subtracted weak lensing convergence $\tilde{\kappa}$. For 3D observables for which the shifted lognormal model would not yield a good match of the PDF and leading order bias, the class of tree models dis- cussed in Appendix A and Bernardeau (2022) offer an alternative with the capability of skewness matching.

### 5.2. Relevance and prediction of super-sample covariance

*Measurements.* The covariance measured from simu- lated boxes with fixed mean density (due to a fixed par- ticle number) $\tilde{\rho}$ cannot capture the super-sample effect created by having different sample densities $\rho_s$ as would arise naturally from computing the covariance of sub- boxes within a much larger simulation box. To estimate this background density driven super-sample covariance effect for the Quijote simulation suite, we make use of their publicly available separate-universe style ‘DC’ runs emulating a background density contrast $\delta_b = \pm 0.035$ through changed cosmological parameters and simula- tion snapshot times from the separate universe approach (Sirko 2005). The super-sample covariance between two data vector entries $D_i$ and $D_j$ can be estimated by

$$ \operatorname{cov}_{\mathrm{SU}}^{\mathrm{SSC}}(D_i, D_j) = \sigma_b^2 \frac{\partial D_i}{\partial \delta_b} \frac{\partial D_j}{\partial \delta_b}, \qquad (81) $$

where $\sigma_b^2$ is the variance of $\delta_b$ (here just $\delta_b^2$), and the other two terms encode the linear response of the data vector, which can be determined from the simulations using fi- nite differences. The impact of super-sample covariance on the 3D matter and halo power spectra, bispectra, void size distributions and the halo mass function has been recently studied in Bayer et al. (2022). As mentioned therein, in addition to the background density effect cap- tured by the separate universe approach, there can be an additional super-sample effect due by tidal fields affect- ing higher-order statistics in a non-trivial way. For the PDF, we expect those tidal effects to be negligible due to its intrinsic symmetry and averaging.

In the lower triangle of the lower panel of Figure 22 we show the covariance matrix of the matter density PDF showing the correlation between different density bins in the central region of the PDF. We see that, as ex- pected, neighbouring bins are positively correlated, while intermediate underdense and overdense bins are anticor- related with each other. The strong positive correlation in the bands along the diagonal is induced by cell over- laps, similarly as discussed previously for the analogous weak lensing case. We illustrate the separate universe super-sample effect on the PDF covariance in the up- per triangle of Figure 22, showing that the super-sample covariance term completely dominates the overall covari-

<footnote n="5">⁵ ensuring the compatibility of measurements in simulations with different particle numbers (changing the mean density of particles in a sphere).</footnote>

<footnote n="6">⁶ The lognormal model is however insufficient to accurately predict the response to cosmological parameters (see Uhlemann et al. 2020).</footnote>
<figure><img src="image_28.png" /><figcaption>FIG. 22.— Comparison of the PDF covariance matrix contributions for $z = 0$, $R = 10$ Mpc/h as measured from the 1(Gpc/h)$^3$ Quijote simulations and the super-sample covariance contribution constructed from the DC mode following equation (81). The upper triangle shows the super-sample covariance contribution while the lower triangle shows the measured covariance from the simulations enhanced by a factor of 10, while the diagonal from the upper panel is masked. The dotted lines indicate the PDF peak location.</figcaption></figure>

ance for a volume of $1(Gpc/h)^3$ and will still lead to a significant effect for a volume of $\mathcal{O}(10)(Gpc/h)^3$.

**Predictions.** Predicting the covariance of the PDF of the normalised matter density $\hat{\rho}$ in a three-dimensional volume is conceptually analogous to the previous two-dimensional case relevant for weak lensing, so we can write the covariance in terms of the two-point PDF as in equation (6) by replacing the weak lensing convergence with the density $\kappa \rightarrow \rho$ and angular by 3D-distances $\theta \rightarrow r$, where $P_d(r)$ indicates the distribution of distances $r$ in a given survey volume. The distance distribution in a cubic box is easily obtained from the one for a 3D unit cube that corresponds to the problem of cube line picking for which a closed-form expression is available (Mathai et al. 1999; Weisstein 2022). Similarly to the weak lensing case, where the mean-subtraction leads to a modification of the relevant bias functions determining the covariance, the relative density definition will result in modified bias functions. This can be predicted by considering the matter PDF in different sample densities $\rho_s$, where the joint PDF of normalised densities $\hat{\rho}_i$ is obtained from the trivariate PDF of physical densities $\rho_i = \rho_s \hat{\rho}_i$ and the background density $\rho_s$ by integrating out the latter

$$ \mathcal{P}(\hat{\rho}_1, \hat{\rho}_2) = \int d\rho_s \rho_s^2 \mathcal{P}(\rho_s, \rho_1 = \rho_s \hat{\rho}_1, \rho_2 = \rho_s \hat{\rho}_2). \quad (82) $$

The super-sample covariance effect (13) for the PDF of the normalised matter density $\hat{\rho}$ is obtained from the leading order effective bias function in the resulting covariance, given by equation (28) in Bernardeau (2022)

$$ \hat{b}_1(\hat{\rho}) = 1 + \frac{\partial \log \mathcal{P}(\hat{\rho})}{\partial \log \hat{\rho}} + b_1(\hat{\rho}). \quad (83) $$

Note that while the expression looks similar to the one

for the mean-subtracted weak lensing case $\tilde{b}_{1,NG}(\tilde{\kappa})$ from equation (68b), this contribution does not vanish even for a Gaussian field. Predictions for the ingredients, the one-point PDF $\mathcal{P}(\rho)$ and the sphere bias $b_1(\rho)$, can be obtained from large-deviation statistics (Codis et al. 2016b; Uhlemann et al. 2017), shifted lognormal models (Section 2.3) and tree models (Appendix A and Bernardeau 2022).

**Comparison.** In Figure 23 we demonstrate that the derivative of the PDF with respect to the background density (green data points) determining the super-sample covariance (81) is well predicted by the effective bias $\hat{b}_1$ from equation (83) using theoretical ingredients (green solid line) or simulated ones (green dotted line). We also show the $\hat{b}_1$ ingredients consisting of the first order sphere bias function $b_1$ as measured from the correlation between neighbouring spheres in one realisation (blue data points) in comparison to the large-deviation theory prediction (blue solid line, following Uhlemann et al. 2017) and the lognormal approximation (blue dashed line), and the logarithmic derivative of the one-point PDF as measured from finite differences (red data points) and predicted from large-deviation theory (red solid line). We conclude that the super-sample covariance effect can be robustly predicted by the theory presented here and that the shifted lognormal model holds promise to compute analytical covariances for the PDF of 3D clustering observables.

**Applications.** The total covariance matrix is a sum of the simulated covariance, $C_{sim}$ and the super-sample covariance term consisting of a dyadic product, $C_{SSC} = v_{SSC}v_{SSC}^\top$ with $v_{SSC} = \sqrt{\xi}(\hat{b}_1\mathcal{P})(\hat{\rho})$. This sum can be inverted using the Sherman-Morrison formula (Sherman & Morrison 1950) to obtain the precision matrix

$$ (C_{\text{sim}} + C_{\text{SSC}})^{-1} = C_{\text{sim}}^{-1} - \frac{C_{\text{sim}}^{-1} C_{\text{SSC}} C_{\text{sim}}^{-1}}{1 + v_{\text{SSC}}^{\top} C_{\text{sim}}^{-1} v_{\text{SSC}}}. \quad (84) $$

The result can be used to predict the super-sample covariance impact on statistical tests like $\chi^2$ and Fisher contours, which generally leads to a degradation of parameter constraints as detailed in (Lacasa & Grain 2019) for the case of Euclid-like photometric galaxy clustering observables.

## 6. CONCLUSION & OUTLOOK

**Conclusion.** In this work, we have shown how to compute covariances for the one-point PDF of the weak lensing convergence on mildly nonlinear scales. Our covariance predictions (6) rely on the two-point PDF, which is integrated over all separations according to the distribution of distances in the survey area. Our formalism includes predictions for the key effects including finite sampling (11) and super-sample covariance (13) which were shown to scale with the inverse survey area. We presented effective models for this joint PDF for a Gaussian field (19), a shifted lognormal field (28) matching the PDF skewness and density-dependence of the two-point clustering (Figure 8), and a field affected by shape noise (42) from galaxy shape measurements. Going beyond (log-)normal models, we laid out a large-separation expansion for the two-point PDF that captures the key effects of the PDF covariance for general non-Gaussian
<figure><img src="image_29.png" /><figcaption>FIG. 23.— Comparison of the theoretically predicted sphere bias (blue solid), the lognormal approximation (blue dashed) and the measurement in one realisation of the fiducial cosmology of Qujote (blue data points) at radius $R = 10$ Mpc/h, distance $d = 20$ Mpc/h and redshift $z = 0$. Additionally shown is the logarithmic derivative of the PDF (red) as predicted (solid) and measured (data points) and the resulting prediction from the effective bias (83) using the theoretical ingredients (green solid) or the simulated ones (green dotted). The predictions agrees very well with the measurements from separate universe simulations (green data points). The vertical dotted line indicates the peak location of the PDF. The $\rho$-range shown equals the one used in the covariance plots in Figure 22.</figcaption></figure>

fields. The validity of our predictions has been established through a set of tests including visual inspection of the covariance, leading order eigenvalues and eigenvectors (Figure 16), the precision matrix (Figure 17), $\chi^2$ tests (Figure 18) and a two-parameter Fisher forecast (Figure 20). With this suite of tests we have demonstrated that our PDF covariance model not only accurately describes large sets of simulated FLASK (Xavier et al. 2016) maps, but also measurements of realistic full sky convergence maps obtained from full N-body simulations (Takahashi et al. 2017).

Finally, we have sketched how the formalism described here for mildly non-Gaussian weak lensing observables is also applicable to more strongly non-Gaussian 3D clustering quantities. Most importantly, we showed how this can yield predictions for the *super-sample covariance effect* (Figures 22 and 23) to complement covariance measurements from finite-box simulations.

**Outlook.** We expect our formalism can be straightforwardly adapted to describe the one-point PDF of related weak-lensing observables such as the CMB weak lensing convergence (Barthelemy et al. 2020b) and aperture mass (Barthelemy et al. 2021), as well as photometric galaxy clustering. Our model can be extended to describe the cross-correlations between weak lensing PDFs at different smoothing scales and/or tomographic redshift bins. Based on our PDF estimator, one can further attempt to model cross-covariances between the one-point PDF and 2-point statistics, which in the mildly non-linear regime are dominated by the correlation between the PDF and the variance (see Figure 12 in Uhlemann et al. 2020). Note that our covariance model and its possible extensions can serve as a starting point for improving covariance estimates from simulations, even in regimes where the model itself becomes inaccurate. The Gaussian and shifted lognormal covariance models can be used as a starting point for a *precision matrix expansion* (Friedrich & Eifler 2017) in which the leading order eigenvectors could iteratively refine this baseline covariance model. Additionally, analytical models can be used to construct priors for a Bayesian covariance estimation with the help of cheap surrogates (Chartier & Wandelt 2022).

## ACKNOWLEDGEMENTS

This work has made use of the Infinity Cluster hosted by Institut d'Astrophysique de Paris. We thank Stéphane Rouberol for smoothly running this cluster for us. We thank the authors of Takahashi et al. (2017) for their publicly available lensing maps, Villaescusa-Navarro et al. (2019) for their publicly available simulation suite products, and Xavier et al. (2016) for their publicly available FLASK code. The figures in this work were created with `MATPLOTLIB` (Hunter et al. 2007) making use of the `NUMPY` (Harris et al. 2020), `SCIPY` (Vitanen et al. 2020), `SCI-KIT IMAGE` (van der Walt et al. 2014) and `ChainConsumer` Python packages.

CU is supported by the STFC Astronomy Theory Consolidated Grant ST/W001020/1 from UK Research & Innovation. CU's work was performed in part at Aspen Center for Physics, which is supported by National Science Foundation grant PHY-1607611 and partially supported by a grant from the Simons Foundation. AG is supported by an EPSRC studentship under Project 2441314 from UK Research & Innovation. SC thanks Fondation Merac and the French *Agence Nationale de la Recherche* (grant ANR-18-CE31-0009) for funding. The authors thank the Euclid work packages for additional galaxy clustering probes and higher-order weak lensing statistics, as well as the Rubin LSST:DESC higher-order statistics group for useful discussions.

## REFERENCES

* Amendola L., et al., 2018, Living Reviews in Relativity, 21

* Baratta P., Bel J., Plaszczynski S., Ealet A., 2020, A&A, 633, A26

* Bardeen J. M., Bond J. R., Kaiser N., Szalay A. S., 1986, ApJ, 304, 15

* Barthelemy A., Codis S., Uhlemann C., Bernardeau F., Gavazzi R., 2020a, MNRAS, 492, 3420

* Barthelemy A., Codis S., Bernardeau F., 2020b, MNRAS, 494, 3368

* Barthelemy A., Codis S., Bernardeau F., 2021, MNRAS, 503, 5204

* Bayer A. E., Liu J., Terasawa R., Barreira A., Zhong Y., Feng Y., 2022, Super-sample covariance of the power spectrum, bispectrum, halos, voids, and their cross-covariances, doi:10.48550/ARXIV.2210.15647,

https://arxiv.org/abs/2210.15647

* Bernardeau F., 1996, A&A, 312, 11

* Bernardeau F., 2022, arXiv e-prints, p. arXiv:2207.02500

* Bernardeau F., Reimberg P., 2016, Phys. Rev. D, 94, 063520

* Bernardeau F., Valageas P., 2000, A&A, 364, 1

* Bernardeau F., Pichon C., Codis S., 2014, Phys. Rev. D, 90, 103519

* Bernardeau F., Codis S., Pichon C., 2015, MNRAS, 449, L105
Bertschinger E., 1987, ApJ, 323, L103

Blas D., Lesgourgues J., Tram T., 2011, Journal of Cosmology and Astroparticle Physics, 2011, 034-034

Boyle A., Uhlemann C., Friedrich O., Barthelemy A., Codis S., Bernardeau F., Giocoli C., Baldi M., 2021, Monthly Notices of the Royal Astronomical Society, 505, 2886

Chartier N., Wandelt B. D., 2022, MNRAS, 515, 1296

Clerkin L., et al., 2016, Monthly Notices of the Royal Astronomical Society, 466, 1444-1461

Codis S., Pichon C., Bernardeau F., Uhlemann C., Prunet S., 2016a, MNRAS, 460, 1549

Codis S., Bernardeau F., Pichon C., 2016b, MNRAS, 460, 1598

Colavincenzo M., et al., 2019, MNRAS, 482, 4883

Coles P., Jones B., 1991, MNRAS, 248, 1

Desjacques V., Jeong D., Schmidt F., 2018, Phys. Rep., 733, 1

Dodelson S., Schneider M. D., 2013a, Phys. Rev. D, 88, 063537

Dodelson S., Schneider M. D., 2013b, Phys. Rev. D, 88, 063537

Friedrich O., Eifler T., 2017, Monthly Notices of the Royal Astronomical Society, 473, 4150

Friedrich O., et al., 2018, Phys. Rev. D, 98, 023508

Friedrich O., Uhlemann C., Villaescusa-Navarro F., Baldauf T., Manera M., Nishimichi T., 2020, MNRAS, 498, 464

Friedrich O., et al., 2021, MNRAS, 508, 3125

Gatti M., et al., 2019, arXiv e-prints, p. arXiv:1911.05568

Giocoli C., Baldi M., Moscardini L., 2018, MNRAS, 481, 2813

Gruen D., et al., 2018, Phys. Rev. D, 98, 023507

Halder A., Friedrich O., Seitz S., Varga T. N., 2021, MNRAS, 506, 2780

Harnois-Déraps J., et al., 2018, MNRAS, 481, 1337

Harris et al. C. R., 2020, Nature, 585, 357

Hartlap J., Simon P., Schneider P., 2006, Astronomy & Astrophysics, 464, 399-404

Hilbert S., Hartlap J., Schneider P., 2011, Astronomy & Astrophysics, 536, A85

Hoffman Y., Ribak E., 1991, ApJ, 380, L5

Hou J., Cahn R. N., Philcox O. H. E., Slepian Z., 2021, arXiv e-prints, p. arXiv:2108.01714

Hunter et al. J. D., 2007, Computing in Science & Engineering, 9, 90

Ivezić Ž., et al., 2019, ApJ, 873, 111

Jamieson D., Loverde M., 2020, Physical Review D, 102

Kaiser N., 1984, ApJ, 284, L9

Kaufman G. M., 1967, Report No. 6710, Center for Operations Research and Econometrics, Catholic University of Louvain, Heverlee, Belgium

Lacasa F., Grain J., 2019, Astronomy & Astrophysics, 624, A61

Laureijs R., et al., 2011, preprint, (arXiv:1110.3193)

Liu J., Bird S., Matilla J. M. Z., Hill J. C., Haiman Z., Madhavacheril M. S., Petri A., Spergel D. N., 2018, Journal of Cosmology and Astroparticle Physics, 2018, 049-049

Martin S., Schneider P., Simon P., 2012, Astronomy & Astrophysics, 540, A9

Mathai A. M., Moschopoulos P., Pederzoli G., 1999, Statistica, 59, 61

Mo H. J., White S. D. M., 1996, MNRAS, 282, 347

Peacock J. A., Smith R. E., 2014, HALOFIT: Nonlinear distribution of cosmological mass and galaxies (ascl:1402.032)

Percival W. J., et al., 2014, MNRAS, 439, 2531

Percival W. J., Friedrich O., Sellentin E., Heavens A., 2022, MNRAS, 510, 3207

Repp A., Szapudi I., 2021a, MNRAS, 500, 3631

Repp A., Szapudi I., 2021b, Monthly Notices of the Royal Astronomical Society, 509, 586

Sarma Boruah S., Rozo E., Fiedorowicz P., 2022, arXiv e-prints, p. arXiv:2204.13216

Sellentin E., Heavens A. F., 2017, MNRAS, 464, 4658

Sherman J., Morrison W. J., 1950, The Annals of Mathematical Statistics, 21, 124

Sheth R. K., Tormen G., 1999, MNRAS, 308, 119

Sirkio E., 2005, ApJ, 634, 728

Szapudi I., Colombi S., 1996, ApJ, 470, 131

Szapudi I., Colombi S., Bernardeau F., 1999, Monthly Notices of the Royal Astronomical Society, 310, 428

Takahashi R., Hamana T., Shirasaki M., Namikawa T., Nishimichi T., Osato K., Shiroyama K., 2017, The Astrophysical Journal, 850, 24

Taylor A., Joachimi B., 2014, MNRAS, 442, 2728

Taylor A., Joachimi B., Kitching T., 2013, MNRAS, 432, 1928

Thiele L., Hill J. C., Smith K. M., 2020, Phys. Rev. D, 102, 123545

Uhlemann C., Codis S., Pichon C., Bernardeau F., Reimberg P., 2016, MNRAS, 460, 1529

Uhlemann C., Codis S., Kim J., Pichon C., Bernardeau F., Pogosyan D., Park C., L'Huillier B., 2017, MNRAS, 466, 2067

Uhlemann C., Pichon C., Codis S., L'Huillier B., Kim J., Bernardeau F., Park C., Prunet S., 2018, MNRAS, 477, 2772

Uhlemann C., Friedrich O., Villaescusa-Navarro F., Banerjee A., Codis S., 2020, MNRAS, 495, 4006

Villaescusa-Navarro F., et al., 2019, arXiv e-prints, p. arXiv:1909.05273

Virtanen et al. P., 2020, Nature Methods, 17, 261

Weisstein E. W., 2022, Cube Line Picking. From MathWorld—A Wolfram Web Resource, http://mathworld.wolfram.com/CubeLinePicking.html

Xavier H. S., Abdalla F. B., Joachimi B., 2016, MNRAS, 459, 3693

van der Walt et al. S., 2014, PeerJ, 2, e453

APPENDIX

A. LARGE-SEPARATION EXPANSION IN THE MINIMAL TREE MODEL

In addition to shifted lognormal models used heavily in the main text, let us discuss one representative from the class of so-called tree models considered in Bernardeau (2022) in the context of PDF covariance modelling. The minimal tree model is the simplest example for which analytical results can be obtained at all orders. Here, we show how the phenomenology of tree models resembles predictions from large-deviation theory for the matter PDF (see the upper panel of Figure 24). Additionally, we present a closed-form expression for the n-th order bias parameter $b_{n,t}$ (A14) that can be used to construct the full large separation series expansion yielding an invertible analytical covariance model.

The one-point cumulant generating function of the minimal tree model is given by

$$ \varphi_t(\lambda) = \frac{\lambda}{1 - \lambda\sigma^2/2} \approx \lambda + \frac{\sigma^2}{2}\lambda^2 + \underbrace{\frac{\sigma^4}{4}}_{\kappa_3/3!} \lambda^3 + \mathcal{O}(\sigma^6), \quad (A1) $$

where the expansion shows that the mean is one, $\sigma^2$ is the variance$^7$ and the reduced skewness is $S_3 = \kappa_3/\sigma^4 = 3/2$. The one-point PDF of the matter density in the minimal tree model can be computed in closed form and is given by a hypergeometric function

$$ \mathcal{P}_t(\rho) = \frac{4}{\sigma^4} \exp \left[ -\frac{2}{\sigma^2} (1 + \rho) \right] {}_0F_1 \left( 2, \frac{4\rho}{\sigma^4} \right) . \quad (A2) $$

The two-point PDF can be obtained from the two-point CGF which can be written in terms of the individual one-point CGFs

$$ \varphi_t(\lambda_1, \lambda_2) = \frac{\varphi(\lambda_1) + \varphi(\lambda_2) + \xi_{12}\varphi(\lambda_1)\varphi(\lambda_2)}{1 - \xi_{12}^2\varphi(\lambda_1)\varphi(\lambda_2)/4} . \quad (A3a) $$

By expanding that expression in the correlation $\xi_{12}$, one can obtain the two leading order bias functions and the additional $q_1$ term in the large-separation covariance ex-

<footnote>$^7$ Note that Bernardeau (2022) uses $\bar{\xi}$ to denote the variance, which should not be confused with our average correlation $\xi$ from equation (13b).</footnote>
pansion (56)

$$b_{1,t}(\rho) = \frac{0F_1\left(1, \frac{4\rho}{\sigma^4}\right)}{0F_1\left(2, \frac{4\rho}{\sigma^4}\right)} - \frac{2}{\sigma^2}, \qquad (A4a)$$

$$\begin{aligned} b_{2,t}(\rho) &= \frac{4[\rho - 1 - \sigma^2 b_1(\rho)]}{\sigma^4} = \frac{4}{\sigma^2} \left[ \frac{\rho - 1}{\sigma^2} - b_1(\rho) \right] \\ &= \frac{4}{\sigma^2} \left[ \frac{\rho + 1}{\sigma^2} - \frac{0F_1\left(1, \frac{4\rho}{\sigma^4}\right)}{0F_1\left(2, \frac{4\rho}{\sigma^4}\right)} \right] \end{aligned} \qquad (A4b)$$

$$q_{1,t}(\rho) = b_{1,t}(\rho)/2. \qquad (A4c)$$

In the limit where $\sigma^2 \to 0$, those expressions tend to the following results, which can be further simplified for small density contrasts $\delta = \rho - 1$ to reproduce the Gaussian result

$$b_{1,t}(\rho) \xrightarrow{\sigma \to 0} \frac{2(\sqrt{\rho} - 1)}{\sigma^2} \xrightarrow{\delta \to 0} \frac{\delta}{\sigma^2} = b_{1,G}(\delta), \qquad (A5a)$$

$$b_{2,t}(\rho) \xrightarrow{\sigma \to 0} \frac{[2(\sqrt{\rho} - 1)]^2 - \sigma^2}{\sigma^4} \xrightarrow{\delta \to 0} \frac{\delta^2 - \sigma^2}{\sigma^4} = b_{2,G}(\delta). \qquad (A5b)$$

The two stages of this limiting behaviour are illustrated in Figure 24 by a sequence of decreasing $\sigma^2$ (from red to green) showing fast approach to the small variance limit (black solid) which reproduces the Gaussian result (black dashed) for small density contrast. The quadratic $q_1$ term in the covariance (56) can be combined with the $b_1$ term (and hence be captured by the first eigenvector), while the $q_1$ term adding to $b_{2,t}$ in (56) does not change much in the case of a small variance $\sigma^2 \ll 1$ and an intermediate density regime.

For the minimal tree model, one can deduce that the covariance at any order is given in terms of the set of bias functions $b_{n,t}(\rho)$. A recursion relation can be obtained as shown in Bernardeau (2022), which can be converted to an explicit closed-form expression as we show here.

*Derivation of closed form bias functions.* Define an auxiliary version of the Laplace transform of the CGF as the integral $I_j(\rho)$:

$$I_j(\rho) = \int \frac{d\lambda}{2\pi i} \exp[-\lambda\rho + j\varphi(\lambda)]. \qquad (A6)$$

The PDF is then the value $I_1(\rho)$. In the minimal tree model, the bias functions are predicted via

$$b_{n,t}(\rho)\mathcal{P}_t(\rho) = \int \frac{d\lambda}{2\pi i} \varphi_{0,t}(\lambda)^n \exp[-\lambda\rho + \varphi_{0,t}(\lambda)] \qquad (A7)$$

where the one-point CGF $\varphi_{0,t}(\lambda)$ is given by equation (A1). The nth bias function equation can also be written as

$$b_{n,t}(\rho)\mathcal{P}_t(\rho) = \frac{\mathrm{d}^n I_j(\rho)}{\mathrm{d}j^n}\bigg|_{j=1}. \qquad (A8)$$

and for the minimal tree model the $I_j(\rho)$ can be obtained as

$$I_j(\rho) = \frac{4j}{\sigma^4} \exp \left[ -\frac{2}{\sigma^2}(j + \rho) \right] {}_0F_1 \left( 2, \frac{4j\rho}{\sigma^4} \right). \qquad (A9)$$

We can therefore obtain a closed form equation for $b_n(\rho)$ by differentiation of equation A9. The nth derivative of

<figure><img src="image_30.png" /><figcaption>FIG. 24.— A comparison of the first two leading order bias functions from the tree model (A4), $b_{1,t}$ (upper panel) and $b_{2,t}$ (lower panel), for different variances (coloured lines) approaching the small-variance limit (A5) (green) and the Gaussian expectation (dashed line). In the upper panel we also show the measurements for $b_1$ from the Quijote simulations (data points) and the large-deviation theory prediction (black line) as in Figure 23, but in linear rather than logarithmic scale.</figcaption></figure>

$I_j$ is given by

$$\frac{\mathrm{d}^n I_j}{\mathrm{d}j^n} = \sum_{k=0}^{n} \sum_{\ell=0}^{k} \binom{n}{k} \binom{k}{\ell} f^{(n-k)}(j) g^{(k-\ell)}(j) h^{(\ell)}(j), \qquad (A10)$$

where we've split $I_j$ into three functions, f, g, h, which we can differentiate separately as follows

$$f(j) = \frac{4j}{\sigma^4}, \qquad g(j) = \exp \left[ -\frac{2}{\sigma^2}(j + \rho) \right], \qquad (A11a)$$

$$\frac{\mathrm{d}^n f}{\mathrm{d}j^n} = \frac{4}{\sigma^4}(j\delta_{n,0} + \delta_{n,1}), \qquad (A11b)$$

$$\frac{\mathrm{d}^n g}{\mathrm{d}j^n} = \left(\frac{-2}{\sigma^2}\right)^n \exp\left[-\frac{2}{\sigma^2}(j + \rho)\right] \qquad (A11c)$$

For the hypergeometric function we make use of the fact that $d/dz {}_0F_1(n, z) = {}_0F_1(n+1, z)/n$, where we defined $z = 4\rho j/\sigma^4$ so that $d/dj = \frac{4\rho}{\sigma^4}d/dz$ and

$$h(j) = {}_0F_1\left(2, \frac{4\rho j}{\sigma^4}\right), \qquad (A12a)$$

$$\frac{\mathrm{d}^n h}{\mathrm{d}j^n} = \left(\frac{4\rho}{\sigma^4}\right)^n \frac{1}{(n+1)!} {}_0F_1\left(2+n, \frac{4\rho j}{\sigma^4}\right). \qquad (A12b)$$
Putting this into the full sum we obtain (using the Kroenecker $\delta$s to collapse the $k$ sum, as the only surviving terms are when $n=k$ and when $n-k=1$).

$$ \frac{d^n I_j}{d j^n} = \frac{4}{\sigma^4} e^{-\frac{2}{\sigma^2}(j+\rho)} \sum_{k=0}^{n} \sum_{\ell=0}^{k} \binom{n}{k} \binom{k}{\ell} (j\delta_{n-k,0} + \delta_{n-k,1}) \\ \times \left(\frac{-2}{\sigma^2}\right)^{k-\ell} \left(\frac{4\rho}{\sigma^4}\right)^k \frac{0F_1\left(k+2, \frac{4j\rho}{\sigma^4}\right)}{(k+1)!} . \quad (A13) $$

Now, we separate this sum into the case where $k=n$ and when $k=n-1$, factor out common terms only depending on $n$ and use the binomial theorem to do these sums as simply $(1-2/\sigma^2)^n$ and $(1-2/\sigma^2)^{n-1}$ respectively. This equation is then the form of $b_{n,t}(\rho)\mathcal{P}_t(\rho)$ when $j=1$. Dividing through by $\mathcal{P}_t(\rho)$ from equation (A2) everywhere and setting $j=1$ gives us the functional form of the $n$-th order bias function for the minimal tree model as

$$ b_{n,t}(\rho) = \frac{1}{n!} \left[ \left( \frac{4\rho}{\sigma^4} \right) \left( 1 - \frac{2}{\sigma^2} \right) \right]^{n-1} \left[ \frac{4\rho \left( 1 - \frac{2}{\sigma^2} \right) \rho_1 \left( n + 2, \frac{4\rho}{\sigma^4} \right)}{\sigma^4 (n+1) \rho_1 \left( 2, \frac{4\rho}{\sigma^4} \right)} + \frac{\rho_1 \left( n + 1, \frac{4\rho}{\sigma^4} \right)}{\sigma^4 \rho_1 \left( 2, \frac{4\rho}{\sigma^4} \right)} \right]. \quad (A14) $$

Note that, the hypergeometrics can always be reexpanded down in terms of $\rho_1$ (1, $z$) and $\rho_1$ (2, $z$) using the identity

$$ \rho_1 \left( n-1, z \right) - \rho_1 \left( n, z \right) = \frac{z}{n(n-1)} \rho_1 \left( n+1, z \right). \quad (A15) $$

## B. REFINED DERIVATION OF COVARIANCE FOR CELLS IN LARGE-SEPARATION LIMIT

In Appendix C of Codis et al. (2016b) the result of the covariance of finding $N_i$ cells with density $\rho_i$ and $N_j$ spheres of density $\rho_j$ is given as (with no summation over repeated indices implied)

$$ \operatorname{cov}(N_i, N_j) = \bar{N}_i \bar{N}_j \bar{\xi} b(\rho_i) b(\rho_j) + \delta_{ij} \bar{N}_i, \quad (B1) $$

where $b(\rho)$ is the cell bias and $\bar{\xi}$ is the mean correlation of spheres defined in equation (B10). We can rewrite the ratio of the number of cells in the bin to the total number in terms of the measured PDF value in the bin and associated density bin size $\Delta_i$ as $N_i \approx N_T \mathcal{P}(\rho_i)\Delta_i$ and similarly $\bar{N}_i \sim N_T \bar{\mathcal{P}}(\rho_i)\Delta_i$. Hence, we have for the covariance of the histogram probability $\mathcal{P}_i$:

$$ \operatorname{cov}(\mathcal{P}(\rho_i), \mathcal{P}(\rho_j)) = \bar{\mathcal{P}}(\rho_i)\bar{\mathcal{P}}(\rho_j)\bar{\xi}b(\rho_i)b(\rho_j) + \delta_{ij} \frac{\bar{\mathcal{P}}(\rho_i)}{\Delta_i N_T}, \quad (B2) $$

Note that the shot-noise term only acts on the diagonal. Let us first focus on the diagonal, where the prediction for the noise to signal ratio is then is

$$ \frac{\operatorname{var}(\mathcal{P}(\rho_i))}{\bar{\mathcal{P}}^2(\rho_i)} = \bar{\xi} b^2(\rho_i) + \frac{1}{\bar{\mathcal{P}}(\rho_i) \Delta_i N_T}, \quad (B3) $$

The correlation matrix can be obtained by dividing the covariance by $\sqrt{\operatorname{var}(\mathcal{P}(\rho_1))\operatorname{var}(\mathcal{P}(\rho_2))}$. In the limit of

infinite volume, where one can place $N_T \to \infty$ well-separated cells, the correlation matrix would read

$$ \operatorname{corr}(\mathcal{P}(\rho_i), \mathcal{P}(\rho_j)) \xrightarrow{N_t \to \infty} \operatorname{sgn}(b(\rho_i)b(\rho_j)) = \pm 1, \quad (B4) $$

which predicts a characteristic 4-tiled pattern.

### B.1. Derivation for autocorrelation of PDF bins

Let us assume that the joint PDF of the density in $N_T$ cells (for simplicity with regular spacing) at large separation$^8$ reads

$$ \mathcal{P}(\hat{\rho}_1, \dots, \hat{\rho}_{N_T}; \{r_{IJ}\}) = \prod_{I=1}^{N_T} \mathcal{P}(\hat{\rho}_I) \left[ 1 + \sum_{I<j} b(\hat{\rho}_I)b(\hat{\rho}_J)\xi(r_{IJ}) \right]. \quad (B5) $$

To derive the result for the biased distribution of finding exactly $N$ cells in a given density bin $\mathcal{P}^b(N)$, we further expand the product in (B5) assuming $b^2\xi \ll 1$ as follows

$$ \mathcal{P}(\hat{\rho}_1, \dots, \hat{\rho}_{N_T}; \{r_{IJ}\}) = \left[ \prod_{I=1}^{N_T} \mathcal{P}(\hat{\rho}_I) \right] \left[ 1 + \sum_{I=1}^{N_T} \sum_{J<i} b(\hat{\rho}_I)b(\hat{\rho}_J)\xi(r_{IJ}) \right]. \quad (B6) $$

The probability to find exactly $N$ cells in a density bin $\hat{\Delta} = [\hat{\rho} - \Delta\rho/2, \hat{\rho} + \Delta\rho/2]$ is given by a) integrating over the probabilities of finding $N$ of the values of $\rho_1, \dots, \rho_{N_T}$ within $\hat{\Delta}$ and $N_T - N$ outside $\neg\hat{\Delta} = \mathbb{R}^+ \setminus \hat{\Delta}$, and b) summing over all different possibilities for the spatial distribution of the $N$ cells. Let us assume that the value $\hat{\rho}_I$ is fixed to a grid point $I$, such that the distance between the densities $\hat{\rho}_I$ and $\hat{\rho}_J$ are given by the distance of grid point $J$ to grid point $I$. Hence, we will drop $\{r_{IJ}\}$ from the argument of the joint PDF in the following. Then considering different spatial distributions of the cells corresponds to selecting a subset $S$ of length $N$ from the set of indices $\{1, \dots, N_T\}$

$$ \mathcal{P}^b(N) = \sum_S \prod_{I \in S} \left( \int_{\hat{\Delta}} d\hat{\rho}_I \right) \prod_{J \notin S} \left( \int_{-\hat{\Delta}} d\hat{\rho}_J \right) \mathcal{P}(\hat{\rho}_1, \dots, \hat{\rho}_{N_T}). \quad (B7) $$

Now we can insert equation (B6) and simplify using

$$ p = \int_{\hat{\Delta}} d\hat{\rho} \mathcal{P}(\hat{\rho}), \quad 1-p = \int_{-\hat{\Delta}} d\hat{\rho} \mathcal{P}(\hat{\rho}) \quad (B8) $$

$$ pb = \int_{\hat{\Delta}} d\hat{\rho} \mathcal{P}(\hat{\rho}) b(\hat{\rho}), \quad -pb = \int_{-\hat{\Delta}} d\hat{\rho} \mathcal{P}(\hat{\rho}) b(\hat{\rho}), \quad (B9) $$

<footnote>$^8$ A generalisation can be obtained by restoring the density-dependent cell two-point correlation $\xi_R(\hat{\rho}_I, \hat{\rho}_J; r_{IJ})$ instead of the factorised version that is valid at large separations.</footnote>
where the normalisation of $\mathcal{P}$ enforces $\int_{\mathbb{R}^+} d\hat{\rho} \mathcal{P}(\hat{\rho})b(\hat{\rho}) = 0$ and we will further use $\bar{b} := pb/p$. This leads to

$$ \begin{aligned} \mathcal{P}^b(N) = & \sum_S \left( p^N (1-p)^{N_T-N} \right. \\ & \qquad \left. + \sum_{S \ni I < J \in S} \xi(r_{IJ}) (pb)^2 p^{N-2} (1-p)^{N_T-N} \right. \\ & \qquad \left. + 2 \sum_{S \ni I < J \notin S} \xi(r_{IJ}) p^{N-1} pb(-pb)(1-p)^{N_T-N-1} \right. \\ & \qquad \left. + \sum_{S \ni I < J \notin S} \xi(r_{IJ}) p^N (-pb)^2 (1-p)^{N_T-N-2} \right) \\ = & p^N (1-p)^{N_T-N} \times \left\{ \begin{pmatrix} N_T \\ N \end{pmatrix} \right. \\ & \qquad \left. + \bar{b}^2 \sum_S \left( \sum_{S \ni I < J \in S} \xi(r_{IJ}) \right. \right. \\ & \qquad \left. \left. + \frac{\bar{N}^2}{(N_T - \bar{N})^2} \sum_{S \ni I < J \notin S} \xi(r_{IJ}) \right. \right. \\ & \qquad \left. \left. - \frac{2\bar{N}}{N_T - \bar{N}} \sum_{S \ni I < J \notin S} \xi(r_{IJ}) \right) \right\}, \end{aligned} $$

where in the second step we have converted $p = \bar{N}/N_T$.

Note that one can rewrite the sums over all $\bar{N}$-component subsets $S$ (of which there are $\binom{N_T}{N}$) in terms of the mean spatial correlation (which agrees with the definition used in Codis et al. (2016b))

$$ \sum_{I < J = 1}^{N_T} \xi(r_{IJ}), \quad \bar{\xi} = \frac{\sum_{I < J = 1}^{N_T} \xi(r_{IJ})}{N_T(N_T - 1)/2} \quad (B10) $$

as follows

$$ \begin{aligned} \sum_S \sum_{S \ni I < J \in S} \xi(r_{IJ}) &= \binom{N_T - 2}{N - 2} \sum_{I} \xi \\ &= \binom{N_T}{N} \frac{N(N-1)}{2} \bar{\xi}, \end{aligned} \quad (B11) $$

$$ \begin{aligned} 2 \sum_S \sum_{S \ni I < J \notin S} \xi(r_{IJ}) &= 2 \binom{N_T - 2}{N - 1} \sum_{I} \xi \\ &= \binom{N_T}{N} N(N_T - N) \bar{\xi}, \end{aligned} \quad (B12) $$

$$ \begin{aligned} \sum_S \sum_{S \ni I < J \notin S} \xi(r_{IJ}) &= \binom{N_T - 2}{N} \sum_{I} \xi \\ &= \binom{N_T}{N} \frac{(N_T - N)(N_T - N - 1)}{2} \bar{\xi}. \end{aligned} \quad (B13) $$

Where the full array of $r_{IJ}$ with $I < J$ is split into three subsets of separations a) among cells in the density bin (length $N(N-1)/2$), b) between cells in the density bin and outside (length $N(N_T-N)$), and c) among cells outside the density bin (length $(N_T - N)(N_T - N - 1)/2$). The combinatorial prefactors come from fixing two indices which corresponds to selecting from a set of $N_T - 2$

remaining indices, $N-2$ when both indices are in $S$, $N-1$ when one index is in $S$ and the other in $\mathcal{I} \setminus S$, and $N$ when both indices are in $\mathcal{I} \setminus S$. With this, we obtain

$$ \mathcal{P}^b(N) = p^N (1-p)^{N_T-N} \binom{N_T}{N} \times \left[ 1 + \bar{b}^2 \bar{\xi} \left( \frac{N(N-1)}{2} - \frac{\bar{N}N(N_T - N)}{N_T - \bar{N}} + \frac{\bar{N}^2(N_T - N)(N_T - N - 1)}{2(N_T - \bar{N})^2} \right) \right], $$

In the Poisson limit this gives<fnref n="9" />

$$ \mathcal{P}^b(N) = P^{\text{Poisss}}(N) \left[ 1 + \bar{b}^2 \bar{\xi} \left( \frac{N(N-1)}{2} - \bar{N}N + \frac{\bar{N}^2}{2} \right) \right], \quad (B14) $$

where $P^{\text{Poisss}}(N) = \bar{N}^N \exp(-\bar{N})/N!$. In Codis et al. (2016a) it was assumed that correlations are the same at all separations, $\xi(r_{IJ}) = \bar{\xi}$ with $\bar{\xi}$ defined in equation (B10). In this case, one can replace all sums, the sum over $r_{IJ}$ by a combinatoric factor for choosing $N$ cells out of $N_T$, and the sums over $I$ and $J$ by the number of pairs. As we showed here, it turns out that this result seems holds true even for the more general case of correlations that depend on separation. From this expression we can verify that the biased PDF is normalised and has mean $\langle N \rangle = \bar{N}$ and variance

$$ \text{var}(N) = \langle N^2 \rangle - \langle N \rangle^2 = \bar{N} + \bar{N}^2 \bar{\xi} \bar{b}^2 \quad (B15) $$

In the continuous limit, one can write the mean correlation as $\bar{\xi} = \int_{r_{\min}}^{r_{\max}} dr P_d(r)\xi(r)$, where $P_d(r)$ is the probability distribution of cell distances $r$. The maximum distance is set by the size of the survey patch/volume, for a $d$-dim box of side length $L$, we have $r_{\max} \approx \sqrt{dL}$. The minimum distance is set by the cell separation under consideration, so $r_{\min} = 2R$ for non-overlapping cells and $r_{\min} \to 0$ for heavily overlapping cells. Note that this is not the same as attempting to measure the average correlation in a periodic simulation box (with fixed particle number) via the following definition

$$ \tilde{\xi} = \frac{\sum_{i \neq j=1}^{N_T} \delta_i \delta_j}{N_T(N_T - 1)} = \frac{\sum_{i=1}^{N} \delta_i \left( \sum_{j=1}^{N} \delta_j - \delta_i \right)}{N_T(N_T - 1)} \quad (B16a) $$

$$ \langle \delta \rangle = 0 \quad \text{and} \quad \bar{\xi} = \frac{\sum_{i=1}^{N} \delta_i^2}{N_T(N_T - 1)} < 0, \quad (B16b) $$

which tends to be small but negative, because by construction the mean density contrast in a box vanishes.

## B.2. Derivation of cross-correlation between bins

Considering different spatial distributions of the cells where a given number are in two density bins corresponds to first selecting a subset of length $N_1$, $S_{N_1}$ from the set of indices $\mathcal{I} = \{1, \dots, N_T\}$ and then selecting another

<footnote n="9">⁹ Note that without the Poisson limit, the variance of the binomial distribution becomes $\bar{N}(1-p)$.</footnote>
subset of length $N_2$, $S_{N_2}$ out of the remaining indices $I \setminus S_{N_1}$. For brevity, let us denote $S_{N_{12}} = S_{N_1} \cup S_{N_2}$

$$ \mathcal{P}^b(N_1, N_2) = \sum_{S_{N_1}, S_{N_2}} \prod_{I \in S_{N_1}} \left( \int_{\hat{\Delta}_1} d\hat{\rho}_I \right) \prod_{J \in S_{N_2}} \left( \int_{\hat{\Delta}_2} d\hat{\rho}_J \right) \prod_{K \notin S_{N_{12}}} \left( \int_{-\hat{\Delta}} d\hat{\rho}_K \right) \mathcal{P}(\hat{\rho}_1, \dots, \hat{\rho}_{N_T}). \quad (B17) $$

As before, the integrals can be simplified

$$ \begin{aligned} & \mathcal{P}^b(N_1, N_2) \\ &= p_1^{N_1} p_2^{N_2} (1-p)^{\Delta N_{T,12}} \times \left\{ \binom{N_T}{N_1} \binom{N_T - N_1}{N_2} \right. \\ & \qquad + \sum_{S_{N_1}, S_{N_2}} \left( \bar{b}_1^2 \sum_{S_{N_1} \ni I < J \in S_{N_1}} \xi(r_{IJ}) + \bar{b}_2^2 \sum_{S_{N_2} \ni I < J \in S_{N_2}} \xi(r_{IJ}) \right. \\ & \qquad \left. + 2\bar{b}_1\bar{b}_2 \sum_{S_{N_1} \ni I < J \in S_{N_2}} \xi(r_{IJ}) \right. \\ & \qquad \left. - 2\frac{p_1\bar{b}_1 + p_2\bar{b}_2}{1-p_1-p_2} \left[ \bar{b}_1 \sum_{S_{N_1} \ni I < J \notin S_{N_{12}}} \xi(r_{IJ}) \right. \right. \\ & \qquad \left. \left. + \bar{b}_2 \sum_{S_{N_2} \ni I < J \notin S_{N_{12}}} \xi(r_{IJ}) \right] \right. \\ & \qquad \left. + \left( \frac{p_1\bar{b}_1 + p_2\bar{b}_2}{1-p_1-p_2} \right)^2 \sum_{S_{N_{12}} \ni I < J \notin S_{N_{12}}} \xi(r_{IJ}) \right\}, \end{aligned} $$

The sums over the different subsets $S_{N_1}$ and $S_{N_2}$ can be simplified using combinatorical properties. We use the average correlation $\bar{\xi}$ and its numerator $\sum \xi$ as defined in equation (B10), $\neg i$ to indicate the opposite index $i = 2$ if $i = 1$ and $\neg i = 1$ if $i = 2$, respectively. Introducing the short hand notation $N_{T,12} := N_T - N_1 - N_2$, we obtain

$$ \begin{aligned} & \sum_{S_{N_1}, S_{N_2}} \sum_{S_{N_i} \ni I < J \in S_{N_i}} \xi(r_{IJ}) \\ &= \binom{N_T - 2}{N_i - 2} \binom{N_T - N_i}{N_T - N_{-i}} \sum \xi \\ &= \binom{N_T}{N_1} \binom{N_T - N_1}{N_2} \bar{\xi} \frac{N_i(N_i - 1)}{2}, \quad (B18) \\ & 2 \sum_{S_{N_1}, S_{N_2}} \sum_{S_{N_1} \ni I < J \in S_{N_2}} \xi(r_{IJ}) \\ &= 2 \binom{N_T - 2}{N_1 - 1} \binom{N_T - N_1 - 1}{N_2 - 1} \sum \xi \\ &= \binom{N_T}{N_1} \binom{N_T - N_1}{N_2} \bar{\xi} N_1 N_2, \quad (B19) \end{aligned} $$

$$ \begin{aligned} & 2 \sum_{S_{N_1}, S_{N_2}} \sum_{S_{N_i} \ni I < J \notin S_{N_{12}}} \xi(r_{IJ}) \\ &= 2 \binom{N_T - 2}{N_i - 1} \binom{N_T - N_i - 1}{N_{-i}} \sum \xi \\ &= \binom{N_T}{N_1} \binom{N_T - N_1}{N_2} \bar{\xi} N_i \Delta N_{T,12}, \quad (B20) \\ & \sum_{S_{N_1}, S_{N_2}} \sum_{S_{N_{12}} \ni I < J \notin S_{N_{12}}} \xi(r_{IJ}) \\ &= \binom{N_T - 2}{N_1} \binom{N_T - N_1 - 2}{N_2} \sum \xi \\ &= \binom{N_T}{N_1} \binom{N_T - N_1}{N_2} \bar{\xi} \frac{\Delta N_{T,12}(\Delta N_{T,12} - 1)}{2}. \quad (B21) \end{aligned} $$

Additionally, one can rewrite

$$ \frac{p_1 \bar{b}_1 + p_2 \bar{b}_2}{1 - p_1 - p_2} = \frac{\bar{N}_1 \bar{b}_1 + \bar{N}_2 \bar{b}_2}{\bar{N}_T - \bar{N}_1 - \bar{N}_2}. \quad (B22) $$

With this one obtains

$$ \begin{aligned} & \mathcal{P}^b(N_1, N_2) \\ &= p_1^{N_1} p_2^{N_2} (1-p)^{\Delta N_{T,12}} \binom{N_T}{N_1} \binom{N_T - N_1}{N_2} \times \left\{ 1 + \right. \\ & \qquad \left. + \bar{\xi} \left[ \frac{N_1(N_1-1)}{2} \bar{b}_1^2 + \frac{N_2(N_2-1)}{2} \bar{b}_2^2 + N_1 N_2 \bar{b}_1 \bar{b}_2 \right. \right. \\ & \qquad \left. - (N_1 \bar{b}_1 + N_2 \bar{b}_2) \Delta N_{T,12} \frac{\bar{N}_1 \bar{b}_1 + \bar{N}_2 \bar{b}_2}{N_T - \bar{N}_1 - \bar{N}_2} \right. \\ & \qquad \left. + \frac{\Delta N_{T,12}(\Delta N_{T,12} - 1)}{2} \left( \frac{\bar{N}_1 \bar{b}_1 + \bar{N}_2 \bar{b}_2}{N_T - \bar{N}_1 - \bar{N}_2} \right)^2 \right] \right\}. \end{aligned} $$

In the Poisson limit this becomes $^{10}$

$$ \mathcal{P}^b(N_1, N_2) = P^{\text{Poisss}}(N_1) P^{\text{Poisss}}(N_2) \times \left[ 1 + \right. \\ \left. + \bar{\xi} \bar{b}_1^2 \left( \frac{N_1(N_1-1)}{2} + \frac{\bar{N}_1^2}{2} - N_1 \bar{N}_1 \right) \right. \\ \left. + \bar{\xi} \bar{b}_2^2 \left( \frac{N_2(N_2-1)}{2} + \frac{\bar{N}_2^2}{2} - N_2 \bar{N}_2 \right) \right. \\ \left. + \bar{\xi} \bar{b}_1 \bar{b}_2 \left( N_1 N_2 + \bar{N}_1 \bar{N}_2 - N_1 \bar{N}_2 - \bar{N}_1 N_2 \right) \right] $$

With this we reproduce the result from Codis et al. (2016a) for the Poisson limit and the associated covariance

$$ \langle N_1 N_2 \rangle = \bar{N}_1 \bar{N}_2 (1 + \bar{\xi} \bar{b}_1 \bar{b}_2). \quad (B23) $$

<footnote>$^{10}$ Note that without the Poisson limit, the covariance of the multinomial distribution becomes $-N_T p_1 p_2 = -\bar{N}_1 \bar{N}_2 / N_T$.</footnote>