© ESO 2018

1. Introduction

Studying the matter distribution of the present large-scale structure reveals a wealth of information about the evolution of the Universe. In particular, its distorting effect on the propagation of light from distant galaxies, known as cosmic shear, can be captured by analysing weak lensing surveys. By comparing the results of cosmological models with the observed signal, one can constrain cosmological parameters (see e.g. Asgari et al. 2021; Abbott et al. 2022; Hamana et al. 2020).

The preferred methods used to infer statistical properties of the matter and galaxy distribution concentrate on second-order statistics, such as the two-point correlation functions or their Fourier counterparts, the power spectra. Although these statistics have an impressive accuracy to describe for instance primordial perturbations visible in the cosmic microwave background (CMB; e.g. Planck Collaboration V 2020) they probe only the Gaussian information present in the density fluctuations. However, these initial conditions developed significant non-Gaussian features by means of non-linear gravitational instability, which can only be investigated with higher-order statistics. Although they are typically more time consuming to model and measure, these higher-order statistics scale differently with cosmological parameters, and are not affected in the same way by residual systematics. Hence, by jointly investigating second- and higher-order statistics, the constraining power on cosmological parameters increases (see e.g. Bergé et al. 2010; Pyne & Joachimi 2021; Pires et al. 2012; Fu et al. 2014; Kilbinger & Schneider 2005).

A large number of analytical models for the two-point statistics exists in the literature (Takahashi et al. 2012; Heitmann et al. 2014; Euclid Collaboration 2021; Mead et al. 2020; Nishimichi et al. 2021); however, the analysis of higher-order statistics is usually based on simulations. Analytical models for higher-order lensing statistics are rare, although they are important not only for scientists to understand physical processes, but also to cross-check simulations, which are usually only tested against Gaussian statistics. For example, Reimberg & Bernardeau (2018) and Barthelemy et al. (2021) used large deviation theory (LDT) to compute the reduced-shear correction to the aperture mass probability distribution function (PDF); Munshi et al. (2020) and Halder et al. (2021) analytically modelled the integrated shear three-point function; the lensing peak count function was modelled in Fan et al. (2010), Lin & Kilbinger (2015) and Shan et al. (2018), while the lensing PDF is modelled in Boyle et al. (2021).

The examples mentioned above all pertain to the analysis of cosmic shear data. However, it has been established in recent analyses that the addition of foreground clustering data, and their cross-correlation with the background source galaxies, yield significantly better constraints (Abbott et al. 2018; Heymans et al. 2021). While the central analyses focused again on two-point statistics, Friedrich et al. (2018, hereafter F18) developed a competitive model based on density split statistics (hereafter DSS). The idea is to measure the mean tangential shear around small sub-areas of the survey, and to stack the signal according to the foreground galaxy density in these sub-areas. We expect the tangential shear to be larger around points with a high density of foreground galaxies, given they correspond to a large matter overdensity on average. The model derived in F18 is based on non-perturbative calculations of the matter density PDF, which predicts the shear profiles and the probability density of galaxy counts in the sub-areas, for a given cosmological model, a redshift distribution for the source and lens galaxy populations, and a mean galaxy density. In Gruen et al. (2018, hereafter G18), the F18 model is used to constrain cosmological parameters from DSS measurements from the Dark Energy Survey (DES) First Year and Sloan Digital Sky Survey (SDSS) data, which yields constraints on the matter density ${\mathrm{\Omega }}_{\text{m}}=0.{26}_{-0.03}^{+0.04}$ that agree and are competitive with the DES analysis of galaxy and shear two-point functions (see Abbott et al. 2018).

One of the motivations of this work is based on Burger et al. (2020, hereafter B20), who use a suite of numerical simulations to show that using matched filter functions for searching peaks and troughs in the galaxy and matter density contrast has clear advantages compared to the top-hat filter used in the F18 model, both in terms of the overall signal S/N and in recovering accurately the galaxy bias term. Another motivation of using compensated filters is that these filters are more confined in Fourier space and are therefore better at smoothing out large ℓ-modes where baryonic effects play an important role (Asgari et al. 2020). Therefore, it is of interest to generalise the DSS to a broader set of filter functions. Smoothing cosmic density fields with filters other than top-hat ones significantly complicates the LDT-like calculations used by F18 and G18 (cf. Barthelemy et al. 2021) because for top-hat filters the Lagrangian to Eulerian mapping inherent in LDT is particularly simple. However, we find here that density split statistics with non-top-hat filters that are sufficiently concentrated around their centres can still be accurately modelled with computationally feasible extensions of approximations made by F18. This paper describes our modifications to the F18 model that will allow us to optimise filtering strategies when applying density split statistics to Stage III weak lensing surveys such as KiDS. Throughout this paper, if not otherwise stated, we assume a spatially flat universe.

This work is structured as follows. In Sect. 2 we review the basics of the aperture statistics; we then detail our changes to the F18 model in Sect. 3. In Sect. 4 we describe the simulations, and the construction of our mock data used to validate the revised model. In Sect. 5 we compare the model predictions with simulations, and establish the model’s limitations. We summarise our work in Sect. 6.

2. Aperture statistics

The lensing convergence κ and shear γ are related via the lensing potential ψ (Schneider et al. 1992) as

$$\begin{array}{c}\hfill \kappa (\mathit{\theta })=\frac{1}{2}({\partial }_1^2+{\partial }_2^2)\psi (\mathit{\theta }),\gamma (\mathit{\theta })=\frac{1}{2}({\partial }_1^2-{\partial }_2^2+2\text{i}{\partial }_1{\partial }_2)\psi (\mathit{\theta }),\end{array}$$(1)

with ${\partial }_i=\frac{\partial }{\partial {\theta }_i}$ and θ the angular position on the sky; we employ the flat-sky approximation. Given a reference point in a Cartesian coordinate system on the sky and a second point whose separation to the first is oriented at an angle ϕ with respect to that coordinate system, we can express the shear at the second point in terms of the tangential and cross-shear with respect to the first point as

$$\begin{array}{c}\hfill {\gamma }_{\mathrm{t}}=-\text{Re}\left(\gamma {\text{e}}^{-2\text{i}\varphi }\right),{\gamma }_{\times }=-\text{Im}\left(\gamma {\text{e}}^{-2\text{i}\varphi }\right),\end{array}$$(2)

where the factor 2 in the exponent is due to the polar nature of the shear. Given a convergence field κ(θ), the aperture mass at position θ is defined as

$$\begin{array}{c}\hfill M_{\text{ap}}\left(\mathit{\theta }\right){\displaystyle \int }{\mathrm{d}}^2{\theta }^{\prime }\kappa (\mathit{\theta }+{\mathit{\theta }}^{\prime })U\left(\left|{\mathit{\theta }}^{\mathbf{\prime }}\right|\right),\end{array}$$(3)

where U(ϑ) is a compensated axisymmetric filter function, such that ∫ϑ U(ϑ) dϑ = 0. As shown in Schneider (1996), if U is compensated, M_ap can also be expressed in terms of the tangential shear γ_t and a related filter function Q as

$$\begin{array}{c}\hfill M_{\text{ap}}(\mathit{\theta })={\displaystyle \int }{\mathrm{d}}^2{\theta }^{\prime }{\gamma }_{\text{t}}(\mathit{\theta }+{\mathit{\theta }}^{\prime })Q\left(\left|{\mathit{\theta }}^{\prime }\right|\right),\end{array}$$(4)

where

$$\begin{array}{c}\hfill Q(\vartheta )=\frac{2}{{\vartheta }^2}{\displaystyle \underset{0}{\overset{\vartheta }{\int }}\mathrm{d}{\vartheta }^{\prime }{\vartheta }^{\prime }U({\vartheta }^{\prime })-U(\vartheta ),}\end{array}$$(5)

which can be inverted, yielding

$$\begin{array}{c}\hfill U(\vartheta )=2{\displaystyle \underset{\vartheta }{\overset{\infty }{\int }}\mathrm{d}{\vartheta }^{\prime }\frac{Q({\vartheta }^{\prime })}{{\vartheta }^{\prime }}-Q(\vartheta ).}\end{array}$$(6)

In analogy to M_ap, we define, as done in B20, the aperture number counts (Schneider 1998), or aperture number, as

$$\begin{array}{c}\hfill N_{\text{ap}}(\mathit{\theta }){\displaystyle \int }{\mathrm{d}}^2{\theta }^{\prime }n(\mathit{\theta }+{\mathit{\theta }}^{\prime })U\left(\left|{\mathit{\theta }}^{\prime }\right|\right),\end{array}$$(7)

where U(ϑ) is the same filter function as in Eq. (3) and n(ϑ) is the (foreground) galaxy number density on the sky. This definition of the aperture number is equivalent to the “Counts-in-Cell” (CiC) from Gruen et al. (2016) if a top-hat filter of the form

$$\begin{array}{c}\hfill U_{\text{th}}(\vartheta )=\frac{1}{A}\mathcal{H}({\vartheta }_{\text{th}}-\vartheta ),\end{array}$$(8)

is used, where ℋ is the Heaviside step function and A is the area of the filter. However, B20 demonstrated that top-hat filters are not optimal, and that a better performance is achieved by an adapted filter in terms of signal-to-noise-ratio (S/N) and in recovering accurately the galaxy bias term. In this paper we compute aperture mass statistics with Eq. (4) using simulated weak lensing catalogues of background source galaxies, notably regarding positions and ellipticities, and aperture number statistics with Eq. (7) from the position of simulated foreground lens galaxies (see Sect. 4).

3. Revised model

In this section we describe our modifications of the original F18 model. Although the derivations shown here are self-contained, we recommend the interested reader to consult the original F18 paper. In particular, it is shown there that the full non-perturbative calculation of the PDF within large deviation theory (LDT) can be well approximated with a log-normal PDF that matches variance and skewness of the LDT result. This allowed F18 and G18 to replace the expensive LDT computation with a faster one, hence making the calculation of full Markov chain Monte Carlo (MCMC) functions feasible. The reason why this approximation works well is that, for top-hat filters, the scaling of variance and higher-order cumulants in LDT is similar to that found in log-normal distributions. This cannot be expected a priori for other filter functions. However, through comparison with N-body simulations we find here (cf. Sect. 5) that either a simple log-normal or a combination of two log-normal distributions still accurately describes the density PDFs required to analyse density split statistics with more general classes of filters. The following section describes these calculations. In order to reduce the mathematical calculations in this section, some derivations are detailed in Appendix A.

We start by defining the line-of-sight projection of the 3D matter density contrast δ_m,3D, weighted by a foreground (lens) galaxy redshift probability distribution n_f(z) as

$$\begin{array}{c}\hfill {\delta }_{\text{m,2D}}(\mathit{\theta })={\displaystyle \int }\mathrm{d}\chi q_{\text{f}}(\chi ){\delta }_{\text{m,3D}}(\chi \mathit{\theta },\chi ),\end{array}$$(9)

where χ is the co-moving distance and the projection kernel q_f(χ) is

$$\begin{array}{c}\hfill q_{\text{f}}(\chi )=n_{\text{f}}\left(z\left[\chi \right]\right)\frac{\mathrm{d}z[\chi ]}{\mathrm{d}\chi }·\end{array}$$(10)

This 2D matter density contrast can then be used together with a linear bias term to represent a tracer density contrast (see Sect. 3.3 or Sect. 4). Following F18, the next step consists of smoothing the results with a filter U of size Θ:

$$\begin{array}{cc}\hfill {\delta }_{\text{m},U}(\mathit{\theta })& \equiv {\displaystyle \underset{|{\mathit{\theta }}^{\prime }|<\mathrm{\Theta }}{\int }{\mathrm{d}}^2{\theta }^{\prime }{\delta }_{\text{m,2D}}(\mathit{\theta }+{\mathit{\theta }}^{\prime })U\left(\left|{\mathit{\theta }}^{\prime }\right|\right).}\hfill \end{array}$$(11)

This simplifies in the case of a top-hat filter of size Θ to

$$\begin{array}{c}\hfill {\delta }_{\text{m},\text{th}}^{\mathrm{\Theta }}(\mathit{\theta })=\frac{1}{A}{\displaystyle \underset{\left|{\mathit{\theta }}^{\prime }\right|<\mathrm{\Theta }}{\int }{\mathrm{d}}^2{\theta }^{\prime }{\delta }_{\text{m,2D}}(\mathit{\theta }+{\mathit{\theta }}^{\prime }).}\end{array}$$(12)

Similar to the 2D density contrast, the convergence, which is needed to describe the DSS signal, is given by

$$\begin{array}{c}\hfill \kappa (\mathit{\theta })={\displaystyle \int }\mathrm{d}\chi W_{\text{s}}(\chi ){\delta }_{\text{m,3D}}(\chi \mathit{\theta },\chi ),\end{array}$$(13)

where W_s(χ) is the lensing efficiency defined as

$$\begin{array}{c}\hfill W_{\text{s}}(\chi )=\frac{3{\mathrm{\Omega }}_{\text{m}}{H}_0^2}{2c^2}{\displaystyle {\int }_{\chi }^{\infty }}\mathrm{d}{\chi }^{\prime }\frac{\chi ({\chi }^{\prime }-\chi )}{{\chi }^{\prime }a(\chi )}{q}_{\text{s}}\left({\chi }^{\prime }\right),\end{array}$$(14)

with $q_{\text{s}}(\chi )=n_{\text{s}}\left(z\left[\chi \right]\right)\frac{\text{d}z[\chi ]}{\text{d}\chi }$ being the line-of-sight probability density of the sources, Ω_m the matter density parameter, H₀ the Hubble parameter, and c the speed of light. The mean convergence inside an angular separation ϑ, κ_{< ϑ}, follows then in analogy to Eq. (12) by substituting δ_m,2D(θ) with κ(θ).

The aim of our model is to predict the tangential shear profiles γ_t given a quantile 𝒬 of the foreground aperture number N_ap, ⟨γ_t|𝒬⟩, where for instance the highest quantile is the set of lines of sight of the sky that have the highest values of N_ap. Therefore, to determine ⟨γ_t|𝒬⟩ the model calculates ⟨γ_t|N_ap⟩ and sums up all that belong to the corresponding quantile 𝒬. The expectation value of ⟨γ_t|N_ap⟩ is computed from the convergence profile as

$$\begin{array}{cc}\hfill \langle {\gamma }_{\text{t}}(\vartheta )|N_{\text{ap}}\rangle & =\langle {\kappa }_{<\vartheta }|N_{\text{ap}}\rangle -\langle {\kappa }_{\vartheta }|N_{\text{ap}}\rangle =-\frac{\vartheta }{2}\frac{\mathrm{d}}{\mathrm{d}\vartheta }\langle {\kappa }_{<\vartheta }|N_{\text{ap}}\rangle ,\hfill \end{array}$$(15)

where κ_ϑ is the azimuthally averaged convergence at angular separation ϑ from the centre of the filter, and κ_{< ϑ} is the average convergence inside that radius. The latter quantity, conditioned on a given N_ap, can be specified by

$$\begin{array}{cc}\hfill \langle {\kappa }_{<\vartheta }|N_{\text{ap}}\rangle & ={\displaystyle \int }\mathrm{d}{\delta }_{\text{m},U}\langle {\kappa }_{<\vartheta }|{\delta }_{\text{m},U},N_{\text{ap}}\rangle p\left({\delta }_{\text{m},U}|N_{\text{ap}}\right)\hfill \end{array}$$(16)

$$\begin{array}{cc}& \approx {\displaystyle \int }\mathrm{d}{\delta }_{\text{m},U}\langle {\kappa }_{<\vartheta }|{\delta }_{\text{m},U}\rangle p\left({\delta }_{\text{m},U}|N_{\text{ap}}\right),\hfill \end{array}$$(17)

where in the second step we assumed that the expected convergence within ϑ only depends on the projected matter density contrast δ_m, U and not on the particular realisation of shot-noise in N_ap found within that fixed matter density contrast¹.

By use of Bayes’ theorem, we can express the conditional PDF as

$$\begin{array}{c}\hfill p\left({\delta }_{\text{m},U}|N_{\text{ap}}\right)=\frac{p\left(N_{\text{ap}}|{\delta }_{\text{m},U}\right)p\left({\delta }_{\text{m},U}\right)}{p\left(N_{\text{ap}}\right)},\end{array}$$(18)

where p(N_ap|δ_m, U) is the probability of finding N_ap given the smoothed density contrast δ_m, U. The normalisation in the denominator of Eq. (18) follows by integrating over the numerator,

$$\begin{array}{c}\hfill p\left(N_{\text{ap}}\right)={\displaystyle \int }\mathrm{d}{\delta }_{\text{m},U}p\left({\delta }_{\text{m},U}\right)p\left(N_{\text{ap}}|{\delta }_{\text{m},U}\right).\end{array}$$(19)

As seen in the derivation above, we are left with three ingredients in order to calculate the tangential shear profiles given a quantile Q of the aperture number ⟨γ_t(ϑ)|N_ap⟩:

1. the PDF of the matter density contrast smoothed with the filter function U (used in Eqs. (18), (19))

   $$\begin{array}{c}\hfill p\left({\delta }_{\text{m},U}\right);\end{array}$$(20)

2. the expectation value of the convergence inside a radius ϑ given the smoothed density contrast (used in Eq. (17))

   $$\begin{array}{c}\hfill \langle {\kappa }_{<\vartheta }|{\delta }_{\text{m},U}\rangle ;\end{array}$$(21)

3. the distribution of N_ap for the given filter function U given the smoothed density contrast (used in Eqs. (18), (19))

   $$\begin{array}{c}\hfill p\left(N_{\text{ap}}|{\delta }_{\text{m},U}\right).\end{array}$$(22)

Since all three ingredients are sensitive to the filter U, we need to adjust all of them coherently with respect to the top-hat case.

3.1. (I):p(δ_m,U)

As shown by F18 the full LDT computation of the matter density PDF can accurately approximated by a shifted log-normal distribution with vanishing mean (Hilbert et al. 2011), which is fully characterised by two parameters, σ and δ₀, as

$$\begin{array}{c}\hfill p\left({\delta }_{\text{m},U}\right)=\frac{1}{\sqrt{2\pi }\sigma ({\delta }_{\text{m},U}+{\delta }_0)}exp(-\frac{{[ln({\delta }_{\text{m},U}/{\delta }_0+1)+{\sigma }^2/2]}^2}{2{\sigma }^2})·\end{array}$$(23)

The two free parameters can be determined by specifying the variance $\langle {\delta }_{\text{m},U}^2\rangle$ and skewness $\langle {\delta }_{\text{m},U}^3\rangle$ of the PDF as (Hilbert et al. 2011)

$$\begin{array}{cc}\hfill \langle {\delta }_{\text{m},U}^2\rangle & ={\delta }_0^2[exp\left({\sigma }^2\right)-1],\hfill \end{array}$$(24)

$$\begin{array}{cc}\hfill \langle {\delta }_{\text{m},U}^3\rangle & =\frac{3}{{\delta }_0}{\langle {\delta }_{\text{m},U}^2\rangle }^2+\frac{1}{{\delta }_0^3}{\langle {\delta }_{\text{m},U}^2\rangle }^3;\hfill \end{array}$$(25)

we derive the expression of $\langle {\delta }_{\text{m},U}^2\rangle$ and $\langle {\delta }_{\text{m},U}^3\rangle$ in Appendix A (see Eq. (A.28)).

As we show later, this approximation works well for non-negative filter functions like top-hat or Gaussian filters. However, the log-normal PDF approximation becomes less accurate for compensated filters that include negative weights. In these cases we instead divide U into U_> (ϑ) = U(ϑ)ℋ(ϑ_ts − ϑ) and U_< (ϑ) = − U(ϑ)ℋ(ϑ − ϑ_ts), where ϑ_ts is the transition scale from positive to negative filter weights. As a consequence, we obtain two correlated log-normal random variables, δ_m, U> and δ_m, U<, whose joint distribution can be represented by a bi-variate log-normal distribution as

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill p({\delta }_{\text{m},U_{>}},{\delta }_{\text{m},U_{<}})& =\frac{1}{2\pi {\sigma }_{>}({\delta }_{\text{m},U_{>}}+{\delta }_{0,>}){\sigma }_{<}({\delta }_{\text{m},U_{<}}+{\delta }_{0,<})\sqrt{1-{\rho }^2}}\hfill \\ \hfill & \times exp(-\frac{1}{2(1-{\rho }^2)}[{\mathrm{\Delta }}_{>}^2+{\mathrm{\Delta }}_{<}^2-2\rho {\mathrm{\Delta }}_{>}{\mathrm{\Delta }}_{<}]),\hfill \end{array}\end{array}$$(26)

where we defined

$$\begin{array}{c}\hfill {\mathrm{\Delta }}_{>}=\frac{ln({\delta }_{\text{m},U_{>}}/{\delta }_{0,>}+1)+{\sigma }_{>}^2/2}{{\sigma }_{>}}\end{array}$$(27)

and similarly for Δ_<. The correlation coefficient ρ is determined by

$$\begin{array}{cc}\hfill \rho & =ln(\frac{\langle {\delta }_{\text{m},U_{>}}{\delta }_{\text{m},U_{<}}\rangle }{{\delta }_{0,>}{\delta }_{0,<}}+1)\frac{1}{{\sigma }_{>}{\sigma }_{<}},\hfill \end{array}$$(28)

and in order to calculate the difference of two independent random variables δ_m, U = δ_m, U> − δ_m, U< we can use the convolution theorem (Arfken & Weber 2008) to get

$$\begin{array}{c}\hfill p\left({\delta }_{\text{m},U}\right)={\displaystyle \underset{-\infty }{\overset{\infty }{\int }}\mathrm{d}{\delta }_{\text{m},U_{>}}p({\delta }_{\text{m},U_{>}},{\delta }_{\text{m},U_{>}}-{\delta }_{\text{m},U}).}\end{array}$$(29)

3.2. (II):⟨κ_{< ϑ}|δ_m, U⟩

In order to calculate the expectation value of the mean convergence inside an angular radius ϑ, κ_{< ϑ}, given the matter density contrast δ_m, U, we assume that both follow a joint log-normal distribution (see e.g. the discussion in Appendix B of G18). In this case, the expectation value can be written as

$$\begin{array}{c}\hfill \frac{\langle {\kappa }_{<\vartheta }|{\delta }_{\text{m},U}\rangle }{{\kappa }_0}=exp\left(\frac{C[2ln({\delta }_{\text{m},U}/{\delta }_0+1)+V-C]}{2V}\right)-1,\end{array}$$(30)

where δ₀ is determined with Eq. (25) and the three variables C, V, and κ₀ can be calculated from the moments $\langle {\delta }_{\text{m},U}^2\rangle$, $\langle {\kappa }_{<\vartheta }{\delta }_{\text{m},U}\rangle$, and $\langle {\kappa }_{<\vartheta }{\delta }_{\text{m},U}^2\rangle$, which follow from the derivation in Appendix A (see Eq. (A.29)):

$$\begin{array}{cc}\hfill V& =ln(1+\frac{\langle {\delta }_{\text{m},U}^2\rangle }{{\delta }_0^2}),\hfill \end{array}$$(31)

$$\begin{array}{cc}\hfill C& =ln(1+\frac{\langle {\kappa }_{<\vartheta }{\delta }_{\text{m},U}\rangle }{{\delta }_0{\kappa }_0}),\hfill \end{array}$$(32)

$$\begin{array}{cc}\hfill {\kappa }_0& =\frac{{\langle {\kappa }_{<\vartheta }{\delta }_{\text{m},U}\rangle }^2{\text{e}}^V}{\langle {\kappa }_{<\vartheta }{\delta }_{\text{m},U}^2\rangle -2\langle {\kappa }_{<\vartheta }{\delta }_{\text{m},U}\rangle \langle {\delta }_{\text{m},U}^2\rangle /{\delta }_0}·\hfill \end{array}$$(33)

We note that the assumption that δ_m, U is log-normal distributed is not well justified for filters with negative weights as we mentioned in the previous section. A possible improvement could be done for instance by assuming again that δ_m, U is made up of two log-normal random variables, and we would need to calculate conditional moments like ⟨κ_{< ϑ}|δ_m, U>−δ_m, U<⟩. This would significantly increase the amount of joint moments needed in our calculation and would render fast modelling unfeasible. However, an improved modelling is also unnecessary at present, given the statistical uncertainties we expect for Stage III weak lensing surveys such as KiDS-1000. We demonstrate this empirically in Sect. 5 by comparison to N-body simulated data, but we also want to give a brief theoretical motivation. The average value of κ_{< ϑ}, given that δ_m, U lies within the range [δ_min, δ_max], is given by

$$\begin{array}{c}\hfill \langle {\kappa }_{<\vartheta }|{\delta }_{\text{m},U}\in [{\delta }_{min},{\delta }_{max}]\rangle =\frac{{\int }_{{\delta }_{min}}^{{\delta }_{max}}\mathrm{d}{\delta }_{\text{m},U}p\left({\delta }_{\text{m},U}\right)\langle {\kappa }_{<\vartheta }|{\delta }_{\text{m},U}\rangle }{{\int }_{{\delta }_{min}}^{{\delta }_{max}}\mathrm{d}{\delta }_{\text{m},U}p\left({\delta }_{\text{m},U}\right)}·\end{array}$$(34)

If κ_{< ϑ} and δ_m, U were joint Gaussian random variables, then p(δ_m, U) would be a Gaussian PDF and we would have $\langle {\kappa }_{<\vartheta }|{\delta }_{\text{m},U}\rangle ={\delta }_{\text{m},U}\langle {\delta }_{\text{m},U}{\kappa }_{<\vartheta }\rangle /\langle {\delta }_{\text{m},U}^2\rangle$. We now argue that the leading-order correction to this Gaussian approximation consists of replacing p(δ_m, U) by our full non-Gaussian model, without changing ⟨κ_{< ϑ}|δ_m, U⟩, since this would be exactly correct in the limit of strong correlation between the two variables. Our log-normal approximation to ⟨κ_{< ϑ}|δ_m, U⟩ is then already a next-to-leading-order correction and a bi-variate log-normal approximation for ⟨κ_{< ϑ}|δ_m, U⟩ would be of even higher order. While this reasoning is admittedly only heuristic, it is proven correct by the accuracy of our model predictions for the lensing signals in Sect. 5.

3.3. (III):p(N_ap|δ_m, U)

The third basic ingredient is the PDF of N_ap given the projected matter density contrast smoothed with the filter U. Assuming a Poisson distribution for N_ap, which is the most straightforward ansatz, is unfortunately not possible because negative values are expected with a compensated filter (i.e. in some of the U_< contributions). We use instead a completely new approach compared to F18, and derive an expression for p(N_ap|δ_m, U) by use of the characteristic function (Papoulis & Pillai 1991, hereafter CF), which is an alternative representation of a probability distribution, similar to the moment generating functions, but based on the Fourier transform of the PDF. Of interest to us, the n th derivative of the CFs can be used to calculate the n th moment of the PDF. The CF corresponding to p(N_ap|δ_m, U) is defined as

$$\begin{array}{c}\hfill \mathrm{\Psi }(t)={\langle {\text{e}}^{i\text{t}{N}_{\text{ap}}}\rangle }_{{\delta }_{\text{m},U}}={\displaystyle {\int }_{\mathbb{R}}\mathrm{d}{N}_{\text{ap}}p\left(N_{\text{ap}}|{\delta }_{\text{m},U}\right){\text{e}}^{\text{i}t{N}_{\text{ap}}},}\end{array}$$(35)

where in our particular case, we derive in Appendix A.4 a closed expression as

$$\begin{array}{c}\hfill \mathrm{\Psi }(t)=exp\left(2\pi n_0{\displaystyle {\int }_0^{\infty }\mathrm{d}\vartheta \vartheta (1+b\langle w_{\vartheta }|{\delta }_{\text{m},U}\rangle )[{\text{e}}^{\text{i}tU(\vartheta )}-1]}\right)\end{array}$$(36)

with n₀ being the mean number density of foreground galaxies on the sky. The assumption of linear galaxy bias enters here by the term b ⟨w_ϑ|δ_m, U⟩, with

$$\begin{array}{c}\hfill w_{\vartheta }=\frac{1}{2\pi }{\displaystyle {\int }_0^{2\pi }\mathrm{d}\varphi {\delta }_{\mathrm{m},2\mathrm{D}}(\vartheta ,\varphi ).}\end{array}$$(37)

Hence, n₀(1+b ⟨w_ϑ|δ_m, U⟩) is the effective number density at ϑ given δ_m, U. The conditional expectation value ⟨w_ϑ|δ_m, U⟩ is given in analogy to Eq. (30), but replacing $\langle {\kappa }_{<\vartheta }{\delta }_{\text{m},U}^k\rangle \to \langle w_{<\vartheta }{\delta }_{\text{m},U}^k\rangle$ in Eqs. (31)–(33) for k = 1, 2 and using that

$$\begin{array}{c}\hfill \langle w_{\vartheta }|{\delta }_{\text{m},U}\rangle =\langle w_{<\vartheta }|{\delta }_{\text{m},U}\rangle +\frac{\vartheta }{2}\frac{\mathrm{d}}{\mathrm{d}\vartheta }\langle w_{<\vartheta }|{\delta }_{\text{m},U}\rangle ,\end{array}$$(38)

where the joint moments $\langle w_{<\vartheta }{\delta }_{\text{m},U}^k\rangle$ are also derived in Eq. (A.30). Next, we re-express Eq. (36) as the product of two terms,

$$\begin{array}{c}\hfill \mathrm{\Psi }(t)=exp[p(t)]exp[\text{i}q(t)],\end{array}$$(39)

where

$$\begin{array}{cc}\hfill p(t)& =2\pi n_0{\displaystyle {\int }_0^{R_{\max }}\mathrm{d}\vartheta \vartheta (1+b\langle w_{\vartheta }|{\delta }_{\text{m},U}\rangle )(cos[tU(\vartheta )]-1),}\hfill \end{array}$$(40)

$$\begin{array}{cc}\hfill q(t)& =2\pi n_0{\displaystyle {\int }_0^{R_{\max }}\mathrm{d}\vartheta \vartheta (1+b\langle w_{\vartheta }|{\delta }_{\text{m},U}\rangle )sin[tU(\vartheta )],}\hfill \end{array}$$(41)

and R_max is the angular radius beyond which U vanishes. We note that G18 and F18 found super-Poisson shot-noise in their work. They interpret these deviations from Poisson noise as having a number ≠1 of galaxies per Poisson halo. This would suggest that we could incorporate non-Poissonian behaviour by replacing n₀ with an effective density of Poisson halos and making this a free parameter of our model. However, more recent investigations (e.g. Friedrich et al., in prep.) cast doubt on the simplified interpretation of F18 and G18. A proper investigation of the problem of non-Poissonian shot-noise is beyond the scope of this work, and we will address it in future investigations.

Finally, the probability density function p(N_ap|δ_m, U) follows from the inverse Fourier transform of the CF

$$\begin{array}{cc}\hfill p\left(N_{\text{ap}}|{\delta }_{\text{m},U}\right)& =\frac{1}{2\pi }{\displaystyle {\int }_{\mathbb{R}}\text{d}texp(-\text{i}t{N}_{\text{ap}})\mathrm{\Psi }(t)}\hfill \\ \hfill & =\frac{1}{2\pi }{\displaystyle {\int }_{\mathbb{R}}\text{d}tcos(q(t)-t{N}_{\text{ap}})exp[p(t)],}\hfill \end{array}$$(42)

where the second step follows from the fact that the imaginary part cancels out.

In Appendix A.4 we discuss a similar approach, where we assume that p(N_ap|δ_m, U) is log-normal distributed. In that case, to specify the PDF, only the first three moments are needed, which follow from derivatives of the CF. As shown in Appendix A.4 both methods yield almost identical results, and since the log-normal approach is significantly faster, we use it hereafter, unless otherwise stated.

To summarise, the major changes compared to the F18 model are the following:

1. To determine p(δ_m, U) we

   • updated the calculation of the variance $\langle {\delta }_{\text{m},U}^2\rangle$ and of the skewness $\langle {\delta }_{\text{m},U}^3\rangle$ in Appendix A to general filter functions;

   • combine in Eqs. (26)–(29) two log-normal random variables for the positive and negative parts for compensated filters to obtain the final expression for any filter shape.

2. To determine p(N_ap|δ_m, U) we

   • calculate the characteristic function of galaxy shot-noise around a given matter density profile via Eq. (35);

   • use log-normal approximation or inverse Fourier transform Eq. (42) to obtain the PDF of shot-noise from its characteristic function.

3. To determine ⟨κ_{< ϑ}|δ_m, U⟩ we

   • updated the calculations of $\langle {\kappa }_{<\vartheta }{\delta }_{\text{m},U}\rangle$ and $\langle {\kappa }_{<\vartheta }{\delta }_{\text{m},U}^2\rangle$ to general filter functions (see Appendix A).

4. Simulation data

Before using our revised model in data analyses, it is mandatory to quantify its precision and range of validity. We use for this validation exercise three simulations suites:

• the full-sky gravitational lensing simulations described in Takahashi et al. (2017, hereafter T17), with which we carry out a detailed investigation of the model in a simple survey configuration;

• the cosmo-SLICS simulations, described in Harnois-Déraps et al. (2019), with which we validate our model on a independent simulation suite;

• the SLICS simulations, described in Harnois-Déraps et al. (2018), with which we construct a KiDS-1000 like covariance matrix.

4.1. T17 simulations

The T17 simulations are constructed from a series of nested cubic boxes with side lengths of L, 2L, 3L, … placed around a fixed vertex representing the observer’s position, with L = 450 Mpc h⁻¹. Each box is replicated eight times and placed around the observer using periodic boundary conditions. The number of particles per box is fixed to 2048³, which results in higher mass and spatial resolutions at lower redshifts. Within each box three spherical lens shells are constructed, each with a width of 150 Mpc h⁻¹, which are then used by the public code GRAYTRIX² to trace the light-ray trajectories from the observer to the last scattering surface³. With the N-body code GADGET2 (Springel et al. 2001) the gravitational evolution of dark matter particles without baryonic processes are followed from the initial conditions, which in turn are determined by use of second-order Lagrangian perturbation theory. The initial linear power spectrum followed from the Code for Anisotropies in the Microwave Background (CAMB; Lewis et al. 2000) with Ω_m = 1 − Ω_Λ = 0.279, Ω_b = 0.046, h = 0.7, σ₈ = 0.82, and n_s = 0.97. The matter power spectrum agrees with theoretical predictions of the revised Halofit (Takahashi et al. 2012) within 5%(10%) for k < 5(6) h Mpc⁻¹ at z < 1. In order to account for the finite shell thickness and angular resolution, T17 provide correction formulae, which we repeat in Appendix B. Although various resolution options are available, for our purpose the realisations with a resolution of NSIDE = 4096 are sufficient.

We use the publicly available matter density contrast maps to create a realistic lens galaxy catalogue that mimics the second and third redshift bins of the luminous red galaxies sample constructed from the KiDS-1000 data (Vakili et al. 2019), as shown by the solid lines in Fig. 1. The reason to mock the LRG sample is that the galaxy bias for this kind of galaxies can be roughly described with a constant linear bias, which is needed for the analytical model. We excluded the lowest-redshift lens bin, first because of its low galaxy number density (n₀ = 0.012 gal arcmin⁻²) in which the shot-noise level is significant, and second because the density field is more non-linear, and hence we expect the log-normal approximation to break down. Since there is a significant overlap between the KiDS-1000 sources and the lenses in the fourth LRG redshift bin, we reject it as well. To create our lens galaxy samples we first project the T17 3D density maps δ_m,3D following the n(z) shown as the step functions in Fig. 1 to get two δ_m,2D maps. For both maps we then distribute galaxies following a Poisson distribution with parameter λ = n(1+b δ_m,2D), where b is a constant linear galaxy bias and n is chosen such that the galaxy number density is n₀ = 0.028 gal arcmin⁻² for the second bin (hereafter the low-redshift bin $z_{\text{l}}^{\text{low}}$) and n₀ = 0.046 gal arcmin⁻² for the third lens bin (hereafter the high-redshift bin $z_{\text{l}}^{\text{high}}$). Since our method requires a constant linear galaxy bias, we specify a bias of 1.72 for lens bin two and 1.74 for lens bin three, similar to those reported in Vakili et al. (2019). F18 found this linear bias assumption to be accurate enough for year 1 data of the Dark Energy Survey, which is similar in constraining power to our target KiDS data (but we note that an investigation of higher-order biasing is underway in Friedrich et al., in prep.).

Fig. 1. Lens galaxy redshift distribution constructed from the T17 simulation given the true n(z) of the second $z_{\text{l}}^{\text{low}}$ and third $z_{\text{l}}^{\text{high}}$ LRG bin (Vakili et al. 2019). The black dashed line shows the redshift of the source galaxies.

Fig. 2. Redshift distributions of the second and third LRG (lens) bins and the last two KiDS-1000 (source) bins of the SLICS simulations. The n(z) are scaled such that a comparison is possible.

Fig. 3. Different filters U used in this work to verify the new model. For all filters we scaled the first bin value to 1/arcmin−2 for comparison. The corresponding Q-filters are shown in Fig. B.1. The wide Mexican filter extends up to 150′.

In our validation test, we use a shear grid at a single source plane located at z = 0.8664, indicated by the black dashed line in Fig. 1. F18 showed that the model works for realistic redshift distributions, and this choice simplifies the generation of our source catalogues. Furthermore, in order to determine a realistic covariance matrix, we transform the shear field into an observed ellipticity field by adding shape noise to the shear grid as

$$\begin{array}{c}\hfill {ϵ}^{\text{obs}}=\frac{{ϵ}^{\text{s}}+g}{1+{ϵ}^{\text{s}}{g}^{\ast }},\end{array}$$(43)

where ϵ^obs, ϵ^s, and γ are complex numbers, and the asterisk ( * ) indicates complex conjugation. The source ellipticities ϵ^s per pixel are generated by drawing random numbers from a Gaussian distribution with width

$$\begin{array}{c}\hfill {\sigma }_{\text{pix}}=\frac{{\sigma }_{ϵ}}{\sqrt{n_{\text{gal}}{A}_{\text{pix}}}}\approx 0.29,\end{array}$$(44)

where A_pix is the pixel area of the shear grid, and the effective number density n_gal and σ_ϵ are chosen such that they are consistent with the KiDS data. While this transformation is valid in terms of the reduced shear g = γ/(1 − κ), we use throughout this paper the approximation γ ≈ g, as the typical values for the convergence are small, |κ|≪1. We neglect the intrinsic alignment of galaxies in this work.

4.2. Extracting the model components from the T17 simulations

In order to validate the different components of our model, we need to extract p(δ_m, U), p(N_ap), and ⟨γ_t|𝒬⟩ from the simulation. The first two follow directly by smoothing the maps of the projected density contrast and the lens galaxy with the corresponding filters. This smoothing can be performed in two different ways. The first is to use the HEALPY function QUERY_DISC, which finds all pixel centres that are located within a given radius, whereas the second approach uses the HEALPY function SMOOTHING, with a given beam window function created by the function BEAM2BL. The two approaches result in PDFs that differ slightly, since the QUERY_DISC does not reproduce an exact top-hat, while the SMOOTHING approach is only over a finite ℓ-range. Nevertheless, we found that both approaches are consistent for NSIDE = 4096 well within the uncertainty we estimate from 48 sub-patches (see discussion below), and hence we use the second approach which is significantly faster.

The tangential shear information ⟨γ_t|𝒬⟩ is measured for each quantile Q by the software TREECORR (Jarvis et al. 2004) in 15 log-spaced bins with angular separation Θ/20 < ϑ < Θ, where Θ is the size of the filter being used. For the top-hat filter we measured the shear profiles from 6′< ϑ < 120′, corresponding to a filter with a size of 120′. We note here that for all measured shear profiles the shear around random points is always subtracted, which ensures that the shear averaged over all quantiles for one realisation vanishes by definition.

In order to have an uncertainty for all three model quantities, we divide the full-sky map into 48 sub-patches, such that each patch has a size of approximately 859.4 deg². For p(δ_m, U) and p(N_ap) we determined for each sub-patch one distribution, such that we were able to calculate a standard deviation from 48 values for each bin in the PDF. For the covariance matrix we use 10 out of the 108 realisations and divide each full-sky map in 48 sub-patches, which then results in a covariance matrix measured from 480 fields. Furthermore, both for the covariance and for the error bars in the plotted shear profiles we use Eq. (43) to create noisy shear profiles for each sub-patch, which are then re-scaled to the effective KiDS-1000 area (see Giblin et al. 2021).

4.3. Cosmo-SLICS

We use the cosmo-SLICS simulations described in Harnois-Déraps et al. (2019) to determine the validity regime of our revised model for different cosmologies. These are a suite of weak lensing simulations sampling 26 points (listed in Table B.1) in a broad cold dark matter (CDM) parameter space, distributed in a Latin hypercube to minimise interpolation errors. Specifically, the matter density Ω_m, the dimensionless Hubble parameter h, the normalisation of the matter power spectrum σ₈, and the time-independent equation-of-state parameter of dark energy w₀ are varied over a range that is large enough to complement the analysis of current weak lensing data (see e.g. Harnois-Déraps et al. 2021). Each simulation follows 1536³ particles inside a cube of co-moving side length L_box = 505 h⁻¹ Mpc and n_c = 3072 grid cells on the side, starting with initial conditions produced with the Zel’dovich approximation. Moreover, the cosmo-SLICS evolve a pair of simulations at each node, designed to suppress the sampling variance (see Harnois-Déraps et al. 2019, for more details). Each cosmological model is ray-traced multiple times to produce 50 pseudo-independent light cones of size 100 deg².

For each realisation, we create KiDS-1000-like sources and KiDS-LRG-like lens catalogues, following the pipeline described in Harnois-Déraps et al. (2018); notably, we reproduce exactly the source galaxy number density and n(z) that is used in Asgari et al. (2021), who report a total number density n_gal = 6.93 arcmin⁻² and a redshift distribution estimated from self-organising maps (see Wright et al. 2020). These mock galaxies are then placed at random angular coordinates on 100 deg² light cones. In contrast to the T17 mocks, we test our model with two source redshift bins, corresponding to the KiDS-1000 fourth and fifth tomographic bins (hereafter $z_{\text{s}}^{\text{low}}$ and $z_{\text{s}}^{\text{high}}$). The source galaxies are assigned a shear signal γ from a series of lensing maps, following the linear interpolation algorithm described in Sect. 2 in Harnois-Déraps et al. (2018). For our lens sample we opted to include the second and third tomographic bin of the LRG galaxies described in Vakili et al. (2019) ($z_{\text{l}}^{\text{low}}$ and $z_{\text{l}}^{\text{high}}$). Compared to the T17 values, the n(z) of the cosmo-SLICS LRG mocks have a coarser redshift resolution of the simulations. Moreover, the n(z) vary slightly for different underlying cosmologies, due to variations in the relation between co-moving distance and redshift. Following Vakili et al. (2019), we generate our LRG catalogues assuming a constant linear galaxy bias of 1.72 and 1.74, with a galaxy number density of n₀ = 0.028 gal arcmin⁻² and n₀ = 0.046 gal arcmin⁻².

4.4. SLICS

In total the SLICS⁴ are a set of over 800 fully independent realisations similar to the fiducial ΛCDM cosmo-SLICS model. The underlying cosmological parameters for each run are the same, fixed to Ω_m = 0.2905, Ω_Λ = 0.7095, Ω_b = 0.0473, h = 0.6898, σ₈ = 0.826 and n_s = 0.969 (see Hinshaw et al. 2013). For Fourier modes k < 2.0 h Mpc⁻¹, the SLICS and cosmo-SLICS three-dimensional dark matter power spectrum P(k) agrees within 2% with the predictions from the Extended Cosmic Emulator (see Heitmann et al. 2014), followed by a progressive deviation for higher k-modes (Harnois-Déraps et al. 2018). We use the SLICS to estimate a reliable covariance matrix, which, combined with the cosmo-SLICS, allows us to test our model for a simulation that is independent of T17. Similar to the T17 simulations, the signal of the SLICS is combined with the randomly oriented intrinsic shapes ϵ^s to create ellipticities, whereas ϵ^s is drawn from a Gaussian distribution with width σ_ϵ directly since the shear information are given here per galaxy. We added an additional layer of realism and used a redshift-dependent shape noise that better reproduces the data properties. Specifically, we used σ_ϵ = 0.25 and 0.27 for the source bins, as reported in Giblin et al. (2021).

4.5. Extracting the SLICS and cosmo-SLICS data vector

The extraction of the data vector for the SLICS and cosmo-SLICS analyses is similar to the T17 case, where shape noise was not included for the cosmo-SLICS data vector to better capture the cosmological signal. Another slight difference is that the light cones are now square, which accentuates the edge effects when the aperture filter overlaps with the light-cone boundaries. In principle, it is possible to weight the outer rims for each N_ap map, so that the whole map can be used. Although this would increase our statistical power, it could also introduce a systematic offset. We opted instead to exclude the outer rim for each realisation resulting in an effective area of (10 − 2Θ)² deg² with Θ the size of the corresponding filter. This procedure also ensures that roughly the same number of background galaxies are used to calculate the shear profile around each pixel.

5. Testing the revised model

We used the simulations described in Sect. 4 to test our revised model and its accuracy in predicting shear profiles. Following the results of F18 we chose a top-hat filter of 20′ as our starting point and we considered a number of more general filters with a similar angular extent, shown in Fig. 3. Our motivation for studying these filters is as follows: We use a Gaussian filter to test whether the model performs well for non-constant but positive filters; the “adapted” filter is the filter that results from B20; the “Mexican” filter removes the local minimum at ϑ ∼ 40′; the “broad Mexican” has a larger width; finally, the “wide Mexican” suppresses the negative tail. In order to lower the amplitude of the negative part while keeping a similar width, we adjusted the upper bound of the wide-Mexican filter to conserve the compensation to 150′, which makes it better suited to large contiguous survey areas.

Before comparing our model to the simulations, we note that we are using here the revised model even for the top-hat filter, for which we could instead use the F18 model directly. Notably, the derivations of $\langle {\delta }_{\text{m},U}^2\rangle$ and $\langle {\delta }_{\text{m},U}^3\rangle$ are identical in the revised model, and we show in the following plots for the top-hat filters that both models yield almost identical results in predicting the shear profiles with a top-hat filter. Therefore, from here on, we only show results from the revised model. In the following three sections, we validate the key model ingredients introduced in Sects. 3.1–3.3.

5.1. Validating p(δ_m, U)

We show in Fig. 4 the PDF of the smoothed two-dimensional density contrast for all six filters, and for the two lens bins. We see by inspecting the different panels that the predictions agree with the simulations for the two lens bins within 1 σ cosmic variance expected for KiDS-1000. We note here that this PDF cannot be measured in real data, and that the real test for the accuracy of our model are the shear signals, with larger uncertainties due to shape noise. Nevertheless, for the top-hat and the Gaussian we have an agreement between model and simulation well within the 1 σ, which indicates that the log-normal approximation for theses filters is good. The other filters show stronger deviations when using a log-normal approximation, but these are weaker when the negative part of the filter approaches zero (wide Mexican) or when the width of the filter increases (broad Mexican), although the negative part of the broad Mexican is stronger than for the Mexican filter. This indicates that probing on larger scales either with a broader or wider filter the log-normal approximation is more accurate. Furthermore, when using the bi-variate log-normal approach discussed in Sect. 3.1, the residuals are even more suppressed, and thus we cannot recognise differences in the match between predicted and measured PDF for all compensated filters. Although the model for the compensated filters is not as good as for the non-negative filters (top-hat and Gaussian), the revised model remains consistent throughout with the T17 simulations.

Fig. 4. PDF of δm, U smoothed with the filters shown in Fig. 3. The orange shaded region is the standard deviation of 48 sub-patches scaled by a 777.4/859.4, where 777.4 deg2 is the effective survey area of KiDS-1000 (see Giblin et al. 2021) and 859.4 deg2 is the area of one sub-patch. The red dashed curve corresponds to a log-normal PDF with the measured moments $\langle {\delta }_{\text{m},U}^2\rangle$ and $\langle {\delta }_{\text{m},U}^3\rangle$ from the smoothed T17 density maps, and indicates the accuracy using a log-normal PDF. The green and the black dashed lines are both from the model; the green corresponds to the PDF of δm, U when using log-normal and the black using the bi-variate approach described in Eq. (26). Lower panels: residuals Δp(δm, U) of all lines with respect to the simulations.

5.2. Validating p(N_ap)

We show in Fig. 5 how well the model can predict p(N_ap) given the galaxy distributions described in Sect. 4. As for p(δ_m, U), the best matches are observed for the non-negative filters, where the simple log-normal PDF is used. For the compensated filters with the bi-variate log-normal p(δ_m, U) we note a slight deviation in the skewness of p(N_ap). These discrepancies are not seen when placing galaxies at random positions regardless of any underlying matter density field as shown in Fig. A.1, which indicates that they must originate either from p(δ_m, U) or from the ⟨w_ϑ|δ_m, U⟩ term (we set the latter to 0 for uniform random fields). It might be that the deviations seen in p(N_ap) are exclusively caused by the deviations in p(δ_m, U), but since they are much smaller, we expect that the assumptions made in computing ⟨w_ϑ|δ_m, U⟩ induce additional inaccuracies. Nevertheless, we show next that these deviations result in shear signals whose residuals are well within the statistical uncertainties of Stage III weak lensing surveys such as KiDS-1000. However the accuracy of the ⟨w_ϑ|δ_m, U⟩ term will likely need to be improved for future surveys like Euclid, as discussed in Sect. 3.2.

Fig. 5. PDF of Nap calculated with the filters U in Fig. 3. The orange lines are determined with the simulations and the orange shaded region is the standard deviation from 48 sub-patches. The black dashed lines correspond to the results from the new model, and for comparison the red dashed line in the upper left panel is from the old model. Lower panels: residuals Δp(Nap) of all lines with respect to the simulations.

5.3. Validating ⟨γ_t|𝒬⟩

Having quantified the accuracy of the basic ingredients of our model, we are now in a position to compare the predicted and measured shear profiles. This is a major result of our paper, which is shown in Fig. 6. Following G18, we used five quantiles and we measured the shear profiles up to 120′ (or 150′ for the wide Mexican case). For the top-hat, Gaussian, and wide Mexican filters we see no significant deviations between the model and the simulations. For the adapted and the smaller Mexicans the shear profiles show minor discrepancies in some quantiles and at large angular scales, but are always consistent within the KiDS-1000 accuracy. The shapes of the signals are affected by the choice of the filter. We can observe shifts in the peak positions and changes in the slope of the signals especially at small scales. This will allow us in the future to select filters that optimise the signal-to-noise ratio of the measurement, while being clean of systematics related to small-scale inaccuracies. Finally, we show in Fig. B.2 that for the compensated filter the difference in using the proposed bi-variate log-normal approach is slightly more accurate than using a plain log-normal. Although the difference does not change the final results noticeably, and although it introduces some inconsistency in the sense that we use a bi-variate approach for p(δ_m,U) but not for ⟨κ_{< θ}|δ_m,U⟩⁵, we decided to stay with the proposed ansatz because it is slightly more accurate, and we plan to use p(N_ap) in future analysis.

Fig. 6. Predicted shear profiles for the two lens samples (dashed black line) and measured shear profiles (in orange) for the new model with filter U. The orange shaded region is the standard deviation on the mean from 48 sub-patches, scaled to the KiDS-1000 area. The residuals between model and simulations were tested to determine whether they can be erased when the PDF of the aperture number is fixed to the measured value from T17, but the same discrepancies were present.

In order to check whether the discrepancies seen for some compensated filters yield biased results, we performed an MCMC analysis. As our data vector we used the T17 shear profiles shown in Fig. 6, where we made a conservative cut and included only scales above 14′, since as shown in F18 the model is not fully accurate for small angular scales. For the comparison we decided to use the adapted filter and the top-hat filter to have one analysis with and one without these discrepancies. Furthermore, since the mean aperture mass summed over all quantiles vanishes per definition, one of the five shear signals is fully determined by the others, and so we discarded for all cases the middle quantile with the lowest signal. Thus, we ended up with data and model vectors of size 88. As explained previously, we measured our covariance matrix from ten T17 simulations, each divided into 48 sub-patches, for a total of 480 sub-patches. We note here that the galaxy number density can slightly deviate between the different realisations due to the Poisson sampling. Given the amplitude of these small fluctuations, these can be safely neglected. Next we de-biased the inverse covariance matrix C⁻¹ following Hartlap et al. (2007),

$$\begin{array}{c}\hfill C^{-1}=\frac{n-p-2}{n-1}{\widehat{C}}^{-1},\end{array}$$(45)

where n is the number of simulations (480) and p the size of the data vector (88). Finally, given our data d measured from only one noise-free T17 realisation, and our model vector m, we measured the χ² statistics as

$$\begin{array}{c}\hfill {\chi }^2={[\mathbf{m}-\mathbf{d}]}^T{C}^{-1}[\mathbf{m}-\mathbf{d}].\end{array}$$(46)

Given this set-up we ran an MCMC varying the matter density parameter Ω_m and normalisation of the power spectrum σ₈ for the adapted and the top-hat filters, where we marginalised over the biases of the lens samples. As shown in Fig. 7 the analysis with the adapted filter results in a biased inference for the Ω_m-σ₈-plane (although still within 1σ); this is not the case for the top-hat filter. We note here that this bias is due to the systematic offset in the slope of the highest quantile, which in turn is is highly sensitive to Ω_m. Since the amplitude of the shear profiles are correct and these are highly correlated with the $S_8={\sigma }_8\sqrt{{\mathrm{\Omega }}_{\text{m}}/0.3}$ parameter, the contours shift to smaller σ₈ values in order to compensate for the bigger Ω_m value⁶. In the next section we calibrate the model to investigate whether this systematic bias can be corrected.

Fig. 7. MCMC results for the top-hat and adapted filter using the model and the T17 simulations as our data vector and a covariance matrix calculated from ten T17 realisations each divided into 48 sub-patches. For the adapted filter a systematic bias for σ8 and Ωm is found, although it cancels out for the $S_8={\sigma }_8\sqrt{{\mathrm{\Omega }}_{\text{m}}/0.3}$ parameter. The contours here are marginalised over the lens galaxy bias parameters.

5.4. Calibrating the model

In this section we calibrate the remaining small inaccuracies of the analytical model seen in Fig. 6 which result in the systematic bias we had observed in the parameter constraints shown in Fig. 7. For this we decided to divide out for each quantile the residuals between the model, γ_MT, at the T17 cosmological parameters, p_T17, and the noiseless shear profiles measured from the T17 simulations, γ_T17, such that the calibrated model at parameters p is defined as

$$\begin{array}{c}\hfill {\gamma }_{\text{M,cal}}(p)={\gamma }_{\text{M}}(p)\frac{{\gamma }_{\text{T17}}}{{\gamma }_{\text{M}}(p_{\text{T17}})}·\end{array}$$(47)

Since we used the n(z) combinations of the fiducial cosmo-SLICS shown in Fig. 2 to validate the calibration, we decided to use the n(z) shown in Fig. 1, and in order to have the source n(z) as close as possible to the one of the cosmo-SLICS we averaged several T17 shear grids at different redshifts for the same realisation weighted by the source n(z) shown in Fig. 2. In Fig. B.4 we show the calibration vectors for the highest and lowest quantile for the top-hat and adapted filter, where it can be seen that the different lens n(z) is more important than the source n(z).

Next, in order to investigate whether the calibration decreases the systematic biases we performed another MCMC analysis on independent simulations, where our data vector is the fiducial cosmology from the cosmo-SLICS shear profiles shown in Fig. 8, with the original model in red and the calibrated one in black. As before we used the adapted filter and the top-hat filter. The match between the predicted and measured shear profiles is slightly degraded compared to the T17 simulations, which could be caused by edge effects (contiguous full-sky vs. 100 deg² patches), the smaller statistic (41 253 deg² vs. 5000 deg²), or differences in the underlying matter power spectrum P(k) that is used in the model. We use the Takahashi et al. (2012) HALOFIT function throughout this paper, which is calibrated on the same N-body code that is used to create the T17 simulations (Springel et al. 2001, GADGET2), and which is known to have an excess power of 5–8% in the mildly non-linear regime (Heitmann et al. 2014). The cosmo-SLICS, in contrast, are produced from cubep³m (Harnois-Déraps et al. 2013), whose P(k) agrees better with the Cosmic Emulator of Heitmann et al. (2014). We tested different choices of power spectra models calculated with the PYCCL package (Chisari et al. 2019), but found differences in the predicted shear profiles that are negligible compared to the expected KiDS-1000 uncertainties.

Fig. 8. Shear profiles for the top-hat filter (left) and for the adapted filter (right) for the fiducial cosmology of cosmo-SLICS. The orange lines are the mean shear profiles and the orange shaded region is the expected KiDS-1000 uncertainty. The red dashed line corresponds to the original model and the black to the calibrated model.

Discarding again the middle quantile with the lowest signal, using four different redshift combinations (Sect. 4.3) and the signal at all scales because the model is calibrated at all scales, we have data and model vectors of size 160. In this scenario we calculated our covariance matrix from 614 SLICS simulations⁷ with shape noise that mimics the KiDS-1000 data. After de-biasing the inverse covariance matrix C⁻¹ with Eq. (45) we calculated the χ² with Eq. (46). Given this set-up we ran multiple MCMC, where we used the original model and the calibrated model. As shown in Fig. 9 the calibrated model for the adapted filter results compared to the original model in less a biased inference. Interestingly, the results for the top-hat filter seen in Fig. B.3 are slightly more biased than the calibrated model for the adapted filter. Since this offset is still inside 1σ, it is likely to be only a statistical fluke due to the remaining residual between model and cosmo-SLICS simulations. The constraining power between top-hat and adapted filter are different because the smoothing scales of the two filters were not adjusted as in Burger et al. (2020), and are here sensitive to different physical scales. Nevertheless, we show in Table 1 the resulting constrains for both filters, where it is seen that the calibration moves the results, also for the top-hat filter, closer to the truth.

Fig. 9. MCMC results for the adapted filter using the original and calibrated model. The data vector is calculated from the fiducial cosmology of cosmo-SLICS and a covariance matrix from 614 SLICS realisations. It is clearly seen that the calibrated model is less biased than the original one. The contours are marginalised over the lens galaxy bias parameters.

In order to compare our results with those of G18, who derived constraints of ${\mathrm{\Omega }}_{\text{m}}=0.{26}_{-0.03}^{+0.04}$ and $S_8=0.{90}_{-0.08}^{+0.10}$ with their fiducial analysis, we need to multiply our uncertainty intervals by $\sqrt{777.4/1321}$ to account for the smaller area of KiDS-1000 (777.4 deg²) compared to the DES Y1 area (1321 deg²). Furthermore, we exclusively used information about the shear profiles, whereas G18 also used the mean aperture number in each quantile. For this work we were a bit sceptical about using the aperture number here for the compensated filters because we have significant residual discrepancies between model and simulation, which would affect our analysis. The match of the shear profiles in turn is very accurate in our simulations, which shows that they are robust against uncertainties in P(N_ap)⁸. For instance, if one monotonically transforms the N_ap values, the predicted P(N_ap) changes, but the segmentation into quantiles is not affected, hence the shear profiles would remain the same. In order to use P(N_ap) in future analysis we need to model shot noise in the galaxy distribution and investigate if the residuals between model and simulations result in systematic biases, but we will keep this for future work. Nevertheless, we see that our constraints from using only information about the shear profiles can be similar to the ones in G18. In addition, due to the calibration method used here, smaller smoothing scales are available than those recommended in F18, where even the top-hat filter has significant deviations. This could allow us to further improve the significance for future DSS analyses or to investigate effects such as baryonic feedback and intrinsic alignments, which are typically relevant on scales < 10 Mpc h⁻¹.

6. Summary and conclusion

In our previous work (Burger et al. 2020) we showed that using compensated filters in the density split statistic (DSS) to quantify over- and underdense regions on the sky have advantages compared to the top-hat filter, both in terms of the overall S/N and of recovering accurately the galaxy bias term. Furthermore, we expect that compensated filters are less influenced by baryonic effect, since they are more confined in Fourier space and therefore are better in smoothing out large ℓ-modes where baryonic effects play an important role. This will be investigated in more detail in a follow-up paper, when we start dealing with real data. Gruen et al. (2016) demonstrated that the DSS is a powerful cosmological tool by constraining cosmological parameters with DSS measurements from the Dark Energy Survey (DES) First Year and Sloan Digital Sky Survey (SDSS) data, using the DSS model derived in Friedrich et al. (2018) which uses a top-hat filter. They found for the matter density parameter ${\mathrm{\Omega }}_{\text{m}}=0.{26}_{-0.03}^{+0.04}$, a constraint that agrees with and is competitive with the DES analysis of galaxy and shear two-point functions (see Abbott et al. 2018).

Following these works, we modify the model of Friedrich et al. (2018) in such a way that it can predict the shear profiles ⟨γ_t|𝒬⟩ for a given quantile Q of the aperture number N_ap for general filters (Gaussian and also compensated filters). This is achieved by recalculating the three basic ingredients, which are the PDF of the projected matter density contrast smoothed with the filter function, p(δ_m, U); the expectation value of the convergence inside a radius ϑ for a fixed smoothed matter density contrast, ⟨κ_{< ϑ}|δ_m, U⟩; and the distribution of N_ap for the given filter function U given the smoothed matter density contrast, p(N_ap|δ_m, U). For ⟨κ_{< ϑ}|δ_m, U⟩ we modified the calculation of the moments for general filters, while we introduced new approaches to calculate p(N_ap|δ_m, U) and p(δ_m, U) for compensated filters. For non-negative filters, δ_m, U is well described by a log-normal PDF, although we found significant deviations for compensated filters. To solve this issue we used a bi-variate log-normal ansatz, where we assumed that δ_m, U can be divided into two log-normal random variables with each separately following a log-normal distribution. For the calculation of p(N_ap|δ_m, U) we derived an expression for the corresponding characteristic function, which can be used either directly to calculate p(N_ap|δ_m, U) by inverse Fourier transformation or by calculating the first three moments, which then specify a log-normal distribution for p(N_ap|δ_m, U). The differences between these two approaches are considerably smaller than the statistical uncertainty, and so we used the latter approach because of its smaller computational time.

In order to validate the revised model, we compared it to the Takahashi et al. (2017) simulations. For non-negative filters like a top-hat or a Gaussian, no significant difference between the model and simulations for the PDF or the tangential shear profiles were detected. For compensated filters, however, we found some discrepancies in the predicted PDF of N_ap and shear signals, which results in a biased inference, although still inside 1σ. To correct this biased result, we calibrated the model to match the noiseless Takahashi et al. (2017) and tested the calibrated model with the independent fiducial cosmology of cosmo-SLICS (Harnois-Déraps et al. 2019). With the calibration applied, all systematic biases are removed, so we are confident that we can apply the model to Stage III surveys such as KiDS-1000. Although this calibration is less important for the top-hat and Gaussian filter, it is still an interesting approach because it allows even smaller scales to be used for both the shear profiles and the filter scales. The use of smaller scales, where the original models fail, makes it possible to increase the constraining power or to study baryonic effects that normally play an important role only at small scales.

After passing all these tests, we are confident that the revised model can be readily applied to Stage III lensing data. We note that a number of systematic effects related to weak lensing analyses will require external simulations, notably regarding the inclusion of secondary signal from the intrinsic alignments of galaxies, or from the impact of baryonic feedback on the matter distribution. However, our model is able to capture the uncertainty on the lens and source redshift distribution, the shape calibration bias, or the galaxy bias at a low computational cost, and is therefore ideally suited to perform competitive weak lensing analyses in the future.

Acknowledgments

We thank the anonymous referee for the very constructive and fruitful comments. This paper went through the whole KiDS review process, where we especially want to thank the KiDS-internal referee Benjamin Joachimi for his fruitful comments to improve this work. Further, we would like to thank Mike Jarvis for maintaining TREECORR and Ryuichi Takahashi for making his simulation suite publicly available. PB acknowledges support by the Deutsche Forschungsgemeinschaft, project SCHN342-13. OF gratefully acknowledges support by the Kavli Foundation and the International Newton Trust through a Newton-Kavli-Junior Fellowship and by Churchill College Cambridge through a postdoctoral By-Fellowship. JHD is supported by a STFC Ernest Rutherford Fellowship (project reference ST/S004858/1). Author contributions: all authors contributed to the development and writing of this paper.

References

Appendix A: Detailed derivations for the new model

In this appendix we show more detailed derivations of the results than in the main text. We start with the calculation of the variances or covariances in the flat-sky approximation, continue with calculations of the third-order moments and finish with the derivation of the PDF of the aperture number given a smoothed density contrast by use of the characteristic function.

A.1. Variance and skewness for general filters at leading order in perturbation theory

Although analytical possible we decided against using the bi-spectrum to calculate third-order moments like the skewness. Instead we use a formalism where we calculate the second- and third-order moments of the smoothed density contrasts within cylinder of physical radius R and physical length L using the flat-sky approximation shown in Appendix B in F18 for a top-hat filter, and apply it to our case with a general filter U. Numerically, this approach is faster since, as we will see below, it is possible to express the third-order moments in terms of second-order moments. Another advantage is that the projection is only along one dimension (radius of the cylinder) compared to the bi-spectrum, where the projection is at least along a 2D grid. Following F18 we start by considering a cylinder of radius R and length L. In Fourier space the top-hat filter for such a cylinder is given by

$$\begin{array}{c}\hfill W_{R,L}(\mathbf{k})=\frac{1}{{(2\pi )}^3}\frac{sin(L{k}_{||}/2)}{L{k}_{||}/2}\frac{2J_1(k_{\perp }R)}{k_{\perp }R}\equiv \frac{1}{{(2\pi )}^3}\frac{sin(L{k}_{||}/2)}{L{k}_{||}/2}{W}_R^{\text{th}}(k_{\perp }),\end{array}$$(A.1)

where J₁ is the first Bessel function, and k_|| and k_⊥ are the components of k parallel and orthogonal to the cylinder, respectively. The variance of the matter contrast within such a cylinder is given at leading order by

$$\begin{array}{c}\hfill \langle {\delta }_{R,L}^2\rangle (\chi )=D_{+}^2{\displaystyle \int }\mathrm{d}{k}_{||}{\mathrm{d}}^2{k}_{\perp }\frac{{sin}^2(L{k}_{||}/2)}{{(L{k}_{||}/2)}^2}{\left[W_R^{\text{th}}(k_{\perp })\right]}^2{P}_{\text{lin},0}(k_{\perp })\approx \frac{2\pi D_{+}^2}{L}{\displaystyle \int }\mathrm{d}kk{\left[W_R^{\text{th}}(k)\right]}^2{P}_{\text{lin},0}(k),\end{array}$$(A.2)

where the last expression follows from L ≫ R, and since the integration depends from now on only on k_⊥ we write the orthogonal component as k. The linear matter power spectrum of P_lin, 0 is calculated using Eisenstein & Hu (1998), and D₊ is the growth factor which depends on the conformal time. We note that the factor 1/L cancels out when projecting the moments in Eq. (A.28-A.30) using the Limber approximation (Limber 1953). According to this derivation for a top-hat filter we get for a general filter U that

$$\begin{array}{cc}\hfill W_{U_{\chi }}(k)& ={\displaystyle \underset{0}{\overset{2\pi }{\int }}\underset{0}{\overset{\infty }{\int }}\mathrm{d}r\mathrm{d}\vartheta U_{\chi }(r){\text{e}}^{-\mathrm{i}krcos\vartheta }=2\pi \underset{0}{\overset{\infty }{\int }}\mathrm{d}r{J}_0(kr)r{U}_{\chi }(r),}\hfill \end{array}$$(A.3)

where U_χ(r) = U(r/χ) = U(ϑ)/χ², with U(ϑ) being a filter measured in angular coordinates (see Fig. 3). Correspondingly, the variance of the matter density contrast for a general filter U in the flat-sky approximation is

$$\begin{array}{c}\hfill \langle {\delta }_{U,L}^2\rangle (\chi )=\frac{2\pi D_{+}^2}{L}{\displaystyle \int }\mathrm{d}kk{W}_{U_{\chi }}^2(k),P_{\text{lin},0}(k).\end{array}$$(A.4)

Following lines similar to those of Appendix B.4 of F18, the leading-order contribution to the skewness of matter density contrast for the general filter U can be calculated as

$$\begin{array}{cc}\hfill \langle {\delta }_{U,L}^3\rangle (\chi )& =3\widehat{c}{\pi }^{-1}{\displaystyle \int \int }\mathrm{d}{q}_1\mathrm{d}{q}_2{q}_1{q}_2{W}_{U_{\chi }}(q_1)W_{U_{\chi }}(q_2)P_{\text{lin},0}(q_1)P_{\text{lin},0}(q_2){\displaystyle \int }\mathrm{d}\varphi W_{U_{\chi }}\left(\sqrt{q_1^2+q_2^2+2q_1{q}_{2}cos\varphi }\right)F_2(q_1,q_2,\varphi )\hfill \\ \hfill & \equiv 3\widehat{c}{\pi }^{-1}{\displaystyle \int \int }\mathrm{d}{q}_1\mathrm{d}{q}_2{q}_1{q}_2{W}_{U_{\chi }}(q_1)W_{U_{\chi }}(q_2)P_{\text{lin},0}(q_1)P_{\text{lin},0}(q_2){\mathrm{\Phi }}_{U_{\chi }}(q_1,q_2),\hfill \end{array}$$(A.5)

where $\widehat{c}=\frac{4{\pi }^2{D}_{+}^4}{L^2}$. The function F₂ in a general ΛCDM universe is given by

$$\begin{array}{cc}\hfill F_2(q_1,q_2,\varphi )& =\frac{1}{2}(2+\frac{q_1}{q_2}cos\varphi +\frac{q_2}{q_1}cos\varphi )+(1+\mu )({cos}^2\varphi -1)=1+\frac{1}{2}cos\varphi (\frac{q_1}{q_2}+\frac{q_2}{q_1})-(1-\mu ){sin}^2\varphi ,\hfill \end{array}$$(A.6)

where μ results from perturbation theory and is a function of the growth factor D₊ (see Appendix B.1 in F18 for more details)⁹, and ϕ is the angle between the vectors with absolute values q₁ and q₂. Given the definition of W_Uχ in Eq. (A.3), Φ_Uχ can be written as

$$\begin{array}{c}\hfill {\mathrm{\Phi }}_{U_{\chi }}(q_1,q_2)=2\pi {\displaystyle \underset{0}{\overset{\infty }{\int }}}\mathrm{d}rr{U}_{\chi }(r){\displaystyle \int }\mathrm{d}\varphi J_0\left(r\sqrt{q_1^2+q_2^2+2q_1{q}_{2}cos\varphi }\right)F_2(q_1,q_2,\varphi ).\end{array}$$(A.7)

Next, we use Graf’s addition theorem (see e.g. Abramowitz & Stegun 1972), which states that

$$\begin{array}{c}\hfill J_0\left(\sqrt{q_1^2+q_2^2+2q_1{q}_{2}cos\varphi }\right)={\displaystyle \sum _{m=-\infty }^{\infty }{(-1)}^m{J}_m(q_1)J_m(q_2){\text{e}}^{\text{i}m\varphi }=J_0(q_1)J_0(q_2)+2\sum _{m=1}^{\infty }{(-1)}^m{J}_m(q_1)J_m(q_2)cos(m\varphi ),}\end{array}$$(A.8)

such that Φ_Uχ(q₁, q₂) becomes

$$\begin{array}{cc}& 2\pi {\displaystyle \underset{0}{\overset{\infty }{\int }}r{U}_{\chi }(r)\mathrm{d}r{\int }_0^{2\pi }\mathrm{d}\varphi [J_0(r{q}_1)J_0(r{q}_2)+2\sum _{m=1}^{\infty }{(-1)}^m{J}_m(r{q}_1)J_m(r{q}_2)cos(m\varphi )][1+\frac{1}{2}cos(\varphi )(\frac{q_1}{q_2}+\frac{q_2}{q_1})-(1-\mu ){sin}^2\varphi ]}\hfill \\ \hfill & =\underset{A}{\underbrace{2{\pi }^2(1+\mu ){\displaystyle \underset{0}{\overset{\infty }{\int }}}r{U}_{\chi }(r)\mathrm{d}r{J}_0(r{q}_1)J_0(r{q}_2)}}-\underset{B}{\underbrace{2{\pi }^2{\displaystyle \underset{0}{\overset{\infty }{\int }}}r{U}_{\chi }(r)\mathrm{d}r{J}_1(r{q}_1)J_1(r{q}_2)(\frac{q_1}{q_2}+\frac{q_2}{q_1})}}+\underset{C}{\underbrace{2{\pi }^2(1-\mu ){\displaystyle \underset{0}{\overset{\infty }{\int }}}r{U}_{\chi }(r)\mathrm{d}r{J}_2(r{q}_1)J_2(r{q}_2)}},\hfill \end{array}$$(A.9)

where we made use of the orthogonality of the trigonometric functions. Plugging Φ_Uχ(q₁, q₂) back into Eq. (A.5) and considering each term separately we get

$$\begin{array}{cc}\hfill A:& 3\widehat{c}{\pi }^{-1}{\displaystyle \int \int }\mathrm{d}{q}_1\mathrm{d}{q}_2{q}_1{q}_2{W}_{U_{\chi }}(q_1)W_{U_{\chi }}(q_2)P_{\text{lin},0}(q_1)P_{\text{lin},0}(q_2)2\pi {\displaystyle \underset{0}{\overset{\infty }{\int }}}\mathrm{d}rr{U}_{\chi }(r)\pi (1+\mu )J_0(r{q}_1)J_0(r{q}_2)\hfill \\ \hfill & =6\pi \widehat{c}(1+\mu ){\displaystyle \underset{0}{\overset{\infty }{\int }}\mathrm{d}rr{U}_{\chi }(r){[\int \mathrm{d}qq{W}_{U_{\chi }}(q)P_{\text{lin},0}(q)J_0(rq)]}^2,}\hfill \end{array}$$(A.10)

and by analogy

$$\begin{array}{cc}\hfill C:& 3\widehat{c}{\pi }^{-1}{\displaystyle \int \int \mathrm{d}{q}_1\mathrm{d}{q}_2{q}_1{q}_2{W}_{U_{\chi }}(q_1)W_{U_{\chi }}(q_2)P_{\text{lin},0}(q_1)P_{\text{lin},0}(q_2)2\pi \underset{0}{\overset{\infty }{\int }}\mathrm{d}rr{U}_{\chi }(r)\pi (1-\mu )J_2(r{q}_1)J_2(r{q}_2)}\hfill \\ \hfill & =6\pi \widehat{c}(1-\mu ){\displaystyle \underset{0}{\overset{\infty }{\int }}\mathrm{d}rr{U}_{\chi }(r){[\int \mathrm{d}qq{W}_{U_{\chi }}(q)P_{\text{lin},0}(q)J_2(rq)]}^2,}\hfill \end{array}$$(A.11)

and finally

$$\begin{array}{cc}\hfill B:& -3\widehat{c}{\pi }^{-1}{\displaystyle \int \int \mathrm{d}{q}_1\mathrm{d}{q}_2{q}_1{q}_2{W}_{U_{\chi }}(q_1)W_{U_{\chi }}(q_2)P_{\text{lin},0}(q_1)P_{\text{lin},0}(q_2)2\pi \underset{0}{\overset{\infty }{\int }}\mathrm{d}rr{U}_{\chi }(r)\pi J_1(r{q}_1)J_1(r{q}_2)[\frac{q_1}{q_2}+\frac{q_2}{q_1}]}\hfill \\ \hfill & =-12\pi \widehat{c}{\displaystyle \underset{0}{\overset{\infty }{\int }}\mathrm{d}rr{U}_{\chi }(r)\int \mathrm{d}{q}_1{q}_1^2{W}_{U_{\chi }}(q_1)P_{\text{lin},0}(q_1)J_1(r{q}_1)\int \mathrm{d}{q}_2{W}_{U_{\chi }}(q_2)P_{\text{lin},0}(q_2)J_1(r{q}_2).}\hfill \end{array}$$(A.12)

The following transformations provide a more compressed expression for $\langle {\delta }_{U_{\chi },L}^3\rangle (\chi )$, which can then be used to verify our derivation by comparing it with the result from F18 for a top-hat filter. For this, we rewrite the expression of Bessel functions in terms of $W_r^{\text{th}}(q)$ as

$$\begin{array}{cc}\hfill J_2(rq)& =\frac{1}{\mathit{rq}}{J}_1(rq)-\frac{1}{q}\frac{\mathrm{d}}{\mathrm{d}r}{J}_1(rq)=\frac{1}{\mathit{rq}}{J}_1(rq)-[r\frac{\mathrm{d}}{\mathrm{d}r}\frac{J_1(rq)}{\mathit{rq}}+\frac{1}{\mathit{rq}}{J}_1(rq)]=-\frac{1}{2}\frac{\mathrm{d}}{\mathrm{d}ln(r)}{W}_r^{\text{th}}(q),\hfill \end{array}$$(A.13)

$$\begin{array}{cc}\hfill J_0(rq)& =\frac{1}{\mathit{rq}}{J}_1(rq)+\frac{1}{q}\frac{\mathrm{d}}{\mathrm{d}r}{J}_1(rq)=\frac{1}{\mathit{rq}}{J}_1(rq)+[r\frac{\mathrm{d}}{\mathrm{d}r}\frac{J_1(rq)}{\mathit{rq}}+\frac{1}{\mathit{rq}}{J}_1(rq)]=W_r^{\text{th}}(q)+\frac{1}{2}\frac{\mathrm{d}}{\mathrm{d}ln(r)}{W}_r^{\text{th}}(q),\hfill \end{array}$$(A.14)

and with

$$\begin{array}{c}\hfill \frac{1}{\mathit{rq}}\frac{{\mathrm{d}}^2}{\mathrm{d}{r}^2}{J}_1(rq)=\frac{1}{2}\frac{{\mathrm{d}}^2}{\mathrm{d}{r}^2}{W}_r^{\text{th}}(q)-\frac{2}{r^2}[\frac{1}{\mathit{rq}}{J}_1(rq)-\frac{1}{q}\frac{\mathrm{d}}{\mathrm{d}r}{J}_1(rq)]=\frac{1}{2}\frac{{\mathrm{d}}^2}{\mathrm{d}{r}^2}{W}_r^{\text{th}}(q)+\frac{1}{r^2}\frac{\mathrm{d}}{\mathrm{d}ln(r)}{W}_r^{\text{th}}(q),\end{array}$$(A.15)