# Density Split Statistics: Cosmological Constraints from Counts and Lensing in Cells in DES Y1 and SDSS Data

D. Gruen, O. Friedrich, E. Krause, J. Derose, R. Cawthon, C. Davis, J.

Elvin-Poole, E.S. Rykoff, R.H. Wechsler, A. Alarcon, et al.

► To cite this version:

D. Gruen, O. Friedrich, E. Krause, J. Derose, R. Cawthon, et al. Density Split Statistics: Cosmological Constraints from Counts and Lensing in Cells in DES Y1 and SDSS Data. Physical Review D, 2018, 98 (2), pp.023507. 10.1103/PhysRevD.98.023507. hal-01846741

HAL Id: hal-01846741

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

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: Cosmological constraints from counts and lensing in cells in DES Y1 and SDSS data

D. Gruen,¹,²* O. Friedrich,³,⁴ E. Krause,¹ J. DeRose,⁵,¹ R. Cawthon,⁶ C. Davis,¹ J. Elvin-Poole,⁷ E. S. Rykoff,¹,²

R. H. Wechsler,⁵,¹,² A. Alarcon,⁸ G. M. Bernstein,⁹ J. Blazek,¹⁰,¹¹ C. Chang,⁶ J. Clampitt,⁹ M. Crocce,⁸ J. De Vicente,¹²

M. Gatti,¹³ M. S. S. Gill,² W. G. Hartley,¹⁴,¹⁵ S. Hilbert,¹⁶,¹⁷ B. Hoyle,⁴,³ B. Jain,⁹ M. Jarvis,⁹ O. Lahav,¹⁴ N. MacCrann,¹⁰,¹⁸

T. McClintock,¹⁹ J. Prat,¹³ R. P. Rollins,⁷ A. J. Ross,¹⁰ E. Rozo,¹⁹ S. Samuroff,⁷ C. Sánchez,¹³ E. Sheldon,²⁰

M. A. Troxel,¹⁰,¹⁸ J. Zuntz,²¹ T. M. C. Abbott,²² F. B. Abdalla,¹⁴,²³ S. Allam,²⁴ J. Annis,²⁴ K. Bechtol,²⁵

A. Benoit-Lévy,²⁶,¹⁴,²⁷ E. Bertin,²⁶,²⁷ S. L. Bridle,⁷ D. Brooks,¹⁴ E. Buckley-Geer,²⁴ A. Carnero Rosell,²⁸,²⁹

M. Carrasco Kind,³⁰,³¹ J. Carretero,¹³ C. E. Cunha,¹ C. B. D'Andrea,⁹ L. N. da Costa,²⁸,²⁹ S. Desai,³² H. T. Diehl,²⁴

J. P. Dietrich,¹⁷,¹⁶ P. Doel,¹⁴ A. Drlica-Wagner,²⁴ E. Fernandez,¹³ B. Flaugher,²⁴ P. Fosalba,⁸ J. Frieman,²⁴,⁶

J. García-Bellido,³³ E. Gaztanaga,⁸ T. Giannantonio,³⁴,³⁵,³ R. A. Gruendl,³⁰,³¹ J. Gschwend,²⁸,²⁹ G. Gutierrez,²⁴

K. Honscheid,¹⁰,¹⁸ D. J. James,³⁶ T. Jeltema,³⁷ K. Kuehn,³⁸ N. Kuropatkin,²⁴ M. Lima,³⁹,²⁸ M. March,⁹ J. L. Marshall,⁴⁰

P. Martini,¹⁰,⁴¹ P. Melchior,⁴² F. Menanteau,³⁰,³¹ R. Miquel,⁴³,¹³ J. J. Mohr,¹⁷,¹⁶,⁴ A. A. Plazas,⁴⁴ A. Roodman,¹,²

E. Sanchez,¹² V. Scarpine,²⁴ M. Schubnell,⁴⁵ I. Sevilla-Noarbe,¹² M. Smith,⁴⁶ R. C. Smith,²²

M. Soares-Santos,²⁴ F. Sobreira,⁴⁷,²⁸ M. E. C. Swanson,³¹ G. Tarle,⁴⁵ D. Thomas,⁴⁸ V. Vikram,⁴⁹

A. R. Walker,²² J. Weller,¹⁷,⁴,³ and Y. Zhang²⁴

(DES Collaboration)

¹Kavli Institute for Particle Astrophysics & Cosmology, P. O. Box 2450, Stanford University, Stanford, California 94305, USA

²SLAC National Accelerator Laboratory, Menlo Park, California 94025, USA

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

⁴Max Planck Institute for Extraterrestrial Physics, Giessenbachstrasse, 85748 Garching, Germany

⁵Department of Physics, Stanford University, 382 Via Pueblo Mall, Stanford, California 94305, USA

⁶Kavli Institute for Cosmological Physics, University of Chicago, Chicago, Illinois 60637, USA

⁷Jodrell Bank Center for Astrophysics, School of Physics and Astronomy, University of Manchester, Oxford Road, Manchester, M13 9PL, United Kingdom

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

⁹Department of Physics and Astronomy, University of Pennsylvania, Philadelphia, Pennsylvania 19104, USA

¹⁰Center for Cosmology and Astro-Particle Physics, The Ohio State University, Columbus, Ohio 43210, USA

¹¹Institute of Physics, Laboratory of Astrophysics, École Polytechnique Fédérale de Lausanne (EPFL), Observatoire de Sauverny, 1290 Versailles, Switzerland

¹²Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas (CIEMAT), Madrid, Spain

¹³Institut de Física d'Altes Energies (IFAE), The Barcelona Institute of Science and Technology, Campus UAB, 08193 Bellaterra (Barcelona) Spain

¹⁴Department of Physics & Astronomy, University College London, Gower Street, London, WC1E 6BT, United Kingdom

¹⁵Department of Physics, ETH Zurich, Wolfgang-Pauli-Strasse 16, CH-8093 Zurich, Switzerland

¹⁶Faculty of Physics, Ludwig-Maximilians-Universität, Scheinerstr. 1, 81679 Munich, Germany

¹⁷Excellence Cluster Universe, Boltzmannstr. 2, 85748 Garching, Germany

¹⁸Department of Physics, The Ohio State University, Columbus, Ohio 43210, USA

¹⁹Department of Physics, University of Arizona, Tucson, Arizona 85721, USA

²⁰Brookhaven National Laboratory, Bldg 510, Upton, New York 11973, USA

²¹Institute for Astronomy, University of Edinburgh, Edinburgh EH9 3HJ, United Kingdom

²²Cerro Tololo Inter-American Observatory, National Optical Astronomy Observatory, Casilla 603, La Serena, Chile

²³Department of Physics and Electronics, Rhodes University, PO Box 94, Grahamstown 6140, South Africa

²⁴Fermi National Accelerator Laboratory, P. O. Box 500, Batavia, Illinois 60510, USA

²⁵LSST, 933 North Cherry Avenue, Tucson, Arizona 85721, USA
<sup>26</sup>CNRS, UMR 7095, Institut d'Astrophysique de Paris, F-75014 Paris, France

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

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

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

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

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

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

<sup>33</sup>Instituto de Física Teórica UAM/CSIC, Universidad Autónoma de Madrid, 28049 Madrid, Spain

<sup>34</sup>Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, United Kingdom

<sup>35</sup>Kavli Institute for Cosmology, University of Cambridge, Madingley Road, Cambridge CB3 0HA, United Kingdom

<sup>36</sup>Astronomy Department, University of Washington, Box 351580, Seattle, Washington 98195, USA

<sup>37</sup>Santa Cruz Institute for Particle Physics, Santa Cruz, California 95064, USA

<sup>38</sup>Australian Astronomical Observatory, North Ryde, New South Wales 2113, Australia

<sup>39</sup>Departamento de Física Matemática, Instituto de Física, Universidade de São Paulo, CP 66318, São Paulo, SP, 05314-970, Brazil

<sup>40</sup>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

<sup>41</sup>Department of Astronomy, The Ohio State University, Columbus, Ohio 43210, USA

<sup>42</sup>Department of Astrophysical Sciences, Princeton University, Peyton Hall, Princeton, New Jersey 08544, USA

<sup>43</sup>Institució Catalana de Recerca i Estudis Avançats, E-08010 Barcelona, Spain

<sup>44</sup>Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, California 91109, USA

<sup>45</sup>Department of Physics, University of Michigan, Ann Arbor, Michigan 48109, USA

<sup>46</sup>School of Physics and Astronomy, University of Southampton, Southampton, SO17 1BJ, United Kingdom

<sup>47</sup>Instituto de Física Gleb Wataghin, Universidade Estadual de Campinas, 13083-859, Campinas, SP, Brazil

<sup>48</sup>Institute of Cosmology & Gravitation, University of Portsmouth, Portsmouth, PO1 3FX, United Kingdom

<sup>49</sup>Argonne National Laboratory, 9700 South Cass Avenue, Lemont, Illinois 60439, USA

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

(Received 17 October 2017; published 13 July 2018)

We derive cosmological constraints from the probability distribution function (PDF) of evolved large-scale matter density fluctuations. We do this by splitting lines of sight by density based on their count of tracer galaxies, and by measuring both gravitational shear around and counts-in-cells in overdense and underdense lines of sight, in Dark Energy Survey (DES) First Year and Sloan Digital Sky Survey (SDSS) data. Our analysis uses a perturbation theory model [O. Friedrich et al., Phys. Rev. D **98**, 023508 (2018)] and is validated using N-body simulation realizations and log-normal mocks. It allows us to constrain cosmology, bias and stochasticity of galaxies with respect to matter density and, in addition, the skewness of the matter density field. From a Bayesian model comparison, we find that the data weakly prefer a connection of galaxies and matter that is stochastic beyond Poisson fluctuations on $\le 20$ arcmin angular smoothing scale. The two stochasticity models we fit yield DES constraints on the matter density $\Omega_m = 0.26^{+0.04}_{-0.03}$ and $\Omega_m = 0.28^{+0.05}_{-0.04}$ that are consistent with each other. These values also agree with the DES analysis of galaxy and shear two-point functions (3x2pt, DES Collaboration et al.) that only uses second moments of the PDF. Constraints on $\sigma_8$ are model dependent ($\sigma_8 = 0.97^{+0.07}_{-0.06}$ and $0.80^{+0.06}_{-0.07}$ for the two stochasticity models), but consistent with each other and with the 3 x 2pt results if stochasticity is at the low end of the posterior range. As an additional test of gravity, counts and lensing in cells allow to compare the skewness $S_3$ of the matter density PDF to its $\Lambda$CDM prediction. We find no evidence of excess skewness in any model or data set, with better than 25 per cent relative precision in the skewness estimate from DES alone.

DOI: 10.1103/PhysRevD.98.023507

<footnote>*Einstein Fellow.
dgruen@stanford.edu</footnote>
# I. INTRODUCTION

Measurements of the two-point correlation function of the evolved matter density field have provided competitive constraints on fundamental cosmological parameters. In combination with cosmic microwave background (CMB) and other geometric data, they are stringent tests of $\Lambda$CDM predictions on the evolution of structure over cosmic time [1–5]. Ostensibly, larger studies of this kind are the primary goal of the upcoming ambitious ground-based and space-based surveys by *Euclid*, LSST, and WFIRST.

On a given smoothing scale, one can describe a field locally by its PDF. An example of this is the PDF of fluctuations of the mean matter density inside spherical or cylindrical volumes. The variance, or second moment, of the PDF is measured by two-point statistics. For Gaussian distributions, this captures all the information in all the moments of the PDF. But when the field is non-Gaussian, the third moment (skewness) can take any value—therefore it contains information that is *not contained* in two-point statistics.

Unlike the primordial CMB, which is extremely close to a Gaussian random field, the density distribution in the evolved Universe has been driven away from Gaussianity by gravitational collapse. Third and higher order moments arise on any scale. Two-point measurements are therefore inherently very incomplete pictures of the matter density field. Even the full hierarchy of N-point correlations ceases to fully describe its statistics [6–8]. This is unfortunate in two ways: (1) A lot of information on cosmology, and (2) a lot of opportunities to test additional, independent $\Lambda$CMB predictions of the growth of structure beyond its variance, are lost by looking at two-point functions alone.

There are other reasons that make two-point correlations a somewhat blunt tool. First, the information to be gained from galaxy auto- or cross-correlations must be related to the clustering of matter by a *bias model*, i.e., a description of how galaxies trace matter density. Yet the information on the bias model that is available from two-point functions alone is limited. The primary reason for the success of joint probes is that they can partially break these degeneracies. For instance, the joint analysis of galaxy clustering and galaxy-galaxy lensing can constrain two combinations of $\sigma_8$, galaxy bias, and galaxy stochasticity. As one pushes to smaller scales where a lot of the cosmological constraining power resides, a linear bias model without stochasticity is not sufficient. The resulting degeneracies thus largely annihilate the information that is gained. Second, the information on the variance of the matter density field can be used only if that variance can be modeled—on nonlinear scales, complex physics that involve baryons and neutrinos begin to influence any moment of the matter density field. Using two-point functions, these complex physics can be constrained [9,10], although (for the same reason as for the bias model) only with limited discriminating power. If we could recover small scale information with models that can be trusted, the ability of presently and

imminently available data sets to confirm or reduce tensions between the CMB and evolved power spectrum would immediately be boosted.

For these reasons, studies of the cosmic density PDF, which address these problems from a different and complementary direction, have gained interest over the last years. The full shape of the joint matter and galaxy density PDF depends on moments of the matter density field and parameters of the bias model that are degenerate in correlation function measurements. Numerical simulations [11–13], tree-level perturbation calculations [14,15], and extensions of theory beyond that [16], have been shown to provide accurate predictions for the matter density PDF. Parameter forecasts show that PDF measurements on data are promising [17–19] due to the complementary information, different degeneracies, and different dependence on observational systematics—factor-of-two improvements in constraining power can be achieved in joint measurements of PDF and two-point functions. While the galaxy count PDF alone can be used to break degeneracies of cosmological and bias parameters [20], gravitational lensing greatly complements this by measuring the actual matter density PDF. Practical application of shear PDF statistics to data has been made with DES [21,22], yet so far with limited use for quantitative constraints.

In this paper, we use the smoothed, joint, projected galaxy count and matter density PDF to constrain cosmological parameters and a galaxy bias model. Our basic concept is to (1) split the sky by the count of tracer galaxies in a top-hat aperture and extended redshift range into quantiles of density, and to (2) measure the gravitational shear around each of the quantiles to reconstruct the matter density PDF. These measurements are a generalization of trough lensing, introduced in Gruen et al. [23] (see also [24–26]). They are also closely related to the galaxy-matter aperture statistics of Simon et al. [27]. We make them on Dark Energy Survey Year 1 (DES Y1) and SDSS DR8 data.

The measurements are analyzed with a tree-level perturbation theory prediction for the joint statistical properties of lensing convergence and density contrast and galaxy bias models of varying complexity (cf. our companion paper Friedrich et al.). We use the analysis not just to provide an independent measurement of cosmological parameters, but also to confront the $\Lambda$CDM prediction for the *skewness* of the matter density field with data, in a model independent test of structure formation. That is, we measure the asymmetry of the low and high density tails of the distribution of matter density in the Universe. Two-point statistics, which only measure the *width* of the matter density distribution, discard this information.

This paper is structured as follows. We describe the data we use and our measurement methodology in Sec. II. Our modeling of these measurements, based on Friedrich et al. [28], is summarized in Sec. III. The covariance matrix we estimate is described in Sec. IV. We combine
measurements, model, and covariance into an inference framework in Sec. V. Results are presented in Sec. VI, and we conclude in Sec. VII. Several tests and technical aspects of this work are detailed in the Appendix.

## II. MEASUREMENT

In the following section, we first describe our method of splitting lines of sight by density based on counts of a tracer galaxy sample (Sec. II A). The REDMAGIC tracer catalogs we use to do this in DES and SDSS and their redshift distribution calibration are presented in Sec. II B 1 and Sec. II C 1. Details on our DES and SDSS source shape and photometric redshift catalogs are given in Sec. II B 2 and Sec. II C 2. The measurement of shear and counts-in-cells signals is described in Sec. II D.

### A. Splitting the sky by density

The basic idea of this study is to split the sky into lines of sight of different density.

To this end, we use a sample of foreground galaxies as tracers of the matter field (the REDMAGIC galaxies at $0.2 < z_T < 0.45$ described in Sec. II B 1 and II C 1). We count these galaxies within circular top-hat apertures with a range of radii $\theta_T = 10', 20', 30', 60'$, centered on a regular HEALPIX [29] grid of $N_{side} = 1024$ (3.4 arcmin grid spacing).

We then assign each line of sight to one of five density quintiles by sorting all lines of sight by galaxy count. The 20 percent of lines of sight with the lowest galaxy count are what we will call the lowest density quintile 1 (or troughs, cf. [23]). The 20 percent of lines of sight with the highest galaxy count (quintile 5) we will denote as overdense lines of sight.

Compared to Gruen et al. [23], we apply a more elaborate scheme of accounting for varying fractions of masked area within the respective survey region. A mask accompanying the REDMAGIC [30] catalog that we will use as our tracers (Sec. II B 1) indicate what fraction of the area inside each pixel in a $N_{side} = 4096$ HEALPIX map is covered by DES Y1 Gold photometry [31] to sufficient depth for detecting REDMAGIC galaxies out to at least $z = 0.45$. For each line of sight, we estimate the fraction of masked area $f_{mask}$ within the corresponding top-hat aperture from the REDMAGIC masks. Centers with more than a fraction $f_{mask}^{max}$ of area within the aperture lost to masking are discarded.

The depth of SDSS is very uniform, with the REDMAGIC sample being complete to $z = 0.45$ everywhere. In this case, we use $f_{mask,SDSS}^{max} = 0.1$. Despite its greater overall depth, DES Y1 [32] is generally more inhomogeneous than the final SDSS imaging data. Where the REDMAGIC sample is not complete to $z = 0.45$, we remove all tracer galaxies and define the area to be fully masked. Due to the larger fraction of masked area, we use $f_{mask,DES}^{max} = 0.2$, above which we discard lines of sight.

To account for residual differences in $f_{mask}$ we apply the following probabilistic scheme of quintile assignment. For each line of sight $i$ with masking fraction $f_{mask,i}$ and raw tracer galaxy count $N_{raw,i}$, we define $N_i$ as a draw from a binomial distribution with $N_{raw,i}$ repetitions and success probability $p_i = 1 - (f_{mask}^{max} - f_{mask,i})$,

$$P(N_i|N_{raw,i}, f_{mask,i}) = \binom{N_{raw,i}}{N_i} p_i^{N_i} (1-p_i)^{N_{raw,i}-N_i} \quad (2.1)$$

This emulates the masking of a fixed fraction $f_{mask}^{max}$ of area within each aperture. It preserves the expectation value of galaxy count in an aperture, regardless of its masking fraction. Under the assumption that galaxies or masked pixels do not cluster, and that galaxy count is not stochastic beyond Poissonian noise, this masking procedure would preserve the full distribution of galaxy counts at fixed matter density (see Appendix A). The latter conditions are not true in practice, which is why the degree and spatial distribution of masking still affects the width of $P(N_i)$ at fixed expectation value. Tests of likelihood runs and the masked $P(N)$ in the Buzzard simulations (see Appendix E and [28]) indicate that this is not a major concern for our analysis.

We assign a line of sight $i$ to a density quintile based on many random realizations of $N_i$. Different realizations of $N_i$ can cause different quintile assignments. To account for this, we define a weight $w_i^q$, proportional to the number of times $N_i$ is in quintile $q$. This weight is assigned to line of sight $i$ when measuring the signal, e.g., the mean tangential shear, of quintile $q$.

Figure 1 shows the result of this quintile assignment procedure for the joint region covered by DES Y1 and SDSS.

### B. Dark Energy Survey Y1 data

The Dark Energy Survey data we use in this work is from the SPT region of the first year of science observations (Y1) performed between 31 August 2013 and 9 February 2014. Details of the data and photometric pipeline are described in Drlica-Wagner et al. [31].

We use catalogs of luminous red galaxies (REDMAGIC galaxies) as tracers of the foreground matter density field and galaxy shape and photometric redshift catalogs for measuring its gravitational shear signal, all of which are described briefly below and in detail in Elvin-Poole et al. [30], Hoyle et al. [33], Zuntz et al. [34].

In all likelihood analyses run on data in this work, we propagate the three most relevant calibration uncertainties of these catalogs:

(i) the multiplicative bias of the shear signal, characterized as $m = \gamma^{obs}/\gamma^{true} - 1$,

(ii) the bias in mean redshift of each source bin $i$, characterized by a $\Delta z_s^i$ which we use to evaluate
<figure><img src="image_3.png" /><figcaption>FIG. 1. Overdense (red) and underdense (cyan) lines of sight in the DES Y1 (south/bottom) and SDSS (north/top) survey regions. Five quintiles of density of the $z_T = 0.2–0.45$ REDMAGIC tracer galaxy sample smoothed within a $\theta_T = 30'$ radius are shown in the same color scheme as in Fig. 3. Pixels are drawn in the color of the quintile with the highest probability/weight $w^q$ (see Sec. II A). Graticule shows lines of $\Delta RA, \Delta dec = \pm 30^\circ$, centered on (RA, dec) = (0, 0).</figcaption></figure>

$n_s^i(z) = n_s^{i,PZ}(z - \Delta z_s^i)$, where $n_s^{i,PZ}$ is the photo-metric estimate of the source redshift distribution, and

(iii) the bias in mean redshift of the tracer galaxy sample, characterized by a $\Delta z_l$ which we use to evaluate $n_l(z) = n_l^{\text{redMaGiC}}(z - \Delta z_l)$.

The derivation of priors on these calibration uncertainties is described or referenced in Sec. V B.

## 1. Tracer catalog

The REDMAGIC [35] algorithm identifies a sample of red galaxies with constant comoving density and fixed luminosity threshold. This is done by fitting the DES photometry of each galaxy in the survey to find its maximum likelihood luminosity and redshift under the assumption of the redMaPPer [36] red sequence template. Galaxies are removed from the REDMAGIC catalog if their fitted luminosity falls below a threshold ($0.5L^*$ for the high density run used in this work). The catalog is further pruned to retain a fixed number density of galaxies per comoving volume element, keeping those that are best fit (in terms of photometric $\chi^2$) by the red sequence template. The resulting galaxy density is $10^{-3}h^3\,\text{Mpc}^{-3}$ in the case of the high density catalog.

This procedure was run on two different photometric measurements of DES Y1 galaxies, one with the SEXTRACTOR MAG_AUTO method and one performing a joint fit to the multiepoch data of multiple overlapping

objects (MOF). Potential correlation of the surface density of REDMAGIC galaxies with observational systematics in DES Y1 have been extensively tested in Elvin-Poole et al. [30] for both versions of the catalog. They found that in the redshift range used for the tracer galaxies, the MAG_AUTO version of the REDMAGIC catalog shows smaller correlations with observational systematics.

We hence adopt MAG_AUTO REDMAGIC with high density as our fiducial tracer catalog. In a trade-off of signal and noise, we choose $z_T = 0.2–0.45$ as the tracer redshift range. We derive weights for the correction of REDMAGIC density for the effect of systematics as in Elvin-Poole et al. [30]. We find significant correlations of REDMAGIC density with r band exposure time and seeing, and with i band sky brightness. In the algorithm described in Sec. II A, we have applied these by dividing the fraction of good area in each pixel by the systematics weight that decorrelates REDMAGIC density with these survey properties.

We do, however, test whether the choice of photometry pipeline (MAG_AUTO or MOF) and the choice of whether we apply the systematics weight in our density splitting procedure makes a difference to our analysis. These tests are detailed in Appendix C and show that the effect on the amplitude of our measured signals is negligible.

The redshift distribution of the tracer galaxy population, estimated by convolving the photometric redshift of each REDMAGIC galaxy with its error estimate $\sigma_z \approx 0.017 \times (1+z)$ [30], is shown as the grey contour in Fig. 2.
<figure><img src="image_4.png" /><figcaption>FIG. 2. The redshift distributions of the DES source redshift bins (violet, red, yellow dashed lines) and SDSS sources (blue, dotted line), and of the REDMAGIC tracer galaxies (grey shaded area) including scatter in their photometric redshift estimates. The dependence of predicted amplitude of the trough lensing signal on source redshift is shown by the black curve. Grey marks on upper axis indicate nominal redshift range of tracers galaxies $z_T = 0.2–0.45$. Lines are normalized to match maxima and the trough signal is evaluated at $\theta = 2\theta_T$, although the dependence of source redshift scaling on angular distance is minor.</figcaption></figure>

Note that due to scatter in photo-z this extends beyond the redshift range $z_T = 0.2–0.45$ inside which these galaxies were selected.

As for other uses of REDMAGIC for cosmology [1,30,37], we limit the catalog to the contiguous DES-SPT area of 1321 deg².

## 2. Lensing source catalogs

Detailed descriptions and tests of the DES Y1 lensing source catalogs are presented in Zuntz et al. [34], Troxel et al. [4] and Prat et al. [37], and the redshift distributions of source galaxies are estimated and calibrated in Davis et al. [38], Gatti et al. [39], Hoyle et al. [33]. We only give a brief summary of the two independent shape catalogs from DES Y1 here.

The fiducial catalog with the larger number of source galaxies is based on the METACALIBRATION method [40,41]. In this scheme, a Gaussian, convolved with the individual exposure point-spread function, is fit jointly to all single-epoch *r*, *i*, and *z*-band images of each galaxy. Galaxies are selected by the size and signal-to-noise ratio of the best fit, and the ellipticity of the Gaussian is used as an estimate of shear. Multiplicative biases in mean shear are caused by both the galaxy selection (selection bias) and the use of a maximum likelihood estimator with a simplified model (noise and model bias). In METACALIBRATION, these are calibrated and removed using a repetition of the Gaussian fit on versions of the galaxy images that have been artificially sheared by a known amount.

As a second catalog, we use IM3SHAPE, which produces a maximum likelihood estimate of shear based on a bulge or

disc fit to all DES Y1 *r* band images of each galaxy. Multiplicative biases in these estimators are calibrated using realistic image simulations of DES Y1 [34,42].

Our estimator of tangential shear around overdense and underdense lines of sight, including the bias corrections, is described in Sec. II D 1. We use the galaxy selection criteria recommended in Zuntz et al. [34]. We split galaxies into redshift bins using the mean *z* of the individual galaxy *p*(*z*) as estimated by the Bayesian Photometric Redshift code (BPZ) [33,43]. We note that for the METACALIBRATION catalog, we run BPZ on METACALIBRATION measurements of galaxy fluxes (both on the original and artificially sheared images) to be able to correct for photo-*z* related shear selection biases (see also Sec. 3.3 of Hoyle et al. [33] and Sec. IV.A.1 of Prat et al. [37]). The three source redshift bins we use are identical to the three highest redshift bins of Hoyle et al. [33], i.e., with sources at mean *z* = 0.43–0.60, 0.60–0.93, 0.93–1.30. Their redshift distributions, as estimated by BPZ using MOF photometry, are shown in Fig. 2.

Uncertainties on residual multiplicative shear bias and on the mean values of the binned redshift distributions [33,34,38,39] are marginalized over in our analysis (see Sec. V B).

## C. SDSS DR8 data

### 1. REDMAGIC tracer catalog

The tracer population in SDSS is the REDMAGIC [35] high density sample, selected by SDSS photometry and cut to the same redshift range $z_T = 0.2–0.45$. Despite this similarity, we will not assume in this work that SDSS and DES REDMAGIC galaxies are the exact same populations.

SDSS has the benefit of an overlapping sample of galaxies with spectroscopic redshifts. We use this to calibrate the mean of the redshift distribution with clustering redshifts, independent of the photometric estimate, in Appendix F 1. We find no significant bias, yet marginalize over the uncertainty in the analyses presented herein (see Sec. V B).

### 2. Lensing source catalogs

We use the shape and photometric redshift $p(z)$ catalog of Sheldon et al. [44] with minor modifications, identical to those in Clampitt and Jain [45]. We refer to these papers for details, but describe our source selection and priors on systematic uncertainties of shears and photometric redshifts below.

Due to the lower observational depth, the SDSS shape catalog peaks at much lower redshift than the one from DES Y1. The source redshift dependence of the trough lensing signal (cf. Fig. 2) and complications arising from significant overlap of sources with the tracer redshift range lead us to only use sources with a mean redshift estimate of $0.45 \le z < 1.0$. We split these into four bins of $z = 0.45–0.5, 0.5–0.55, 0.55–0.6$ and $0.6–1.0$. Within
each of these bins, each individual source is assigned a minimum-variance relative weight (cf. [45]) of

$$w_i = [\sigma_{i,shape,meas}^2 + 0.32^2]^{-1}. \quad (2.2)$$

Due to the moderate signal-to-noise ratio and redshift range of sources, we combine the four source redshift bins into one for the purpose of our final data vector. In this, we apply an optimal relative weighting of the bins as follows.

The predicted amplitude of shear around our troughs at $z_T = 0.2-0.45$ (see black line in Fig. 2) scales with source redshift approximately as the amplitude of gravitational shear $\Sigma_{crit}^{-1}$ [see Eq. (2.4)] due to a lens at $z_d = 0.36$. We use the value of $\langle\Sigma_{crit}^{-1}\rangle$ estimated for $z_d = 0.36$ and the stacked $p(z)$ of each of these four bins to apply a relative weight of $W_{bin} = 1, 1.30, 1.56$, and $1.81$ to each of them. Because the number density of sources is steeply falling with source redshift in this range, the effective total relative weights of the four bins (equal to this $W_{bin}$ times the sum of all source $w_i$) are 1,0.697,0.453, and 0.202. We use these effective weights to combine the measured shear signal and the $n_s(z)$ from each of the four bins into a single source sample.

As a calibration of the photometric estimate, the mean redshift of the sources is constrained by their angular cross-correlation with galaxies with known spectroscopic redshift (Appendix F 2).

## D. Measured signals

Our data vector in this work contains two components, the modeling of which was extensively tested in Friedrich et al. [28]. In Sec. II D 1, we describe the measurement of gravitational shear signals around overdense and underdense lines of sight. Section II D 2 details our measurement of mean counts-in-cells in each density quintile in the presence of masking.

All measurements are made in jackknife resamplings of the survey. The covariance model constructed in Sec. IV can therefore be compared to an jackknife covariance. These were made based on 100 and 200 patches in the DES and SDSS footprint, respectively, defined by k-means clustering<fnref n="1" /> of the tracer galaxies, an algorithm that splits the tracer galaxies into spatially compact subsets by their distance to the nearest among a set of centers, optimizing the center positions to minimize these distances.

### 1. Shear

The ellipticity of a galaxy is a pseudovector with two components, $e_1$ and $e_2$ that, for any lens position, can equivalently be described by a component tangential to a circle around a lens ($e_t$) and by a component rotated by $\pi/4$ relative to that ($e_\times$).

Gravitational shear due to any single lens only affects the mean component of $e_t$ for an ensemble of sources sampling a full annulus around the lens. As a function of angular separation $\theta$ from the lens, this effect is described by the tangential shear profile $\gamma_t(\theta)$. For a single lens, the tangential shear profile is directly related to the azimuthally averaged, projected surface mass density $\Sigma(\theta)$ of the lens, i.e., the projected mass per physical area, as

$$\gamma_t(\theta) = [\langle\Sigma\rangle(<\theta) - \Sigma(\theta)] \times \Sigma_{crit}^{-1} \equiv \langle\kappa\rangle(<\theta) - \kappa(\theta), \quad (2.3)$$

where, in a flat universe,

$$\Sigma_{crit}^{-1} = \frac{4\pi G \chi_d (\chi_s - \chi_d)}{c^2 \chi_s (1 + z_d)} \quad (2.4)$$

is the inverse of the critical surface mass density and $\chi_{d,s}$ is the comoving distances to the deflector at redshift $z_d$ and the lensed source, respectively. For a set of lenses along the line of sight, the signal on any source is close to the sum of the effects of all lenses. One can still define a convergence $\kappa$ related to mean gravitational shear as in Eq. (2.3), although it is no longer relatable to a uniformly weighted surface mass density [see Eq. (3.7)].

The relation between tangential shear and measured tangential ellipticity is less straightforward and depends on the implementation of the selection and measurement of source ellipticities. For the two schemes used on DES Y1, the responsivity $R = d\langle e_t \rangle/d\gamma_t$ of observed ellipticity to applied tangential shear is calibrated very differently: for METACALIBRATION, it is estimated from versions of the actual galaxy images sheared with image manipulation algorithms; for IM3SHAPE, it is estimated from realistic simulations of DES imaging data (and usually defined as $m = R - 1$). Both types of calibration contain an explicit or implicit correction for selection biases, i.e., the shear dependence of the choice of whether to include a galaxy in the source sample.

For METACALIBRATION, we define the estimator $\hat{\gamma}_t^q$ of mean tangential shear around lines of sight $i$ with probability $w_i^q$ to be in a given density quintile $q$ as

$$\hat{\gamma}_t^q = \hat{\gamma}_t^{\text{q,signal}} - \hat{\gamma}_t^{\text{random}} = \frac{\sum_{i,j} w_i^q e_{ij,t}}{R \sum_{i,j} w_i^q} - \frac{\sum_{i,j} e_{ij,t}}{R \sum_{i,j} 1}, \quad (2.5)$$

where $e_{ij,t}$ is the ellipticity of source $j$ in the tangential direction around line of sight $i$, the sums run over all lines of sight $i$ in the mask of the density-split sky and all sources $j$ in an angular bin around each line of sight. The second term subtracts shear around random lines of sight—for our statistic, these are all healpix pixels around which the masked fraction of area $f_{mask} < f_{mask}^{\text{max}}$ (cf. Sec. II A). $R$ is the sum of shear and selection responsivity,

<footnote n="1">¹https://github.com/esheldon/kmeans-radec.</footnote>
$$R = R_{\gamma} + R_{S} = \frac{1}{2} \left\langle \frac{e_{1}^{+} - e_{1}^{-}}{2\Delta\gamma_{1}} + \frac{e_{2}^{+} - e_{2}^{-}}{2\Delta\gamma_{2}} \right\rangle + \frac{1}{2} \left( \frac{\langle e_{1} \rangle^{+} - \langle e_{1} \rangle^{-}}{2\Delta\gamma_{1}} + \frac{\langle e_{2} \rangle^{+} - \langle e_{2} \rangle^{-}}{2\Delta\gamma_{2}} \right), \quad (2.6)$$

where superscripts ± on e indicate an ellipticity measured on an image artificially sheared by $\Delta\gamma$ in the same component and superscripts ± on ... indicate an average taken on an ensemble of source selected by quantities measured on an image artificially sheared by $\Delta\gamma$.

We note that this is identical to the methodology for DES Y1 galaxy-galaxy lensing employed in Prat et al. [37], except that we estimate the responsivity separately for the source galaxies in each radial bin around the cluster, rather than as a global scalar. The scale dependence, however, is negligible—the METACALIBRATION R is equal to within 0.5 per-cent for any two angular bins. As in other DES Y1 lensing analyses [4,37], we weight all sources in a bin uniformly—using the inverse variance of the shape measurement underlying the METACALIBRATION scheme as a weight would require a rederivation of the redshift calibration [33] and additional bookkeeping for selection bias correction, yet increases signal-to-noise ratio only mildly.

For IM3SHAPE, we use the source weights $W_j$ defined in Zuntz et al. [34] to first measure the weighted mean R of the source sample, then define the estimator for tangential shear as

$$\hat{\gamma}_t^q = \frac{\sum_{i,j} w_i^q W_j e_{ij,t}}{R \sum_{i,j} w_i^q W_j} - \frac{\sum_{i,j} W_j e_{ij,t}}{R \sum_{i,j} W_j}. \quad (2.7)$$

The IM3SHAPE e are defined with the calibration correction for additive bias already applied.

For SDSS, multiplicative bias is already corrected in the source catalog. We therefore measure tangential shear with the above equation by setting $R = 1$, and use weights $W_j = w_j \times W_{\text{bin}}$ [see Eq. (2.2) and subsequent description].

The measured shear signals are shown for DES in Fig. 3, slicing the data by density percentile, source redshift, and smoothing scale of the tracer galaxy field. SDSS signals are in Fig. 4.

## 2. Counts-in-cells

The discriminating power of density split lensing signals for cosmological parameters and parameters describing the connection of galaxies and matter is greatly improved by adding some degree of information of galaxy clustering or bias. Here, we use a very basic statistic, the mean tracer galaxy overdensity in our density quintiles, that was extensively tested in Friedrich et al. [28]—other signals could significantly improve the constraining power in the future.

Operationally, we define the mean tracer galaxy overdensity in all quintiles q as follows. We convert the raw tracer galaxy count $N_i^{\text{raw}}$ within the aperture radius around

each line of sight i to a stochastically masked count $N_i$ with fixed masking fraction by a Bernoulli draw (Sec. II A and Appendix A). We then order lines of sight by $N_i$ and take the mean of $N_i$ in each quintile q of that list as $\langle N_i \rangle^q$. The mean tracer galaxy overdensity in quintile q is

$$C^q = \frac{\langle N_i \rangle^q}{\langle N_i \rangle}, \quad (2.8)$$

where the average in the denominator runs over all lines of sight.

We note that this does account for the fact, in a stochastic fashion, that a given line of sight can end up in different density quintiles depending on the realization on masking that decides the galaxy count.

Figure 5 shows the full $P(N)$ distribution in both DES and SDSS, alongside a model evaluated at the maximum likelihood parameter values fit to the shear signal and mean tracer galaxy overdensity in quintiles. The model not only fits these mean overdensities, but also the full $P(N)$ extremely well: absolute differences in probabilities of finding N galaxies in a random line of sight, $|P^{\text{model}} - P^{\text{data}}|$, are below $10^{-3}$ and $3 \times 10^{-4}$ for any N in DES and SDSS, respectively. The bias model used for the plot is one with two-parametric stochasticity (called $b, \alpha_0, \alpha_1$ in Sec. III C), although even a simpler model can reproduce the $P(N)$ well.

## III. MODEL

In order to describe our signal as a function of

(i) cosmological parameters,

(ii) parameters that connect galaxy counts to the matter (over)density, and

(iii) nuisance parameters,

we use the model developed and tested in Friedrich et al. [28]. We only briefly summarize it here, with an emphasis on required extensions for the use on observational data, and refer the reader to that paper for details.

Let $\hat{\mathbf{n}}$ be a unit vector on the sky. The signal we have to predict in this work is the shear profile around lines of sight that fall into a certain quintile of foreground tracer density. Also, our data vector includes the average tracer density contrast in each of those density quintiles. To model these two parts of our data vector, we have to consider the following fields on the sky:

(i) $\delta_{m,2D}(\hat{\mathbf{n}})$: the line-of-sight density contrast underlying our tracer galaxies. Given the redshift distribution $n_l(z)$ of our tracer sample, this is given by

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

where $\chi$ is comoving distance and the projection kernel $q_l(\chi)$ is given in terms of $n_l(z)$ as

$$q_l(\chi) = n_l(z[\chi]) \frac{dz[\chi]}{d\chi}. \quad (3.2)$$
<figure><img src="image_5.png" /><figcaption>FIG. 3. Shear signal around overdense and underdense lines of sight in DES Y1, split by line-of-sight density quintile (top row), source redshift $z_s$ (central row) from BPZ run on METACALIBRATION (left) and MOF (right) photometry, and aperture radius (bottom panel). Our fiducial data vector is the shear around the most underdense and most overdense quintile of $\theta_T = 20'$ lines of sight as seen by sources in $z_s = 0.63-0.90$. Left-hand panels: measurements with METACALIBRATION shears. Right-hand panels: measurements with IM3SHAPE. Error bars are from 100 jackknife resamplings of the survey, consistent with our model covariance (Sec. IV). Dotted lines indicate model prediction at maximum likelihood parameters (Sec. VI) and are a good fit to the data ($\chi^2 = 171$ (METACALIBRATION) and 201 (IM3SHAPE) for $\approx 200$ d.o.f., as determined only after unblinding).</figcaption></figure>

(ii) $\delta_{m,T}(\hat{\mathbf{n}})$: the result of smoothing the field $\delta_{m,2D}(\hat{\mathbf{n}})$ with a circular top-hat aperture $T$.

(iii) $N_T(\hat{\mathbf{n}})$: the number of tracer galaxies in the aperture $T$ around the line-of-sight $\hat{\mathbf{n}}$

(iv) $\kappa_{<\theta}(\hat{\mathbf{n}})$: the convergence inside an angular radius $\theta$ around the line-of-sight $\hat{\mathbf{n}}$.

Because the Universe is isotropic, we will omit the dependence on $\hat{\mathbf{n}}$, i.e., only consider a single line of sight.
<figure><img src="image_6.png" /><figcaption>FIG. 4. Shear signal around overdense and underdense lines of sight in SDSS, shown for all quintiles, with $\theta_T = 20'$ and a single source bin of $z_s = 0.45-1$. Error bars are from 200 jackknife resamplings of the survey. Dotted lines indicate model prediction at maximum likelihood parameters (Sec. VI) and are a good fit to the data ($\chi^2 = 81$ for $\approx 70$ d.o.f., as determined only after unblinding). For comparison with Fig. 3, note the changed $\gamma_t$ axis scale.</figcaption></figure>

As detailed in Friedrich et al. [28], the density split lensing signal can be calculated from the convergence profile around lines of sight with a fixed value of $N_T$. This profile can be computed as

$$ \langle \kappa_{<\theta} | N_T \rangle = \int d\delta_{m,T} \langle \kappa_{<\theta} | \delta_{m,T} \rangle p(\delta_{m,T} | N_T) \quad (3.3) $$

where Bayes' theorem can be used to express the PDF of $\delta_{m,T}$ at fixed $N_T$ as

<figure><img src="image_7.png" /><figcaption>FIG. 5. Counts-in-cells distribution for REDMAGIC $z_T = 0.2-0.45$ galaxies in circular $\theta_T = 20'$ top-hat apertures for DES with 20 per cent masking fraction (blue) and SDSS with 10 per cent masking fraction (red). Lines indicate prediction of $P(N)$ for the maximum likelihood model fit to the lensing and counts-in-cells data in the $(b, \alpha_0, \alpha_1)$ bias model. The color of the line denotes quintiles, in the color scheme of Fig. 4, i.e., the integral under each colored segment is 0.2. The mean count in each quintile is part of the data vector and consistent with our best-fit model, as checked via its $\chi^2$ after unblinding.</figcaption></figure>

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

Here $p(\delta_{m,T})$ is the overall PDF of $\delta_{m,T}$, $P(N_T|\delta_{m,T})$ is the probability of finding $N_T$ in a line-of-sight with fixed $\delta_{m,T}$ and

$$ P(N_T) = \int d\delta_{m,T} P(N_T | \delta_{m,T}) p(\delta_{m,T}). \quad (3.5) $$

From the convergence profile $\langle \kappa_{<\theta} | N_T \rangle$ the corresponding shear profile can be computed as (cf. Friedrich et al.)

$$ \langle \gamma_t(\theta) | N_T \rangle = \frac{\cos \theta - 1}{\sin \theta} \frac{d}{d\theta} \langle \kappa_{<\theta} | N_T \rangle \quad (3.6) $$

and the shear profile around a certain quintile of tracer density is given by the average of $\langle \gamma_t(\theta) | N_T \rangle$ over the values $N_T$ occurring in that quintile (cf. Friedrich et al. for further details). Fundamentally, we therefore have to model

(i) $p(\delta_{m,T})$, the PDF of matter density smoothed inside our aperture,

(ii) $\langle \kappa(<\theta) | \delta_{m,T} \rangle$, the expectation value of convergence inside an angular radius $\theta$ around a line-of-sight with given density contrast $\delta_{m,T}$ inside our aperture, and

(iii) $P(N_T|\delta_{m,T})$, the probability of finding $N_T$ galaxies in an aperture, given its density contrast is $\delta_{m,T}$.

We describe our approaches on each of these ingredients in the following subsections, and close with a description of how we account for biases in source redshift and shear estimates, and overlap of the source redshift distribution with the tracer redshift distribution.

In all these steps, in order to predict the nonlinear 3D matter power spectrum, we use the Takahashi et al. [46] halofit approximation with the Eisenstein and Hu [47] transfer function with baryonic features, which is sufficiently accurate given our large scale binning.

## A. PDF of matter density contrast

The PDF $p(\delta_{m,T})$ can be computed from its cumulant generating function (CGF). This function can be derived at tree-level in perturbation theory with the help of the cylindrical collapse model ([28], see also pioneering work on the computation of the CGF in [14,15,48]).

The computations are numerically involved and, at least in our implementation, too slow for application in a likelihood analysis. We however show in Friedrich et al. that, on the scales used in this work, the perturbation theory computation of $p(\delta_{m,T})$ is well approximated by a log-normal distribution that matches the second and third moments $\langle \delta_{m,T}^2 \rangle$ and $\langle \delta_{m,T}^3 \rangle$ of the perturbation theory approach. We use this log-normal model [Sec. 4.1.1 in [28]] for the smoothed, projected matter density field in this work.
photon and non-Poissonian electron shot noise in CCD detectors, described by a gain factor $\alpha$. It could similarly accommodate non-Poissonianity in halo counts [55]. We allow for $\alpha$ to be different in higher and lower density regions, e.g., because more massive haloes might be more common in the former, by means of a linear dependence of $\alpha$ on $\delta_{m,T}$ as

$$ \alpha(\delta_{m,T}) = \alpha_0 + \delta_{m,T}\alpha_1. \qquad (3.17) $$

We note that a bias model without stochasticity is a common assumption made for the galaxy distribution on large scales [e.g., [1]]. MacCrann et al. [56] show that in the *Buzzard* simulations, large scale stochasticity is present. From the combination of probes with different sensitivity to b and r, such as the three galaxy and convergence auto- and cross-correlation functions, the two parameters could be disentangled. Density split statistics, in addition, are sensitive to differences in higher moments of the galaxy and matter density field, and can test and, potentially, constrain, more complex models such as $b, \alpha_0, \alpha_1$.

## D. Nuisance effects on data

In all runs on data, biases $\Delta z$ in the means of redshift distributions in DES and SDSS are accounted for at the level of the model: we marginalize over lens redshift and (multiple, in the case of a tomographic analysis) source redshift bias parameters $\Delta z$ by shifting the tracer galaxy and source galaxy redshift distributions accordingly before computing predictions for the signals. Likewise, we scale the predicted shear signal by $(1+m)$ to account for multiplicative shear biases $m$.

A more complex issue arises from the clustering of sources with the overdense and anticorrelation of sources with the underdense lines of sight. This is a common problem in cluster lensing or galaxy-galaxy lensing, accounted for by so-called *boost factors* [57,58].

In the case of density split lensing, we apply the assumption of linear bias to predict the radius dependence of boost factors and their effect, given the non-thin lenses, on our model predictions. For a given tracer redshift distribution and a the matter field at redshift $z_s$, the angular clustering $w_q(\theta, z_s)$ of quintile $q$ with matter can be calculated with the same formalism as the convergence in Sec. III B. Assuming a linear bias of source galaxies $b_s$, their redshift distribution at separation $\theta$ from quintile $q$ changes due to clustering to

$$ n_s(z) \rightarrow n_{s,q}(z, \theta) = [1 + b_s w_q(\theta, z)] n_s(z) \times \left[ \int dz_s [1 + b_s w_q(\theta, z_s)] n_s(z_s) \right]^{-1}. \qquad (3.18) $$

The lowest redshift bin in DES Y1 or the *Buzzard* simulations and the sources in SDSS have sufficiently strong overlap with the lens redshift distribution that we include this effect in the modeling and marginalize over $b_s$ in the analysis. This means that we use a different source redshift distribution for predicting each point of the density split, radially binned shear signal data vector. While $b_s$ is in reality a function of $z_s$, one can very accurately describe the deboosting of the lensing signals by an effective $b_s$ because the radial profile shape of the shear signal is almost independent of source redshift.

We note that in this derivation we neglect a second, but likely subdominant effect: the source redshift dependence of the probability of failing to include a source in the DES shape catalogs due to blending, that might cause a similar density dependence of source $n_s(z)$.

A potential spurious signal is due to intrinsic alignment of physical source galaxy shapes with the underdense or overdense lines of sight due to gravitational interactions (see [59,60] for a review). For cross-correlations between the positions of object and gravitational shear, such as counts and lensing in cells, intrinsic alignments affect only the signal from source galaxies that are physically associated with the lensing objects, i.e., if redshift distributions of source galaxies and lensing objects overlap. This is the case primarily in the lowest redshift bin for DES Y1. Hence test (4) in Sec. V D, which demonstrates the robustness of the results to removing the lowest source redshift bin from the data vector, indicates that the current analysis is at most weakly affected by intrinsic alignments. This is in agreement with our expectation that the tidal alignment of galaxies with the comparatively small mean over- and underdensities of our density quintiles is small at $\theta > 20'$ separation.

On small scales, baryonic effects can modify the matter power spectrum from its dark matter only prediction, primarily by affecting overdense regions. For our statistic, this could be absorbed by the bias model on scales smaller than the top-hat aperture $\theta_T$. The shear signal is used on scales larger than $\theta_T = 20'$ only, and parameter constraints are robust to a more conservative scale cut of $\theta > 40'$ (Sec. V D). We hence do not expect a significant impact of baryonic effects on the parameter constraints at the accuracy level of the current analysis, but note that these effects require further study for future, more constraining analyses.

## IV. COVARIANCE

In order to interpret our measurements, we need an accurate description of their covariance. We construct this covariance from a large number of mock realizations of our data vectors. In that, we make use of the fact that the noise in our measurements can be separated into two components: a contribution from shape noise and a contribution from large scale structure and shot noise in the galaxy
catalog. This approach is similar to the one of Murata
et al. [61].

In the following, we describe how we measure these
contributions, and how we combine them into a covariance
matrix.

We assume in all following analyses that the signals
measured in SDSS and DES Y1 are uncorrelated, justified
by the fact that the survey footprints (using only the
contiguous SPT region of DES) are well separated.

A. Shape noise

The primary contribution to the shape we measure for
any individual galaxy in our survey is the sum of its
intrinsic shape and measurement noise, not the weak
gravitational shear that distorts the galaxy image.

Because of this dominance of the noise over the signal,
and because the intrinsic shapes of neighboring galaxies are
almost uncorrelated, we can measure shape noise by
rotating each galaxy in our shape catalog by an independent
random angle. The shear signal around our actual under-
dense and overdense lines of sight as measured from these
rotated source catalogs represents a random realization of
the shape noise (cf., e.g., [62], for a similar technique for
shear peak statistics, [63] for void lensing, and [61] for
cluster lensing).

In measuring the signals on the rotated catalogs, we take
care to use the same methodology as for the measurements
on data. That means we use each randomly rotated source
catalog for cross-correlation with all of the maps contrib-
uting to our data vector. As on the data, we subtract the
mean shears measured around random points, for which we
simply use the centers of all HEALPIX pixels that are used as
lines of sight in any density quintile. The subtraction of
shear around random points considerably reduces shape
noise on large scales (see Fig. 6, and refer to [64] for a
detailed study of the effect).

<figure><img src="image_8.png" /><figcaption>FIG. 6. Variance of the shear signal around troughs due to shape noise in DES Y1 METACALIBRATION $z_s = 0.63–0.9$ galaxies. Subtracting the shape noise around random points (cross symbols), as we also do in our data, lowers the variance considerably on scales above 300 arcmin.</figcaption></figure>

B. Cosmic variance and shot noise

Two additional effects cause our signal to deviate from
its expectation value:

(i) The cosmic density field present in our survey
volume is a random realization. This is true both
for the volume in which our tracer galaxies are
located (and in which the signal of troughs and
overdense lines of sight originates) and for the
redshift range along the line of sight in between
us and the source galaxies that is not contained in
that volume. This causes there to be noise in the true
convergence around the lines of sight we identify,
and in the counts-in-cells distribution.

(ii) In a given realization of the matter density field, tracer galaxies could be placed differently (for instance, according to Poisson noise around their expectation value in any given volume). Which of these possible galaxy catalogs is realized causes there to be a different true shear signal around what we identify as troughs and overdense lines of sight, and a different counts-in-cells distribution.

On the scales we care about in this work, we can measure
the sum of both contributions to the covariance, to good
approximation, from log-normal simulations of the related
matter and convergence fields and Poissonian realizations
of the tracer galaxy catalog. We do this by generating a
large number of realizations of these fields and catalogs
with FLASK [49].

We note that this part of the covariance is dependent on cosmology and the parameters describing the connection of galaxies and matter. For the covariance in this work, we will assume the settings of the *Buzzard* simulations, namely a fiducial flat $\Lambda$CDM cosmology with $\Omega_{m,0}=0.286$, $\sigma_8 = 0.82$, $\Omega_b = 0.047$, $n_s = 0.96$ and $H_0 = h \times 100 \text{ km s}^{-1}$ with $h = 0.7$.

For the matter and associated galaxy field in the tracer
redshift range, we use the power spectrum with a linear bias
of $b = 1.54$, a redshift distribution, and a mean density of
the tracer galaxy population as in the *Buzzard-v1.1* suite of
simulations. We assume Poissonianity of the galaxy count
at fixed density, i.e., the *b* model (Sec. III C). Note that the
relation of galaxies and matter in the *Buzzard* simulations
[56] and, potentially, the actual Universe is more complex
than that. We ensure, using mock likelihood runs on
*Buzzard*, that this does not mean our covariance from
the log-normal mocks is significantly underestimated (see
Appendix E). To set the log-normal parameter of the
projected matter field (i.e., the minimum allowed value
of $\delta_m \ge -\delta_0$), we use the methodology of Friedrich *et al.*
[28]. In the *Buzzard* cosmology and at a top-hat smoothing
scale of 20 arcmin, this yields $\delta_0 = 0.669$. Details are given
in Appendix B.

For the source redshift distributions of the simulated
convergence fields, we use those estimated for the source
samples in our data.
We separate the convergence field into two parts: a component correlated with the matter field that our tracer galaxies populate, and an uncorrelated component (mostly comprised of the parts of the lensing kernel in front and behind our tracer galaxies). The correlated component is modeled as a log-normal field with cross-power spectrum and $\kappa_0$ set to match the perturbation theory predictions for $\langle\delta\kappa\rangle$ and $\langle\delta_T^2\kappa\rangle$ at a fiducial smoothing scale. This constrains the auto power spectrum of this component to be only a fraction of the total convergence power spectrum. The uncorrelated component is then simulated as a Gaussian random field that is uncorrelated to all other fields and whose power spectrum is chosen to give the correct total convergence power spectrum (see Appendix B for the details of the procedure).

We apply the same mask to the tracer galaxies as in our data (or in our simulations, for the mock analysis described in Appendix E), and the same prescription for splitting the survey into lines of sight of different density.

We then measure tangential shear signals, as in our data, yet on the noiseless shear maps output by FLASK at $N_{side} = 4096$ resolution. In order to sample the density fields as in our data, we use the sum of weights of sources in our actual shear catalogs situated in a pixel as the weight of the shear signal in that pixel. We do this both for the correlated and the uncorrelated part of the convergence field (see above) and coadd the two signals. In addition, as in our data, we measure the mean tracer galaxy overdensity in our density quintiles.

On scales much smaller than the aperture radius $\theta_T = 20'$, a checkerboard pattern in the off-diagonal shape noise covariance is apparent (see Fig. 7). We find that this is due to an interference of the HEALPIX grid we use to sample the density field and the angular binning scheme for our shear signal—for adjacent HEALPIX pixels, sources move from one angular bin to the next and their intrinsic shape orientation changes relative to the pixel centers. Since these effects are present in the data as well (as seen from the jackknife covariance) and only significant on angular scales below our scale cut, we do not attempt to address them further.

### C. Constructing the covariance matrix

We create 1000 realizations of both the shape noise and the large scale structure and shot noise contributions to the covariance. Despite this relatively large number, there is noise in our estimated covariance matrix. When inverting the covariance matrix to calculate $\chi^2$ values and run a likelihood analysis, this noise has two consequences.

First, the inverse of a noisy estimate of the covariance matrix is a biased estimate of the inverse covariance matrix. We follow the correction described in Hartlap *et al.* [65] to correct for this effect, i.e., we multiply the $\chi^2$ calculated

<figure><img src="image_9.png" /><figcaption>FIG. 7. Covariance of shear signal around SDSS troughs from shape noise (top), large scale structure and shot noise (second from top, note the different scale). The model (third from top) is the sum of these two and closely matches the jackknife covariance (bottom). White lines indicate aperture radius $\theta_T = 20'$—only data above that is used in the likelihood.</figcaption></figure>
with the inverse of our estimated covariance matrix by a factor

$$f^{\text{AH}} = \frac{N^{\text{cov}} - N^{\text{data}} - 2}{N^{\text{cov}} - 1}. \qquad (4.1)$$

The number of entries in our data vector is at most $N^{\text{data}} = 208$ in the fiducial DES analysis and we use $N^{\text{cov}} = 960$ realizations of the log-normal field to estimate the covariance, which means $f^{\text{AH}}$ is 0.78 or larger for all our likelihood runs. We confirm, using independent log-normal mock realizations of our data vector, that the inverse covariance matrix rescaled such does lead to a consistent $\chi^2$ distribution of residuals (Appendix D).

Second, the noise in the inverse variance leads to additional scatter in the best fit we find [66–68]. Under the assumption that the model is linear in all parameters within the range probed, this can be compensated by multiplying $\chi^2$ with a factor

$$f^{\text{DS}} = \left[ 1 + \frac{(N^{\text{data}} - N^{\text{par}})(N^{\text{cov}} - N^{\text{data}} - 2)}{(N^{\text{cov}} - N^{\text{data}} - 1)(N^{\text{cov}} - N^{\text{data}} - 4)} \right]^{-1}, \qquad (4.2)$$

where $N^{\text{par}}$ is the number of free parameters in the model.

These corrections are only appropriate for a monolithic covariance estimated from a fixed number of independent realizations. From the previous subsections, we can get independent, unbiased estimates of the two contributions, Cov<sup>shape</sup> and Cov<sup>LSS</sup>. The sum of Cov<sup>shape</sup> + Cov<sup>LSS</sup> would be an unbiased and less noisy estimate of the total covariance. But to apply the above corrections, we need to resort to coadding shape noise and cosmic variance realizations before estimating the full covariance matrix.

In addition to the 960 realizations used to estimate the covariance, we use 40 independent realizations to confirm that our prediction indeed matches the mean signal measured from the log-normal simulations at the expected $\chi^2 \approx N_{\text{data}}$. This is a test of both the numerical scheme employed by FLASK and the implementation of the analytical calculations of Friedrich et al. [28]. We find that the two are in excellent agreement, except for a small offset of the predicted and measured counts-in-cells statistic. The mean tracer galaxy overdensities (Sec. II D 2) we measure in log-normal mocks are offset from the predictions at most at the $10^{-3}$ level. We hypothesize that this is due to resolution effects of the simulations, but cannot exclude that similar effects could also present in the data.<fnref n="2" /> To compensate for this, we boost the variance of each of the four counts-in-cells entries in our data vector by 0.002².

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

<figure><img src="image_11.png" /><figcaption>FIG. 8. Covariance of tomographic shear signal around DES density quintiles (20', lowest to highest density quintile, then lowest to highest redshift source bin, division indicated by white lines) and counts in cells (last 5 bins) from shape noise (top), large scale structure and shot noise (second from top, note the different scale), and the full model from the sum of these two (third from top). Bottom panel: diagonal of shear around troughs in intermediate source redshift bin. Counts in cells residuals have been rescaled by factor 1/50 to match a common color scale.</figcaption></figure>

<footnote n="2">²We confirm, however, that the mean tracer galaxy overdensities in our data are well fit by the model at its maximum likelihood parameters.</footnote>
Using this covariance and $f^{\text{AH}}$ (but not $f^{\text{DS}}$) as defined above, the mean of all realizations with no shape noise matches the predicted signal at the true input parameters at total $\chi^2 = 0.19$ with 208 degrees of freedom, proving the numerical accuracy of the prediction code at a sufficient level. Additional tests of our likelihood pipeline run on the 40 log-normal realizations are shown in Appendix D.

### D. Comparison with jackknife covariances

From jackknife resamplings of our data, we can internally estimate the covariance matrix. While more care would have to be taken for applying this estimate of the covariance matrix in a likelihood analysis [69], it does provide confirmation of our scheme to compare the jackknife estimate to the covariance estimated above.

Figure 7 shows the shape noise and cosmic variance + shot noise components of the covariance matrix and compares their sum to the jackknife covariance, for the shear signal around underdense lines of sight in SDSS. The same for the full density and source-redshift tomographic covariance matrix of DES, including counts-in-cells, is displayed in Fig. 8.

### V. LIKELIHOOD

We compare our data to model predictions $M$ in a Bayesian fashion, i.e., we sample the posterior distribution of model parameters $p$ with a Monte Carlo Markov Chain (MCMC) run on the likelihood

$$-2 \ln \mathcal{L} = f^{\text{AH}} f^{\text{DS}}[D - M(p)]^{\text{T}} C^{-1} [D - M(p)] + \text{Prior}(p). \quad (5.1)$$

For our fiducial likelihood analysis, we remove the following parts of the full data vector:

(i) lensing and counts-in-cells signal for any aperture radii other than $\theta_T = 20'$—on smaller smoothing scales, small but significant deviations of our model and measurements in $N$-body simulations appear [28]. Smoothing on larger scales than 20' yields signals with errors that are highly correlated to the 20' measurements, thus adding little independent information.

(ii) lensing signal on scales smaller than $\theta_T = 20'$—small but significant deviations of our model and measurements in $N$-body simulations are present on scales smaller than the aperture radius $\theta_T = 20'$. The lensing signal on these scales has low signal-to-noise ratio. In addition, shape noise in adjacent small-scale bins is anticorrelated, visible as the checkerboard pattern in the lower left of Fig. 7. This is due to interference of the radial binning scheme with the HEALPIX grid of lines of sight: when we measure the contribution of a source galaxy to the shear signal around two adjacent lines of sight, its intrinsic orientation relative to a line of sight and its distance

from the line of sight change coherently. While the effect is consistently seen in jackknife and model covariance, it makes these small scale lensing signals numerically redundant. This leaves 17 angular bins in each shear profile.

(iii) signal for quintile 3—the signals we use are not linearly independent between all quintiles; we therefore discard the signal in the median quintile, which is close to zero by construction anyway.

Therefore, in all of the following, unless otherwise noted, $D$ contains the shear signals measured at $\theta = 20-600'$ and the relative overdensity of tracer galaxy count for the lower two and upper two quintiles of galaxy count, measured in $\theta_T = 20'$ apertures. For the source tomographic DES Y1 analysis, these are 208 entries (72 for SDSS).

The precision matrix $C^{-1}$ is estimated as detailed in Sec. IV.

In the following subsections, we describe our choice of parametrization, the nuisance parameters and associated priors, and the consistency tests we perform before unblinding the estimated cosmological parameters.

### A. Cosmological parameters

Since this is our first cosmological analysis of counts and lensing in cells, we choose to only vary a minimal set of cosmological parameters, adopting fixed priors for ones that the density split lensing and counts signal is not very sensitive to. For the fiducial run of our likelihood, we validate this approach by marginalizing over these parameters with informative external priors.

All likelihoods assume a flat $\Lambda$CDM cosmology. The main parameters we wish to constrain are the matter density in units of the critical density $\Omega_m$, and the amplitude of structure in the present day Universe, parametrized as the RMS of overdensity fluctuations on $8h^{-1}$ Mpc scale, $\sigma_8$.

In an alternate run of our likelihoods, we will also leave free the parameter $S_3$ that describes the skewness of the matter density field when smoothed over the given aperture and redshift range,

$$S_3 = \langle \delta^3 \rangle / \langle \delta^2 \rangle^2. \quad (5.2)$$

$S_3$ was first defined by Peebles [70], who derived a perturbation theory prediction $S_3 = 34/7$ for the unsmoothed matter density field, and later generalized to top-hat smoothed fields and higher orders [48,71,72]. Perturbation theory predicts also the smoothed $S_3$ to be almost independent of $\Omega_m$ and $\sigma_8$ and to only vary slowly with redshift or scale. A skewness $S_3$ that is inconsistent with these predictions could be caused either by non-Gaussian initial density fluctuations (although CMB limits set tight constraints on these [73]) or by physics beyond gravity that affect collapse either in the overdense or underdense regime.

We assume wide, flat priors for $\Omega_m$, $\sigma_8$ and, in the likelihood runs that vary it, $S_3$ that do not limit the range
<table><caption>TABLE I. Priors for likelihood runs.</caption><thead><tr><th>Parameter</th><th>Prior</th></tr></thead><tbody><tr><td colspan="2"><em>Cosmology</em></td></tr><tr><td>Ω<sub>m</sub></td><td>flat (0.1, 0.9)</td></tr><tr><td>σ<sub>8</sub></td><td>flat (0.2, 1.6)</td></tr><tr><td>S<sub>3</sub></td><td>fixed to PT/flat</td></tr><tr><td>Ω<sub>b</sub></td><td>fixed (0.047)</td></tr><tr><td>h</td><td>fixed (0.70)</td></tr><tr><td>Ω<sub>r</sub>h<sup>2</sup></td><td>fixed (4.15 × 10<sup>-5</sup>)</td></tr><tr><td colspan="2"><em>Tracer galaxies</em></td></tr><tr><td>b</td><td>flat (0.8, 2.5)</td></tr><tr><td>r</td><td>flat (0, 1)</td></tr><tr><td>α<sub>0</sub></td><td>flat (0.1, 3.0)</td></tr><tr><td>α<sub>1</sub></td><td>flat (-1.0, 4.0)</td></tr><tr><td colspan="2"><em>Tracer galaxy photo-z shift</em></td></tr><tr><td>Δz<sub>1,DES</sub></td><td>Gauss (0.003, 0.008)</td></tr><tr><td>Δz<sub>1,SDSS</sub></td><td>Gauss (0.002, 0.006)</td></tr><tr><td colspan="2"><em>Source photo-z shift</em></td></tr><tr><td>Δz<sub>s,METACALIBRATION</sub><sup>2</sup></td><td>Gauss (-0.019, 0.018)</td></tr><tr><td>Δz<sub>s,METACALIBRATION</sub><sup>3</sup></td><td>Gauss (+0.009, 0.016)</td></tr><tr><td>Δz<sub>s,METACALIBRATION</sub><sup>4</sup></td><td>Gauss (-0.018, 0.031)</td></tr><tr><td>Δz<sub>s,IM3SHAPE</sub><sup>2</sup></td><td>Gauss (-0.024, 0.018)</td></tr><tr><td>Δz<sub>s,IM3SHAPE</sub><sup>3</sup></td><td>Gauss (-0.003, 0.016)</td></tr><tr><td>Δz<sub>s,IM3SHAPE</sub><sup>4</sup></td><td>Gauss (-0.057, 0.031)</td></tr><tr><td>Δz<sub>s,SDSS</sub><sup>2</sup></td><td>Gauss (-0.014, 0.011)</td></tr><tr><td colspan="2"><em>Shear calibration</em></td></tr><tr><td>m<sup>i</sup><sub>METACALIBRATION</sub>(i = 2, 3, 4)</td><td>Gauss (0.012, 0.023)</td></tr><tr><td>m<sup>i</sup><sub>IM3SHAPE</sub>(i = 2, 3, 4)</td><td>Gauss (0.0, 0.029)</td></tr><tr><td>m<sub>SDSS</sub></td><td>Gauss(0.0, 0.05)</td></tr></tbody></table>

sampled by the likelihoods. We fix the Baryon density $\Omega_b = 0.047$, the spectral index of primordial density fluctuations $n_s = 0.96$, and a dimensionless Hubble parameter $h = 0.7$, equal to the values used in the *Buzzard* simulations and consistent with best constraints. For the transfer function of primordial to initial matter power spectrum, we assume a radiation density $\Omega_r h^2 = 4.15 \times 10^{-5}$. The evolution of expansion and growth of structure in the late universe assumes only matter and a cosmological constant.

An overview of these choices is given in Table I.

## B. Nuisance parameters

In our likelihoods, we apply three different models to describe the distribution of REDMAGIC galaxy count $N_T$ inside an aperture at given mean matter overdensity $\delta_{m,T}$, $P(N_T|\delta_{m,T})$. Details of this are described in Sec. III C, and sampling ranges for the parameters $b$, $(b, r)$, or $(b, \alpha_0, \alpha_1)$, designed to span any physically sensible configurations, are listed in Table I.

Similar to previous cosmological constraints derived from DES Y1 data, we assume and always marginalize over nuisance parameters describing photometric redshift

and shear biases in our measurements. The nuisance parameter for the redshift bias of REDMAGIC sources in the $z_T = 0.2-0.45$ redshift range that is constrained from cross-correlations with a sample of galaxies with spectroscopic redshifts as in Cawthon et al. [74]. Specifics of this are described in Appendix F 1.

The photometric redshift biases and multiplicative shear biases of source galaxies are described by two parameters in each redshift bin. The three bins we use, i.e., all but the lowest redshift bin of DES Collaboration et al. [1], are labeled as $i = 2, 3, 4$ in Table I. Priors on the redshift biases are taken from the combination of the redshift distributions of a matched sample of galaxies in the COSMOS survey and angular cross-correlation with REDMAGIC galaxies [38,39,75] as described in detail in Hoyle et al. [33]. The priors on multiplicative shear bias in DES Y1 are described in detail in Zuntz et al. [34]. Both of these priors are widened in our analysis to account for their potential correlation between bins [see Appendices of [33,34]], conservatively assuming comparable signal-to-noise ratio in each bin.

Multiplicative bias in an independent SDSS shear catalog that is consistent with the one we use [76] was investigated in detail in Mandelbaum et al. [3]. The authors in that paper find a Gaussian uncertainty related to multiplicative shear calibration of $\sigma_m = 0.037$, in addition to photometric redshift biases over which we marginalize separately. We assume a slightly more conservative Gaussian uncertainty of $\sigma_m = 0.05$ for the multiplicative shear bias in the SDSS catalog used in this work.

These priors are also summarized in Table I.

## C. Sampling and evidence

To sample the posterior likelihoods efficiently, we employ both the EMCEE [77] and the MULTIEST [78] algorithm. The latter has the advantage of also estimating Bayesian evidences $E$,

$$E \propto p(D|\text{model}) = \int d\mu p(D|\text{model}, \mu) p(\mu|\text{model}), \qquad (5.3)$$

where $\mu$ are the parameters of the model.

Knowing the evidence of two models 1 and 2 allows comparing them with the Bayes factor, $E_1/E_2$. If the latter ratio exceeds 3.2 or 10, the evidence for model 1 over model 2 can be considered substantial or strong in the nomenclature of Jeffreys [79].

We confirm in the tests performed in Appendices D and E that both sampling algorithms and the analysis of their outputs with CHAINCONSUMER [80] and custom codes yield reliable parameter constraints and test results.

## D. Blinding and tests

Since most of the tests in this paper were performed after the scaling factors of the initial, blinded shear catalogs had
been revealed [34], we primarily rely on parameter level blinding. This means that we do not compare measurements on data to predictions in a known cosmology before the following tests are passed:

(1) Log-normal simulations show $\chi^2$ values of data vs. model at the input set of parameters consistent with a fiducial $\chi^2$ distribution with the appropriate number of d.o.f. Likelihood runs on mock data have a coverage within expectations (i.e., the input cosmology lies within the confidence interval the expected fraction of times). Results are unremarkable and described in Appendix D.

(2) Bayesian model comparisons run on log-normal simulations without stochasticity do not provide evidence for more complex models. We find that this requirement is met, both for the *b*, *r* model of stochasticity and models with free skewness $S_3$, in Appendix D.

(3) 21 independent *Buzzard N*-body realizations of our data vector give consistent $\chi^2$ relative to the model evaluated at the input cosmology and independently measured nuisance parameters. Their coverage in likelihood runs, i.e., the number of times the input cosmology is within derived confidence limits, is within expectations only for the *b*, $\alpha_0$, $\alpha_1$ model of bias and stochasticity (Appendix E). The fact that the most complex bias model is required may be particular to these mock galaxy catalogs, which may have different relations to matter density than real REDMAGIC galaxies. We still take this as evidence that the most general stochasticity model, unless disfavored by the data, needs to be considered in our analysis.

(4) Likelihood runs on *N*-body realizations are insensitive to replacing true source redshift distributions with source redshift distributions estimated from BPZ and marginalizing over $\Delta z$ uncertainties. Results: we find that the mean shifts in cosmological parameters are at or below the ten percent level of their statistical uncertainty, and that the statistical uncertainty increases by less than five per cent due to marginalization over $\Delta z$, both tested with the *b*, $\alpha_0$, $\alpha_1$ model for the galaxy-matter connection.

Once these tests are successful, we continue to make tests on likelihood analyses run on the data itself. To ensure that these do not introduce experimenter bias, before looking at any chains we shift all cosmological parameters by a constant unknown vector, uniformly distributed between +2 and -2 standard deviations of the parameters as found from *N*-body simulations. We then proceed with the following tests, the results of which are shown in Fig. 9:

(5) Cosmological constraints from the data are consistent between the fiducial METACALIBRATION and additional IM3SHAPE measurements. For models including galaxy stochasticity (lower two panels of Fig. 9) this is the case to a fraction of the

<figure><img src="image_12.png" /><figcaption>FIG. 9. Likelihoods run on variants of the data to test robustness to analysis choices. Use of IM3SHAPE shape catalogs (turquoise), removal of the highest (orange) or lowest (violet) source redshift bin, use of a direct estimate of source redshift distributions from COSMOS instead of BPZ (magenta), removal of small (green) or large scales (yellow) in the lensing data or neglect of the clustering of sources in the lowest redshift bin with the lenses (brown) do not have an effect on recovered parameters beyond their statistical uncertainty. Contours are centered on the fiducial result in the *b* model. The same likelihood run on fully independent SDSS data vectors (grey) yields consistent results. Panels show different models for the connection of galaxies and matter (top: *b* only, center: *b*, *r*, bottom: *b*, $\alpha_0$, $\alpha_1$).</figcaption></figure>
statistical uncertainty. For the *b* only model (top panel), IM3SHAPE constraints on $\Omega_m$ are offset by $\approx\sigma/2$. Accounting for the fact that shape noise is largely uncorrelated between the two catalogs, is possible that this is simply a statistical fluctuation. We note, however, that DES Collaboration et al. [1] found a similar discrepancy, likely attributed to the multiplicative bias or source redshift calibration of the IM3SHAPE highest source redshift bin.

(6) Cosmological constraints are robust to removing the lowest or highest source redshift bin from the data vector. This is the case for all models, indicating that the calibration of METACALIBRATION catalogs is consistent between bins.

(7) Cosmological constraints are robust to replacing the source redshift distributions estimated by BPZ by ones directly estimated from COSMOS. Again, this is the case for all models, with a noticeable but insignificant offset in the *b* only model.

(8) Cosmological constraints are robust to cutting scales smaller than $2 \times \theta_T$ or larger than 250' from the shear signal. Removal of small scale shear information shifts $\sigma_8$ in the $b, \alpha_0, \alpha_1$ model by approximately $1\sigma$. Given the cosmic variance in large-scale modes and the unremarkable result of all other variants of the scale cut test, this does not pose a significant issue.

(9) Cosmological constraints are robust to not correcting for clustering of the overlapping source redshift bins with the matter distributions around overdense and underdense lines of sight. While this is not necessary a null test—it could be possible that we need to account for the effect—we find that marginalizing over $b_s$, the source bias in the lowest redshift source bin, neither significantly widens nor shifts the contours in either model.

(10) Cosmological constraints are consistent between DES and SDSS. We note that these are completely independent data sets, i.e., have no cross-covariance, and thus we *a priori* expect larger offsets between the two than in the other tests. We find that constraints on $\Omega_m$ are very similar and $\sigma_8$ is offset by $\approx 1\sigma$, both consistent with these expectations.

In addition, we confirm that the central value of the nuisance parameter priors (multiplicative shear bias, tracer, and source galaxy redshift biases as defined in Hoyle et al. [33], Zuntz et al. [34] Appendices F1 and F2) is within the $1\sigma$ confidence interval of the posteriors for both DES and SDSS.

Only after unblinding do we test whether the model at its maximum likelihood parameters is a good fit to the data. For the tomographic data vector of DES Y1, there are 208 elements fit with 13 parameters in the ($b, \alpha_0, \alpha_1$) model. Because the model is not linear in the parameters, the number of degrees of freedom and expectation value for the $\chi^2$ distribution is not known precisely [81], but likely between 208 – 13 and 208. Its standard deviation is

$\sigma_{\chi^2} = \sqrt{2N_{\text{d.o.f.}}} \approx 20$. The data vectors for the two shear pipelines have a $\chi^2_{\text{METACALIBRATION}} = 171$ and $\chi^2_{\text{IM3SHAPE}} = 201$, respectively, both consistent with expectations for multivariate Gaussian noise around a signal that is correctly described by our model. The *b* only and ($b, r$) models give equally acceptable fits. For the SDSS single source bin data vector with 72 entries and 9 parameters in the ($b, \alpha_0, \alpha_1$) model, we find $\chi^2 = 81$, and equally acceptable results for the other models.

We also perform a run of the fiducial DES Y1 data vector with a full cosmological model that marginalizes, in addition, over baryon density $\Omega_b$, spectral index of primordial density fluctuations $n_s$ and Hubble parameter $h$ with the flat priors also used in DES Collaboration et al. [1]. We find that this does not shift or increase the uncertainty on the reported parameters at a discernible level in any of the models for the connection of galaxies and matter.

## VI. COSMOLOGICAL CONSTRAINTS

We perform likelihood analyses, i.e., we determine the probability of finding our fiducial data vectors (Sec. II) as a

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

FIG. 10. Constraints on matter density $\Omega_m$, amplitude of late-time structure formation $\sigma_8$, and galaxy bias of the REDMAGIC tracer galaxies. Shown are results for three different models for the connection of matter density and galaxy density on the scales of our apertures: linearly biased tracers with Poissonian noise (black, *b*); biased, stochastic tracers (green, *b*, *r*); and biased tracers with density-dependent stochasticity (orange, *b*, $\alpha_0$, $\alpha_1$). The Bayes factors of these different models indicate substantial evidence that stochasticity of galaxy count is required to model the measurements. Constraints from DES galaxy and shear auto- and cross-two-point correlation functions for fixed neutrino mass are shown as blue, dotted contours [1].
<table><caption>TABLE II. Constraints from counts and lensing in cells likelihood runs. Bayes factors are quoted relative to the *b*, *r* model. Results for the *b* only likelihoods are not shown, since this simpler model is disfavored by the data. Constraints on cosmological parameters and REDMAGIC bias in *z* ≈ 0.2–0.45 from two-point functions are reproduced from DES Collaboration *et al.* [1] for comparison.</caption><thead><tr><th>Data</th><th>Model</th><th>Bayes factor</th><th>S<sub>8</sub></th><th>Ω<sub>m</sub></th><th>σ<sub>8</sub></th><th>b</th><th>r</th><th>α<sub>0</sub></th><th>α<sub>1</sub></th></tr></thead><tbody><tr><td>DES</td><td>b, r</td><td>≡1.0</td><td>0.90<sup>+0.11</sup><sub>-0.08</sub></td><td>0.26<sup>+0.04</sup><sub>-0.03</sub></td><td>0.97<sup>+0.07</sup><sub>-0.06</sub></td><td>1.45<sup>+0.10</sup><sub>-0.11</sub></td><td>0.77<sup>+0.10</sup><sub>-0.13</sub></td><td>...</td><td>...</td></tr><tr><td>SDSS</td><td>b, r</td><td>≡1.0</td><td>0.78<sup>+0.13</sup><sub>-0.08</sub></td><td>0.25<sup>+0.05</sup><sub>-0.04</sub></td><td>0.86<sup>+0.06</sup><sub>-0.05</sub></td><td>1.48<sup>+0.09</sup><sub>-0.09</sub></td><td>0.70<sup>+0.16</sup><sub>-0.14</sub></td><td>...</td><td>...</td></tr><tr><td>DES</td><td>b, α<sub>0</sub>, α<sub>1</sub></td><td>0.7</td><td>0.78<sup>+0.05</sup><sub>-0.04</sub></td><td>0.28<sup>+0.05</sup><sub>-0.04</sub></td><td>0.80<sup>+0.06</sup><sub>-0.07</sub></td><td>1.75<sup>+0.22</sup><sub>-0.26</sub></td><td>...</td><td>1.5<sup>+0.4</sup><sub>-0.6</sub></td><td>1.7<sup>+1.1</sup><sub>-0.9</sub></td></tr><tr><td>SDSS</td><td>b, α<sub>0</sub>, α<sub>1</sub></td><td>1.6</td><td>0.76<sup>+0.08</sup><sub>-0.07</sub></td><td>0.28<sup>+0.07</sup><sub>-0.05</sub></td><td>0.80<sup>+0.08</sup><sub>-0.11</sub></td><td>1.18<sup>+0.37</sup><sub>-0.23</sub></td><td>...</td><td>2.3<sup>+0.3</sup><sub>-0.5</sub></td><td>2.9<sup>+1.1</sup><sub>-1.0</sub></td></tr><tr><td>3 × 2pt, fixed ν</td><td></td><td></td><td>0.80<sup>+0.02</sup><sub>-0.02</sub></td><td>0.26<sup>+0.02</sup><sub>-0.03</sub></td><td>0.85<sup>+0.06</sup><sub>-0.05</sub></td><td>1.54<sup>+0.09</sup><sub>-0.10</sub></td><td>DES Collaboration *et al.* [1]</td><td></td><td></td></tr></tbody></table>

function of the parameters of our model (Sec. III) and given their covariance (Sec. IV). We use models of different complexity for the connection of galaxies and matter density—one with linear bias only (*b*), one adding stochasticity (*b*, *r*), and one allowing for density dependence of stochasticity (*b*, *α*<sub>0</sub>, *α*<sub>1</sub>) (for details see Sec. III C). Our philosophy, decided with parameters still blinded, will be to compare these models via their Bayesian evidence (Sec. V C) and report the results for models that are supported by the data.

Figure 10 shows constraints on the matter density $\Omega_m$, the amplitude of late-time structure formation $\sigma_8$, and galaxy bias of the REDMAGIC tracer galaxies, when marginalizing over all remaining model parameters. Confidence limits are summarized in Table II. The stochastic models, favored by the data (see next subsection), constrain the matter density consistently as $\Omega_m = 0.26_{-0.03}^{+0.04}$ for the (*b*, *r*) model and $\Omega_m = 0.28_{-0.04}^{+0.05}$ for the (*b*, $\alpha_0$, $\alpha_1$) model from DES data. The degeneracy directions of $\Omega_m - \sigma_8$ for the two stochasticity models are different, thus leading to a higher central value of $\sigma_8 = 0.97_{-0.06}^{+0.07}$ for (*b*, *r*) and a lower $\sigma_8 = 0.80_{-0.07}^{+0.06}$ for (*b*, $\alpha_0$, $\alpha_1$).

Bias models and cosmology are thus interdependent: a prior, even a weak one, on the values of stochasticity parameters would significantly improve these cosmological constraints—if galaxies have less stochasticity, the relevant regions of the green and orange contours in Fig. 10 are closer to the black region. Likewise, external information on cosmological parameters allows to constrain bias parameters and, potentially, even choose between the bias models. If the true cosmology is $\Omega_m = 0.3$ and $\sigma_8 = 0.8$, both the *b* and the *b*, *r* model are ≈2σ off, while the *b*, $\alpha_0$, $\alpha_1$ model is consistent.

For a sense of how these results compare to two-point function measurements, Fig. 10 contains constraints from the three tomographic auto- and cross-correlation functions of DES REDMAGIC galaxy positions and source galaxy shapes [3x2pt, [1]].<fnref n="3" /> The 3x2pt contours are

tighter than the constraints that counts and lensing in cells yield. This is due to the wide freedom on stochasticity parameters and models that we have allowed: if we could fix the stochasticity (such as in the black contour with *r* = 1), the smaller scale density PDF measurements would yield highly competitive cosmological constraints. It is clear from this and the different degeneracies that a joint analysis would result in improved constraints—yet we are lacking a covariance matrix and inference pipeline to perform this at this point. Prima facie, the 3x2pt constraint is consistent with any of the bias models, and indicates a relatively small stochasticity, i.e., a point in parameter space close to where the black, green, and orange contours intersect.

Finally, we compare the results from our DES Y1 and SDSS analysis. Within their mutual uncertainty, the two independent data sets provide consistent measurements of cosmological parameters. It is less clear whether the bias model of REDMAGIC galaxies in SDSS and DES is identical, a question we turn to in the following subsection.

## A. Results on bias and stochasticity

The Bayes factors for the stochastic models, i.e., the ratio of their evidence over the evidence of the *b* only model, are 3.6 (*b*, *r*) and 2.5 (*b*, $\alpha_0$, $\alpha_1$). This means that there is substantial evidence, as defined by the Jeffreys scale, for stochasticity in the count of REDMAGIC galaxies at fixed projected matter density within 20' apertures and with a redshift range of $z_T = 0.2–0.45$. Similar observations are made in SDSS, with Bayes factors 2.8 and 4.5 for the introduction of the stochastic models. The data thus prefers a model with stochasticity, but at an odds ratio of ≈3:1, the preference is not very conclusive.

The DES constraint on *r* is $r = 0.77_{-0.13}^{+0.10}$. In likelihood runs of the (*b*, *r*) model on log-normal mocks with no stochasticity (Appendix D), we find smaller central values for *r* than this in 3 out of 40 independent realizations.

We note that this finding is *not* in conflict with the nondetection of stochasticity in the DES Collaboration *et al.* [1] 3x2pt analysis, and the associated explicit tests

<footnote n="3">³We use the version of the likelihood that does not vary neutrino mass, as in our counts and lensing in cells analysis. The galaxy bias parameter plotted is the mean bias of REDMAGIC galaxies in the first two bins ($z = 0.15–0.3$ and $0.3–0.45$), weighted 1:2, which is not quite the same as the bias of our single $z_T = 0.2–0.45$ lens bin.</footnote>
<figure><img src="image_14.png" /><figcaption>FIG. 11. Constraints on $S_8 = \sigma_8\sqrt{\Omega_m/0.3}$, bias and stochasticity of the redmagic tracer galaxies from density split lensing and counts-in-cells in DES Y1 (blue) and SDSS (red).</figcaption></figure>

for consistency of the clustering and galaxy-galaxy lensing constraints on bias [30,37]. Those analyses use significantly larger scales ($> 27'$ and $> 45'$ in the lowest lens redshift bins, corresponding to $> 8h^{-1}$ comoving Mpc for clustering, and $> 12h^{-1}$ comoving Mpc for galaxy-galaxy lensing), on which stochasticity, if present, is expected to be small. Our statistic in sensitive to stochasticity on scales smaller or equal to the radius $\theta_T = 20'$ of the apertures inside which we count tracer galaxies. Physically, this corresponds to $\lesssim 3.5–7h^{-1}$ comoving Mpc in the $z_T = 0.2–0.45$ redshift range. Uncertainty as to whether the nonstochastic bias model would be sufficient on these smaller scales was a primary reason for the conservative 3x2pt scale cuts [82].

Figure 11 shows constraints on $S_8 \equiv \sigma_8\sqrt{\Omega_m/0.3}$, bias and stochasticity of the tracer galaxies, in both DES and SDSS. The deviation of $r$ from unity is at the $\approx 2\sigma$ level. Individual parameter constraints are consistent, while there is a hint for a lower value of $r$ in SDSS at fixed cosmology. Note that the primary degeneracy of $S_8$ is not with bias, but with stochasticity—even a mildly informative prior on $r$ would significantly lower uncertainty on $S_8$. For cosmological constraints close to those of the 3x2pt analysis ($S_8 = 0.80 \pm 0.02$ for the run with fixed $\nu$ mass), the counts and lensing in cells data is well fit by small stochasticity.

For the more complex $\alpha_0, \alpha_1$ model of density dependent stochasticity, which is not significantly preferred to $r$ by the data anywhere, the situation is qualitatively similar. The SDSS data does not constrain these parameters very well, especially $\alpha_1$ (cf. Table II),

but there is an indication of super-Poissonian scatter in galaxy count at fixed matter density, the amplitude of which increases with density, broadly consistent with the effect of a single stochasticity parameter $r$ ([28] their Fig. 6).

It is difficult to compare this tentative detection of stochasticity on the $\lesssim 7h^{-1}$ comoving Mpc aperture scale to the literature. Various works have found levels of stochasticity that are broadly consistent, using a range of samples and scales in numerical simulations [e.g., [83,84]] and data [e.g., [85–90]]. The comparison of low-z galaxy clustering and galaxy-galaxy lensing in DES SV on scales above $4h^{-1}$ comoving Mpc provided similar hints of $r < 1$ ([91,92], see also [93]). Even those studies that found no evidence for $r \neq 1$ do not exclude a mild stochasticity on the relevant scales within their uncertainties [94,95]. Most of these studies use two-point correlations, which means their results on stochasticity would have to be transformed to aperture statistics using a numerical model or simulations.

Note that we do not attempt to combine the DES and SDSS results because, without more detailed study, it is not certain that the redmagic samples trace the exact same galaxy populations. A larger stochasticity of redmagic galaxies in SDSS, if at all significant, could also be due to correlations of the redmagic density with SDSS observational systematics that, unlike in the case of DES [30], has not been removed.

## B. Test for excess skewness of matter density

As described in Sec. VA, we can allow for the skewness of the projected, smoothed matter density field, $S_3$, to be a free parameter in our likelihood, rather than predicting it from perturbation theory. Parameter constraints from these likelihood runs are shown in Table III.

We first test whether the introduction of this additional parameter to our model is justified by the data. The Bayes factor of the extended models with $\Delta S_3$ as a free parameter, relative to any of the three models for the connection of galaxies and matter with fixed $S_3$, is smaller than unity, both on DES and on SDSS runs. This indicates no evidence that such an extension is required.

If we still perform a likelihood analysis of the extended models despite of this, we can find constraints on $S_3$. For the $(b, r)$ model, these are shown in Fig. 12. DES Y1 and SDSS provide independent constraints, both of which are consistent with $\Delta S_3 = 0$. The DES constraint, $\Delta S_3/S_3 = -0.08_{-0.20}^{+0.25}$ is significantly tighter, primarily due to the fact that the lensing signal that breaks the degeneracy of bias and skewness is measured with higher signal-to-noise ratio.

Generalizing the likelihood to a two-parametric $\alpha_0, \alpha_1$ model for stochasticity and leaving $S_3$ free yields similarly tight constraints on from DES data, again consistent with no excess skewness at $\Delta S_3 = -0.18_{-0.22}^{+0.25}$. In SDSS, the
<figure><img src="image_15.png" /><figcaption>FIG. 12. Constraints on excess skewness of the projected matter density field $\Delta S_3$, amplitude of structure $S_8 = \sigma_8 \sqrt{\Omega_m/0.3}$, bias and stochasticity of the REDMAGIC tracer galaxies from density split lensing and counts-in-cells in DES Y1 (blue) and SDSS (red). Bayesian model comparison indicates that the introduction of $S_3$ as a free parameter, rather than fixing it from perturbation theory calculations (the $\Delta S_3 = 0$ indicated by dashed lines), is not necessary to describe the data.</figcaption></figure>

posterior distribution of $\alpha_1$ is not constrained in this model within our sampling range.

The joint interpretation of these results is that we find no hints for an excess or deficit in skewness of the matter density relative to our $\Lambda$CDM perturbation theory prediction. This conclusion is largely independent of the bias model we choose, and tested at the 20 per cent level. Future analyses with larger data sets or joint constraints from counts and lensing in cells and additional probes could provide much tighter constraints on $S_3$.

## VII. CONCLUSIONS

We perform the first cosmological analysis using counts and lensing in cells, a method that constrains the matter density PDF with the combination of counts-in-cells and gravitational lensing signals around low and high density

lines of sight. We do this by creating quintiles based on the galaxy counts in apertures and evaluating the stacked lensing for each quintile.

This analysis is tested extensively, using the perturbation theory model of Friedrich et al. [28], by applying it to log-normal density fields and realistic N-body mock catalogs from the *Buzzard* simulations. Robustness to systematics in the data and choices in the analysis is confirmed by a series of tests, performed while blind to the cosmological parameter values constrained by our data.

Applying the analysis to data vectors from DES and SDSS, we find that

(i) The data prefer stochasticity (beyond Poisson sampling) in galaxy count at fixed matter density, on the < 20' or < 3.5–7h$^{-1}$ comoving Mpc smoothing scale of our aperture. This is indicated
by a Bayesian model comparison of bias parame-
trizations, whose odds ratios of ≈3–4:1 in
favor of stochastic models represent substantial
evidence on the Jeffreys scale. Our data does not
discriminate between the two different models of
stochasticity we apply, one with the correlation
coefficient of galaxy and matter density r as a free
parameter, and one with two parameters α₀, α₁ that
allow for a density dependent deviation from
Poissonian noise.

(ii) Either of these models yields consistent DES constraints on the cosmic matter density, $\Omega_m = 0.26_{-0.03}^{+0.04}$ for the (b, r) model and $\Omega_m = 0.28_{-0.04}^{+0.05}$ for the (b, $\alpha_0$, $\alpha_1$) model. These are consistent with, and not much less constraining than, results from the three tomographic auto- and cross-two-point correlation functions of galaxy counts and gravitational shear (3x2pt) in DES Collaboration et al. [1], $\Omega_m = 0.26_{-0.03}^{+0.02}$.

(iii) The degeneracy direction of $\Omega_m$ and $\sigma_8$ and best-fit values of $\sigma_8$ depend on the stochasticity model—for (b, r) we find a higher amplitude $\sigma_8 = 0.97_{-0.06}^{+0.07}$ and $\sigma_8 = 0.86_{-0.05}^{+0.06}$ from DES and SDSS than for the (b, $\alpha_0$, $\alpha_1$) model. The latter has a central value of $\sigma_8 = 0.80$ and somewhat larger uncertainty. In turn, this means that external constraints on cosmology would yield an improved model selection and posterior on stochasticity and vice versa. For small stochasticity, both models are consistent with each other and with the 3x2pt constraint. The assumption of linear bias with no stochasticity, which is however mildly disfavored by the data, would allow constraints from the matter PDF that are competitive with 3x2pt. Thus if we could use prior knowledge to select a particular model or a narrower range of possible bias parameter values, this would greatly improve the constraining power on cosmological parameters.

(iv) Because counts and lensing in cells measure the width and skewness of the matter PDF independ-
ently, they can be used to constrain $S_3 = \langle\delta^2\rangle^2/\langle\delta^3\rangle$, for which $\Lambda$CDM and gravitational collapse make a
very stable prediction. We find that the deviation
from this prediction, $\Delta S_3/S_3$, is consistent with zero,
within 20-30 per cent uncertainty, in any bias and
stochasticity model.

This analysis of the density PDF opens several avenues
for future research:

(i) Towards tightly constraining cosmological parame-
ters and performing model independent tests on
higher moments of the matter density PDF generated
by gravity—this could be achieved best with ex-
ternal priors or data on stochasticity or bias, or in the

regime of larger scales or larger tracer density, where
predicted signals are less sensitive to the bias model.

(ii) Towards discriminating between and constraining parameters of bias and stochasticity models—optimally by analyzing smaller scales, which will require improved models of the matter density PDF and baryonic effects on it.

(iii) Towards joint analyses, e.g., with two-point func-
tions or CMB lensing data, that can break cosmo-
logical, galaxy bias parameter and nuisance
parameter degeneracies.

Given a suitable model and analysis framework, counts and lensing in cells with present and imminently available data will allow tight constraints on cosmological parameters and the hierarchy of moments of the matter density field.

ACKNOWLEDGMENTS

D. G. thanks Yao-Yuan Mao, Cora Uhlemann, Zvonimir Vlah, and numerous members of the DES WL, LSS and Theory working groups for helpful discussions. 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 NAS8-03060. O. F. acknowledges funding by SFB-Transregio 33 ‘The Dark Universe’ by the Deutsche Forschungsgemeinschaft (DFG) and the DFG Cluster of Excellence ‘Origin and Structure of the Universe.’ 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 Collaborating Institutions are Argonne National Laboratory, the University of California at Santa Cruz, the University of Cambridge, Centro de Investigaciones Energéticas, 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, Texas A&M University, and the OzDES Membership Consortium. Based in part on observations at Cerro Tololo Inter-American Observatory, National Optical Astronomy Observatory, which is operated by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation. The DES data management system is supported by the National Science Foundation under Grant No. AST-1138766 and No. AST-1536171. The DES participants from Spanish institutions are partially supported by MINECO under Grants No. AYA2015-71825, No. ESP2015-88861, No. FPA2015-68048, No. SEV-2012-0234, No. SEV-2016-0597, and No. MDM-2015-0509, some of which include ERDF funds from the European Union. I. F. A. E. is partially funded by the CERCA program of the Generalitat de Catalunya. Research leading to these results has received funding from the European Research Council under the European Union's Seventh Framework Program (FP7/2007-2013) including ERC Grant agreements No. 240672, No. 291329, and No. 306478. We acknowledge support from the Australian Research Council Centre of Excellence for All-sky Astrophysics (CAASTRO), through Project No. CE110001020. S. H. acknowledges support by the DFG cluster of excellence 'Origin and Structure of the Universe' (http://www.universe-cluster.de). This manuscript has been authored by Fermi Research Alliance, LLC under Contract No. DE-AC02-07CH11359 with the U.S. Department of Energy, Office of Science, Office of High Energy Physics. The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a nonexclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this manuscript, or allow others to do so, for United States Government purposes. This paper has gone through internal review by the DES collaboration.

## APPENDIX A: BERNOULLI MASKING OF A COUNT WITH POISSONIAN NOISE

Assume that the true number of galaxies N in a randomly selected volume follows a Poisson distribution around the

expectation value $\bar{N}$ (which is conditional on the matter density inside the volume),

$$P(N) = \frac{e^{-\bar{N}}}{N!} \times \bar{N}^N. \qquad (A1)$$

If incompleteness and masking act on each galaxy with some detection probability p, then the observed number of galaxies in that volume $N_{obs}$ is related to N by a Bernoulli process,

$$P(N_{\text{obs}}|N) = \frac{N!}{N_{\text{obs}}!(N - N_{\text{obs}})!} p^{N_{\text{obs}}} (1-p)^{N - N_{\text{obs}}}. \qquad (A2)$$

The resulting probability distribution for $N_{obs}$,

$$\begin{aligned} P(N_{\text{obs}}) &= \sum_{N} P(N_{\text{obs}}|N)P(N) \\ &= \frac{e^{-\bar{N}}}{N_{\text{obs}}!} (p\bar{N})^{N_{\text{obs}}} \sum_{N=N_{\text{obs}}}^{\infty} \frac{[(1-p)\bar{N}]^{N-N_{\text{obs}}}}{(N-N_{\text{obs}})!} \\ &= \frac{e^{-p\bar{N}}}{N_{\text{obs}}!} (p\bar{N})^{N_{\text{obs}}}, \end{aligned} \qquad (A3)$$

is again Poisson-distributed around the expectation value $p\bar{N}$. Masking-induced scatter in $N_{obs}$ is therefore correctly described by Poisson noise. This justifies the accounting for masking used in Sec. II A. Note that we have implicitly assumed intrinsic Poisson noise and independent random removal of galaxies due to masking, when in fact clustering implies a more complex form of Eqs. (A1) and (A2), yet at a level not relevant at first order.

## APPENDIX B: CHOICE OF LOG-NORMAL PARAMETERS FOR THE SIMULATED DENSITY AND CONVERGENCE FIELDS

In this Appendix, we describe how we chose the log-normal parameters (i.e., the minimum allowed values of the log-normal PDFs) and power spectra for generating simulated convergence and density fields that closely match the 2-point and 3-point auto- and cross-correlation statistics we expect from our fiducial model.

### 1. Configuration of FLASK maps of simulated density contrast

As we have shown in Friedrich et al. [28], at a scale of $\theta_A = 20'$ the PDF of the smoothed matter density contrast $\delta_{m,T}$ is well described by a zero-mean shifted log-normal distribution, when the parameters of the log-normal PDF are chosen such as to match the variance and skewness of $\delta_{m,T}$. For the cosmological parameters and redshift distributions of REDMAGIC galaxies in the Buzzard
simulations, the variance and perturbation theory prediction for skewness result in a log-normal shift parameter of $\delta_0 = 0.669$.

We input this parameter and an angular power spectrum computed by Limber's approximation from our fiducial matter power spectrum to the FLASK tool to generate maps of projected density contrast. FLASK will generate HEALPIX maps of $\delta_{m,2D}$ such that on the pixel scale of these maps $\delta_{m,2D}$ is a log-normal random variable with $\delta_0 = 0.669$. In our case, the pixel scale is much smaller than the smoothing scale $\theta = 20'$ for which we determined the log-normal parameter. Fortunately, a limit theorem derived by Szyszkowicz and Yanikomeroglu [96] ensures that a smoothed version of a log-normal random field, while not formally log-normal, is still well described by a log-normal field with the same shift parameter. We have verified this approximation to be accurate in our situation.

## 2. Configuration of FLASK maps of simulated convergence fields

For the convergence field $\kappa$, defined in Eq. (3.7), we again need to fix the parameters of the log-normal simulations. The expectation value of $\kappa_\theta$ around overdense or underdense lines of sight with given matter contrast $\delta_{m,T}$ is fully determined by the moments $\langle \delta_{m,T}^n \rangle$ and joint moments $\langle \delta_{m,T}^n \kappa_\theta \rangle$, $n \ge 1$.

Hence, the expectation value of the density split lensing signal obtains contributions only from those redshifts where the distribution of tracers and the lensing efficiency kernel $W_s(w)$ overlap. The covariance of the signal, however, also has contributions from foreground and background structures at distances where only $W_s(w)$ is non-zero.

As pointed out by, e.g., Xavier et al. [49], 2D projections of the 3D density contrast such as $\delta_{m,T}$ and $\kappa_\theta$ are not well described by a joint log-normal distribution if the kernels of their projection have a very different width along the line-of-sight. In our situation, the lensing kernel $W_s(w)$ is significantly broader than the distribution in co-moving distance of our tracer galaxies. In order to still accurately match the higher-order statistics of $\delta_{m,\theta}$ and $\kappa_\theta$ with the FLASK log-normal simulations, we split $\kappa$ into a contribution from the overlap of tracers and $W_s(w)$ and a contribution from foreground and background structures, i.e.,

$$ \kappa(\hat{n}) = \kappa_{\text{overlap}}(\hat{n}) + \kappa_{\text{nonoverlap}}(\hat{n}) \quad (B1) $$

with

$$ \kappa_{\text{overlap}}(\hat{n}) = \int_{w_{\min}}^{w_{\max}} dw W_s(w) \delta_{m,3D}(w\hat{n}, w) \quad (B2) $$

and

$$ \kappa_{\text{nonoverlap}}(\hat{n}) = \int_{0}^{w_{\min}} dw W_s(w) \delta_{m,3D}(w\hat{n}, w) + \int_{w_{\max}}^{\infty} dw W_s(w) \delta_{m,3D}(w\hat{n}, w). \quad (B3) $$

Here, $w_{\min}$ and $w_{\max}$ are the minimum and maximum comoving distances of our tracer population. We separately compute the power spectra of $\kappa_{\text{overlap}}$ and $\kappa_{\text{nonoverlap}}$ using Limber's approximation and our fiducial matter power spectrum.

Instead of approximating the distribution of $\kappa_{\text{overlap}}$ and $\delta_{m,T}$ with a joint log-normal distribution, we point out in Friedrich et al. [28] that it is better to further split $\kappa_{\text{overlap}}$ into two contributions as

$$ \kappa_{\text{overlap}} = \kappa_{\text{log-normal}} + \kappa_{\text{uncorr}}, \quad (B4) $$

where only $\kappa_{\text{log-normal}}$ is a log-normal variable and $\kappa_{\text{uncorr}}$ is assumed to be completely uncorrelated with $\delta_{m,T}$. The reason for this is the following: The density split lensing signal mainly depends on the moments $\langle \delta_{m,T}^2 \rangle$ and $\langle \delta_{m,T}^3 \rangle$ as well as $\langle \delta_{m,T} \kappa_{\text{overlap}} \rangle$ and $\langle \delta_{m,T}^2 \kappa_{\text{overlap}} \rangle$. Requiring our simulated convergence fields to obey our analytic predictions of these moments would already fix the log-normal PDF for $\kappa_{\text{overlap}}$. Importantly, it would also fix the variance $\langle \kappa_{\text{overlap}}^2 \rangle$. And in general, this variance will disagree with the variance of $\kappa_{\text{overlap}}$ as defined in Eq. (B2) that is predicted from our power spectrum. Splitting $\kappa_{\text{overlap}}$ into $\kappa_{\text{log-normal}}$ and $\kappa_{\text{uncorr}}$ solves this disagreement, since we can use the above moments to fix the log-normal PDF of $\kappa_{\text{log-normal}}$ and attribute part of the variance of $\kappa_{\text{overlap}}$ to $\kappa_{\text{uncorr}}$ to keep the total variance in agreement with our power spectrum. We then use FLASK to simulate the contributions to $\kappa_{\text{overlap}}$ as two distinct random fields. We assume that $\kappa_{\text{log-normal}}$ and $\kappa_{\text{uncorr}}$ are uncorrelated and that their power spectra are simply proportional to that of $\kappa_{\text{overlap}}$. The proportionality factors are determined such that the two power spectra sum up to the total power spectrum of $\kappa_{\text{overlap}}$ and also such that the variance of $\kappa_{\text{log-normal}}$ is indeed the one predicted by its log-normal PDF.

Finally, while in Friedrich et al. [28] we allow the log-normal PDF to vary depending on the scale $\theta$ of the convergence field, we have to choose a fixed scale in order to generate log-normal random fields with FLASK. We consider $\theta = 20'$ as a reasonable choice. At larger scales the log-normal PDF of our formalism quickly transitions to a Gaussian PDF anyway, and we do not consider smaller scales in our analysis. The log-normal shift parameters we get this way are $\kappa_0 = 0.0088, 0.0150,$
<table><caption>TABLE III. Constraints from counts and lensing in cells likelihood runs with skewness of the matter density field $S_3$ as a free parameter. Bayes factors are quoted relative to the b, r model without free $S_3$. SDSS does not constrain the two-parametric stochasticity model $\alpha_0, \alpha_1$ within our sampling range.</caption><thead><tr><th>Data</th><th>Model</th><th>Bayes factor</th><th>S<sub>8</sub></th><th>Ω<sub>m</sub></th><th>σ<sub>8</sub></th><th>b</th><th>r</th><th>α<sub>0</sub></th><th>α<sub>1</sub></th><th>ΔS<sub>3</sub>/S<sub>3</sub></th></tr></thead><tbody><tr><td>DES</td><td>b, r, ΔS<sub>3</sub></td><td>0.3</td><td>0.91<sup>+0.10</sup><sub>-0.10</sub></td><td>0.26<sup>+0.07</sup><sub>-0.05</sub></td><td>0.96<sup>+0.17</sup><sub>-0.13</sub></td><td>1.37<sup>+0.32</sup><sub>-0.27</sub></td><td>0.72<sup>+0.14</sup><sub>-0.10</sub></td><td>...</td><td>...</td><td>-0.08<sup>+0.25</sup><sub>-0.20</sub></td></tr><tr><td>SDSS</td><td>b, r, ΔS<sub>3</sub></td><td>0.4</td><td>0.76<sup>+0.12</sup><sub>-0.09</sub></td><td>0.28<sup>+0.09</sup><sub>-0.07</sub></td><td>0.72<sup>+0.17</sup><sub>-0.13</sub></td><td>1.64<sup>+0.44</sup><sub>-0.46</sub></td><td>0.73<sup>+0.15</sup><sub>-0.15</sub></td><td>...</td><td>...</td><td>0.06<sup>+0.40</sup><sub>-0.27</sub></td></tr><tr><td>DES</td><td>b, α<sub>0</sub>, α<sub>1</sub>, ΔS<sub>3</sub></td><td>0.3</td><td>0.80<sup>+0.06</sup><sub>-0.05</sub></td><td>0.26<sup>+0.07</sup><sub>-0.05</sub></td><td>0.86<sup>+0.10</sup><sub>-0.13</sub></td><td>1.48<sup>+0.41</sup><sub>-0.32</sub></td><td>...</td><td>1.7<sup>+0.4</sup><sub>-0.6</sub></td><td>2.0<sup>+1.1</sup><sub>-0.9</sub></td><td>-0.18<sup>+0.25</sup><sub>-0.22</sub></td></tr><tr><td>SDSS</td><td>b, α<sub>0</sub>, α<sub>1</sub>, ΔS<sub>3</sub></td><td>0.6</td><td>0.78<sup>+0.07</sup><sub>-0.08</sub></td><td>0.27<sup>+0.09</sup><sub>-0.05</sub></td><td>0.80<sup>+0.14</sup><sub>-0.12</sub></td><td>0.98<sup>+0.49</sup><sub>-0.17</sub></td><td>...</td><td>2.5<sup>+0.3</sup><sub>-0.5</sub></td><td>>2.1 (68% c.l.)</td><td>-0.14<sup>+0.44</sup><sub>-0.39</sub></td></tr></tbody></table>

0.0181 for our three DES source redshift bins and 0.0094 for the SDSS sources. The remaining ingredient needed by FLASK to generate the field $\kappa_{\log-normal}$ is its cross power spectrum with $\delta_{m,2D}$. Since we assume all other contributions to the convergence to be uncorrelated to the tracer density contrast, this is just the cross power spectrum of the total convergence and $\delta_{m,2D}$. Also, the theorem by Szyszkowicz and Yanikomeroglu [96] ensures that generating log-normal fields with the above values of $\kappa_0$ at the pixel scale is sufficient to obtain the same log-normal properties also on larger smoothing scales.

The contributions $\kappa_{uncorr}$ and $\kappa_{nonoverlap}$ enter the covariance of our signal via 2-point statistics only. Hence, we just use the sum of their power spectra to generate a single Gaussian random field, i.e., no log-normal shift parameters have to be determined for these components. Since we include 3 source bins in our analysis, we also have to include their cross-power spectra in our FLASK configuration. The cross-power spectra of the nonoverlap contributions to the convergence can be computed straightforwardly using Limber's approximation. The cross-power spectra of $\kappa_{overlap}$ between different source bins are split to cross-power spectra of the fields $\kappa_{\log-normal}$ and the fields $\kappa_{uncorr}$ in a similar way as the auto-power spectra for identical source bins: this time we assume that $\delta_{m,2D}$ and the two convergence fields $\kappa_{\log-normal,i}$ and $\kappa_{\log-normal,j}$ corresponding to source bins i and j have a joint log-normal distribution. We then compute their combined third order moment with the method presented in Friedrich *et al.* [28],

which fixes this PDF and allows us to compute the covariance of $\kappa_{\log-normal,i}$ and $\kappa_{\log-normal,j}$. The cross power spectra are then split between the $\kappa_{\log-normal}$ and $\kappa_{uncorr}$ contributions such that they add up to the total cross power spectrum of $\kappa_{overlap,i}$ and $\kappa_{overlap,j}$ while also giving the correct covariance between the $\kappa_{\log-normal}$ contributions in each source bin. No cross-correlation is assumed for $\kappa_{\log-normal}$ and $\kappa_{uncorr}$ between any combination of source bins.

## APPENDIX C: DEPENDENCE OF SIGNAL ON REDMAGIC VARIANT

As discussed in Sec. II B 1, there are two variants of REDMAGIC catalogs in DES Y1 (based on either MAG_AUTO or MOF photometry), and either can optionally be corrected with a set of weights that removes the correlation of galaxy density with observational systematics. Our fiducial choice, as in Elvin-Poole *et al.* [30] for the redshift range we use, is the MAG_AUTO catalog with these weights applied. In this appendix, we repeat our measurements with the remaining three variants of the REDMAGIC catalogs to see whether there are appreciable differences in the recovered signals.

Results from this are shown in Table IV and, for selected samples, in Fig. 13. We find that the ratios of shear signals between the variants are consistent with unity, with a hint of lower signals in the MOF variant when not applying the weights to correct for systematics related density variations. This is consistent with the analysis of Elvin-Poole *et al.*

<table><caption>TABLE IV. Shear ratios of REDMAGIC variants.</caption><thead><tr><th>z<sub>s</sub> bin</th><th>Variant</th><th>γ<sub>t</sub> ratio-1 [10<sup>-2</sup>]</th><th>γ<sub>t,0</sub> ratio-1 [10<sup>-2</sup>]</th><th>γ<sub>t,4</sub> ratio-1 [10<sup>-2</sup>]</th></tr></thead><tbody><tr><td>0.43–0.63</td><td>COADD, no corr.</td><td>-0.2 ± 2.7</td><td>0.7</td><td>-1.1</td></tr><tr><td>0.63–0.90</td><td>COADD, no corr.</td><td>0.6 ± 1.2</td><td>2.0</td><td>0.9</td></tr><tr><td>0.90–1.30</td><td>COADD, no corr.</td><td>-1.2 ± 1.3</td><td>-0.2</td><td>-2.2</td></tr><tr><td>0.43–0.63</td><td>MOF, corr.</td><td>-4.1 ± 7.1</td><td>7.5</td><td>-15.7</td></tr><tr><td>0.63–0.90</td><td>MOF, corr.</td><td>-2.2 ± 3.9</td><td>-5.7</td><td>-1.2</td></tr><tr><td>0.90–1.30</td><td>MOF, corr.</td><td>-1.2 ± 4.3</td><td>5.5</td><td>-7.8</td></tr><tr><td>0.43–0.63</td><td>MOF, no corr.</td><td>-6.3 ± 7.4</td><td>3.6</td><td>-16.2</td></tr><tr><td>0.63–0.90</td><td>MOF, no corr.</td><td>-3.5 ± 4.1</td><td>-7.3</td><td>0.4</td></tr><tr><td>0.90–1.30</td><td>MOF, no corr.</td><td>-2.7 ± 4.7</td><td>3.6</td><td>-9.0</td></tr></tbody></table>
<figure><img src="image_16.png" /><figcaption>FIG. 13. Ratio of tangential shear measured around troughs (cyan) and overdense lines of sight (red), identified using different variants of the REDMAGIC catalog. The fiducial variant is based on co-add MAG_AUTO photometry and corrects tracer density for its correlation with observational systematics. The top and bottom panel compare this to REDMAGIC runs on the same photometry and multi-epoch MOF photometry, both without correcting for the correlation with systematics.</figcaption></figure>

[30], who found larger, significant correlations of galaxy density with observational systematics in this catalog.

## APPENDIX D: MOCK ANALYSIS OF LOG-NORMAL SIMULATION

We perform a set of tests of our model prediction code, our covariance estimation and likelihood pipeline using log-normal realizations of our data vector.

We first check that the $\chi^2$ distributions of residuals of the individual log-normal realizations follow expectations. To this end, we add realizations of shape noise from random rotations of our actual source catalog to the log-normal data vectors. We have used $N_{\text{real}} = 960$ out of 1000 realizations of the sum of cosmic variance and shape noise to estimate the covariance matrix, and for these the mean $\chi^2$ when evaluated with the inverse of the estimated covariance matrix is equal, within uncertainties, to the number of degrees of freedom. For the additional 40 realizations that

we have reserved for the purpose of testing, the mean $\chi^2$ is larger by a factor consistent with the inverse of Eq. (4.1) [65,97]. We apply this factor to our estimated covariance matrix for all following analyses. We also confirm that this correction is *not* appropriate when estimating the cosmic variance and shape noise parts of our covariance independently and coadding them in the end, which is why we do not apply this (statistically desirable) procedure [but cf. [67]] for a possible way out].

The mean data vector of all 1000 log-normal realizations is described well by the analytical prediction of Friedrich *et al.* [28]. The total $\chi^2$ of the residual of the mean vs. the model, at our fiducial scale cuts and for the DES analysis with three source redshift bins with 208 d.o.f., evaluated at the input cosmology and galaxy bias and the covariance matrix of Sec. IV, is 0.19, without the increase in the diagonal of the counts-in-cells covariance discussed in Sec. IV C. This is well below the statistical uncertainty ($\sqrt{\text{Var}[\chi^2_{208 \text{ d.o.f.}}]} > 20$) and confirms that the analytical calculations used in Friedrich *et al.* [28] to get from input power spectra and skewness parameters for the matter and convergence fields to a prediction of our signals are derived and implemented correctly and to sufficient precision in our codes.

Next, we run a full likelihood analysis on the 40 reserved realizations to check the coverage, i.e., whether the 68 per cent confidence contour contains the input cosmology in a sufficient fraction of cases. Since nuisance parameters for measurement systematics are fixed to zero in these data vectors, we only vary $b, \alpha_0, \alpha_1, \Omega_m, \sigma_8$, and in an additional run $\Delta S_3$ in these analyses.

Figure 14 shows the resulting confidence contours. The fraction of times that these contain the true input cosmology (horizontal/vertical lines) is consistent with 68 percent. This statement is also true when adding the additional parameter $\Delta S_3$. The mean half-width of the marginalized confidence intervals for $(\Omega_m, \sigma_8)$ are (0.05, 0.07) and for $(\Omega_m, \sigma_8, \Delta S_3)$ are (0.07, 0.10, 0.28).

We repeat the test for the $(b, r)$ model, finding consistent coverage for $\sigma_8$ and marginally low coverage for $\Omega_m$, potentially related to the asymmetry of the flat prior around $r$, which is correlated with $\Omega_m$ (and which cannot physically take a value larger than 1, its true value in the log-normal simulations).

Finally, to test our methodology of model comparison with Bayes factors, we sample these likelihoods again with the MULTINEST algorithm. We first confirm that parameter constraints derived from MULTINEST and EMCEE closely match each other. We then determine ratio of Bayesian evidence between complex models and the most simple run (linear bias with no stochasticity, $\Delta S_3 = 0$). Since the log-normal simulations should be fully described by this simple model, there should be no evidence for any more complex model. We show in Fig. 15 that this is indeed the case.
<figure><img src="image_17.png" /><figcaption>FIG. 14. Realizations of likelihoods on 40 log-normal mocks, sampled with EMCEE, with the input parameters given by the dashed black lines. The coverage of these mock likelihoods is within expectations, i.e., any input parameter lies within its marginalized $1\sigma$ confidence interval ≈68 per cent of the times.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>FIG. 15. Model comparison based on the Bayes factor, i.e., the ratio of evidences of the more complex over the less complex model, from likelihoods run on log-normal mocks. These mocks contain no stochasticity and no $\Delta S_3$ and hence should not favor the more complex models. Dotted and solid vertical lines indicate the Jeffreys scale for substantial and strong evidence, showing that this is indeed the case.</figcaption></figure>

## APPENDIX E: MOCK ANALYSIS OF BUZZARD SIMULATIONS

We repeat the steps of Appendix D, but instead of log-normal realizations using density split lensing and counts-in-cells signals measured on the *Buzzard* suite of N-body simulations.

The *Buzzard* simulations are a suite of mock galaxy catalogs built on top of dark matter only N-body simulations. We describe only the essential details here and refer the reader to more complete descriptions in DeRose *et al.* [98]. Each set of 6 *Buzzard* DES Y1 catalogs is generated from a combination of 3 N-body lightcones runs using L-Gadget2, a version of Gadget2 [99] optimized for memory efficiency. 2nd order Lagrangian perturbation theory initial conditions were generated using 2LPTIC [100]. The lightcones were produced on the fly as the simulations ran from boxes of volumes $1050^3$, $2600^3$, and $4000^3$ ($h^{-1}$ Mpc)$^3$ and mass resolutions of $2.7 \times 10^{10}$, $1.3 \times 10^{11}$, $4.8 \times 10^{11} h^{-1} M_{\odot}$ respectively. The light cones are joined at redshifts 0.34
and 0.9, and arranged such that the highest resolution
simulations are used at lower redshifts.

Galaxy catalogs are produced from the dark matter
lightcones using the ADDGALS algorithm [101].
ADDGALS uses the relation between large scale density
and r-band absolute magnitude determined from a subhalo
abundance matching [102–104] algorithm run on a high
resolution N-body simulation to place galaxies with magni-
tudes into a low resolution density field. Galaxies are then
assigned SEDs from SDSS DR7 [105] based on the distance
to their fifth nearest neighbor. DES griz fluxes are generated
from these SEDs and photometric noise is added to them
using the DES Y1 depth map. These galaxies are then lensed
using the multiple plane raytracing algorithm CALCLENS
[106] which uses a spherical harmonic transform Poisson
solver allowing for curved sky boundary conditions. Two
different versions of ADDGALS were used to create the
catalogs referred to as *Buzzard-v1.1* and *Buzzard-v1.6*
below. The main differences come from changes to the
assumed luminosity function of galaxies, the evolution of the
red fraction of galaxies with redshift and the use of different
depth maps to produce photometric noise.

In order to obtain a simulated REDMAGIC galaxy
sample, we run the REDMAGIC algorithm on the simu-
lations with the same configuration as the data, yielding
very similar photometric redshift and clustering properties
to those found in DES Y1. A METACALIBRATION like
sample is created by making signal to noise cuts in order
to approximate the source density found in the data.

While the true input cosmology is known and identical
between versions 1.1 and 1.6 of the *Buzzard* mock catalogs,
their model for early-type galaxy SEDs vary, and thus the
REDMAGIC selection is not the same. Consequently, we
cannot use the values for bias and stochasticity parameters
determined in Friedrich *et al.* [[28], their Secs. 4.3.1 and
4.3.2] for version 1.1 as truth inputs for version 1.6. We
repeat the analysis performed there to find that REDMAGIC
galaxies in *z* = 0.2–0.45 have a somewhat larger bias and
stochasticity in version 1.6 (*b*, *α*₀, *α*₁ = 1.72, 1.36, 0.29)
than in version 1.1 (*b*, *α*₀, *α*₁ = 1.54, 1.26, 0.29). The same
is true in the simpler parametrization with (*b*, *r* = 1.84, 0.96
vs. 1.62,0.96). We note that the finding that these galaxies in
the *Buzzard-v1.6* simulation show large-scale stochasticity is
in line with the analysis of MacCrann *et al.* [56].

The masks of the *Buzzard-v1.1* simulations cover a
somewhat smaller sky area (cf. Friedrich *et al.*, their Fig. 1)
than the full DES-SPT footprint (cf. Fig. 1). We therefore
generate a separate covariance matrix for *Buzzard-v1.1* by
cutting the mask of the FLASK realizations accordingly and
repeating the procedure of Sec. IV B. The larger galaxy bias in
*Buzzard-v1.6* causes a larger cosmic variance than the one
derived in Sec. IV for *Buzzard-v1.1* parameters. We account
for this by simply re-scaling the *Buzzard-v1.6* covariance
matrix by a factor 1.05 that brings the mean $\chi^2$ of the *Buzzard-
v.1.6* and v1.1 simulations to agreement.

We split the survey area by density and measure counts-
in-cells and lensing signals in *Buzzard-v1.1* and *Buzzard-
v1.6* mock catalogs as described in Sec. II. As source
galaxies, we use the approximated lensing source sample
described above. For source redshift tomography, we
define three bins based on BPZ redshift expectation
values at $z_s = 0.43-0.70, 0.70-0.78, 0.92-1.30$ in v1.1
and $z_s = 0.43-0.66, 0.66-0.76, 0.76-1.00$. These limits
were defined such that the mean true redshift matches that
of the three highest redshift source bins in Hoyle *et al.* [33].
In measuring tangential shear profiles, we use true shear
information from raytracing in *Buzzard* [107]. Shape noise
is added at the level of the data vector using measurements
of density split lensing made with randomly rotated source
galaxy catalogs on DES Y1 data itself that were not used
for estimating the covariance matrix (Sec. IVA). We
generate 62 data vectors in total, 28 (7 different shape
noise realizations each) based on the 4 *Buzzard-v1.1*

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

FIG. 16. Residuals of measurements of lensing and counts-in-cells signal in *Buzzard* relative to model evaluated at true parameters of the simulations. Top panel: $\chi^2$ of each realization with 208 entries in the data vector. Bottom panel: Residuals of individual data points in units of their expected standard deviation. Black lines indicate expected distributions. The figures use 62 realizations based on 21 independent N-body mocks from the *Buzzard-v1.1* and *v1.6* runs, to which different realizations of Y1 METACALIBRATION shape noise were added.
$N$-body mocks and 34 (2 different shape noise realizations each) from 17 *Buzzard-v1.6* mocks. Caution must be taken in interpreting results, especially on coverage, due to the fact that these realizations are not completely independent.

For a test of our model and the covariance matrix, we predict the signal, using the model of Friedrich *et al.* [28], with the true input parameters for cosmology and the directly measured parameters for the $b, \alpha_0, \alpha_1$ model of galaxy biasing. In this, we use the mean true redshift distributions of REDMAGIC and lensing source galaxies, and the mean REDMAGIC counts in the v1.1 or v1.6 sets of *Buzzard* simulations. For the source bias of the lowest redshift bin, we assume $b_{s1} = 0.5$.

From the residuals of this model relative to the 62 realizations with 208 data points each, we generate the statistics shown in Fig. 16. The distribution of $\chi^2$ is consistent with expectations (left panel). The root-mean-square of individual data point residuals in units of the

expected standard deviation according to the covariance matrix is slightly larger than unity, but by less than 5 per cent (right panel). We have confirmed that in the case of the *Buzzard-v1.1* simulations it is indeed consistent with unity, and suspect that the small increase in *Buzzard-v1.6* might be due to the larger cosmic variance caused by the increased REDMAGIC galaxy bias, which we have not fully accounted for in the covariance matrix (see discussion above). There is thus no apparent non-Gaussianity in the distribution of residuals [108]. These statistics look very similar when using the $b, r$ instead of the $b, \alpha_0, \alpha_1$ model, i.e., the prediction evaluated at the true cosmology and directly measured stochasticity is a good fit to the simulations.

For a test of our inference methodology, we run mock likelihood chains on these *Buzzard* realizations, with $\Omega_m, \sigma_8, b$, stochasticity parameters, and galaxy bias of sources in the lowest redshift bin $b_{s1}$ as free parameters. All other parameters are fixed to their input values, and we use true redshift distributions for the predictions.

<figure><img src="image_20.png" /><figcaption>FIG. 17. Realizations of likelihoods on 34 mock realizations based on the *Buzzard-v1.6* catalogs with the true parameters given by the dashed black lines. The coverage of these mock likelihoods is within expectations, i.e., any input parameter lies within its marginalized $1\sigma$ confidence interval ≈68 per cent of the times. Likelihoods are sampled with MULTINEST, which is a significant improvement in speed. While the appearance is more patchy, results are consistent with long EMCEE chains, as we have checked from a subset of dual runs on identical simulations.</figcaption></figure>
<figure><img src="image_21.png" /></figure>

FIG. 18. Model comparison based on the Bayes factor, i.e., the ratio of evidences of the more complex over the less complex model, from likelihoods run on *Buzzard-v1.6* mocks. These mocks contain no $\Delta S_3$ but do show a complex relation of galaxy count and matter density and hence may favor the more complex models. Dotted and solid vertical lines indicate the Jeffreys scale for substantial and strong evidence, showing that this is indeed the case.

Figure 17 shows the most relevant run, of the fiducial $b, \alpha_0, \alpha_1$ biasing model on the 17 independent *Buzzard-v1.6* realizations with 2 independent versions of shape noise each. Coverage of cosmological parameters is within expectations for 17 independent realizations (58 per cent for $\Omega_m$, 62 per cent for $\sigma_8$). We find a mean source bias of $b_{s1} = 0.62$.

For the $b, r$ model, in contrast, coverage of $\sigma_8$ is low (only 18 percent) because $\sigma_8$ is biased high (with the best fit being 0.93 on average, significantly above the true input 0.82). We confirm that this is not only caused by the asymmetry of the allowed parameter range around the directly measured $r = 0.96$: allowing linear extrapolation of the model to $r > 1$ still yields a coverage of only 33 percent in $\sigma_8$. We also confirm that the same low coverage is found regardless of whether we use EMCEE or MULTINEST as a sampler.

Some understanding might be gained from the Bayesian model comparison we perform in Fig. 18. While the *Buzzard* simulations should not prefer a model with free $\Delta S_3$, it is up to tests like this to determine what level of complexity is needed to describe the stochasticity in the galaxy count distribution. Notably, in some cases the $b, \alpha_0, \alpha_1$ model is strongly preferred to the $b, r$ model, and in one case the $b, r$ model prefers the introduction of $\Delta S_3$ as a free parameter, potentially to compensate for the excess skewness in a galaxy count distribution with density-dependent stochasticity. Both could be an indication of density dependent stochasticity actually being present in *Buzzard-v1.6*. In this case, the high bias of $\sigma_8$ in the $b, r$ model seen in the coverage tests could be a partial compensation of the missing variance.

We conclude that the $b, \alpha_0, \alpha_1$ model is required for describing the galaxy distribution in *Buzzard-v1.6*

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

FIG. 19. $\chi^2$ of comparing *Buzzard* realizations of the data vector to a model of the signal in which the skewness of the matter density field was set to 0, i.e., the smoothed matter density field was assumed to be Gaussian. In all realizations, the Gaussian model is excluded at more than $3\sigma$.

sufficiently well. In this framework, the coverage test and model comparison test for adding $\Delta S_3$ as a free parameter have results that are within expectations for a reliable inference scheme. Whether or not the complexity in *Buzzard* is realistic or caused, in part, by peculiarities of the simulations (such as the stitching of simulation boxes or the placement of galaxies) is not clear at this point [56].

Finally, we test how well the *Buzzard* realizations of our data vector are fit in a model in which the smoothed matter field $\delta_{m,T}$ is assumed to be Gaussian with the variance predicted by the power spectrum, i.e., to have no skewness. This is a model that provided a reasonable fit to the data in Gruen et al. [23]. Figure 19 shows the $\chi^2$ of the residuals of data and Gaussian model, where we have used true cosmology and directly measured galaxy bias ($b, \alpha_0, \alpha_1$) parameters for the latter. Even for the least $\chi^2$ among the 62 realizations, the large $\chi^2$ allows us to exclude the Gaussian model at $p < 0.001$.

## APPENDIX F: CLUSTERING CONSTRAINTS ON SOURCE AND LENS REDSHIFT DISTRIBUTIONS

For both our REDMAGIC tracer population and our lensing source galaxies, we calibrate the mean value of redshift distributions using clustering redshifts. In these, excess angular correlations with thinly-sliced spectroscopic or spectroscopic-like samples are used to determine the redshift distributions [38,39,74,75,109,110].

For the case of DES Y1, this calibration and its systematic uncertainties are described in Cawthon et al. [74], Davis et al. [38], Gatti et al. [39]. Here, we give the missing details relevant to the specific tracer galaxy selection in DES and the tracer and source samples in SDSS not described in those papers.
<figure><img src="image_23.png" /><figcaption>FIG. 20. Comparison of the distribution of DES (left panel) and SDSS (right) REDMAGIC photometric redshifts with the estimated redshift distribution from clustering using SDSS DR12 LOWZ and CMASS as a reference sample. By shifting the photo-z distribution to fit the mean of the clustering estimate, we estimate the photometric bias to be $\Delta z = 0.003 \pm 0.008$ (DES) and $\Delta z = 0.002 \pm 0.006$ (SDSS).</figcaption></figure>

## 1. REDMAGIC galaxies

We perform a clustering redshift analysis to assess the accuracy of the REDMAGIC photometric redshift algorithm. In this, we follow the same techniques as described in Cawthon *et al.* [74], yet with a different redshift cut on the REDMAGIC sample.

We cross-correlate the SDSS and DES REDMAGIC high density sample, selected by the photometric estimate $z_{\text{red}} = 0.2-0.45$, with the BOSS spectroscopic samples LOWZ and CMASS [111]. The clustering redshift measurement follows the method described in Schmidt *et al.* [112], using physical scales of 0.5 to 1.5 Mpc. Statistical error estimates are from jackknife resampling.

As described in Cawthon *et al.* [74], the main systematic to this measurement is the redshift evolution of the galaxy bias of the two samples across the redshift bin. We can measure the evolution of the spectroscopic samples using autocorrelations. The amplitude of the autocorrelations of REDMAGIC though will be impacted by photo-z errors. Results in Cawthon *et al.* [74] suggest that the quantity $b_l\sqrt{w_{\text{mm}}}$, where $b_l$ is the galaxy bias of REDMAGIC and $w_{\text{mm}}$ is the autocorrelation of the matter density on the scales we measure, shows very little evolution with redshift for REDMAGIC. We parameterize the galaxy bias evolution with $b_l\sqrt{w_{\text{mm}}} \propto (1+z)^\gamma$. In Cawthon *et al.* [74], we assumed that $\gamma$ is in the range $0.0 \pm 2.0$, but found that it shows less spread around 0 for wider redshift bins. Since the redshift range $z_T = 0.2-0.45$ used in this work is significantly wider than that used in the bins of Cawthon *et al.* [74], we assume $\gamma = 0 \pm 1.5$ here. These choices of $\gamma$ range broadly fit the various estimates of $\gamma$ from the auto-correlations of REDMAGIC on the full DES sample, the Stripe 82 sample which contains the galaxies that overlap with

BOSS, and simulations as measured in Cawthon *et al.* [74].

After the galaxy bias correction is selected, the clustering measurement is also narrowed to $\pm 2.5\sigma$ around the mean of the clustering redshift distribution estimate, with $\sigma$ being the standard deviation of that estimate. This cut is indicated in Figs. 20 and 21. This is done since the clustering measurement can be noisy and biased where the signal is low. We then fit for a single photometric bias parameter, $\Delta z = z - z_{\text{phot}}$, where $z$ is the clustering redshift mean and $z_{\text{phot}}$ is the

<figure><img src="image_24.png" /><figcaption>FIG. 21. Clustering measurements on a subset of the SDSS lens sample that has spectroscopic measurements. (79, 583, 4.9% of the full sample.) The clustering measurement contains a correction for the evolution of galaxy bias, $\gamma = -1.4$ which best fits the true redshift distribution mean. This subsample gives a photo-z bias of $-0.0015$.</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>FIG. 22. Redshift distribution of combined SDSS samples as measured by photo-z (solid line) and clustering (points) with SDSS DR12 spectra Alam et al. [113]. The photo-z is calibrated with a mean redshift offset of $\Delta_z = -0.014 \pm 0.011$, while the clustering estimates are calibrated with a bias evolution model of $b_{\text{PZ}}b_{\text{spectra}} \propto (1+z)^\gamma$ with $\gamma = -2.0 \pm 0.6$.</figcaption></figure>

photometric redshift mean over the redshift range selected by the $\pm 2.5\sigma$ cut.

The results of this analysis are shown in Table V and Figs. 20. The measurement on DES REDMAGIC can only be done on a subsample of 20,347 galaxies in Stripe 82, which has overlap with the BOSS spectroscopic samples. For this reason, the statistical uncertainty on the DES measurement is significantly larger than on SDSS REDMAGIC. For both samples, the dominant systematic error is the uncertainty in the REDMAGIC galaxy bias evolution parameterized by $\gamma$. The uncertainty in $\gamma$ of $\pm 1.5$ leads to approximately an uncertainty of $\Delta z$ of $\pm 0.006$. Figure 21 shows the cross correlation redshift estimate using just a subsample of SDSS REDMAGIC that has spectra. This is a biased sample that is brighter and likely has a different photo-z bias and galaxy bias, though it appears to confirm some similar trends seen in Fig. 20, such as a lower galaxy density around $z = 0.4$ and higher density around $z = 0.3$ compared to the photo-z code.

<table><caption>TABLE V. Cross-correlation estimates of photo-z bias on REDMAGIC.</caption><thead><tr><th>Sample</th><th>Δz</th><th>δΔz (syst)</th><th>δΔz (stat)</th></tr></thead><tbody><tr><td>DES</td><td>0.003 ± 0.008</td><td>0.006</td><td>0.005</td></tr><tr><td>SDSS</td><td>0.002 ± 0.006</td><td>0.006</td><td>0.001</td></tr></tbody></table>

## 2. SDSS source galaxies

We measure the excess angular clustering from 500 to 1500 kpc between SDSS source galaxies selected and weighted as in Sec. II C 2 and the SDSS DR12 spectroscopic galaxy sample. This method is similar to that detailed in Ménard et al. [109]. Statistical errors in clustering-z are estimated by jackknife resampling. We find that reliable clustering signal can only be recovered over a subset of the sample, from $0.32 < z < 0.80$, but this is sufficient to calibrate the photo-z distribution with our method, where we assume the shape of the source $n(z)$ is correctly estimated by photo-z and only determine its shift. More details for the calibration procedure may be found in Davis et al. [38] as well as in Gatti et al. [39] and Cawthon et al. [74].

We fit a relative redshift bias, $z \rightarrow z - \Delta_z$ in the photo-z distribution, and parametrize the clustering bias evolution in the clustering-z as a power-law with free exponent, $b_{\text{PZ}}b_{\text{spectra}} \propto (1+z)^\gamma$. As we found in Gatti et al. [39] we expect that systematic uncertainties in modeling the underlying bias evolution to dominate over our statistical uncertainties. Systematic uncertainties in Gatti et al. [39] are of the order $\sigma_{\Delta z} \approx 0.01$, while our statistical uncertainty is $\sigma_{\Delta z} = 0.002$. In quadrature, these combine to be an uncertainty of $\sigma_{\Delta z} = 0.011$. Clustering constraints on the redshift distribution are shown in Fig. 22. We find $\Delta_z = -0.014 \pm 0.011$ and $\gamma = -2.0 \pm 0.6$, indicating only a marginal preference for a mean offset in redshift, and moderate combined bias evolution. We note that the bias evolution measured is the product of the bias evolutions of the lens and source samples, and so the bias evolution measured here will likely differ from the calibration of the lens bin even though both use the same reference galaxies and scales.

[1] DES Collaboration et al., arXiv:1708.01530.

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

[3] 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).

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

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

[6] J. Carron, Phys. Rev. Lett. **108**, 071301 (2012).

[7] J. Carron and M. C. Neyrinck, Astrophys. J. **750**, 28 (2012).

[8] P. Coles and B. Jones, Mon. Not. R. Astron. Soc. **248**, 1 (1991).

[9] S. Foreman, M. R. Becker, and R. H. Wechsler, Mon. Not. R. Astron. Soc. **463**, 3326 (2016).

[10] N. MacCrann et al., Mon. Not. R. Astron. Soc. **465**, 2567 (2017).
[11] A. Klypin, F. Prada, J. Betancort-Rijo, and F. D. Albareti, arXiv:1706.01909.

[12] B. Pandey, S. D. M. White, V. Springel, and R. E. Angulo, Mon. Not. R. Astron. Soc. **435**, 2968 (2013).

[13] R. Takahashi, M. Oguri, M. Sato, and T. Hamana, Astrophys. J. **742**, 15 (2011).

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

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

[16] C. Uhlemann, S. Codis, C. Pichon, F. Bernardeau, and P. Reimberg, Mon. Not. R. Astron. Soc. **460**, 1529 (2016).

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

[18] J. Liu, J. C. Hill, B. D. Sherwin, A. Petri, V. Böhm, and Z. Haiman, Phys. Rev. D **94**, 103501 (2016).

[19] K. Patton, J. Blazek, K. Honscheid, E. Huff, P. Melchior, A. J. Ross, and E. Suchyta, Mon. Not. R. Astron. Soc. **472**, 439 (2017).

[20] A. J. Ross, R. J. Brunner, and A. D. Myers, Astrophys. J. **682**, 737 (2008).

[21] C. Chang *et al.*, Mon. Not. R. Astron. Soc. **475**, 3165 (2018).

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

[23] D. Gruen *et al.*, Mon. Not. R. Astron. Soc. **455**, 3367 (2016).

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

[25] M. Cautun, E. Paillas, Y.-C. Cai, S. Bose, J. Armijo, B. Li, and N. Padilla, Mon. Not. R. Astron. Soc. **476**, 3195 (2018).

[26] Y. Higuchi and Shirasaki, Mon. Not. R. Astron. Soc. **459**, 2762 (2016).

[27] P. Simon *et al.*, Mon. Not. R. Astron. Soc. **430**, 2476 (2013).

[28] O. Friedrich *et al.*, following paper, Phys. Rev. D **98**, 023508 (2018).

[29] 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).

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

[31] A. Drlica-Wagner *et al.*, The Astrophysical Journal **235**, 33 (2018).

[32] H. T. Diehl *et al.*, in Observatory Operations: Strategies, Processes, and Systems V, 91490V (2014).

[33] B. Hoyle *et al.*, arXiv:1708.01532.

[34] J. Zuntz *et al.*, arXiv:1708.01533.

[35] E. Rozo *et al.*, Mon. Not. R. Astron. Soc. **461**, 1431 (2016).

[36] E. S. Rykoff *et al.*, Astrophys. J. **785**, 104 (2014).

[37] J. Prat *et al.*, arXiv:1708.01537.

[38] C. Davis *et al.*, arXiv:1710.02517.

[39] M. Gatti *et al.*, Mon. Not. R. Astron. Soc. **477**, 1664 (2018).

[40] E. Huff and R. Mandelbaum, arXiv:1702.02600.

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

[42] S. Samuroff *et al.*, Mon. Not. R. Astron. Soc. **475**, 4524 (2018).

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

[44] E. S. Sheldon *et al.*, Astrophys. J. **703**, 2217 (2009).

[45] J. Clampitt and B. Jain, Mon. Not. R. Astron. Soc. **454**, 3357 (2015).

[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] H. S. Xavier, F. B. Abdalla, and B. Joachimi, Mon. Not. R. Astron. Soc. **459**, 3693 (2016).

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

[51] M. Boylan-Kolchin, V. Springel, S. D. M. White, and A. Jenkins, Mon. Not. R. Astron. Soc. **406**, 896 (2010).

[52] F. Jiang and F. C. van den Bosch, Mon. Not. R. Astron. Soc. **472**, 657 (2017).

[53] Y.-Y. Mao, M. Williamson, and R. H. Wechsler, Astrophys. J. **810**, 21 (2015).

[54] J. N. Fry, S. Colombi, P. Fosalba, A. Balaraman, I. Szapudi, and R. Teyssier, Mon. Not. R. Astron. Soc. **415**, 153 (2011).

[55] M. C. Neyrinck, M. A. Aragón-Calvo, D. Jeong, and X. Wang, Mon. Not. R. Astron. Soc. **441**, 646 (2014).

[56] N. MacCrann, J. DeRose *et al.*, arXiv:1803.09795.

[57] R. Mandelbaum, C. M. Hirata, U. Seljak, J. Guzik, N. Padmanabhan, C. Blake, M. R. Blanton, R. Lupton, and J. Brinkmann, Mon. Not. R. Astron. Soc. **361**, 1287 (2005).

[58] E. S. Sheldon *et al.*, Astron. J. **127**, 2544 (2004).

[59] B. Joachimi, M. Cacciato, T. D. Kitching, A. Leonard, R. Mandelbaum, B. M. Schäfer, C. Sifón, H. Hoekstra, A. Kiessling, D. Kirk, and A. Rassat, Space Sci. Rev. **193**, 1 (2015).

[60] M. Troxel and M. Ishak, Phys. Rep. **558**, 1 (2015).

[61] R. Murata, T. Nishimichi, M. Takada, H. Miyatake, M. Shirasaki, S. More, R. Takahashi, and K. Osato, arXiv: 1707.01907.

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

[63] C. Sánchez *et al.*, Mon. Not. R. Astron. Soc. **465**, 746 (2017).

[64] S. Singh, R. Mandelbaum, U. Seljak, A. Slosar, and J. V. Gonzalez, Mon. Not. R. Astron. Soc. **471**, 3827 (2017).

[65] J. Hartlap, P. Simon, and P. Schneider, Astron. Astrophys. **464**, 399 (2007).

[66] S. Dodelson and M. D. Schneider, Phys. Rev. D **88**, 063537 (2013).

[67] O. Friedrich and T. Eifler, Mon. Not. R. Astron. Soc. **476**, 3195 (2018).

[68] E. Sellentin and A. F. Heavens, Mon. Not. R. Astron. Soc. **456**, L132 (2016).

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

[70] P. J. E. Peebles, *The Large-Scale Structure of the Universe* (Princeton University Press, Princeton, New Jersey, 1980).

[71] J. N. Fry, Astrophys. J. **279**, 499 (1984).

[72] R. Juszkiewicz, F. R. Bouchet, and S. Colombi, Astrophys. J. **412**, L9 (1993).

[73] Planck Collaboration *et al.*, Astron. Astrophys. **594**, A17 (2016).

[74] R. Cawthon *et al.*, arXiv:1712.07298.

[75] C. Davis *et al.*, Mon. Not. R. Astron. Soc. **477**, 2196 (2018).
[76] M. Simet, T. McClintock, R. Mandelbaum, E. Rozo, E. Rykoff, E. Sheldon, and R. H. Wechsler, Mon. Not. R. Astron. Soc. **466**, 3103 (2017).

[77] D. Foreman-Mackey, D. W. Hogg, D. Lang, and J. Goodman, Publ. Astron. Soc. Pac. **125**, 306 (2013).

[78] F. Feroz, M. P. Hobson, and M. Bridges, Mon. Not. R. Astron. Soc. **398**, 1601 (2009).

[79] H. Jeffreys, *Theory of Probability*, 3rd ed. (Oxford University Press, New York, 1961).

[80] S. Hinton, The Journal of Open Source Software **1**, (2016).

[81] R. Andrae, T. Schulze-Hartung, and P. Melchior, arXiv: 1012.3754.

[82] E. Krause et al., arXiv:1706.09359.

[83] S. Bonoli and U. L. Pen, Mon. Not. R. Astron. Soc. **396**, 1610 (2009).

[84] R. S. Somerville, G. Lemson, Y. Sigad, A. Dekel, G. Kauffmann, and S. D. M. White, Mon. Not. R. Astron. Soc. **320**, 289 (2001).

[85] H. Hoekstra, L. van Waebeke, M. D. Gladders, Y. Mellier, and H. K. C. Yee, Astrophys. J. **577**, 604 (2002).

[86] A. Leauthaud et al., Mon. Not. R. Astron. Soc. **467**, 3024 (2017).

[87] P. Simon, M. Hetterscheidt, M. Schirmer, T. Erben, P. Schneider, C. Wolf, and K. Meisenheimer, Astron. Astrophys. **461**, 861 (2007).

[88] M. E. C. Swanson, M. Tegmark, M. Blanton, and I. Zehavi, Mon. Not. R. Astron. Soc. **385**, 1635 (2008).

[89] Y. Wang, X. Yang, and H. J. Mo, and F. C. van den Bosch, Astrophys. J. **664**, 608 (2007).

[90] V. Wild et al., Mon. Not. R. Astron. Soc. **356**, 247 (2005).

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

[92] J. Prat et al., Mon. Not. R. Astron. Soc. **473**, 1667 (2018).

[93] T. Giannantonio et al., Mon. Not. R. Astron. Soc. **456**, 3213 (2016).

[94] J. Comparat et al., Mon. Not. R. Astron. Soc. **433**, 1146 (2013).

[95] E. Jullo, J. Rhodes, A. Kiessling, J. E. Taylor, R. Massey, J. Berge, C. Schimd, J.-P. Kneib, and N. Scoville, Astrophys. J. **750**, 37 (2012).

[96] S. S. Szyszkowicz and H. Yanikomeroglu, IEEE Trans. Commun. **57**, 3538 (2009).

[97] T. W. Anderson, *An Introduction to Multivariate Statistical Analysis*, 3rd ed. (Wiley, New York, 2003).

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

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

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

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

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

[103] B. V. Lehmann, Y.-Y. Mao, M. R. Becker, S. W. Skillman, and R. H. Wechsler, Astrophys. J. **834**, 37 (2017).

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

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

[106] M. R. Becker, Mon. Not. R. Astron. Soc. **435**, 115 (2013).

[107] M. R. Becker, Ph.D. thesis, The University of Chicago, 2013.

[108] E. Sellentin and A. F. Heavens, Mon. Not. R. Astron. Soc. **473**, 2355 (2018).

[109] B. Ménard, R. Scranton, S. Schmidt, C. Morrison, D. Jeong, T. Budavari, and M. Rahman, arXiv:1303.4722.

[110] J. A. Newman, Astrophys. J. **684**, 88 (2008).

[111] K. S. Dawson et al., Astron. J. **145**, 10 (2013).

[112] S. J. Schmidt, B. Ménard, R. Scranton, C. Morrison, and C. K. McBride, Mon. Not. R. Astron. Soc. **431**, 3307 (2013).

[113] S. Alam et al., Astrophys. J. Suppl. Ser. **219**, 12 (2015).