# The Splashback Feature around DES Galaxy Clusters: Galaxy Density and Weak Lensing Profiles

Chihway Chang, E. Baxter, B. Jain, C. Sánchez, S. Adhikari, T.N. Varga, Y. Fang, Eduardo Rozo, E.S. Rykoff, Andrey Kravtsov, et al.

► To cite this version:

Chihway Chang, E. Baxter, B. Jain, C. Sánchez, S. Adhikari, et al. The Splashback Feature around DES Galaxy Clusters: Galaxy Density and Weak Lensing Profiles. The Astrophysical Journal, 2018, 864 (1), pp.83. 10.3847/1538-4357/aad5e7. hal-01872309

HAL Id: hal-01872309

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

Submitted on 19 Aug 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
<figure><img src="image_2.png" /></figure>

# The Splashback Feature around DES Galaxy Clusters: Galaxy Density and Weak Lensing Profiles

C. Chang¹, E. Baxter², B. Jain², C. Sánchez²,³, S. Adhikari⁴,⁵, T. N. Varga⁶,⁷, Y. Fang², E. Rozo⁸, E. S. Rykoff⁵,⁹, A. Kravtsov¹⁰,¹¹,¹², D. Gruen⁵,⁹, W. Hartley¹³, E. M. Huff¹⁴, M. Jarvis², A. G. Kim¹⁵, J. Prat³, N. MacCrann¹⁶,¹⁷, T. McClintock⁸, A. Palmese¹³, D. Rapetti¹⁸,¹⁹, R. P. Rollins²⁰, S. Samuroff²⁰, E. Sheldon²¹, M. A. Troxel¹⁶,¹⁷, R. H. Wechsler⁵,⁹,²², Y. Zhang²³, J. Zuntz²⁴, T. M. C. Abbott²⁵, F. B. Abdalla¹³,²⁶, S. Allam²³, J. Annis²³, K. Bechtol²⁷, A. Benoit-Lévy¹³,²⁸,²⁹, G. M. Bernstein², D. Brooks¹³, E. Buckley-Geer²³, A. Carnero Rosell³⁰,³¹, M. Carrasco Kind³²,³³, J. Carretero³, C. B. D'Andrea², L. N. da Costa³⁰,³¹, C. Davis⁵, S. Desai³⁴, H. T. Diehl²³, J. P. Dietrich³⁵,³⁶, A. Drlica-Wagner²³, T. F. Eifler¹⁴,³⁷, B. Flaugher²³, P. Fosalba³⁸, J. Frieman¹,²³, J. García-Bellido³⁹, E. Gaztanaga³⁸, D. W. Gerdes⁴⁰,⁴¹, R. A. Gruendl³²,³³, J. Gschwend³⁰,³¹, G. Gutierrez²³, K. Honscheid¹⁶,¹⁷, D. J. James⁴², T. Jeltema⁴³, E. Krause⁵, K. Kuehn⁴⁴, O. Lahav¹³, M. Lima³⁰,⁴⁵, M. March², J. L. Marshall⁴⁶, P. Martini¹⁶,⁴⁷, P. Melchior⁴⁸, F. Menanteau³²,³³, R. Miquel³,⁴⁹, I. J. J. Mohr⁷,³⁵,³⁶, B. Nord²³, R. L. C. Ogando³⁰,³¹, A. A. Plazas¹⁴, E. Sanchez⁵⁰, V. Scarpine²³, R. Schindler⁹, M. Schubnell⁴¹, I. Sevilla-Noarbe⁵⁰, M. Smith⁵¹, R. C. Smith²⁵, M. Soares-Santos²³, F. Sobreira³⁰,⁵², E. Suchyta⁵³, M. E. C. Swanson³³, G. Tarle⁴¹, and J. Weller⁶,⁷,³⁵

(DES Collaboration)

¹ Kavli Institute for Cosmological Physics, University of Chicago, Chicago, IL 60637, USA; chihway@kicp.uchicago.edu

² Department of Physics and Astronomy, University of Pennsylvania, Philadelphia, PA 19104, USA

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

⁴ Department of Astronomy, University of Illinois at Urbana-Champaign, Champaign, IL 61801, USA

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

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

⁷ Max Planck Institute for Extraterrestrial Physics, Giessenbachstrasse, D-85748 Garching, Germany

⁸ Department of Physics, University of Arizona, Tucson, AZ 85721, USA

⁹ SLAC National Accelerator Laboratory, Menlo Park, CA 94025, USA

¹⁰ Kavli Institute for Cosmological Physics, The University of Chicago, Chicago, IL 60637, USA

¹¹ Department of Astronomy and Astrophysics, The University of Chicago, Chicago, IL 60637, USA

¹² Enrico Fermi Institute, The University of Chicago, Chicago, IL 60637, USA

¹³ Department of Physics & Astronomy, University College London, Gower Street, London, WC1E 6BT, UK

¹⁴ Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Dr., Pasadena, CA 91109, USA

¹⁵ Lawrence Berkeley National Laboratory, 1 Cyclotron Road, Berkeley, CA 94720, USA

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

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

¹⁸ Center for Astrophysics and Space Astronomy, Department of Astrophysical and Planetary Science, University of Colorado, Boulder, CO 80309, USA

¹⁹ NASA Ames Research Center, Moffett Field, CA 94035, USA

²⁰ Jodrell Bank Center for Astrophysics, School of Physics and Astronomy, University of Manchester, Oxford Road, Manchester, M13 9PL, UK

²¹ Brookhaven National Laboratory, Bldg. 510, Upton, NY 11973, USA

²² Department of Physics, Stanford University, 382 Via Pueblo Mall, Stanford, CA 94305, USA

²³ Fermi National Accelerator Laboratory, P.O. Box 500, Batavia, IL 60510, USA

²⁴ Institute for Astronomy, University of Edinburgh, Edinburgh EH9 3HJ, UK

²⁵ 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

²⁷ LSST, 933 North Cherry Avenue, Tucson, AZ 85721, USA

²⁸ CNRS, UMR 7095, Institut d'Astrophysique de Paris, F-75014, Paris, France

²⁹ Sorbonne Universités, UPMC Univ Paris 06, UMR 7095, Institut d'Astrophysique de Paris, F-75014, Paris, France

³⁰ Laboratório Interinstitucional de e-Astronomia—LIneA, Rua Gal. José Cristino 77, Rio de Janeiro, RJ—20921-400, Brazil

³¹ Observatório Nacional, Rua Gal. José Cristino 77, Rio de Janeiro, RJ—20921-400, Brazil

³² Department of Astronomy, University of Illinois, 1002 W. Green Street, Urbana, IL 61801, USA

³³ National Center for Supercomputing Applications, 1205 West Clark St., Urbana, IL 61801, USA

³⁴ Department of Physics, IIT Hyderabad, Kandi, Telangana 502285, India

³⁵ Excellence Cluster Universe, Boltzmannstr. 2, D-85748 Garching, Germany

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

³⁷ Department of Physics, California Institute of Technology, Pasadena, CA 91125, USA

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

³⁹ Instituto de Física Teorica UAM/CSIC, Universidad Autónoma de Madrid, E-28049 Madrid, Spain

⁴⁰ Department of Astronomy, University of Michigan, Ann Arbor, MI 48109, USA

⁴¹ Department of Physics, University of Michigan, Ann Arbor, MI 48109, USA

⁴² Astronomy Department, University of Washington, Box 351580, Seattle, WA 98195, USA

⁴³ Santa Cruz Institute for Particle Physics, Santa Cruz, CA 95064, USA

⁴⁴ Australian Astronomical Observatory, North Ryde, NSW 2113, Australia

⁴⁵ Departamento de Física Matemática, Instituto de Física, Universidade de São Paulo, CP 66318, São Paulo, SP, 05314-970, Brazil

⁴⁶ George P. and Cynthia Woods Mitchell Institute for Fundamental Physics and Astronomy, and Department of Physics and Astronomy, Texas A&M University, College Station, TX 77843, USA

⁴⁷ Department of Astronomy, The Ohio State University, Columbus, OH 43210, USA

⁴⁸ Department of Astrophysical Sciences, Princeton University, Peyton Hall, Princeton, NJ 08544, USA

⁴⁹ Institució Catalana de Recerca i Estudis Avançats, E-08010 Barcelona, Spain

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

⁵¹ School of Physics and Astronomy, University of Southampton, Southampton, SO17 1BJ, UK
**Abstract**

Splashback refers to the process of matter that is accreting onto a dark matter halo reaching its first orbital apocenter and turning around in its orbit. The clustercentric radius at which this process occurs, $r_{sp}$, defines a halo boundary that is connected to the dynamics of the cluster. A rapid decline in the halo profile is expected near $r_{sp}$. We measure the galaxy number density and weak lensing mass profiles around REDMAPPER galaxy clusters in the first-year Dark Energy Survey (DES) data. For a cluster sample with mean $M_{200m}$ mass $\approx 2.5 \times 10^{14} M_{\odot}$, we find strong evidence of a splashback-like steepening of the galaxy density profile and measure $r_{sp} = 1.13 \pm 0.07 h^{-1}$ Mpc, consistent with the earlier Sloan Digital Sky Survey measurements of More et al. and Baxter et al. Moreover, our weak lensing measurement demonstrates for the first time the existence of a splashback-like steepening of the matter profile of galaxy clusters. We measure $r_{sp} = 1.34 \pm 0.21 h^{-1}$ Mpc from the weak lensing data, in good agreement with our galaxy density measurements. For different cluster and galaxy samples, we find that, consistent with ΛCDM simulations, $r_{sp}$ scales with $R_{200m}$ and does not evolve with redshift over the redshift range of 0.3–0.6. We also find that potential systematic effects associated with the REDMAPPER algorithm may impact the location of $r_{sp}$. We discuss the progress needed to understand the systematic uncertainties and fully exploit forthcoming data from DES and future surveys, emphasizing the importance of more realistic mock catalogs and independent cluster samples.

**Key words:** cosmology: observations – galaxies: clusters: general – galaxies: halos – gravitational lensing: weak

## 1. Introduction

The density profiles of dark matter halos in N-body simulations exhibit a steepening at radii comparable to the halo virial radius (Diemer & Kravtsov 2014, hereafter DK14). Such a feature was predicted by the analytical collapse models of Gunn & Gott (1972), Fillmore & Goldreich (1984), and Bertschinger (1985). The sharp decline in the profile can be understood as resulting from an absence of particles orbiting beyond the radius of second turnaround.<fnref n="54" /> In simulations, a phase-space caustic cleanly separates matter that is experiencing second turnaround from matter that is on first infall, leading to a very sharp steepening in the halo profile (DK14; Adhikari et al. 2014; Diemer et al. 2017). As measurements from individual halos are noisy, to detect this sharp steepening, one needs to “stack,” or average, over a large number of halos. This makes the caustic structure less clear, since the dark matter halos are oftentimes nonspherical. Nevertheless, DK14 showed that some of the steepening can, in principle, be detected based on simulations. Mansfield et al. (2017) later found that one can improve on the stacking procedure by accounting for the effects of subhalos, which sharpens the steepening of the profile even more. This feature—which appears as a narrow minimum in the logarithmic derivative of the halo density profile—has been termed *splashback*.

The splashback feature is potentially interesting for several reasons. First, it defines a physical boundary of a dark matter halo that is motivated by dynamics (More et al. 2015). This is different from other common halo boundary definitions, such as $R_{200m}$ (the radius within which the mean density is 200 times the mean density of the universe at that redshift), which need not be associated with any change in physical properties across the boundary. Furthermore, the location of the splashback feature has been shown in simulations to correlate with the halo

accretion rate (DK14; Diemer et al. 2017). Since the feature is, in principle, straightforward to measure in data, it could potentially be used to constrain halo accretion rates of clusters, which are otherwise challenging to measure. Finally, the sharpness of the feature and the relatively simple dynamics that are responsible for its generation make it a potentially powerful probe of new physics, such as dark matter self-interaction.

The first measurement of the splashback feature in data was performed by More et al. (2016, hereafter M16) using DR8 data from the Sloan Digital Sky Survey (SDSS; Aihara et al. 2011). M16 measured the projected galaxy density profiles, $\Sigma_g(R)$, around galaxy clusters in the REDMAPPER catalog of Rykoff et al. (2014), finding evidence for a sharp minimum in the logarithmic derivatives of these profiles. Note, however, that a minimum in the logarithmic derivative does not by itself constitute evidence for a splashback feature. Indeed, if the matter profile of the halo is described by a Navarro–Frenk–White (NFW; Navarro et al. 1996) profile at small scales and the halo–matter correlation function at large scales, there will necessarily be a minimum in the logarithmic derivative in the transition regime. In the language of the halo model (for a review, see Cooray & Sheth 2002), a minimum of the logarithmic derivative is naturally associated with the transition regime between the one-halo and two-halo terms. Defining the splashback feature is a way to isolate this dynamical feature that is not explicitly described by either the one-halo or two-halo term. The splashback process produces a profile that is significantly steeper at this transition region than what is expected from the above naive picture of an NFW profile plus the halo–matter correlation function.

By fitting different models to the measured $\Sigma_g$ profiles, M16 determined that the data show strong evidence of the existence of a splashback feature. However, M16 determined that the location of the splashback feature (henceforth the splashback radius) measured in SDSS data appears to be smaller than that predicted by dark matter–only N-body simulations. The explanation for this discrepancy remains unclear; M16 considered several possibilities, including dark matter

<footnote n="52">⁵² Instituto de Física Gleb Wataghin, Universidade Estadual de Campinas, 13083-859, Campinas, SP, Brazil</footnote>

<footnote n="53">⁵³ Computer Science and Mathematics Division, Oak Ridge National Laboratory, Oak Ridge, TN 37831, USA

*Received 2017 October 23; revised 2018 July 23; accepted 2018 July 23; published 2018 August 31*</footnote>

<footnote n="54">⁵⁴ The radius of first turnaround is the radius at which a particle first separates from the Hubble flow and begins to fall toward an overdensity. The radius of second turnaround is the radius at which a particle that has passed by the halo once turns around in its orbit.</footnote>
self-interaction. Two follow-up studies (Zu et al. 2016; Busch & White 2017) examined potential systematic effects in the estimation of the splashback radius and showed that projection along the line of sight could affect the estimated splashback radius, especially when employing a selection based on $\langle R_{\text{mem}} \rangle$ (the weighted average member distance to the cluster center; see more discussion in Section 7), as was done in M16.

On the other hand, Baxter et al. (2017, hereafter B17) pointed out a difficulty associated with quantifying the evidence for a splashback feature using the model parameterizations of DK14 and M16. In particular, these parameterizations rely on a truncated Einasto profile (Einasto 1965), with the truncation term representing the splashback feature. However, the Einasto model is sufficiently flexible that even without such a truncation, it can still reproduce a splashback-like steepening in the outer halo profile. Consequently, the evidence for splashback quantified either with a $\Delta\chi^2$ (as in More et al. 2016) or with a Bayesian evidence ratio (as in B17) can be misleadingly low, even when there is significant steeping of the outer halo profile. This problem becomes more severe when additional flexibility is introduced to the model to account for halo mis-centering, for example.

A more robust approach proposed by B17 is to instead use the model fits to separate the contributions to the total profile from infalling and collapsed material. The logarithmic derivative of the profile of the collapsed material can then be used to identify the presence of a splashback feature in the density profile. The SDSS clusters show a dramatic steepening of the collapsed material profile slightly outside the virial radius, consistent with the presence of a splashback feature. B17 also measured the profiles of red and blue galaxies around REDMAPPER clusters, showing that the fraction consisting of red galaxies exhibits a sharp transition at the splashback radius. This is consistent with the interpretation of star formation being quenched in galaxies that have orbited through the cluster and adds additional support to the picture of a physically motivated halo boundary.

Measurement of the splashback radius has also been used recently by Adhikari et al. (2016, hereafter A16) to measure dynamical friction in galaxy clusters. Dynamical friction refers to an effective drag force induced on a massive object via gravitational interaction with nearby matter (Chandrasekhar 1949; Binney & Tremaine 2008). As pointed out by M16, dynamical friction will act to reduce the splashback radii of galaxy clusters, since a subhalo that has experienced dynamical friction will turn around after first infall sooner compared to a subhalo not experiencing dynamical friction. A16 tested this hypothesis using a sample of lower-richness clusters from SDSS, since the impact of dynamical friction on a fixed galaxy sample is expected to be larger for low-mass parent halos. As the effect of dynamical friction increases with subhalo mass, one expects more massive subhalos to have smaller splashback radii. Indeed, A16 found that the splashback radius identified using a bright galaxy sample (which is expected to live in more massive subhalos) was smaller than the splashback radius identified using a fainter galaxy sample. Two caveats to this simple picture are that galaxies of different magnitudes may have different orbits through the cluster and that quenching may result in changes to galaxy magnitudes within the cluster.

In this work, we measure the galaxy density and weak lensing mass profiles around galaxy clusters in data from the

first year of the Dark Energy Survey (DES Y1). We carry out analyses based on the methodology developed in B17 to characterize the splashback feature. This new data set provides several advances over previous measurements. First, the DES footprint maps a different part of the sky from that of SDSS, providing an independent measurement of a nonoverlapping sample. Second, the DES data extend to higher redshift ranges and fainter galaxies than SDSS. Finally, the DES Y1 weak lensing measurements have significantly higher signal-to-noise compared to the previous lensing measurements in Umetsu & Diemer (2017), which were only able to place a lower bound on the location of the splashback radius. While the signal-to-noise of the lensing measurements is still lower than that of the galaxy density measurements, lensing has the advantage of directly probing the mass profile of the halos, in contrast to the galaxy density profile, which makes it cleaner to compare with dark matter simulations.

The paper is organized as follows. In Section 2, we introduce the data set used in this work and the selection criteria of the samples. In Section 3, we outline the model that is used to describe the observed galaxy and weak lensing profiles around clusters. In Section 4, the measurement methods for both the galaxy and the lensing profiles are described, followed by the model-fitting procedure and a summary of the model priors. In Section 5, we present the fiducial measurements of the galaxy and lensing profiles, as well as a comparison with the simulations. We then investigate in Section 6 the topic of dynamical friction by looking at the splashback feature for galaxy samples of different luminosities. We discuss the potential systematic effects in the measurements in Section 7 and conclude in Section 8.

If not otherwise specified, throughout the paper, we assume a flat $\Lambda$CDM cosmology with $h = 0.7$ and $\Omega_m = 0.3$. In addition, all calculations and plots use comoving coordinates.

## 2. Data

The measurements in this work are based on the DES Y1 data (Diehl et al. 2014). Here we briefly describe the relevant catalogs used, including the REDMAPPER galaxy cluster catalog, the photometric galaxy catalog, the weak lensing shear catalogs, and the photometric redshift (photo-z) catalog.

### 2.1. The REDMAPPER Galaxy Cluster Catalog

We use a galaxy cluster catalog constructed using the REDMAPPER algorithm described in Rykoff et al. (2014, 2016). REDMAPPER is a red-sequence cluster finder that is optimized for large-scale optical surveys such as DES. The same algorithm was employed to construct the SDSS-based cluster catalog used in M16 and B17.

The “fiducial” cluster catalog used in this work is constructed from a volume-limited sample (similar to that described in Rykoff et al. 2016) with a redshift selection of $0.2 < z < 0.55$ and a richness selection of $20 < \lambda < 100$. We only use the Y1 region at decl. $< -35^\circ$, where most of the tests for the weak lensing and photo-z catalogs were conducted. The cluster redshifts used for selection are determined by the REDMAPPER algorithm and are expected to have mean uncertainties $\sigma_z \sim 0.01(1+z)$ (Rykoff et al. 2016). We also consider different subsets of the fiducial sample, as well as a “high-z” sample that is not contained in the fiducial sample. The characteristics of all the samples used in this paper are
<table><caption>**Table 1**
Selection Criteria and Sample Sizes for the Cluster Samples Used in This Work</caption><thead><tr><th>Sample</th><th>z Selection</th><th>λ Selection</th><th>No. of Clusters</th><th>⟨z⟩</th><th>⟨λ⟩</th><th>⟨R<sub>λ</sub>⟩ (h<sup>-1</sup> Mpc)</th><th>⟨M<sub>200m</sub>⟩ (10<sup>14</sup> M<sub>c</sub>)</th></tr></thead><tbody><tr><td>Fiducial</td><td>0.2 &lt; z &lt; 0.55</td><td>20 &lt; λ &lt; 100</td><td>3684</td><td>0.41</td><td>31.6</td><td>0.79</td><td>2.5</td></tr><tr><td>Low-z</td><td>0.2 &lt; z &lt; 0.4</td><td>20 &lt; λ &lt; 100</td><td>1588</td><td>0.32</td><td>32.2</td><td>0.79</td><td>2.5</td></tr><tr><td>Mid-z</td><td>0.4 &lt; z &lt; 0.55</td><td>20 &lt; λ &lt; 100</td><td>2096</td><td>0.48</td><td>31.1</td><td>0.78</td><td>2.5</td></tr><tr><td>High-z</td><td>0.55 &lt; z &lt; 0.7</td><td>20 &lt; λ &lt; 100</td><td>1518</td><td>0.61</td><td>30.3</td><td>0.78</td><td>2.4</td></tr><tr><td>Low-λ</td><td>0.2 &lt; z &lt; 0.55</td><td>20 &lt; λ &lt; 28</td><td>1964</td><td>0.41</td><td>23.3</td><td>0.75</td><td>1.8</td></tr><tr><td>High-λ</td><td>0.2 &lt; z &lt; 0.55</td><td>28 &lt; λ &lt; 100</td><td>1720</td><td>0.40</td><td>41.1</td><td>0.83</td><td>3.3</td></tr></tbody></table>

Note. The mean $M_{200m}$ mass is derived via the mass–richness relation of Melchior et al. (2017).

summarized in Table 1, while the redshift–richness distribution of the fiducial sample is shown in Figure 1. The mean mass for each sample listed in Table 1 is calculated from the mass–richness relation derived in Melchior et al. (2017) for RED-MAPPER clusters identified in DES science verification (SV) data. The REDMAPPER algorithm is expected to be approximately survey-independent. However, small differences in the Melchior et al. (2017) mass–richness relation for SV data and the mass–richness relation of Y1 data may exist due to, e.g., differences in data quality or statistical fluctuations. We rely on the SV mass calibration here because it was derived using DES data; the mass calibration of Y1 REDMAPPER clusters using galaxy lensing is forthcoming. As further support for our use of the SV mass–richness relation from Melchior et al. (2017), we note that Baxter et al. (2018) performed a mass calibration of DES Y1 REDMAPPER clusters using gravitational lensing of the cosmic background radiation, finding excellent consistency with the Melchior et al. (2017) results. We also list in Table 1 the mean $R_λ = (λ/100)^{0.2} h^{-1}$ Mpc values in physical units, which are used later in modeling the cluster mis-centering (Section 3.1).

We see that the number of clusters falls steeply with richness and increases by a factor of ~2 over the redshift range. The structure in the redshift distribution is associated with the DES filter transition and the 4000 Å break, where the photo-z redshifts are less certain. We overlay the histogram with the theoretical expectation of the number of halos given mass and redshift using the halo mass function from Tinker et al. (2008) and the mass–richness relation from Melchior et al. (2017). The data roughly follow the expectation, with higher discrepancies in the low-richness bins. We also make use of the REDMAPPER random catalogs, which uniformly sample the volume over which a real cluster could have been observed. As described in Section 4, the random catalog is used to estimate the background mean galaxy distribution in the absence of galaxy clusters.

Finally, uncertainty in the cluster center position is important in this analysis. According to Rykoff et al. (2016), roughly 22% of the clusters are mis-centered at about 0.3 $R_λ$. We discuss in Section 3.1 how this is incorporated into our model.

## 2.2. The Photometric Galaxy Catalog

We use photometric galaxies from the DES Y1 Gold catalog described in Drlica-Wagner et al. (2018). Our galaxy selection begins with a flux-limited sample of $i < 21.5$ with the following flag cuts: `flags_badregion = flags.gold = 0;` the following color cuts: `-1 < mag_auto_g-mag_auto_r < 3`, `-1 < mag_auto_r-mag_auto_i < 2.5`, `-1 < mag_auto_i-mag_auto_z < 2`; and the following star-galaxy separation cut: `spread_model_i + (5./3.) spreaderr_model_i > 0.007`. The flag cuts are DES-specific and

<figure><img src="image_3.png" /><figcaption>**Figure 1.** $λ$ and $z$ distributions of the fiducial REDMAPPER sample used in this work. The lower left panel shows the 2D histogram in the $z-λ$ plane, while the upper left and lower right panels show the individual 1D histograms of $z$ and $λ$. The color bars in the 2D histogram are shown in log scale. The naive theoretical expectation of the redshift and richness distribution for our sample based on Tinker et al. (2008) and Melchior et al. (2017) is shown by the black dashed curves.</figcaption></figure>

described in Drlica-Wagner et al. (2018), while the other cuts are based on SourceExtractor (Bertin & Arnouts 1996) columns. We further require the errors on the galaxy magnitudes to be less than 0.1. After applying the depth mask, as well as the REDMAPPER mask, the total number of galaxies in this sample is 11,263,383. A random catalog that uniformly samples the galaxy catalog mask is generated from the intersection of the $i > 21.5$ depth mask and the REDMAPPER mask. The final area used is ~1297 deg$^2$.

## 2.3. The Weak Lensing Shear and Photo-z Catalogs

For the lensing measurements performed in this work, we use the two DES Y1 shear catalogs: METACALIBRATION (Huff & Mandelbaum 2017) and IM3SHAPE (Zuntz et al. 2013). Both catalogs are tested and described in detail in Zuntz et al. (2017). The two catalogs were generated using completely independent pipelines; performing the measurements using both catalogs is therefore a powerful test of weak lensing systematics, as shown in, e.g., Troxel et al. (2017). The MetaCalibration catalog
contains 34.8 million galaxies, roughly 60% more than the IM3SHAPE catalog, due to the fact that METACALIBRATION uses the combined information of the r-, i-, and z-band images while IM3SHAPE only uses r-band images. We present our main results using the METACALIBRATION catalog but have checked that the IM3SHAPE measurements show consistent results.

For the weak lensing measurement in this work, redshift information is needed for each source galaxy. We use the photo-z catalog described in Hoyle et al. (2018), which is based on the template-based Bayesian photometric redshifts (BPZ) algorithm (Benítez 2000). Following T. McClintock et al. (2018, in preparation), we use both the mean of the probability density function (PDF) and a random draw from the full PDF for each galaxy when estimating the weak lensing mass profile. We describe the procedure in detail in Section 4.2.

Both the shear calibration biases associated with the shear catalogs and the biases in the photo-z catalog are well characterized in Hoyle et al. (2018) and Zuntz et al. (2017). We do not account for these in the modeling, since it mainly contributes to a scale-independent multiplicative factor at ~2% and does not impact the inference of the splashback feature.

### 3. Formalism

To model the 3D density profile around the clusters, we use the analytical model profile of DK14, which was found to be a good description of dark matter halos in simulations across a wide range of mass, redshift, and accretion rate. The model includes two components: “collapsed” matter that has passed through at least one orbital pericenter and is in orbit around the halo and “infalling” material that is falling toward the halo but has not experienced an orbital pericenter. The profile of the collapsed matter is modeled by a truncated Einasto profile (Einasto 1965), while the infalling material is modeled by a power law. The truncation of the Einasto profile accounts for the splashback feature and is modeled using the $f_{\text{trans}}(r)$ term below.

The complete model for the 3D density, $\rho(r)$, is

$$ \rho(r) = \rho^{\text{coll}}(r) + \rho^{\text{infall}}(r), \quad (1) $$

$$ \rho^{\text{coll}} = \rho^{\text{Ein}}(r) f_{\text{trans}}(r), \quad (2) $$

$$ \rho^{\text{Ein}} = \rho_s \exp \left( -\frac{2}{\alpha} \left[ \left( \frac{r}{r_s} \right)^{\alpha} - 1 \right] \right), \quad (3) $$

$$ f_{\text{trans}}(r) = \left[ 1 + \left( \frac{r}{r_t} \right)^{\beta} \right]^{-\gamma/\beta}, \quad (4) $$

$$ \rho^{\text{infall}} = \rho_0 \left( \frac{r}{r_0} \right)^{-s_e}. \quad (5) $$

Since $r_0$ is completely degenerate with $\rho_0$, we fix $r_0 = 1.5 h^{-1} \text{ Mpc}$ throughout. Also, Equation (5) differs from the formalism in DK14 slightly in that we model the mean-subtracted density profile, so there is no term corresponding to the mean density. Using the mean-subtracted profile allows us to model the average profiles of clusters at different redshifts more easily.

In practice, we measure the 2D projected profile instead of the 3D profile, so it is useful to compute the projected density,

$\Sigma(R)$, which is related to the 3D density by

$$ \Sigma(R) = \int_{-h_{\text{max}}}^{h_{\text{max}}} dh \rho(\sqrt{R^2 + h^2}), \quad (6) $$

where R is the projected distance to the halo center and $h_{\text{max}}$ is the maximum scale of integration. We set $h_{\text{max}} = 40 h^{-1} \text{ Mpc}$ and test in Appendix C that changing $h_{\text{max}}$ does not significantly affect the inferred value of the splashback radius. However, as also shown in Appendix C, changing $h_{\text{max}}$ does impact the inferred slope of the density profile, especially at large distances. Our choice of $h_{\text{max}} = 40 h^{-1} \text{ Mpc}$ is sufficiently large so that increasing $h_{\text{max}}$ by 50% only changes the large-scale slope by <1σ.

The model formulated above was intended for fitting the distribution of mass around halos in simulations. Following M16 and B17, we nevertheless apply the model above to the measured galaxy distribution (replacing $\Sigma(R)$ with $\Sigma_g(R)$, the galaxy density profile), as well as the mass profile. In this approach, the unknown details of galaxy bias are absorbed into the model parameters. In addition, by adopting this model, we have assumed that the average profile of clusters with a range of mass, richness, redshift, and mis-centering parameters (see Section 3.1) can be described by one effective $\Sigma(R)$ or $\Sigma_g(R)$ profile for the whole sample.

From weak lensing shear measurements, we derive the differential mass profile $\Delta\Sigma(R)$,

$$ \Delta\Sigma(R) = \bar{\Sigma}(R) - \Sigma(R), \quad (7) $$

where $\Sigma(R)$ is the projected surface mass density, and $\bar{\Sigma}(R)$ is the average of $\Sigma(R)$ within the circle of radius R, i.e.,

$$ \bar{\Sigma}(R) = \frac{2\pi \int_{0}^{R} dR' R' \Sigma(R')}{\pi R^2}. \quad (8) $$

Equations (7) and (8) make it clear that the lensing profile $\Delta\Sigma(R)$ depends on the density profile of the cluster all the way down to $R = 0$. This is problematic, since the lensing measurements on small scales may be affected by systematics, and the halo density profile may depart from the simple Einasto model at small scales as a result of baryonic effects. It is therefore convenient to introduce a new parameter, $\mu$, defined by

$$ \mu = \int_{0}^{R_{\text{min}}} dR' R' \Sigma(R'), \quad (9) $$

where $R_{\text{min}}$ can be set to the minimum scale at which $\Delta\Sigma$ is measured; we set $R_{\text{min}} = 0.2 h^{-1} \text{ Mpc}$. The expression for $\Delta\Sigma$ can then be written as

$$ \Delta\Sigma(R) = 2 \frac{\mu + \int_{R_{\text{min}}}^{R} dR' R' \Sigma(R')}{R^2} - \Sigma(R). \quad (10) $$

Treating $\mu$ as a free parameter effectively removes any sensitivity of the $\Delta\Sigma$ profile to $\Sigma(R < R_{\text{min}})$, where there are no measurements.

Equations (1)–(5) and (10) use many free parameters to fit functions that are very smooth as a function of radius (see Figures 2 and 3). Consequently, there may be significant degeneracies between the various parameters. We emphasize, however, that our intention here is not to extract robust constraints on the model parameters themselves. Instead, the
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Top: stacked surface density of DES Y1 galaxies around REDMAPPER clusters with $0.2 < z < 0.55$ and $20 < \lambda < 100$ (black points with error bars). The red line shows the model fit to the measurements. The inferred 3D $r_{sp}$ is shown as the vertical orange band, with the width of the band indicating the $1\sigma$ uncertainty; the vertical gray band shows the inferred 2D $R_{sp}$. We note that since the measurement is in projection, the gray band (instead of the orange band) indicates the point of steepest slope of the red line. Bottom: difference in the model and the measurements divided by the uncertainty in the measurement.</figcaption></figure>

goal of the model fitting is mainly to smoothly interpolate between the data points to enable the computation of the logarithmic derivatives of the 3D density profile.

We describe below two additional pieces of modeling that we incorporate into the above formalism in order to capture two important observational complications: cluster mis-centering and weak lensing boost factor.

### 3.1. Cluster Mis-centering

We model the effects of mis-centering following the approach of Melchior et al. (2017) and Baxter et al. (2017). The mis-centered density profile, $\Sigma$, can be related to the profile in the absence of mis-centering, $\Sigma_0$, via

$$ \Sigma = (1 - f_{\text{mis}})\Sigma_0 + f_{\text{mis}}\Sigma_{\text{mis}}, \quad (11) $$

where $f_{\text{mis}}$ is the fraction of clusters that are mis-centered, and $\Sigma_{\text{mis}}$ is the density profile of the mis-centered clusters. For clusters that are mis-centered by $R_{\text{mis}}$ from the true halo center, the corresponding azimuthally averaged density profile is (Yang et al. 2006; Johnston et al. 2007)

$$ \begin{aligned} \Sigma_{\text{mis}}(R|R_{\text{mis}}) \\ = \int_{0}^{2\pi} \frac{d\theta}{2\pi} \Sigma_{0}(\sqrt{R^{2} + R_{\text{mis}}^{2} + 2RR_{\text{mis}}\cos\theta}). \end{aligned} \quad (12) $$

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Top: lensing stacked excess surface mass density around DES REDMAPPER clusters with $0.2 < z < 0.55$ and $20 < \lambda < 100$ (black points with error bars). The red line shows the model fit to the measurements. The inferred 3D $r_{sp}$ is shown as the vertical orange band, with the width of the band indicating the $1\sigma$ uncertainty. Bottom: difference in the model and the measurements divided by the uncertainty in the measurement.</figcaption></figure>

The profile averaged across the distribution of $R_{\text{mis}}$ values is then

$$ \Sigma_{\text{mis}}(R) = \int dR_{\text{mis}} P(R_{\text{mis}}) \Sigma_{\text{mis}}(R|R_{\text{mis}}), \quad (13) $$

where $P(R_{\text{mis}})$ is the probability that a cluster is mis-centered by a (comoving) distance $R_{\text{mis}}$. Following Rykoff et al. (2016), we assume that $P(R_{\text{mis}})$ results from a mis-centering distribution that is a 2D Gaussian on the sky. The 1D $P(R_{\text{mis}})$ is then given by a Rayleigh distribution,

$$ P(R_{\text{mis}}) = \frac{R_{\text{mis}}}{\sigma_{R}^{2}} \exp\left[-\frac{R_{\text{mis}}^{2}}{2\sigma_{R}^{2}}\right], \quad (14) $$

where $\sigma_R$ controls the width of the distribution. Rykoff et al. (2016) assumed $\sigma_R = c_{\text{mis}}R_\lambda$, where $R_\lambda = (\lambda/100)^{0.2} h^{-1} \text{Mpc}$, and used a combination of X-ray and SZ data to measure $\ln(c_{\text{mis}}) = -1.13 \pm 0.22$ and $f_{\text{mis}} = 0.22 \pm 0.11$. We introduce $c_{\text{mis}}$ and $f_{\text{mis}}$ as free parameters in our analysis of both the galaxy density and lensing profiles, imposing priors corresponding to the Rykoff et al. (2016) constraints. We adopt the mean value of $R_\lambda$ for our sample when computing mis-centering corrections, as listed in Table 1.

### 3.2. Weak Lensing Boost Factor

Galaxies that are included in the shear catalog but are not behind the galaxy clusters of our sample will not be lensed and will therefore dilute the inferred shear. Since clusters contain
many galaxies, the odds of such an occurrence increases toward the cluster center, resulting in systematic underestimation of the true $\Delta\Sigma$ profile. One typically calculates a boost factor to correct for this systematic (Sheldon et al. 2004). Our boost factor model is derived the same way as in T. McClintock et al. (2018, in preparation). We calculate the lensing-weighted average $p_{\text{phot}}(z|R_{\text{clust}})$ of source galaxies as a function of clustercentric radius and compare it with the corresponding reference $p_{\text{phot}}(z|\text{field})$ of field galaxies. The excess probability represents the member contamination in the source catalog. We then decompose the $p_{\text{phot}}(z|R_{\text{clust}})$ into two components: the reference distribution of field galaxies $p_{\text{phot}}(z|\text{field})$ and a Gaussian $p_{\text{phot}}(z|\text{Gauss})$ for the cluster member component. The decomposition is done jointly for all radial scales, such that consistency is enforced for the position and width of the Gaussian $p_{\text{phot}}(z|\text{Gauss})$ components, and only the mixing amplitude at each scale, $A(R)$, is allowed to vary. The $A(R)$ is then related to the traditional boost factor, $B(R)$, via $B(R) = 1/(1 - A(R))$. The observed lensing signal, $\Delta\Sigma_{\text{measured}}(R)$, is then related to $B(R)$ via

$$ \Delta\Sigma_{\text{measured}}(R) = \Delta\Sigma(R)/B(R). \quad (15) $$

At the minimum lensing scale considered in this work, $R = 0.2 h^{-1}$ Mpc, we find $B(R = 0.2 h^{-1}$ Mpc) – 1 ~ 0.26; for scales greater than $1 h^{-1}$ Mpc, we find $B(R) – 1 < 0.05$.

# 4. Measurement and Analysis

## 4.1. Galaxy Density Profile, $\Sigma_g$

We first measure the distribution of galaxies around the REDMAPPER clusters. The galaxy distribution is expected to roughly trace the matter distribution and has higher signal-to-noise than the weak lensing measurements in our data. The density profile of galaxies around a cluster is directly related to the galaxy cluster correlation function, $w(R)$, where $R$ is the projected comoving distance to the cluster center. We work in comoving distances so that $R = (1 + z)R_{\text{phys}}$, where $R_{\text{phys}}$ is the projected physical distance.

As shown in DK14, after averaging over the distribution of accretion rates at fixed halo mass, the location of the splashback feature is expected to scale with physical $R_{200m}$. Since the physical $R_{200m}$ is proportional to $(1 + z)^{-1}$ for fixed $M_{200m}$, measuring clustercentric radii in comoving units implicitly accounts for the redshift dependence of physical $R_{200m}$ when stacking clusters of fixed $\lambda$ at different redshifts. However, such scaling does not account for potential systematic evolution of the mean halo accretion rate with redshift.

Because our fiducial cluster sample includes a broad range of cluster richnesses ($20 < \lambda < 100$), there will be some smearing of the stacked signal due to variation in $r_{\text{sp}}$ across the bin. In Appendix A, we investigate the improvement in signal-to-noise when we approximately scale each measurement in the radial direction by the expected $R_{200m}$. Because the improvement in signal-to-noise is modest and scatter in the mass-richness relation will complicate the relationship between the scaled measurements and the simulations, our fiducial analysis does not employ this scaling.

To measure the mean-subtracted galaxy density around the REDMAPPER clusters, we divide the clusters into redshift bins of $\Delta z = 0.05$ and measure the mean cluster–galaxy angular correlation function for each bin $i$, $w(\theta, z_i)$, using the

Landy–Szalay estimator (Landy & Szalay 1993). The angular correlation function $w(\theta, z_i)$ is then converted into $w(R, z_i)$, where $R$ is the projected comoving distance for the angular separation $\theta$ at $z_i$, the center of the redshift bin. The measurements for all the redshift bins are combined by weighting $w(R, z_i)$ with the number of cluster random-galaxy random pairs $P_{\text{ran}}^i$ in each bin. The $P_{\text{ran}}^i$ values are calculated by counting the cluster–galaxy pairs in each angular bin using the cluster and galaxy random catalogs and then normalizing the pair counts by the number of clusters (galaxies) over the number of cluster randoms (galaxy randoms). That is, we calculate

$$ w(R) = \frac{\sum_{i} w(R, z_i) P_{\text{ran}}^i}{\sum_{i} P_{\text{ran}}^i}. \quad (16) $$

Finally, to convert the measured correlation function into the mean-subtracted density profile, $\Sigma_g(R)$, we multiply $w(R)$ by the mean density of galaxies around clusters, $\bar{\Sigma}_g$. This is calculated via the weighted mean of the mean galaxy density in each redshift bin $\Sigma_g^i$, where the weight is the number of clusters $N_c^i$ in that bin. We have

$$ \Sigma_g(R) = \bar{\Sigma}_g w(R), \quad (17) $$

where

$$ \bar{\Sigma}_g = \frac{\sum_{i} \Sigma_g^i N_c^i}{\sum_{i} N_c^i}. \quad (18) $$

Since we only have photo-z information for each galaxy, we cannot select galaxies that are close to each cluster in redshift very accurately. To avoid mixing galaxies with very different luminosities across the full redshift range, we create an absolute magnitude–limited sample following the approach of M16. That is, before calculating the cluster–galaxy cross-correlation in each redshift bin, we apply a luminosity cut on the galaxy sample where the absolute magnitudes are calculated assuming all galaxies are at the same redshift as the clusters. In our fiducial sample, this luminosity cut is $M^* \equiv M - 5 \log(h) < -20.23$, where the upper limit is set to be the absolute magnitude for galaxies with apparent magnitude $i = 21.5$ at $z = 0.55$.

The covariance between data points of different $R$ bins is derived using 100 jackknife samples, where the jackknife regions are derived using the “k-means” method (MacQueen 1967). The k-means method splits the data points into groups, where the groups are divided so that the spatial coordinate of all the data points in each group is closest to the mean of them. With 100 jackknife samples, each jackknife region is approximately $3.7 \times 3.7 \text{ deg}^2$, which means we can reliably measure effects up to scales ~20 Mpc at the lowest redshift of interest, $z = 0.2$.

## 4.2. Lensing Profile, $\Delta\Sigma$

The tangential shear $\gamma_t$ of a background source galaxy around a cluster is given by

$$ \gamma_t = -\gamma_1 \cos(2\phi) - \gamma_2 \sin(2\phi), \quad (19) $$

where $\phi$ is the position angle of the source galaxy with respect to the horizontal axis of a Cartesian coordinate system centered on the cluster and $\gamma_1$ and $\gamma_2$ are the two components of shear
measured with respect to the same coordinate system. For a given lens redshift $z_l$ and source redshift $z_s$, the excess surface mass density (see also Equation (7)) is related to the tangential shear according to

$$ \Delta\Sigma(z_l, z_s; R) = \langle\gamma_l(z_l, z_s; R)\rangle \Sigma_{\text{crit}}(z_l, z_s), \quad (20) $$

where $\langle\gamma_l(z_s; R)\rangle$ is the mean tangential shear for all lens-source pairs at these redshifts and $\Sigma_{\text{crit}}$ is the critical surface density in comoving units defined through

$$ \Sigma_{\text{crit}}^{-1}(z_l, z_s) = 4\pi G(1 + z_l)\chi(z_l) \left[ 1 - \frac{\chi(z_l)}{\chi(z_s)} \right], \quad (21) $$

where $G$ is the gravitational constant and $\chi(z)$ is the comoving distance to redshift $z$. To combine the full source and lens redshift distributions, we follow the same approach used in T. McClintock et al. (2018, in preparation), where we measure

$$ \Delta\Sigma = \frac{\sum_{i,j} s_{i,j} \gamma_i^t}{\sum_{i,j} s_{i,j} \Sigma_{\text{crit}}^{-1} R_{\gamma_i, i} + \left(\sum_{i,j} s_{i,j} \Sigma_{\text{crit}}^{-1}\right) \cdot \langle R_s^t \rangle} \quad (22) $$

The $i$ and $j$ subscripts denote the source and lens galaxies, respectively. The weight for each source–lens pair, $s_{i,j}$, is calculated via

$$ s_{i,j} = \Sigma_{\text{crit}}^{-1}(z_l, \langle z_s \rangle)_{i,j}, \quad (23) $$

where $z_l$ and $z_s$ are the mean redshift point-estimate provided by BPZ, and $\Sigma_{\text{crit}}^t$ is the critical density calculated using a random redshift value drawn from the source $p(z)$. Here $R_\gamma$ and $R_s$ are the METACALIBRATION responses to correct for the biased estimator and the selection bias on the ellipticity measurements.

For the covariance between data points in the different $R$ bins, we use the same 100 jackknife covariance as in the galaxy measurements (Section 4.1).

### 4.3. Model Fitting

We fit the models developed above to the data using a Bayesian approach. We define a Gaussian likelihood, $\mathcal{L}$,

$$ \mathcal{L}(d|\theta) = [d - m(\theta)]^T C^{-1} [d - m(\theta)], \quad (24) $$

where $d$ is the data vector (either $\Sigma_g$ or $\Delta\Sigma$), $m(\theta)$ is the model vector (again for either $\Sigma_g$ or $\Delta\Sigma$) evaluated at parameter values $\theta$, and $C$ is the covariance matrix of the data. The free parameters of the model are $\rho_0$, $\rho_s$, $r_t$, $r_s$, $\alpha$, $\beta$, $\gamma$, and $s_e$ and the mis-centering parameters $f_{\text{mis}}$ and $\ln c_{\text{mis}}$. Additionally, when fitting the $\Delta\Sigma$ data, we also fit for $\mu$ as defined in Equation (9).

Throughout this analysis, we restrict the range of scales we use to fit the data to 0.1–10 $h^{-1}$ Mpc for galaxies and 0.2–10 $h^{-1}$ Mpc for lensing. We do not use the 0.1–0.2 $h^{-1}$ Mpc scales for lensing in order to minimize sensitivity to the boost factor. Below ~0.1 $h^{-1}$ Mpc, the galaxy density measurements are not well understood, as the crowded nature of cluster fields and the existence of the brightest cluster galaxy renders various complications in the analysis on small scales related to detection incompleteness, photometry inaccuracy, and blending (Melchior et al. 2015, 2017). We restrict our fits to data at scales less than

<caption>**Table 2**</caption>

Priors Used for the Model Fits of Galaxy Density and Weak Lensing Mass Profiles

<table><thead><tr><th>Parameter</th><th>Priors</th></tr></thead><tbody><tr><td>ρ<sub>0</sub></td><td>[0, 10] g cm<sup>-3</sup></td></tr><tr><td>ρ<sub>s</sub></td><td>[0, 10] g cm<sup>-3</sup></td></tr><tr><td>r<sub>t</sub></td><td>[0.1, 5.0] h<sup>-1</sup> Mpc</td></tr><tr><td>r<sub>s</sub></td><td>[0.1, 5.0] h<sup>-1</sup> Mpc</td></tr><tr><td>log(α)</td><td>log(0.19) ± 0.2 (log(0.19) ± 0.1)</td></tr><tr><td>log(β)</td><td>log(6.0) ± 0.2</td></tr><tr><td>log(γ)</td><td>log(4.0) ± 0.2</td></tr><tr><td>s<sub>e</sub></td><td>[1, 10]</td></tr><tr><td>f<sub>mis</sub></td><td>0.22 ± 0.11</td></tr><tr><td>ln(c<sub>mis</sub>)</td><td>-1.13 ± 0.22</td></tr></tbody></table>

**Notes.** The value in parentheses for $\log(\alpha)$ is only applied to lensing measurements. The ranges specified in brackets are for uniform priors, while for the others, we quote the mean and standard deviation of the Gaussian priors.

10 $h^{-1}$ Mpc, since the model introduced in Section 3 breaks down at ~9 $R_{\text{vir}}$, where the power law of Equation (5) is no longer a good description of the infalling term (DK14). Since the mean virial radius for the clusters in our sample is roughly 1 $h^{-1}$ Mpc, we set the upper limit to be 10 $h^{-1}$ Mpc.

We constrain the model parameters using a Markov chain Monte Carlo (MCMC) algorithm. Following M16 and B17, we use the set of priors listed in Table 2. For $\alpha$, we use a tighter prior than that used in M16 and B17, although these priors are still very wide compared to the $\alpha$ values seen in simulations (Gao et al. 2008). We have also checked that widening the $\alpha$ priors does not affect the resulting constraints on the splashback feature. For $\beta$ and $\gamma$, the priors are informative: widening the priors leads to less constraining model fits. We discuss this point later in Section 5.2 but note that these priors still allow for a large range of profiles with and without the splashback feature. As such, the data (and not the prior) are still the main drivers that determine the slope of the profile around the splashback radius. We also note that there was a typographical error in M16. The actual Gaussian priors used were $\log(\beta) = \log(6.0) \pm 0.2$ and $\log(\gamma) = \log(4.0) \pm 0.2$, which is what is used in this work. For $f_{\text{mis}}$ and $c_{\text{mis}}$, we use values estimated in Melchior et al. (2017). Otherwise, we implement the same priors as in M16. We sample the posterior of the parameters using the `emcee` code (Foreman-Mackey et al. 2013). Convergence of the MCMC is assessed using trace plots and the Geweke statistic.

For the remainder of the paper, we define the splashback radius $r_{\text{sp}}$ to be the radius at which the logarithmic derivative of the 3D density profile $\rho(r)$ is at its minimum. To facilitate comparison with previous literature, we also define $R_{\text{sp}}$ to be the location where the logarithmic derivative of the projected galaxy density profile ($\Sigma_g$) has a minimum.

### 5. Mass Profiles of REDMAPPER Clusters

In this section, we first present the measurements of the galaxy density and lensing profiles around the fiducial cluster sample in Section 5.1. In Section 5.2, we compare the galaxy and lensing measurements and discuss the implications. These measurements are then compared to the measurements from dark matter simulations in Section 5.3. We follow that by investigating the redshift and richness dependencies of the splashback feature in Section 5.4.
## 5.1. Galaxy and Lensing Profiles

The measurement of $\Sigma_g$ around REDMAPPER clusters in our fiducial sample of $20 < \lambda < 100$ and $0.2 < z < 0.55$ is shown in the top panel of Figure 2. The red curve shows the model fit of Equations (1)–(6) to the data points with the inclusion of the mis-centering prescription of Equations (11)–(14). The bottom panel of Figure 2 shows the residuals to the fit divided by the uncertainty in the measurements. The residuals are consistent with the uncertainties, indicating that the model is a good fit to the data. The gray vertical band marks the steepest slope of the $\Sigma_g$ model, which corresponds to the 2D splashback radius $R_{sp}$. The orange vertical band marks the steepest slope of the inferred 3D density profile, $\rho(r)$, which corresponds to the splashback radius $r_{sp}$. Our measurement of the splashback radius for the fiducial sample is consistent with that from M16 within $1\sigma$ measurement uncertainty, which provides a good confirmation of their results using an independent data set and analysis pipeline.<fnref n="55" />

In Figure 3, we show the weak lensing measurement of $\Delta\Sigma$ around our fiducial REDMAPPER clusters sample with $20 < \lambda < 100$ and $0.2 < z < 0.55$. The top panel shows the data and the model fit, while the bottom panel shows the residuals of the fit divided by the uncertainty of the measurements. Again, the model provides an excellent fit to the data. The corresponding $r_{sp}$ is marked by the orange band. Although the uncertainty in $r_{sp}$ here is larger than that derived from the galaxy density profile, we note that $r_{sp}$ is very well constrained (compared to, e.g., Umetsu & Diemer 2017). The high signal-to-noise of this measurement is a result of the combination of a large number of clusters and background source galaxies.

## 5.2. Splashback Feature of REDMAPPER Clusters

In Figure 4, we present the results of the model fits to the galaxy density and weak lensing measurements. Throughout the three panels, the vertical lines mark the mean $r_{sp}$ derived from the galaxy density (gray) and the weak lensing (red) profiles, whereas the horizontal bars in the middle panel indicate the uncertainties (standard deviation of the $r_{sp}$ distribution) of the two $r_{sp}$ values. The gray and red bands show the 16th–84th percentile confidence range for each profile. The top panel shows the ratio $\rho^{\text{coll}}(r)/\rho(r)$, i.e., the fraction of the total density profile that is part of the collapsed material profile. The difference in normalization between the galaxy density and lensing measurements cancels in this ratio. We find that the galaxy density and lensing measurements yield very consistent collapsed fractions, with the lensing measurements being slightly higher. One might worry that we are drawing conclusions about the collapsed material in a regime where it is completely dominated by the infalling term. The top panel of Figure 4 makes it clear that this is not the case: near $r_{sp}$, the collapsed profile term makes up 40%–50% of the total profile. Our inferences about the collapsed profile in the transition regime are therefore robust as long as the infalling material has a relatively smooth profile, which is a good assumption here. Alternatively, one can introduce an additional “limiting density” term in the denominator as in Equation (42).

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Comparison of model-fit results from galaxy density $\Sigma_g$ (gray) and weak lensing $\Delta\Sigma$ (red). Top: fraction of the density profile for the collapsed material over the total density profile. Middle: logarithmic derivative of the total density profile compared to the logarithmic derivative of an NFW profile (dashed curve). Bottom: logarithmic derivative of the profile for the collapsed material compared to the logarithmic derivative of an NFW profile. The vertical lines mark the mean $r_{sp}$ inferred from the model fits for both galaxy and lensing measurements, while the horizontal bars in the middle panel indicate the uncertainties on $r_{sp}$.</figcaption></figure>

of Diemer (2017) to avoid spurious contribution on small scales from the infalling term.

The middle panel of Figure 4 shows the logarithmic derivative of the total density profile inferred from the galaxy density and lensing measurements. The locations of the steepest slopes in the two profiles are consistent with the lensing measurement: the galaxy profile gives $r_{sp} = 1.13 \pm 0.07 \, h^{-1} \, \text{Mpc}$, and the weak lensing profile gives $r_{sp} = 1.34 \pm 0.21 \, h^{-1} \, \text{Mpc}$. The amplitudes and shapes of the logarithmic derivative profiles are quite consistent, with the galaxy density profile slightly steeper at large radii. We find the total profile of both our galaxy and lensing measurements to be steeper than an NFW profile of similar mass at $r_{sp}$ (as we discuss in more detail below). This is consistent with the expectation for a splashback feature.

An alternative is to look at the logarithmic slope of the collapsed profile, which is also the approach taken by B17. This approach includes our model for the profile of the infalling material, which is assumed to be a power law. In the bottom panel of Figure 4, we show the logarithmic slope of the collapsed profile inferred from the galaxy density and lensing measurements. We find that at $r_{sp}$, the inferred collapsed profiles from both galaxy and lensing profiles exhibit rapid

<footnote n="55">⁵⁵ The redshift ranges of the clusters are somewhat different in the two analyses, but as we show in Section 5.4, there does not appear to be significant redshift evolution in this range of redshifts.</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Top: posterior distributions of $r_{sp}$ for galaxy density (gray) and lensing (red) data. Middle: posterior distribution of the slope of the total matter profile at $r_{sp}$. Bottom: posterior distribution of the slope of the collapsed matter at $r_{sp}$. All the distributions are marginalized over all nuisance parameters.</figcaption></figure>

steepening, achieving values much steeper than the slope of an NFW profile at scales around $r_{sp}$ and beyond. This again is consistent with the picture that a splashback feature exists at the outskirts of these clusters.

The posterior distributions of $r_{sp}$ and the slope of the total and collapsed profiles in Figure 4 are shown in Figure 5. Here we clearly see that the galaxy and lensing measurements of $r_{sp}$ and the slopes of the profiles are consistent with each other, with the lensing measurements having larger uncertainties. The measured logarithmic slope of the total profile at $r_{sp}$ is -3.6 ± 0.3 and -3.5 ± 0.4 for the galaxy density and lensing profiles, respectively. The measured logarithmic slope of the collapsed profile is -5.9 ± 0.7 and -5.3 ± 0.9 for the galaxy density and lensing profiles, respectively. These measured slopes can be compared to the expectation for an NFW profile. For the NFW profile predicted by the mass-richness relation of Melchior et al. (2017), the logarithmic slope at $r_{sp}$ is ~-2.7, while the maximum possible slope is -3. The slope of the total profile is therefore steeper than the NFW at roughly 3.0σ for the galaxy density measurements and 2.0σ for the lensing measurements. However, the NFW profile does not fully capture the contribution from infalling material near the cluster, which generically makes the profile less steep at $r_{sp}$. Comparing the slope of only the collapsed component to that of the NFW profile, we find that it is steeper than the NFW by 4.6σ for the galaxy density profile and 2.9σ for the lensing

profile. The values of $r_{sp}$ derived from the MCMC, as well as the model parameters, are listed in Table 3.

As discussed in Section 4.3, the parameters $\beta$ and $\gamma$ are important for determining the behavior of the profile around the splashback feature. These parameters are degenerate, and the priors that we place on them are informative. To test how relaxing these priors would affect the splashback measurement from lensing, we completely relax the $\gamma$ priors and examine the constraints on the slope of the profiles. We find the slope of the total (collapsed) profile at $r_{sp}$ to be -3.7 ± 0.6 (-6.2 ± 2.0) for the lensing measurement. This corresponds to a roughly 1.6$\sigma$ (1.8$\sigma$) steeper profile compared to the NFW profile at $r_{sp}$. We also perform an additional check to see whether the priors are wide enough to span a range of profiles with and without a splashback feature “detection.” That is, we check that the priors are not driving us to falsely detect a splashback-like steepening. To check this, we sample the priors of $\alpha, \beta, \gamma$, and $r_t$ (the most relevant parameters for the splashback feature), generate model profiles, and measure the slope of the profile at $r_{sp}$. The resulting slope distribution is shown in Figure 6. Noting that the minimum logarithmic slope achieved by an NFW profile is -3, we see that the priors allow profiles with slopes both shallower and steeper than the NFW.

### 5.3. Comparison to N-body Simulations

We now compare the galaxy density and lensing measurements around REDMAPPER clusters to similar measurements made using dark matter-only N-body simulations. For this purpose, we use the MultiDark Planck 2 simulation from the CosmoSim database (Prada et al. 2012; Riebe et al. 2013; also see www.cosmosim.org).

Using the ROCKSTAR (Behroozi et al. 2013) halo catalogs made available by CosmoSim, we identify a set of halos that is matched to the REDMAPPER cluster catalog used in this work. We match the REDMAPPER clusters to the simulated dark matter halos on the basis of halo mass and redshift. Using the best-fit mass-richness relation from Melchior et al. (2017), we calculate the mean $M_{200m}$ halo mass of our fiducial sample to be $2.5 \times 10^{14} M_{\odot}$. We then determine in the simulations a mass threshold, $M_{min}$, such that the mean $M_{200m}$ mass of simulated halos between $M_{min}$ and $10^{15} h^{-1} M_{\odot}$ is equal to $\bar{M}$. We find $M_{min} = 1.0 \times 10^{14} M_{\odot}$. The upper mass limit here has little impact on our results but ensures that a very small number of extremely massive halos is not skewing our predictions. A total of 11,745 clusters are used in the simulations, and all are at a single redshift $z = 0.52$.<fnref n="56" /> We note that, as we are employing the mass-richness relation from a different sample (DES SV) and there is approximately a 10% scatter on the mass-richness relation, we can imagine that the masses in Table 1 could be over/underestimated. If the mass estimates were off by $1\sigma$, the inferred $r_{sp}$ would move by ~3%.

The matching of our galaxy sample to objects in the MultiDark simulations is more complicated, since the simulations only contain dark matter. We use both dark matter subhalos and dark matter particles to perform the comparison with the galaxy and lensing measurements. The connection between galaxies and subhalos depends on a combination of environmental parameters (Reddick et al. 2013;

<footnote n="56">⁵⁶ The mean redshift of our cluster sample is $z = 0.41$. We have chosen the closest redshift slice at $z = 0.52$ in our simulations for an approximate comparison with the data. However, we do not expect the subhalo profiles to vary significantly over this redshift range.</footnote>
<table><caption>**Table 3**
Model Parameters for $\Sigma_g$ and $\Delta\Sigma$</caption><thead><tr><th></th><th>Sample</th><th>log(r<sub>s</sub>)</th><th>log(r<sub>t</sub>)</th><th>log(α)</th><th>log(β)</th><th>log(γ)</th><th>s<sub>e</sub></th><th>r<sub>sp</sub></th><th>r<sub>sp</sub>/r<sub>200m</sub></th></tr></thead><tbody><tr><td rowspan="6">Σ<sub>g</sub></td><td>Fiducial</td><td>-0.82 ± 0.10</td><td>0.03 ± 0.05</td><td>-0.83 ± 0.12</td><td>0.92 ± 0.14</td><td>0.70 ± 0.15</td><td>1.57 ± 0.07</td><td>1.13 ± 0.07</td><td>0.82 ± 0.05</td></tr><tr><td>Low-z</td><td>-0.65 ± 0.12</td><td>-0.01 ± 0.05</td><td>-0.84 ± 0.14</td><td>0.92 ± 0.17</td><td>0.64 ± 0.16</td><td>1.61 ± 0.09</td><td>1.07 ± 0.09</td><td>0.73 ± 0.06</td></tr><tr><td>Mid-z</td><td>-0.78 ± 0.09</td><td>0.07 ± 0.09</td><td>-0.72 ± 0.14</td><td>0.82 ± 0.16</td><td>0.66 ± 0.17</td><td>1.62 ± 0.09</td><td>1.12 ± 0.14</td><td>0.85 ± 0.11</td></tr><tr><td>High-z</td><td>-0.66 ± 0.12</td><td>0.04 ± 0.11</td><td>-0.73 ± 0.15</td><td>0.86 ± 0.18</td><td>0.60 ± 0.20</td><td>1.46 ± 0.11</td><td>1.14 ± 0.15</td><td>0.96 ± 0.13</td></tr><tr><td>Low-λ</td><td>-0.75 ± 0.09</td><td>0.01 ± 0.06</td><td>-0.73 ± 0.15</td><td>0.91 ± 0.16</td><td>0.73 ± 0.17</td><td>1.50 ± 0.08</td><td>1.05 ± 0.09</td><td>0.85 ± 0.07</td></tr><tr><td>High-λ</td><td>-0.81 ± 0.10</td><td>0.10 ± 0.08</td><td>-0.83 ± 0.12</td><td>0.79 ± 0.14</td><td>0.62 ± 0.17</td><td>1.53 ± 0.08</td><td>1.27 ± 0.14</td><td>0.83 ± 0.09</td></tr><tr><td>ΔΣ</td><td>Fiducial</td><td>-0.62 ± 0.15</td><td>0.15 ± 0.11</td><td>-0.71 ± 0.10</td><td>0.75 ± 0.17</td><td>0.72 ± 0.18</td><td>1.23 ± 0.24</td><td>1.34 ± 0.21</td><td>0.97 ± 0.15</td></tr></tbody></table>

Note. Here $r_s$, $r_t$, and $r_{sp}$ are in units of $h^{-1}$ Mpc, and $r_{200m}$ is calculated using the mean mass and redshift for each cluster sample as listed in Table 1.

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Distribution of the logarithmic slope at $r_{sp}$ when sampling over the prior distribution of several model parameters ($\alpha$, $\beta$, $\gamma$, $r_s$). The black dashed line indicates a rough indicator of the slope for an NFW profile.</figcaption></figure>

Chaves-Montero et al. 2016). Furthermore, subhalo density profiles around massive halos are known to be flatter on small scales compared to galaxy density profiles tidally stripped of mass near the cluster center (plus, there can be numerical artifacts on these scales in the simulations). Galaxies, on the other hand, tend to live at the centers of their dark matter halos and are therefore less likely to be tidally stripped (e.g., Nagai & Kravtsov 2005; Budzynski et al. 2012). These effects result in differences between the matter, subhalo, and galaxy profiles but mainly in the inner regions of the parent halo, well below the splashback radius.

To construct a subhalo sample corresponding to our galaxy sample, we select ROCKSTAR (Behroozi et al. 2013)-identified subhalos using $v_p$, the largest circular velocity attained by the subhalo over its history, which corresponds roughly to a mass cut on the subhalos at the time of accretion (Reddick et al. 2013). Unlike selection based on subhalo mass, selection on $v_p$ is expected to lead to a sample that more closely approximates real galaxies, since it is unaffected by tidal stripping of mass from the subhalo. The subhalo density profiles, $\Sigma_{sub}$, for the simulated subhalos around dark matter halos of $M_{min} < M < 10^{15} h^{-1} M_{\odot}$ in the MultiDark simulations are shown in Figure 7, overlaid with the data $\Sigma_g$ measurements. At large scales, the amplitude of the $\Sigma_g$ curve scales with the abundance and can be compared with the galaxy density profile to find the approximate subhalo mass corresponding to a given galaxy sample. We consider three $v_p^{min}$ values: 135, 178, and 280 km s$^{-1}$. In the simulations, there are, on average, 28, 17, and 7 subhalos per cluster within $1.5 h^{-1}$ Mpc of the halo center for the three $v_p$ cuts, respectively. As seen in Figure 7, the galaxy sample in our data lies between the two subhalo samples $v_p^{min} = 135$ and

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Subhalo density profiles measured in simulations around halos with mass similar to that of our fiducial cluster sample. Different colors correspond to different choices of subhalo $v_p$. The data points are the galaxy profile measured with our fiducial sample, which lie between the two lower-mass subhalo samples. The lighter curves indicate the range excluded from the model fits described in Section 5.3.</figcaption></figure>

$178$ km s$^{-1}$. That is, we can identify our galaxy sample with subhalos that are less massive than the $v_p^{min} = 178$ km s$^{-1}$ sample and more massive than the $v_p^{min} = 135$ km s$^{-1}$ sample. Note that we do not employ a rigorous abundance-matching procedure similar to M16. As a result, the amplitudes of our data points in Figure 7 and later in Figure 10 do not match the subhalo profiles exactly.

We fit the model described in Section 3 to the subhalo profiles, excluding scales below $0.5 h^{-1}$ Mpc to minimize bias induced by the tidal-stripping effect on small scales mentioned above. The model describes the subhalo profiles well after excluding the small scales. In the top panel of Figure 8, we compare the logarithmic derivative of the model profile from our fiducial sample and from the two lower-mass subhalo bins (since these bins bracket our galaxy sample). The inferred $r_{sp}$ and uncertainty for each of the curves shown in the top panel of Figure 8 are marked by horizontal bars on the top of the panel. As seen in the figure, the two lowest-mass subhalo bins have essentially the same $r_{sp}$, indicating that these subhalos are sufficiently small that they are not affected by dynamical friction. Since these two subhalo samples have masses that bracket that of our galaxy sample, we conclude that our measurements of $r_{sp}$ from the galaxy density profile are not affected by dynamical friction. We will present a more thorough analysis of dynamical friction in Section 6.
<figure><img src="image_10.png" /></figure>

**Figure 8.** Comparison of measurements from dark matter simulations and data. Top: log derivatives of the model fit to the galaxy profiles in data and the subhalo profiles in simulations. The horizontal bars in each panel indicate the inferred location and uncertainty of $r_{sp}$. Note that $r_{sp}$ in the data is smaller than in the subhalo cases that are best matched to our galaxies. The lighter sections of the green and red curves indicate the regime where we expect differences between the data and simulations, as we do not fit the subhalo profiles on small scales. Bottom: same as top panel but comparing the slope of the profile of the dark matter particles with the lensing measurements.

The $r_{sp}$ inferred from our galaxy density profile (1.13 ± 0.07 $h^{-1}$ Mpc) is significantly smaller than the corresponding subhalo measurements (1.46 ± 0.05 $h^{-1}$ Mpc for the $v_p^{\min} = 178$ km s$^{-1}$ subhalo sample), as seen in Figure 8. However, the steepest slope inferred from the simulations and data appears to be consistent, suggesting that we are seeing a level of steepening in the galaxy profile that is consistent with the splashback feature in simulations. The overall shape of the galaxy profile in the data differs somewhat from that of subhalos in the simulations, where the small-scale differences have been addressed above. These findings are consistent with those of M16.

In the bottom panel of Figure 8, we compare the lensing measurements with the dark matter particles. When fitting to the particle measurements, we do not include the effects of miscentering. We find that the particles give consistent $r_{sp}$ values as the two lower-mass subhalo samples in the middle panel and are larger than the lensing measurements by about 18%. We note that the seemingly better agreement between the measurements and the simulations (about $1\sigma$) is mainly driven by the fact that the lensing measurements have larger uncertainties. The slope of the lensing profile at large radii is shallower than the simulation particles; the same trend is seen in the galaxy versus subhalo profiles. We have not investigated possible sources of this ≈2$\sigma$ discrepancy.

## 5.4. Richness and Redshift Dependences of $r_{sp}$

We now consider the richness dependence of the splashback feature. According to simulation tests in DK14 and A14, one would expect the splashback feature to be shallower and appear at smaller scales for lower-mass (or richness) clusters. We measure the richness dependence of the splashback location by

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

**Figure 9.** Top: logarithmic derivative of the model fits to the $\Sigma_g$ measurements with different richnesses. Bottom: similar to the top panel, but for different redshift bins. The horizontal bars in each panel indicate the inferred location and uncertainty of $r_{sp}$ in the different subsamples.

dividing the fiducial cluster sample into two richness subsamples: $-20 < \lambda < 28$ and $28 < \lambda < 100$. The bins are chosen so that the number of clusters is approximately equal in both bins. The mean richness in the two bins is 23.3 and 41.1, respectively. In the top panel of Figure 9, we show the log derivatives of the model fits to the galaxy density profiles of these two subsamples. We find that $r_{sp}$ is 1.05 ± 0.09 $h^{-1}$ Mpc and 1.27 ± 0.14 $h^{-1}$ Mpc for the low- and high-richness samples, respectively. The dependence of the mean $r_{sp}$ on the mean $\lambda$ is roughly $r_{sp} \propto \lambda^{0.33 \pm 0.24}$, which is consistent with the expectation<fnref n="57" /> from the slope of the mass–richness relation of REDMAPPER clusters measured in Melchior et al. (2017), $r_{sp} \propto \lambda^{0.37}$. We note, however, that the detailed shapes of the logarithmic derivatives measured from the data exhibit some puzzling differences from simulations. In particular, we find that the high-richness cluster sample has a shallower splashback feature than the low-richness cluster sample. In the simulations of Diemer & Kravtsov (2014), on the other hand, higher-mass halos tend to have sharper splashback features.

In principle, our measurement of the richness dependence of the splashback radius could be impacted by dynamical friction. As discussed in Section 1, dynamical friction will result in a decrease in the observed splashback radius measured via the galaxy density profile. This effect is expected to be weaker for larger host halos, which could result in an increase in the observed scaling of the splashback radius with mass relative to the expectation from particles in simulations (the particle profile is not impacted by dynamical friction). However, as we show in Section 6, for our fiducial galaxy sample, dynamical friction does not appear to have a significant impact on the inferred splashback radius. Consequently, our measurement of the richness dependence of the splashback radius can be

<footnote n="57">⁵⁷ Since $r_{sp} \propto R_{200m}$, we expect $r_{sp} \propto M_{200m}^{1/3}$. Melchior et al. (2017) found $M_{200m} \propto \lambda^{1.12}$, suggesting $r_{sp} \propto \lambda^{0.37}$.</footnote>
compared directly to the expectation from particles in
simulations.

We next consider the redshift dependence of the splashback feature. A14 looked at the redshift dependence of the splashback feature in simulations, finding that for a given accretion rate, $r_{sp}$ becomes larger at higher redshift, which results from a simple scaling with the background cosmology (specifically $\Omega_m$). When averaged over a distribution of accretion rates, however, DK14 found that the results are consistent with no redshift evolution. We test this by performing the same $\Sigma_g$ measurement in three redshift bins: $0.2 < z < 0.4$, $0.4 < z < 0.55$, and $0.55 < z < 0.75$. The lowest-redshift bin is similar to that used in M16, whereas the highest-redshift sample is not strictly volume-limited. In the bottom panel of Figure 9, we show the log derivative of the model fit to the measurements for the three redshift bins, with the inferred $r_{sp}$ marked on the plot and listed in Table 3. We find no evidence of redshift evolution of $r_{sp}$ over this redshift range. Given that we do not select the clusters in accretion rate, our finding of no redshift evolution is consistent with that found in DK14. One might worry that the mass-richness relation also evolves with redshift, which could complicate the comparison. However, in our sample, we do not find a significant evolution of mass over the three redshift samples (see Table 1), which means we indeed do not see a redshift evolution of $r_{sp}$ for fixed halo mass.

## 6. Effect of Dynamical Friction

As discussed in Section 1, measuring the splashback radius provides an avenue for detecting the effects of dynamical friction in galaxy clusters. The rate of deceleration due to dynamical friction for a subhalo traveling through a cluster is proportional to the mass of the subhalo. Consequently, more massive (brighter) galaxies are expected to splashback at smaller radii. We first test this expectation in simulations by looking at the log derivative of the model fits to the three $\Sigma_{sub}$ curves in Figure 7. The corresponding log-derivative profiles are shown in the middle panel of Figure 10, together with the inferred $r_{sp}$. It is clear that subhalos with $v_p > 280$ km s$^{-1}$ have a significantly smaller splashback radius than lower-mass subhalos, the expected consequence of dynamical friction. For $v_p > 280$ km s$^{-1}$ subhalos, we find $r_{sp} = 1.21 h^{-1}$ Mpc, while for the other two subhalo samples, we find $r_{sp} = 1.47 h^{-1}$ Mpc. This ~20% difference is consistent with that found in Diemer et al. (2017).

Since we cannot directly measure the masses of the galaxies in our sample, we divide the galaxies based on luminosity, which correlates with mass. We define three luminosity bins from our galaxy sample ($M_* < -19.4, -20.4$, and $-21.4$) and measure the resultant $\Sigma_g$ profiles around the low-$z$ cluster sample ($20 < \lambda < 100$, $0.2 < z < 0.4$), as shown in the top panel of Figure 10. We use the low-$z$ sample so that we can lower the luminosity cut on the galaxies and have higher signal-to-noise measurements. Overlaying the same subhalo profiles from the dark matter simulation as in Figure 7, we find that the brightest galaxy bin ($M_* < -21.4$) roughly corresponds to the most massive subhalo bin ($v_p > 280$ km s$^{-1}$), which is also the sample that has showed signs of dynamical friction. The two fainter galaxy bins roughly correspond to the two lower-mass subhalo samples. We fit all three galaxy measurements to the same model used in Section 5.1 and show the log-derivative profile of the models in the bottom panel of

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Effects of changing the galaxy luminosity cut on the inferred $r_{sp}$ around $20 < \lambda < 100$, $0.2 < z < 0.4$ clusters. Top: measurement of $\Sigma_g$ profiles for the three luminosity bins (data points) and the three subhalo profiles in Figure 7 (solid lines). The lighter curves indicate the range excluded from the model fits. The subhalo samples are not abundance-matched to the galaxies; therefore, we do not expect the amplitudes of the data points to agree with the solid lines. Middle: log derivatives of model fits to the subhalo density profiles measured in simulations from the top panel. The faded section of the curves indicates the regime where we expect differences between the data and simulations, as we do not fit the subhalo profiles on small scales. Bottom: log derivative of the model fits to the three galaxy density profiles. The horizontal bars indicate the inferred location and uncertainty of the 3D $r_{sp}$ for each galaxy sample.</figcaption></figure>

**Figure 10.** The galaxies show similar behaviors to what was observed with the subhalos in the dark matter simulations: the two fainter galaxy bins have consistent $r_{sp}$ measurements, while the brightest galaxy bin has a slightly smaller $r_{sp}$. However, the difference between the brightest galaxy bin and the other two bins is smaller than what is expected from the simulations and well within the measurement uncertainties. Furthermore, as we show in Appendix B, this measurement is
sensitive to the choice of the REDMAPPER parameter $R_0$ (see Section 7). This test does, on the other hand, confirm that our fiducial galaxy sample used in Section 5 is not affected by dynamical friction.

Comparing in more detail the bottom two panels of Figure 10, we also find other qualitative differences in the profiles: the most massive galaxy sample shows a shallower log derivative compared to the other two galaxy bins, which is in the opposite direction of what is expected from the subhalo simulations. To further investigate these subtle differences and systematics effects would require more realistic simulations that capture the baryonic physics on small scales. We defer this study to future work.

One can imagine further increasing the effect of the dynamical friction by going to lower-mass clusters, an approach taken by Adhikari et al. (2016). This is because one expects the effect of dynamical friction to be larger for smaller host halos (the effect of dynamical friction scales with $M_{sub}/M_{host}$, where $M_{sub}$ is the subhalo mass and $M_{host}$ is the host halo mass). However, we note that the mass estimates for REDMAPPER clusters below $\lambda = 20$ are less reliable, as shown in Melchior et al. (2017). We therefore do not perform further measurements using the low-richness clusters. We also tested that our conclusion of this analysis does not change when using nonoverlapping magnitude bins, which could, in principle, enhance the effect of dynamical friction.

## 7. Potential Biases Due to the REDMAPPER Algorithm

One potential concern for splashback measurements relying on REDMAPPER clusters was pointed out in recent work by Busch & White (2017). REDMAPPER identifies clusters based on overdensities of red galaxies on the sky. Selection effects in REDMAPPER could therefore result in changes to the measured galaxy density profile relative to the true galaxy density profile around REDMAPPER clusters, which could potentially result in biases to splashback measurements that use the galaxy density profile around these clusters. We review and investigate this issue below.

In the M16 measurement of splashback with SDSS, REDMAPPER-identified clusters were split into two subsamples based on $\langle R_{mem} \rangle$, defined as the average clustercentric distance of the cluster members weighted by membership probability. These two samples were found to have significantly different $r_{sp}$, as well as large-scale clustering amplitudes. Zu et al. (2016) and Busch & White (2017) later pointed out that $\langle R_{mem} \rangle$ can be strongly affected by projection effects due to the way REDMAPPER assigns members to clusters. That is, clusters that live in dense environments are likely to have a large number of spurious members from line-of-sight projections that have low membership probability but contribute to large $\langle R_{mem} \rangle$ values. These clusters will have a higher large-scale clustering amplitude as a result of their association with projected structures along the line of sight. Selecting on $\langle R_{mem} \rangle$ can therefore result in spurious assembly bias signals (Miyatake et al. 2016; Zu et al. 2016). Given the sensitivity of $\langle R_{mem} \rangle$ to projection effects, we have not employed this quantity in our analysis.

Using the Millennium simulation (Springel et al. 2005), Busch & White (2017) also argued that the aperture radius, $R_c$, used by REDMAPPER to define cluster richness could impact the cluster density profile and therefore the splashback feature. REDMAPPER computes the cluster richness as a weighted sum

over galaxies within $R_c$ of the assumed cluster center, where

$$R_c(\lambda) = R_0(\lambda/\lambda_0)^\beta. \qquad (25)$$

Values of $R_0 = 1.0 h^{-1}$ Mpc, $\beta = 0.2$, and $\lambda_0 = 100$ were chosen to minimize scatter in the mass-richness relation (Rozo et al. 2009). Busch & White (2017) considered the effects of changing $R_0$ on the results of their simulated REDMAPPER measurement.<fnref n="58" /> For catalogs generated with $R_0 = 0.67, 1$, and $1.5 h^{-1}$ Mpc, Busch & White (2017) found that the inferred $r_{sp}$ and profile shapes were altered. Note that, as mentioned above, these analyses were carried out using the Millennium simulation (Springel et al. 2005) and a simplified procedure that approximates the REDMAPPER algorithm. Consequently, the exact quantitative effect on $r_{sp}$ from the $R_c$ selection may not be directly applicable to our data measurements.

To test the impact of $R_c$ on our splashback measurements, we rerun REDMAPPER, setting $R_0 = 0.75$ and $1.25 h^{-1}$ Mpc. The resultant cluster catalogs will have a new richness estimate $\lambda'$ for each cluster. We rank in descending order the old cluster catalog by $\lambda$ and the new cluster catalog by $\lambda'$, then select the clusters in the new catalog that have the same ranking as the fiducial sample in our original cluster catalog with $20 < \lambda' < 100$. We find that $16.5 < \lambda' < 75.3$ ($22.8 < \lambda' < 117.3$) gives roughly the same number of clusters with the same ranking for $R_0 = 0.75$ ($1.25 h^{-1}$ Mpc). We measure $\Sigma_g$ and $\Delta\Sigma$ for these two new cluster catalogs and fit them to our model. We have checked that the amplitudes of the $\Delta\Sigma$ measurements are nearly identical for the different $R_0$ settings, suggesting that the mean mass of the samples did not change significantly when we changed $R_0$. We note that the choices of $R_0$ here are rather extreme, and the REDMAPPER code is not well tested at these $R_0$ values. For instance, we expect the scatter in the mass-richness relation to be much larger at these extreme $R_0$ values, which could have an effect on the resulting stacked profiles. As a result, the tests below should be treated as bounds for the potential systematic effects introduced by the $R_c$ settings.

The top panel of Figure 11 shows the log derivative of the model fits to $\Sigma_g$ for these two cases together with the fiducial setting of $R_0 = 1 h^{-1}$ Mpc. We find that the profiles do indeed change as a function of $R_0$, similar to what was seen in Busch & White (2017). The quantitative change in our measurements is, however, smaller than that seen in Busch & White (2017), likely because Busch & White (2017) employed a simplified REDMAPPER-like cluster finder.

For $\Delta\Sigma$, we show in the bottom panel of Figure 11 the resulting log derivative of the model fits for the two alternative $R_0$ values. Similar to what is seen in the $\Sigma_g$, the location of $r_{sp}$ moves outward as the $R_0$ value increases. However, for each of the $R_0$ settings, the lensing-inferred $r_{sp}$ remains consistent with the galaxy measurements at better than $1\sigma$. Also, compared to the NFW profile, shown by the dashed black curve, the slope at $r_{sp}$ for the lensing remains steeper than the NFW by about $1\sigma$. We note, however, that the lensing-inferred $r_{sp}$ appears to be more robust to the change in $R_0$ than the galaxies.

The variation of $r_{sp}$ with $R_0$ is not necessarily indicative of a systematic error in the measured mass profile or the inferred values of $r_{sp}$. Instead, it could suggest another source of

<footnote n="58">⁵⁸ Since increasing $R_0$ necessarily means that clusters will have larger richness, Busch & White (2017) simultaneously varied $R_0$ and $\lambda_0$.</footnote>
the less massive galaxies, which would result in a smaller $r_{sp}$. We measure the profiles of galaxies in different luminosity bins around clusters and find that the highest-luminosity galaxies indeed exhibit a slightly smaller $r_{sp}$, a behavior that matches the corresponding subhalo profiles of the dark matter simulations. However, the difference is smaller than expected from simulations and within measurement uncertainties. We also tested that this measurement is sensitive to the change in $R_0$ mentioned above.

Looking toward the next DES data set, which covers the full footprint of 5000 deg², we can expect significant improvement in the statistical uncertainties in both the galaxy and the lensing measurements, as well as the redshift coverage. However, interpreting the subtle systematic effects in the cluster-finding algorithm and measurement process will be the crucial next step for a deeper understanding of the connection between the true splashback feature and the observed cluster profiles. One important step is to develop more realistic simulations that can reproduce the observables. In parallel, exploring the splashback feature for cluster samples selected in other wavelengths (in particular, SZ- and X-ray-selected samples) would be a good test for potential systematics in the optical cluster finder. A lensing mass-selected cluster sample must await surveys that are deep enough to provide high-significance detections of individual clusters. Improvements in simulations and cluster selection will enable us to control for systematic effects and pursue the effects of standard and new physics associated with the splashback feature.

C.C. and A.K. were supported in part by the Kavli Institute for Cosmological Physics at the University of Chicago through grant NSF PHY-1125897 and an endowment from the Kavli Foundation and its founder, Fred Kavli. E.B. and B.J. are partially supported by US Department of Energy grant DE-SC0007901. T.N.V. was supported by the SFB-Transregio 33 The Dark Universe by the Deutsche Forschungsgemeinschaft (DFG) and the DFG Cluster of Excellence “Origin and Structure of the Universe.” The weak lensing boost factors were calculated and calibrated using the computing facilities of the Computational Center for Particle and Astrophysics (C2PAP). D.R. is supported by a NASA Postdoctoral Program Senior Fellowship at the NASA Ames Research Center, administered by the Universities Space Research Association under contract with NASA.

The CosmoSim database used in this paper is a service by the Leibniz-Institute for Astrophysics Potsdam (AIP). The MultiDark database was developed in cooperation with the Spanish MultiDark Consolider Project CSD2009-00064. The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) and the Partnership for Advanced Supercomputing in Europe (PRACE; www.prace-ri.eu) for funding the MultiDark simulation project by providing computing time on the GCS Supercomputer SuperMUC at the Leibniz Supercomputing Centre (LRZ; www.lrz.de). The Bolshoi simulations have been performed within the Bolshoi project of the University of California High-Performance AstroComputing Center (UC-HiPACC) and were run at the NASA Ames Research Center.

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 InterAmerican 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 Numbers AST-1138766 and AST-1536171. The DES participants from Spanish institutions are partially supported by MINECO under grants AYA2015-71825, ESP2015-66861, FPA2015-68048, SEV-2016-0588, SEV-2016-0597, and MDM-2015-0509, some of which include ERDF funds from the European Union. The IFAE 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 240672, 291329, and 306478. We acknowledge support from the Australian Research Council Centre of Excellence for All-sky Astrophysics (CAASTRO) through project number CE110001020.

This manuscript has been authored by the 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, worldwide 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.
<figure><img src="image_13.png" /><figcaption>**Figure 12.** Log derivative for fiducial measurements with and without the λ-scaling that accounts for the range of λ inside the bin. The inferred $r_{sp}$ values are marked as horizontal bars at the top of the figure. The x-axis for the scaled case is $r/R_{200m}\bar{R}_{200m}$, where $\bar{R}_{200m}$ is the $R_{200m}$ at the mean richness.</figcaption></figure>

## Appendix A
Scaling Measurements with Richness λ

As shown in DK14, the splashback radius scales with physical $R_{200m}$ and the accretion rate. We discussed in Section 4.1 that our measurements are performed in comoving distances that take into account the redshift evolution of $R_{200m}$. Here we investigate the improvement in the measurements if we were to also take into account the richness dependence of $R_{200m}$ within the sample. We do not attempt to correct for the accretion rate dependence in the same way, as the estimation of accretion rates for clusters is nontrivial.

We subdivide our fiducial cluster sample (with selection $20 < \lambda < 100$ and $0.2 < z < 0.55$) into 10 logarithmic $\lambda$ bins and repeat the measurement in Section 4.1. The measured distances of galaxies from the cluster centers in each $\lambda$ bin are then scaled by $(\bar{\lambda}_i/\bar{\lambda}_{\text{full}})^{F/3}$, where $\bar{\lambda}_i$ is the mean richness in the bin, $\bar{\lambda}_{\text{full}}$ is the mean richness in the full sample, and $F = 1.12$ is the exponent of the mass–richness relation derived in Melchior et al. (2017).

Figure 12 shows the log derivative of the model fit to the fiducial measurement (Figure 2) with and without the λ-scaling. We find that after taking into account the λ dependence of $r_{sp}$, the splashback feature does not change significantly. This is somewhat counterintuitive given the results in Figure 9. Taking a closer look at the measurements, we find that the improvement in the λ-scaling is mostly washed out by the slightly increased error bars in the measurements, which is likely a result of the large scatter in the mass–richness relation.

## Appendix B
Additional Tests of the Impact of $R_0$ on the Splashback Feature

As discussed in Section 7, we have seen that the choice of $R_0$ in REDMAPPER affects the inferred splashback radius, likely a result of selection effects. In this appendix, we carry out a few more tests to see the impact of $R_0$ on other measurements in this paper.

### B.1. Effect of $R_0$ on $r_{sp}$ in High-richness Clusters

Recently, using a combination of weak lensing and abundance measurements, Murata et al. (2018) found that the REDMAPPER mass–richness relation exhibits unexpectedly large scatter at low richness. A nonnegligible fraction of the clusters with richness ~20 comes from halos of mass

<figure><img src="image_14.png" /><figcaption>**Figure 13.** Same as the top panel of Figure 11, but for clusters of richness $28 < \lambda < 100$.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 14.** Log derivative for galaxy profiles of different luminosity, with cluster samples of different $R_0$ settings: the top (bottom) panel shows the same measurements as the bottom (top) panel of Figure 10, but with $R_0 = 0.75$ ($R_0 = 1.25$) $h^{-1}$ Mpc.</figcaption></figure>

$\approx 10^{13} M_{\odot}$. One of the hypotheses in Murata et al. (2018) is that the low-richness clusters are affected by projection effects and thus less reliable. Inspired by this finding, we perform the test on the high-richness clusters in Section 5.4 to see whether they are more or less sensitive to the choice of $R_0$. The results are shown in Figure 13. We find that the high-richness clusters give lower signal-to-noise results and are similarly affected by the $R_0$ settings. We also compared the $r_{sp}$ inferred from the high-richness clusters with what is expected from the dark matter simulations and did not see significantly improved agreement. These results show that with the statistical uncertainties in our data set, we do not gain by switching to a higher-richness sample.

### B.2. Effect of $R_0$ on $r_{sp}$ in Dynamical Friction Measurements

Here we test how the $R_0$ settings affect our measurements in Section 6. To do this, we repeat the measurements in Section 6 using the two cluster catalogs described in Section 7, which were derived using different $R_0$ values. The resulting measurements are shown in Figure 14. We find that with $R_0 = 1.25 h^{-1}$ Mpc, the three galaxy luminosity bins show
<figure><img src="image_16.png" /><figcaption>**Figure 15.** Log derivative of galaxy density profiles calculated from model fits for varying choices of $h_{max}$, the maximum distance along the line of sight to integrate when converting the 3D profile into a projected profile.</figcaption></figure>

similar trends as our fiducial case of $R_0 = 1 h^{-1} Mpc$, where there is a hint of dynamical friction but at lower significance. For the $R_0 = 0.75 h^{-1} Mpc$ case, the three-galaxy sample gives consistent $r_{sp}$ values, and no sign of dynamical friction is seen. These findings again show that the measurement of $r_{sp}$ is sensitive to the choice of $R_0$.

## Appendix C
### Effect of Varying $h_{max}$

As discussed in Section 3, our model for the projected density profile, $\Sigma(R)$, is obtained by integrating the 3D profile, $\rho(r)$, along the line of sight. Throughout this analysis, we impose a maximum line-of-sight integration distance of $h_{max} = 40 h^{-1} Mpc$. Figure 15 shows the effect on our results of varying $h_{max}$. In general, we find that our inferences about the 3D profile are quite insensitive to the choice of $h_{max}$. There is a small change in the inferred slope of the outer density profile, but we note that most of our main results are not sensitive to the precise value of this outer density profile.

## ORCID iDs

C. Chang https://orcid.org/0000-0002-7887-0896

E. Baxter https://orcid.org/0000-0002-6836-3196

E. S. Rykoff https://orcid.org/0000-0001-9376-3135

A. Kravtsov https://orcid.org/0000-0003-4307-634X

D. Gruen https://orcid.org/0000-0003-3270-7644

E. M. Huff https://orcid.org/0000-0002-9378-3424

A. Palmese https://orcid.org/0000-0002-6011-0530

D. Rapetti https://orcid.org/0000-0003-2196-6675

R. H. Wechsler https://orcid.org/0000-0003-2229-011X

J. Annis https://orcid.org/0000-0002-0609-3987

A. Carnero Rosell https://orcid.org/0000-0003-3044-5150

A. Drlica-Wagner https://orcid.org/0000-0001-8251-933X

D. W. Gerdes https://orcid.org/0000-0001-6942-2736

R. A. Gruendl https://orcid.org/0000-0002-4588-6517

P. Martini https://orcid.org/0000-0002-4279-4182

P. Melchior https://orcid.org/0000-0002-8873-5065

R. Miquel https://orcid.org/0000-0002-6610-4836

B. Nord https://orcid.org/0000-0001-6706-8972

R. C. Smith https://orcid.org/0000-0001-5303-6830

M. Soares-Santos https://orcid.org/0000-0001-6082-8529

G. Tarle https://orcid.org/0000-0003-1704-0781

## References

Adhikari, S., Dalal, N., & Chamberlain, R. T. 2014, JCAP, 11, 019

Adhikari, S., Dalal, N., & Clampitt, J. 2016, JCAP, 7, 022

Aihara, H., Allende Prieto, C., An, D., et al. 2011, ApJS, 193, 29

Baxter, E., Chang, C., Jain, B., et al. 2017, ApJ, 841, 18

Baxter, E. J., Raghunathan, S., Crawford, T. M., et al. 2018, MNRAS, 476, 2674

Behroozi, P. S., Wechsler, R. H., & Wu, H.-Y. 2013, ApJ, 762, 109

Benítez, N. 2000, ApJ, 536, 571

Bertin, E., & Arnouts, S. 1996, A&AS, 117, 393

Bertschinger, E. 1985, ApJS, 58, 39

Binney, J., & Tremaine, S. 2008, Galactic Dynamics (2nd ed.; Princeton, NJ: Princeton Univ. Press)

Budzynski, J. M., Koposov, S. E., McCarthy, I. G., McGee, S. L., & Belokurov, V. 2012, MNRAS, 423, 104

Busch, P., & White, S. D. M. 2017, MNRAS, 470, 4767

Chandrasekhar, S. 1949, RvMP, 21, 383

Chaves-Montero, J., Angulo, R. E., Schaye, J., et al. 2016, MNRAS, 460, 3100

Cooray, A., & Sheth, R. 2002, PhR, 372, 1

Diehl, H. T., Abbott, T. M. C., Annis, J., et al. 2014, Proc. SPIE, 9149, 91490V

Diemer, B. 2017, arXiv:1712.04512

Diemer, B., & Kravtsov, A. V. 2014, ApJ, 789, 1

Diemer, B., Mansfield, P., Kravtsov, A. V., & More, S. 2017, ApJ, 843, 140

Drlica-Wagner, A., Sevilla-Noarbe, I., Rykoff, E. S., et al. 2018, ApJS, 235, 33

Einasto, J. 1965, TrAlm, 5, 87

Fillmore, J. A., & Goldreich, P. 1984, ApJ, 281, 1

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Gao, L., Navarro, J. F., Cole, S., et al. 2008, MNRAS, 387, 536

Gunn, J. E., & Gott, J. R., III 1972, ApJ, 176, 1

Hoyle, B., Gruen, D., Bernstein, G. M., et al. 2018, MNRAS, 478, 592

Huff, E., & Mandelbaum, R. 2017, arXiv:1702.02600

Johnston, D. E., Sheldon, E. S., Wechsler, R. H., et al. 2007, arXiv:0709.1159

Landy, S. D., & Szalay, A. S. 1993, ApJ, 412, 64

MacQueen, J. 1967, in Proc. Fifth Berkeley Symp. Mathematical Statistics and Probability, Vol. 1 (Berkeley, CA: Univ. California Press), 281, http://projecteuclid.org/euclid.bsmsp/1200512992

Mansfield, P., Kravtsov, A. V., & Diemer, B. 2017, ApJ, 841, 34

Melchior, P., Gruen, D., McClintock, T., et al. 2017, MNRAS, 469, 4899

Melchior, P., Suchyta, E., Huff, E., et al. 2015, MNRAS, 449, 2219

Miyatake, H., More, S., Takada, M., et al. 2016, PhRvL, 116, 041301

More, S., Diemer, B., & Kravtsov, A. V. 2015, ApJ, 810, 36

More, S., Miyatake, H., Takada, M., et al. 2016, ApJ, 825, 39

Murata, R., Nishimichi, T., Takada, M., et al. 2018, ApJ, 854, 120

Nagai, D., & Kravtsov, A. V. 2005, ApJ, 618, 557

Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563

Prada, F., Klypin, A. A., Cuesta, A. J., Betancort-Rijo, J. E., & Primack, J. 2012, MNRAS, 423, 3018

Reddick, R. M., Wechsler, R. H., Tinker, J. L., & Behroozi, P. S. 2013, ApJ, 771, 30

Riebe, K., Partl, A. M., Enke, H., et al. 2013, AN, 334, 691

Rozo, E., Rykoff, E. S., Koester, B. P., et al. 2009, ApJ, 703, 601

Rykoff, E. S., Rozo, E., Busha, M. T., et al. 2014, ApJ, 785, 104

Rykoff, E. S., Rozo, E., Hollowood, D., et al. 2016, ApJS, 224, 1

Sheldon, E. S., Johnston, D. E., Frieman, J. A., et al. 2004, AJ, 127, 2544

Springel, V., White, S. D. M., Jenkins, A., et al. 2005, Natur, 435, 629

Tinker, J., Kravtsov, A. V., Klypin, A., et al. 2008, ApJ, 688, 709

Troxel, M. A., MacCrann, N., Zuntz, J., et al. 2017, arXiv:1708.01538

Umetsu, K., & Diemer, B. 2017, ApJ, 836, 231

Yang, X., Mo, H. J., van den Bosch, F. C., et al. 2006, MNRAS, 373, 1159

Zu, Y., Mandelbaum, R., Simet, M., Rozo, E., & Rykoff, E. S. 2016, arXiv:1611.00366

Zuntz, J., Kacprzak, T., Voigt, L., et al. 2013, MNRAS, 434, 1604

Zuntz, J., Sheldon, E., Samuroff, S., et al. 2017, arXiv:1708.01533