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

# Probing Galaxy Evolution in Massive Clusters Using ACT and DES: Splashback as a Cosmic Clock

Susmita Adhikari¹,² ID, Tae-hyeon Shin³, Bhuvnesh Jain³, Matt Hilton⁴,⁵ ID, Eric Baxter⁶ ID, Chihway Chang⁷,⁸ ID,
Risa H. Wechsler¹,² ID, Nick Battaglia⁹, J. Richard Bond¹⁰ ID, Sebastian Bocquet¹¹ ID, Steve K. Choi⁹,¹² ID,
Joseph DeRose¹³,¹⁴,¹⁵ ID, Mark Devlin³ ID, Jo Dunkley¹⁶,¹⁷, August E. Evrard¹⁸, Simone Ferraro¹⁹ ID, J. Colin Hill²⁰,²¹,
John P. Hughes²² ID, Patricio A. Gallardo⁹, Martine Lokken¹⁰,²³,²⁴, Amanda MacInnis²⁵, Mathew S. Madhavacheril²⁶ ID,
Jeffrey McMahon⁷,⁸,²⁷,²⁸, Frederico Nati²⁹, Laura B. Newburgh³⁰, Michael D. Niemack⁹,¹², Lyman A. Page³¹ ID,
Antonella Palmese⁷,³² ID, Bruce Partridge³³, Eduardo Rozo³⁴, Eli Rykoff¹,² ID, Maria Salatino¹, Alessandro Schillaci³⁵,
Neelima Sehgal²⁵ ID, Cristóbal Sifón³⁶ ID, Chun-Hao To¹,², Ed Wollack³⁷ ID, Hao-Yi Wu³⁸,³⁹, Zhilei Xu³ ID, Michel Aguena⁴⁰,⁴¹,
Sahar Allam³², Alexandra Amon¹,², James Annis³² ID, Santiago Avila⁴² ID, David Bacon⁴³, Emmanuel Bertin⁴⁴,⁴⁵ ID,
Sunayana Bhargava⁴⁶, David Brooks⁴⁷ ID, David L. Burke¹,² ID, Aurelio C. Rosell⁴⁸,⁴⁹, Matias Carrasco Kind⁵⁰,⁵¹ ID,
Jorge Carretero⁵² ID, Francisco Javier Castander⁵³,⁵⁴ ID, Ami Choi³⁸, Matteo Costanzi⁵⁵,⁵⁶ ID, Luiz N. da Costa⁴¹,⁵⁷ ID,
Juan De Vicente⁵⁸, Shantanu Desai⁵⁹ ID, Thomas H. Diehl³², Peter Doel⁴⁷, Spencer Everett⁶⁰, Ismael Ferrero⁶¹, Agnès Ferté⁶²,
Brenna Flaugher³², Pablo Fosalba⁵³,⁵⁴ ID, Josh Frieman⁷,³², Juan García-Bellido⁴² ID, Enrique Gaztanaga⁵³,⁵⁴ ID, Daniel Gruen¹,² ID,
Robert A. Gruendl⁵⁰,⁵¹ ID, Julia Gschwend⁴¹,⁵⁷ ID, Gaston Gutierrez³², Will G. Hartley⁴⁷,⁶³,⁶⁴ ID, Samuel R. Hinton⁶⁵ ID,
Devon L. Hollowood⁶⁰ ID, Klaus Honscheid³⁸,⁶⁶, David J. James⁶⁷, Tesla Jeltema⁶⁰, Kyler Kuehn⁶⁸,⁶⁹ ID, Nikolay Kuropatkin³² ID,
Ofer Lahav⁴⁷, Marcos Lima⁴⁰,⁴¹, Marcio A. G. Maia⁴¹,⁵⁷, Jennifer L. Marshall⁷⁰ ID, Paul Martini³⁸,⁷¹, Peter Melchior⁷² ID,
Felipe Menanteau⁵⁰,⁵¹ ID, Ramon Miquel⁵²,⁷³ ID, Robert Morgan⁷⁴, Ricardo L. C. Ogando⁴¹,⁵⁷, Francisco Paz-Chinchón⁷⁵ ID,
Andrés Plazas Malagón⁷², Eusebio Sanchez⁵⁸ ID, Basilio Santiago⁴¹,⁵¹,⁷⁶, Vic Scarpine³², Santiago Serrano⁵³,⁵⁴ ID,
Ignacio Sevilla-Noarbe⁵⁸ ID, Mathew Smith⁷⁷ ID, Marcelle Soares-Santos⁷⁸ ID, Eric Suchyta⁷⁹, Molly E. C. Swanson⁵¹ ID,
Tamas N. Varga⁸⁰,⁸¹, Reese D. Wilkinson⁴⁶, Yuanyuan Zhang³², Jason E. Austermann⁸², James A. Beall⁸², Daniel T. Becker⁸²,
Edward V. Denison⁸², Shannon M. Duff⁸², Gene C. Hilton⁸², Johannes Hubmayr⁸², Joel N. Ullom⁸², Jeff Van Lanen⁸²,
Leila R. Vale⁸²
(DES Collaboration)
and
(ACT Collaboration)

¹ Kavli Institute for Particle Astrophysics and Cosmology and Department of Physics, Stanford University, Stanford, CA 94305, USA; susmita@stanford.edu

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

³ Department of Physics and Astronomy, University of Pennsylvania, 209 S. 33rd Street, Philadelphia, PA 19104, USA; taeshin@sas.upenn.edu

⁴ Astrophysics Research Centre, University of KwaZulu-Natal, Westville Campus, Durban 4041, South Africa

⁵ School of Mathematics, Statistics & Computer Science, University of KwaZulu-Natal, Westville Campus, Durban 4041, South Africa

⁶ Institute for Astronomy, University of Hawaii, Manoa, Hawaii, USA

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

⁸ Department of Physics, The University of Chicago, Chicago, IL 60637, USA

⁹ Department of Astronomy, Cornell University, Ithaca, NY 14853, USA

¹⁰ Canadian Institute for Theoretical Astrophysics, 60 St. George Street, University of Toronto, Toronto, ON, M5S 3H8, Canada

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

¹² Department of Physics, Cornell University, Ithaca, NY, 14853 USA

¹³ Department of Astronomy and Astrophysics, University of California, Santa Cruz, 1156 High Street, Santa Cruz, CA 95064 USA

¹⁴ Berkeley Center for Cosmological Physics, Berkeley, CA 94720, USA

¹⁵ Department of Physics, University of California, Berkeley, CA 94720, USA

¹⁶ Department of Astrophysical Sciences, Princeton University, Princeton, NJ, 08544, USA

¹⁷ Department of Physics, Princeton University, Princeton, NJ 08544, USA

¹⁸ Department of Astronomy, University of Michigan, MI, USA

¹⁹ Lawrence Berkeley National Laboratory, One Cyclotron Road, Berkeley, CA 94720, USA

²⁰ Department of Physics, Columbia University, New York, NY, 10027 USA

²¹ Center for Computational Astrophysics, Flatiron Institute, New York, NY, 10010 USA

²² Department of Physics and Astronomy, Rutgers, the State University of New Jersey, 136 Frelinghuysen Road, Piscataway, NJ 08854-8019, USA

²³ Dunlap Institute for Astronomy and Astrophysics, University of Toronto, 50 St. George Street, Toronto, ON M5S 3H4, Canada

²⁴ Department of Astronomy and Astrophysics, University of Toronto, 50 St. George Street, Toronto, ON M5S 3H4, Canada

²⁵ Physics and Astronomy Department, Stony Brook University, Stony Brook, NY 11794, USA

²⁶ Centre for the Universe, Perimeter Institute, Waterloo, ON N2L 2Y5, Canada

²⁷ Department of Astronomy and Astrophysics, University of Chicago, Chicago, IL 60637, USA

²⁸ Enrico Fermi Institute, University of Chicago, Chicago, IL 60637, USA

²⁹ Department of Physics, University of Milano-Bicocca, Piazza della Scienza 3, I-20126 Milano, Italy

³⁰ Department of Physics, Yale University, 217 Prospect Street, New Haven, CT 06511, USA

³¹ Joseph Henry Laboratories of Physics, Jadwin Hall, Princeton University, Princeton, NJ 08544, USA

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

³³ Department of Physics and Astronomy, Haverford College, Haverford, PA, 19041, USA

³⁴ Department of Physics, University of Arizona, 1118 E 4th Street, Tucson, AZ 85721 USA

³⁵ Department of Physics, California Institute of Technology, Pasadena, CA, USA
³⁶ Instituto de Física, Pontificia Universidad Católica de Valparaíso, Casilla 4059, Valparaíso, Chile

³⁷ Code 665, NASA/Goddard Space Flight Center, Greenbelt, MD 20771, USA

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

³⁹ Department of Physics, Boise State University, Boise, ID 83725, USA

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

⁴¹ Laboratório Interinstitucional de e-Astronomia—LineA, Rua Gal. José Cristino 77, Rio de Janeiro, RJ—20921-400, Brazil

⁴² Instituto de Física Teorica UAM/CSIC, Universidad Autonoma de Madrid, E-28049 Madrid, Spain

⁴³ Institute of Cosmology and Gravitation, University of Portsmouth, Portsmouth, PO1 3FX, UK

⁴⁴ 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

⁴⁶ Department of Physics and Astronomy, Pevensey Building, University of Sussex, Brighton, BN1 9QH, UK

⁴⁷ Department of Physics & Astronomy, University College London, Gower Street, London, WC1E 6BT, UK

⁴⁸ Instituto de Astrofísica de Canarias, E-38205 La Laguna, Tenerife, Spain

⁴⁹ Universidad de La Laguna, Dpto. Astrofísica, E-38206 La Laguna, Tenerife, Spain

⁵⁰ Department of Astronomy, University of Illinois at Urbana-Champaign, 1002 W. Green Street, Urbana, IL 61801, USA

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

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

⁵³ Institut d'Estudis Espacials de Catalunya (IEEC), 08034 Barcelona, Spain

⁵⁴ Institute of Space Sciences (ICE, CSIC), Campus UAB, Carrer de Can Magrans, s/n, E-08193 Barcelona, Spain

⁵⁵ INAF-Osservatorio Astronomico di Trieste, via G.B. Tiepolo 11, I-34143 Trieste, Italy

⁵⁶ Institute for Fundamental Physics of the Universe, Via Beirut 2, I-34014 Trieste, Italy

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

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

⁵⁹ Department of Physics, IIT Hyderabad, Kandi, Telangana 502285, India

⁶⁰ Santa Cruz Institute for Particle Physics, Santa Cruz, CA 95064, USA

⁶¹ Institute of Theoretical Astrophysics, University of Oslo, P.O. Box 1029 Blindern, NO-0315 Oslo, Norway

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

⁶³ Département de Physique Théorique and Center for Astroparticle Physics, Université de Genève, 24 quai Ernest Ansermet, CH-1211 Geneva, Switzerland

⁶⁴ Department of Physics, ETH Zurich, Wolfgang-Pauli-Strasse 16, CH-8093 Zurich, Switzerland

⁶⁵ School of Mathematics and Physics, University of Queensland, Brisbane, QLD 4072, Australia

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

⁶⁷ Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA

⁶⁸ Australian Astronomical Optics, Macquarie University, North Ryde, NSW 2113, Australia

⁶⁹ Lowell Observatory, 1400 Mars Hill Road, Flagstaff, AZ 86001, USA

⁷⁰ 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

⁷⁴ Physics Department, 2320 Chamberlin Hall, University of Wisconsin-Madison, 1150 University Avenue, Madison, WI 53706-1390, USA

⁷⁵ Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK

⁷⁶ Instituto de Física, UFRGS, Caixa Postal 15051, Porto Alegre, RS—91501-970, Brazil

⁷⁷ School of Physics and Astronomy, University of Southampton, Southampton, SO17 1BJ, UK

⁷⁸ Department of Physics, University of Michigan, Ann Arbor, MI 48109, USA

⁷⁹ Computer Science and Mathematics Division, Oak Ridge National Laboratory, Oak Ridge, TN 37831, USA

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

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

⁸² National Institute of Standards and Technology, Boulder, CO, USA

Received 2020 September 10; revised 2021 May 30; accepted 2021 June 10; published 2021 December 8

## Abstract

We measure the projected number density profiles of galaxies and the splashback feature in clusters selected by the Sunyaev–Zel’dovich effect from the Advanced Atacama Cosmology Telescope (AdvACT) survey using galaxies observed by the Dark Energy Survey (DES). The splashback radius is consistent with CDM-only simulations and is located at $2.4_{-0.4}^{+0.3}$ Mpc $h^{-1}$. We split the galaxies on color and find significant differences in their profile shapes. Red and green-valley galaxies show a splashback-like minimum in their slope profile consistent with theory, while the bluest galaxies show a weak feature at a smaller radius. We develop a mapping of galaxies to subhalos in simulations and assign colors based on infall time onto their hosts. We find that the shift in location of the steepest slope and different profile shapes can be mapped to the average time of infall of galaxies of different colors. The steepest slope traces a discontinuity in the phase space of dark matter halos. By relating spatial profiles to infall time, we can use splashback as a clock to understand galaxy quenching. We find that red galaxies have on average been in clusters over 3.2 Gyr, green galaxies about 2.2 Gyr, while blue galaxies have been accreted most recently and have not reached apocenter. Using the full radial profiles, we fit a simple quenching model and find that the onset of galaxy quenching occurs after a delay of about a gigayear and that galaxies quench rapidly thereafter with an exponential timescale of 0.6 Gyr.

Unified Astronomy Thesaurus concepts: Galaxy clusters (584); Cold dark matter (265); Galaxy dark matter halos (1880); Galaxy kinematics (602); Galaxy evolution (594); Surveys (1671)
# 1. Introduction

Galaxy clusters, the largest bound objects in the universe, are a unique laboratory for studying the nonlinear evolution of the universe. These objects can be observed across the electromagnetic spectrum in X-ray, submillimeter, and optical, with each band illuminating different components of the cluster. Large-scale surveys of the sky have provided us with statistical samples of clusters that give us the opportunity to study the coevolution of galaxies and matter in these extreme environments in great detail. Understanding the evolution of the galaxy populations that make up a cluster and the underlying dark matter potential can provide fundamental insights into the connection between dark matter and galaxies.

The properties of galaxies in a cluster are significantly different from those in the field. One striking feature of clusters is that they appear to be dominated by red galaxies, with little or no star formation (Oemler 1974; Dressler 1980; Dressler & Gunn 1983; Balogh et al. 1997; Poggianti et al. 1999). Intracluster processes like ram pressure stripping (Gunn & Gott 1972; Abadi et al. 1999), strangulation (Larson et al. 1980), harassment, or tidal disruption are known to quench star formation in galaxies as they orbit within the cluster. Understanding the response of baryonic matter to the dark matter halo environment and the evolution of star formation has been a long-standing question in galaxy evolution and cosmology. The assembly and evolution of galaxies are expected to be closely related to the evolution of its parent dark matter halo (Bullock et al. 2002; Wechsler et al. 2002; Moustakas & Somerville 2002; Conroy et al. 2006; Cooray 2006; Conroy & Wechsler 2009; Kravtsov & Borgani 2012; Somerville & Davé 2015; Wechsler & Tinker 2018). While the total spatial distribution of galaxies in a cluster halo is determined by the dark matter potential, studying spatial correlations between various galaxy properties is a powerful way to infer the coevolution of dark matter and galaxies.

In recent years, detailed studies of the galaxy distribution within clusters have emerged as a robust way to gain insight into dynamic properties of the halo. In particular, examining the logarithmic slope of the cluster halo profile highlights a novel feature in the density profile termed the splashback radius, where the slope reaches a minimum in a narrow localized region. The splashback radius relates to the splashback surface, which can be thought of as the physical boundary of a halo (Diemer & Kravtsov 2014; Adhikari et al. 2014; More et al. 2015; Shi 2016; Mansfield et al. 2017). The splashback surface traces the location of the apocenters in the orbits of the most recently accreted matter onto a dark matter halo (Adhikari et al. 2014; Diemer 2017; Diemer et al. 2017). This boundary therefore separates the multistreaming or “virialized” region of a halo, where particles or satellite galaxies are in orbits, from the region where there are only infalling galaxies or particles (see also Aung et al. 2021). The splashback surface denotes a sharp boundary in phase space, and it has been found that the slope of the spherically averaged, logarithmic matter density profile reaches a minimum at this boundary. This can be explained by the fact that the density in the virialized region of the halo follows a Navarro, Frenk, and White (NFW)-like profile that approaches a slope of -3 in the outer regions (Navarro et al. 1996), while the density of particles in the infall region follows a power law with index ~-1.5 (Baxter et al. 2017); a sharp transition between the virialized and infalling regimes at splashback therefore causes the slope to reach a

minimum in a narrow, localized region before it asymptotes back to the background value.

The splashback radius, as a feature in the spherically averaged density profile of both matter and galaxies around dark matter halos, can be accessed observationally. In several recent studies (More et al. 2016; Baxter et al. 2017; Chang et al. 2018; Shin et al. 2019; Zürcher & More 2019; Murata et al. 2020), the splashback radius was measured in the projected number density of galaxies around massive clusters selected both optically and using the Sunyaev–Zel’dovich (SZ; Sunyaev & Zel’dovich 1972) effect. The splashback radius has also been measured in the dark matter distribution itself, using weak lensing of galaxies around massive clusters; while Umetsu & Diemer (2017) and Contigiani et al. (2019) have attempted to measure the feature using a small sample of massive X-ray selected clusters, Chang et al. (2018) have used a large, statistical sample of optically selected, RedMaPPer clusters (Rykoff et al. 2014) in the Dark Energy Survey (DES; Abbott et al. 2005) data. The latter finds evidence for a significant steepening of the density profile at the edge of galaxy clusters. Xhakaj et al. (2020) find that future surveys like LSST (Abell et al. 2009) and Euclid (Laureijs et al. 2011) will measure the splashback radius through weak lensing with nearly 10% precision. In Tomooka et al. (2020) a sharp boundary was also measured in the line-of-sight velocity dispersion of spectroscopic galaxies around clusters in SDSS (Abazajian et al. 2009). Okumura et al. (2017, 2018) have also explored theoretical predictions for the splashback radius in the observed velocity field. All of these studies provide evidence for the existence of a sharp transition around cluster boundaries, a feature that must exist in the presence of dark matter. Furthermore, studies have shown that the location of this feature can be sensitive to the nature of gravity and dark matter itself (Adhikari et al. 2018; Banerjee et al. 2020).

While the splashback radius is observed as a feature in the spatial distribution of matter and galaxies (and their velocities), it also has an inherent timescale associated with it: the time for a particle or galaxy to reach the apocenter of its first orbit. In other words, dark matter particles or galaxies that form the splashback region have been inside the halo for approximately one orbital time, from accretion to first apocenter. Unlike dark matter particles, galaxy properties like star formation rates (SFRs) and morphology can evolve significantly on similar timescales if the cluster environment plays a role in galaxy evolution. Because the splashback feature is closely related to the orbital time, it is possible to use it as a clock to study different populations of galaxies and their time-evolving properties.

A galaxy within a cluster evolves from being blue and star-forming to red over the course of its multiple orbits through the halo. Depending on the specific mechanism of quenching, the quenching timescale can vary. The longer a galaxy orbits within a halo, the more likely it is to be affected by the intracluster medium and evolve in color, which in turn changes the radial profiles of galaxies binned in color. Thus by using the shapes of the profiles of galaxies of different colors within the splashback radius, we can learn about the quenching timescale. In the scenario that a population of galaxies has not reached apocenter (for example, if we assume that all blue galaxies quench at or right after pericenter passage), we should not expect to see a true splashback feature in its density profile. We
show below that this applies to the bluest galaxies in our sample.

In this paper we investigate the relation between galaxy evolution and density distributions within dark matter halos, focusing on information contained in the splashback feature. In particular, we study whether the location of the minimum in the profile's slope can provide additional insight into the properties of galaxies with different SFRs. We study how this information can help us relate the distribution of observed galaxies to the dark matter phase space.

Our analysis relies on a catalog of galaxy clusters (Hilton et al. 2021) that have been observed via their SZ signal in submillimeter-wave maps from the Atacama Cosmology Telescope (ACT; Aiola et al. 2020; Naess et al. 2020). SZ selection is essentially redshift independent and is relatively closer to a mass-limited sample because of the small scatter in the relationship between cluster mass and SZ observable (Nagai et al. 2007). Moreover, clusters selected with the SZ effect are known to suffer less from systematic effects such as line-of-sight projections and triaxiality than optically selected clusters (e.g., Shin et al. 2019). We study the distribution of galaxies around the SZ-selected clusters from ACT by cross-correlating the cluster positions with galaxies from DES data (Abbott et al. 2005).

This paper is organized as follows. In Section 2 we study and review the relation between infall time, density profiles, and splashback radius from N-body simulations. In Section 3 we describe observations of splashback radius using galaxy clusters selected in SZ data, as a function of galaxy color; in Section 4 we model the quenching of star-forming galaxies as a function of infall time to infer relevant timescales for cluster-mass halos; and we summarize our conclusions in Section 5.

## 2. Phase Space and Infall Time of Subhalos and Dark Matter in Simulations

The collisionless dark matter particles that collapse gravitationally to form dark matter halos occupy a very specific region in the phase space of halocentric radius and radial velocity. Particles, as they fall into the halo, follow a trajectory in this space as a function of time, forming what is known as the “multisreaming” region of the halo. In other words, the multisreaming region is where particles are orbiting in the potential of the halo, and at any point in space, there are multiple streams of particles, at different velocities, that have entered the halo at different times. The spherically averaged 3D density profile is the integral of the phase-space distribution over the radial velocity. The dynamics of galaxies that fall into cluster halos largely follow the collisionless dark matter particles and can also be traced in this space. In this section, we study the evolution of particles and substructure around cluster halos in phase space and infer how observed density profiles are composed of components accreted at different times.

### 2.1. Simulations

We use two separate simulations to study halo evolution. To study dark matter particles we use a cosmological N-body simulation of cold dark matter (CDM) in a 1 Gpc $h^{-1}$ box, with 1024³ particles, run using Gadget2 (Springel 2005). We study a sample of cluster halos selected based on their $M_{500c}$, the mass within the radius that encloses 500 times the critical density of the universe. We apply a minimum mass threshold such that the

<table><caption>**Table 1**
Summary of Important Length Scales and Timescales Used in The Paper</caption><thead><tr><td>Variable</td><td>Definition</td></tr></thead><tbody><tr><td>r<sub>sp</sub></td><td>Location of the 3D splashback radius</td></tr><tr><td>r<sub>min</sub></td><td>Location of the minimum in the slope profile</td></tr><tr><td>t<sub>in</sub></td><td>Time since infall/accretion onto the halo boundary</td></tr><tr><td>z<sub>acc</sub></td><td>Redshift at accretion</td></tr><tr><td>t<sub>max</sub></td><td>Maximum infall time threshold</td></tr><tr><td>t<sub>d</sub></td><td>Delay time in quenching model</td></tr><tr><td>t<sub>q</sub></td><td>Exponential quenching timescale in quenching model</td></tr></tbody></table>

mean mass of our sample $\langle M_{500c} \rangle = 3.1 \times 10^{14} M_{\odot} h^{-1}$. We extract all particles around the cluster halos at $z=0$ and find the orbits of subsampled sets of 1000 particles, randomly selected, within 5 Mpc $h^{-1}$ of each object. To study the distribution of subhalos and halos around clusters, we require simulations with higher resolution; for this we use the publicly available Multidark Planck (MDPL2, Riebe et al. 2013) simulation. This is a CDM-only simulation of a 1 Gpc³ $h^{-3}$ volume with 3840³ particles. The simulation assumes the best-fitting flat $\Lambda$CDM Planck cosmological model (Ade et al. 2014), with $\Omega_m = 0.307$ and $h = 0.677$. The properties of the halos and their subhalos were obtained using the ROCKSTAR halo finder (Behroozi et al. 2013) in both simulations.

Henceforth, we will refer to the host, cluster-mass dark matter halos as “clusters” to avoid confusion while referring to the other halos around these cluster-mass objects. This is primarily because we study the overall distribution of collapsed structures around clusters in the subsequent sections, which include both subhalos, within the virial radius, and other halos in the vicinity of the cluster extending beyond the virial radius. We use several timescales and length scales in this paper; a summary of variables is provided in Table 1 for reference. We describe them in the text where they appear.

### 2.2. Particles

In this section we study the distribution of particles around clusters in N-body simulations. In Figure 1 we show the location of dark matter particles around clusters in the $r-v_{rad}$ plane, where r is the distance from the center of the cluster and $v_{rad}$ is the radial velocity with respect to the center. The three top panels correspond to particle populations that were accreted onto the cluster at different times. The particles are separated based on when, in their orbital history, they first crossed within 4 Mpc $h^{-1}$ of the host cluster.<fnref n="83" /> The figure shows the stacked phase-space distribution of particles from cluster halos at $z=0$ and with $\langle M_{500c} \rangle = 3.1 \times 10^{14} M_{\odot} h^{-1}$. This mass cut is chosen to match the mean mass of the observed galaxy cluster sample described in Section 3.1.1. Here we only focus on clusters at a single redshift for simplicity.

In the top-left panel we see particles that have been accreted onto the cluster halo most recently, with $z_{acc} < 0.25$, where $z_{acc}$ is the redshift at which the particle crossed the boundary at 4 Mpc $h^{-1}$. It appears that these particles mainly occupy the first “infall” stream in the region with $v_{rad} < 0$. This population corresponds to all particles that have infall times, $t_{in} < t_{max}$, where $t_{max}$ is the lookback time corresponding to upper-limit on

<footnote n="83">⁸³ We choose 4 Mpc $h^{-1}$ to be safely outside the virial radius of the halo, which is typically ~1.5–2 Mpc $h^{-1}$ for the halo masses considered in this work.</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** (Top) Stacked phase space of particles accreted with different maximum time thresholds. Moving from large to small clustercentric radius, the leftmost top panel corresponds to particles that have been within 4 Mpc $h^{-1}$ of the host for less than 3 Gyr, the second and third panels correspond to particles that have been inside for less than 5.6 and 9 Gyr, respectively. The particles that have been inside for the least amount of time are present only in the infall stream (left panel). The middle panel shows a discontinuity (transition from multistream to infall stream), shown by the white vertical line, at a distance smaller than the traditional splashback radius. The actual splashback radius (right panel) is traced only by particles that have been inside the cluster long enough to reach apocenter. (Bottom) The left panel shows the density profiles for the three cases; the right panel shows the logarithmic slope. The location of the minimum in the slope profile, $r_{smin}$, coincides with the location of the phase-space discontinuity shown by the vertical lines.</figcaption></figure>

$z_{acc}$, 3.1 Gyr. The infall or accretion time, $t_{in}$, is also the time since crossing the 4 Mpc $h^{-1}$ boundary. The crossing time to the center of the halo from the boundary is ~3 Gyr, indicating that these particles have mostly not crossed pericenter. On the other hand, if we increase the maximum time threshold, $t_{max}$, to 5.6 Gyr (top row, middle panel), such that we now have a population of particles that have been within the boundary for, at most, 5.6 Gyr, we find that they live in a stream that has wrapped around in phase space once. This population occupies the $v_{rad} < 0$ infall stream and also the region past pericenter with positive radial velocity, $v_{rad} > 0$. We note that these particles have not reached apocenter or splashback in their orbits (at which $v_{rad} \sim 0$ km s$^{-1}$). If the maximum time threshold, $t_{max}$, is increased further to show the phase-space distribution of all particles that have been inside the cluster boundary for more than 9 Gyr, for example, we find that particles have wrapped around in multiple streams (top-right panel), forming the complete virialized region of the halo.

The bottom-row panels of Figure 1 show the spherically averaged 3D density, and its logarithmic slope as a function of radius, for the particles in the three panels of the top row. The slopes have been measured by smoothing the density profile using a Savitsky–Golay smoothing filter (Diemer & Kravtsov 2014). We note, in the left panel, that the shape of the density profile changes as a function of the maximum infall time threshold. From the slope profiles, shown in the right panel, we infer that the density of particles in the infall stream (lightest gray curve) is well approximated by a single power law, but the same is not true for the other two cases. In fact, we find that the slopes show a distinct minimum in $r$ for particles that have been inside the halo longer than the time required for pericentric passage. The location of the minimum of the slope shifts for particles that have been accreted at different time thresholds. The white dashed vertical lines in the top three panels show the location of the minimum evaluated from the slope profiles in phase space. In particular, we find that the minimum appears to trace the boundary between the multistream
and single-stream region. This is most clearly demonstrated by comparing the middle and right panels in the top row; for the population of particles with $t_{in} < 5.6$ Gyr, the white dashed line separates the two-stream region on its left at small $r$ from the single infall stream on its right, whereas for particles with $t_{in} < 9$ Gyr, the minimum is at the location of the conventionally described splashback radius, i.e., at the boundary of the halo. The three particular values for the time thresholds were chosen to demonstrate the three distinct regimes.

We conclude that the minimum of the slope traces a phase-space discontinuity that shifts for tracer populations accreted at different times. If we can separate populations of objects that have accreted onto halos at different times, we may expect to be able to measure this shift. Henceforth, we define a new generalized quantity, $r_{smin}$, to refer to the location of the steepest slope. This may be distinct from the traditional splashback radius for particle or galaxy populations that have not reached the apocenter of their orbits. We note that this result was discussed in the context of subhalos in Shin et al. (2019); here we verify that the trend holds for particle distributions as well.

Here we use cumulative time bins or the maximum infall time thresholds, rather than differential bins, to connect with galaxy populations that enter the clusters at different times. For example, if we study the density of red galaxies around a cluster, we expect them to be present both in the infall stream and in the virialized region. In the next section, we continue to study the distribution of the halos of such galaxies (subhalos in cluster halos) in simulations.

## 2.3. Subhalos

Dark matter halos form hierarchically, in the sense that small objects merge to form massive structures. A cluster halo is a relatively young object and contains a lot of existing substructure or subhalos that have not been destroyed through multiple orbits within its potential. The overall dynamics of halos in a cluster environment is also dictated by the cluster potential and is similar to that of the dark matter particles. All halos, if they are sufficiently massive, are expected to host galaxies at the center of their potential wells, and therefore studying the evolution and distribution of substructure around clusters in simulations can help us understand the evolution of galaxies in observed clusters. As mentioned before, halos in our simulations are found using the code ROCKSTAR (Behroozi et al. 2013). The merger history of the halos are generated using the CONSISTENT-TREES algorithm (Behroozi et al. 2012). Among other properties, the final halo catalog provides the redshift at the time of accretion, $z_{acc}$, of a subhalo onto its parent host halo. This time is recorded when a subhalo crosses the virial radius of the parent cluster. In general, we study the entire field of halos around each cluster; this includes objects that have been tagged as subhalos of the cluster in the halo catalog and also others in the infall region. Halos that have not crossed the virial boundary have $z_{acc} = z_{cluster}$.

We study the distribution of substructure/halos around clusters in the MDPL2 simulation. We select a sample of clusters that matches the distribution of mass and redshift of our observed SZ sample (3.1.1). All clusters lie in the redshift range $0.15 < z < 0.7$ and have a mean mass $\langle M_{500c} \rangle = 3.1 \times 10^{14} M_{\odot} h^{-1}$. We extract all halos that have peak maximum circular velocity $v_{peak} > 150$ km s$^{-1}$ within a spherical volume of radius $5$ Mpc $h^{-1}$ around the cluster centers at every redshift.

The property $v_{peak}$ is a proxy for halo mass and is known to correlate best with the luminosity of the galaxy that it hosts (Reddick et al. 2013).

Halos that cross the virial radius of the cluster at different times occupy distinct regions in phase space in a manner similar to dark matter particles discussed in the previous section. We study the density profiles of these populations; in particular, we study how the 3D number density profile changes if we vary the maximum infall time, $t_{max}$, of a given population of halos around a cluster. We compute the stacked 3D number density of halos in every case as a function of clustercentric radius and measure the logarithmic slope of the profile using the Savitsky-Golay smoothing filter. We summarize our results in Figure 2. In the top panel we show the logarithmic slope of the number density profile of a population of halos around clusters that have infall times $t_{in} < t_{max}$; each curve corresponds to a different $t_{max}$. The radius $r_{smin}$ is the location of the minimum in the slope profile. The middle panel maps the movement of $r_{smin}$ with $t_{max}$. Halo populations that contain objects that have been accreted earlier than ~4 Gyr show a minimum in their profiles at the location of the cluster splashback radius, while $r_{smin}$ falls off to smaller values for more recently accreted populations. This shows that the location of the slope minimum, as is the case for particles, traces the location of discontinuity in phase space and contains information about the average infall time of a galaxy population. The dashed curve in the middle panel corresponds to the relation between $t_{max}$ and $r_{smin}$ for a sample of lower-mass halos with $\langle M_{500c} \rangle = 1.2 \times 10^{14} M_{\odot} h^{-1}$. The shape of the curve appears to be similar; the offset between the two curves arises from the overall smaller physical sizes of low-mass objects that scale approximately as $M^{1/3}$.

Dark matter particles are indistinguishable from each other, and probes of the matter distribution, like weak lensing, are sensitive to the total enclosed mass as opposed to populations of particles with different infall times. On the other hand, the observed properties of galaxies change with time after they are accreted onto clusters. Measuring the rate of change is often difficult because we do not know how long a particular population has been inside cluster. The methods outlined above allow us to use the number density profile of galaxies and its slope to probe the infall history of galaxies and provide a way to determine how long galaxies have been in a cluster. In the bottom panel of Figure 2, we show the spread of allowed curves corresponding to the $r_{smin}$ measured from three different populations of galaxies with varying colors. The following discussion in the paper describes the connection between these two pictures. In the next section we use insights from the evolution of subhalos in simulations to understand galaxy evolution in observed galaxy clusters using ACT and DES data.

## 3. Observations: Galaxy Profiles and Splashback Radius

In this section we describe the data and methods used for measuring the density distribution of galaxies and their splashback radius in massive galaxy clusters. Our procedure closely follows the method implemented in Shin et al. (2019); we refer readers to that work for details. Here we briefly summarize the method to measure profiles, describe our color selection criteria, and present our results.
<figure><img src="image_3.png" /><figcaption>**Figure 2.** (Top) The slope of the density profile as a function of radius for populations of subhalos that have different maximum infall time thresholds. Each curve corresponds to the profiles of subhalos that have been inside the virial radius for at most $t_{\text{max}}$ Gyr, i.e., $t_{\text{in}} < t_{\text{max}}$ (objects that have not entered the host at all have $t_{\text{in}} = 0$). (Middle) $t_{\text{max}}$ of a population of subhalos as a function of the location of the minimum of the slope, $r_{\text{smin}}$. The two separate curves correspond to different cluster-mass samples. The width of the blue, green, and red bands correspond to the $1\sigma$ measurements of $r_{\text{smin}}$ for the observed blue, green, and red galaxies (Section 4.1). (Bottom) Slope profiles for subhalo populations that have $r_{\text{smin}}$ in the range of values obtained from data (shaded curves) compared to the estimated timescales shown in the middle panel.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Distribution of ACT DR5 clusters in signal to noise (S/N) and redshift. Clusters with S/N > 4 and $0.15 < z < 0.7$ within the DES Y3 footprint, 908 in total, are shown in the figure and used in the analysis. The mean mass of the cluster sample is $3.1 \times 10^{14} M_{\odot}h^{-1}$.</figcaption></figure>

## 3.1. Data

### 3.1.1. ACT Cluster Catalog

We measure cluster–galaxy cross-correlations to estimate the splashback radius around massive galaxy clusters detected by the SZ effect from the Advanced Atacama Cosmology Telescope (AdvACT) survey. The cluster catalog is part of the fifth data release from ACT (ACT DR5; Hilton et al. 2021). The cluster catalog is derived from applying a multifrequency matched filter (e.g., Melin et al. 2006; Williamson et al. 2011) to 98 and 150 GHz ACT day- and night-time observations taken from 2008 to 2018. For this work, the ACT maps for each observing season and detector array were coadded using the procedure described in Naess et al. (2020).

The cluster signal is modeled using the universal pressure profile (UPP; Arnaud et al. 2010), and masses are inferred from the SZ signal in a method similar to Hilton et al. (2018). Here we make use of masses that have been rescaled according to a richness-based weak-lensing mass calibration procedure by a factor of $0.69 \pm 0.07$ (see Hilton et al. 2018). The survey selection function is defined using maps filtered at a single reference filter scale (2.4; equivalent to a cluster with $M_{500c} = 2 \times 10^{14} M_{\odot}$ at $z=0.4$). We apply a cut on the signal-to-noise ratio, S/N > 4, and redshift, $0.15 < z < 0.7$, which gives a total number of 908 clusters in the 4552 deg$^2$ overlapping area between DES and ACT. The mean cluster mass for this sample is $3.1 \times 10^{14} h^{-1} M_{\odot}$. The redshift and S/N distribution of the SZ clusters are shown in Figure 3. The redshifts for the clusters have been obtained by compiling spectroscopic redshifts, where available, and preexisting catalogs that have overlap with ACT (see Hilton et al. 2021, for a detailed description).

We also require a mock cluster catalog that closely follows the survey geometry of ACT, located at random positions in the sky to reliably measure cluster–galaxy cross-correlation functions. The mock “cluster randomness” were generated from the SZ-signal noise map generated by the cluster finder by
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Color selection method in DES data. The plot shows the excess number density of galaxies in the color–color plane around 1 Mpc h⁻¹ of the ACT DR5 clusters (with respect to random locations on the sky). The dotted–dashed lines are used to assign galaxy colors, red/green/blue, from upper right to lower left, respectively. See Section 3.2.1 for details. For comparison, the black dashed contour marks the 68% range of the red-sequence galaxies drawn from the optically selected clusters (RedMaPPer) in DES data.</figcaption></figure>

sampling from the Tinker et al. (2008) mass function and applying an SZ-signal–mass relation adjusted to reproduce the observed number of SZ S/N > 6 clusters found in the DES footprint (Hilton et al. 2021).

### 3.1.2. DES Galaxy Catalog

To measure the galaxy distribution around the selected clusters, we use galaxies observed in DES. DES is a 5000 deg² survey that covers the Southern Galactic cap. The survey used the 570 megapixel Dark Energy Camera (Flaugher et al. 2015), mounted on the 4 m Blanco Telescope in Chile, to image the sky in grizY filters. In this work we use data from the first three years of observation, in particular the DES Y3 gold catalog (Sevilla-Noarbe et al. 2021) similarly generated as in Drlica-Wagner et al. (2018). The image-processing pipeline used in DES is described in Morganson et al. (2018). We use the galaxy magnitudes that have been corrected for differential reddening across the DES footprint. After excluding galaxies with extreme colors (outside the range $-1 < (g-r) < 3$, $-1 < (r-i) < 2.5$, $-1 < (i-z) < 2$), we take all galaxies with i-band magnitude, $m_i$, smaller than 22.5 and only use parts of the footprint for which the depth of the survey in $m_i$ is larger than 22.5 to ensure the completeness; this leaves ~90% (~4500 deg²) of the entire DES footprint. We also require all galaxies to have uncertainties on the magnitude smaller than 0.1. When calculating the cross-correlation function between the clusters and the galaxies, we apply a further limit on the absolute magnitude, $M_i < -19.87$, which corresponds to the apparent magnitude limit of $m_i < 22.5$ at the maximum redshift of $z=0.7$. This is to ensure the same luminosity cut on the galaxies regardless of the redshift of the clusters. The photometric redshifts for the galaxies are estimated with the Directional Neighbourhood Fitting (DNF) algorithm (De Vicente et al. 2016). The galaxy sample used in this paper is

almost identical to the one used in Shin et al. (2019); we refer the readers there for further details.

## 3.2. Method

### 3.2.1. Color Selection in Data

Figure 4 shows how DES galaxies are assigned to different color bins. We measure the density of galaxies in the $(g-r)$–$(r-z)$ color space in redshift bins of $\Delta z = 0.05$. Specifically, we measure the density of galaxies in color space around our cluster sample, within 1 Mpc h⁻¹ of the center, and also around random points on the sky; we then subtract the latter from the former. The plot shows the resultant map of “overdensity” in the color space around our cluster sample, with respect to the global density, in the redshift bin of $z = [0.45, 0.50]$ as an example. One can see an excess of red galaxies and a deficit of blue galaxies around the clusters, due to the quenching of galaxy star formation inside the clusters. This tendency is prevalent in every redshift bin.

We define color bins for the galaxies as follows. First, we identify the “red peak” as the average overdensity-weighted location of the five points with the largest values of overdensity (red point). We identify the “blue peak” similarly, with the five points of the smallest overdensity values (blue point). We finally define the green valley as the location where the absolute value of overdensity is minimized between the red peak and the blue peak (green point). We then draw two lines that pass through (1) the midpoint between the red peak and the green valley and (2) the midpoint between the green valley and the blue peak. These lines are also perpendicular to the line adjoining the red peak and the green valley and the line adjoining the green valley and the blue peak, although the specific choice of the slope of the lines does not affect our color selection significantly because most of the cluster member galaxies are located around the narrow path in the color–color space, as shown in the figure. These two lines then separate the red/green/blue galaxies (three regions divided by the two lines). The black contour in the figure includes 68% of the red-sequence galaxies drawn from the member catalog of the optically selected (RedMaPPer; Rykoff et al. 2014) clusters in the DES, which identifies the clusters by detecting red-sequence overdensities at each redshift. The agreement between the locations of the red-sequence galaxies and the calculated red peak is evident.

### 3.2.2. Modeling Galaxy Number Density Profiles

We model the profiles of galaxies using the fitting formula for galaxy number density profiles described in Diemer & Kravtsov (2014) and used in subsequent work on splashback in different cluster and galaxy samples (More et al. 2016; Baxter et al. 2017; Chang et al. 2018; Shin et al. 2019; Zürcher & More 2019). The 3D number density profile is defined by an inner, virialized Einasto component (Einasto 1965) and an outer, infall term along with a transition region:<fnref n="84" />

$$ \rho(r) = \rho_{\text{inner}}(r)f_{\text{trans}}(r) + \rho_{\text{out}}(r) \quad (1) $$

<footnote n="84">⁸⁴ See Section 3.3 in Diemer & Kravtsov (2014).</footnote>
where

$$ \rho_{\text{inner}}(r) = \rho_s \exp\left(-\frac{2}{\alpha}\left[\left(\frac{r}{r_s}\right)^{\alpha} - 1\right]\right) \quad (2) $$

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

$$ \rho_{\text{outer}}(r) = \rho_0 \left(\frac{r}{r_0}\right)^{-s_e}. \quad (4) $$

The 3D density is integrated to obtain the expected projected density, $\Sigma(R)$, as a function of projected radius, $R$, with the maximum projection length of 40 Mpc $h^{-1}$.

We also account for cluster miscentering. Due to the finite beam size in the CMB survey and other systematic effects, the calculated center of the clusters could differ from the true cluster center of mass.

The azimuthally averaged profile of a cluster miscentered by a distance $R_{\text{mis}}$ is

$$ \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}), \quad (5) $$

where $\Sigma_0$ is the profile without miscentering. We average the profiles over the distribution in $R_{\text{mis}}$ as

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

where $P(R_{\text{mis}})$ is the probability distribution of a cluster to be miscentered by a distance $R_{\text{mis}}$ from the true center, which we model as a Rayleigh distribution (Saro et al. 2015):

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

where $\sigma_R$ characterizes the width of the $R_{\text{mis}}$ distribution. Miscentering is generally measured as a fraction of the FWHM of the ACT beam, which has a fixed angular size; therefore the radius of miscentering is a function of the cluster redshift. In this paper we treat $\sigma_R$ as a free parameter initialized at $\sigma_R = 0.13 \text{ Mpc } h^{-1}$ in radial units, with a wide prior range. We have verified that our choice of prior is a conservative one and is consistent with an angular width of $0.3'$ for S/N = 4 clusters at the mean redshift of the sample (Hilton et al. 2021).

Our model has a total of nine parameters, eight from the halo model and one from the miscentering model. The parameters $r_0$ and $\rho_0$ are degenerate with each other; therefore, we fix $r_0=1.5 \text{ Mpc } h^{-1}$. We fit the projected surface number density profile, described in Section 3.2.3, with the nine-parameter model above using a Markov Chain Monte Carlo (MCMC) method implemented in the `emcee` package (Foreman-Mackey et al. 2013), with priors described in Table 2. In comparison to the priors assumed in Shin et al. (2019), we use (1) the same prior on the Einasto slope parameter $\alpha$ because the mean mass is the same $(3.1 \times 10^{14} \text{ h}^{-1}\text{M}_\odot)$, (2) three times wider prior ranges on the slope parameters of $f_{\text{trans}}(r)$, $(\beta, \gamma)$ to allow more flexibility on the fitting especially for the blue and green galaxies, and (3) three times wider prior on the miscentering parameters.

<table><caption>**Table 2**
Prior Range of Each Model Parameter</caption><thead><tr><th>Parameter</th><th>Prior</th></tr></thead><tbody><tr><td>log ρ<sub>s</sub></td><td>[-∞, ∞]</td></tr><tr><td>log α</td><td>N(log(0.22), 0.6<sup>2</sup>)</td></tr><tr><td>log r<sub>s</sub></td><td>[log(0.01), log(5.0)]</td></tr><tr><td>log r<sub>i</sub></td><td>[log(0.1), log(5.0)]</td></tr><tr><td>log β</td><td>N(log(6.0), 0.6<sup>2</sup>)</td></tr><tr><td>log γ</td><td>N(log(4.0), 0.6<sup>2</sup>)</td></tr><tr><td>log ρ<sub>0</sub></td><td>[-∞, ∞]</td></tr><tr><td>s<sub>e</sub></td><td>[0.1, 10.0]</td></tr><tr><td>ln σ<sub>R</sub></td><td>N(-2.0, 1.2<sup>2</sup>)</td></tr></tbody></table>

**Note.** $N(m, \sigma^2)$ represents a Gaussian prior with mean $m$ and standard deviation $\sigma$. The densities, $\rho_s$ and $\rho_i$, are in units of Mpc$^3 h^{-1}$ and the radii, $r_s$ and $r_0$, are in units of Mpc $h^{-3}$. The miscentering parameter $\sigma_R$ is also in units of Mpc $h^{-1}$.

### 3.2.3. Measurement and MCMC Fitting of Galaxy Profiles

We adopt the same method to measure the galaxy surface number density profile, $\Sigma_g$, as implemented in Shin et al. (2019). We briefly summarize the procedure here and refer readers to Shin et al. (2019) for details.

The galaxy surface density profile, $\Sigma_g$, around the clusters is related to the cluster–galaxy cross-correlation function, $\omega(R)$, as

$$ \Sigma_g = \bar{\Sigma}_g \omega(R), \quad (8) $$

where $\bar{\Sigma}_g$ is the mean galaxy number density. We first divide the clusters into redshift bins of width $dz = 0.025$. In each redshift bin, we measure the cluster–galaxy cross-correlation function, $\omega(\theta, z_i)$, using the Landy–Szalay estimator (Landy & Szalay 1993). The angular cross-correlation function, $\omega(\theta, z_i)$, is then converted to $\omega(R, z_i)$ assuming the midpoint redshift value of the bin. We finally average $\omega(R, z_i)$ over the redshift bins weighted by the number of clusters in each bin to obtain $\omega(R)$. Then we multiply it with the mean galaxy density to generate the final estimate of $\Sigma_g$. We use 15 bins between $0.1 < R < 20 \text{ Mpc } h^{-1}$ spaced equally on a logarithmic scale.

The covariance matrix of the galaxy surface density profile is derived by the jackknife resampling method (Norberg et al. 2009) with 100 patches of similar size. Each jackknife patch retains $\sim 4.4^2 \text{ deg}^2$ of area. The length of each patch corresponds to $\sim 100 \text{ Mpc } h^{-1}$ at $z = 0.5$, significantly larger than the maximum distance scale of interest in this paper.

With the covariance matrix estimated by the jackknife method, $C$, we assume the likelihood ($\mathcal{L}$) to be Gaussian. Given the data, $d$, and the model parameters, $\vec{\theta}$ (Table 2), the likelihood is written as

$$ \ln \mathcal{L}[d|m(\vec{\theta})] = -\frac{1}{2}[d - m(\vec{\theta})]^T C^{-1} [d - m(\vec{\theta})], \quad (9) $$

where $m(\vec{\theta})$ is the model evaluated at the parameter $\vec{\theta}$. The posterior on the model parameters is then expressed as

$$ \ln \mathcal{P}(\vec{\theta}|d) = \ln[\mathcal{L}(d|m(\vec{\theta}))\Pr(\vec{\theta})], \quad (10) $$

where $\Pr(\vec{\theta})$ are the priors imposed on $\vec{\theta}$.
<figure><img src="image_6.png" /><figcaption>**Figure 5.** (Top) Measured surface density profiles of galaxies around clusters from ACT DR5 as a function of 2D projected radius. The shaded regions in red, green, and blue show the $1\sigma$ ranges of the fitted profiles of galaxies of the corresponding color from the MCMC runs (see Section 3.2.3). (Bottom) Fitted 3D slope of the density profile as a function of galaxy color. The vertical line specifies the location of the splashback radius measured from subhalos in the simulation (see Section 3.3.2). The colored segments (around the vertical line) show the estimated splashback radius from the data and its $1\sigma$ uncertainty.</figcaption></figure>

## 3.3. Results

### 3.3.1. Splashback Radius in Clusters from ACT DR5

The top panel of Figure 5 shows the measurement of the projected galaxy number density as a function of radius. Measurement errors are shown as solid curves. The shaded regions show the model fit to data obtained using the method described in the previous section; this region encompasses all curves within $1\sigma$ of the best-fit model parameters. The bottom panel shows the logarithmic slope of the corresponding 3D density profile obtained from the fitted parameters. The errors here are also derived from the $1\sigma$ values of the model parameters.

The black curve and the gray shaded region in the top panel corresponds to the projected number density and the model fit for “all” galaxies around the cluster. The gray shaded region in the bottom panel shows the corresponding, fitted, 3D slope profile. The steepest slope measured from the profile is located at the radius $r_{smin} = 2.4^{+0.3}_{-0.4}$ Mpc h$^{-1}$. The black vertical line denotes the theoretical prediction for the location of the steepest slope measured from the density profile of halos with

<table><caption>**Table 3**
Comparison of the Measurement of the Location of Splashback Radius from All Galaxies with Earlier Work</caption><thead><tr><th>Sample</th><th>[⟨M<sub>200m</sub>⟩, ⟨z⟩]</th><th>r<sub>sp</sub> (Mpc h<sup>-1</sup>)</th><th>Reference</th></tr></thead><tbody><tr><td>SDSS RM</td><td>[1.9, 0.24]</td><td>1.18 ± 0.08</td><td>More et al. (2016)</td></tr><tr><td>DES RM</td><td>[1.8, 0.41]</td><td>1.13 ± 0.07</td><td>Baxter et al. (2017)</td></tr><tr><td>DES RM*</td><td>[1.8, 0.41]</td><td>1.4 ± 0.2</td><td>Chang et al. (2018)</td></tr><tr><td>SPT SZ</td><td>[5.3, 0.49]</td><td>2.4<sup>+0.5</sup><sub>-0.5</sub></td><td>Shin et al. (2019)</td></tr><tr><td>ACT SZ</td><td>[5.8, 0.49]</td><td>2.2<sup>+0.7</sup><sub>-0.6</sub></td><td>Shin et al. (2019)</td></tr><tr><td>Planck SZ</td><td>[6.2, 0.177]</td><td>1.9<sup>+0.4</sup><sub>-0.3</sub></td><td>Zürcher & More (2019)</td></tr><tr><td>ACT DR5 SZ</td><td>[5.8, 0.49]</td><td>2.4<sup>+0.3</sup><sub>-0.4</sub></td><td>This work</td></tr></tbody></table>

**Note.** RM corresponds to the RedMaPPer cluster sample and SZ to the Sunyaev–Ze’ldovich samples. Each row designates the splashback radius measured using the galaxy number density profiles except row 3, DES RM*, which specifies the splashback measured using weak lensing. The mean mass of the sample is quoted in units of $10^{14} M_{\odot} h^{-1}$.

$v_{peak} > 150$ km s$^{-1}$ around matched clusters in the MDPL2 simulation. The $v_{peak}$ threshold has been derived by matching the magnitude of the number density profile at large radius (More et al. 2016).

When all halos are considered without splitting on accretion time, the transition from the infall region to the multistreaming region occurs at the traditional splashback radius. Therefore, the location of the steepest slope measured using all galaxies around the clusters corresponds to the splashback radius, $r_{sp}$, for the sample. We find that the observed value is statistically consistent with theoretical predictions from simulations. We also confirm, in agreement with previous results from Shin et al. (2019) and Zürcher & More (2019), that the splashback radius measured in SZ clusters does not show any significant discrepancy between theory and observations (see the black line and gray shaded curve in Figure 5). This is unlike the splashback radius measured using optically selected RedMaPPer clusters (More et al. 2016; Baxter et al. 2017; Chang et al. 2018; Shin et al. 2019), where the measured value is smaller than that expected from theory. A summary of the measurements of the splashback radius from different cluster samples with a range of masses is shown in Table 3.

### 3.3.2. Galaxy Density and Slope Profile as a Function of Galaxy Color

In this section we study the distribution of galaxies of different colors around the clusters from the ACT DR5 catalog. The red, green, and blue curves in Figure 5 show the measurements and their best-fit model curves for galaxies assigned to the corresponding color bins using the method described in Section 3.2.1.

We find that the overall shape of the projected number density profiles, and consequently the slope profiles, is different for the three populations. The apparent location of the splashback radius, traditionally defined as the location of the minimum of the slope profile, appears to show shifts as a function of galaxy color. The locations of the slope minimum, $r_{smin}$, for the red, green, and blue galaxies are constrained to $2.2^{+0.2}_{-0.1}$, $1.9^{+0.2}_{-0.2}$, and $0.95^{+0.2}_{-0.1}$ Mpc h$^{-1}$, respectively. The slope minimum for green galaxies is shallower and at a smaller radius compared to red ones, but the difference is not statistically significant; however, the blue galaxies show a weak feature, with the minimum in the slope located at a smaller radius
compared to both green and red galaxies. We also note that the red galaxies show a significantly steeper inner profile compared to green and blue galaxies.

Previous work has also measured the slope of the density profile as a function of galaxy color. Baxter et al. (2017) measured the splashback radius in the bluest and reddest quartiles of the $(g-r)$ galaxy colors for RedMaPPer clusters and found that the splashback feature was more prominent for galaxies in the red quartile and slightly larger than the full galaxy sample, while the galaxies in the blue quartile show a weak splashback-like feature. Shin et al. (2019) developed an improved color selection and measured the splashback locations of galaxies of different colors around SZ-selected clusters from the South Pole Telescope (SPT) and ACT. Shin et al. (2019) looked at not only galaxies in the blue and red quartiles but also galaxies that lie in the green valley. While the bluest galaxies did not show a significant minimum in the transition region, the green and red galaxies both show a distinct splashback-like minimum feature. Due to the small number of clusters (300 for the SPT sample and about 100 for the ACT sample), small shifts of the minimum between galaxies of different colors, if they exist, could not be detected. In the current sample of ACT DR5 clusters, with three times the number of clusters compared to the SPT sample, we can study more subtle changes in the density profiles with color.

Considering that galaxies evolve in clusters due to astrophysical processes, with their SFRs decreasing with time, we expect galaxies to evolve from being blue and star-forming to red and quenched over the course of their orbits, migrating between the different color bins. The density profiles that we measure here are a snapshot in time of this process. In Section 2 we studied the signatures of the net time spent within a cluster on the density profile of a population of halos. In particular we found that the location of the minimum of the slope of the density profile traces a phase-space discontinuity and encodes information about how long a population of subhalos has been inside a halo. Therefore, one possible explanation for the movement of the location of the steepest slope for galaxies of different colors observed in data can be the movement of the phase-space boundary for the different populations as they evolve through the phase space of the halo. In the next section, we elaborate on the possible connection between galaxy colors and infall time and attempt to model the process of galaxy quenching within clusters by mapping galaxies to subhalos in N-body simulations.

# 4. Galaxy Quenching in Clusters

## 4.1. Infall Time of Observed Galaxies

If intracluster processes are indeed responsible for quenching the SFRs of galaxies that fall into clusters, we may expect that an infalling population of galaxies will evolve to a redder population over time. Considering that quenching is a continuous process, galaxies that are originally blue will convert to green galaxies and eventually to red ones, while green galaxies will turn red with time.<fnref n="85" /> This implies, first, that the number of galaxies as a function of color is not conserved at each point in the phase space of clusters. Second, if the quenching timescales are short enough, blue galaxies will not exist in the phase-space locations populated by galaxies that

have been within the cluster over multiple orbits, therefore affecting the density profiles in ways discussed in Section 2. For example, if all blue galaxies are quenched before pericentric passage or at pericenter, they will only exist in the single infall stream and will not show a minimum in their logarithmic slopes. Similarly, if all green galaxies convert to red ones before they reach the apocenter of their first orbits (i.e., splashback), $r_{smin}$ will be smaller compared to $r_{sp}$. In this scenario, the galaxies that reach splashback are only red and the steepest slope in their profile $r_{smin,red} \equiv r_{sp}$.

We compare the $r_{smin}$ observed for different color bins in the observed DES galaxy population with $r_{smin}$ for halo populations in simulations that have been accreted onto their host clusters at different times on average. These comparisons are shown in Figure 2. The middle panel of the figure demonstrates that different values of $r_{smin}$ correspond to populations of halos that have a different maximum infall time threshold, $t_{max}$. Each $r_{smin}$ is derived from the logarithmic slope of the number density profile of halos by varying the quantity $t_{max}$, such that each population comprises all halos around clusters with $t_{in} < t_{max}$ (top panel). The horizontal width of the red, green, and blue shaded regions in the middle panel of Figure 2 covers the $1\sigma$ region of $r_{smin}$ for the galaxies with different colors measured from data. Comparing with the simulations, we find that red galaxies in data show an $r_{smin}$ consistent with halos in the simulation that have been inside the halo for longer than 3.2 Gyr. The $r_{smin}$ for green galaxies corresponds to the $r_{smin}$ for halos in simulations that have been in the cluster for at least $t_{in} > 2.2$ Gyr; they most likely do not survive beyond 3.2 Gyr. The blue galaxies are consistent with a population of objects that have been inside the cluster for less than 1.5 Gyr. The bottom panel of the figure shows the simulation curves for the slope profiles for halo populations that have their $r_{smin}$ in the $1\sigma$ region of $r_{smin}$ for observed galaxies.

We conclude that the shift in the minimum of the slope of the density profile for observed galaxies can be inferred to arise from the phenomenon that different populations of galaxies survive in the halo for different net amounts of time. We emphasize that this is not in fact a movement of the “splashback radius” or the boundary of the halo, traced by the apocenter of the first orbit for different galaxy colors but instead a movement of the phase-space boundary for populations of different colors. In other words while the “splashback” radius traces the boundary of the whole multistreaming region of the halo, different color galaxies, by virtue of the fact that galaxies change color over their orbits, can have their streams end at different locations in phase space.

## 4.2. Modeling Galaxy Quenching

In this section we model intracluster galaxy quenching to obtain the relevant timescales involved in star formation quenching within clusters. We adopt the commonly used exponential quenching model (Papovich et al. 2001; Shapley et al. 2005; Lee et al. 2010; Förster Schreiber et al. 2009); the details of our model follow Wetzel et al. (2013). In this model the quenching of star formation in a galaxy falling into a cluster is described by two relevant timescales, a delay time, $t_d$, which is the duration of time that the galaxy SFR remains unaffected by intracluster processes, and an exponential decay timescale, $t_q$, after the period of delay has passed. The outskirts of clusters have low densities of both gas and dark matter, it can be conjectured therefore that infalling galaxies start quenching

<footnote n="85">⁸⁵ We reiterate that by green galaxies we mean galaxies that are in the green valley, between red galaxies with low SFRs and blue ones with high SFRs.</footnote>
strongly only when they reach the central regions of the host; the delay time, $t_d$, can therefore be related to the pericenter-crossing time of the galaxy. While $t_d$ and $t_q$ capture star formation quenching within the cluster, we note that galaxies also have an intrinsic decay time for their SFR in the field (Noeske et al. 2007), i.e., isolated or central galaxies (galaxies in the field) eventually stop forming stars with time. We call this intrinsic decay timescale, $t_{q,iso}$, the quenching timescale for isolated (field) galaxies; we assume that this timescale applies to all galaxies before they are impacted by cluster processes.

Adopting the model from Wetzel et al. (2013), the SFR of a satellite galaxy that falls into the cluster can be defined as

$$SFR_{\text{sat}}(t) = \begin{cases} SFR_{\text{iso}}(t), & \text{if } t < t_d \\ SFR_{\text{iso}}(t) \exp\left(-\frac{t-t_d}{t_q}\right), & \text{if } t > t_d, \end{cases} \quad (11)$$

where $SFR_{\text{sat}}(t)$ represents the SFR of galaxies that are satellites of clusters at time $t$, $SFR_{\text{iso}}(t)$ is the SFR of isolated galaxies, $t$ is the time since infall of the galaxy into the cluster virial radius.

The intrinsic evolution of SFR of an isolated galaxy is modeled as

$$SFR_{\text{iso}}(t) = SFR(t_i) \exp(-t/t_{q,\text{iso}}), \quad (12)$$

where $t_{q,iso}$ is the intrinsic timescale for the decay of star formation, independent of the cluster environment, and $SFR(t_i)$ is the initial SFR when the galaxies are formed at time $t_i$.

We obtain constraints on our model for quenching in two separate ways. First, we use only the information about the location of the splashback radius, and second, we use the profile of the ratio of galaxy number densities of different colors. In the following sections, we describe our method.

### 4.3. Constraints on Galaxy Quenching Timescales

Our goal is to constrain the galaxy star formation evolution timescales, $t_d$ and $t_q$, in our model. The time $t_d$ is the delay before quenching begins inside a cluster, and $t_q$ is the exponential decay rate of star formation after the onset of quenching. To constrain the quenching timescales, we make mock galaxy catalogs from the MDPL2 simulations by assigning galaxies to halos around cluster-mass objects and evolve every mock galaxy in our sample with the model described above to find the parameters that best describe our data. We match the redshift distribution and mean mass of the simulation cluster sample to the observed sample of SZ clusters. We reiterate that we refer to the main host cluster dark matter halo as a “cluster” and study the entire halo field around it, including both subhalos and halos outside the virial radius of the cluster.

First, we select a minimum $v_{\text{peak}}$ threshold for the halos in the simulation to correspond to the magnitude limit of the observed galaxies. We use abundance matching to determine this threshold, i.e., the cumulative distribution of observed galaxy magnitudes is matched to the $v_{\text{peak}}$ distribution of halos in simulation.<fnref n="86" /> We find that our sample of galaxies with $M_r < -19.87$ matches with a halo $v_{\text{peak}}$ threshold of 150 km s⁻¹. We extract all subhalos with $v_{\text{peak}} > 150$ km s⁻¹.

in a spherical volume with radius $10 \text{ Mpc h}^{-1}$ around each cluster. Further, we assign SFRs to the simulation halos based on the quenching model described in Equation (11). The sample of clusters and the halos in their neighborhood is the same as the one used in previous sections on simulations (Section 2.3), extended out to $10 \text{ Mpc h}^{-1}$.

#### 4.3.1. Initial Distribution of SFR

To constrain the quenching timescales of galaxies within the cluster halos, it is essential to correctly model the initial distribution of SFRs before they were accreted and became satellites. To model the distribution of the SFR of field/central galaxies (i.e., galaxies that have not fallen into larger host halos) we use the Universe Machine simulations (Behroozi et al. 2019). Universe Machine is an empirical model that populates CDM N-body simulations with galaxies using an extensive set of observational constraints. It provides galaxy properties like stellar masses and SFRs that have been assigned based on a detailed parameterization of the galaxy-halo connection based on the growth history of each halo.

In general, the SFR distribution in the universe at a given redshift is a bimodal function, with its two peaks at high and low SFRs corresponding to star-forming and quiescent galaxies. The fraction of quenched or star-forming galaxies, in turn, is a function of the galaxy stellar mass or alternatively the virial mass of its host dark matter halo. Initial SFRs are therefore assigned to halos based on redshift and maximum circular velocity, $v_{\text{cmax}}$, which is a proxy for halo mass. SFRs are assigned to subhalos of the cluster based on their $v_{\text{acc}}$, i.e., the $v_{\text{cmax}}$ at the time of accretion, and from the distribution of the SFR of central galaxies at the accretion redshift, $z_{\text{acc}}$. In other words we start the clock for subhalos at $z_{\text{acc}}$ and draw their initial SFRs from the distribution of colors for isolated galaxies at that redshift, inferred from Universe Machine. In detail, we model the distribution of the SFR of galaxies as a double Gaussian in bins of $v_{\text{cmax}}$ for each redshift,

$$f(SFR_i) = c_1 G(\mu_1, \sigma_1) + c_2 G(\mu_2, \sigma_2), \quad (13)$$

where $G(\mu, \sigma)$ is a Gaussian with mean $\mu$ and width $\sigma$, and $c_1$ and $c_2$ are parameters that control the relative fraction of galaxies in the two Gaussians. The variables $\mu$, $\sigma$, and $c$ all depend on redshift and the $v_{\text{cmax}}$ of the halo. We use Universe Machine to calibrate these parameters for both central and satellite galaxies at all redshift snapshots between $z=0-4$ available in Universe Machine in three bins of $v_{\text{cmax}}$: $[150 < v_{\text{cmax}} < 170, 170 < v_{\text{cmax}} < 200, v_{\text{cmax}} > 200]$ in units of km s⁻¹. Subhalos of the cluster are assigned initial SFRs by randomly drawing from the Gaussian distribution of SFRs described by Equation (13) at $z_{\text{acc}}$ for the corresponding $v_{\text{acc}}$ bin, where the parameters of the Gaussian are calibrated using central galaxies from Universe Machine.

SFRs are also assigned to all halos outside the virial radius of the cluster out to $10 \text{ Mpc h}^{-1}$ based on their $v_{\text{cmax}}$ or $v_{\text{acc}}$ (if they are subhalos).<fnref n="87" />

Once the initial distributions of SFRs are assigned, with galaxies within subhalos of the cluster set to initial states before

<footnote n="86">⁸⁶ We use peak quantities for subhalos under the assumption that the galaxy stellar mass or luminosity traces the mass of the original unstripped subhalo before it falls into the cluster’s tidal field.</footnote>

<footnote n="87">⁸⁷ We find that the region around clusters can have a large number of galaxies that are satellites within larger halos, these satellites tend to have lower SFRs than centrals, and their distribution should therefore be drawn from satellite SFR distributions.</footnote>
infall, we evolve the SFRs of all halos based on the quenching model, Equation (11).

We summarize our main steps here:

1. Galaxies are assigned to all $v_{\text{peak}} > 150 \text{ km s}^{-1}$ halos around cluster-mass halos.

2. All galaxies are assigned initial SFR based on their $v_{\text{cmax}}$ and redshift by drawing from distributions of $v_{\text{cmax}}$-SFR (z) calibrated using Universe Machine.

3. Galaxies in subhalos are assigned SFRs based on the $v_{\text{cmax}}$ at their accretion redshift, $z_{\text{acc}}$.

4. Galaxies are evolved with the quenching model from their initial SFRs.

In the following sections, we describe the constraints on the quenching model obtained from the location of the steepest slope $r_{\text{smin}}$ and the number density profiles of galaxies.

### 4.3.2. Constraints from Slope Minimum, $r_{\text{smin}}$

Here we describe the constraints we obtain on the galaxy quenching model from the location of the slope minimum, $r_{\text{smin}}$, for the three populations of galaxies with different colors. The free parameters of our model are the delay time, $t_d$, and the quenching timescale, $t_q$. We divide our parameter space into 30 bins in delay time between $0.1 < t_d < 2$ Gyr and 100 bins in quenching time between $0 < t_q < 5$ Gyr. We evaluate the location of the steepest slopes for mock galaxy populations of different colors at each point in the parameter space and compare it to data.

The procedure in detail is as follows. Using the time since infall, $t$, provided by the ROCKSTAR catalog, we evolve each mock galaxy's initial SFR to the current redshift according to Equation (11), for every pair of $t_q$ and $t_d$ in the parameter space. During the initial period of delay ($t_d$), before the onset of exponential quenching, they are evolved according to Equation (12), following the evolution of isolated galaxies calibrated using Universe Machine. The evolution timescale, $t_{q,iso}$, varies linearly between 1 and 3 Gyr for the logarithmic velocity range of the subhalos used in this paper. Once the galaxy SFRs have been evolved to the current time (the redshift at which their host cluster is observed), we divide them into red, green, and blue based on the new SFR. The splits in the SFR space to assign simulation galaxies red, green, or blue colors are adjusted to match the color fraction profile outside the splashback radius by fitting the color fraction profiles from the radial bin at $r = 2.64 \text{ Mpc h}^{-1}$ outwards (see Section 4.3.3 for a detailed description). Following the color split, we compute the slope profile for mock galaxy populations in each color bin and find the location of the steepest slope, $r_{\text{smin}}$. We compare the location of the steepest slope at each point in the plane of $t_d$ and $t_q$ to the values obtained from the data.

The constraints on the model from the location of $r_{\text{smin}}$ are shown in Figure 6. The gray contours in the top panel show the allowed region of the parameter space constrained from the ACT DR5 cluster sample. The dark gray region corresponds to the 68% confidence interval, which represents the statistical uncertainty obtained using simulations as described above. The likelihood is calculated assuming a Gaussian form. As expected, there is a degeneracy between the delay time and the quenching timescale. The direction of degeneracy is denoted by the black dashed line in the top panel; short delay times with long $t_q$ and long delay times with short $t_q$ are both allowed by the data. By extrapolating the black line to $t_d = 0$

<figure><img src="image_7.png" /><figcaption>**Figure 6.** (Top) Constraints on the quenching model parameters from $r_{\text{smin}}$. The darker and the lighter gray contours represent the $1\sigma$ and the $2\sigma$ confidence regions, respectively (see Section 4.3.2). The black dashed line is the best-constrained linear combination of the timescales. (Bottom) Constraints on the best-constrained linear combination of $t_q$ and $t_d$ from $r_{\text{smin}}$. The blue shaded region shows the $1\sigma$ constraint in the linear combination, $1.85_{-0.15}^{+0.12}$ Gyr.</figcaption></figure>

Gyr to obtain its intercept on the y-axis, it appears that a galaxy takes at least $1.15 \pm 0.3$ Gyr to transition from a star-forming phase to a quenched phase in our model, in the sense that the sum of $t_q$ and $t_d$ is always greater than $\sim 1.15$ Gyr.

We find that the location of $r_{\text{smin}}$ can be used to best constrain the linear combination $t_d + 1.65 t_q$ (bottom panel Figure 6). The best-fit value for this combination is $1.85_{-0.15}^{+0.12}$ Gyr. Note that $r_{\text{smin}}$ is sensitive to the total time that a galaxy population spends inside the cluster, particularly before apocenter crossing; therefore, it constrains the linear combination comprised of the sum of the total delay time and a multiple of the e-folding time of quenching $t_q$, considering the number galaxies with a given SFR reduces by $\sim 65\%$ after first e-folding.
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Color fraction—the measured profiles of the fraction of red/green/ blue galaxies (solid lines with error bars) as a function of three-dimensional clustercentric radius, as defined in Section 4. The shaded regions in the red, green, and blue colors are the 1σ confidence intervals from fitting the quenching model to the color ratio. The χ² per degree of freedom of the fit is 1.46. The inner (gray shaded) and outer cluster regions are fitted separately. The vertical lines denote the 1σ interval of the measured splashback radius for all galaxies in the SZ cluster sample.</figcaption></figure>

While the location of the minimum gives us an intuitive understanding of the phase-space picture of galaxies that is easy to interpret, there is, theoretically, more information contained in the entire density profile measured from this rich data set. Referring to Figure 1, $r_{smin}$ gives us the location of the discontinuity between multistream and single stream in phase space, therefore it can tell us about the maximum age of a galaxy population within a halo. On the other hand, the density profile is the integration of the number of galaxies in different streams at each radius; its shape therefore contains information about the radial location of the transition from delay to quenching phases, helping to break the degeneracy between the two timescales. In the next section we use the color ratios to constrain the quenching parameters.

### 4.3.3. Constraints from the Color Fraction Profile

To break the degeneracy between $t_q$ and $t_d$, we explore constraints from the complete number density profile of galaxies in this section. In particular we use the color fraction profile, the ratio of the number of galaxies of a particular color at a given radius to the total number of galaxies at that radius. The solid lines with error bars in Figure 7 show the fraction of galaxies in each bin of color around the ACT DR5 cluster sample as a function of 3D clustercentric radius. The 3D color fraction profiles and their covariance matrices are retrieved from the MCMC chains of the model fitting described in Section 3.3.2.

We use the best-fit color fraction curves with 1σ error bars as input data to constrain the quenching model. We fit our quenching model to nine radial bins between a radius of 0.6 and 10 Mpc $h^{-1}$. Note that we do not use any radial bins

smaller than 0.6 Mpc $h^{-1}$ as the subhalo distributions in the simulations are known to deviate significantly from observed galaxy distribution in this region. In particular, subhalos can be stripped to masses below the halo-finder resolution limit due to tidal stripping or artificial disruption (e.g., van den Bosch et al. 2018), an effect that is not necessarily reflected in the galaxies. Moreover, blending near cluster centers may lead to an underestimation of the number of galaxies in the innermost bins (Zhang et al. 2015). This effect is particularly important for our analysis as blue galaxies are systematically fainter than red ones (Bell et al. 2004; Strateva et al. 2001), creating a possible systematic bias in the color ratio due to incompleteness near the center. Although we expect our absolute magnitude cut to mitigate some of this effect, it provides us with additional motivation to use only the radial bins larger than 0.6 Mpc $h^{-1}$. We do not expect systematic effects for the innermost galaxies to affect the slope or ratio measurements in the outskirts.

We constrain our quenching model by comparing the measured fraction of galaxy colors (Figure 7) to that calculated around the clusters in the MDPL2 simulations using the halo number density profiles around them. The color fraction profiles shown in Figure 7 span a large range that encompasses the virialized region and the region where halos are still falling into the cluster. Bearing in mind that the infalling halos outside the clusters could themselves be part of larger structures like group-mass halos whose quenching timescales can, in principle, differ from the timescales within clusters, we treat the inner and outer regions around the cluster separately. We denote the radial bins outside the splashback boundary of the cluster sample with $r \ge 2.64$ Mpc $h^{-1}$ as the “outer” region of the cluster and the radial bins within 2.64 Mpc $h^{-1}$ as the inner region. We allow the quenching timescale of satellites outside the cluster to be different from that inside it.

For clarity, we discuss our results from the outer and inner regions separately below.

**Constraints from the outer regions of clusters:** To map from the observed color space to the SFR space in simulations, we only use the color fraction in the radial bins outside the clusters. We define our boundary between red, green, and blue galaxies in SFR space based on the relative density of red, green, and blue galaxies between 2.64–10 Mpc $h^{-1}$, as this region does not evolve with the quenching timescales within the cluster.

We use two parameters, $f_{rg}$ and $f_{gb}$, to specify the splits in the SFR space to assign simulation halos to red, green, and blue color bins. The value of SFR that separates red from green and green from blue are described as follows:

$$ 
\begin{aligned}
SFR_{rg} &= f_{rg} SFR_r + (1 - f_{rg}) SFR_g \\
SFR_{gb} &= f_{gb} SFR_g + (1 - f_{gb}) SFR_g,
\end{aligned}
\quad (14) $$

where $SFR_r$ and $SFR_g$ correspond to the location of the two peaks at low and high SFR in the bimodal SFR distribution, respectively. $SFR_g$ is the location of the minimum separating the two peaks, often referred to as the green valley. The splitting points are described with respect to peaks in each individual $v_{cmax}$ bin.

Reviewing Section 4.3.1, the SFRs of central halos at $r > 2.64$ Mpc $h^{-1}$ are assigned based on their current $v_{cmax}$. However, subhalos in this region that are within the virial radius of more massive halos have an accretion time associated with them that corresponds to the redshift at which they crossed
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Constraints on quenching timescales from 3D color fraction profiles. (Red shaded) The $1\sigma$ and the $2\sigma$ contours from fitting the color fraction profile between 0.6 and $2.64 \text{ Mpc h}^{-1}$ ("inner" profile). The red shaded region corresponds to quenching timescales within the clusters from ACT DR5. (Blue lines) The constraints on parameters using the color fraction profile of the "outer" region beyond $2.64 \text{ Mpc h}^{-1}$. The darker and the lighter blue contours correspond to the 68% and 95% confidence regions, respectively. These correspond to quenching timescales within halos outside the virial radius of the cluster.</figcaption></figure>

into the virial radius of their current host. We assign initial SFRs to these objects based on the $v_{cmax}$ at the time of accretion onto their host halos and use the quenching model to evolve them to the observed redshift. This allows us to infer quenching timescales associated with groups outside and around the ACT DR5 clusters.

The profile in the outer region therefore allows us to fix the color-SFR mapping and additionally constrain any quenching timescales associated with groups outside the cluster. We refer to the quenching parameters in the outer region as $t_{d,out}$ and $t_{q,out}$.

We fit the outer color fraction profile by running MCMC chains assuming a Gaussian likelihood (see Section 3.2.3). We vary the four model parameters, $f_{rg}$, $f_{gb}$, $t_{q,out}$, and $t_{d,out}$; split the galaxy SFRs in simulations into red, green, and blue based on $f_{rg}$ and $f_{gb}$; and compute the color fraction profile as a function of the clustercentric radius. The color fraction profiles in the simulation are then compared to those measured from data (Figure 7) to derive the constraints on the model parameters.

The blue contours in Figure 8 show the resultant constraints on the quenching parameters from the outer profiles. The 68% confidence interval of the delay time, $t_{d,out}$, is 0.8–2.7 Gyr and $t_{q,out} < 0.7$ Gyr at 95% confidence level (blue contours). We emphasize that these quenching parameters correspond to timescales within subhalo hosts outside the cluster, in a sense, they correspond to preprocessing timescales within groups before infall. The typical mass of a subhalo-host in the infall region is $10^{12.5} \text{ M}_\odot \text{ h}^{-1}$. This result is consistent with previous studies for quenching timescales around group-sized halos, including Wetzel et al. (2013), who report a delay time of 2–4 Gyr and an upper limit for the quenching timescale of 0.8 Gyr.

**Constraints from the inner profiles of clusters:** The analysis of the inner region is performed in a similar manner as the outer region. We adopt the splitting locations $f_{rg}$ and $f_{gb}$, constrained by the outer region as priors for the color splits. As before, at every step in the MCMC chain, we vary four parameters, the two splitting locations, and the quenching parameters in the inner cluster, $t_{q,in}$ and $t_{d,in}$, where "in" denotes the inner region.

The final constraint on the quenching parameters for the interior of the cluster is shown by the red contours in Figure 8. The best-fit value for the exponential quenching timescale, $t_{q,in}$, is $0.6_{-0.1}^{+0.1}$ and for the delay time, $t_{d,in}$, is $1.0_{-0.3}^{+0.1}$. While the quenching timescale, $t_{q,in}$, is consistent with that in Wetzel et al. (2013), we find a significantly smaller delay time, $t_{d,in}$, than Wetzel et al. (2013). It appears that the onset of exponential quenching happens at an earlier time after infall in the massive clusters that we observe compared to the lower-mass, group-sized objects explored in Wetzel et al. (2013). This can be attributed to the fact that the density of gas and dark matter in lower-mass objects can be significantly smaller compared to massive clusters used in this study, allowing them to survive over longer times and more orbits before quenching begins. We also note that our sample is at a higher mean redshift compared to Wetzel et al. (2013). The quenching timescale of satellite galaxies can in principle evolve with redshift, we compare our best-fit results for the sum of the quenching timescales to previous work that have studied this evolution in Figure 11. The gray band shows the expected decay of dynamical friction timescales of satellites within halos; we find that our results are consistent with the overall trend of decaying quenching timescales at higher redshifts.

In Figure 9 we illustrate the likely distribution of the galaxies in the phase space of dark matter halos using the best-fit quenching parameters. We generate the distribution of galaxies in the $r-v_{rad}$ plane using the simulated galaxies and separating them in color space. The color indicates the fraction of galaxies in each sample that are in a given pixel. Our quenching parameters indicate that the blue galaxies are primarily on first infall, having only recently crossed pericenter, whereas the red and green galaxies are in the virialized region of the cluster. We note that previous work like Adhikari et al. (2019), Oman & Hudson (2016), and Orsi & Angulo (2018) that have used spectroscopic data to infer dynamical properties of galaxies also conclude, similar to our findings, that blue galaxies have been accreted onto their host clusters recently and live mostly in the infall stream. Our method, combined with the inference drawn from the location $r_{smin}$, further implies that they are past pericenter in their orbits. We do not compare the values of the timescales in these works here due to different cluster and color selection choices made in the papers mentioned above.

### 4.3.4. Discussion on Quenching Constraints

Figure 10 shows the comparison between the constraints obtained on the quenching timescales within clusters from the two different methods. The gray and red shaded regions show the constraints from $r_{smin}$ and the 3D color ratio profiles, respectively. The two results are consistent with each other statistically. The color fraction profiles give relatively tighter constraints due to the greater amount of information contained in them. While the relative locations of $r_{smin}$ allow for longer quenching times and shorter delays, the color fraction profile breaks the degeneracy to a large extent, preferring a delay time that is close to 1 Gyr. Comparing with subhalo accretion scales
<figure><img src="image_10.png" /><figcaption>**Figure 9.** The inferred distribution of galaxies of different colors in phase space from the best-fit parameters of the quenching timescales obtained from data. The phase-space diagram has been generated using simulations halos that have been assigned colors based on the quenching model. The three panels correspond to red, green, and blue galaxies from left to right. The red galaxies appear to be mostly concentrated in the virialized, multistreaming region while the blue galaxies are still infalling. The green galaxy population is a combination of infalling galaxies and galaxies that have crossed over to a second stream in phase space but not reached apocenter. (see also Figure 1).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 10.** Comparison between constraints on $t_q$ and $t_d$ from two techniques: the location of $r_{smin}$ and the 3D color fraction profiles. The gray bands correspond to the $1\sigma$ and $2\sigma$ constraints from $r_{smin}$. The red shaded contour corresponds to constraints obtained using the entire inner profile (same as Figure 8).</figcaption></figure>

from simulations, this result implies that the onset of exponential quenching begins close to the pericenter of a cluster, as it takes a subhalo approximately a gigayear to reach pericentric passage from the virial radius of a cluster.

It is important also to consider that a large fraction of the current subhalos of a cluster, enter as parts of bound groups, i.e., they are often subhalos of groups even before infall into the observed cluster. The fraction of subhalos that enter as parts of groups can be as large as 50% in cluster-mass objects (Wetzel et al. 2013). As the quenching timescales in this work are calculated with respect to the “first accretion time” provided by the ROCKSTAR catalog, i.e., the time when a halo crosses the virial radius of another halo, more massive than itself, for the first time, we note that some of the subhalos in our current cluster have actually been partly preprocessed before infall. We

defer a detailed study of the two separate populations of subhalos to future study.

With regard to our constraints from the minimum of the slope, it has been known that the location of the splashback radius can be affected by various halo properties like the mass and the halo accretion rate (Adhikari et al. 2014; Diemer & Kravtsov 2014). However, we do not expect that the shift in the slope minimum is being caused solely by galaxy colors selecting significantly different types of halos. For example, we can rule out a large scatter in halo mass given the fact that red galaxies show a relatively narrow and steep feature at splashback. Another possibility to consider is that faster accreting halos in our sample have a higher fraction of blue galaxies and therefore appear to have a smaller splashback radius. However, the splashback radius for fast accreting systems are usually deeper and the trend with accretion rate and depth does not hold with our different color samples. Nonetheless, it will be useful to constrain the accretion history distribution of our sample of clusters in the future and incorporate it in our modeling.

We note some caveats that must be considered with regard to the methods used in this paper. Our simple quenching model with two free parameters does not include all modeling uncertainties that may be present in the galaxy-halo connection. For example, we do not include scatter in the abundance matching relationship, and we also do not include scatter in the relation between the SFR and $v_{cmax}$; both of these can broaden the error contours shown in Figure 5. Another theoretical uncertainty not included in this work is the modeling of orphan galaxies in simulations; these are galaxies that are expected to live within destroyed subhalos in the simulation. Subhalos suffer enhanced disruption compared to galaxies as they are more extended and feel stronger tidal forces, and also suffer from artificial disruption due to resolution effects; not accounting for galaxies within such subhalos can significantly change the cluster profiles. Therefore, we caution against a direct comparison of the slope and number density profile from the innermost region of clusters in simulations to the slopes observed in galaxies of real clusters. Both of these effects are stronger near cluster centers and we have mitigated them by considering galaxy profiles outside the very center of the cluster by using only the region, $r > 0.6$ Mpc $h^{-1}$ in our analysis; we also use relative color fractions for comparison
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Comparison between the total quenching timescale inferred in this paper at $z = 0.484$ with previous work based on Tinker et al. (2013). The blue point shows the best-fit sum of the quenching and delay time from our color fraction analysis and its $2\sigma$ error bar. The gray band corresponds to a theoretical estimate of the relationship between the dynamical friction timescale within a halo as a function of redshift calibrated using the quenching timescale obtained from low-redshift SDSS groups in Wetzel et al. (2013). The yellow data point was obtained from Tinker et al. (2010) by analyzing the clustering of high-redshift red galaxies, and the green points have been obtained from clustering studies in Tinker & Wetzel (2010); both show $2\sigma$ error bars.</figcaption></figure>

instead of the absolute values of the number density at each radius. We have also confirmed that our basic physical picture —that of the phase-space discontinuity tracing the infall time of a population—is robust, directly using particles in simulation (Section 2.2), which do not suffer from artifacts like subhalo disruption. We also note that while the difference in the maximum time spent within the cluster by different galaxy populations can explain the relative location of splashback between galaxies with different colors, it does not necessarily explain the complete shape of the density profile. A comparison between Figures 2 and 5 clearly demonstrates this point. While the location of splashback for green galaxies agrees with simulations, these galaxies show significantly shallower slopes in the inner region in observations compared to simulations.

Furthermore, in this work we do not explicitly account for differences in quenching timescales between galaxies that are on different types of orbits. In principle, galaxies on radial orbits that form the inner regions of the cluster can have different quenching timescales from those that are on tangential orbits, as they probe different cluster environments. We note that while our model does not explicitly include orbital parameters, infall times themselves correlate with pericenter distances. We defer a detailed quantitative analysis of this effect to future studies. With regard to the location of $r_{smin}$, we note that while the radial profile can be significantly different for galaxies in tangential orbits compared to those on radial orbits, we find that the splashback boundary for these populations does not move significantly; therefore, differences in orbital angular momentum cannot alone explain the behavior of $r_{smin}$ observed in data (see the Appendix).

In this paper we have used CDM-only simulations to map the distribution of observed galaxies around clusters to satellites in dark matter halos. However, in principle a more detailed study of the specific effects of ram pressure or harassment that lead to quenching can be conducted using controlled hydrodynamical simulations. We note that some of the issues like the disruption of satellites near cluster centers are less important in hydrodynamical simulations as they simulate galaxies directly; however, the exact timescales involved in star formation evolution can have significant uncertainties and depend strongly on the modeling of feedback mechanisms. The timescales obtained from this analysis appear to be consistent with that found in hydrodynamical simulations for massive clusters in Bahe & McCarthy (2015); however, the delay time in our work is smaller compared to that in Rhee et al. (2020). We note that the cluster sample in our work is significantly more massive and at a higher mean redshift compared to both the studies quoted above, and a detailed comparison with hydrodynamical simulations with the mass and redshift distribution matched to the cluster sample is a useful direction for future work.

## 5. Conclusion

We study the distribution of galaxies around SZ-selected galaxy clusters as a function of galaxy color. The galaxy sample is derived from the DES Y3 gold galaxy catalog, and the SZ-selected clusters are identified by AdvACT and will be published as a part of ACT DR5 Hilton et al. (2021). The number density profiles encode important information about the dynamical history of galaxies in these massive systems; in this paper, we study them in the context of the splashback radius $r_{sp}$. We motivate the use of a more general parameter $r_{smin}$ along with $r_{sp}$, as it represents the location of the steepest slope for galaxies that have not reached splashback. While the slope minimum itself has been explored as the splashback radius before in numerous studies, we point out here that it does not necessarily correspond to the outer edge of a halo when galaxies are split on color. This occurs most distinctly for blue galaxies that are rapidly quenched and become red before reaching apocenter.

We demonstrate that even in the absence of any velocity information, the profiles can be used to extract time evolution information for segments of the galaxy population. The time elapsed between crossing into the cluster and reaching $r_{smin}$ traces the maximum time members of a galaxy population have been inside a cluster, in a model-independent way. Populations of galaxies with different SFRs can be mapped to halos that have entered clusters at different times. This allows us to estimate the timescale(s) associated with galaxy quenching.

Using only photometric data and the spatial distribution of galaxies, we are thus able to map galaxies onto different regions of the 3D phase space of dark matter halos. Just like the splashback radius traces the boundary between the virialized, multistreaming region and the infall region, $r_{smin}$ traces the discontinuity in the phase-space distributions even for galaxies that do not reach splashback. While the galaxy distribution in the inner regions of massive clusters can be significantly sensitive to baryonic physics, the location of $r_{smin}$ is more robust. Further, the information in the full density profiles, excluding the innermost regions but including the slope minimum, provides more detailed constraints on the parameters for quenching.
We use CDM-only simulations to create mock galaxy catalogs that reproduce the observed galaxy distributions of different colors and study their density distribution and slope profiles. We use a simple quenching model to estimate the quenching timescales for galaxies in cluster-mass halos corresponding to our SZ sample. Our principal empirical findings are:

1. The location of the splashback radius for the SZ-selected clusters in ACT DR5, measured using the complete galaxy sample, is $2.4_{-0.4}^{+0.3}$ Mpc h$^{-1}$. This is consistent with theoretical expectations from N-body simulations.

2. The shape of the density profile of galaxies differ significantly as a function of galaxy color. The observed density profiles show evidence of shifts of $r_{smin}$, the location of the steepest slope; in particular, the blue galaxies appear to show a weak minimum at a smaller radius compared to red and green galaxies (see Figure 5).

3. By comparing the $r_{smin}$ estimated for galaxies with simulated subhalos accreted at different times, we find that blue galaxies largely live in the infall stream, have not reached splashback, and are likely to have been accreted between 1.1 and 1.5 Gyr ago. Red galaxies show a sharp slope minimum at $2.2 \text{ Mpc h}^{-1}$ and have been in the halo for more than 3.2 Gyr on average, while green galaxies have been in the halo for at least 2.3 Gyr on average and have nearly reached their apocenter (see Figure 2).

4. We use the entire color fraction profile of galaxies between $0.5 \text{ Mpc h}^{-1}$ to $2.6 \text{ Mpc h}^{-1}$ to obtain constraints on our quenching parameters (see Figure 8). We find that the delay time is $t_d = 1.0_{-0.3}^{+0.1}$ Gyr and the exponential quenching timescale is $t_q = 0.6_{-0.1}^{+0.1}$ Gyr. The constraints obtained using the profiles agree with the timescales obtained from the location of the slope minimum. The color fraction profiles thus imply a short quenching time and a longer delay time that is comparable to the pericenter-crossing time of the cluster.

5. Using color ratio profiles in the region outside the virial radius of the cluster, we find that the total time required for quenching of satellite galaxies in this region is longer than inside the cluster, as expected because these galaxies are likely to occupy group-sized halos. The delay time is constrained to be in the range $0.8 < t_d < 2.7$ Gyr, and the exponential quenching time is constrained to $t_q < 0.7$ Gyr (see Figure 8).

We note that the error bars quoted in this paper are statistical and derived from a comparison with $\Lambda$CDM simulations without baryons. The simulations can introduce systematic biases particularly in the innermost radial bins, as discussed in Section 4.3.4. Because we focus mainly on the region near the outskirts of galaxy clusters, these effects are not expected to be severe. We do not include any baryonic effects on subhalo orbits and evolution. These can impact the galaxy-halo connection, and a detailed comparison of these results with the density and slope profiles of galaxies in hydrodynamic simulations is a promising direction to explore in the future.

The splashback radius is a robust length scale within dark matter halos that can be used to understand the evolution of nonlinear, virialized structures in the universe and the evolution of galaxies within them. We emphasize that while the location of the minimum of the density slope appears to be a function of galaxy

color, the splashback radius, i.e., the location of the first apocenter of galaxy orbits, is not. The movement of the location of the minimum is due to the movement of the phase-space boundary between a single-stream and multistream region.

We note that the location of the steepest slope for red galaxies is consistent with the location of $r_{smin}$ for the entire galaxy sample, implying that the total galaxy sample provides a robust estimate of the location of splashback. This location, we expect, should also be consistent with the splashback radius measured from dark matter particles itself for the samples in our data—specifically galaxies that are not experiencing significant dynamical friction (Adhikari et al. 2016; More et al. 2016; Baxter et al. 2017; Chang et al. 2018). We also note that while the location of $r_{smin}$ for red galaxies is the same as that of the total galaxy sample, the feature itself is much sharper and deeper in the latter case. This is further evidence that the slope minimum (splashback radius in this case) traces the cluster boundary, where a caustic-like feature is expected from galaxies/particles on turnaround. The sharpening of the feature is similar to the sharpening observed in Adhikari et al. (2014), Figure 1, wherein splashback was measured by subselecting particles with low radial velocities, tracing the caustic at turnaround or apocenter more clearly.

In this paper we have mostly focused on the most massive clusters and on the mean trend for clusters in a range of redshifts. We find that blue galaxies trace recent accretion. In particular, they infer that blue galaxies should be present mostly in the infall stream. In theory, halo accretion history can differ as a function of redshift, for example a similar-mass halo at two different redshifts will have different growth rates; studying star formation evolution as a function of redshift using the distribution of galaxies in phase space is an intriguing possibility in the future. High-redshift spectroscopic studies of X-ray-detected clusters like Willis et al. (2020) and Miller et al. (2018) also provide the opportunity to study clusters in their early formation stages.

The analysis in this paper can naturally also be extended to lower-mass clusters and groups that can be found using optical selection or using low-mass X-ray and SZ clusters from upcoming surveys like eRosita (Predehl et al. 2010; Pillepich et al. 2012), Simons Observatory (Ade et al. 2019), and CMB-S4 (Abazajian et al. 2019). Further, we can study quenching timescales as a function of cluster properties; existing DES data can already be used to study quenching as a function of properties of the host cluster, e.g., stellar mass (Palmese et al. 2020), redshift, and other tracers of history. We defer this natural extension to future work. In the immediate future, DES Year 6 data will allow us to probe fainter galaxies with lower stellar mass than our current sample around clusters that are within the DES footprint. This will enable a study of quenching jointly as a function of galaxy stellar mass and SFR. Ongoing and future galaxy surveys should provide an unprecedented wealth of data that can be used in novel ways, like studying the galaxy distribution in light of the splashback radius, to understand more deeply the connection between galaxies and their dark matter halos.

We thank Arka Banerjee, Peter Behroozi, Benedikt Diemer, and Andrey Kravtsov for illuminating discussions and reviewing early drafts of the paper. We thank the anonymous referee for helpful suggestions. This work was supported in part by the US Department of Energy contract to SLAC no. DE-AC02-
76SF00515 and made use of computational resources at SLAC National Accelerator Laboratory, a US Department of Energy Office, and the Sherlock cluster at the Stanford Research Computing Center (SRCC); the authors are thankful for the support of the SLAC and SRCC computing teams. The CosmoSim database used in this paper is a service by the Leibniz-Institute for Astrophysics Potsdam (AIP). The Multi-Dark database was developed in cooperation with the Spanish MultiDark Consolider Project CSD2009-00064. J.P.H. acknowledges funding for SZ cluster studies from NSF grant number AST-1615657. N.S. acknowledges support from NSF grant numbers AST-1513618 and AST-1907657. This paper has gone through internal review by the DES collaboration and the ACT collaboration.

Funding for the DES Projects has been provided by the DOE and NSF(USA), MEC/MICINN/MINECO (Spain), STFC (UK), HEFCE (UK). NCSA (UIUC), KICP (U. Chicago), CCAPP (Ohio State), MIFPA (Texas A&M), CNPQ, FAPERJ, FINEP (Brazil), DFG (Germany), and the Collaborating Institutions in the Dark Energy Survey.

The Collaborating Institutions are Argonne Lab, UC Santa Cruz, University of Cambridge, CIEMAT-Madrid, University of Chicago, University College London, DES-Brazil Consortium, University of Edinburgh, ETH Zürich, Fermilab, University of Illinois, ICE (IEEC-CSIC), IFAE Barcelona, Lawrence Berkeley Lab, LMU München and the associated Excellence Cluster Universe, University of Michigan, NFS's NOIRLab, University of Nottingham, Ohio State University, University of Pennsylvania, University of Portsmouth, SLAC National Lab, Stanford University, University of Sussex, Texas A&M University, and the OzDES Membership Consortium.

Based in part on observations at Cerro Tololo Inter-American Observatory at NSF's NOIRLab (NOIRLab Prop. ID 2012B-0001; PI: J. Frieman), which is managed 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 NSF under grant Nos. AST-1138766 and AST-1536171. The DES participants from Spanish institutions are partially supported by MICINN under grants ESP2017-89838, PGC2018-094773, PGC2018-102021, SEV-2016-0588, SEV-2016-0597, and MDM-2015-0509, some of which include ERDF funds from the European Union. 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 Brazilian Instituto Nacional de Ciência e Tecnologia (INCT) do e-Universo (CNPq grant 465376/2014-2).

This manuscript has been authored by Fermi Research Alliance, LLC under Contract No. DE-AC02-07CH11359 with the US Department of Energy, Office of Science, and Office of High Energy Physics.

The ACT project is supported by the US National Science Foundation through awards AST-1440226, AST-0965625, and AST-0408698, as well as awards PHY-1214379 and PHY-0855887. Funding was also provided by Princeton University, the University of Pennsylvania, and a Canada Foundation for Innovation (CFI) award to UBC. ACT operates in the Parque Astronómico Atacama in northern Chile under the auspices of

the Comisión Nacional de Investigación Científica y Tecnológica de Chile (CONICYT). Computations were performed on the GPC supercomputer at the SciNet HPC Consortium and on the hippo cluster at the University of KwaZulu-Natal. SciNet is funded by the CFI under the auspices of Compute Canada, the Government of Ontario, the Ontario Research Fund—Research Excellence, and the University of Toronto. The development of multichroic detectors and lenses was supported by NASA grants NNX13AE56G and NNX14AB58G.

# Appendix

## Splashback as a Function of the Orbital Angular Momentum of Particle Orbits

One possible explanation for the movement of the measured minimum of the density profile may be related to the distribution of orbits of galaxies. It is possible that galaxies that are on tangential orbits quench more slowly than galaxies that are on radial orbits. If tangential and radial orbits reach apocenters at different locations on average, we may expect that the movement of splashback between red, green, and blue galaxies may be explained by the difference in orbital histories. To check this hypothesis, we use simulations to study the splashback radius of particles with different pericentersg.

We track the particles orbits from the time they cross into the $4\ \text{Mpc}\ h^{-1}$ comoving distance from the cluster center and find the time and location of their first pericentric passage. Orbits with low pericenters correspond to particles that come in on radial orbits, while those with high pericenters come in with high angular momentum in tangential orbits. We randomly assign colors to all particles within $4\ \text{Mpc}\ h^{-1}$. We assign particles that have pericenter below $0.6\ \text{Mpc}\ h^{-1}$ random colors with a red fraction of 0.8 and particles that have pericenters larger than $0.6\ \text{Mpc}\ h^{-1}$ a red fraction of 0.4. In Figure 12 we show the logarithmic slope of the density profiles of particles

<figure><img src="image_13.png" /><figcaption>**Figure 12.** The slope of the density as a function of radius for particles that have been assigned colors based on their different orbital pericenters. The blue curve corresponds to particles that have on average high pericenters and red curve to particles that have low pericenters inside the halo. The location of the minimum of the slope remains unchanged in this toy model.</figcaption></figure>
that have been assigned red colors and those that have been assigned blue colors. We find that the location of the outer minimum does not change for the two samples. We find the same result if we alter the pericenter limit of 0.6 Mpc h⁻¹ to a different value.

Therefore, we conclude that the orbital differences in particles do not affect the splashback radius to a large radius; however, the inner density and slopes of particles can be significantly affected by it.

## ORCID iDs

Susmita Adhikari https://orcid.org/0000-0002-0298-4432

Matt Hilton https://orcid.org/0000-0002-8490-8117

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

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

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

J. Richard Bond https://orcid.org/0000-0003-2358-9949

Sebastian Bocquet https://orcid.org/0000-0002-4900-805X

Steve K. Choi https://orcid.org/0000-0002-9113-7058

Joseph DeRose https://orcid.org/0000-0002-0728-0960

Mark Devlin https://orcid.org/0000-0002-3169-9761

Simone Ferraro https://orcid.org/0000-0003-4992-7854

John P. Hughes https://orcid.org/0000-0002-8816-6800

Mathew S. Madhavacheril https://orcid.org/0000-0001-6740-5350

Lyman A. Page https://orcid.org/0000-0002-9828-3525

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

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

Neelima Sehgal https://orcid.org/0000-0002-9674-4527

Cristóbal Sifón https://orcid.org/0000-0002-8149-1352

Ed Wollack https://orcid.org/0000-0002-7567-4451

Zhilei Xu https://orcid.org/0000-0001-5112-2567

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

Santiago Avila https://orcid.org/0000-0001-5043-3662

Emmanuel Bertin https://orcid.org/0000-0002-3602-3664

David Brooks https://orcid.org/0000-0002-8458-5047

David L. Burke https://orcid.org/0000-0003-1866-1950

Jorge Carretero https://orcid.org/0000-0002-3130-0204

Francisco Javier Castander https://orcid.org/0000-0001-7316-4573

Matteo Costanzi https://orcid.org/0000-0001-8158-1449

Shantanu Desai https://orcid.org/0000-0002-0466-3288

Pablo Fosalba https://orcid.org/0000-0002-1510-5214

Juan García-Bellido https://orcid.org/0000-0002-9370-8360

Enrique Gaztanaga https://orcid.org/0000-0001-9632-0815

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

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

Julia Gschwend https://orcid.org/0000-0003-3023-8362

Samuel R. Hinton https://orcid.org/0000-0003-2071-9349

Devon L. Hollowood https://orcid.org/0000-0002-9369-4157

Kyler Kuehn https://orcid.org/0000-0003-0120-0808

Nikolay Kuropatkin https://orcid.org/0000-0003-2511-0946

Jennifer L. Marshall https://orcid.org/0000-0003-0710-9474

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

Felipe Menanteau https://orcid.org/0000-0002-1372-2534

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

Eusebio Sanchez https://orcid.org/0000-0002-9646-8198

Ignacio Sevilla-Noarbe https://orcid.org/0000-0002-1831-1953

Mathew Smith https://orcid.org/0000-0002-3321-1432

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

## References

Abadi, M. G., Moore, B., & Bower, R. G. 1999, MNRAS, 308, 947

Abazajian, K. N., Adelman-McCarthy, J. K., Agüeros, M. A., et al. 2009, ApJS, 182, 543

Abazajian, K., Addison, G., Adshead, P., et al. 2019, arXiv:1907.04473

Abbott, T., et al. 2005, arXiv:astro-ph/0510346

Abell, P. A., Allison, J., Anderson, S. F., et al. 2009, arXiv:0912.0201

Ade, P. A. R., Aghanim, N., Armitage-Caplan, C., et al. 2014, A&A, 571, A16

Ade, P., Aguirre, J., Ahmed, Z., et al. 2019, JCAP, 02, 056

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

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

Adhikari, S., Dalal, N., More, S., & Wetzel, A. 2019, ApJ, 878, 9

Adhikari, S., Sakstein, J., Jain, B., Dalal, N., & Li, B. 2018, JCAP, 1811, 033

Aiola, S., Palabrese, E., Maurin, L., et al. 2020, JCAP, 12, 047

Arnaud, M., Pratt, G. W., Piffaretti, R., et al. 2010, A&A, 517, A92

Aung, H., Nagai, D., Rozo, E., & Garcia, R. 2021, MNRAS, 502, 1041

Bahe, Y. M., & McCarthy, I. G. 2015, MNRAS, 447, 969

Balogh, M., Morris, S. L., Yee, H. K. C., Carlberg, R. G., & Ellingson, E., 1997, ApJ, 488, L75

Banerjee, A., Adhikari, S., Dalal, N., More, S., & Kravtsov, A. 2020, JCAP, 02, 024

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

Behroozi, P., Wechsler, R. H., Hearin, A. P., & Conroy, C. 2019, MNRAS, 488, 3143

Behroozi, P. S., Wechsler, R. H., Wu, H.-Y., et al. 2012, ApJ, 763, 18

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

Bell, E. F., Wolf, C., Meisenheimer, K., et al. 2004, ApJ, 608, 752

Bullock, J. S., Wechsler, R. H., & Somerville, R. S. 2002, MNRAS, 329, 246

Chang, C., Baxter, E., Jain, B., et al. 2018, ApJ, 864, 83

Conroy, C., & Wechsler, R. H. 2009, ApJ, 696, 620

Conroy, C., Wechsler, R. H., & Kravtsov, A. V. 2006, ApJ, 647, 201

Contigiani, O., Hoekstra, H., & Bahé, Y. 2019, MNRAS, 485, 408

Cooray, A. 2006, MNRAS, 365, 842

De Vicente, J., Sánchez, E., & Sevilla-Noarbe, I. 2016, MNRAS, 459, 3078

Diemer, B. 2017, ApJS, 231, 5

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

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

Dressler, A. 1980, ApJ, 236, 351

Dressler, A., & Gunn, J. E. 1983, ApJ, 270, 7

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

Einasto, J. 1965, TrAlm, 5, 87, 1965

Flaugher, B., Diehl, H. T., Honscheid, K., et al. 2015, AJ, 150, 150

Foreman-Mackey, D., Conley, A., Meierjurgen Farr, W., et al. 2013, emce: The MCMC Hammer, Astrophysics Source Code Library, ascl:1303.002

Förster Schreiber, N. M., Genzel, R., Bouché, N., et al. 2009, ApJ, 706, 1364

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

Hilton, M., Hasselfield, M., Sifón, C., et al. 2018, ApJS, 235, 20

Hilton, M., Sifón, C., Naess, S., et al. 2021, ApJS, 253, 3

Kravtsov, A., & Borgani, S. 2012, ARA&A, 50, 353

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

Larson, R. B., Tinsley, B. M., & Caldwell, C. N. 1980, ApJ, 237, 692

Laureijs, R., Amiaux, J., Arduini, S., et al. 2011, arXiv:1110.3193

Lee, S.-K., Ferguson, H. C., Somerville, R. S., Wiklind, T., & Giavalisco, M. 2010, ApJ, 725, 1644

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

Melin, J. B., Bartlett, J. G., & Delabrouille, J. 2006, A&A, 459, 341

Miller, T., Chapman, S. C., Aravena, M., et al. 2018, Natur, 556, 469

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

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

Morganson, E., Gruendl, R. A., Menanteau, F., et al. 2018, PASP, 130, 074501

Moustakas, L. A., & Somerville, R. S. 2002, ApJ, 577, 1

Murata, R., Sunayama, T., Oguri, M., et al. 2020, PASJ, 72, 64

Naess, S., Aiola, S., Austermann, J. E., et al. 2020, JCAP, 12, 046

Nagai, D., Kravtsov, A. V., & Vikhlinin, A. 2007, ApJ, 668, 1

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

Noeske, K. G., Faber, S. M., Weiner, B. J., et al. 2007, ApJL, 66, L47

Norberg, P., Baugh, C. M., Gaztañaga, E., & Croton, D. J. 2009, MNRAS, 396, 19

Oemler, A., Jr. 1974, ApJ, 194, 1

Okumura, T., Nishimichi, T., Umetsu, K., & Osato, K. 2017, arXiv:1706.08860
Okumura, T., Nishimichi, T., Umetsu, K., & Osato, K. 2018, PhRvD, 98, 023523

Oman, K. A., & Hudson, M. J. 2016, MNRAS, 463, 3083

Orsi, A. A., & Angulo, R. E. 2018, MNRAS, 475, 2530

Palmese, A., Annis, J., Burgad, J., et al. 2020, MNRAS, 493, 4591

Papovich, C., Dickinson, M., & Ferguson, H. C. 2001, ApJ, 559, 620

Pillepich, A., Porciani, C., & Reiprich, T. H. 2012, MNRAS, 422, 44

Poggianti, B., Smail, I., Dressler, A., et al. 1999, ApJ, 518, 576

Predehl, P., Andritschke, R., Böhringer, H., et al. 2010, Proc. SPIE, 7732, 77320U

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

Rhee, J., Smith, R., Choi, H., et al. 2020, ApJS, 247, 45

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

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

Saro, A., Bocquet, S., Rozo, E., et al. 2015, MNRAS, 454, 2305

Sevilla-Noarbe, I., Bechtol, K., Carrasco Kind, M., et al. 2021, ApJS, 254, 24

Shapley, A. E., Steidel, C. C., Erb, D. K., et al. 2005, ApJ, 626, 698

Shi, X. 2016, MNRAS, 459, 3711

Shin, T., Adhikari, S., Baxter, E. J., et al. 2019, MNRAS, 487, 2900

Somerville, R. S., & Davé, R. 2015, Ann. Rev. Astron. Astrophys, 53, 51

Springel, V. 2005, MNRAS, 364, 1105

Strateva, I., Ivezić, Ž., Knapp, G. R., et al. 2001, AJ, 122, 1861

Sunyaev, R., & Zel'dovich, Y. 1972, Astrophys. Space Phys., 4, 173

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

Tinker, J. L., Leauthaud, A., Bundy, K., et al. 2013, ApJ, 778, 93

Tinker, J. L., Wechsler, R. H., & Zheng, Z. 2010, ApJ, 709, 67

Tinker, J. L., & Wetzel, A. R. 2010, ApJ, 719, 88

Tomooka, P., Rozo, E., Wagoner, E. L., et al. 2020, MNRAS, 499, 1291

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

van den Bosch, F. C., Ogiya, G., Hahn, O., & Burkert, A. 2018, Mon. Not. Roy. Astron. Soc, 474, 3043

Wechsler, R. H., Bullock, J. S., Primack, J. R., Kravtsov, A. V., & Dekel, A. 2002, ApJ, 568, 52

Wechsler, R. H., & Tinker, J. L. 2018, ARA&A, 56, 435

Wetzel, A. R., Tinker, J. L., Conroy, C., & Bosch, F. C. d. 2013, MNRAS, 432, 336

Williamson, R., Benson, B. A., High, F. W., et al. 2011, ApJ, 738, 139

Willis, J., Canning, R., Noordeh, E., et al. 2020, Natur, 577, 39

Xhakaj, E., Diemer, B., Leauthaud, A., et al. 2020, MNRAS, 499, 3534

Zhang, Y., McKay, T. A., Bertin, E., et al. 2015, PASP, 127, 1183

Zürcher, D., & More, S. 2019, ApJ, 874, 184