# Density Split Statistics: Joint Model of Counts and Lensing in Cells

O. Friedrich, D. Gruen, J. Derose, D. Kirk, E. Krause, T. Mcclintock, E.S.

Rykoff, S. Seitz, R.H. Wechsler, G.M. Bernstein, et al.

► To cite this version:

O. Friedrich, D. Gruen, J. Derose, D. Kirk, E. Krause, et al., Density Split Statistics: Joint
Model of Counts and Lensing in Cells. Physical Review D, 2018, 98 (2), pp.023508. 10.1103/Phys-
RevD.98.023508. hal-01846850

HAL Id: hal-01846850

https://hal.science/hal-01846850v1

Submitted on 2 Oct 2025

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
# Density split statistics: Joint model of counts and lensing in cells

O. Friedrich,<sup>1,2,*</sup> D. Gruen,<sup>3,4</sup> J. DeRose,<sup>5,3</sup> D. Kirk,<sup>6</sup> E. Krause,<sup>3</sup> T. McClintock,<sup>7</sup> E. S. Rykoff,<sup>3,4</sup> S. Seitz,<sup>2,1</sup>

R. H. Wechsler,<sup>5,3,4</sup> G. M. Bernstein,<sup>8</sup> J. Blazek,<sup>9,10</sup> C. Chang,<sup>11</sup> S. Hilbert,<sup>12,13</sup> B. Jain,<sup>8</sup> A. Kovacs,<sup>14</sup> O. Lahav,<sup>6</sup>

F. B. Abdalla,<sup>6,15</sup> S. Allam,<sup>16</sup> J. Annis,<sup>16</sup> K. Bechtol,<sup>17</sup> A. Benoit-Lévy,<sup>18,6,19</sup> E. Bertin,<sup>18,19</sup> D. Brooks,<sup>6</sup>

A. Carnero Rosell,<sup>20,21</sup> M. Carrasco Kind,<sup>22,23</sup> J. Carretero,<sup>14</sup> C. E. Cunha,<sup>3</sup> C. B. D'Andrea,<sup>8</sup> L. N. da Costa,<sup>20,21</sup> C. Davis,<sup>3</sup>

S. Desai,<sup>24</sup> H. T. Diehl,<sup>16</sup> J. P. Dietrich,<sup>12,13</sup> A. Drlica-Wagner,<sup>16</sup> T. F. Eifler,<sup>25,26</sup> P. Fosalba,<sup>27</sup> J. Frieman,<sup>16,11</sup>

J. García-Bellido,<sup>28</sup> E. Gaztanaga,<sup>27</sup> D. W. Gerdes,<sup>29,30</sup> T. Giannantonio,<sup>31,32,1</sup> R. A. Gruendl,<sup>22,23</sup> J. Gschwend,<sup>20,21</sup>

G. Gutierrez,<sup>16</sup> K. Honscheid,<sup>9,33</sup> D. J. James,<sup>34</sup> M. Jarvis,<sup>8</sup> K. Kuehn,<sup>35</sup> N. Kuropatkin,<sup>16</sup> M. Lima,<sup>36,20</sup> M. March,<sup>8</sup>

J. L. Marshall,<sup>37</sup> P. Melchior,<sup>38</sup> F. Menanteau,<sup>22,23</sup> R. Miquel,<sup>39,14</sup> J. J. Mohr,<sup>13,12,2</sup> B. Nord,<sup>16</sup> A. A. Plazas,<sup>26</sup> E. Sanchez,<sup>40</sup>

V. Scarpine,<sup>16</sup> R. Schindler,<sup>4</sup> M. Schubnell,<sup>30</sup> I. Sevilla-Noarbe,<sup>40</sup> E. Sheldon,<sup>41</sup> M. Smith,<sup>42</sup> M. Soares-Santos,<sup>16</sup>

F. Sobreira,<sup>43,20</sup> E. Suchyta,<sup>44</sup> M. E. C. Swanson,<sup>23</sup> G. Tarle,<sup>30</sup> D. Thomas,<sup>45</sup> M. A. Troxel,<sup>9,33</sup>

V. Vikram,<sup>46</sup> and J. Weller<sup>13,2,1</sup>

(DES Collaboration)

<sup>1</sup>Universitäts-Sternwarte, Fakultät für Physik, Ludwig-Maximilians Universität München,
Scheinerstr. 1, 81679 München, Germany

<sup>2</sup>Max Planck Institute for Extraterrestrial Physics, Giessenbachstrasse, 85748 Garching, Germany

<sup>3</sup>Kavli Institute for Particle Astrophysics and Cosmology, P. O. Box 2450, Stanford University,
Stanford, California 94305, USA

<sup>4</sup>SLAC National Accelerator Laboratory, Menlo Park, California 94025, USA

<sup>5</sup>Department of Physics, Stanford University, 382 Via Pueblo Mall, Stanford, California 94305, USA

<sup>6</sup>Department of Physics and Astronomy, University College London,
Gower Street, London WC1E 6BT, United Kingdom

<sup>7</sup>Department of Physics, University of Arizona, Tucson, Arizona 85721, USA

<sup>8</sup>Department of Physics and Astronomy, University of Pennsylvania,
Philadelphia, Pennsylvania 19104, USA

<sup>9</sup>Center for Cosmology and Astro-Particle Physics, The Ohio State University,
Columbus, Ohio 43210, USA

<sup>10</sup>Institute of Physics, Laboratory of Astrophysics, École Polytechnique Fédérale de Lausanne (EPFL),
Observatoire de Sauvernay, 1290 Versoix, Switzerland

<sup>11</sup>Kavli Institute for Cosmological Physics, University of Chicago, Chicago, Illinois 60637, USA

<sup>12</sup>Faculty of Physics, Ludwig-Maximilians-Universität, Scheinerstr. 1, 81679 Munich, Germany

<sup>13</sup>Excellence Cluster Universe, Boltzmannstr. 2, 85748 Garching, Germany

<sup>14</sup>Institut de Física d'Altes Energies (IFAE), The Barcelona Institute of Science and Technology, Campus
UAB, 08193 Bellaterra (Barcelona) Spain

<sup>15</sup>Department of Physics and Electronics, Rhodes University,
PO Box 94, Grahamstown 6140, South Africa

<sup>16</sup>Fermi National Accelerator Laboratory, P. O. Box 500, Batavia, Illinois 60510, USA

<sup>17</sup>LSST, 933 North Cherry Avenue, Tucson, Arizona 85721, USA

<sup>18</sup>CNRS, UMR 7095, Institut d'Astrophysique de Paris, F-75014 Paris, France

<sup>19</sup>Sorbonne Universités, UPMC Univ Paris 06, UMR 7095, Institut d'Astrophysique de Paris,
F-75014 Paris, France

<sup>20</sup>Laboratório Interinstitucional de e-Astronomia—LIneA, Rua Gal. José Cristino 77,
Rio de Janeiro RJ—20921-400, Brazil

<sup>21</sup>Observatório Nacional, Rua Gal. José Cristino 77, Rio de Janeiro RJ—20921-400, Brazil

<sup>22</sup>Department of Astronomy, University of Illinois, 1002 W. Green Street, Urbana, Illinois 61801, USA

<sup>23</sup>National Center for Supercomputing Applications, 1205 West Clark Street, Urbana, Illinois 61801, USA

<sup>24</sup>Department of Physics, IIT Hyderabad, Kandi, Telangana 502285, India

<sup>25</sup>Department of Physics, California Institute of Technology, Pasadena, California 91125, USA

<sup>26</sup>Jet Propulsion Laboratory, California Institute of Technology,
4800 Oak Grove Dr., Pasadena, California 91109, USA

<sup>27</sup>Institute of Space Sciences, IEEC-CSIC, Campus UAB, Carrer de Can Magrans,
s/n, 08193 Barcelona, Spain

<sup>28</sup>Instituto de Física Teórica UAM/CSIC, Universidad Autónoma de Madrid, 28049 Madrid, Spain
$^{29}$Department of Astronomy, University of Michigan, Ann Arbor, Michigan 48109, USA

$^{30}$Department of Physics, University of Michigan, Ann Arbor, Michigan 48109, USA

$^{31}$Institute of Astronomy, University of Cambridge, Madingley Road,
Cambridge CB3 0HA, United Kingdom

$^{32}$Kavli Institute for Cosmology, University of Cambridge, Madingley Road,
Cambridge CB3 0HA, United Kingdom

$^{33}$Department of Physics, The Ohio State University, Columbus, Ohio 43210, USA

$^{34}$Astronomy Department, University of Washington, Box 351580, Seattle, Washington 98195, USA

$^{35}$Australian Astronomical Observatory, North Ryde, NSW 2113, Australia

$^{36}$Departamento de Física Matemática, Instituto de Física, Universidade de São Paulo,
CP 66318, São Paulo, São Paulo 05314-970, Brazil

$^{37}$George P. and Cynthia Woods Mitchell Institute for Fundamental Physics and Astronomy,
and Department of Physics and Astronomy, Texas A&M University, College Station, Texas 77843, USA

$^{38}$Department of Astrophysical Sciences, Princeton University,
Peyton Hall, Princeton, New Jersey 08544, USA

$^{39}$Institució Catalana de Recerca i Estudis Avançats, E-08010 Barcelona, Spain

$^{40}$Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas (CIEMAT), 28040 Madrid, Spain

$^{41}$Brookhaven National Laboratory, Bldg 510, Upton, New York 11973, USA

$^{42}$School of Physics and Astronomy, University of Southampton, Southampton SO17 1BJ, United Kingdom

$^{43}$Instituto de Física Gleb Wataghin, Universidade Estadual de Campinas,
Campinas, São Paulo 13083-859, Brazil

$^{44}$Computer Science and Mathematics Division, Oak Ridge National Laboratory,
Oak Ridge, Tennessee 37831, USA

$^{45}$Institute of Cosmology and Gravitation, University of Portsmouth,
Portsmouth PO1 3FX, United Kingdom

$^{46}$Argonne National Laboratory, 9700 South Cass Avenue, Lemont, Illinois 60439, USA

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

(Received 24 October 2017; published 13 July 2018)

We present density split statistics, a framework that studies lensing and counts-in-cells as a function of foreground galaxy density, thereby providing a large-scale measurement of both 2-point and 3-point statistics. Our method extends our earlier work on trough lensing and is summarized as follows: given a foreground (low redshift) population of galaxies, we divide the sky into subareas of equal size but distinct galaxy density. We then measure lensing around uniformly spaced points separately in each of these subareas, as well as counts-in-cells statistics (CiC). The lensing signals trace the matter density contrast around regions of fixed galaxy density. Through the CiC measurements this can be related to the density profile around regions of fixed matter density. Together, these measurements constitute a powerful probe of cosmology, the skewness of the density field and the connection of galaxies and matter. In this paper we show how to model both the density split lensing signal and CiC from basic ingredients: a non-linear power spectrum, clustering hierarchy coefficients from perturbation theory and a parametric model for galaxy bias and shot-noise. Using N-body simulations, we demonstrate that this model is sufficiently accurate for a cosmological analysis on year 1 data from the Dark Energy Survey.

DOI: 10.1103/PhysRevD.98.023508

# I. INTRODUCTION

The large-scale structure (LSS) observed today is thought to originate from almost perfectly Gaussian density perturbations in the early Universe. This means that there was a complete symmetry in the abundance and amplitude of underdense and overdense regions in very early times. Gravitational attraction then caused initial overdensities to collapse to small but highly overdense structures such as

galaxy clusters, while initial underdensities expanded but stayed moderately underdense and, e.g., became voids. As a consequence the majority of the volume in the late-time Universe is underdense, compensated by the presence of few highly overdense spots. Or, in other words, a positive skewness in the distribution of density fluctuations emerges due to gravitational collapse.

A variety of probes have been used to study the statistical properties of the late-time density field and to thereby understand the physics of gravitational collapse as well as the processes responsible for the properties of the initial density fluctuations. So far, the most extensive studies have

<footnote>* Corresponding author.
oliverf@usm.uni-muenchen.de</footnote>
been carried out on the 2-point statistics of density fluctuations, i.e., on measuring the variance of density fluctuations as a function of scale. This has, e.g., been done through measurements of cosmic shear 2-point correlation functions [e.g., [1–10]], galaxy clustering [e.g., [11–13]] and galaxy-galaxy lensing [e.g., [14–20]] as well as combined measurements thereof [e.g., [21–23]].

While 2-point statistics are only sensitive to the overall amplitude of density fluctuations, higher-order statistics also know about the skewness arising from the different behaviour of underdense and overdense regions. This does not necessarily mean that higher-order statistics are better than 2-point statistics in discriminating between particular choices of cosmological parameters [24]. But they scale differently with parameters such as $\Omega_m$, $\sigma_8$, galaxy bias and galaxy stochasticity than their 2-point counterparts. Hence, in a cosmological analysis that varies a large number of parameters, probes that are sensitive to both 2-point and higher order statistics have the power to break degeneracies between these parameters [25–28].

Observations of higher-order statistical features of the density field include measurements of three point correlation functions [29], shear peak statistics [30–32], and the cluster mass function [33]. Also, a number of probes have been suggested (and in some cases measured in data) that study the correlation of 2-point statistics and background density. Chiang et al. [34] have measured this by means of the integrated bispectrum. Simpson et al. [35–37] have proposed a clipped power spectrum approach, where 2-point statistics are measured on the sky after excluding high density regions. They have shown that these measurements contain information complementary to the corresponding measurements on the full sky.

A similar direction was investigated by Gruen et al. [38] who separately measured the lensing power spectrum in underdense and overdense lines of sight. The framework presented in this paper is based on their concept of trough lensing. We will call it *density split lensing* when only lensing measurements are involved and *density split statistics* when it is combined with counts-in-cells measurements. This method can be summarized as follows: we consider a foreground (low redshift) population of galaxies and smooth their position field with a circular top-hat aperture. This smoothed density field is then used to divide the sky into sub-areas of equal size but distinct galaxy density. In this paper we consider in particular 5 sub-areas and call them quintiles of galaxy density. As a next step, we use a background (high redshift) population of galaxies to measure the tangential shear of these galaxies around a set of uniformly spaced points within the area of each density quintile. The resulting lensing signals trace the matter density contrast around regions of fixed foreground galaxy density. This data vector is then complemented by the histogram of counts-in-cells of the foreground galaxies to pin down their bias and stochasticity. As we show in this

paper, a cosmological analysis based on this density split data vector has a number of desirable features:

(i) it allows for an accurate analytic modeling with the help of cosmological perturbation theory and a nonlinear power spectrum,

(ii) it yields high signal-to-noise measurement,

(iii) it avoids systematics common to cosmic shear such as additive shear biases or intrinsic alignments (as long as tracer sample and source sample do not overlap in redshift),

(iv) it has a very intuitive interpretation.

This paper is a companion paper of Gruen et al. [39], where we present our actual data analysis, including tests for systematic effects as well as a description of how we estimate the covariance of our signal. This paper is presenting the modeling framework used in that analysis. Our Sec. II gives a general overview of density split statistics: we describe our data vector, explain how it can be modeled and also present forecasts on cosmological parameter constraints, both for a $\Lambda$CDM model and an extended model that allows gravitational collapse to behave differently than within general relativity. In Sec. III we describe the simulated data used in this work. Section IV explains details of the model presented in Sec. II. There we also compare individual components of this model directly to measurements in N-body simulations. In Sec. V we show that our model for a data vector combining density split lensing and counts-in-cells statistics is accurate enough to recover the cosmology underlying our N-body simulations. Any possible deviation between our model and the simulations is shown to be well within statistical uncertainties of year 1 data of the Dark Energy Survey (DES Y1).

In Appendix A, we review a number of differential equations that govern gravitational collapse. In Appendix B, we review the leading order perturbative calculation of the 3-point statistics of the cosmic density field for a general $\Lambda$CDM model. Appendix C qualitatively compares our model of the cosmic density PDF to a second set of N-body simulations as a complement to the comparison carried out in the main text. Appendix D derives properties of joint log-normal random fields and Appendix E repeats the validation of our model for an alternative shot-noise parametrization.

## II. DENSITY SPLIT STATISTICS: DATA VECTOR, MODELING AND FORECASTS

This section provides an introduction to the program of density split statistics. In Sec. II A we describe how we obtain the density split lensing signal and how this signal can be further complemented with information on galaxy bias and stochasticity from counts-in-cells. In Sec. II B we outline our modeling of this signal (but postponing technical details of this model to Sec. IV). In Sec. II C we provide forecasts on the cosmological information that can
be obtained with a measurement of density split statistics in year-1 data of the Dark Energy Survey (DES Y1).

## A. Measuring density split statistics

Density split lensing is a generalization of trough lensing [38] and can be described in three steps:

### 1. Splitting the sky into quantiles of different foreground galaxy density

Consider a sample of low-redshift galaxies that are tracing the line-of-sight density of matter with some redshift distribution $n_l(z)$. We will call these galaxies the foreground sample. For an angular radius $\theta_T$, which we will call the top-hat aperture radius, we define $N_T(\hat{n})$ to be the number of galaxies found within a radius $\theta_T$ around the point on the sky specified by the vector $\hat{n}$ on the unit sphere. The field $N_T(\hat{n})$ can be used to divide the sky into regions of different galaxy density. Gruen et al. [38] have done this by discretizing the sky with a HEALPIX.<fnref n="1" /> grid and sorting the pixels according to their value of $N_T(\hat{n})$. Then they considered the 20% of the pixels with the lowest values of $N_T$, calling them troughs. In the limit of a fine pixelization these pixels can be considered the most underdense quintile of the sky area. This can be generalized to the second most underdense quintile, the third most underdense quintile etc. or even to finer splits using more than just 5 quantiles.

We stick to a division into 5 quantiles (quintiles) throughout this paper. The upper panel of Fig. 1 illustrates such a subdivision on a patch of a simulated sky (from the Buzzard flock, see Sec. III and especially DeRose et al. [41] for details). There we use a top-hat aperture radius $\theta_T = 20'$ and the tracer galaxies have the redshift distribution that is displayed by the solid line in Fig. 2. Figure 1 shows the most underdense quintile of the simulated patch in cyan, the most overdense quintile in red, and the three intermediate quintiles in blue, green, and orange.

Note that the sum of the 5 lensing signals will vanish on average (we subtracted 1/5 times the shear around random points from each signal, but due to boundary effects their sum will not vanish exactly). This means that roughly 4 of the 5 signals contain independent information. We have not investigated, whether our choice of 5 quantiles is in any way optimal. Choosing 3 quantiles would leave us with 2 independent signals and would hence suffice to be sensitive to both the variance and skewness of the density field. 5 quantiles enable a sentitivity beyond the 3rd moment of the density field. And they also allow us to explicitly show, that the median universe is underdense (which we could not do with 4 quantiles). In Sec. IV we investigate different radii of our top-hat aperture and find that $\theta_T = 20'$ is the smallest

radius at which our model is reliable (given the redshift distribution we use in Gruen et al. [39]).

### 2. Tracing the mean dark matter density in each sky quintile with gravitational lensing

Now consider a second sample of galaxies at higher redshifts than the foreground sample (the source sample, see, e.g., the dashed and dotted redshift distributions in Fig. 2). As the light of these galaxies passes the large-scale structure of the foreground density distribution it undergoes gravitational lensing effects such as gravitational shear [see, e.g., [42]]. The density split lensing signal around each quintile of the sky is obtained by measuring the stacked radial profile of tangential shear around random points located within that quintile. These points are constrained to lie within the part of the sky covered by a certain quintile of galaxy density but are otherwise random in their location. Because these random points are split according to the density quintile they are located in, their stacked shear signals trace the average profile of density contrast around each quintile.

In the lower right panel of Fig. 1 we show the signals measured for each density quintile in our mock data. The points show the average measurement from 4 Buzzard realisations of DES year-1 data (using the highest redshift source population shown in Fig. 2) and the solid lines show predictions by the model presented in this paper. The error bars are derived from a set of log-normal realizations (using the FLASK tool by Xavier et al. [43]; in Gruen et al. [39] we describe in detail how we configured FLASK to generate our mock catalogs). Two main features of the density split lensing signals are apparent: first, the amplitude of the radial shear around the 20% most underdense pixels is lower than the amplitude of the tangential shear around the 20% most overdense pixels. This is reflecting the skewness of the cosmic density PDF. Secondly, the signal around points in the third quintile is still significantly negative, which reflects the fact that the median universe is underdense. A more subtle feature is the fact that the underdense signals fall off less rapidly with increasing scale than the overdense signals. This is because on large scales the density field becomes Gaussian and hence recovers its initial symmetry between overdensities and underdensities.

### 3. Measuring the average counts-in-cells in each density quintile to obtain additional information on galaxy bias and stochasticity

If galaxy counts and the matter density field were perfectly correlated, then a split of the sky by galaxy density would be identical to a split by matter density. Hence, in this limit the density split lensing signals would be independent of the bias of the tracer galaxies. In a realistic scenario however, shot-noise of the galaxies smears out our attempts to divide the sky into areas of

<footnote n="1">¹See Górski et al. [40] for details on HEALPIX.</footnote>
<figure><img src="image_3.png" /><figcaption>FIG. 1. Top panel: splitting the lines of sight in one DES-Y1 like Buzzard simulation into 5 quantiles of galaxy density (color coding from cyan, most underdense, to red, most overdense). The map uses a 20 arcmin top-hat radius and REDMAGIC galaxies with a redshift range of $0.2 \le z \le 0.45$. Bottom left: histogram of REDMAGIC galaxy counts in 20 arcmin radii (counts-in-cells). We show the mean histogram from 4 Buzzard realizations of DES-Y1 (black points), our model based on perturbation theory and cylindrical collapse (solid line) and a model that assumes the projected density contrast to be a Gaussian random field (dotted line). The color coding corresponds exactly to the density quantiles in the top panel. Bottom right: Lensing signals around random points split by the density quantile in which these points are located. We show the mean measurement from 4 Buzzard realisations (black points), our perturbation theory model (solid line) and a model that assumes projected density contrast and lensing convergence to be joint Gaussian random variables (dotted line). Color coding is the same as in the other panels. The asymmetry between the lensing signals around the most underdense and most overdense lines-of-sight indicates the skewness of the cosmic density PDF.</figcaption></figure>

different matter density. Hence the density split lensing signals obtain a dependence on galaxy bias, but also on galaxy stochasticity. Increasing the linear bias of galaxy clustering will sharpen the tracers' ability to distinguish between overdensities and underdensities. Thus, increasing this bias will increase the amplitudes of the signals. This means that linear bias is to some degree degenerate with the amplitude of density fluctuations, $\sigma_8$. But $\sigma_8$ and bias influence the third moments of the density field in different ways and their degeneracy is not complete. As a consequence, it is possible to obtain constraints on cosmological parameters from the lensing signals alone (cf. Sec. II C and the blue contour in the left panel of Fig. 3).

But additional information on bias and stochasticity nevertheless helps to tighten these constraints. In this paper we decided to add that information in the form of normalized quantiles of the counts-in-cells (CiC) histogram of the tracer galaxies: we measure the histogram of tracer counts within the same aperture that was used to identify our density quintiles. Then we identify the parts of this histogram that correspond to these quintiles (cf. lower left panel of Fig. 1). For each quintile $q$ we then determine the mean galaxy count in that quintile, $N_q$, and normalize it by the overall mean galaxy count in our aperture, $\bar{N}$. I.e., for each quintile we add $N_q/\bar{N}$ to our data vector. This indeed helps to tighten constraints on cosmological parameters
<figure><img src="image_4.png" /><figcaption>FIG. 2. Redshift distributions of the tracer galaxy sample and the source samples of our N-body realizations of DES-Y1.</figcaption></figure>

(cf. Sec. IIC and the green contour in the left panel of Fig. 3).

## B. Modeling density split statistics

We now outline a general framework for modeling the data vector described above, leaving details of this framework to Sec. IV. Unless stated differently, we will assume a flat $\Lambda$CDM universe throughout this paper.

Let us start by introducing the quantities whose relations need to be modeled. First, we denote with $\delta_{m,2D}$ the line-of-sight projection of the 3D density contrast according to

the redshift distribution $n_l(z)$ of our foreground galaxy sample, i.e.,

$$ \delta_{m,2D}(\hat{\mathbf{n}}) = \int dw q_l(w) \delta_{m,3D}(w\hat{\mathbf{n}}, w) \quad (2.1) $$

where $\hat{\mathbf{n}}$ denotes a unit vector on the sky, $w$ is comoving distance and the projection kernel $q_l(w)$ is given in terms of $n_l(z)$ as

$$ q_l(w) = n_l(z[w]) \frac{dz[w]}{dw}. \quad (2.2) $$

We furthermore define $\delta_{m,T}$ to be the average of $\delta_{m,2D}$ over top-hat filters with aperture radius $\theta_T$, i.e.

$$ \delta_{m,T}(\hat{\mathbf{n}}) = \int_{|\hat{\mathbf{n}}, \hat{\mathbf{n}}'| < \theta_T} d\Omega' \frac{\delta_{m,2D}(\hat{\mathbf{n}}')}{2\pi(1 - \cos\theta_T)}. \quad (2.3) $$

Here $|\cdot, \cdot|$ denotes the angular distance between two points on the sky.

We identify regions of different density by means of our foreground galaxy sample. When smoothed with a top-hat filter of radius $\theta_T$, these galaxies are biased and possibly stochastic tracers of $\delta_{m,T}(\hat{\mathbf{n}})$. Hence our model also needs to include a description of how $N_T(\hat{\mathbf{n}})$, the number of tracer galaxies found within an angular radius $\theta_T$ around the line-of-sight $\hat{\mathbf{n}}$, relates to $\delta_{m,T}(\hat{\mathbf{n}})$.

Finally, in order to describe the density split lensing signal, we need to consider the lensing convergence field

<figure><img src="image_5.png" /><figcaption>FIG. 3. Left panel: Forecast of $1\sigma$ and $2\sigma$ constraints on $\Omega_m$ and $\sigma_8$ achievable with density split lensing and counts-in-cells in DES Y1 data. The constraints are marginalized over $\Omega_b$, $n_s$, $h_{100}$, REDMAGIC galaxy bias $b$, and galaxy-matter correlation coefficient $r$. For the parameters $\Omega_b$, $n_s$, $h_{100}$ we have assumed the same flat priors as used in the DES Y1 combined probes analysis presented in DES Collaboration [21]. Right panel: $\Delta S_3/S_3$ measures relative deviations from our fiducial perturbation theory prediction of the scaling coefficient $S_3 \equiv (\delta^3)^2 / (\delta^2)^2$ [cf. Eq. (2.30) and Sec. IV for details]. It can hence be thought of as the Bispectrum amplitude. We show $1\sigma$ constraints on this parameter achievable with density split lensing and counts-in-cells in DES data alone (solid lines) and using additional information on cosmology from Planck (dashed lines, no lensing). The sharp cutoff of the contours for low values of $\Sigma_8$ is caused by the requirement that matter density and a shot-noise free galaxy density field must have a correlation coefficient $r \le 1$. All likelihoods are centered around our fiducial cosmology, i.e., the parameters describing the Buzzard simulations.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>FIG. 4. Forecast of posterior constraints on galaxy bias b and galaxy-matter correlation coefficient r achievable with density split lensing and counts-in-cells in DES Y1 data. The constraints are marginalized over $\Omega_m$, $\sigma_8$, $\Omega_h$, $n_s$, and $h_{100}$, where for the last three parameters we have assumed the same flat priors as used in the DES Y1 combined probes analysis presented in DES Collaboration [21].</figcaption></figure>

for our population of source galaxies. Given the source redshift distribution $n_s(z)$, the convergence $\kappa$ is given by the line-of-sight projection

$$ \kappa(\hat{\mathbf{n}}) = \int dw W_s(w) \delta_{m,3D}(w\hat{\mathbf{n}}, w), \quad (2.4) $$

where $W_s$ is the lensing efficiency, which is defined by

$$ W_s(w) = \frac{3\Omega_m H_0^2}{2c^2} \int_w^\infty dw' \frac{w(w' - w)}{w'a(w)} q_s(w'), \quad (2.5) $$

and $q_s(w) = n_s(z[w]) \frac{dz[w]}{dw}$ is the line-of-sight density of the sources. Smoothing the convergence field with a circular aperture of radius $\theta$ results in a field which we will denote by $\kappa_{<\theta}(\hat{\mathbf{n}})$.

Because of the isotropy of the universe, we will now omit the dependence of the above quantities on $\hat{\mathbf{n}}$. To model the density split lensing signal one needs to answer the following questions:

(i) Given the number of galaxies $N_T$ found around a line-of-sight $\hat{\mathbf{n}}$, what distribution can be inferred for the matter density contrast $\delta_{m,T}$ in that line-of-sight? i.e., what is the expectation value $\langle \delta_{m,T} | N_T \rangle$?

(ii) Given the matter density contrast $\delta_{m,T}$ in the line-of-sight $\hat{\mathbf{n}}$, what lensing convergence $\kappa_{<\theta}$ is expected inside an angular distance $\theta$ from that line-of-sight? i.e. what is the expectation value $\langle \kappa_{<\theta} | \delta_{m,T} \rangle$? The tangential shear profile around that line-of-sight can then be inferred from the convergence profile as

$$ \begin{aligned} \langle \gamma_t(\theta) | \delta_{m,T} \rangle &= \langle \kappa_{<\theta} | \delta_{m,T} \rangle - \langle \kappa_{\theta} | \delta_{m,T} \rangle \\ &= \frac{\cos \theta - 1}{\sin \theta} \frac{d}{d\theta} \langle \kappa_{<\theta} | \delta_{m,T} \rangle, \end{aligned} \quad (2.6) $$

where $\kappa_{\theta}$ is the average convergence at the radius $\theta$.

The first of the above questions can be answered in the form of a conditional PDF of $\delta_{m,T}$ given a certain value of $N_T$, i.e. $p(\delta_{m,T} | N_T)$. Using Bayes' theorem this can be written as

$$ p(\delta_{m,T} | N_T) = \frac{P(N_T | \delta_{m,T}) p(\delta_{m,T})}{P(N_T)}, \quad (2.7) $$

where $P(N_T | \delta_{m,T})$ is the probability of finding a number of galaxies $N_T$ given that the density contrast is $\delta_{m,T}$ and where $p(\delta_{m,T})$ and $P(N_T)$ are the total PDF of $\delta_{m,T}$ and the total probability of finding $N_T$ tracer galaxies. The average convergence profile around a circle with $N_T$ galaxies is then given by

$$ \begin{aligned} \langle \kappa_{<\theta} | N_T \rangle &= \int d\delta_{m,T} \langle \kappa_{<\theta} | \delta_{m,T}, N_T \rangle p(\delta_{m,T} | N_T) \\ &\approx \int d\delta_{m,T} \langle \kappa_{<\theta} | \delta_{m,T} \rangle p(\delta_{m,T} | N_T), \end{aligned} \quad (2.8) $$

where in the second step we have assumed that the expected convergence within $\theta$ only depends on the total matter density contrast within $\theta_T$.

We now divide the sky into different quintiles of tracer galaxy density. Let us denote with $Q[0.0, 0.2]$ the 20% of the lines-of-sight on the sky that have the lowest value of $N_T$. There will be a maximal value $N_T \le N_{max}$ in this quintile and the stacked convergence profile around these lines-of-sight is given by

$$ \begin{aligned} \langle \kappa_{<\theta} | Q[0.0, 0.2] \rangle &= \frac{1}{0.2} \left( \sum_{N < N_{\text{max}}} P(N) \langle \kappa_{<\theta} | N_T = N \rangle + \alpha \langle \kappa_{<\theta} | N_T = N_{\text{max}} \rangle \right). \end{aligned} \quad (2.9) $$
Here the factor $\alpha$ in the second term accounts for the fact that the lines-of-sight with exactly $N_T = N_{max}$ might have to be split between the quintile $Q[0.0, 0.2]$ and the quintile $Q[0.2, 0.4]$. It is given by

$$ \alpha = 0.2 - \sum_{N<N_{max}} P(N). \qquad (2.10) $$

This can be easily generalized to the other quintiles $Q[q_{min}, q_{max}]$ and also to the case of dividing the sky into more than 5 density regimes. Finally, we also add the average of the counts-in-cells in each quintile normalized by the mean galaxy count $\bar{N}$ to our data vector. For the quintile $Q[0.0, 0.2]$ this is given by

$$ \frac{\langle N_T | Q[0.0, 0.2] \rangle}{\bar{N}} = \frac{1}{0.2\bar{N}} \left( \sum_{N<N_{max}} P(N)N + \alpha N_{max} \right), \qquad (2.11) $$

which is also straightforward to generalize to other quan- tiles $Q[q_{min}, q_{max}]$.

The probabilities $P(N)$ can be computed from the normalization of Eq. (2.7). Hence, our model for the density split lensing signal needs the following three ingredients:

(i) The PDF of matter density contrast, smoothed with a top-hat filter of radius $\theta_T$,

$$ p(\delta_{m,T}). \qquad (2.12) $$

(ii) The expectation value of convergence inside a distance $\theta$ given the density contrast inside $\theta_T$,

$$ \langle \kappa_{<\theta} | \delta_{m,T} \rangle. \qquad (2.13) $$

(iii) The distribution of galaxy counts inside the top-hat radius $\theta_T$ given the density contrast within that radius,

$$ P(N_T | \delta_{m,T}). \qquad (2.14) $$

Gruen et al. [38] assumed $\delta_{m,T}$ and $\kappa_{<\theta}$ to have a joint Gaussian distribution. This allowed them to compute (i) and (ii) solely from the dark matter clustering power spectrum. To compute (iii) they assumed a linear galaxy bias and Poissonian shot-noise of the tracer galaxies. These assumptions allowed a sufficient model for their measurements made on DES Science Verification data. But as can be seen from the dotted lines in the lower panels of Fig. 1, a Gaussian model for the density PDF is not sufficient within the much smaller uncertainty of DES-Y1. Also, in Sec. IV we demonstrate that the shot-noise of the tracer galaxies cannot necessarily be assumed to be Poissonian. In this work we hence want to revise their model.

For each of the model components (i) and (ii) we investigate two different modeling approaches—a baseline approach and an approach with increased complexity used to check the validity of the baseline model. In the following we briefly outline each approach. The reader interested in details of each modeling ansatz is referred to Sec. IV. Readers who are not interested in this technical part of the paper should feel free to skip Sec. IV.

(i) Baseline model for $p(\delta_{m,T})$:

In our fiducial model we assume $\delta_{m,T}$ to be a log-normal random field [44]. The PDF of such a variable can, e.g., be fixed by specifying the variance $\langle \delta_{m,T}^2 \rangle$ and skewness $\langle \delta_{m,T}^3 \rangle$. We predict the variance of $\delta_{m,T}$ from the nonlinear matter power spectrum [cf. [38]]. The latter is computed using HALOFIT [45,46] and an analytic transfer function [47]. We then use leading order perturbation theory to compute a scaling relation between the bispectrum and the power spectrum of the density field. Together with our power spectrum this fixes the skewness of $\delta_{m,T}$.

Alternative model for $p(\delta_{m,T})$:

As an alternative we compute the PDF $p(\delta_{m,T})$ from its cumulant generating function (see Sec. IV for a definition and further details). To model this function we use a cylindrical collapse approach based on the work of Bernardeau [48], Bernardeau and Valageas [49] and Valageas [50].

(ii) Baseline model for $\langle \kappa_{<\theta} | \delta_{m,T} \rangle$:

In our fiducial model we assume that $\kappa_{<\theta}$ can be expressed as the sum of two random variables,

$$ \kappa_{<\theta} = \kappa_{<\theta,\text{corr}} + \kappa_{<\theta,\text{uncorr}}, \qquad (2.15) $$

where $\kappa_{<\theta,\text{uncorr}}$ is assumed to be completely un-correlated to $\delta_{m,T}$ and hence does not contribute to the density split lensing signal. As a result we have

$$ \langle \kappa_{<\theta} | \delta_{m,T} \rangle = \langle \kappa_{<\theta,\text{corr}} | \delta_{m,T} \rangle. \qquad (2.16) $$

We assume a joint log-normal PDF for the two random variables $\delta_{m,T}$ and $\kappa_{<\theta,\text{corr}}$. The expectation value $\langle \kappa_{<\theta,\text{corr}} | \delta_{m,T} \rangle$ is then fixed by specifying the moments

$$ \langle \delta_{m,T}^2 \rangle, \langle \delta_{m,T}^3 \rangle \qquad (2.17) $$

as well as

$$ \langle \kappa_{<\theta} \delta_{m,T} \rangle = \langle \kappa_{<\theta,\text{corr}} \delta_{m,T} \rangle \qquad (2.18) $$

and

$$ \langle \kappa_{<\theta} \delta_{m,T}^2 \rangle = \langle \kappa_{<\theta,\text{corr}} \delta_{m,T}^2 \rangle. \qquad (2.19) $$
Second order moments are again computed from our nonlinear power spectrum while third order moments are inferred from perturbation theory. (The introduction of $\kappa_{<\theta,uncorr}$ is only necessary for consistency reasons: a joint log-normal PDF of $\delta_{m,T}$ and $\kappa_{<\theta}$ characterized by the above moments would be inconsistent with the variance $\langle \kappa_{<\theta}^2 \rangle$ derived from our power spectrum.)

Alternative model for $\langle \kappa_{<\theta} | \delta_{m,T} \rangle$:

As an alternative we compute $\langle \kappa_{<\theta} | \delta_{m,T} \rangle$ from the joint cumulant generating function of the variables $\delta_{m,T}$ and $\kappa_{<\theta}$. This function can also be modelled in a cylindrical collapse approach.

For model component (iii) we also investigate two different modeling approaches—one ansatz introducing 2 free parameters and one ansatz introducing 3 free parameters. We find that our simulated tracer catalogs are well described by the 2-parametric model. But—anticipating real galaxies to behave more complicated than simulated ones—we do not consider either of these models as our baseline model and instead apply both approaches to DES data in Gruen et al. [39]. We summarize both ansatzes here, but the interested reader is again referred to Sec. IV for details of each ansatz.

(iii) Model 1 for $P(N_T|\delta_{m,T})$:

In our fiducial model we introduce an auxilliary field $\delta_{g,T}$ such that our foreground galaxies are Poissonian tracers of that field, i.e.,

$$P(N_T|\delta_{g,T}) = \frac{[\bar{N}(1 + \delta_{g,T})]^{N_T}}{N_T!} e^{-\bar{N}(1+\delta_{g,T})}. \quad (2.20)$$

$\delta_{g,T}$ can be thought of as a smooth (shot noise free) galaxy density contrast. We then assume that $\delta_{g,T}$ and $\delta_{m,T}$ are joint log-normal random variables with

$$\langle \delta_{g,T}^2 \rangle = b^2 \langle \delta_{m,T}^2 \rangle, \quad \langle \delta_{g,T}^3 \rangle = b^3 \langle \delta_{m,T}^3 \rangle \quad (2.21)$$

and

$$\langle \delta_{g,T} \delta_{m,T} \rangle = br \langle \delta_{m,T}^2 \rangle. \quad (2.22)$$

The parameters $b$ and $r$ will be called galaxy bias and galaxy stochasticity and are free parameters of the model.

Model 2 for $P(N_T|\delta_{m,T})$:

As an alternative we assume $P(N_T|\delta_{m,T})$ to be a generalization of the Poisson distribution, that allows for

$$\langle N_T^2 | \delta_{m,T} \rangle \neq \langle N_T | \delta_{m,T} \rangle + \langle N_T | \delta_{m,T} \rangle^2, \quad (2.23)$$

i.e. for a shot-noise that is either enhanced or suppressed with respect to the Poisson case. The

enhancement of shot-noise is also allowed to be a function of $\delta_{m,T}$ of (approximately) the form

$$\frac{\langle N_T^2 | \delta_{m,T} \rangle - \langle N_T | \delta_{m,T} \rangle^2}{\langle N_T | \delta_{m,T} \rangle} \approx \alpha_0 + \alpha_1 \delta_{m,T}. \quad (2.24)$$

This model introduces an alternative bias parameter $\tilde{b}$ such that

$$\langle N_T | \delta_{m,T} \rangle = \bar{N}[1 + \tilde{b} \delta_{m,T}]. \quad (2.25)$$

For the model components (i) and (ii) and on the scales considered in this paper, the baseline and alternative modeling approaches yield almost indistinguishable predictions (cf. Sec. IV). For component (iii) the modeling ansatzes 1 and 2 are not necessarily identical, because they introduce a different number of degrees of freedom. Figure 1 as well as all parameter contours shown in this paper are using the baseline model for components (i) and (ii) and ansatz 1 for component (iii). The predictions derived from different modeling approaches are however almost indistinguishable (cf. Fig. 5).

### C. Data vector and forecasts on parameter constraints

#### 1. Binning and scales

Throughout this paper, we assume that one sample of tracer galaxies is used to identify density quintiles and that the lensing profiles around these quintiles are measured with two source redshift bins (cf. redshift distributions in Fig. 2). To identify the different density quintiles, we use a top-hat filter with fiducial smoothing radius of $\theta_T = 20'$.

We measure the density split lensing signal in 24 log-spaced angular bins with

$$5' < \theta < 600'. \quad (2.26)$$

But in the parameter forecasts and likelihood contours shown in the following, we exclude bins with scales $\le \theta_T$. This reduces all lensing signals to 17 log-spaced angular bins with

$$20' \le \theta < 600'. \quad (2.27)$$

The scales with $\theta < \theta_T$ are excluded from our analysis for a number of reasons: first, the signal-to-noise ratio of our lensing profiles drops quickly below $\theta_T$ (cf. Fig. 1). Second, we chose our fiducial value of $\theta_T = 20'$ because we still trust the modeling of the density PDF described in Sec. IV on this scale, and we do not want smaller angular scales to contribute to our fiducial data vector. Third, the approximation made in Eq. (2.8) might fail at scales smaller than our aperture.
<figure><img src="image_7.png" /><figcaption>FIG. 5. The PDF of projected density contrast in Buzzard compared to several models for various smoothing scales ($\delta_{m,T}$ is the projected density contrast $\delta_m$ smoothed by our top-hat radius $\theta_T$). In the upper panel of each plot the black line shows a histogram of $\delta_{m,T}$ measured in an all-sky map from Buzzard. The grey lines show histograms measured in 14 DES year1 shaped patches of that all-sky map. The blue lines show the PDF predicted by our tree-level computation of the cumulant generating function, the green lines show the PT-motivated log-normal model and the red lines show a Gaussian PDF with the same variance as the other two models. The bottom panels of each plot are showing the ratio of each PDF to the one measured in the Buzzard all-sky. For all aperture radii our HALOFIT power spectrum is predicting a standard deviation of $\delta_{m,T}$ that is $\gtrsim 2\%$ higher than what we find in Buzzard (cf. Fig. 6). For $\theta_T = 20'$ and $30'$ this is the dominant source of mismatch.</figcaption></figure>
<table><caption>TABLE I. Model parameters of the forecast described in Sec. IIC 2. The second column shows our fiducial values (the values describing the Buzzard simulations). The third column shows the parameter priors used to cut our prediction for the posterior distribution of best-fit parameters. The priors on $\Omega_b$, $n_s$, and $h_{100}$ are informative and chosen to be the same as used by the DES Collaboration [21]. The prior on $r$ is needed for mathematical consistency. And the other priors only have to be introduced formally since we are approximating our analytic posterior by an MCMC. The 4th column shows the standard deviation of each parameter (as computed from the MCMC) after marginalization over all other parameters. The 5th column shows the same standard deviations but for the case where also $\Delta S_3/S_3$ is introduced as a free parameter of our model. In column 6, we again fix $\Delta S_3/S_3$ but assume that only the lensing part of the data vector is used. These forecasts can be compared to the actual errors we find in Gruen et al. [[39], their tables 2 and 3] which are close despite marginalizing over systematic uncertainties.</caption><thead><tr><th>Parameter</th><th>Fiducial value (in Buzzard)</th><th>Prior in likelihood analysis</th><th>DES-Y1 constraints without ΔS<sub>3</sub>/S<sub>3</sub></th><th>DES-Y1 constraints with ΔS<sub>3</sub>/S<sub>3</sub></th><th>DES-Y1 constraints without ΔS<sub>3</sub>/S<sub>3</sub> (lensing only)</th></tr></thead><tbody><tr><td>σ<sub>8</sub></td><td>0.82</td><td>[0.2, 1.6]</td><td>±0.05</td><td>±0.10</td><td>±0.11</td></tr><tr><td>Ω<sub>m</sub></td><td>0.286</td><td>[0.1, 0.9]</td><td>±0.04</td><td>±0.06</td><td>±0.06</td></tr><tr><td>Ω<sub>b</sub></td><td>0.047</td><td>[0.01, 0.07]</td><td>...</td><td>...</td><td>...</td></tr><tr><td>n<sub>s</sub></td><td>0.96</td><td>[0.7, 1.3]</td><td>...</td><td>...</td><td>...</td></tr><tr><td>h<sub>100</sub></td><td>0.7</td><td>[0.55, 0.91]</td><td>...</td><td>...</td><td>...</td></tr><tr><td>b</td><td>1.618</td><td>[0.8, 2.5]<br>(lensing only: [0.0, 4.5])</td><td>±0.11</td><td>±0.27</td><td>±0.57</td></tr><tr><td>r</td><td>0.956</td><td>[0.0, 1.0]</td><td>±0.10</td><td>±0.11</td><td>±0.18</td></tr><tr><td>ΔS<sub>3</sub>/S<sub>3</sub></td><td>0.0</td><td>[-1.0, 2.0]</td><td>...</td><td>±0.20</td><td>...</td></tr></tbody></table>

The sum of all 5 lensing signals will be very close to zero (though not exactly zero because of masking effects), so that they are not an independent set of observations. Hence, we only include the 2 most overdense and the 2 most underdense quintiles in our analysis, i.e., the cyan, blue, orange, and red lines in Fig. 1. For the same reason, we also use the normalized mean galaxy count only in four of the five quintiles of the counts-in-cells histogram to complement the lensing measurement. The total number of data points in our data vector is thus

$$ 
\begin{aligned}
N_{\text{data}} &= N_{\text{lens}} + N_{\text{CiC}} \\
&= [N_{\text{quant}} - 1] \cdot N_{\text{source}} \cdot N_{\text{ang}} + [N_{\text{quant}} - 1] \\
&= 4 \cdot 2 \cdot 17 + 4 \\
&= 140.
\end{aligned} 
\qquad (2.28) 
$$

## 2. Model parameters and forecast on constraining power

We now investigate what constraints on model parameters can be expected from the above data vector when measured in DES-Y1 data. For this, we assume the optimistic case that component (iii) of our model is sufficiently described by two parameters, i.e., by modeling approach 1 in the previous section.<fnref n="2" /> This means that our model is determined by the following 7 parameters:

(1) $\Omega_m$: present total matter density in units the critical density of the universe,

(2) $\sigma_8$: amplitude of present day density fluctuations in spheres of 8 Mpc/h radius as predicted by the linear power spectrum,

(3) $\Omega_b$: present density of baryonic matter in units of the critical density of the universe,

(4) $n_s$: the spectral index of the primordial power spectrum,

(5) $h_{100}$: the present day Hubble parameter in units of 100 Mpc/(km/s),

(6) $b$: linear bias of the tracers with respect to matter density [cf. Eq. (2.21)],

(7) $r$: correlation coefficient between $\delta_{m,T}$ and $\delta_{g,T}$ [cf. Eq. (2.22)].

We summarize these parameters and our fiducial values for them in Table I. Throughout this paper, we assume the universe to be flat.

If $\pi_\alpha$ and $\pi_\beta$ are any two of the above parameters then let $\hat{\pi}_{\alpha,ML}$ and $\hat{\pi}_{\beta,ML}$ be maximum likelihood estimates of these parameters based on a measurement of density split statistics. The covariance of $\hat{\pi}_{\alpha,ML}$ and $\hat{\pi}_{\beta,ML}$ can be estimated by

$$ \text{Cov}[\hat{\pi}_{\alpha,ML}, \hat{\pi}_{\beta,ML}]^{-1} \approx \frac{\partial \mathbf{d}_{\text{th}}^T}{\partial \pi_{\alpha}} \cdot \mathbf{C}_d^{-1} \cdot \frac{\partial \mathbf{d}_{\text{th}}}{\partial \pi_{\beta}}, \quad (2.29) $$

where $\mathbf{d}_{\text{th}}$ is our model of the density split data vector and $\mathbf{C}_d$ is the covariance matrix of a measurement of this signal in DES-Y1. We will in the following use an estimate of $\mathbf{C}_d$ from log-normal mocks and real DES Y1 shape noise (see Sec. III B for a brief summary of our covariance estimation and Gruen et al. [39] for details). The parameter covariance computed with Eq. (2.29) can then be used to approximate the expected distribution of our best-fit parameter estimates with a multivariate Gaussian distribution.

Since the three parameters $\Omega_b$, $n_s$, and $h_{100}$ are only poorly constrained by our data vector we are forced to assume prior knowledge on them. To do so, we cut the

<footnote n="2">²Gruen et al. [39] will also consider modeling approach 2.</footnote>
Gaussian posterior predicted from the parameter covariance with flat informative priors. These priors are chosen to be the same used by the DES Collaboration [21]. For reasons of mathematical consistency we furthermore have to demand that $r \in [0, 1]$. These hard cuts of our originally Gaussian approximation to the posterior distribution of best-fit parameters make it difficult to marginalize over individual parameters. We hence numerically evaluate our analytic posterior with a Monte-Carlo Markov-Chain (MCMC). This chain is used in the following visualizations. Since we are using an MCMC to trace our analytic posterior, we have to formally define priors for all other model parameters. These are chosen to be flat and to extend well beyond the single-parameter standard deviations of the posterior. In the third column of Table I, we summarize the priors chosen for each model parameter.

The left panel of Fig. 3 shows the $1\sigma$ and $2\sigma$ constraints achievable in the $\Omega_m - \sigma_8$ plane. These contours are marginalized over the other model parameters, using the priors mentioned above. The blue contours assume that only the density split lensing signal has been used while the green contours allow for complementary information from the tracer counts-in-cells histogram. In the fourth column of Table I we show the standard deviation of each parameter as found in our approximation to the posterior (and assuming the full data vector, including lensing and counts-in-cell). The 6th column shows the same constraints but using only the lensing part of the data vector. Note that even without the CiC part of the data vector density split statistics can constrain the tracer bias $b$ and the galaxy-matter correlation coefficient $r$ (see Fig. 4 for the full posterior distributions of those parameters).

Density split statistics is complementary to an analysis based on 2-point statistics not just because it has a different dependence on the connection of galaxies and matter, but also since it is sensitive to higher order moments of the density field. We demonstrate this by introducing an additional degree of freedom in our model, described by an additional parameter:

(8) $\Delta S_3/S_3$: a factor multiplying all third order statistics in our predictions. The notation for this parameter is based on the ratio

$$S_3 = \frac{\langle \delta^3 \rangle}{\langle \delta^2 \rangle^2} \qquad (2.30)$$

which connects third and second moments of the density contrast and hence characterizes the amplitude of the density bispectrum (see Sec. IVA 2 for details). In our fiducial setup we compute it from leading order perturbation theory and $\Delta S_3/S_3$ hence describes a relative deviation from that result.

Within the $\Lambda$CDM model and at leading order in perturbation theory, the scaling between 2-point and 3-point statistics of the density field is almost independent of the cosmological parameters $\Omega_m$ and $\sigma_8$ [51]. Hence, a value of

$\Delta S_3/S_3 \neq 0$ would allow for deviations from the leading order result that cannot be compensated by changing $\Omega_m$ or $\sigma_8$. Such deviations could be caused nonstandard physics of dark matter and dark energy that affect overdensities and underdensities differently (see, e.g., Lue et al. [52,53]; though $f(R)$ modified gravity theories have been shown to largely preserve the $\Lambda$CDM scaling, cf. Borisov and Jain [54], Jain and Zhang [55]). Alternatively, they could indicate a break down the perturbative scaling relations due to highly nonlinear evolution of the density field or any small scale baryonic physics that do not follow the scaling relations of perturbation theory [cf. [28,51,56]].

In the right panel of Fig. 3 we show how density split statistics including lensing and counts-in-cells can simultaneously constrain $\Delta S_3/S_3$ and the parameter

$$\Sigma_8 = \sigma_8 \sqrt{\Omega_m}, \qquad (2.31)$$

even after marginalization over the other model parameters. We also project how these constraints will improve when moving to year 5 (Y5) data of DES or when adding cosmological information from the cosmic microwave background. For the latter we estimated the parameter covariance in a Planck chain<fnref n="3" /> and added this covariance as an additional Gaussian prior around our fiducial cosmology.

With DES Y1 alone, we will be able to constrain the amplitude of third order statistics of the density field to about 20% accuracy (cf. last column of Table I). Combining DES Y5 and Planck, this improves to about 5%. And this is even underestimating the power of DES-Y5: to project our constraints onto year-5 we simply divided our covariance by a factor of 4 in order to account for the larger area of the final DES survey. But this does not take into account the fact that DES Y5 will also be deeper than DES-Y1, which reduces shape noise and opens up the possibility of analyzing a larger number of redshift bins.

### III. SIMULATED DATA AND COVARIANCE MATRIX

In this work we use two sets of simulated data: the Buzzard galaxy catalogs which are constructed from high-resolution N-body simulations [41,58,59] and simulated random fields on the sky generated by the FLASK tool [43]. We briefly describe these data sets in the following sections.

#### A. Buzzard mock galaxy catalogs

Here we describe the key aspects of the Buzzard galaxy catalogs for the purposes of this work and refer the reader to more detailed descriptions elsewhere [41,58,59]. We use four independent realizations of a DES Y1-like survey in version Buzzard-v1.1. These catalogs were constructed

<footnote n="3">³plikHM_TT_lowTEB in COM_CosmoParams_base-plikHM_R2.00.tar.gz from the Planck legacy archive https://pla.esac.esa.int/pla/, no lensing, cf. Planck Collaboration et al. [57].</footnote>
from N-body simulations run using L-GADGET2 [60], a version of GADGET2 modified for memory efficiency. Second-order Lagrangian perturbation theory initial conditions [61] were employed using 2LPTIC [61], and light-cones were output on the fly. Each galaxy catalog is built from a set of three nested lightcones using progressively larger volume and lower resolution at higher redshifts. The force resolution in each box is 20, 35, and 53 kpc/h with the boundaries between the lightcones falling at redshifts $z = 0.34$ and $z = 0.9$.

The galaxy catalogs are constructed from the lightcones using the ADDGALS algorithm [59] which assigns galaxy luminosities and positions based on the relation between redshift, r-band absolute magnitude, and large-scale density, $p(\delta|M_r, z)$, found in a subhalo abundance matching (SHAM) model [[62,63], e.g.], in a high resolution N-body simulation. Spectral energy distributions (SEDs) are given to each simulated galaxy by finding a SDSS DR7 galaxy [64] that has a close match in $M_r$ and distance to its fifth nearest neighbor galaxy and assigning the SDSS galaxy's SED to the simulated galaxy. Galaxy sizes and ellipticities were assigned by drawing from distributions fit to high resolution SuprimeCam i'-band data [65]. Observed magnitudes in griz were generated by redshifting the SEDs to the observer frame and integrating over the DES passbands. Photometric errors were added using the DES Y1 multi object fitting (MOF) depth maps.

The effects of weak gravitational lensing are calculated using the multiple-plane raytracing algorithm CALCLENS [66]. The raytracing is done on a $n_{side} = 4096$ HEALPIX [40] grid yielding an effective angular resolution of 0.85'. At each lens plane, the inverse magnification matrix of the ray closest to every galaxy is interpolated to the galaxy's 3D position and used to shear and magnify the galaxy.

With galaxy catalogs with magnitudes, sizes and lensing effects in hand, REDMAGIC and METACALIBRATION [67–69] samples of galaxies are selected from the simulations in an effort to approximate the selections done in the data. In the case of REDMAGIC, the same algorithm which is used for selection in the data is applied to the simulations, resulting in a tracer galaxy catalog of equivalent volume density and at least comparable bias. For the METACALIBRATION sample, as we do not run full image simulations, we must make approximate cuts on signal to noise of the galaxies to create a facsimile of the source sample in the data. For in depth comparisons of these simulated samples to their data counterparts see DeRose et al. [41].

For the density split analysis in this work, we use REDMAGIC high density tracer galaxies ($L > 0.5L_{star}$, $\rho = 10^{-3}$ Mpc$^{-3}$ comoving density) selected at a REDMAGIC photometric redshift estimate of $0.2 < z < 0.45$. For the source redshift split lensing signals, we bin source galaxies by the expectation value of their $p(z)$ as estimated with BPZ [70,71] from the Buzzard mock photometry. Bin limits are chosen such that the true mean

redshifts of the bins match the mean redshifts of the two highest redshift source samples defined in Hoyle et al. [71].

## B. Simulated density and convergence fields from FLASK and Covariance matrix

For testing the numerical implementation of the model described in the following section and for estimating a covariance matrix of the density split lensing and counts-in-cells signals in the Buzzard simulations and the DES data, we use log-normal realizations of matter density and convergence fields. We summarize their properties here, with details given in Appendix A of Gruen et al. [39].

We generate these fields as all-sky HEALPIX maps of matter density and convergence using the FLASK software [43]. For the matter field, we choose the true redshift distribution of REDMAGIC galaxies in the Buzzard simulations, selected as described in Sec. III A. The matter field is sampled by a tracer population with REDMAGIC density, bias of $b = 1.54$, and Poissonian noise, from which lines of sight of different density are identified by the same algorithm as in Buzzard or in the data. For the convergence fields, we choose the estimated redshift distributions of DES source galaxies in the two highest redshift bins of Hoyle et al. [71]. Log-normal parameters of the density and convergence fields are set by the perturbation theory formalism described in the following section.

960 of these realizations are used to estimate large-scale structure and shot noise contributions to the covariance matrix of the signals modeled herein. The contribution of shape noise is estimated by measuring the lensing signals in actual DES Y1 data with 960 realizations of the METACALIBRATION shape catalog [69] in which each galaxy ellipticity was rotated by a random angle.

# IV. MODELING DETAILS AND COMPARISON TO SIMULATIONS

In this section we present the approximations used to compute the ingredients (i), (ii), and (iii) of our model that are listed at the end of Sec. II B. We also test the model ingredients (i) and (iii) directly against our N-body simulations.

Section IVA describes our model for the PDF of projected density contrast within the top-hat smoothing radius $\theta_T$, $p(\delta_{m,T})$. Section IV B describes how we model the convergence profile $\langle K_{<\theta} | \delta_{m,T} \rangle$ around apertures of fixed density contrast $\delta_{m,T}$. And in Sec. IV C we describe our modeling of the probability $P(N_T|\delta_{m,T})$ of finding $N_T$ tracer galaxies in an aperture with density contrast $\delta_{m,T}$.

Section IV D summarizes our fiducial model and the approximations used therein.

## A. Projected density PDF

The computation of the PDF of the density field when smoothed by top-hat filters has e.g. been addressed by Bernardeau [48], Bernardeau and Valageas [49], and
Valageas [50] (see also more recent developments in [28,72–75] which however do not yet enter our formalism). Bernardeau and Valageas [49] demonstrated how to extend methods for the computation of the 3D density PDF to the weak lensing aperture mass which is a projected quantity. In the following we show how to modify their formalism in order to compute the line-of-sight density PDF in angular circles of radius $\theta_T$. To do so, we have to consider the cumulant generating function (CGF) of the field $\delta_{m,T}(\hat{\mathbf{n}})$. The moment generating function (MGF) of $\delta_{m,T}(\hat{\mathbf{n}})$ is defined as

$$ \psi(y) = \sum_k \frac{\langle \delta_{m,T}(\hat{\mathbf{n}})^k \rangle}{k!} y^k. \quad (4.1) $$

Due to the isotropy of the universe it does not depend on $\hat{\mathbf{n}}$. The CGF $\varphi(y)$ is given in terms of the MGF $\psi(y)$ as

$$ \begin{aligned} \varphi(y) &= \ln \psi(y) \\ &\equiv \sum_k \frac{\langle \delta_{m,T}(\hat{\mathbf{n}})^k \rangle_c}{k!} y^k, \end{aligned} \quad (4.2) $$

where in the last line we have defined the connected moments or cumulants $\langle \delta_{m,T}(\hat{\mathbf{n}})^k \rangle_c$ of $\delta_{m,T}(\hat{\mathbf{n}})$. The CGF of $\delta_{m,T}(\hat{\mathbf{n}})$ is related to its PDF via

$$ e^{\varphi(y)} = \int d\delta_{m,T} e^{y\delta_{m,T}} p(\delta_{m,T}), \quad (4.3) $$

which is the Laplace transform of the PDF. Hence, if $\varphi(y)$ is a known analytic function, then $p(\delta_{m,T})$ can be computed by the inverse Laplace transform

$$ p(\delta_{m,T}) = \int_{-\infty}^{\infty} \frac{dy}{2\pi} e^{-iy\delta_{m,T} + \varphi(iy)}. \quad (4.4) $$

This integral in the complex plain is most efficiently evaluated along the path of steepest descent (cf. [72], especially their Appendix B). The cumulants of $\delta_{m,T}(\hat{\mathbf{n}})$ can be approximated as (cf. Bernardeau and Valageas [49])

$$ \langle \delta_{m,T}(\hat{\mathbf{n}})^k \rangle_c \approx \int dw \frac{\langle \{\delta_{cy,\theta,w,L}(w\hat{\mathbf{n}}, w)q_l(w)L\}^k \rangle_c}{L}. \quad (4.5) $$

Here $q_l$ is the line-of-sight density of tracer galaxies defined in Eq. (2.2) and $\delta_{cy,R,L}$ is the average of $\delta_{m,3D}$ over a cylinder of length $L$ and radius $R$, where $L$ has to be chosen such that correlations of $\delta_{m,3D}$ over a distance $L$ vanish for all practical purposes. Equation (4.5) employs a small angle approximation and a Limber-like approximation [76]. This means it assumes that any n-point correlation function between density contrast at different redshifts $z_i$, $i = 1...n$, varies much more quickly with the redshift differences $\Delta z_{ij}$ than the projection

kernel $q_l$. Comparing (4.2) and (4.5) we see that the CGF of $\delta_{m,T}$ can be computed in terms of the CGF of $\delta_{cy,R,L}$ as

$$ \varphi(y) \approx \int dw \frac{\varphi_{cy,\theta,w,L}(q_l(w)Ly, w)}{L}. \quad (4.6) $$

Hence, we have reduced the task of computing $p(\delta_{m,T})$ to the task of computing the connected moments of matter contrast in a long 3D cylinder.

To proceed we consider two different ansatzes. The first is to approximate $p(\delta_{m,T})$ by a log-normal distribution which is fixed by the first three connected moments of $\delta_{m,T}(\hat{\mathbf{n}})$. The second ansatz is to compute $\varphi_{cy,\theta,w,L}(y)$ as a whole in a way similar to the one of Bernardeau [48] and Valageas [50] for the matter contrast in a 3-dimensional sphere. Using (4.6) we can then attempt to solve the integral in (4.4) directly. We present details of both approaches in the following subsections.

## 1. Log-normal approximation for the PDF

Instead of computing the complete cumulant generating function of $\delta_{m,T}$ via Eq. (4.6) we start with an approach that only requires knowledge of the second and third cumulant, i.e., the moments $\langle \delta_{m,T}^2 \rangle_c$ and $\langle \delta_{m,T}^3 \rangle_c$ (implicitly we also assume $\langle \delta_{m,T} \rangle_c = 0$ for the first cumulant). To do so, we approximate $\delta_{m,T}$ as a shifted log-normal random variable [43,44]. In this case the PDF of $\delta_{m,T}$ is given by

$$ p(\delta_{m,T}) = \frac{\exp\left(-\frac{\left[\ln(\delta_{m,T}/\delta_0 - 1) + \sigma^2/2\right]^2}{2\sigma^2}\right)}{\sqrt{2\pi}\sigma(\delta_{m,T} - \delta_0)} \quad (4.7) $$

if $\delta_{m,T} > \delta_0$ and $p(\delta_{m,T}) = 0$ otherwise. The expectation value of this PDF is zero, as is appropriate. The variance and skewness of $\delta_{m,T}$ are given in terms of the parameters $\delta_0$ and $\sigma$ by [44]

$$ \begin{aligned} \langle \delta_{m,T}^2 \rangle_c &= \delta_0^2 (e^{\sigma^2} - 1) \\ \langle \delta_{m,T}^3 \rangle_c &= \frac{3}{\delta_0} \langle \delta_{m,T}^2 \rangle_c^2 + \frac{1}{\delta_0^3} \langle \delta_{m,T}^2 \rangle_c^3. \end{aligned} \quad (4.8) $$

The ansatz of a log-normal PDF has been shown to be consistent with early DES data by Clerkin et al. [77]. For the 3-dimensional density contrast $\delta_{m,3D}$ this can be reasonably motivated from theory by observing that at leading order in perturbation theory the skewness of $\delta_{m,3D}$ scales as

$$ \langle \hat{\delta}_{m,3D}^3 \rangle_c \sim \langle \hat{\delta}_{m,3D}^2 \rangle_c^2 \quad (4.9) $$

with corrections of the order $\langle \hat{\delta}_{m,3D}^2 \rangle_c^3$. This is exactly the scaling obeyed by the log-normal distribution and choosing the parameter $\delta_0$ appropriately allows one to exactly reproduce the scaling coefficients predicted by perturbation theory.
At least for a power law power spectrum, the same kind of scaling is observed also for 2-dimensional, projected versions of the density field. Hence, one of the ansatzes used in this paper to compute the PDF of $\delta_{m,T}$ is to derive its variance and skewness as described in Appendix B and fix $\delta_0$ and $\sigma$ by demanding that the PDF in (4.7) has the same second and third moments.

## 2. Tree level computation of $\varphi_{cy,\theta,w,L}(y,w)$ in the cylindrical collapse model

Let us first consider the cumulant generating function $\varphi_{3D}(y, \tau)$ of the 3-dimensional density contrast $\delta_{3D}(\mathbf{x}, \tau)$. To compute $\varphi_{3D}$ at tree-level in perturbation theory it is sufficient to assume spherical symmetry around a particular point $\mathbf{x}$ [see, e.g., [50]]. Doing so, we can directly compute $\delta_{3D}(\mathbf{x}, \tau)$ as a function of the linear density contrast $\delta_{3D,lin}(\mathbf{x}, \tau)$ by means of the spherical collapse model [50,78], i.e.

$$ \delta_{3D}(\mathbf{x}, \tau) = F[\delta_{3D,lin}(\mathbf{x}, \tau), \tau] \quad (4.10) $$

where $F$ is determined by one of the differential equations presented in appendix A. Hence, using the assumption that the linear density contrast has a Gaussian distribution with variance $\sigma_{3D,lin}^2(\tau)$ we can express the cumulant generating function as [cf. Eq. (4.3)]

$$ \begin{aligned} \exp\{\varphi_{3D}(y, \tau)\} &= \int d\delta_{3D} e^{y\delta_{3D}} p(\delta_{3D}, \tau) \\ &= \int \frac{d\delta_{3D,lin}}{\sqrt{2\pi\sigma_{3D,lin}^2(\tau)}} \\ & \quad \times \exp\left(yF[\delta_{3D,lin}, \tau] - \frac{\delta_{3D,lin}^2}{2\sigma_{3D,lin}^2(\tau)}\right), \end{aligned} \quad (4.11) $$

where in the last step we simply performed a change of variables from $\delta_{3D}$ to $\delta_{3D,lin}$. We now employ Laplace's method, which states that a function $f(x)$ which strongly peaks around $x_0$ fulfils

$$ \int dx e^{f(x)} \approx \sqrt{\frac{2\pi}{|f''(x_0)|}} e^{f(x_0)}. \quad (4.12) $$

This way we can approximate the last line of (4.11) as

$$ e^{\varphi_{3D}(y, \tau)} \approx \sqrt{\frac{1}{|1 - yF''[\delta^*, \tau]\sigma_{3D,lin}^2(\tau)|}} \\ \times \exp\left(yF[\delta^*, \tau] - \frac{\delta^{*2}}{2\sigma_{3D,lin}^2(\tau)}\right), \quad (4.13) $$

where $\delta^*$ is the linear density contrast that maximizes the exponent in (4.11) and ' denotes derivation with respect to $\delta$. In the quasilinear limit of $\sigma_{3D,lin}^2 \to 0$ this gives

$$ \varphi_{3D}(y, \tau) \approx yF[\delta^*, \tau] - \frac{\delta^{*2}}{2\sigma_{3D,lin}^2(\tau)}, \quad (4.14) $$

where $\delta^*$ has to be determined by the implicit equation

$$ \frac{\delta^*}{\sigma_{3D,lin}^2(\tau)} = yF'[\delta^*, \tau]. \quad (4.15) $$

It should be noted that Eqs. (4.14) and (4.15) reproduce exactly the tree-level results for the cumulant generating function [cf. [48,50,51,73]]. As described in Bernardeau et al. [51] the coefficients

$$ S_n = \frac{\langle \delta^n \rangle_c}{\langle \delta^2 \rangle_c^{n-1}} \quad (4.16) $$

are given quite accurately by the lowest order of perturbation theory. Hence, using HALOFIT [45,46] and an analytic transfer function [47] to compute the nonlinear matter power spectrum, we can compute the nonlinear variance $\langle \delta^2 \rangle_{c,non lin}$ and then rescale the leading order CGF to its nonlinear version via

$$ \varphi_{non lin}(y, \tau) = \frac{\langle \delta^2 \rangle_{c,lin}}{\langle \delta^2 \rangle_{c,non lin}} \varphi_{lead} \left( y \frac{\langle \delta^2 \rangle_{c,non lin}}{\langle \delta^2 \rangle_{c,lin}}, \tau \right). \quad (4.17) $$

To perform the projection integral in Eq. (4.6) we need to know the cumulant generating function of the density contrast in a long cylinder of radius R and length $L \gg R$, $\delta_{cy,R,L}$. Bernardeau [48] [1994; see also [50], and the other references above] have generalized Eq. (4.14) to the case of matter density in a spherical aperture. Their results can easily be transferred to cylindrical apertures, yielding

$$ \varphi_{cy,R,L}(y, \tau) \approx yF_{cyl}[\delta^*, \tau] - \frac{\delta^{*2}}{2\sigma_R^2 \sqrt{1+F_{cyl}[\delta^*], L, lin}(\tau)}, \quad (4.18) $$

where $\sigma_{R,L,lin}^2$ is the variance of linear density contrast in a cylinder, $F_{cyl}$ is now determined by cylindrical collapse and $\delta^*$ has to be determined by the implicit equation

$$ \frac{1}{2} \frac{d}{d\delta^*} \frac{\delta^{*2}}{\sigma_R^2 \sqrt{1+F_{cyl}[\delta^*], L, lin}(\tau)} = yF'[\delta^*, \tau]. \quad (4.19) $$

Using Eq. (4.17) we can again rescale this leading order result for the generating function to its nonlinear counterpart.

The validity of Eq. (4.17) is ultimately limiting the accuracy of our model for the distribution of the density contrast inside the aperture radius $\theta_T$, $p(\delta_{m,T})$. In Fig. 5 we compare the PDF measured in the Buzzard simulations to both the log-normal model and the full CGF computation of the PDF for aperture radii $\theta_T = 10', 20', 30'$. For both
$\theta_T = 20'$ and $30'$ our model PDF's and the Buzzard simulations agree within DES-Y1 cosmic variance as can be seen in the residuals shown in the lower panels of each plot. Also, the difference between log-normal approximation and full CGF computation is completely negligible. To investigate the agreement of Buzzard and our models more quantitatively, we also compare the variance and skewness of each PDF. In Fig. 6 we show the ratio of these moments as found in Buzzard to our predictions. For $\theta_T = 10'$, the density field in Buzzard seems to have a significantly higher skewness than predicted by our model. We attribute this indeed to the failing of the quasilinear rescaling, Eq. (4.17). For the other radii the skewness values agree to within 2%–3%. A similar relative agreement is achieved for the variance of the distributions. At our fiducial radius $\theta_T = 20'$ the variances of Buzzard and our model differ by about 2.4%. This corresponds to a disagreement in the fluctuation amplitude $\sigma_8$ of about 1.2%.

For comparison, we also show a Gaussian model for the density PDF in Fig. 5, using the same variance as for the other PDF models. It can clearly be seen, that $p(\delta_{m,T})$ cannot be well described by a Gaussian distribution for any of the considered smoothing radii.

## B. Convergence profile around lines-of-sight of fix density contrast $\delta_{m,T}$

We now want to know the average lensing convergence $\kappa_{<\theta}$ inside a radius $\theta$ around a line-of-sight with a given value of $\delta_{m,T}$. By means of Eq. (2.6) this can be turned into a prediction of the density split lensing signal. We start by looking at the joint moment generating function of $\kappa_{<\theta}$ and $\delta_{m,T}$,

<figure><img src="image_8.png" /><figcaption>FIG. 6. This figure shows the ratio of moments of $\delta_{m,T}$ measured in a Buzzard all-sky density map to the moments predicted by our model as a function of $\theta_T$. The error bars represent the standard deviations of the moments found in a set of 14 DES Y1-sized patches in the Buzzard map. At our fiducial radius $\theta_T = 20'$ we find a ~2.4% deviation between the variance of $\delta_{m,T}$ measured in Buzzard and that predicted by our model. This would correspond to a ~1.2% deviation in $\sigma_8$.</figcaption></figure>

$$ \psi_{\theta}(y, z) = \sum_{k,l} \frac{\langle \delta_{m,T}^k \kappa_{<\theta}^l \rangle}{k!l!} y^k z^l, \quad (4.20) $$

and their joint cumulant generating function defined by

$$ \begin{aligned} \varphi_{\theta}(y, z) &= \ln \psi_{\theta}(y, z) \\ &\equiv \sum_{k,l} \frac{\langle \delta_{m,T}^k \kappa_{<\theta}^l \rangle_c}{k!l!} y^k z^l. \end{aligned} \quad (4.21) $$

Using a Limber-like approximation similar to the one employed in Eq. (4.6), one can write $\varphi_{\theta}(y, z)$ as a line-of-sight projection of the CGF of matter density contrasts that are smoothed over concentrical cylindrical apertures with length $L$ and radii $R_1, R_2, \varphi_{cyl, R_1, R_2, L}(y, z, w)$:

$$ \varphi_{\theta}(y, z) \approx \int dw \frac{\varphi_{cyl, \theta_T, w, \theta_w, L}(q_l(w)Ly, W_s(w)Lz, w)}{L}. \quad (4.22) $$

Here, $q_l(w)$ is again the line-of-sight density of lens galaxies and $W_s(w)$ is the lensing efficiency defined in Eq. (2.5).

The joint PDF of $\kappa_{<\theta}$ and $\delta_{m,T}$ is related to $\varphi_{\theta}$ via

$$ \begin{aligned} p(\delta_{m,T} = s, \kappa_{<\theta} = r) \\ &= \int_{-\infty}^{\infty} \frac{dy dz}{(2\pi)^2} e^{-iys - izr + \varphi_{\theta}(iy, iz)}. \end{aligned} \quad (4.23) $$

The expectation value of $\kappa_{<\theta}$ given a certain value of $\delta_{m,T}$ is then given by

$$ \begin{aligned} \langle \kappa_{<\theta} | \delta_{m,T} = s \rangle &= \int dr p(\delta_{m,T} = s, \kappa_{<\theta} = r) \\ &= \int_{-\infty}^{\infty} \frac{dy dz}{(2\pi)^2} e^{-iys + \varphi_{\theta}(iy, iz)} \int dr r e^{-izr} \\ &= \int_{-\infty}^{\infty} \frac{dy dz}{2\pi} e^{-iys + \varphi_{\theta}(iy, iz)} i \frac{d}{dz} \delta_D(z) \\ &= \int_{-\infty}^{\infty} \frac{dy}{2\pi} e^{-iys + \varphi_{m,T}(iy)} G_{\theta}(iy), \end{aligned} \quad (4.24) $$

with

$$ \begin{aligned} G_{\theta}(y) &= \frac{d}{dz} \varphi_{\theta}(y, z) \Big|_{z=0} \\ &= \sum_k \frac{\langle \delta_{m,T}^k \kappa_{<\theta} \rangle_c}{k!} y^k. \end{aligned} \quad (4.25) $$

Using equation (4.22) we can express $G_{\theta}(y)$ by the corresponding cylindrical, 3-dimensional quantity:
$$G_{\theta}(y) \approx \int dw W_{s}(w) G_{\text{cyl}, \theta_{T} w, \theta_{w} L}(q_{l}(w) L y, w). \quad (4.26)$$

We again pursue two ansatzes for computing $G_{\theta}(y)$: one involving a log-normal model for the joint cumulants of $\kappa_{<\theta}$ and $\delta_{m,T}$ and one involving the model of cylindrical collapse to compute a leading order perturbation theory prediction for $G_{\text{cyl},\theta_T w, \theta_w L}(y, w)$. We are detailing these ansatzes in the following subsections.

## 1. Log-normal model for the joint cumulants of $\kappa_{<\theta}$ and $\delta_{m,T}$

From Eq. (4.25) one can see that only joint cumulants of the form

$$\langle \delta_{m,T}^k \kappa_{<\theta} \rangle_c \quad (4.27)$$

enter the convergence profile around a given density contrast. Hence, in a spirit similar to Sec. IVA 1 we only compute the leading order cumulants $\langle \delta_{m,T} \kappa_{<\theta} \rangle_c$ and $\langle \delta_{m,T}^2 \kappa_{<\theta} \rangle_c$ as described in Appendix B and use these moments to fix a joint log-normal PDF for $\kappa_{<\theta}$ and $\delta_{m,T}$. Note that this is indeed not assuming, that $\kappa_{<\theta}$ is a log-normal random variable. It rather assumes that

$$\kappa_{<\theta} = \kappa_{\text{log-normal}} + \kappa_{\text{uncorr}}, \quad (4.28)$$

where $\kappa_{\text{log-normal}}$ is log-normal and $\kappa_{\text{uncorr}}$ is an unspecified random variable that is uncorrelated with $\delta_{m,T}$. Only $\kappa_{\text{log-normal}}$ will actually contribute to the density split lensing signal and we can model the expectation value $\langle \kappa_{<\theta} | \delta_{m,T} = s \rangle$ by the following relation holding for two joint log-normal variables:

$$\frac{\langle \kappa_{<\theta} | \delta_{m,T} = s \rangle}{\kappa_0} = \exp \left( \frac{C(2 \log(1 + s/\delta_0) + V - C)}{2V} \right) - 1 \quad (4.29)$$

with

$$\begin{align*}
V &= \log \left( 1.0 + \frac{\langle \delta_{m,T}^2 \rangle_c}{\delta_0^2} \right) \\
C &= \log \left( 1.0 + \frac{\langle \delta_{m,T} \kappa_{<\theta} \rangle_c}{\delta_0 \kappa_0} \right) \\
\kappa_0 &= \frac{\langle \delta_{m,T} \kappa_{<\theta} \rangle_c^2 e^V}{\langle \delta_{m,T}^2 \kappa_{<\theta} \rangle_c - 2 \langle \delta_{m,T} \kappa_{<\theta} \rangle_c \langle \delta_{m,T}^2 \rangle_c / \delta_0} \quad (4.30)
\end{align*}$$

and $\delta_0$ determined as described in Sec. IVA 1.

It should be noted that the log-normal parameter $\kappa_0$ which we use to approximate the contribution of $\kappa_{<\theta}$ to the lensing signal is dependent on the smoothing scale $\theta$.

This indicates even further, that we do indeed not approximate $\kappa$ as a log-normal field.

## 2. Tree level computation of $G_{\text{cyl},R_1,R_2,L}(y,w)$ in the cylindrical collapse model

For convenience we will shorten the notation of Sec. IVA 2 by defining

$$\begin{align*}
G(y) &= G_{\text{cyl},R_1,R_2,L}(y, w) \\
\varphi(y, z) &= \varphi_{\text{cyl},R_1,R_2,L}(y, z, w) \\
F[\delta] &= F_{\text{cyl}}[\delta, \tau] \\
\sigma_R &= \sigma_{R,L,\text{lin}} \\
R_{i,L}(\delta) &= R_i \sqrt{1 + F[\delta]}, \quad i = 1, 2. \quad (4.31)
\end{align*}$$

The joint cumulant generating function of density contrast in concentric cylinders is then [in complete analogy to Eq. (4.18); see also [49], who present very similar calculations] given by

$$\varphi(y, z) \approx yF[\delta_1^*] + zF[\delta_2^*] - \frac{1}{2} \sum_{i,j} \delta_i^* \delta_j^* (C^{-1})_{ij}, \quad (4.32)$$

where the elements of the matrix $C$ are given by $C_{11} = \sigma_{R_1,L}^2(\delta_1^*)$, $C = \sigma_{R_2,L}^2(\delta_2^*)$ and $C_{12} = C_{21}$ is the linear covariance of density contrasts in concentric cylinders of radii $R_{1,L}(\delta_1^*)$ and $R_{2,L}(\delta_2^*)$. This time the critical linear density contrasts $\delta_1^*$ and $\delta_2^*$ are given by the implicit equations

$$\frac{1}{2} \frac{\partial}{\partial \delta_1^*} \sum_{i,j} \delta_i^* \delta_j^* (C^{-1})_{ij} = yF'[\delta_1^*] \quad (4.33)$$

$$\frac{1}{2} \frac{\partial}{\partial \delta_2^*} \sum_{i,j} \delta_i^* \delta_j^* (C^{-1})_{ij} = zF'[\delta_2^*]. \quad (4.34)$$

Note that these conditions force each $\delta_i^*$ to be a function of both $y$ and $z$, i.e., $\delta_i^* = \delta_i^*(y, z)$.

To predict the convergence profile around apertures of a given density contrast $\delta_{m,T}$ by means of Eqs. (4.24) and (4.26) we are interested in computing the function

$$G(y) = \left. \frac{\partial}{\partial z} \varphi(y, z) \right|_{z=0}. \quad (4.35)$$

Using the conditions (4.33) and (4.34) one can see right away that

$$G(y) = F[\delta_2^*(y, 0)]. \quad (4.36)$$

Furthermore, for $z=0$ Eqs. (4.33) and (4.34) can be simplified to
$$ \frac{1}{2} \frac{d}{d\delta_1^*} \frac{\delta_1^{*2}}{C_{11}} = yF'[\delta_1^*] \quad (4.37) $$

$$ \delta_2^* = \frac{C_{12}}{C_{11}} \delta_1^*. \quad (4.38) $$

This way we obtain a solution for $G(y)$ at leading order in perturbation theory. In Appendix B we argue that for $R_2 \ge R_1$ the cumulants $\langle \delta_{R_1}^k \delta_{R_2} \rangle_c$ approximately follow the scaling relation

$$ \langle \delta_{R_1}^k \delta_{R_2} \rangle_c \sim \langle \delta_{R_1} \delta_{R_2} \rangle_c \langle \delta_{R_1}^2 \rangle_c^{k-1}. \quad (4.39) $$

This can be used to correct the tree-level approximation of $G(y)$ for the nonlinear evolution of the power spectrum. To do so, we first determine the proportionality factors of the relation (4.39) at leading order by fitting a polynomial in $y$ to the function $G(y)$ and extracting the cumulants $\langle \delta_{R_1}^k \delta_{R_2} \rangle_c$ from the polynomial coefficients. In practice, we do this with a polynomial of degree 10, but already a polynomial of degree 5 gives almost identical results.<fnref n="4" /> Then, we use the revised HALOFIT of Takahashi et al. [46] to compute a late-time version of the right-hand side of (4.39). This, together with the tree-level proportionality factors determined before, yields a nonlinear approximation of the polynomial coefficients representing $G(y)$. We use those to re-compute $G(y)$ and then carry out the projection integral in Eq. (4.26).

Our rescaling of the coefficient corresponding to the cumulant $\langle \delta_{R_1}^2 \delta_{R_2} \rangle_c$ is in fact more complicated than described here, cf. Appendix B 5. But we find our prediction of the density split lensing signal to be insensitive to the details of the rescaling procedure.

### C. Shot-noise, stochasticity, and counts-in-cells

We now want to model the conditional probability $P(N_T|\delta_{m,T})$ of finding $N_T$ galaxies in an angular radius of $\theta_T$, when the projected density contrast in that radius is $\delta_{m,T}$. This is the third ingredient of the framework described in Sec. II and completes our modeling of the density split lensing signal as well as the counts-in-cells histogram.

To analyze the relation of $N_T$ and $\delta_{m,T}$ in a systematic way, let us introduce the auxiliary field $\delta_{g,T}$. We assume that $\delta_{g,T}(\hat{n})$ is a smooth field in the sky and that $N_T$ is a Poissonian tracer of this field. This means we will assume that

$$ P(N_T = N | \delta_{g,T}) = \frac{[\bar{N}(1 + \delta_{g,T})]^N}{N!} e^{-\bar{N}(1 + \delta_{g,T})}, \quad (4.40) $$

where $\bar{N} \equiv \langle N_T \rangle$. A consequence of this assumption is that the expectation value of $N_T$ for fixed $\delta_{g,T}$ is given by

$$ \langle N_T | \delta_{g,T} \rangle = \bar{N}(1 + \delta_{g,T}) \quad (4.41) $$

and that the variance of $N_T$ for fixed $\delta_{g,T}$ fulfils

$$ \frac{\mathrm{Var}[N_T | \delta_{g,T}]}{\langle N_T | \delta_{g,T} \rangle} = 1. \quad (4.42) $$

To connect the galaxy field to the lensing convergence we however need to know the relation between $N_T$ and $\delta_{m,T}$. Assuming a generic joint PDF of $\delta_{m,T}$ and $\delta_{g,T}$ we can write the expectation values of $N_T$ for fixed $\delta_{m,T}$ as

$$ \begin{aligned} \langle N_T | \delta_{m,T} \rangle &= \int d\delta_{g,T} p(\delta_{g,T} | \delta_{m,T}) \langle N_T | \delta_{g,T} \rangle \\ &= \bar{N}(1 + \langle \delta_{g,T} | \delta_{m,T} \rangle). \end{aligned} \quad (4.43) $$

Also, it can be shown that the variance of $N_T$ for a fixed value of $\delta_{m,T}$ is given by

$$ \mathrm{Var}[N_T | \delta_{m,T}] = \langle N_T | \delta_{m,T} \rangle + \bar{N}^2 \mathrm{Var}[\delta_{g,T} | \delta_{m,T}]. \quad (4.44) $$

From Eq. (4.44) we can see that the distribution of $N_T$ given $\delta_{m,T}$ can only be a Poisson distribution if $\mathrm{Var}[\delta_{g,T} | \delta_{m,T}] \equiv 0$. This is the simplest model for the connection of $N_T$ and $\delta_{m,T}$ that we test in this work and in Gruen et al. [39]. If $\mathrm{Var}[\delta_{g,T} | \delta_{m,T}] \ne 0$ we will say that there is a stochasticity between the galaxy field and the matter density field, and we cannot assume a Poisson distribution for $P(N_T|\delta_{m,T})$. We note that “stochasticity” in this context could arise from a nonlinear biasing relationship between $\delta_{g,T}$ and $\delta_{m,T}$, including, e.g., a dependence on higher powers of $\delta_{m,T}$ or effects from halo exclusion [79], or from physical stochasticity in galaxy formation.

We explore two ways to account for a possible stochasticity (see also Dekel and Lahav [80], who have discussed similar concepts). In our first approach we introduce a free parameter to our model—a Pearson correlation coefficient $r \ne 1$ between the random fields $\delta_{g,T}$ and $\delta_{m,T}$. Within our log-normal framework we show that this automatically leads to a $\delta_{m,T}$-dependence of the ratio in Eq. (4.46). We explain the details of this in Sec. IV C 1).

In our second approach we employ a generalized Poisson distribution for $P(N_T|\delta_{m,T})$ that allows for

$$ \frac{\mathrm{Var}[N_T | \delta_{m,T}]}{\langle N_T | \delta_{m,T} \rangle} \ne 1. \quad (4.45) $$

In this approach we introduce 2 parameters, $\alpha_0$ and $\alpha_1$, to our model such that

<footnote n="4">⁴The coefficients of linear and quadratic order in $y$ are always obtained from the exact perturbation theory computation of Appendix B.</footnote>
$$ \frac{\mathrm{Var}[N_T | \delta_{m,T}]}{\langle N_T | \delta_{m,T} \rangle} \approx \alpha_0 + \alpha_1 \delta_{m,T}. \quad (4.46) $$

The details of this are explained in Sec. IV C 2.

Both of our approaches match our simulated data equally well (cf. Fig. 8). This means that, for the galaxies in these realizations, the model based on the correlation coefficient *r* is a sufficient description. It will thus be the fiducial model in this paper, used in all figures unless otherwise noted. In Gruen et al. [39] we will nevertheless apply both this and the two-parametric model to account for the possibility that the shot-noise of real galaxies behaves in a more complicated way than that of our simulated galaxies.

# 1. Shot-noise model 1: Correlation *r* ≠ 1 between galaxy density and matter density

In our fiducial model of $P(N_T|\delta_{m,T})$ we approximate the joint distribution of both $\delta_{m,T}$ and $\delta_{g,T}$ with a joint log-normal distribution (cf. Eq. (4.7) and Hilbert et al. [44] for properties of joint log-normal distributions). The joint PDF of two log-normal random variables is characterized by 5 parameters, e.g., by the variance and skewness of each variable and the covariance between the two variables.

In our case, we compute the variance and skewness of $\delta_{m,T}$ as described in Sec. IVA 1 and set the variance and skewness of $\delta_{g,T}$ to

$$ \begin{aligned} \langle \delta_{g,T}^2 \rangle_c &= b^2 \langle \delta_{m,T}^2 \rangle_c \\ \langle \delta_{g,T}^3 \rangle_c &= b^3 \langle \delta_{m,T}^3 \rangle_c \end{aligned} \quad (4.47) $$

where the galaxy bias *b* is a free parameter. The covariance of $\delta_{m,T}$ and $\delta_{g,T}$ is parametrized by their correlation coefficient

$$ r = \frac{\langle \delta_{g,T} \delta_{m,T} \rangle_c}{\sqrt{\langle \delta_{g,T}^2 \rangle_c \langle \delta_{m,T}^2 \rangle_c}} = \frac{\langle \delta_{g,T} \delta_{m,T} \rangle_c}{b \langle \delta_{m,T}^2 \rangle_c}, \quad (4.48) $$

i.e.,

$$ \langle \delta_{g,T} \delta_{m,T} \rangle_c = rb \langle \delta_{m,T}^2 \rangle_c. \quad (4.49) $$

The log-normal model for the joint PDF of $\delta_{m,T}$ and $\delta_{g,T}$ now allows us to compute the variance of galaxy counts as a function of $\delta_{m,T}$ and more generally to compute $P(N_T|\delta_{m,T})$. We present the necessary derivations in detail in Appendix D.

In our data analysis we consider *b* and *r* as free parameters. The only restrictions we impose on them are

$$ 0 < b, \quad 0 \le r \le 1. \quad (4.50) $$

To test how accurately this model describes the behavior of our mock REDMAGIC catalogs based on the Buzzard

N-body simulations we nevertheless want to determine what values of *b* and *r* are underlying our simulations. To this end, we generate HEALPIX maps of $\delta_{m,T}$ with different top-hat aperture radii $\theta_T = 10, 20, 30$ arcmin, based on particle count maps at resolution $N_{side} = 8192$ in slices of comoving $50 h^{-1}$ Mpc thickness. We co-add these maps to reproduce a redshift range close to that of our fiducial analysis, $z = 0.2100...0.4453$. We then select REDMAGIC galaxies with true redshift in this range and determine their counts around the same HEALPIX pixel centers and within the same aperture radii. The REDMAGIC mock catalogs have a complex mask similar to that of real DES data, which adds complication because the fraction of masked area in each aperture must be equal in order to meaningfully sort lines of sight by galaxy count. To this end, we convert all counts to a masking fraction of 20 percent of area within the aperture radius using the stochastic method of Gruen et al. [[39], their Sec. II.A]. This leaves us with simulated 2D maps of $\delta_{m,T}$ and $N_T$ within a DES-Y1 shaped mask.

We can then measure the variances of these maps, $\mathrm{Var}(\delta_{m,T})$ and $\mathrm{Var}(N_T)$, as well as their covariance $\mathrm{Cov}(N_T, \delta_{m,T})$. These fulfill the relations

$$ \mathrm{Var}(N_T) = \bar{N} + \bar{N}^2 b^2 \mathrm{Var}(\delta_{m,T}). \quad (4.51) $$

and

$$ \mathrm{Cov}(N_T, \delta_{m,T}) = \bar{N}br \mathrm{Var}(\delta_{m,T}) \quad (4.52) $$

which then fixes *b* and *r*. The values determined in this way are shown in table II.

We now need to check whether these value for *b* and *r* together with our assumption of a log-normal PDF for $\delta_{m,T}$ and $\delta_{g,T}$ describe the properties of our tracer galaxies well. Using our simulated maps of $\delta_{m,T}$ and $N_T$ we can measure the expectation value

$$ \langle \delta_{g,T} | \delta_{m,T} \rangle = \frac{\langle N_T | \delta_{m,T} \rangle}{\bar{N}} - 1 \quad (4.53) $$

as a function of $\delta_{m,T}$. Within the log-normal model (cf. Appendix D for the relevant formulas) this is very well approximated by

<table><caption>TABLE II. Best-fit values galaxy bias and correlation coefficients of our simulated tracer galaxies within the model presented in Sec. IV C 1. Error bars are estimated from a jackknife approach.</caption><thead><tr><th>Smoothing<br>Scale [arcmin]</th><th>b</th><th>r</th></tr></thead><tbody><tr><td>10</td><td>1.644 ± 0.008</td><td>0.938 ± 0.001</td></tr><tr><td>20</td><td>1.618 ± 0.008</td><td>0.956 ± 0.001</td></tr><tr><td>30</td><td>1.603 ± 0.008</td><td>0.961 ± 0.001</td></tr></tbody></table>
<figure><img src="image_9.png" /><figcaption>FIG. 7. Average galaxy overdensity $\delta_{g,T}$ as a function of matter overdensity, $\delta_{m,T}$, for our simulated maps at different smoothing scales: 10 arcmin [left], 20 arcmin [middle], and 30 arcmin [right]. Solid lines show a linear bias model with the bias parameters obtained from maximizing the likelihood in Eq. (4.63) and the residual between the two are shown in the bottom panels. Note that the coefficient of linearity found with (4.63) ($\approx 1.54$) is almost identical to the value of the product $b \cdot r$ determined with equation (4.54). To indicate the range of $\delta_{m,T}$ that is relevant to our computation, we also show the density PDFs of Fig. 5 as shaded regions. The error bars were derived from a jackknife approach.</figcaption></figure>

$$ \langle \delta_{g,T} | \delta_{m,T} \rangle \approx rb\delta_{m,T} \qquad (4.54) $$

which becomes exact for Gaussian random variables. In Fig. 7 we compare measurements of $\langle \delta_{g,T} | \delta_{m,T} \rangle$ with the different smoothing radii $\theta_T = 10, 20, 30$ arcmin to the prediction of the log-normal model. We find that in our simulations $\langle \delta_{g,T} | \delta_{m,T} \rangle$ is consistent with a linear relation in $\delta_{m,T}$. Interestingly, the scale dependence of $b$ and $r$ we find in Table II almost perfectly cancels to give a scale independent proportionality coefficient

$$ rb \approx 1.54. \qquad (4.55) $$

Next, we also measure the variance of galaxy counts $N_T$ as a function of $\delta_{m,T}$ in our simulated maps and compare to the prediction of the log-normal model. In Fig. 8 we indeed find that

$$ \frac{\mathrm{Var}[N_T | \delta_{m,T}]}{\langle N_T | \delta_{m,T} \rangle} \neq 1 \qquad (4.56) $$

and that the $\delta_{m,T}$-dependence is very well described by the log-normal model and the values of $b$ and $r$ we determined before. Finally, in Fig. 9 we show the residuals between our baseline prediction of the counts-in-cells histogram with $\theta_T = 20$ arcmin<fnref n="5" /> and the average of measurements in 4 Buzzard realizations of DES Y1 data (cf. Fig. 1). The residuals are well contained within DES Y1 errorbars.

## 2. Shot-noise model 2: Parametric model for super-Poissonianity

Our model of shot-noise based on galaxy bias $b$ and galaxy matter correlation coefficient $r$ describes our simulated tracer

<figure><img src="image_10.png" /><figcaption>FIG. 8. Ratio of the variance in the galaxy count distribution to average galaxy counts in our simulated maps as a function of matter overdensity. Differently coloured solid lines show the result for each smoothing scale. The dashed black lines show the predictions of the 2-parametric shot-noise model described in Sec. IV C 2. The blue dotted lines show the corresponding prediction of the alternative, 1-parametric model described in Sec. IV C 1. The horizontal solid line shows the expectation if shot noise was purely Poissonian. The coloured regions show the 95% confidence limits derived from jackknife resampling.</figcaption></figure>

<footnote n="5">⁵This is the smoothing radius used in our data analysis [39].</footnote>
<figure><img src="image_11.png" /><figcaption>FIG. 9. For the fiducial smoothing radius of our data analysis presented in Gruen et al. [39], $\theta_T = 20$ arcmin, we show the residuals between our baseline prediction of the counts-in-cells histogram and the average of measurements in 4 Buzzard realizations of DES Y1 data. Blue error bars represent the uncertainties we expect for DES Y1 while green error bars show the actual uncertainties of the mean measurement from our simulations.</figcaption></figure>

catalog well. But it contains the arbitrary assumption that both the variance and the skewness of $\delta_{m,T}$ and $\delta_{g,T}$ are related through the bias parameter $b$ (cf. Eq. (4.47)). To account for the possibility that real galaxies might behave in a more complicated way, we also consider a more flexible model of the conditional distribution $P(N_T|\delta_{m,T})$.

Gruen et al. [38] assumed that there is no stochasticity in the relation of $\delta_{m,T}$ and $\delta_{g,T}$ and that galaxies trace the matter density with a linear bias and Poissonian shot noise. This means they set

$$P(N_T = N | \delta_{m,T} = s) = \exp(N \ln[\bar{N}(1 + \tilde{b}s)] - [\bar{N}(1 + \tilde{b}s)] - \ln \Gamma(N + 1)), \quad (4.57)$$

where $\tilde{b}$ is the galaxy bias, and where we now use a generalizable definition of the Poisson distribution based on the Gamma function $\Gamma$, for reasons that will appear later. The galaxy bias $\tilde{b}$ defined this way is not identical to the definition in our fiducial model. We now rather have

$$\delta_{g,T} \equiv \tilde{b}\delta_{m,T}. \quad (4.58)$$

We test this in our simulated maps of $N_T$ and $\delta_{m,T}$ by fitting a linear relation to the mean smoothed galaxy contrast as a function of dark matter contrast that was shown in Fig. 7. We indeed find that this linear biasing model describes the

simulations very well and that $\tilde{b} \approx br$ as expected from our arguments in Sec. IV C 1.

The model used by Gruen et al. [38] however predicts that

$$\frac{\mathrm{Var}[N_T | \delta_{m,T}]}{\langle N_T | \delta_{m,T} \rangle} = 1 \quad (4.59)$$

which is not what we find in Fig. 8. To account for the deviations we observe from pure Poissonian shot-noise, we hence model the distribution of $N_T$ given $\delta_{m,T}$ as

$$P(N_T = N | \delta_{m,T} = s) \sim \mathcal{N} \exp \left\{ \frac{N}{\alpha} \ln \left[ \frac{\bar{N}_T}{\alpha} (1 + b s) \right] - \ln \Gamma \left( \frac{N+1}{\alpha} \right) - \frac{\bar{N}_T}{\alpha} (1 + b s) \right\}, \quad (4.60)$$

where the parameter $\alpha > 0$ generalizes the distribution to one where groups of $\alpha$ galaxies appear with Poissonian noise and where the normalization coefficient $\mathcal{N}$ is needed to ensure that $\int P(N_T = N)dN = 1$. We find $\mathcal{N}$ to be very close to $\alpha^{-1}$ and identical to $\alpha^{-1}$ in the case where $\alpha$ is an integer value.

To account for the observed increase of super-Poissonianity with density, we allow $\alpha$ to depend on $\delta_{m,T}$,

$$\alpha(\delta_{m,T}) = \alpha_0 + \alpha_1 \times \delta_{m,T}. \quad (4.61)$$

This indeed leads to a $\delta_{m,T}$-dependence of the variance of galaxy counts that is close to the relation mentioned in Eq. (4.46). In our analysis we treat $\alpha_0$ and $\alpha_1$ as free parameters within the ranges

$$\begin{aligned} \alpha_0 &\in [0.1, 3.0] \\ \alpha_1 &\in [-1.0, 4.0]. \end{aligned} \quad (4.62)$$

In principle we could allow any value $\alpha_0 > 0$ but we choose the boundary $0.1 < \alpha_0$ because it is numerically difficult (and slow) to predict the CiC histogram for values close to $\alpha_0 = 0.0$. The other boundaries roughly enclose the $2-\sigma$ region of the posterior distribution of $\alpha_0$ and $\alpha_1$ we infer with DES Y1 like errors around the mean signal measured in Buzzard (after marginalizing over our other model parameters, cf. Appendix E). Also, the constraints on $\alpha_0$ and $\alpha_1$ we derive in [39] on DES data are well contained within our prior distributions.

Nevertheless, these priors must be considered mildly informative. We expect that even stronger priors can be motivated theoretically. [79] find that for their most massive halos shot-noise is reduced with respect to Poisson expectation by a factor of $\approx 2$, indicating that $\alpha_0 \gtrsim 0.5$, while for halo masses comparable to redmaGiC halo masses
(cf. Clampitt et al. [81]) they find shot-noise to be close to Poissonian. Also, there is evidence that the fraction of red galaxies that are satellites (resp. the fraction of satellite galaxies that are red) increases with environment density (see e.g. [82,83]). According to [79] this will cause an increase of galaxy stochasticity with environment density, corresponding to $\alpha_1 > 0.0$. We intend to investigate implications of models for halo occupation distributions (HOD) on our shot-noise parametrizations in future studies (see, e.g., the work by [84,85] on connecting HOD models and parametric models of galaxy bias and stochasticity).

To compare this parametric shot-noise model to our simulations, we are nevertheless interested in the particular value of $\alpha_0$ and $\alpha_1$ that describe these simulations. From the tuples of $(N_T, \delta_{m,T})_i$ measured in our simulated maps, we can constrain bias and the $\alpha_{0/1}$ parameters with a likelihood $L$ that is simply the product of the probabilities of the individual tuples from Eq. (4.60),

$$ \ln L = \sum_i \left[ \begin{aligned} & N_i / \alpha(\delta_{m,T}^i) \ln \left[ \begin{aligned} & (\bar{N}_i / \alpha(\delta_{m,T}^i))(1 + b \times \delta_{m,T}^i) \\ & - (\bar{N}_i / \alpha(\delta_{m,T}^i))(1 + b \times \delta_{m,T}^i) \\ & - \ln G[N_i / \alpha(\delta_{m,T}^i) + 1] \end{aligned} \right] \end{aligned} \right] - \ln \alpha(\delta_{m,T}^i). \quad (4.63) $$

Because the tuples have correlated counts and densities, this is not an exact expression for the likelihood of our measurements, but it should be sufficient to obtain reasonable best fit values for $b$, $\alpha_0$ and $\alpha_1$. We determine the uncertainties of these best-fit values by finding the maximum of Eq. (4.63) on jackknife resamplings of the simulations. The resulting parameter values are shown in Table III and displayed in Figs. 7 and 8. We find that our simulated REDMAGIC galaxies are indeed well described as linearly biased tracers of the density field with a small, but significant, scale and density dependent super-Poissonian shot-noise.

## D. Summary of fiducial model and approximations therein

For each ingredient (i) to (iii) of the framework described in Sec. II we have introduced at least two different modeling ansatzes. We want to once more describe our baseline model built from these ansatzes (cf. also Sec. II B).

<table><caption>TABLE III. Best-fit values galaxy bias and shot-noise parameters of our simulated tracer galaxies within the model presented in Sec. IV C 2. Error bars are again estimated from a jackknife approach.</caption><thead><tr><th>Smoothing scale [arcmin]</th><th>β̃</th><th>α<sub>0</sub></th><th>α<sub>1</sub></th></tr></thead><tbody><tr><td>10</td><td>1.54 ± 0.001</td><td>1.15 ± 0.001</td><td>0.22 ± 0.003</td></tr><tr><td>20</td><td>1.54 ± 0.002</td><td>1.26 ± 0.002</td><td>0.29 ± 0.010</td></tr><tr><td>30</td><td>1.54 ± 0.002</td><td>1.39 ± 0.003</td><td>0.45 ± 0.020</td></tr></tbody></table>

This is the model we consider in Sec. V and which we use in the data analysis presented in Gruen et al. [39].

(i) $p(\delta_{m,T})$: We find that the log-normal model (section IVA 1) and our model based in cylindrical collapse (Sec. IVA 2) describe the PDF of projected density contrast equally well. The computations based on the log-normal model are however significantly faster. Hence in our fiducial analysis we employ the log-normal model.

(ii) $\langle \kappa_{<\theta} | \delta_{m,T} \rangle$: We also introduced a log-normal model (Sec. IV B 1) and a model based on cylindrical collapse (Sec. IV B 2) for the convergence profile around lines-of-sight with fixed density contrast $\delta_{m,T}$. Both models lead to almost identical predictions for the density split lensing signal. Hence we again choose the log-normal model for our fiducial analysis, because of the shorter computation time.

(iii) $P(N_T | \delta_{m,T})$: We introduced two models for the distribution of tracer counts $N_T$ in lines-of-sight of matter density $\delta_{m,T}$. The first was based on linear galaxy bias $b$ and galaxy-matter-correlation coefficient $r$ (Sec. IV C 1). The second was based on an alternative definition of galaxy bias and on two parameters $\alpha_0$ and $\alpha_1$ describing density dependent deviations from Poissonian shot-noise (Sec. IV C 2). Both models describe the behavior of our simulated tracer galaxies in Buzzard-v1.1 similarly well. But anticipating that real galaxies might behave in a more complicated way, we will consider both ansatzes in our fiducial analysis.

In the following list, we are summarizing the approximations that went into the derivation of our baseline model.

(1) We assumed, that for fixed value of $\delta_{m,T}$ the convergence within angular radius $\theta$ is not dependent on $N_T$ [cf. Eq. (2.8)].

(2) All second order moments in our formalism are computed with a HALOFIT power spectrum [46] using an analytic approximation for the transfer function [47].

(3) Equations (4.5), (4.6), (4.22) and (4.26) employ a small angle and a Limber-like approximation [following [49]].

(4) We compute the cumulant generating function of density contrast in long cylinders by means of the cylindrical collapse approximation (cf. Sec. IVA 2).

(5) We assume that the tree-level result of the cumulant generating function can be corrected for the full nonlinear evolution of the density field by means of Eqs. (4.16) and (4.17).

(6) We approximate the PDF of $\delta_{m,T}$ resulting from a cylindrical collapse approximation by a log-normal PDF (cf. Sec. IVA 1).

(7) We employed approximations similar to (4), (5), and (6) for the joint distribution of $\delta_{m,T}$ and $\kappa_{<\theta}$ [cf. Sec. IV B 2, Eq. (4.39) and Sec. IV B 1].
(8) We assume that galaxies are linearly biased tracers of the density field. We consider two different models for shot-noise (resp. stochasticity), assuming that the full distribution $P(N_T|\delta_{m,T})$ is well described by either two parameters (b, r) or three parameters ($\tilde{b}$, $\alpha_0$, $\alpha_1$).

Despite this long list of approximations, this baseline model describes our measurements in the Buzzard simulations well within DES Y1 error bars (cf. Fig. 1). As shown in the next section, the model is also accurate enough to recover the true cosmology of our simulation within DES Y1 uncertainties in a simulated likelihood analysis. In Gruen et al. [39] (and using an extended set of simulations) we furthermore show that the values of $\chi^2$ found between our fiducial model and individual simulation measurements are consistent with the $\chi^2$-distribution expected from our number of data points, and that the coverage (i.e., the fraction of times the true simulation cosmology is within the confidence interval) matches expectations.

## V. RECOVERING COSMOLOGY IN N-BODY SIMULATIONS

In this section we want to test whether the modeling that was described in Secs. II and IV is sufficient to recover the cosmology underlying a density split data vector measured in N-body simulations. The simulations we use are described in Sec, III A. They are the same simulations against which we tested the ingredients of our model in the previous section. A likelihood analysis based on a density split data vector measured in these simulations is presented in Sec. VA. We only run a cosmological analysis on the mean data vector measured on 4 DES-Y1 realisations. The goal of this is to show that any possible systematic deviations between our modeling of density split statistics and the behaviour of our N-body simulations is smaller than the statistical uncertainties of DES-Y1. A more extensive validation of our likelihood pipeline is presented in Gruen et al. [39].

### A. Simulated likelihood analysis

We now measure the data vector that was described in Sec. II C in 4N-body realizations of DES-Y1. We always use the mean of these 4 data vectors. In order to further reduce the noise of this measurement, we turn off the shape noise in our simulated source catalogs, i.e. we measure our signal directly from the gravitational shear acting each galaxy. We then run Monte-Carlo Markov Chains of our model around this data vector. For this we assume a Gaussian likelihood function with the full covariance (i.e. including shape noise) that was estimated by Gruen et al. [39] for a DES-Y1 data set. The goal of this analysis is to test whether the fiducial cosmology of the Buzzard simulations is well contained within the $1\sigma$ contraints derived from this likelihood analysis. A more extensive

validation of our likelihood pipeline is presented in Gruen et al. [39].

In the top panel of Fig. 10 we show the $1\sigma$ and $2\sigma$ constraints obtained from our simulated likelihood analysis in the $\Omega_m-\sigma_8$ plane after marginalizing over different sets of parameters. First, we only vary $\Omega_m$ and $\sigma_8$, setting other cosmological parameters to the inputs of the Buzzard

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

FIG. 10. To test our model for possible systematic deviations from N-body simulations, we try to recover the Buzzard cosmology in a simulated likelihood analysis. Top panel: $1\sigma$ and $2\sigma$ contours in the $\Omega_m-\sigma_8$ plane from a likelihood computed around the mean of 4 shape noise free realisations of DES Y1 (but assuming the full covariance matrix for a single DES Y1). The green contours are marginalized over $\Omega_b, n_s, h_{100}$, redMaGC galaxy bias $b$, and galaxy-matter correlation coefficient $r$. For the parameters $\Omega_b, n_s, h_{100}$ we have assumed the same flat priors as used in the DES Y1 combined probes analysis presented in DES Collaboration [21]. The red contours are marginalized only over $b$ and $r$ and the blue contours only vary $\Omega_m$ and $\sigma_8$. Even when going to this small parameter space, our model agrees with Buzzard within $1\sigma$ errors of DES Y1. Bottom panel: Same contours but in the $\Sigma_8 - \Delta S_3/S_3$ plane and varying one additional parameter, $\Delta S_3/S_3$.
simulations and the parameters connecting galaxy count and matter density to the values we found from the Buzzard galaxy and density maps. Note that those values would be inaccessible in a real measurement. The corresponding constraints are very tight, but the fiducial values of our parameters are still well contained in the $1\sigma$ contour. Then, we also marginalize over the galaxy bias $b$ and the galaxy stochasticity $r$, demanding that $0 < r \le 1$. The contours now widen, and the fiducial values of $\Omega_m$ and $\sigma_8$ are still located well within the corresponding $1\sigma$ contour. Finally, we also marginalize over $\Omega_b$, $n_s$, and $h_{100}$, using the same informative priors that have been used in the DES-Y1 combined probes analysis [21]. The contours widen further, but our model and our simulations still agree well within $1\sigma$ uncertainty.

In the bottom panel of Fig. 10 we repeat this analysis, but now also vary the parameter $\Delta S_3/S_3$ that was introduced in Sec. II C. This parameter allows for deviations of the 3-point statistics of the density field from our fiducial model. Within our statistical uncertainties we find that the scaling between 3-point and 2-point statistics in our simulations is well described by our fiducial assumptions ($\Delta S_3/S_3 = 0$).

We repeat this analysis with our alternative shot-noise model in appendix E.

## VI. DISCUSSION AND CONCLUSIONS

In this work we introduced density split statistics, a technique to separately measure contributions to weak lensing and counts-in-cells from regions of different foreground galaxy density. Based on the pioneering work of Bernardeau [48], Bernardeau and Valageas [49] and Valageas [50] (see also references therein) on modeling the cosmic density PDF we were able to model the density split lensing signal as well as the counts-in-cells histogram from basic principles. With the help of this model, we then showed that density split statistics has two features that make it a potentially powerful cosmological probe:

(i) it is able to constrain the cosmological parameters $\Omega_m$ and $\sigma_8$ even if the relation of galaxy density and matter density is assumed to have 2 degrees of freedom: galaxy bias and galaxy-matter-correlation coefficient,

(ii) it is able to constrain the amplitude of 3-point statistics of the density field with almost no degeneracy to constraints on the amplitude of 2-point statistics.

In our fiducial model we predict 3-point statistics from cosmological perturbation theory. Deviations from that fiducial prediction may hint to nonstandard physics, that affect overdense and underdense parts of the matter field differently, or to any non-linear dynamics or small scale physics that break the scaling relations of $\Lambda$CDM perturbation theory. We showed that a DES-Y5 data set combined

with data form the cosmic microwave background can measure the amplitude of 3-points statistics to a $1\sigma$ accuracy of $\lesssim 5\%$. This is a conservative estimate since our projections neglect the fact that DES-Y5 will be a deeper data set than DES-Y1. Also, we so far neglected the possibility of a combined analysis including density split statistics and measurements of 2-point correlation functions.

Using measurements in high-resolution N-body simulations we showed that our model of the density split lensing signal and the counts-in-cells histogram is accurate to well within the statistical uncertainties of the DES-Y1 data set. Especially, in a mock likelihood analysis we were able to recover the input cosmology of our simulations to well within DES-Y1 parameter errors. Cosmological constraints from DES-Y1 data based on density split statistics are presented in Gruen *et al.* [39].

## ACKNOWLEDGMENTS

O. F. was supported by SFB-Transregio 33 ‘The Dark Universe’ by the Deutsche Forschungsgemeinschaft (DFG). Support for D. G. was provided by NASA through Einstein Postdoctoral Fellowship Grant No. PF5-160138 awarded by the Chandra X-ray Center, which is operated by the Smithsonian Astrophysical Observatory for NASA under contract No. NAS8-03060. O. F. and S. H. acknowledge support by the DFG cluster of excellence ‘Origin and Structure of the Universe’ (www.universe-cluster.de). Part of our computations have been carried out on the computing facilities of the Computational Center for Particle and Astrophysics (C2PAP). This paper has gone through internal review by the DES collaboration. We want to thank all the members of the DES WL, LSS and Theory working groups that have contributed with helpful comments and discussions. We also want to thank the anonymous journal referee for very helpful comments. Funding for the DES Projects has been provided by the U.S. Department of Energy, the U.S. National Science Foundation, the Ministry of Science and Education of Spain, the Science and Technology Facilities Council of the United Kingdom, the Higher Education Funding Council for England, the National Center for Supercomputing Applications at the University of Illinois at Urbana-Champaign, the Kavli Institute of Cosmological Physics at the University of Chicago, the Center for Cosmology and Astro-Particle Physics at the Ohio State University, the Mitchell Institute for Fundamental Physics and Astronomy at Texas A&M University, Financiadora de Estudos e Projetos, Fundação Carlos Chagas Filho de Amparo à Pesquisa do Estado do Rio de Janeiro, Conselho Nacional de Desenvolvimento Científico e Tecnológico and the Ministério da Ciência, Tecnologia e Inovação, the Deutsche Forschungsgemeinschaft and the Collaborating Institutions in the Dark Energy Survey. The DES data management system is supported by the National Science Foundation under Grant No. AST-1138766.
The Collaborating Institutions are Argonne National Laboratory, the University of California at Santa Cruz, the University of Cambridge, Centro de Investigaciones Enérgicas, Medioambientales y Tecnológicas-Madrid, the University of Chicago, University College London, the DES-Brazil Consortium, the University of Edinburgh, the Eidgenössische Technische Hochschule (ETH) Zürich, Fermi National Accelerator Laboratory, the University of Illinois at Urbana-Champaign, the Institut de Ciències de l'Espai (IEEC/CSIC), the Institut de Física d'Altes Energies, Lawrence Berkeley National Laboratory, the Ludwig-Maximilians Universität München and the associated Excellence Cluster Universe, the University of Michigan, the National Optical Astronomy Observatory, the University of Nottingham, The Ohio State University, the University of Pennsylvania, the University of Portsmouth, SLAC National Accelerator Laboratory, Stanford University, the University of Sussex, and Texas A&M University. The DES participants from Spanish institutions are partially supported by MINECO under Grants No. AYA2012-39559, No. ESP2013-48274, No. FPA2013-47986, and Centro de Excelencia Severo Ochoa SEV-2012-0234. Research leading to these results has received funding from the European Research Council under the European Union's Seventh Framework Programme (FP7/2007-2013) including ERC Grant Agreements No. 240672, No. 291329, and No. 306478. D. G. is an Einstein Fellow.

## APPENDIX A: FRIEDMANN EQUATIONS, LINEAR GROWTH, SPHERICAL COLLAPSE, AND CYLINDRICAL COLLAPSE

Throughout this section we set $G = 1 = c$ and we assume a flat $\Lambda$CDM universe. In proper comoving time $t$ the Friedmann equations take the form

$$H^2 = \frac{8\pi}{3}(\bar{\rho}_m + \bar{\rho}_{\Lambda}) \qquad (A1)$$

$$\frac{dH}{dt} + H^2 = -\frac{4\pi}{3}(\bar{\rho}_m - 2\bar{\rho}_{\Lambda}), \qquad (A2)$$

where $H = \frac{d}{dt} \ln a$. In conformal time, defined by $dt = a d\tau$, this changes to

$$\mathcal{H}^2 = \frac{8\pi}{3}a^2(\bar{\rho}_m + \bar{\rho}_{\Lambda}) \qquad (A3)$$

$$\frac{d\mathcal{H}}{d\tau} = -\frac{4\pi}{3}a^2(\bar{\rho}_m - 2\bar{\rho}_{\Lambda}), \qquad (A4)$$

where $\mathcal{H} = \frac{d}{d\tau} \ln a$. We will from now put $\frac{d}{d\tau} \equiv$.

In the Newtonian approximation, i.e., on scales much smaller than the curvature horizon of the universe, the evolution of a spherical, cylindrical or planar perturbation $\delta$ is given by the equation

$$\frac{d^2\delta}{dt^2} + 2H\frac{d\delta}{dt} - \frac{N+1}{N}\frac{1}{1+\delta}\left(\frac{d\delta}{dt}\right)^2 = 4\pi\bar{\rho}_m\delta(1+\delta), \qquad (A5)$$

where $N = 3$ for a spherical perturbation, $N = 2$ for a cylindrical perturbation and $N = 1$ for a planar perturbation (see Mukhanov [86] where this is demonstrated for $N = 3$ and $N = 1$).

In conformal time this equation reads

$$\ddot{\delta} + \mathcal{H}\dot{\delta} - \frac{N+1}{N}\frac{\dot{\delta}^2}{1+\delta} = 4\pi\bar{\rho}_m a^2 \delta(1+\delta). \qquad (A6)$$

To linear order in $\delta$ this becomes

$$\ddot{\delta} + \mathcal{H}\dot{\delta} = 4\pi\bar{\rho}_m a^2 \delta, \qquad (A7)$$

which is indeed independent of the particular shape of the perturbation.

## APPENDIX B: $\Lambda$CDM PERTURBATION THEORY

Consider the matter density contrast $\delta$ and the divergence of the velocity field $\theta = \nabla v$. In the Newtonian approximation the Fourier space equations of motion of $\delta$ and $\theta$ are [cf. [51]]

$$\begin{aligned} & \frac{\partial \tilde{\delta}(\mathbf{k}, \tau)}{\partial \tau} + \tilde{\theta}(\mathbf{k}, \tau) \\ & = -\int d^3 k_1 d^3 k_2 \delta_D(\mathbf{k} - \mathbf{k}_{12}) \alpha(\mathbf{k}_1, \mathbf{k}_2) \tilde{\delta}(\mathbf{k}_1, \tau) \tilde{\theta}(\mathbf{k}_2, \tau) \\ & \frac{\partial \tilde{\theta}(\mathbf{k}, \tau)}{\partial \tau} + \mathcal{H} \tilde{\theta}(\mathbf{k}, \tau) + \frac{3\Omega_m^0 H_0^2}{2a} \tilde{\delta}(\mathbf{k}, \tau) \\ & = -\int d^3 k_1 d^3 k_2 \delta_D(\mathbf{k} - \mathbf{k}_{12}) \beta(\mathbf{k}_1, \mathbf{k}_2) \tilde{\theta}(\mathbf{k}_1, \tau) \tilde{\theta}(\mathbf{k}_2, \tau), \end{aligned} \qquad (B1)$$

where $\mathbf{k}_{12} = \mathbf{k}_1 + \mathbf{k}_2$ and $\alpha$ and $\beta$ are given by

$$\begin{aligned} \alpha(\mathbf{k}_1, \mathbf{k}_2) &= 1 + \frac{1}{2} \frac{\mathbf{k}_1 \cdot \mathbf{k}_2}{k_1 k_2} \left( \frac{k_1}{k_2} + \frac{k_2}{k_1} \right) \\ \beta(\mathbf{k}_1, \mathbf{k}_2) &= \frac{1}{2} \frac{\mathbf{k}_1 \cdot \mathbf{k}_2}{k_1 k_2} \left( \frac{k_1}{k_2} + \frac{k_2}{k_1} \right) + \frac{(\mathbf{k}_1 \cdot \mathbf{k}_2)^2}{k_1^2 k_2^2}. \end{aligned} \qquad (B2)$$

In the following we will abbreviate the integrals involving $\alpha$ and $\beta$ as $\alpha[\tilde{\delta}, \tilde{\theta}, \mathbf{k}]$ and $\beta[\tilde{\theta}, \tilde{\theta}, \mathbf{k}]$. In perturbation theory we write $\tilde{\delta}$ and $\tilde{\theta}$ as

$$\begin{aligned} \tilde{\delta}(\mathbf{k}, \tau) &= \sum_{n=1}^{\infty} \delta_n(\mathbf{k}, \tau) \\ \text{and} \quad \tilde{\theta}(\mathbf{k}, \tau) &= -\frac{\partial \ln D_+(\tau)}{\partial \tau} \sum_{n=1}^{\infty} \theta_n(\mathbf{k}, \tau), \end{aligned} \qquad (B3)$$
where $\delta_n$ and $\theta_n$ are of order $n$ in the linearly approximated fields $\delta_1$ and $\theta_1$ and $D_+(\tau)$ is the linear growth factor. (We will ignore the decaying mode of linear growth here.) At linear order we have

$$ \delta_1(\mathbf{k}, \tau) = \theta_1(\mathbf{k}, \tau) = \frac{D_+(\tau)}{D_+(\tau_0)} \delta_1(\mathbf{k}, \tau_0) \equiv D_+(\tau) \delta_{1,1}(\mathbf{k}), \quad (B4) $$

where we have assumed that $D_+(\tau_0) = 1$ at present time $\tau_0$ and introduced the notation $\delta_{1,1}(\mathbf{k}) = \delta_1(\mathbf{k}, \tau_0)$ whose purpose will become clear at the end of this section. To get $\tilde{\delta}$ at second order we first note that

$$ \frac{\partial \tilde{\theta}(\mathbf{k}, \tau)}{\partial \tau} + \mathcal{H} \tilde{\theta}(\mathbf{k}, \tau) = \frac{1}{a} \frac{\partial}{\partial \tau} (a \tilde{\theta}(\mathbf{k}, \tau)). \quad (B5) $$

Hence, multiplying the first of Eq. (B1) with $a$ and differentiating with respect to $\tau$ and then multiplying with $1/a$ we get

$$ \begin{aligned} \frac{\partial^2 \delta_2(\mathbf{k}, \tau)}{\partial^2 \tau} + \mathcal{H} \frac{\partial \delta_2(\mathbf{k}, \tau)}{\partial \tau} - \frac{3\Omega_m^0 H_0^2}{2a} \delta_2(\mathbf{k}, \tau) &= \left(\frac{\partial D_+}{\partial \tau}\right)^2 \beta[\delta_{1,1}, \delta_{1,1}, \mathbf{k}] + \left(D \frac{\partial^2 D_+}{\partial \tau^2} + D \mathcal{H} \frac{\partial D_+}{\partial \tau} + \left(\frac{\partial D_+}{\partial \tau}\right)^2\right) \alpha[\delta_{1,1}, \delta_{1,1}, \mathbf{k}] \\ &= \left(\frac{\partial D_+}{\partial \tau}\right)^2 \beta[\delta_{1,1}, \delta_{1,1}, \mathbf{k}] + \left(\frac{3\Omega_m^0 H_0^2}{2a} D^2 + \left(\frac{\partial D_+}{\partial \tau}\right)^2\right) \alpha[\delta_{1,1}, \delta_{1,1}, \mathbf{k}] \\ &= \alpha[\delta_{1,1}, \delta_{1,1}, \mathbf{k}] \left(\frac{3\Omega_m^0 H_0^2}{2a} D^2 + 2 \left(\frac{\partial D_+}{\partial \tau}\right)^2\right) \\ &\quad + (\beta[\delta_{1,1}, \delta_{1,1}, \mathbf{k}] - \alpha[\delta_{1,1}, \delta_{1,1}, \mathbf{k}]) \left(\frac{\partial D_+}{\partial \tau}\right)^2. \end{aligned} \quad (B8) $$

This is solved by

$$ \delta_2(\mathbf{k}, \tau) = D_{2,1}(\tau) \delta_{2,1}(\mathbf{k}) + D_{2,2}(\tau) \delta_{2,2}(\mathbf{k}) \quad (B9) $$

where

$$ \begin{aligned} D_{2,1}(\tau) &\equiv D_+^2(\tau), & \delta_{2,1}(\mathbf{k}) &= \alpha[\delta_{1,1}, \delta_{1,1}, \mathbf{k}], \\ \delta_{2,2}(\mathbf{k}) &= \beta[\delta_{1,1}, \delta_{1,1}, \mathbf{k}] - \alpha[\delta_{1,1}, \delta_{1,1}, \mathbf{k}] \end{aligned} \quad (B10) $$

and $D_{2,2}$ is given by the differential equation

$$ \frac{\partial^2 D_{2,2}(\tau)}{\partial^2 \tau} + \mathcal{H} \frac{\partial D_{2,2}(\tau)}{\partial \tau} - \frac{3\Omega_m^0 H_0^2}{2a} D_{2,2}(\tau) = \left(\frac{\partial D_+}{\partial \tau}\right)^2. \quad (B11) $$

## 1. Second order of $\delta$ in Einstein-de Sitter universe

Let us define $1 - \mu \equiv D_{2,2}/D_+^2$. Then the general solution to (B8) is given by

$$ \begin{aligned} & \frac{\partial^2 \tilde{\delta}(\mathbf{k}, \tau)}{\partial^2 \tau} + \mathcal{H} \frac{\partial \tilde{\delta}(\mathbf{k}, \tau)}{\partial \tau} + \frac{\partial \tilde{\theta}(\mathbf{k}, \tau)}{\partial \tau} + \mathcal{H} \tilde{\theta}(\mathbf{k}, \tau) \\ &= -\alpha \left[ \frac{\partial \tilde{\delta}}{\partial \tau}, \tilde{\theta}, \mathbf{k} \right] - \alpha \left[ \tilde{\delta}, \frac{\partial \tilde{\theta}}{\partial \tau}, \mathbf{k} \right] - \mathcal{H} \alpha[\tilde{\delta}, \tilde{\theta}, \mathbf{k}]. \end{aligned} \quad (B6) $$

Now the second of Eq. (B1) can be used to eliminate $\tilde{\theta}$ from the right-hand side, giving

$$ \begin{aligned} & \frac{\partial^2 \tilde{\delta}(\mathbf{k}, \tau)}{\partial^2 \tau} + \mathcal{H} \frac{\partial \tilde{\delta}(\mathbf{k}, \tau)}{\partial \tau} - \frac{3\Omega_m^0 H_0^2}{2a} \tilde{\delta}(\mathbf{k}, \tau) \\ &= \beta[\tilde{\theta}, \tilde{\theta}, \mathbf{k}] - \alpha \left[ \frac{\partial \tilde{\delta}}{\partial \tau}, \tilde{\theta}, \mathbf{k} \right] - \alpha \left[ \tilde{\delta}, \frac{\partial \tilde{\theta}}{\partial \tau}, \mathbf{k} \right] - \mathcal{H} \alpha[\tilde{\delta}, \tilde{\theta}, \mathbf{k}]. \end{aligned} \quad (B7) $$

At second order in perturbation theory this equation becomes

$$ \delta_2(\mathbf{k}, \tau) = D_+^2([1 - \mu]\beta[\delta_{1,1}, \delta_{1,1}, \mathbf{k}] + \mu\alpha[\delta_{1,1}, \delta_{1,1}, \mathbf{k}]). \quad (B12) $$

In an Einstein-de Sitter universe where $\Omega_m^0 = 1$ and $D \equiv a$ we have

$$ D_{2,2} = \frac{2}{7} D_+^2, \quad \mu = \frac{5}{7} \quad (B13) $$

and $\delta_2$ is hence given by

$$ \begin{aligned} \delta_2(\mathbf{k}, \tau) &= D_+^2 \left( \frac{2}{7} \beta[\delta_{1,1}, \delta_{1,1}, \mathbf{k}] + \frac{5}{7} \alpha[\delta_{1,1}, \delta_{1,1}, \mathbf{k}] \right) \\ &= \int d^3 k_1 d^3 k_2 \delta_D(\mathbf{k} - \mathbf{k}_{12}) F_2(\mathbf{k}_1, \mathbf{k}_2) \\ &\quad \times \delta_{1,1}(\mathbf{k}_1) \delta_{1,1}(\mathbf{k}_2) \end{aligned} \quad (B14) $$

with
$$ F_2(\mathbf{k}_1, \mathbf{k}_2) = \frac{5}{7}\alpha(\mathbf{k}_1, \mathbf{k}_2) + \frac{2}{7}\beta(\mathbf{k}_1, \mathbf{k}_2) \\ = \frac{5}{7} + \frac{1}{2}\frac{\mathbf{k}_1 \cdot \mathbf{k}_2}{k_1 k_2} \left(\frac{k_1}{k_2} + \frac{k_2}{k_1}\right) + \frac{2}{7}\frac{(\mathbf{k}_1 \cdot \mathbf{k}_2)^2}{k_1^2 k_2^2}. \qquad (B15) $$

## 2. Second order of $\delta$ in $\Lambda$CDM universe

In a general $\Lambda$CDM universe the function $F_2$ becomes time dependent. It is given by

$$ F_2(\mathbf{k}_1, \mathbf{k}_2, \tau) = \mu(\tau) + \frac{1}{2} \frac{\mathbf{k}_1 \cdot \mathbf{k}_2}{k_1 k_2} \left( \frac{k_1}{k_2} + \frac{k_2}{k_1} \right) + [1 - \mu(\tau)] \frac{(\mathbf{k}_1 \cdot \mathbf{k}_2)^2}{k_1^2 k_2^2}. \qquad (B16) $$

A useful property or this kernel is that

$$ F_2(\mathbf{k}, -\mathbf{k}, \tau) = \mu(\tau) + \frac{1-k^2}{2k^2}(1+1) + [1-\mu(\tau)] \frac{k^4}{k^4} \\ = \mu(\tau) - 1 + 1 - \mu(\tau) = 0. \qquad (B17) $$

Denoting the angle between $\mathbf{k}_1$ and $\mathbf{k}_2$ with $\phi$ one can also arrive at the following form of $F_2(\mathbf{k}_1, \mathbf{k}_2, \tau)$ which will be useful when computing the skewness of matter inside a long cylinder:

$$ \begin{aligned} \langle \delta_{1,1}(\mathbf{k}_1) \delta_{1,1}(\mathbf{k}_2) \delta_{1,1}(\mathbf{q}_1) \delta_{1,1}(\mathbf{q}_2) \rangle &= \langle \delta_{1,1}(\mathbf{k}_1) \delta_{1,1}(\mathbf{k}_2) \rangle \langle \delta_{1,1}(\mathbf{q}_1) \delta_{1,1}(\mathbf{q}_2) \rangle + \langle \delta_{1,1}(\mathbf{k}_1) \delta_{1,1}(\mathbf{q}_1) \rangle \langle \delta_{1,1}(\mathbf{k}_2) \delta_{1,1}(\mathbf{q}_2) \rangle \\ &\quad + \langle \delta_{1,1}(\mathbf{k}_1) \delta_{1,1}(\mathbf{q}_2) \rangle \langle \delta_{1,1}(\mathbf{k}_2) \delta_{1,1}(\mathbf{q}_1) \rangle \\ &= \delta_D(\mathbf{k}_1 + \mathbf{k}_2) \delta_D(\mathbf{q}_1 + \mathbf{q}_2) P_{\text{lin},0}(k_1) P_{\text{lin},0}(q_1) \\ &\quad + \delta_D(\mathbf{k}_1 + \mathbf{q}_1) \delta_D(\mathbf{k}_2 + \mathbf{q}_2) P_{\text{lin},0}(k_1) P_{\text{lin},0}(q_2) \\ &\quad + \delta_D(\mathbf{k}_1 + \mathbf{q}_2) \delta_D(\mathbf{k}_2 + \mathbf{q}_1) P_{\text{lin},0}(k_1) P_{\text{lin},0}(q_1) \end{aligned} \qquad (B21) $$

Because of Eq. (B17) the contribution of the first term to the bispectrum is zero. Using the symmetry $1 \leftrightarrow 2$ between the second and third term we hence get

$$ \begin{aligned} \langle \tilde{\delta}(\mathbf{k}_1, \tau) \tilde{\delta}(\mathbf{k}_2, \tau) \tilde{\delta}(\mathbf{k}_3, \tau) \rangle_{\text{2nd}} &= 2D_{\perp}^2 \delta_D(\mathbf{k}_1 + \mathbf{k}_2 + \mathbf{k}_3) F_2(\mathbf{k}_1, \mathbf{k}_2, \tau) P_{\text{lin},0}(k_1) P_{\text{lin},0}(k_2) \\ &\quad + \text{cycl.} \end{aligned} \qquad (B22) $$

## 4. Variance and skewness of long cylinder at leading order in perturbation theory

Consider a cylinder with radius $R$ and length $L$. In Fourier space the tophat filter for this cylinder is given by

$$ F_2(\mathbf{k}_1, \mathbf{k}_2, \tau) = \frac{1}{2} \left\{ \left( 1 + \frac{k_1}{k_2} \cos \phi \right) + \left( 1 + \frac{k_2}{k_1} \cos \phi \right) \right\} \\ + [1 - \mu(\tau)](\cos^2 \phi - 1). \qquad (B18) $$

## 3. Bispectrum and 3-point function at leading order

The bispectrum $B(k_1, k_2, k_3, \tau)$ of $\delta$ is defined by

$$ \begin{aligned} \langle \tilde{\delta}(\mathbf{k}_1, \tau) \tilde{\delta}(\mathbf{k}_2, \tau) \tilde{\delta}(\mathbf{k}_3, \tau) \rangle &= \delta_D(\mathbf{k}_1 + \mathbf{k}_2 + \mathbf{k}_3) B(k_1, k_2, k_3, \tau). \end{aligned} \qquad (B19) $$

At leading order in perturbation theory this can be calculated as

$$ \begin{aligned} \langle \tilde{\delta}(\mathbf{k}_1, \tau) \tilde{\delta}(\mathbf{k}_2, \tau) \tilde{\delta}(\mathbf{k}_3, \tau) \rangle_{\text{2nd}} &= D_{\perp}^2 \int d^3 q_1 d^3 q_2 \delta_D(\mathbf{k}_3 - \mathbf{q}_{12}) \\ &\quad \times F_2(\mathbf{q}_1, \mathbf{q}_2, \tau) \langle \delta_{1,1}(\mathbf{k}_1) \delta_{1,1}(\mathbf{k}_2) \delta_{1,1}(\mathbf{q}_1) \delta_{1,1}(\mathbf{q}_2) \rangle \\ &\quad + \text{cycl}, \end{aligned} \qquad (B20) $$

where cycl indicates that the integral on the right-hand side should appear for all possible permutations of $\mathbf{k}_1, \mathbf{k}_2$, and $\mathbf{k}_3$. Since we assume the linear density field to be a Gaussian random field, the expectation value on the left-hand side factorizes as

$$ W_{R,L}(\mathbf{k}) = \frac{1}{(2\pi)^3} W_L(k_{\parallel}) W_R(\mathbf{k}_{\perp}) \qquad (B23) $$

where we denote the component of $\mathbf{k}$ parallel to the cylinder axis with $k_{\parallel}$ and the components orthogonal to it are represented by the two-dimensional vector $\mathbf{k}_{\perp}$ and $W_L$ and $W_R$ given by

$$ W_L(k_{\parallel}) = \frac{\sin(Lk_{\parallel}/2)}{Lk_{\parallel}/2}, \qquad W_R(\mathbf{k}_{\perp}) = \frac{2J_1(Rk_{\perp})}{Rk_{\perp}}. \qquad (B24) $$

Here $k_{\perp} = |\mathbf{k}_{\perp}|$ and $J_{\nu}$ are the cylindrical bessel functions. At leading order or tree level in perturbation theory the variance of matter contrast within the cylinder is then given by
$$
\begin{align}
\langle \delta_{R,L}^2 \rangle_{\text{tree}}(\tau) &= D_+^2(\tau) \int dk_{\parallel,1} dk_{\parallel,2} d^2k_{\perp,1} d^2k_{\perp,2} W_L(k_{\parallel,1}) \nonumber \\
&\quad \times W_L(k_{\parallel,2}) W_R(\mathbf{k}_{\perp,1}) W_R(\mathbf{k}_{\perp,2}) \nonumber \\
&\quad \times \langle \delta_{1,1}(\mathbf{k}_1) \delta_{1,1}(\mathbf{k}_2) \rangle \nonumber \\
&= D_+^2(\tau) \int dk_{\parallel} d^2k_{\perp} W_L(k_{\parallel})^2 W_R(\mathbf{k}_{\perp})^2 \nonumber \\
&\quad \times P_{\text{lin},0}(\mathbf{k}). \tag{B25}
\end{align}
$$
$$
W_L(k_{\parallel})^2 \approx \frac{2\pi}{L} \delta_D(k_{\parallel}) \tag{B26}
$$
$$
\langle \delta_{R,L}^2 \rangle_{\text{tree}}(\tau) \approx \frac{2\pi D_+^2}{L} \int dk k W_R(k)^2 P_{\text{lin},0}(k). \tag{B27}
$$

Here $P_{\text{lin},0}(\mathbf{k})$ is today's linear power spectrum. For $L \gg R$
we can actually approximate $W_L$ by

The third moment at tree level is given by

$$
\begin{align}
\langle \delta_{R,L}^3 \rangle_{\text{tree}}(\tau) &= 3D_+^2 \int dk_{\parallel,1} dk_{\parallel,2} dq_{\parallel} d^2k_{\perp,1} d^2k_{\perp,2} d^2q_{\perp} W_L(k_{\parallel,1}) W_L(k_{\parallel,2}) W_L(q_{\parallel}) W_R(\mathbf{k}_{\perp,1}) W_R(\mathbf{k}_{\perp,2}) W_R(\mathbf{q}_{\perp}) \nonumber \\
&\quad \times \langle \delta_{1,1}(\mathbf{k}_1) \delta_{1,1}(\mathbf{k}_2) \delta_{2}(\mathbf{q}, \tau) \rangle \nonumber \\
&= 3D_+^4 \int dk_{\parallel,1} dk_{\parallel,2} dq_{\parallel} d^2k_{\perp,1} d^2k_{\perp,2} d^2q_{\perp} d^3q_1 d^3q_2 W_L(k_{\parallel,1}) W_L(k_{\parallel,2}) W_L(q_{\parallel}) W_R(\mathbf{k}_{\perp,1}) W_R(\mathbf{k}_{\perp,2}) W_R(\mathbf{q}_{\perp}) \nonumber \\
&\quad \times \delta_D(\mathbf{q} - \mathbf{q}_1 - \mathbf{q}_2) F_2(\mathbf{q}_1, \mathbf{q}_2, \tau) \langle \delta_{1,1}(\mathbf{k}_1) \delta_{1,1}(\mathbf{k}_2) \delta_{1,1}(\mathbf{q}_1) \delta_{1,1}(\mathbf{q}_2) \rangle \nonumber \\
&= 3D_+^4 \int dk_{\parallel,1} dk_{\parallel,2} d^2k_{\perp,1} d^2k_{\perp,2} d^3q_1 d^3q_2 W_L(k_{\parallel,1}) W_L(k_{\parallel,2}) W_R(\mathbf{k}_{\perp,1}) W_R(\mathbf{k}_{\perp,2}) \nonumber \\
&\quad \times W_R(\mathbf{q}_{\perp,1} + \mathbf{q}_{\perp,2}) W_L(q_{\parallel,1} + \mathbf{q}_{\parallel,2}) F_2(\mathbf{q}_1, \mathbf{q}_2, \tau) \langle \delta_{1,1}(\mathbf{k}_1) \delta_{1,1}(\mathbf{k}_2) \delta_{1,1}(\mathbf{q}_1) \delta_{1,1}(\mathbf{q}_2) \rangle. \tag{B28}
\end{align}
$$

Since we assume the linear density field to be a Gaussian random field, the expectation value on the left-hand-side factorizes as

$$
\begin{align}
\langle \delta_{1,1}(\mathbf{k}_1) \delta_{1,1}(\mathbf{k}_2) \delta_{1,1}(\mathbf{q}_1) \delta_{1,1}(\mathbf{q}_2) \rangle &= \langle \delta_{1,1}(\mathbf{k}_1) \delta_{1,1}(\mathbf{k}_2) \rangle \langle \delta_{1,1}(\mathbf{q}_1) \delta_{1,1}(\mathbf{q}_2) \rangle + \langle \delta_{1,1}(\mathbf{k}_1) \delta_{1,1}(\mathbf{q}_1) \rangle \langle \delta_{1,1}(\mathbf{k}_2) \delta_{1,1}(\mathbf{q}_2) \rangle \nonumber \\
&\quad + \langle \delta_{1,1}(\mathbf{k}_1) \delta_{1,1}(\mathbf{q}_2) \rangle \langle \delta_{1,1}(\mathbf{k}_2) \delta_{1,1}(\mathbf{q}_1) \rangle \nonumber \\
&= \delta_D(\mathbf{k}_1 + \mathbf{k}_2) \delta_D(\mathbf{q}_1 + \mathbf{q}_2) P_{\text{lin},0}(\mathbf{k}_1) P_{\text{lin},0}(\mathbf{q}_1) \nonumber \\
&\quad + \delta_D(\mathbf{k}_1 + \mathbf{q}_1) \delta_D(\mathbf{k}_2 + \mathbf{q}_2) P_{\text{lin},0}(\mathbf{k}_1) P_{\text{lin},0}(\mathbf{q}_2) \nonumber \\
&\quad + \delta_D(\mathbf{k}_1 + \mathbf{q}_2) \delta_D(\mathbf{k}_2 + \mathbf{q}_1) P_{\text{lin},0}(\mathbf{k}_1) P_{\text{lin},0}(\mathbf{q}_1) \tag{B29}
\end{align}
$$

Because of Eq. (B17) the contribution of the first term to the skewness is zero. Using the symmetry 1 ↔ 2 between the second and third term we hence get

$$
\begin{equation}
\begin{split}
\langle \delta_{R,L}^3 \rangle_{\text{tree}}(\tau) = {}& 6D_+^4 \int dq_{\parallel,1} dq_{\parallel,2} d^2q_{\perp,1} d^2q_{\perp,2} W_L(q_{\parallel,1}) W_L(q_{\parallel,2}) W_L(q_{\parallel,1} + q_{\parallel,2}) W_R(\mathbf{q}_1) W_R(\mathbf{q}_2) W_R(\mathbf{q}_1 + \mathbf{q}_2) \\
& \times P_{\text{lin},0}(q_1) P_{\text{lin},0}(q_2) F_2(\mathbf{q}_1, \mathbf{q}_2, \tau).
\end{split}
\tag{B30}
\end{equation}
$$

For $L \gg R$ we can use the approximation

$$
W_L(q_{\parallel,1}) W_L(q_{\parallel,2}) W_L(q_{\parallel,1} + q_{\parallel,2}) \approx \frac{1}{L^2} \delta_D^2(q_{\parallel,1}, q_{\parallel,2}). \tag{B31}
$$

This gives

$$
\langle \delta_{R,L}^3 \rangle_{\text{tree}}(\tau) = \frac{6D_+^4}{L^2} \int d^2 q_1 d^2 q_2 W_R(\mathbf{q}_1) W_R(\mathbf{q}_2) W_R(\mathbf{q}_1 + \mathbf{q}_2) P_{\text{lin},0}(q_1) P_{\text{lin},0}(q_2) F_2(\mathbf{q}_1, \mathbf{q}_2, \tau), \quad (B32)
$$

where we will consider all vectors to be 2-dimensional from now on. Using Eq. (B18) to express $F_2$ in terms of $q_1$, $q_2$, and $\phi$
we can simplify this to
$$ \langle \delta_{R,L}^3 \rangle_{\text{tree}}(\tau) = \frac{12\pi D_+^4}{L^2} \int dq_1 dq_2 q_1 W_R(q_1) q_2 W_R(q_2) P_{\text{lin},0}(q_1) P_{\text{lin},0}(q_2) \\ \times \int d\phi W_R \left[ \sqrt{q_1^2 + q_2^2 + 2q_1 q_2 \cos \phi} \right] F_2(q_1, q_2, \phi, \tau). \qquad (B33) $$

Using relations given in Bernardeau (1995) or Buchalter et al. (2000) one can simplify the integral over $\phi$ as

$$ \begin{aligned} & \int d\phi W_R \left[ \sqrt{q_1^2 + q_2^2 + 2q_1 q_2 \cos \phi} \right] F_2(q_1, q_2, \phi, \tau) \\ &= \frac{1}{2} \int d\phi W_R \left[ \sqrt{q_1^2 + q_2^2 + 2q_1 q_2 \cos \phi} \right] \left\{ \left( 1 + \frac{k_1}{k_2} \cos \phi \right) + \left( 1 + \frac{k_2}{k_1} \cos \phi \right) \right\} \\ & \quad + [1 - \mu(\tau)] \int d\phi W_R \left[ \sqrt{q_1^2 + q_2^2 + 2q_1 q_2 \cos \phi} \right] (\cos^2 \phi - 1) \\ &= 2\pi W_R(q_1) W_R(q_2) + \frac{\pi}{2} \left\{ W_R(q_1) R q_2 \left. \frac{\partial W_R(x)}{\partial x} \right|_{x=Rq_2} + W_R(q_2) R q_1 \left. \frac{\partial W_R(x)}{\partial x} \right|_{x=Rq_1} \right\} - \pi [1 - \mu(\tau)] W_R(q_1) W_R(q_2) \\ &= \pi [1 + \mu(\tau)] W_R(q_1) W_R(q_2) + \frac{\pi}{2} \frac{\partial}{\partial \ln R} \{ W_R(q_1) W_R(q_2) \}. \end{aligned} \qquad (B34) $$

For the third moment of $\delta_{R,L}$ this gives

$$ \begin{aligned} \langle \delta_{R,L}^3 \rangle_{\text{tree}}(\tau) &= [1 + \mu(\tau)] \frac{12\pi^2 D_+^4(\tau)}{L^2} \int dq_1 dq_2 q_1 q_2 W_R(q_1)^2 W_R(q_2)^2 P_{\text{lin},0}(q_1) P_{\text{lin},0}(q_2) \\ & \quad + \frac{6\pi^2 D_+^4(\tau)}{L^2} \int dq_1 dq_2 q_1 q_2 W_R(q_1) W_R(q_2) \frac{\partial}{\partial \ln R} \{ W_R(q_1) W_R(q_2) \} P_{\text{lin},0}(q_1) P_{\text{lin},0}(q_2) \\ &= 3[1 + \mu(\tau)] \left( \frac{2\pi D_+^2(\tau)}{L} \int dq_1 q W_R(q)^2 P_{\text{lin},0}(q) \right)^2 \\ & \quad + \frac{3\pi^2 D_+^4(\tau)}{L^2} \frac{\partial}{\partial \ln R} \int dq_1 dq_2 q_1 q_2 W_R(q_1)^2 W_R(q_2)^2 P_{\text{lin},0}(q_1) P_{\text{lin},0}(q_2) \\ &= 3[1 + \mu(\tau)] \left( \frac{2\pi D_+^2(\tau)}{L} \int dq_1 q W_R(q)^2 P_{\text{lin},0}(q) \right)^2 + \frac{3}{4} \frac{\partial}{\partial \ln R} \left( \frac{2\pi D_+^2(\tau)}{L^2} \int dq_1 q W_R(q)^2 P_{\text{lin},0}(q) \right)^2 \\ &= 3[1 + \mu(\tau)] \left( \langle \delta_{R,L}^2 \rangle_{\text{tree}}(\tau) \right)^2 + \frac{3}{4} \frac{\partial}{\partial \ln R} \left( \langle \delta_{R,L}^2 \rangle_{\text{tree}}(\tau) \right)^2. \end{aligned} \qquad (B35) $$

Especially we have

$$ S_3 = \frac{\langle \delta_{R,L}^3 \rangle_{\text{tree}}(\tau)}{\langle \delta_{R,L}^2 \rangle_{\text{tree}}(\tau)^2} = 3[1 + \mu(\tau)] + \frac{3}{2} \frac{\partial \ln \langle \delta_{R,L}^2 \rangle_{\text{tree}}(\tau)}{\partial \ln R}. \qquad (B36) $$

For an Einstein-de Sitter universe and a power law power spectrum $P(k) \sim k^n$ this gives $S_3 = 36/7 - 3/2(n + 2)$. The leading order prediction for $S_3$ is surprisingly good, even in the mildly nonlinear regime [see, [51] and references therein]. Hence in order to predict the nonlinear skewness, we simply employ the approximation

$$ \langle \delta_{R,L}^3 \rangle_{\text{non-lin}}(\tau) \approx S_3 \langle \delta_{R,L}^2 \rangle_{\text{non-lin}}(\tau)^2, \qquad (B37) $$

where we compute the nonlinear variance with the use of HALOFIT as detailed in Takahashi et al. [87] which is a revised version of Smith et al. [45].

## 5. The moment $\langle \delta_{R,A,L}^2 \delta_{R,B,L} \rangle_{\text{tree}}$

For predicting the density split lensing signal we are also interested in the moment $\langle \delta_{R_A,L}^2 \delta_{R_B,L} \rangle_{\text{tree}}$, where $R_A$ and $R_B$ are two different Radii. The above derivations can be generalized to give
$$ \langle \delta_{R_A,L}^2 \delta_{R_B,L} \rangle_{\text{tree}}(\tau) = \text{Var}(R_A) \text{Cov}(R_A, R_B) \left\{ \begin{aligned} & 2[1 + \mu(\tau)] + \frac{1}{2} \frac{\partial \ln \text{Var}(R_A)}{\partial \ln R_A} \\ & + \frac{\partial \ln \text{Cov}(R_A, R_B)}{\partial \ln R_A} \\ & + \text{Cov}(R_A, R_B)^2 \left\{ [1 + \mu(\tau)] + \frac{\partial \ln \text{Cov}(R_A, R_B)}{\partial \ln R_B} \right\} \end{aligned} \right\} \quad (B38) $$

where we defined

$$ \begin{aligned} \text{Var}(R_A) &= \langle \delta_{R_A,L}^2 \rangle_{\text{tree}}, \\ \text{Cov}(R_A, R_B) &= \langle \delta_{R_A,L} \delta_{R_B,L} \rangle_{\text{tree}}. \end{aligned} \quad (B39) $$

To correct this expression for the nonlinear evolution of the power spectrum, we compute Var($R_A$) and Cov($R_A, R_B$) with our HALOFIT power spectrum whenever they appear outside of the logarithmic derivatives. This is a generalization of the rescaling of Var($R_A$) by means of $S_3$.

For $R_B \gg R_A$ this rescaling is dominated by the first term on the right-hand side of Eq. (B38). For $R_B \approx R_A$ it reduces to Eq. (B37). As a consequence, using the procedure described around (4.39) to correct for the nonlinear power spectrum evolution yields a prediction for the density split lensing signal that is almost identical to the procedure described here. Also, it can be considered accurate to the extend that Eq. (B37) is accurate. We nevertheless rescale the 3rd order moments in the more elaborate way described here.

## 6. The moment $\langle \delta_{R_A,L}^n \delta_{R_B,L} \rangle_{\text{tree}}$

Using a diagrammatic representation of perturbation theory (see, e.g., [51]) one can see that the tree-level result for the moment $\langle \delta_{R_A,L}^n \delta_{R_B,L} \rangle_c$ will consist of terms that scale as

$$ \sim \text{Cov}(R_A, R_B)^k \text{Var}(R_A)^{n-k}, \quad 1 \le k \le n. \quad (B40) $$

For $R_B \approx R_A$ each of these scalings reduces to $\sim \text{Var}(R_A)^n$ [cf. (B37) and the definition of $S_{n+1}$ in (4.16)]. On the other hand, for $R_B \gg R_A$ the terms scaling as $\sim \text{Cov}(R_A, R_B) \text{Var}(R_A)^{n-1}$ are the dominant contributions (cf. the last section for the case $n=2$). This is why we use (4.39) when rescaling moments with $n > 2$ in $G_{\text{cyl},\theta_T,\theta_W,\theta_L}(q_l(w)Ly, w)$ [see also (4.26)].

# APPENDIX C: COMPARISON WITH MILLENIUM SIMULATION

In Fig. 11, we compare our model for the PDF of projected density contrast to another set of N-body simulations, the Millennium Run [MR [88]]. The MR has a smaller simulation volume of only $(500 \text{ h}^{-1} \text{ Mpc})^3$ co-moving, but a force resolution of $5 \text{ h}^{-1} \text{ kpc}$ that is 4–10

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

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

<figure><img src="image_15.png" /><figcaption>FIG. 11. The PDF of projected density contrast $\delta_{m,T}$ in the Millennium Run (MR) compared to our model. In each plot, the black line shows a histogram of $\delta_{m,T}$ measured from 64 patches of $4 \times 4 \text{ deg}^2$ made from the MR by projecting the 3D density contrast with a constant selection function $q_l$ between 0.19 $\lesssim$ $z \lesssim$ 0.43, i.e. with a constant comoving density between those redshifts. The blue lines display the PDF predicted by our PT-motivated log-normal model, and the red lines show a Gaussian PDF with the same variance. The grey band is using the subsample covariance to estimate the error on the mean of all patches [89].</figcaption></figure>
times higher than that of the Buzzard simulations. The fiducial model and the log-normal model describe the distribution of $\delta_{m,T}$ measured from the MR well considering the large statistical uncertainty on $p(\delta_{m,T})$ due to the limited simulated sky area.

## APPENDIX D: GALAXY STOCHASTICITY

Consider the field of galaxy density contrast $\delta_{g,T}$ and the field of matter density contrast $\delta_{m,T}$, where both fields are assumed to be smoothed over a fix circular aperture. The number of galaxies found inside such an aperture is assumed to be a Poissonian random variable with first and second moments for a given value of $\delta_{g,T}$ are given by

$$ \langle \hat{N} | \delta_{g,T} \rangle = \bar{N}(1 + \delta_{g,T}) \quad (D1) $$

and

$$ \langle \hat{N}^2 | \delta_{g,T} \rangle = \bar{N}(1 + \delta_{g,T}) + \bar{N}^2(1 + \delta_{g,T})^2. \quad (D2) $$

Let $\text{Var}_m$ be the variance of $\delta_{m,T}$ and $\text{Var}_g = b^2 \text{Var}_m$ the variance of $\delta_{g,T}$, where b is the galaxy bias. Then the galaxy stochasticity r is defined by $\text{Cov}_{mg} = rb\text{Var}_m$, i.e. it is the correlation coefficient of $\delta_{g,T}$ and $\delta_{m,T}$.

We will now assume both $\delta_{g,T}$ and $\delta_{m,T}$ to be joint log-normal random variables, i.e.,

$$ \begin{aligned} \delta_{m,T} &= [e^{n_m} - 1]\delta_{m,0} \\ \delta_{g,T} &= [e^{n_g} - 1]\delta_{g,0}, \end{aligned} \quad (D3) $$

where $n_m$ and $n_g$ have a joint Gaussian distribution and $\delta_{g,0} = b\delta_{m,0}$. The variances of $n_m$ and $n_g$ are given by

$$ \begin{aligned} \sigma_m^2 &= \ln \left\{ 1 + \frac{\text{Var}_m}{\delta_{m,0}^2} \right\} \\ \sigma_g^2 &= \ln \left\{ 1 + \frac{\text{Var}_g}{\delta_{g,0}^2} \right\} = \sigma_m^2 \end{aligned} \quad (D4) $$

and their covariance is given by

$$ \xi_{mg} = \ln \left\{ 1 + \frac{\text{Cov}_{mg}}{\delta_{m,0}\delta_{g,0}} \right\} = \ln \left\{ 1 + r \frac{\text{Var}_m}{\delta_{m,0}^2} \right\}. \quad (D5) $$

Let us denote the correlation coefficient of the Gaussian field by

$$ \rho = \frac{\xi_{mg}}{\sigma_m^2} = \frac{\ln \left\{ 1 + r \frac{\text{Var}_m}{\delta_{m,0}^2} \right\}}{\ln \left\{ 1 + \frac{\text{Var}_m}{\delta_{m,0}^2} \right\}}. \quad (D6) $$

Note that $\rho$ will depend on scale even of b and r do not.

Now we want to compute the conditional moments $\langle \delta_{g,T} | \delta_{m,T} \rangle$ and $\langle \delta_{g,T}^2 | \delta_{m,T} \rangle$. First,

$$ \begin{aligned} \langle e^{n_g} | n_m \rangle &= e^{\langle n_g | n_m \rangle + \sigma_g^2(1-\rho^2)/2} \\ &= e^{\rho(n_m + \sigma_m^2/2) - \sigma_g^2\rho^2/2} \\ &= e^{\sigma_g^2(\rho-\rho^2)/2} e^{\rho n_m}. \end{aligned} \quad (D7) $$

Second,

$$ \begin{aligned} \text{Var}(e^{n_g} | n_m) &= (e^{\sigma_g^2(1-\rho^2)} - 1)e^{2\langle n_g | n_m \rangle + \sigma_g^2(1-\rho^2)} \\ &= (e^{\sigma_g^2(1-\rho^2)} - 1)e^{\sigma_g^2(\rho-\rho^2)} e^{2\rho n_m}. \end{aligned} \quad (D8) $$

Now what is $\text{Var}(\hat{N}|\delta_{m,T})$?

$$ \begin{aligned} \langle \hat{N}^2 | \delta_{g,T} \rangle &= \langle \hat{N}^2 | n_m \rangle \\ &= \int d\delta_{g,T} p(\delta_{g,T} | n_m) (\bar{N}[1 + \delta_{g,T}] + \bar{N}^2[1 + \delta_{g,T}]^2) \\ &= \bar{N} + \bar{N}^2 + \int d\delta_{g,T} p(\delta_{g,T} | n_m) (\delta_{g,T}[\bar{N} + 2\bar{N}^2] + \bar{N}^2\delta_{g,T}^2) \\ &= \bar{N} + \bar{N}^2 + [\bar{N} + 2\bar{N}^2]\langle \delta_{g,T} | n_m \rangle + \bar{N}^2(\text{Var}(\delta_{g,T} | n_m) + \langle \delta_{g,T} | n_m \rangle^2) \\ \Rightarrow \text{Var}(\hat{N} | \delta_{m,T}) &= \bar{N}(1 + \langle \delta_{g,T} | n_m \rangle) + \bar{N}^2\text{Var}(\delta_{g,T} | n_m) \\ &= \bar{N}(1 + \langle \delta_{g,T} | n_m \rangle) + \bar{N}^2\delta_{g,0}^2\text{Var}(e^{n_g} | n_m). \end{aligned} \quad (D9) $$

The probability $P(N_A|\delta_{m,T})$ can be computed in a similar way, by numerically evaluating

$$ P(N_A|\delta_{m,T}) = \int d\delta_{g,T} p(\delta_{g,T}|\delta_{m,T}) P(N_A|\delta_{g,T}), \quad (D10) $$

where $p(\delta_{g,T}|\delta_{m,T})$ can be computed from basic relations for joint log-normal random variables.
# APPENDIX E: VALIDATION OF ALTERNATIVE SHOT-NOISE MODEL

In our data analysis [39] we investigate both shot-noise parametrizations introduced in Sec. IV C. We hence check whether our alternative shot-noise parametrization, i.e., the one that uses three parameters to describe the relation between matter and galaxies (b, $\alpha_0$ and $\alpha_1$, cf. Sec. IV C 2), recovers the true cosmology and shot-noise parameters of our mock data.

In Fig. 12 we show the posterior constraints derived for the two shot-noise parameters $\alpha_0$ and $\alpha_1$, when marginalizing over different sets of model parameters. Our priors $0.1 < \alpha_0 < 3.0$ and $-1.0 < \alpha_1 < 4.0$ are mildly informative. We however expect that even stronger priors can be motivated (cf. our discussion in Sec. IV C 2) and will investigate this in future work. Figure 13 shows that our alternative shot-noise parametrization also recovers the correct Buzzard cosmology (cf. Fig. 10, which presents the same test for our baseline model).

<figure><img src="image_16.png" /><figcaption>FIG. 12. 1$\sigma$ and 2$\sigma$ contours in the $\alpha_0$-$\alpha_1$ plane from a likelihood computed around the mean of 4 shape noise free realisations of DES Y1 (but assuming the full covariance matrix for a single DES Y1). The blue contour only varies $\alpha_0$ and $\alpha_1$. The red contour marginalizes over $\Omega_m$, $\sigma_8$ and galaxy bias b. The green contour additionally marginalizes over $\Omega_b$, $n_s$, $h_{100}$, assuming the priors used by [21]. And the black contour also allows variation of the parameter $\Delta S_3/S_3$. Dotted lines show the values of $\alpha_0$ and $\alpha_1$ that were found to describe our mock data best in Sec. IV C 2.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>FIG. 13. In analogy to figure 10, we test whether our alternative shot-noise parametrization can recover the Buzzard cosmology in a simulated likelihood analysis. Top panel: 1$\sigma$ and 2$\sigma$ contours in the $\Omega_m-\sigma_8$ plane from a likelihood computed around the mean of 4 shape noise free realisations of DES Y1 (but assuming the full covariance matrix for a single DES Y1). The green contours are marginalized over $\Omega_b$, $n_s$, $h_{100}$, redMaGC galaxy bias b as well as the shot-noise parameters $\alpha_0$ and $\alpha_1$. For the parameters $\Omega_b$, $n_s$, $h_{100}$ we have assumed the same flat priors as used in the DES Y1 combined probes analysis presented in Abbott et al. (in prep.). The red contours are marginalised only over bias and shot-noise parameters and the blue contours only vary $\Omega_m$ and $\sigma_8$. Even when going to this small parameter space, our model agrees with Buzzard within 1$\sigma$ errors of DES Y1. Bottom panel: Same contours but in the $\Sigma_8 - \Delta S_3/S_3$ plane and varying one additional parameter, $\Delta S_3/S_3$. Dotted lines show the true Buzzard cosmology and our fiducial value of $\Delta S_3/S_3 = 0$.</figcaption></figure>
[1] M. R. Becker et al. (Dark Energy Survey Collaboration), Phys. Rev. D **94**, 022002 (2016).

[2] J. Benjamin, C. Heymans, E. Semboloni, L. Van Waerbeke, H. Hoekstra, T. Erben, M. D. Gladders, M. Hetterscheidt, Y. Mellier, and H. K. C. Yee, Mon. Not. R. Astron. Soc. **381**, 702 (2007).

[3] L. Fu et al., Astron. Astrophys. **479**, 9 (2008).

[4] H. Hildebrandt et al., Mon. Not. R. Astron. Soc. **465**, 1454 (2017).

[5] M. J. Jee, J. A. Tyson, S. Hilbert, M. D. Schneider, S. Schmidt, and D. Wittman, Astrophys. J. **824**, 77 (2016).

[6] M. Kilbinger et al., Mon. Not. R. Astron. Soc. **430**, 2200 (2013).

[7] T. Schrabback et al., Astron. Astrophys. **516**, A63 (2010).

[8] M. A. Troxel et al., arXiv:1708.01538.

[9] L. Van Waerbeke et al., Astron. Astrophys. **358**, 30 (2000).

[10] D. M. Wittman, J. A. Tyson, D. Kirkman, I. Dell'Antonio, and G. Bernstein, Nature (London) **405**, 143 (2000).

[11] S. Alam et al., Mon. Not. R. Astron. Soc. **470**, 2617 (2017).

[12] M. Crocce et al., Mon. Not. R. Astron. Soc. **455**, 4301 (2016).

[13] J. Elvin-Poole et al., arXiv:1708.01536.

[14] T. G. Brainerd, R. D. Blandford, and I. Smail, Astrophys. J. **466**, 623 (1996).

[15] F. Brimioulle, S. Seitz, M. Lerchster, R. Bender, and J. Snigula, Mon. Not. R. Astron. Soc. **432**, 1046 (2013).

[16] J. Clampitt et al., Mon. Not. R. Astron. Soc. **465**, 4204 (2017).

[17] M. J. Hudson, S. D. J. Gwyn, H. Dahle, and N. Kaiser, Astrophys. J. **503**, 531 (1998).

[18] J. Prat et al., arXiv:1708.01537.

[19] G. Wilson, N. Kaiser, G. A. Luppino, and L. L. Cowie, Astrophys. J. **555**, 572 (2001).

[20] E. van Uitert, H. Hoekstra, M. Velander, D. G. Gilbank, M. D. Gladders, and H. K. C. Yee, Astron. Astrophys. **534**, A14 (2011).

[21] DES Collaboration, arXiv:1708.01530.

[22] R. Mandelbaum, A. Slosar, T. Baldauf, U. Seljak, C. M. Hirata, R. Nakajima, R. Reyes, and R. E. Smith, Mon. Not. R. Astron. Soc. **432**, 1544 (2013).

[23] E. van Uitert et al., Mon. Not. R. Astron. Soc. **476**, 4662 (2018).

[24] A. Barreira, S. Bose, B. Li, and C. Linares, J. Cosmol. Astropart. Phys. 02 (2017) 031.

[25] F. Bernardeau, L. van Waerbeke, and Y. Mellier, Astron. Astrophys. **322**, 1 (1997).

[26] S. Pires, A. Leonard, and J.-L. Starck, Mon. Not. R. Astron. Soc. **423**, 983 (2012).

[27] M. Takada and B. Jain, Mon. Not. R. Astron. Soc. **337**, 875 (2002).

[28] C. Uhlemann et al., Mon. Not. R. Astron. Soc. **473**, 5098 (2018).

[29] E. Semboloni, T. Schrabback, L. van Waerbeke, S. Vafaei, J. Hartlap, and S. Hilbert, Mon. Not. R. Astron. Soc. **410**, 143 (2011).

[30] T. Kacprzak et al., Mon. Not. R. Astron. Soc. **463**, 3653 (2016).

[31] C.-A. Lin and M. Kilbinger, Astron. Astrophys. **576**, A24 (2015).

[32] X. Liu et al., Mon. Not. R. Astron. Soc. **450**, 2888 (2015).

[33] A. B. Mantz, S. W. Allen, R. Glenn Morris, A. von der Linden, D. E. Applegate, P. L. Kelly, D. L. Burke, D. Donovan, and H. Ebeling, Mon. Not. R. Astron. Soc. **463**, 3582 (2016).

[34] C.-T. Chiang, C. Wagner, A. G. Sánchez, F. Schmidt, and E. Komatsu, J. Cosmol. Astropart. Phys. 09 (2015) 028.

[35] F. Simpson, J. B. James, A. F. Heavens, and C. Heymans, Phys. Rev. Lett. **107**, 271301 (2011).

[36] F. Simpson, A. F. Heavens, and C. Heymans, Phys. Rev. D **88**, 083510 (2013).

[37] F. Simpson, J. Harnois-Déraps, C. Heymans, R. Jimenez, B. Joachimi, and L. Verde, Mon. Not. R. Astron. Soc. **456**, 278 (2016).

[38] D. Gruen and O. Friedrich (The DES Collaboration), Mon. Not. R. Astron. Soc. **455**, 3367 (2016).

[39] D. Gruen, O. Friedrich et al., preceding paper, Phys. Rev. D **98**, 023507 (2018).

[40] K. M. Górski, E. Hivon, A. J. Banday, B. D. Wandelt, F. K. Hansen, M. Reinecke, and M. Bartelmann, Astrophys. J. **622**, 759 (2005).

[41] J. DeRose, R. Wechsler, E. Rykoff et al. (to be published).

[42] M. Bartelmann and P. Schneider, Phys. Rep. **340**, 291 (2001).

[43] H. S. Xavier, F. B. Abdalla, and B. Joachimi, Mon. Not. R. Astron. Soc. **459**, 3693 (2016).

[44] S. Hilbert, J. Hartlap, and P. Schneider, Astron. Astrophys. **536**, A85 (2011).

[45] R. E. Smith, J. A. Peacock, A. Jenkins, S. D. M. White, C. S. Frenk, F. R. Pearce, P. A. Thomas, G. Efstathiou, and H. M. P. Couchman, Mon. Not. R. Astron. Soc. **341**, 1311 (2003).

[46] R. Takahashi, M. Sato, T. Nishimichi, A. Taruya, and M. Oguri, Astrophys. J. **761**, 152 (2012).

[47] D. J. Eisenstein and W. Hu, Astrophys. J. **496**, 605 (1998).

[48] F. Bernardeau, Astron. Astrophys. **291**, 697 (1994).

[49] F. Bernardeau and P. Valageas, Astron. Astrophys. **364**, 1 (2000).

[50] P. Valageas, Astron. Astrophys. **382**, 412 (2002).

[51] F. Bernardeau, S. Colombi, E. Gaztañaga, and R. Scoccimarro, Phys. Rep. **367**, 1 (2002).

[52] A. Lue, R. Scoccimarro, and G. Starkman, Phys. Rev. D **69**, 044005 (2004).

[53] T. Multamäki, E. Gaztañaga, and M. Manera, Mon. Not. R. Astron. Soc. **344**, 761 (2003).

[54] B. Jain and P. Zhang, Phys. Rev. D **78**, 063503 (2008).

[55] A. Borisov and B. Jain, Phys. Rev. D **79**, 103506 (2009).

[56] D. Jeong and E. Komatsu, Astrophys. J. **703**, 1230 (2009).

[57] Planck Collaboration et al., Astron. Astrophys. **594**, A13 (2016).

[58] N. MacCrann, J. DeRose, R. Wechsler et al. (to be published).

[59] R. Wechsler, J. DeRose, M. Busha et al. (to be published).

[60] V. Springel, Mon. Not. R. Astron. Soc. **364**, 1105 (2005).

[61] M. Crocce, S. Pueblas, and R. Scoccimarro, Mon. Not. R. Astron. Soc. **373**, 369 (2006).

[62] C. Conroy, R. H. Wechsler, and A. V. Kravtsov, Astrophys. J. **647**, 201 (2006).

[63] R. M. Reddick, R. H. Wechsler, J. L. Tinker, and P. S. Behroozi, Astrophys. J. **771**, 30 (2013).

[64] M. C. Cooper et al., Astrophys. J. Suppl. Ser. **193**, 14 (2011).

[65] S. Miyazaki et al., Publ. Astron. Soc. Jpn. **54**, 833 (2002).

[66] M. R. Becker, Mon. Not. R. Astron. Soc. **435**, 115 (2013).
[67] E. Huff and R. Mandelbaum, arXiv:1702.02600.

[68] E. S. Sheldon and E. M. Huff, Astrophys. J. **841**, 24 (2017).

[69] J. Zuntz *et al.* (to be published).

[70] N. Benítez, Astrophys. J. **536**, 571 (2000).

[71] B. Hoyle *et al.*, Mon. Not. R. Astron. Soc. (2018), DOI: 10.1093/mnras/sty957.

[72] F. Bernardeau, C. Pichon, and S. Codis, Phys. Rev. D **90**, 103519 (2014).

[73] F. Bernardeau, S. Codis, and C. Pichon, Mon. Not. R. Astron. Soc. **449**, L105 (2015).

[74] S. Codis, C. Pichon, F. Bernardeau, C. Uhlemann, and S. Prunet, Mon. Not. R. Astron. Soc. **460**, 1549 (2016).

[75] S. Codis, F. Bernardeau, and C. Pichon, Mon. Not. R. Astron. Soc. **460**, 1598 (2016).

[76] D. N. Limber, Astrophys. J. **117**, 134 (1953).

[77] L. Clerkin *et al.*, Mon. Not. R. Astron. Soc. **466**, 1444 (2017).

[78] P. Fosalba and E. Gaztanaga, Mon. Not. R. Astron. Soc. **301**, 503 (1998).

[79] T. Baldauf, U. Seljak, R. E. Smith, N. Hamaus, and V. Desjacques, Phys. Rev. D **88**, 083507 (2013).

[80] A. Dekel and O. Lahav, Astrophys. J. **520**, 24 (1999).

[81] J. Clampitt *et al.*, Mon. Not. R. Astron. Soc. **465**, 4204 (2017).

[82] R. Mandelbaum, U. Seljak, G. Kauffmann, C. M. Hirata, and J. Brinkmann, Mon. Not. R. Astron. Soc. **368**, 715 (2006).

[83] Y.-j. Peng, S. J. Lilly, A. Renzini, and M. Carollo, Astrophys. J. **757**, 4 (2012).

[84] M. Cacciato, O. Lahav, F. C. van den Bosch, H. Hoekstra, and A. Dekel, Mon. Not. R. Astron. Soc. **426**, 566 (2012).

[85] A. Dvornik *et al.*, arXiv:1802.00734.

[86] V. Mukhanov, Physical Foundations of Cosmology **1**, 279 (2005).

[87] R. Takahashi, S. Soma, M. Takada, and I. Kayo, Mon. Not. R. Astron. Soc. **444**, 3473 (2014).

[88] V. Springel *et al.*, Nature (London) **435**, 629 (2005).

[89] O. Friedrich, S. Seitz, T. F. Eifler, and D. Gruen, Mon. Not. R. Astron. Soc. **456**, 2662 (2016).