# The Observed Evolution of the Stellar Mass-Halo Mass Relation for Brightest Central Galaxies

Jesse B. Golden-Marx, C. J. Miller, Y. Zhang, R. L. C. Ogando, A. Palmese, T. M. C. Abbott, M. Aguena, S. Allam, F. Andrade-Oliveira, J. Annis, et al.

► To cite this version:

Jesse B. Golden-Marx, C. J. Miller, Y. Zhang, R. L. C. Ogando, A. Palmese, et al. The Observed Evolution of the Stellar Mass-Halo Mass Relation for Brightest Central Galaxies. The Astrophysical Journal, 2022, 928, 10.3847/1538-4357/ac4cb4. insu-03748291

HAL Id: insu-03748291

https://insu.hal.science/insu-03748291

Submitted on 9 Aug 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
OPEN ACCESS

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

CrossMark

# The Observed Evolution of the Stellar Mass–Halo Mass Relation for Brightest Central Galaxies

Jesse B. Golden-Marx<sup>1,2</sup> ◊, C. J. Miller<sup>2,3</sup>, Y. Zhang<sup>4</sup> ◊, R. L. C. Ogando<sup>5,6</sup> ◊, A. Palmese<sup>4,7</sup> ◊, T. M. C. Abbott<sup>8</sup> ◊, M. Aguena<sup>5</sup> ◊, S. Allam<sup>4</sup>, F. Andrade-Oliveira<sup>5,9</sup>, J. Annis<sup>4</sup> ◊, D. Bacon<sup>10</sup> ◊, E. Bertin<sup>11,12</sup> ◊, D. Brooks<sup>13</sup> ◊, E. Buckley-Geer<sup>4,14</sup>, A. Carnero Rosell<sup>5,15,16</sup> ◊, M. Carrasco Kind<sup>17,18</sup> ◊, F. J. Castander<sup>19,20</sup> ◊, M. Costanzi<sup>21,22,23</sup> ◊, M. Crocce<sup>19,20</sup> ◊, L. N. da Costa<sup>5,6</sup> ◊, M. E. S. Pereira<sup>3</sup>, J. De Vicente<sup>24</sup> ◊, S. Desai<sup>25</sup> ◊, H. T. Diehl<sup>4</sup> ◊, P. Doel<sup>13</sup>, A. Drlica-Wagner<sup>4,7,14</sup> ◊, S. Everett<sup>26</sup>, A. E. Evrard<sup>2,3</sup> ◊, I. Ferrero<sup>27</sup> ◊, B. Flaugher<sup>4</sup>, P. Fosalba<sup>19,20</sup> ◊, J. Frieman<sup>4,7</sup>, J. García-Bellido<sup>28</sup> ◊, E. Gaztanaga<sup>19,20</sup> ◊, D. W. Gerdes<sup>2,3</sup> ◊, D. Gruen<sup>29,30,31</sup> ◊, R. A. Gruendl<sup>17,18</sup> ◊, J. Gschwend<sup>5,6</sup> ◊, G. Gutierrez<sup>4</sup>, W. G. Hartley<sup>32</sup>, S. R. Hinton<sup>33</sup> ◊, D. L. Hollowood<sup>26</sup> ◊, K. Honscheid<sup>34,35</sup>, B. Hoyle<sup>36</sup>, D. J. James<sup>37</sup>, T. Jeltema<sup>26</sup> ◊, A. G. Kim<sup>38</sup>, E. Krause<sup>39</sup>, K. Kuehn<sup>40,41</sup> ◊, N. Kuropatkin<sup>4</sup> ◊, O. Lahav<sup>13</sup>, M. Lima<sup>5,42</sup>, M. A. G. Maia<sup>5,6</sup>, J. L. Marshall<sup>43</sup> ◊, P. Melchior<sup>44</sup> ◊, F. Menanteau<sup>17,18</sup> ◊, R. Miquel<sup>45,46</sup> ◊, J. J. Mohr<sup>36,47</sup>, R. Morgan<sup>48</sup> ◊, F. Paz-Chinchón<sup>17,49</sup> ◊, D. Petravick<sup>17</sup>, A. Pieres<sup>5,6</sup> ◊, A. A. Plazas Malagón<sup>44</sup>, J. Prat<sup>7,14</sup>, A. K. Romer<sup>50</sup> ◊, E. Sanchez<sup>24</sup> ◊, B. Santiago<sup>5,51</sup>, V. Scarpine<sup>4</sup>, M. Schubnell<sup>3</sup>, S. Serrano<sup>19,20</sup> ◊, I. Sevilla-Noarbe<sup>24</sup> ◊, M. Smith<sup>52</sup>, M. Soares-Santos<sup>3</sup> ◊, E. Suchyta<sup>53</sup>, G. Tarle<sup>3</sup> ◊, and T. N. Varga<sup>47,54</sup>

(DES Collaboration)

<sup>1</sup> Department of Astronomy, Shanghai Jiao Tong University, Shanghai 200240, People's Republic of China; jessegm@sjtu.edu.cn

<sup>2</sup> Department of Astronomy, University of Michigan, Ann Arbor, MI 48109, USA

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

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

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

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

<sup>7</sup> Kavli Institute for Cosmological Physics, University of Chicago, Chicago, IL 60637, USA

<sup>8</sup> Cerro Tololo Inter-American Observatory, NSF's National Optical-Infrared Astronomy Research Laboratory, Casilla 603, La Serena, Chile

<sup>9</sup> Instituto de Física Teórica, Universidade Estadual Paulista, São Paulo, Brazil

<sup>10</sup> Institute of Cosmology and Gravitation, University of Portsmouth, Portsmouth, PO1 3FX, UK

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

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

<sup>13</sup> Department of Physics & Astronomy, University College London, Gower Street, London, WC1E 6BT, UK

<sup>14</sup> Department of Astronomy and Astrophysics, University of Chicago, Chicago, IL 60637, USA

<sup>15</sup> Instituto de Astrofísica de Canarias, E-38205 La Laguna, Tenerife, Spain

<sup>16</sup> Universidad de La Laguna, Dpto. Astrofísica, E-38206 La Laguna, Tenerife, Spain

<sup>17</sup> Center for Astrophysical Surveys, National Center for Supercomputing Applications, 1205 West Clark St., Urbana, IL 61801, USA

<sup>18</sup> Department of Astronomy, University of Illinois at Urbana-Champaign, 1002 W. Green Street, Urbana, IL 61801, USA

<sup>19</sup> Institut d'Estudis Espacials de Catalunya (IEEC), E-08034 Barcelona, Spain

<sup>20</sup> Institute of Space Sciences (ICE, CSIC), Campus UAB, Carrer de Can Magrans, s/n, E-08193 Barcelona, Spain

<sup>21</sup> Astronomy Unit, Department of Physics, University of Trieste, via Tiepolo 11, I-34131 Trieste, Italy

<sup>22</sup> INAF-Osservatorio Astronomico di Trieste, via G. B. Tiepolo 11, I-34143 Trieste, Italy

<sup>23</sup> Institute for Fundamental Physics of the Universe, Via Beirut 2, I-34014 Trieste, Italy

<sup>24</sup> Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas (CIEMAT), Madrid, Spain

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

<sup>26</sup> Santa Cruz Institute for Particle Physics, Santa Cruz, CA 95064, USA

<sup>27</sup> Institute of Theoretical Astrophysics, University of Oslo. P.O. Box 1029 Blindern, NO-0315 Oslo, Norway

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

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

<sup>30</sup> Kavli Institute for Particle Astrophysics & Cosmology, Stanford University, P. O. Box 2450, Stanford, CA 94305, USA

<sup>31</sup> SLAC National Accelerator Laboratory, Menlo Park, CA 94025, USA

<sup>32</sup> Department of Astronomy, University of Geneva, ch. d'Écogia 16, CH-1290 Versoix, Switzerland

<sup>33</sup> School of Mathematics and Physics, University of Queensland, Brisbane, QLD 4072, Australia

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

<sup>35</sup> Department of Physics, The Ohio State University, Columbus, OH 43210, USA

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

<sup>37</sup> Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA

<sup>38</sup> Lawrence Berkeley National Laboratory, 1 Cyclotron Road, Berkeley, CA 94720, USA

<sup>39</sup> Department of Astronomy/Steward Observatory, University of Arizona, 933 North Cherry Avenue, Tucson, AZ 85721-0065, USA

<sup>40</sup> Australian Astronomical Optics, Macquarie University, North Ryde, NSW 2113, Australia

<sup>41</sup> Lowell Observatory, 1400 Mars Hill Rd, Flagstaff, AZ 86001, USA

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

<sup>43</sup> George P. and Cynthia Woods Mitchell Institute for Fundamental Physics and Astronomy, and Department of Physics and Astronomy, Texas A&M University, College Station, TX 77843, USA

<sup>44</sup> Department of Astrophysical Sciences, Princeton University, Peyton Hall, Princeton, NJ 08544, USA

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

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

<sup>47</sup> Max Planck Institute for Extraterrestrial Physics, Giessenbachstrasse, D-85748 Garching, Germany

<sup>48</sup> Physics Department, 2320 Chamberlin Hall, University of Wisconsin-Madison, 1150 University Avenue Madison, WI 53706-1390, USA

<sup>49</sup> Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK

<sup>50</sup> Department of Physics and Astronomy, Pevensey Building, University of Sussex, Brighton, BN1 9QH, UK

<sup>51</sup> Instituto de Física, UFRGS, Caixa Postal 15051, Porto Alegre, RS—91501-970, Brazil

<sup>52</sup> School of Physics and Astronomy, University of Southampton, Southampton, SO17 1BJ, UK
**Abstract**

We quantify evolution in the cluster-scale stellar mass–halo mass (SMHM) relation’s parameters using 2323 clusters and brightest central galaxies (BCGs) over the redshift range 0.03 ≲ z ≲ 0.60. The precision on the inferred SMHM parameters is improved by including the magnitude gap ($m_{gap}$) between the BCG and fourth-brightest cluster member (M14) as a third parameter in the SMHM relation. At fixed halo mass, accounting for $m_{gap}$, through a stretch parameter, reduces the SMHM relation’s intrinsic scatter. To explore this redshift range, we use clusters, BCGs, and cluster members identified using the Sloan Digital Sky Survey C4 and redMaPPer cluster catalogs and the Dark Energy Survey redMaPPer catalog. Through this joint analysis, we detect no systematic differences in BCG stellar mass, $m_{gap}$, and cluster mass (inferred from richness) between the data sets. We utilize the Pareto function to quantify each parameter’s evolution. We confirm prior findings of negative evolution in the SMHM relation’s slope (3.5$\sigma$), and detect negative evolution in the stretch parameter (4.0$\sigma$) and positive evolution in the offset parameter (5.8$\sigma$). This observed evolution, combined with the absence of BCG growth, when stellar mass is measured within 50 kpc, suggests that this evolution results from changes in the cluster’s $m_{gap}$. For this to occur, late-term growth must be in the intracluster light surrounding the BCG. We also compare the observed results to IllustrisTNG 300-1 cosmological hydrodynamic simulations and find modest qualitative agreement. However, the simulations lack the evolutionary features detected in the real data.

*Unified Astronomy Thesaurus concepts: Galaxy clusters (584); Galaxy evolution (594); Brightest cluster galaxies (181)*

**1. Introduction**

The stellar mass–halo mass (SMHM) relation is the primary mechanism used to quantify and characterize the galaxy–dark matter halo connection. Since multiple versions of the SMHM relation exist, we note that, for this analysis, we study the brightest central galaxy (BCG) SMHM relation for galaxy clusters ($\log_{10}(M_{halo}/(M_{\odot}/h)) \ge 14.0$), the linear correlation that compares the stellar mass of the BCG to the total halo cluster mass, which includes the dark matter. We do not account for the stellar mass contained within the satellites in this analysis. The parameters measured as part of the SMHM relation can constrain galaxy formation models, including the amount of active galactic nucleus feedback in central galaxies (Kravtsov et al. 2018). The intrinsic scatter in stellar mass at fixed halo mass ($\sigma_{int}$) can constrain the processes responsible for quenching star formation in central galaxies (Tinker 2017), as well as characterize dark matter halo assembly (Gu et al. 2016). Additionally, the redshift evolutions of the slope and scatter provide insights into how BCGs grow and evolve over cosmic time (Gu et al. 2016; Golden-Marx & Miller 2019).

BCGs, which solely make up the stellar mass portion of the cluster-scale SMHM relation, are massive, radially extended, elliptical galaxies that emit a significant fraction of the total light within their host cluster (Schombert 1986; Jones et al. 2000; Lin & Mohr 2004; Bernardi et al. 2007; Lauer et al. 2007; von der Linden et al. 2007; Aguerri et al. 2011; Brough et al. 2011; Proctor et al. 2011; Harrison et al. 2012). BCGs are located near the cluster’s X-ray center. This location, along with their hierarchical formation (e.g., De Lucia & Blaizot 2007; Oser et al. 2010; van Dokkum et al. 2010), lead to correlations between the BCG properties and those of their host cluster (Jones & Forman 1984; Rhee & Latour 1991; Lin & Mohr 2004; Lauer et al. 2014). Additionally, BCGs are

surrounded by diffuse halos of intracluster light (ICL; Zwicky 1933, 1951), which are observed to extend radially as far out as ≈1 Mpc from the center of the BCG (Zhang et al. 2019a), and mostly result from the BCG’s hierarchical assembly (Murante et al. 2007).

BCGs grow “inside-out” (van Dokkum et al. 2010), following a two-phase formation scenario (Oser et al. 2010); at high redshifts ($z > 2$), the dense core ($r < ≈10$ kpc) forms via in situ star formation, and at lower redshifts ($z < 2$), the outer envelope grows hierarchically via major/minor mergers. The two-phase formation scenario is supported both by observations (van Dokkum et al. 2010; Huang et al. 2018) and dark matter–only cosmological simulations that use empirical or semi-analytic models to quantify central galaxy stellar mass growth (e.g., Croton et al. 2006; De Lucia & Blaizot 2007; Guo et al. 2011; Tonini et al. 2012; Shankar et al. 2015).

As a result of this “inside-out” growth, all of the information about the BCG’s recent stellar mass growth is contained within the BCG’s outer envelope, which extends to the ICL (Oser et al. 2010; van Dokkum et al. 2010). Moreover, recent observations suggest that the majority of the BCG’s stellar mass may be contained within a radial aperture of 100 kpc, centered on the BCG (Huang et al. 2018), and that the 100 kpc boundary may represent a transitional regime between the BCG’s outer envelope and the ICL (Zhang et al. 2019a). Therefore, when characterizing the BCG evolution associated with the parameters of the SMHM relation, it is vital to measure the BCG photometry within large radii, as opposed to the more commonly used 20–30 kpc aperture radius (e.g., Lin et al. 2013; Zhang et al. 2016; Lin et al. 2017), as discussed in Golden-Marx & Miller (2019), hereafter referred to as GM&M19. More specifically, including stellar mass within large radii strengthens the correlation between the BCG’s stellar mass and halo mass (Moster et al. 2018; GM&M19). However, to yield a stronger correlation, one must also incorporate a third parameter related to BCG growth.

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

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

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

<footnote n="54">⁵⁴ Universitäts-Sternwarte, Fakultät für Physik, Ludwig-Maximilians Universität München, Scheinerstr. 1, D-81679 München, Germany
*Received 2021 July 1; revised 2021 December 20; accepted 2022 January 17; published 2022 March 24*</footnote>
One observational measurement inherently tied to BCG hierarchical growth is $m_{gap}$, the difference in r-band magnitude between the BCG and the fourth-brightest cluster member within half the radius enclosing 200 times the critical density of the universe ($R_{200}$) (Dariush et al. 2010). Throughout this paper, we refer to the $m_{gap}$ between the BCG and the fourth-brightest member as M14. Using N-body simulations, Solanes et al. (2016) find that BCG stellar mass linearly increases with the number of progenitor galaxies. Since the BCG's central location leads to faster merger growth than in noncentral galaxies, as BCGs grow hierarchically, their stellar mass and magnitude increase, while those same parameters for the fourth-brightest member remain the same (unless that galaxy is involved in the BCG merger). Therefore, BCG growth results in a corresponding increase in the magnitude gap ($m_{gap}$), and yields the correlation between $m_{gap}$ and BCG stellar mass (Harrison et al. 2012; Golden-Marx & Miller 2018). Thus, it follows that $m_{gap}$ can be thought of as a statistical latent parameter within the cluster SMHM relation, as first presented by Golden-Marx & Miller (2018), hereafter referred to as GM&M18. Additionally, the correlation between $m_{gap}$ and stellar mass, which results from hierarchical growth, suggests that $m_{gap}$ may be a tracer of formation redshift (GM&M18). Given that, we use M14 in this analysis, as opposed to alternative $m_{gap}$ measures, because Dariush et al. (2010) find that systems with large M14 measurements are more efficiently identified as earlier-forming systems than those with large values of M12, the $m_{gap}$ between the BCG and the second-brightest cluster member.

GM&M18 incorporate $m_{gap}$ into the cluster SMHM relation as a linear stretch parameter, which acts to spread the observed range of stellar masses at fixed halo mass. This is just clarifying one of the primary components of the intrinsic scatter in the classic SMHM relation (i.e., without using $m_{gap}$ as a third parameter). In other words, the intrinsic scatter measured using the standard two-parameter SMHM relation is larger than the intrinsic scatter in the SMHM relation after accounting for this third parameter (e.g., akin to a fundamental plane). The inferred intrinsic scatter in the SMHM relation found by GM&M18 is less than 0.1 dex, which is smaller than previous studies by as much as a factor of two (e.g., Zu & Mandelbaum 2015; Gu et al. 2016; Tinker et al. 2017; Pillepich et al. 2018a; Kravtsov et al. 2018). Since the scatter in the SMHM relation is quite small, the other parameters can be more precisely constrained, but only after incorporating the stretch parameter, which measures the strength of the correlation between $m_{gap}$ and stellar mass at fixed halo mass.

Next, consider the evolution of the SMHM relation over cosmic time. This evolution can inform us about how BCGs grow over time, as well as the fraction of stellar material ejected into the ICL as a result of major/minor mergers. Using empirical models with abundance-matching techniques to infer halo masses, Behroozi et al. (2013) and Moster et al. (2013) find that the slope of the SMHM relation increases by a factor of 1.5–2.0 from $z=1.0$ to $z=0.0$, which would suggest that BCGs continue to grow significantly via mergers over this redshift range. Moreover, Moster et al. (2013) detect moderate evolution from $z=0.5$ to $z=0$. In contrast, Pillepich et al. (2018a) and Engler et al. (2021), using the IllustrisTNG 300-1 cosmological hydrodynamic simulation, measure little change in the slope between $z=1.0$ and $z=0.0$. This absence of

evolution over this same redshift range was also found in Rodríguez-Puebla et al. (2017), who use semi-empirical modeling to study the evolution of the SMHM relation. In addition to the slope, the expected redshift evolution of the intrinsic scatter, $\sigma_{int}$, has also been investigated in models and simulations (Gu et al. 2016; Matthee et al. 2017; Pillepich et al. 2018a). However, as was the case for the slope, there is no consensus between these studies.

Since our analysis accounts for the satellite population via $m_{gap}$, it is also worth highlighting two recent results looking at the evolution within the total stellar mass of Dark Energy Survey (DES) clusters using the parameter $\mu_*$, the sum of the individual galaxy stellar masses weighted by their membership probability. Palmese et al. (2020) find no evolution in the correlation between $\mu_*$ and richness ($\lambda$). In contrast, Pereira et al. (2020) explicitly account for redshift evolution in their relation, and find weak redshift evolution. However, Pereira et al. (2020) note that their evolution is within the accepted uncertainty of the total stellar mass ($\approx 0.1$ dex). Therefore, as for the SMHM relation for BCGs, it is currently unclear how the total stellar mass of the cluster is evolving.

Using observational data, prior studies have been unable to constrain the SMHM relation's late-time redshift evolution (Oliva-Altamirano et al. 2014; Gozaliasl et al. 2016; Erfanianfar et al. 2019). However, by incorporating $m_{gap}$, GM&M19 placed the first statistically significant observational constraints on the redshift evolution of the slope of the SMHM relation. Over the redshift range $0.03 < z < 0.30$, the slope of the SMHM relation decreases by $\approx 0.20$ dex or 40%. To expand upon those results, the primary goal of this paper is to characterize the evolution of the cluster SMHM relation out to $z \sim 0.6$.

To constrain the evolution in the SMHM relation to higher redshifts, we combine the lower-redshift Sloan Digital Sky Survey (SDSS) data with data from the Dark Energy Survey Year 3 (DESY3) release (The Dark Energy Survey Collaboration 2005; Flaugher et al. 2015; Abbott et al. 2018). We choose DES data for a few key reasons. First, tens of thousands of galaxy clusters are identified in DES, and the survey is complete out to $z \approx 0.6$, significantly deeper than the redshift range probed by SDSS (e.g., Rykoff et al. 2014; Alam et al. 2015). While surveys such as the Hyper Suprime-Cam Subaru Strategic Program (Aihara et al. 2018) or the Atcama Cosmology Telescope (ACT; Hilton et al. 2021) may offer a similarly deep (or deeper) redshift coverage, those surveys do not provide a large enough sample of clusters to reduce the statistical uncertainty needed to provide tight constraints on the parameters associated with the SMHM relation using our Bayesian model. Additionally, the deep DES photometry allows us to accurately measure large-aperture photometry for our BCGs. DES provides a wide field of view around each BCG, allowing us to easily determine $m_{gap}$ as well. Finally, one further goal of this analysis is to create a homogeneous data set to study redshift evolution. Since the redMaPPer algorithm (Rykoff et al. 2014, 2016) has been applied to DES data, and there exists a set of clusters observed by both SDSS and DES, this makes the processes of creating such a homogeneous sample and determining the associated uncertainty simpler, since the membership restrictions and measurement methods applied to SDSS data can be similarly applied to DES data.
The outline for the remainder of this paper is as follows. In Section 2, we discuss the observational and simulated data (IllustrisTNG 300-1) used to measure the stellar masses, halo masses, and $m_{gap}$ values for our SMHM relation. In Section 3, we describe the hierarchical Bayesian Markov Chain Monte Carlo (MCMC) model used to evaluate the redshift evolution of the SMHM relation. In Section 4, we present our results. In Section 5, we discuss our findings and conclude.

Except for the case of the IllustrisTNG 300-1 simulated data, for which the cosmological parameters were previously defined ($\Omega_M = 0.3089$, $\Omega_\Lambda = 0.6911$, and $H_0 = 100 \text{ h km s}^{-1} \text{ Mpc}^{-1}$ with $h = 0.6774$), we assume a flat $\Lambda$CDM universe for our analysis, with $\Omega_M = 0.30$, $\Omega_\Lambda = 0.70$, and $H_0 = 100 \text{ h km s}^{-1} \text{ Mpc}^{-1}$ with $h = 0.7$.

## 2. Data

To characterize the SMHM relation and its evolution, we require mass measurements of the central galaxies, as well as enough satellite galaxies to infer $m_{gap}$ and the richness of the halo, the latter of which allows for an estimate of the halo's mass. To obtain these measurements over the desired redshift range, we utilize two survey data sets: SDSS Data Release 12 (SDSS DR12; Alam et al. 2015) and DESY3 (Flaugher et al. 2015; Abbott et al. 2018; Morganson et al. 2018), which are briefly summarized below.

DES is an optical to near-infrared photometric survey covering 5000 deg² in the South Galactic Cap in the DES *grizY* bands (for the purpose of this analysis, only the *g*-, *r*-, *i*-, and *z*-bands are used). In total, over 575 nights of observations were taken over a 6 yr period, beginning in 2013. The observations were taken at the Cerro Tololo International Observatory in Chile, using the ≈3 deg² CCD Dark Energy Camera (Flaugher et al. 2015) on the Blanco 4 m telescope. The data used in this analysis were taken over the first 3 yr of observations.

SDSS is a photometric survey with overlapping spectroscopic data collected by the Baryon Oscillation Spectroscopic Survey that covers a footprint in the northern sky of 14,055 deg² in the SDSS *griz* bands (for the purpose of this analysis, only the *g*-, *r*-, and *i*-bands are used), observed between 1998 and 2009. The observations were taken using the Sloan Foundation 2.5 m telescope at Apache Point Observatory in New Mexico. As in GM&M18 and GM&M19, the data used in this analysis come from SDSS DR12. The only difference between the data used in this analysis and in the prior studies is in the radii within which the SDSS BCG magnitudes are measured, as discussed in Section 2.2.

The galaxy clusters that are used in this analysis come from the low-redshift SDSS-C4 (Miller et al. 2005), SDSS-redMaPPer v6.3 (Rykoff et al. 2014), and DES-redMaPPer v6.4.22 ($\lambda > 20$ volume-limited (Rykoff et al. 2016) catalogs, where $\lambda$ is the redMaPPer richness measurement. These redMaPPer cluster catalogs have both high purity and high completeness over the redshift and $\lambda$ ranges that we are studying (Rykoff et al. 2014, 2016). However, no single cluster catalog individually covers the entire redshift range that we aim to study: $0.03 < z < 0.60$. The SDSS-C4 sample covers $0.03 < z < 0.15$, with $z_{med} = 0.08$; the SDSS-redMaPPer sample covers $0.08 < z < 0.30$; and the DES-redMaPPer sample covers $0.20 < z < 0.60$. Therefore, we combine the cluster catalogs to create one parent sample.

<figure><img src="image_4.png" /><figcaption>**Figure 1.** The redshift distribution of the combined sample of SDSS and DES data for this analysis. The dashed lines represent the edges of the eight bins used in this analysis.</figcaption></figure>

GM&M19 used clusters in the redshift range $0.08 \le z \le 0.12$ to characterize any differences in the halo mass, richness, central galaxy magnitude, stellar mass, and magnitude gaps between the SDSS-C4 and SDSS-redMaPPer data. By conducting a direct comparison of individual clusters in both data sets, they ruled out systematic differences in the mean observables (e.g., biases) between the two samples that could mimic real evolutionary trends in the SMHM relation. We conduct a similar analysis on an overlapping redshift region for the SDSS-redMaPPer and DES-redMaPPer clusters in this work, described in Section 2.6.

In our SMHM relation analysis, we constrain the evolution of the parameters with and without redshift binning to emphasize consistency in our statistical analysis. For the redshift-binned analysis, the parent sample of SDSS and DES clusters is divided into eight redshift bins, as shown in Figure 1 and given in the Appendix.

In the following subsections, we describe the measurements of our observables: BCG stellar mass, cluster mass (via richness), and $m_{gap}$. We will specifically highlight important differences compared to the measurement approaches used in GM&M19. In addition to the observables, our Bayesian analysis also requires priors on the measurement uncertainties. We estimate these uncertainties from a bootstrapped Bayesian analysis, described in Section 2.6.

### 2.1. BCG Identification

The BCG in every cluster is identified using a combination of visual identification, magnitude, color, and redshift. While the BCG is nominally the brightest galaxy at the center of the cluster's gravitational potential well, difficulties in BCG identification arise due to cluster centering accuracy, foreground/background contamination, photometric accuracy, etc. For the low-redshift SDSS-C4 clusters, the BCGs were visually confirmed. For most of the redMaPPer clusters, we use the statistically most probable BCG from the redMaPPer algorithm. However, in the overlap samples (i.e., the clusters that appear in two or more of the three catalogs), we visually confirm the BCG (a thorough discussion of miscentering in the redMaPPer algorithm is provided in, e.g., Hoshino et al. 2015 and Zhang et al. 2019b). We note that the BCG must lie within $0.5 \times R_{200,crit}$ of the cluster center (see Section 2.5 for details). The photometric algorithms and visual confirmations ensure
<figure><img src="image_5.png" /><figcaption>**Figure 2.** The left image is a coadded DES cluster, centered on the BCG. The right image is the SDSS image of the same BCG. Each image is a 1.5 by 1.5 box centered on the BCG. The insets in both panels are representations of the BCG after masking. The black circle is the 50 kpc radius, within which we measure the BCG stellar mass. Note that while the DES postage stamp appears more spherical than the SDSS postage stamp, we do not use shape information in this work.</figcaption></figure>

that the BCGs have similar colors to the rest of the galaxies within the so-called E/S0 ridgeline, and that the BCG morphologies exclude disk, disturbed, or merging galaxies.

## 2.2. BCG Light Profiles

GM&M19 found that the slope of the SMHM relation is dependent on the radius within which the BCG's stellar mass is measured. The SMHM relation's slope reaches an asymptote when the projected aperture used to estimate the stellar mass within a galaxy is between 60 and 100 kpc. Therefore, to homogeneously infer the slope, we use fixed physical aperture BCG magnitudes, as opposed to alternatives such as the Petrosian or Kron magnitudes (Petrosian 1976; Kron 1980).

For the SDSS BCGs, we use the SDSS pipeline-processed radial light profiles to measure the fixed aperture magnitudes. The DES pipeline does not provide radial aperture photometry. Therefore, we follow the procedure described in Zhang et al. (2019a) to measure the DES BCG light profiles.

We coadd and stack the DES individual image frames out to 0.15 from the BCG locations. In Figure 2, we compare a DES coadded and sky-subtracted image to an SDSS pipeline-processed image of the central region for a specific cluster. Unsurprisingly, the DES data reach a much lower surface brightness than the SDSS photometry. In the SDSS pipeline, the radial light profiles are constructed after nearby objects are masked. Masking removes the majority of the excess light associated with the neighboring galaxies, yielding a clean measurement of the radial light profile centered on the BCG, which includes the ICL. For the DES data, we mask all objects brighter than 30th mag in the i-band out to a radius of 2.5 $R_{\text{kron}}$ for each detected object. In the inset images of Figure 2, we compare the masked and sky-subtracted DES BCG to the SDSS Atlas image, the SDSS equivalent. In the DES masked, recovered image, the masked pixels are replaced with the radially averaged flux level.

After the mask is applied to the DES coadded image, we measure the BCG's radial light profile in annuli centered on the BCG. We subtract a background determined from the median flux at radii beyond 500 kpc from the BCG. In Figure 3, we compare the r-band light profiles for a single BCG in both SDSS and DES photometry, which shows that while the light measured within the central aperture is very similar, there is more light in the DES photometry when compared to the SDSS photometry at larger radii.

We identify 48 BCGs within $0.20 < z < 0.35$ observed in both DES and SDSS that have SDSS light profiles that are above the background to 100 kpc. We use this subset of data to quantitatively characterize the radially dependent magnitude differences between DES and SDSS BCG photometry. Beyond 50 kpc, the SDSS photometric measurements are consistently fainter than the DES photometry, as shown in Figure 4. The differences begin to grow beyond 50 kpc, with the DES magnitudes being nearly 0.5 mag brighter in the g-band compared to SDSS when measured at 100 kpc. The differences are less pronounced in the r- and i-bands, but are large enough to cause concern about using the 100 kpc aperture magnitude for BCGs. Based on this analysis, we choose 50 kpc as the BCG aperture magnitude for the remaining analyses. We use this aperture for all SDSS and DES BCGs. We note that since we do not use the SDSS z-band photometry, no comparison is made between the SDSS and DES z-bands.

As a final test, we compare the colors of the 48 BCGs that exist in both the SDSS and DES data. We convert the SDSS magnitudes to DES magnitudes using the available filter curves for each survey (e.g., Alam et al. 2015; Li et al. 2016; Burke et al. 2018). We then calculate the mean and standard deviation of the difference between the SDSS and DES BCG colors. We find that the mean (error) is 0.035 mag, with a standard deviation of 0.137 mag. Therefore, we find no bias between the two data sets for the BCG colors.
<figure><img src="image_6.png" /><figcaption>**Figure 3.** The SDSS and DES r-band light profiles for the BCG in Figure 2. For such a comparison, we scaled the SDSS photometry to match that of the DES photometry, since there are slight differences in the wave bands used for each survey. The DES profile matches the SDSS within the inner 50 kpc region, but then becomes brighter as more light above the background is observed in the deeper DES observations.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 4.** The difference between the SDSS and DES photometric measurements for 48 BCGs identified in both the SDSS- and DES-redMaPPer catalogs. We show the cumulative magnitude difference for the g-, r-, and i-bands as a function of the radial aperture. The SDSS and DES photometry begin to diverge at radii >50 kpc, particularly in the g-band, where differences between the two surveys are larger than the average magnitude measurement uncertainty on the BCG magnitude.</figcaption></figure>

## 2.3. BCG Stellar Masses

We use the observed radial light profiles to measure the projected BCG luminosity/magnitude within a 50 kpc aperture. The magnitudes in different bands allow us to estimate the stellar mass. We follow the same procedure as outlined in GM&M19, which is summarized here. For each cluster, redMaPPer assigns every galaxy a membership probability, $P_{\text{mem}}$, which is dependent on the cluster's richness, density profile, and background density (Rykoff et al. 2014, 2016). The $P_{\text{mem}} > 0.7$ members are then used to estimate the cluster's photometric redshift, which we use for the BCG.

Using the BCG apparent magnitudes, colors, and redshifts, our goal is to estimate a consistent stellar mass for a single and homogeneous population of galaxies. To do this using our three-color photometry, we use a passively evolving spectral model using the spectral energy distribution (SED) modeling software package EzGal (Mancone & Gonzalez 2012), because our available photometry limits us from statistically being able to constrain additional parameters such as burst times and

formation epochs. This model fitting allows us to infer a color-based stellar mass. We assume a Bruzual & Charlot (2003) stellar population synthesis model, a Salpeter (1955) initial mass function (IMF), a formation redshift of $z=4.9$, and a metallicity of 0.008, 0.02, or 0.05. We note that our stellar mass estimates are well within $1\sigma$ regardless of our choice of formation redshift. Furthermore, the choice of metallicity for each DES BCG is determined based on which model yielded the lowest chi-squared statistic between the measured and modeled photometry. We note that GM&M19 found that >99% of lower-redshift SDSS clusters are best constrained by the model when a metallicity of 0.008 is used; however, that fraction decreases to 87% for the high-redshift DES data. Based on the results from McDermid et al. (2015), we suspect that our choice of a subsolar metallicity likely results from our inclusion of light from the outskirts of the BCG envelope. McDermid et al. (2015) find that the metallicity of early-type galaxies decreases from the core region out to the outskirts measured at the effective radius. This selection is further supported by Emsellem et al. (2007), who show that massive early-type galaxies are predominantly slow rotators, and Krajnović et al. (1919), who find that such slow rotators are characterized by a subsolar metallicity. Moreover, though farther out into the ICL than we explore, Montes & Trujillo (2018) find that at radii between 50 and 100 kpc, the BCG+ICL system is characterized by a low metallicity. Therefore, it is unsurprising that, as a result of our larger aperture radii, we select a subsolar metallicity.

We determine the best-fit SED using a Bayesian MCMC approach, using emcee (Foreman-Mackey et al. 2013), where we treat the absolute magnitude (the EzGal normalization parameter) as a free parameter with a uniform prior. The colors generated by the EzGal model are then compared against the $g-r$ and $r-i$ colors for SDSS and either the $g-r$ and $r-i$ or $r-i$ and $i-z$ colors for DES to determine the absolute magnitude that yields the best fit for our observations. We use different colors depending on the redshift of the data because model degeneracies can become a problem in the $g-r$ colors at $z > 0.35$. Additionally, we note that had we chosen either a different IMF or formation redshift, the only impact on our results would be a uniform shift in the stellar mass values, which would only impact the value, and not the evolution in the $\alpha$ parameter.

Based on a comparison of 61 clusters between $0.20 \le z \le 0.35$ in both SDSS and DES (with light profiles out to 50kpc), we find that the mean difference between our stellar mass measurements is 0.04 dex with a standard deviation of 0.07. Thus, the difference between the 50 kpc SED-inferred BCG stellar masses from two independent imaging surveys is statistically consistent with zero. Therefore, we find excellent agreement between the stellar masses estimated for BCGs with both DES and SDSS photometry. We also test to confirm that the *gri* and *riz* photometrically determined stellar masses are consistent, except for increased scatter when color degeneracies appear, and find a median difference of -0.001 with a standard deviation of 0.013.

## 2.4. DES Cluster Richnesses and Masses

For the low-redshift SDSS-C4 sample, we use the preliminary mass-richness relation from GM&M19, which was shown to have masses that agree with the SDSS-redMaPPer clusters to within 0.1 dex. For the SDSS-redMaPPer clusters,
we use the Simet et al. (2017) mass–richness relation, which is given by Equation (1),

$$M_{\text{halo}}/(h^{-1}M_{\odot}) = 10^{14.344}(\lambda/40)^{1.33}, \quad (1)$$

and for the DES clusters, we use the mass–richness relation from McClintock et al. (2019), which is calibrated for the DES Year 1 redMaPPer data, given by Equation (2):

$$M_{\text{halo}}/(h^{-1}M_{\odot}) = 10^{14.344}(\lambda/40)^{1.356} \left( \frac{1+z_{\text{red}}}{1+0.35} \right)^{-0.30} . \quad (2)$$

For both equations, $M_{\text{halo}}$ is $M_{200m}$, $z_{\text{red}}$ is the redMaPPer photometric redshift, and $\lambda$ is the redMaPPer-defined richness. We note that both the Simet et al. (2017) and McClintock et al. (2019) mass–richness relations have intrinsic scatters associated with the halo mass at fixed richness. While not shown in Equations (1) or (2), we account for this scatter in our Bayesian MCMC analysis, as discussed in Sections 2.6 and 3. The primary difference between the two redMaPPer mass–richness relations is the redshift evolution parameter incorporated into the McClintock et al. (2019) version using the DES-redMaPPer Year 1 data. We note that we are actually using DESY3 data, and that a preliminary analysis of the DES-redMaPPer Y3 richness estimate is consistent when compared to the Y1 analysis. We also compare the redMaPPer richnesses for the 61 clusters that are in both SDSS and DES and find excellent agreement between the two measurements, with a difference of $0.9 \pm 10.1$. This translates to an offset of $0.01 \pm 0.11$ dex in units of $\log_{10} M_{\odot}$ in halo mass.

### 2.5. $m_{gap}$

For the low-redshift SDSS-C4 clusters, we use all available spectroscopic information to identify the four brightest cluster galaxies and measure $m_{gap}$. We use a radius of $0.5 \times R_{200}$, where $R_{200}$ is the radius where the mass density reaches $200 \times$ the critical density. For the SDSS-C4 clusters, $R_{200}$ is determined from the cluster masses (see GM&M18 for details).

For the redMaPPer clusters in SDSS and DES, spectroscopic completeness is too low to be of any value for membership criteria. Therefore, we use the redMaPPer red sequence–based galaxy membership criteria to define the $m_{gap}$ values. As discussed in GM&M19, we use all galaxies with a redMaPPer $P_{\text{mem}} \ge 0.984$ to define the fourth-brightest galaxy within half the virial radius. This high $P_{\text{mem}}$ value was chosen because it yielded a match in the red-sequence color parameter space density between a sample of clusters identified in both SDSS-C4 and SDSS-redMaPPer. Additionally, our final sample requires all clusters to have four or more members (including the BCG) within $0.5 R_{200}$, which we approximate using Equation (3) (Rykoff et al. 2014; McClintock et al. 2019):

$$R_{200} \approx 1.5R_c(\lambda), \quad (3)$$

where $\lambda$ is the redMaPPer richness and $R_c$ is the redMaPPer cutoff radius, given by Equation (4):

$$R_c(\lambda) = 1.0h^{-1} Mpc(\lambda/100)^{0.2}. \quad (4)$$

Using this estimate for $R_{200}$, we define M14 as the difference between the r-band apparent model magnitude of the fourth-brightest cluster member with $P_{\text{mem}} \ge 0.984$ within $0.5 R_{200}$ and the BCG’s 50 kpc r-band apparent magnitude. We note that

this differs from our definition in GM&M19, which used the apparent model magnitude of the BCG. However, the choice of the 50 kpc mag ensures consistency between the DES and SDSS BCG magnitudes, as previously shown in Figure 4, and the $m_{gap}$ values.

We use the same sample of 61 clusters as before to compare the M14 measurements between the DES and SDSS data, and find $\Delta M_{14} = -0.08 \pm 0.5$, which is consistent with zero. We note that there is significant scatter in the data when comparing the $m_{gap}$ measurements between DES and SDSS, which likely results from different galaxies being identified as the fourth-brightest cluster member.

### 2.6. Statistical Uncertainties

For each of the observables used in the analysis, we show that there is no statistically significant difference in the measurements between the observables for the different catalogs and the survey data. Recall that we use two surveys (SDSS and DES) and three cluster catalogs to make a combined sample that can be analyzed with and without redshift binning over the redshift range $0.03 \le z \le 0.6$.

We will also incorporate estimates for the uncertainties on the BCG stellar mass, cluster (or halo) mass, and M14 in the Bayesian inference of the SMHM relation. Therefore, just as we need to ensure that differences in the measurements do not introduce systematic evolution in the SMHM relation, we also need to ensure that no such biases arise from the uncertainties on the measurements. To address these uncertainties, we take a similar approach as in GM&M19, where we analyze a redshift bin that has data in both the SDSS and the DES surveys.

GM&M19 used a combined analysis of the SDSS-C4 and SDSS-redMaPPer data to ensure that the statistical uncertainties were similar in the redshift bin $0.08 \le z \le 0.12$, where both catalogs have data. To infer the uncertainties on the observables, they conducted a constrained Bayesian analysis by subsampling the SDSS-redMaPPer clusters to have the same mass distribution as the SDSS-C4 sample (and over the same redshift) and using strong priors on the four main parameters that describe the richness-dependent SMHM relation, the offset $\alpha$, the slope $\beta$, the stretch $\gamma$, and the intrinsic scatter $\sigma_{\text{int}}$ in the multivariate linear relation:

$$\log_{10}(M_{*,BCG}/(M_{\odot}/h^2)) = \alpha + \beta \times \log_{10}(M_{\text{halo}}/(M_{\odot}/h)) + \gamma \times M_{14}. \quad (5)$$

Each BCG stellar mass is treated as a draw from a normal distribution, with the mean defined by the above equation and an intrinsic scatter (standard deviation) defined as $\sigma_{\text{int}}$. The priors used for $\alpha$, $\beta$, $\gamma$, and $\sigma_{\text{int}}$ were also defined as normal distributions, with means and variances from the analysis on the SDSS-C4 clusters (GM&M18). Before running the Bayesian analysis, we shift the data by the difference between the minimum and maximum of the stellar mass and halo mass ($x_{\text{pivot}} = 14.65$ and $y_{\text{pivot}} = 11.50$). Doing this subtraction removes covariance between $\alpha$ and $\beta$. This is a well-known and established technique for constraining scaling law parameters (e.g., Rozo et al. 2014; Simet et al. 2017; McClintock et al. 2019). We then treat the uncertainties on the SDSS-redMaPPer observables as free parameters and regress for their values.
There are numerous advantages to using the SDSS-C4 clusters as the initial rung in the redshift ladder. First, the low-redshift SDSS-C4 data include cluster masses that can be inferred both from dynamics (caustic halo masses) and from a mass-richness relation that provides a self-consistent estimate of the mass uncertainties on the richness-based masses (see GM&M19 for more details). Second, the high spectroscopic completeness of the low-redshift clusters minimizes (or eliminates) foreground/background contamination in the $m_{gap}$ measurement. Third, there exist multiple simulation-based mock galaxy samples that mimic the SDSS main galaxy sample. These mock galaxies allow for alternative estimates on the uncertainties of the BCG stellar masses, membership, $m_{gap}$, and cluster masses. Moreover, by using the results from the SDSS-C4 catalog as the initial rung, we ensure consistency between our measured values, in particular $\sigma_{int}$ and those of prior studies. In GM&M18, we found excellent agreement between our $\sigma_{int}$ when $m_{gap}$ was not accounted for and other prior results (e.g., Tinker et al. 2017; Pillepich et al. 2018a; Kravtsov et al. 2018). This consistency was further reproduced in GM&M19 (see Table 2), which leads to our choice here.

We further remind the reader that, in this analysis, we are calibrating the uncertainty associated with our measurements, in particular the stellar mass, based on the spectroscopic SDSS-C4 sample (Miller et al. 2005). We note that any changes to the uncertainty associated with these stellar mass estimates would lead to a disagreement with prior estimates for $\sigma_{int}$, as described above. Moreover, we reiterate that the primary goal of this analysis is not to present absolute values of $\alpha, \beta, \gamma$, and $\sigma_{int}$, but rather to determine how these parameters change over cosmic time.

We follow the same procedure described above for the SDSS data to calibrate the uncertainties for the DES BCG stellar masses and $m_{gap}$ values. The SDSS-redMaPPer and DES-redMaPPer samples overlap within ~100 square degrees and in the redshift range $0.206 < z < 0.30$. We then create a subsample of DES-redMaPPer clusters that matches the redshift range of the SDSS-redMaPPer clusters. We label these two cluster samples SDSS-Calibration and DES-Calibration in the Appendix.

We note that we have remeasured the SDSS-redMaPPer BCG stellar masses using a 50 kpc aperture. Therefore, we first conduct a Bayesian regression analysis on the SDSS-Calibration data to infer the SMHM parameters in Equation (5). We use the same uncertainties as given in GM&M19 and follow the same algorithm described there. However, there are differences in the data, which is why we conduct this analysis again. Besides the smaller radius used to estimate the stellar mass, the SDSS-redMaPPer BCG sample is larger than the one that was used in GM&M19, since more BCGs have light profiles measured out to 50 kpc than to 100 kpc. We find that the fitted parameters of the SDSS-Calibration sample are nearly identical to those from GM&M19, and provide the results of this analysis in the Appendix.

We expect differences in the uncertainties between the SDSS and DES data for two reasons. First, the DES data are of much higher quality and depth, leading to a higher signal-to-noise ratio at a fixed aperture in the light profiles (on average). Second, the deeper DES data should make it more difficult to identify the fourth-brightest galaxy (on average), which is often located in the central region of these higher-redshift clusters. The former should lead to a reduced uncertainty on the stellar

masses (relative to SDSS), and the latter should lead to a larger scatter in the $m_{gap}$ measurements.

We conduct the simplified (no evolution) Bayesian analysis on the DES-Calibration sample and use the SDSS-Calibration posterior distributions for $\alpha, \beta, \gamma$, and $\sigma_{int}$ as strong priors. We regress for the mean errors associated with $m_{gap}$ and the BCG stellar masses in the DES data. We find that the DES stellar mass uncertainty is best fit to a value of 0.06 dex, which is smaller than the uncertainty on the SDSS-redMaPPer BCGs (0.08 dex), consistent with expectations. We also find that the uncertainties associated with $m_{gap}$ have gone up compared to the SDSS clusters, to 0.31 (from 0.15 mag). This is likely due to the DES photometry identifying more objects in the cluster core (see Figure 2) than SDSS-redMaPPer, as well as the lower spectroscopic completeness of the DES training set. We treat the halo mass errors as a fixed value identical to what was found for the SDSS BCGs in GM&M19 (0.087 dex) because the halo masses are determined from identical mass-richness relations and because we find no differences in the richness measures for the overlap sample (see Section 2.4). We note that, as discussed in GM&M19, our chosen value, 0.087 dex, was determined as the result of a joint analysis, where the parameters for the SMHM relation were simultaneously determined for a sample with halo masses estimated by both richness and the caustic phase-space technique. This value in scatter in halo mass at fixed richness corresponds to 0.20 when a natural log scale is used instead of a $\log_{10}$ scale. Therefore, our measured uncertainty has excellent agreement with the results presented in Rozo et al. (2015), who find the value to be between 0.17 and 0.21.

We do one final test to ensure that our Bayesian-inferred uncertainties on $m_{gap}$ and the BCG stellar masses in the DES-redMaPPer sample are sensible. We fix those uncertainties to the values inferred by the prior analysis, and we let the SMHM relation parameters $\beta, \gamma$, and $\sigma_{int}$ be free in the Bayesian regression. We note that we fix $\alpha$ to the value measured from the SDSS-Calibration, our remeasurement of the results from GM&M19. This is because of the strong degeneracy between $\alpha$ and $\gamma$. The results are presented in rows 2 and 3 of the Appendix. We also report the results from GM&M19 in row 1 of the Appendix. We note that the original results from GM&M19 use a slightly different pivot point for the data. However, despite this difference, we find $1\sigma$ agreement between the GM&M19 measurements and our SDSS and DES-Calibration measurements for all parameters except $\alpha$, which is impacted most significantly by the offset and is a $1.5\sigma$ difference.

In summary, our goal in this subsection has been to calibrate the uncertainties on the observables in a fixed redshift bin. By ensuring this, and also ensuring that the mean values of the observables are identical where we expect them to be (e.g., Sections 2.2, 2.3, 2.4, and 2.5), we avoid the possibility of inferring evolution where there is none. After quantifying these uncertainties, we test the model to ensure that our conclusions remain unchanged from GM&M19, which is the case. Next, we combine the SDSS and DES data into a single catalog to search for evolution in the parameters of the SMHM relation.

## 2.7. DES Final Sample

We analyze how the SMHM relation evolves with redshift using two approaches. First, we divide our data, incorporating the SDSS-C4, SDSS-redMaPPer, and DES-redMaPPer data, into
eight bins, sorted by redshift, and measure our Bayesian MCMC posteriors for each bin with the redshift evolution parameters set to 0.0. The redshift range of those bins is given in the Appendix and shown visually in Figure 1. Second, we incorporate redshift evolution using four additional parameters, described in Section 3, and fit them over the SDSS and DES clusters. We note that we include all of the clusters considered in GM&M19, not just those of the final sample (i.e., we include those clusters that were previously removed as a result of our completeness analyses, since the low-redshift analysis is redone here). In total, this homogeneous multisurvey data set covers the redshift range $0.03 \le z \le 0.60$. We note that there is spatial overlap between each of these surveys. To account for this, we remove galaxies that appear in multiple surveys. We keep the data from DES over the data from SDSS since the photometry is deeper, and we keep SDSS-C4 over SDSS-redMaPPer because of the more stringent red-sequence cluster membership.

Following this initial selection criterion, our sample consists of 1172 SDSS clusters and 1564 DES clusters with halo masses greater than $10^{14} M_{\odot} h^{-1}$, which yields a total sample of 2736 clusters. This number does not account for further halo mass limits; however, as in GM&M19, we expect this total data sample to have differing halo mass lower limits as we move to higher redshifts. We also check for $m_{gap}$ incompleteness, since both SDSS and DES are flux-limited surveys.

For each bin, as in GM&M18 and GM&M19, we apply an $m_{gap}$ completeness criterion, based on the binning of the BCG and fourth-brightest galaxy's absolute magnitudes against the BCG's apparent magnitude and $m_{gap}$, to determine the apparent magnitude limit of the sample (a redshift-dependent limit; Colless 1989; Garilli et al. 1999; La Barbera et al. 2010; Trevisan et al. 2017; GM&M18). Additionally, since the halo mass distribution can be approximated as Gaussian, the peak indicates the mass where the sample becomes incomplete. However, we apply a lower halo mass cut, located at the halo mass where the amplitude of the binned halo mass distribution decreases to 70% of the peak value, to ensure high completeness out to higher redshifts. This halo mass criterion is conservative, and results in a redMaPPer richness threshold of ≈22, well above the detection limit. However, when combined with the $m_{gap}$ completeness analysis, these cuts reduce our available sample down to 2323 clusters, a reduction of ~15.1%. Slightly more restrictive (higher) halo mass lower limits do not impact our final results. Of the 2323 clusters in our final sample, 1062 come from SDSS and 1261 come from DES.

Figure 5 visualizes our final sample and shows the 50 kpc stellar masses versus the halo masses, color coded by M14. We also show the error bars on each set of survey data (see Section 2.6). Figure 5 includes the entire final sample (following all halo mass and $m_{gap}$ completeness cuts), as described, and spans the redshift range $0.03 < z < 0.60$, as shown in Figure 1. Since the dependence on M14 is evident, it is clear that the previously detected stellar mass–M14 stratification continues to exist at higher redshifts.

## 2.8. Simulated Data

For the simulated analysis of this study, we examine the evolution of the SMHM-$m_{gap}$ relation using the magnetohydrodynamical cosmological galaxy formation simulations from the IllustrisTNG<fnref n="55" /> suite (Pillepich et al. 2018b; Springel et al. 2018;

<figure><img src="image_8.png" /><figcaption>**Figure 5.** SMHM–M14 relation for the SDSS-C4, SDSS-redMaPPer, and DES-redMaPPer samples, colored via M14. We see that a stellar mass-$m_{gap}$ stratification continues to persist when measured out to high redshifts. The black cross represents the error in halo mass, 0.087 dex, and stellar mass, 0.08 dex, for the SDSS data; and the red cross represents the error in halo mass, 0.087 dex, and stellar mass, 0.06 dex, for the DES data. The black square represents the pivot point that is used in our Bayesian analysis.</figcaption></figure>

Nelson et al. 2019). Specifically, we use the IllustrisTNG 300-1 simulation (Pillepich et al. 2018a), with snapshots analyzed at redshifts of 0.08, 0.11, 0.15, 0.18, 0.24, 0.31, 0.40, and 0.52 (snapshots 92, 90, 87, 85, 81, 77, 72, and 66), the redshifts that best match the median of our eight binned samples given in the Appendix. For this comparison, we identify the 260 clusters with $log10(M_{200m}/(M_{\odot}/h)) > 14.0$ in the redshift 0.0 snapshot, and use these same clusters in the higher-redshift bins.

For each simulation box, we use the 3D information provided directly from each snapshot of the TNG 300-1 simulation, including the $M_{200m}$ halo masses, measured within $R_{200m}$; the galaxy positions, x, y, and z; $R_{200m}$; and the magnitudes. Cluster membership is determined using positional information (x, y, z) and a fit to the red sequence, such that cluster member candidates within $0.5 R_{200m}$, and within $2\sigma$ from the red sequence, are identified as members.

For the observed data, we use a 50 kpc aperture on the deblended BCG (see Sections 2.2 and 2.3). In the simulation, the parent halo and subhalo have been identified using Subfind (Springel et al. 2001). We therefore use a fixed 50 kpc physical aperture for the BCG subhalo alone to calculate the BCG stellar masses. In other words, we use the Subfind deblending algorithm to separate the stellar components of the BCG from the halo satellite galaxies. An example of one such projected image of a TNG 300-1 BCG is shown in Figure 6.

The use of a simulation allows us to quantify the impact of deblending in a controlled fashion, since the 3D positions of the satellites are known. In Figure 7, we show that all but four BCGs have a close companion within the 50 kpc aperture. The peak of the distribution is around three or four satellites per BCG, with some BCGs being in very crowded but very localized environments. Proper deblending is therefore necessary to study the BCG itself. However, we note that as discussed in Zhang et al. (2019a), these satellites are not included in the ICL that is present and included in our measurements of the BCG stellar masses, in both the data and the simulations. However, we note that stellar particles that may be observationally identified as part of the ICL, but in simulations are associated with the satellite galaxies, are thus not incorporated into our stellar mass measurement.

<footnote n="55">⁵⁵ http://www.tng-project.org/</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 6.** An example of a TNG 300-1 image showing the stellar particle information of the halo, which shows that there are many satellite galaxies near to the BCG. The primary image is the inner 500 kpc centered on the BCG. The insert shows the stellar particle information for just the inner 100 kpc of the subhalo containing the BCG. The circle represents the 50 kpc aperture.</figcaption></figure>

We use the snapshot information for each BCG subhalo in each of the eight redshift snapshots. For M14, we follow a similar procedure to what is done in GM&M19, where M14 is the difference between the r-band model magnitude of the fourth-brightest member and the r-band model magnitude of the TNG 300-1 BCG. We do not apply a completeness criterion to our simulated data. However, to make our approaches between the data and the simulation as homogeneous as possible, we apply a halo mass completeness limit, in the same manner as described in Section 2.7, which accounts for the fact that the 260 most massive clusters at $z = 0.0$ are not guaranteed to be the most massive at each of the higher-redshift snapshots, because each dark matter halo follows a unique accretion history.

## 3. Bayesian MCMC Model

We use a similar hierarchical Bayesian MCMC approach to what is described in GM&M19 to determine the values of $\alpha, \beta, \gamma$, and $\sigma_{\text{int}}$ given in Equation (5). Any changes to the underlying model are designed to improve the efficiency of our analysis. The Bayesian formalism works by convolving the prior information for a selected model with the likelihood of the observations given that model, to yield the probability of observing the data for the model—the posterior distribution up to a normalization constant called the Bayesian evidence.

To determine the posterior distributions for each parameter in the SMHM relation, our MCMC model generates values for the observed stellar mass, halo mass, and $m_{\text{gap}}$ values at each step of our likelihood analysis, which are directly compared to our observed measurements.

### 3.1. The Observed Quantities

We model the $\log_{10}$ BCG stellar masses (y), $\log_{10}$ halo masses (x), and M14 values (z) as normal distributions with

<figure><img src="image_10.png" /><figcaption>**Figure 7.** The distribution of the number of satellite galaxies within 50 kpc of the TNG 300-1 BCGs, which highlights that either masking or using the particle information directly is necessary to estimate the stellar masses of the TNG 300-1 BCGs.</figcaption></figure>

mean values (locations) taken from our observed/measured results. The standard deviations associated with each point are taken from the uncertainties on each measurement, which are determined using the sample of clusters observed by SDSS and DES, and include an estimate of the observational uncertainty ($\sigma_{x_0}, \sigma_{y_0}, \sigma_{z_0}$), as well as a stochastic component from a beta function, $\beta(0.5, 100)$ (GM&M18), which allows for additional uncertainty on the observational errors. These errors are treated statistically in the Bayesian model as free nuisance parameters $\sigma_x, \sigma_y$, and $\sigma_z$.

### 3.2. The Unobserved Quantities

Our aim is to constrain the parameters of the SMHM relation: the offset, slope, stretch, and intrinsic scatter ($\alpha, \beta, \gamma$, and $\sigma_{\text{int}}$) as a function of redshift. In GM&M19, we modeled the evolution of these parameters as power laws $(1+z)^{n_i}$, where $n_i$ defined the amount and shape of the redshift evolution for each of the four model parameters, $\alpha, \beta, \gamma$, and $\sigma_{\text{int}}$.

In this work, we extend the data from the limit of $z \sim 0.3$ to $z \sim 0.6$ using the deeper DES data. We initially tried the same power-law parameterization as in GM&M19. However, we found that a simple power law could not accurately characterize the flattening of the slope at $z > 0.3$ to a single value, while simultaneously having the sharp increase in the slope at low redshift. We explored numerous functional forms and found that the Pareto function, given in Equation (6), has the correct shape over the range of data explored in this work:

$$f(a) = \begin{cases} \text{Constant} + n_i (a_m^{n_i} / a^{n_i+1}) & a \ge a_m \\ f(a_m) & a < a_m \end{cases} \quad (6)$$

This is a Pareto Type I distribution, which is characterized by a scale parameter $a_m$ and a shape parameter $n_i$, which in our case defines the strength of the evolution for $\alpha, \beta, \gamma$, and $\sigma_{\text{int}}$.

A second change from GM&M19 is that we evolve against lookback time as opposed to each cluster's photometric redshift, given that from an astrophysical perspective, stellar evolution is characterized by time as opposed to the universe's scale factor. We choose $a_m = 0.4$ Gyr, which is a redshift of $\sim 0.03$, the lowest redshift our data probe. Below this lower limit, the Pareto distribution is a constant fixed to the value at its lowest data point. The Pareto distribution also asymptotes to
<table><caption>**Table 1**

Bayesian Analysis Parameters for the Combined SDSS-C4, SDSS-redMaPPer, and DES-redMaPPer Sample</caption><thead><tr><th>Symbol</th><th>Description</th><th>Prior</th></tr></thead><tbody><tr><td>α<sub>0</sub></td><td>The offset of the SMHM relation</td><td>U(-20,20)</td></tr><tr><td>β<sub>0</sub></td><td>The high-mass power-law slope</td><td>Linear Regression Prior</td></tr><tr><td>γ<sub>0</sub></td><td>The stretch parameter, which describes the stellar mass–M14 stratification</td><td>Linear Regression Prior</td></tr><tr><td>σ<sub>int0</sub></td><td>The uncertainty in the intrinsic stellar mass at fixed halo mass</td><td>U(0.0, 0.5)</td></tr><tr><td>y<sub>i</sub></td><td>The underlying distribution in stellar mass</td><td>Equation (7)</td></tr><tr><td>x<sub>i</sub></td><td>The underlying halo mass distribution</td><td>N(14.23,0.18<sup>2</sup>)</td></tr><tr><td>z<sub>i</sub></td><td>The underlying m<sub>gap</sub> distribution</td><td>N(2.51,0.62<sup>2</sup>)</td></tr><tr><td>n<sub>α</sub></td><td>The shape parameter associated with the evolution of α</td><td>U(-5.0, 5.0)</td></tr><tr><td>n<sub>β</sub></td><td>The shape parameter associated with the evolution of β</td><td>U(-5.0, 5.0)</td></tr><tr><td>n<sub>γ</sub></td><td>The shape parameter associated with the evolution of γ</td><td>U(-5.0, 5.0)</td></tr><tr><td>n<sub>σ</sub></td><td>The power law associated with the evolution of σ<sub>int</sub></td><td>U(-5.0, 5.0)</td></tr><tr><td>σ<sub>y0i</sub></td><td>The uncertainty between the observed stellar mass and intrinsic stellar mass distribution</td><td>0.08 or 0.06 dex</td></tr><tr><td>σ<sub>x0i</sub></td><td>The uncertainty associated with the mass–richness relation</td><td>0.087 dex</td></tr><tr><td>σ<sub>z0i</sub></td><td>The uncertainty between the underlying and observed m<sub>gap</sub> distribution</td><td>0.15 or 0.31</td></tr></tbody></table>

Note. U(a, b) refers to a uniform distribution, where a and b are the upper and lower limits. The linear regression prior is of the form $-1.5 \times \log(1 + \text{value}^2)$. N(a, b) refers to a normal distribution with mean and variance of a and b. Additionally, we note that for $x_i$ and $z_i$, the means and widths given in this table are example values belonging to the lowest redshift bin.

a constant at large a (high redshift or lookback time), which we treat as a free nuisance parameter in the analysis. In other words, this is the constant in Equation (6).

Using the Pareto distribution, we model the cluster portion of the SMHM relation linearly as:

$$y_i = \frac{\alpha_0 + ((n_\alpha) \times \frac{(0.4^{n_\alpha})}{(t)^{1+n_\alpha}})}{\underbrace{\text{evolving offset, } \alpha(t)}_{\text{evolving slope, } \beta(t)}} + \frac{(\beta_0 + ((n_\beta) \times \frac{(0.4^{n_\beta})}{(t)^{1+n_\beta}})}{\underbrace{\text{evolving offset, } \alpha(t)}_{\text{evolving slope, } \beta(t)}} \times x_i + \frac{(\gamma_0 + ((n_\gamma) \times \frac{(0.4^{n_\gamma})}{(t)^{1+n_\gamma}})}{\underbrace{\text{evolving offset, } \alpha(t)}_{\text{evolving slope, } \beta(t)}} \times z_i, \quad (7)$$

where x, y, and z are the observed halo masses, BCG stellar masses, and m<sub>gap</sub> values, and t is the lookback time, calculated

Equation (7), the values of α<sub>0</sub>, β<sub>0</sub>, γ<sub>0</sub>, and σ<sub>int0</sub> represent the parameter values at large t, when the shape is flat—the maximum redshift of the sample. This is different than the α, β, γ, and σ<sub>int</sub> parameters in GM&M19, which represent the values at z = 0.0. Therefore, these two sets of parameters cannot be directly compared unless the same model is used in both (either a power-law or a Pareto function), which is done and discussed in Section 4.

This Bayesian model regresses the generated values against the observed stellar mass, halo mass, and m<sub>gap</sub> values simultaneously and self-consistently. The parameters that model the underlying distributions and their uncertainties are nuisance parameters, and thus are marginalized over when we present the posterior distributions. Each parameter in the Bayesian analysis, along with its prior information, is presented in Table 1. We express the entire posterior as:

$$p(\alpha_0, \beta_0, \gamma_0, \sigma_{\text{int0}}, n_\alpha, n_\beta, n_\gamma, n_\sigma, x_i, z_i, \sigma_{y_i}, \sigma_{x_i}, \sigma_{z_i} | x, y, z) \propto \frac{\sum_i P(y_{0i} | \alpha_0, \beta_0, \gamma_0, \sigma_{y_i}, n_\alpha, n_\beta, n_\gamma, n_\sigma, \sigma_{\text{int}}, x_i, z_i) P(x_{0i} | x_i, \sigma_{x_i}) P(z_{0i} | z_i, \sigma_{z_i})}{p(x_i) p(z_i) p(\sigma_{x_i}) p(\sigma_{y_i}) p(\sigma_{z_i}) p(\alpha) p(\beta) p(\gamma) p(\sigma_{\text{int}}) p(n_\alpha, n_\beta, n_\gamma, n_\sigma)}, \quad (8)$$

using the photometric redshift. α<sub>0</sub>, β<sub>0</sub>, and γ<sub>0</sub> are free parameter offsets that are asymptotic at high lookback time (high redshift). Note that for zero evolution, Equation (7) reverts to Equation (5). We also assume a Gaussian likelihood form, with σ<sub>int</sub>(t) evolving with redshift: σ<sub>int0</sub> + ((n<sub>σ</sub>) × (0.4<sup>n<sub>σ</sub></sup>)/(t)<sup>1+n<sub>σ</sub></sup>). n<sub>α</sub>, n<sub>β</sub>, n<sub>γ</sub>, and n<sub>σ</sub> measure the redshift evolution of α, β, γ, and σ<sub>int</sub>, respectively.

This model is nested. Thus, for the redshift-binned samples, these n parameters are set to 0.0, and, as in GM&M19, the zero-redshift evolution model from GM&M18 is returned. This approach allows us to interpret how much better a given model is (e.g., with redshift evolution versus without) using only the posterior distribution. We also note that in

where each ith cluster is a component in the summed log likelihood.

This is a hierarchical Bayes model because the priors on the true halo masses (x<sub>i</sub>) and M14 values (z<sub>i</sub>) depend on the models themselves (the observed halo mass and M14 distributions).

## 4. Results

We evaluate the strength of the redshift evolution in the SMHM–M14 relation using our previously described MCMC model (Section 3), Bayesian formalism, and linear SMHM relation (Equation (7)). For this analysis, we have run the MCMC chains to convergence by examining the parameter
<figure><img src="image_11.png" /><figcaption>**Figure 8.** SMHM–M14 parameter posteriors from Equation (7). The posterior distribution for $\alpha_0$, $\beta_0$, $\gamma_0$, $\sigma_{\text{int0}}$, $n_\alpha$, $n_\beta$, $n_\gamma$, and $n_\sigma$. As in GM&M18, we see that $\gamma$ is significantly nonzero and $\sigma_{\text{int0}}$ is less than 0.1 dex. However, we note that as a result of the modified redshift evolution form given by Equation (7), the values are not directly comparable to the results from GM&M19. Instead, the posteriors for $\alpha_0$, $\beta_0$, $\gamma_0$, and $\sigma_{\text{int0}}$ represent the values that these parameters asymptote to. To see the values at the redshifts measured in our study, see Figure 9. Using this model, $n_\alpha$ is 5.8$\sigma$ from 0.0, $n_\beta$ is 3.5$\sigma$ from 0.0, and $n_\gamma$ is 4.0$\sigma$ from 0.0.</figcaption></figure>

autocorrelation functions. We run our analysis for 10 million steps with a burn-in of 2 million steps. The triangle plot, Figure 8, shows the 1D and 2D posterior distributions for the eight SMHM relation parameters, $\alpha_0$, $\beta_0$, $\gamma_0$, $\sigma_{\text{int0}}$, $n_\alpha$, $n_\beta$, $n_\gamma$, and $n_\sigma$.

A negative (positive) value for the evolution parameters (the “n”s) indicates that the parameter itself ($\alpha(t)$, $\beta(t)$, $\gamma(t)$, or $\sigma_{\text{int}}(t)$) is growing (shrinking) with increasing lookback time. The 1D marginalized posteriors in Figure 8, given in row 3 of

Table 1, illustrate the evidence for evolution in the SMHM relation in the offset ($n_\alpha$), slope ($n_\beta$), and stretch ($n_\gamma$). We find no evidence for evolution in the intrinsic scatter ($n_{\sigma_{\text{int}}}$), which is small, well within $2\sigma$ of 0.0, and consistent with the prior results presented in GM&M18 and GM&M19.

The parameter fits and their errors are provided in Table 2. That table starts with constraints from a revised analysis of the data used in GM&M19. Recall that GM&M19 used a simple power-law evolution model and 100 kpc apertures for the BCG
<table><caption>**Table 2**

Posterior Distribution Results with Lookback Time Evolution</caption><thead><tr><th>Sample</th><th>α<sub>0</sub></th><th>β<sub>0</sub></th><th>γ<sub>0</sub></th><th>σ<sub>int0</sub></th><th>n<sub>α</sub></th><th>n<sub>β</sub></th><th>n<sub>γ</sub></th><th>n<sub>σ</sub></th></tr></thead><tbody><tr><td>GM&M19<sup>a</sup></td><td>-0.325<sup>+0.030</sup><sub>-0.029</sub></td><td>0.224<sup>+0.020</sup><sub>-0.021</sub></td><td>0.142<sup>+0.016</sup><sub>-0.015</sub></td><td>0.078<sup>+0.010</sup><sub>-0.009</sub></td><td>0.094<sup>+0.068</sup><sub>-0.061</sub></td><td>0.573<sup>+0.152</sup><sub>-0.141</sub></td><td>0.018<sup>+0.026</sup><sub>-0.025</sub></td><td>0.015<sup>+0.018</sup><sub>-0.017</sub></td></tr><tr><td>GM&M19<sup>b</sup></td><td>-0.267<sup>+0.032</sup><sub>-0.028</sub></td><td>0.268<sup>+0.024</sup><sub>-0.023</sub></td><td>0.155 ± 0.016</td><td>0.081 ± 0.009</td><td>0.037<sup>+0.056</sup><sub>-0.057</sub></td><td>0.307<sup>+0.105</sup><sub>-0.095</sub></td><td>0.003<sup>+0.027</sup><sub>-0.025</sub></td><td>0.002<sup>+0.016</sup><sub>-0.016</sub></td></tr><tr><td>This paper</td><td>-0.069<sup>+0.028</sup><sub>-0.027</sub></td><td>0.281<sup>+0.018</sup><sub>-0.017</sub></td><td>0.095<sup>+0.008</sup><sub>-0.007</sub></td><td>0.085 ± 0.005</td><td>-0.196<sup>+0.034</sup><sub>-0.032</sub></td><td>0.263<sup>+0.086</sup><sub>-0.075</sub></td><td>0.084 ± 0.021</td><td>-0.016 ± 0.011</td></tr><tr><td>w/o M14</td><td>0.154<sup>+0.015</sup><sub>-0.013</sub></td><td>0.241<sup>+0.030</sup><sub>-0.027</sub></td><td>...</td><td>0.098 ± 0.005</td><td>-0.199 ± 0.029</td><td>-0.026<sup>+0.057</sup><sub>-0.055</sub></td><td>...</td><td>0.013<sup>+0.013</sup><sub>-0.012</sub></td></tr><tr><td>z > 0.09<sup>c</sup></td><td>0.284<sup>+0.086</sup><sub>-0.113</sub></td><td>0.267<sup>+0.025</sup><sub>-0.023</sub></td><td>0.076<sup>+0.07</sup><sub>-0.006</sub></td><td>0.084 ± 0.006</td><td>-0.705<sup>+0.118</sup><sub>-0.077</sub></td><td>0.260<sup>+0.080</sup><sub>-0.079</sub></td><td>0.149<sup>+0.016</sup><sub>-0.021</sub></td><td>-0.014 ± 0.016</td></tr></tbody></table>

Notes. The Equation (7) parameter fits for the SDSS data in GM&M19 ($z < 0.3$) compared to the fits in this paper that use DES data to extend the analysis to $z = 0.6$.

a The same data from GM&M19 was reanalyzed using the model from this paper (e.g., Equation (8)) and the original 100 kpc apertures.

b The same data from GM&M19 was reanalyzed using the model from this paper (e.g., Equation (8)) and 50 kpc apertures for a fair comparison to the results from the data in this paper.

c All data in this analysis, except the lowest redshift bin. We note that while the posterior results differ, when plotted in the redshift range of interest, we find $1\sigma$ agreement, as shown in Figure 9.

stellar masses. In addition, the GM&M19 data only extend to $z = 0.3$. Notable differences between GM&M19 and this work include the use of the Pareto function to describe the evolution, smaller 50 kpc aperture stellar masses, and higher redshift data out to $z = 0.6$. Therefore, we reanalyze the GM&M19 data using this new model (Equation (7)) for that original data set, as well as with the new model and 50 kpc apertures for the SDSS BCG stellar masses.

The first row can be compared to the original GM&M19 discovery of evolution in the slope of the SMHM relation, which was reported at the $3.5\sigma$ statistical level. Switching from the power-law fitting function to the Pareto function, we find that the evolution in the slope is also significantly nonzero, $n_β = 0.573^{+0.152}_{-0.141}$ or $\sim 3.8\sigma$. We find that by using the Pareto function, the confidence in the detection of the slope evolution has gone up, likely because the Pareto function more closely matches the shape of the data as a function of lookback time. In the second row, we compare to what happens as we take a smaller physical radius to measure the stellar mass, and find no statistically significant differences between the parameter values as a result of using a smaller aperture.

The third row in Table 2 shows the parameter fits for the data described in this paper, which uses DES to extend the analysis to $z = 0.6$. The slope evolution is now detected at $3.5\sigma$. Most of the error bars on the parameters have also decreased compared to the analysis on the SDSS data alone. Although we have nearly doubled the sample size from GM&M19, the evolution in $\beta(t)$ is in the late universe, where DES does not provide new data. However, the DES data is useful for pinning down the amplitude of the flattening of the tail of the Pareto function for the parameters. We note that we have a similar amount of DES data in our higher-redshift bins compared to what is present in the SDSS data in GM&M19 in their highest redshift bin (see Figure 1). Therefore, we do not expect a significant drop in the error bars on the inferred parameters when moving from the SDSS-only data set to the combined SDSS and DES data set.

If the parameters $\alpha(t)$, $\beta(t)$, $\gamma(t)$, and $\sigma_{\text{int}}(t)$ are evolving between $z = 0.3$ and $z = 0.6$, we would expect to see differences in the zero-points $\alpha_0$, $\beta_0$, $\gamma_0$, and $\sigma_{\text{int0}}$ between the second and third rows, because these parameters represent the values after the Pareto function flattens out to a constant at the upper limit of the redshift traced by the data, which is deeper for DES than for SDSS. We do not detect any changes in $\beta_0$ or $\sigma_{\text{int0}}$ after we extend the analysis to $z = 0.6$ using the DES data. However, we do find that $\alpha_0$ is significantly higher and $\gamma_0$ is significantly lower as we extend the data from $z = 0.3$

to $z = 0.6$. In fact, we detect evolution in the offset $\alpha$ at $5.8\sigma$ and evolution in the stretch $\gamma$ at $4.0\sigma$.

We also consider what would happen if we excluded $m_{\text{gap}}$ from our model by dropping $z_i$ in Equation (7). We find that the significance of the evolution of the slope drops from $n_β = 0.263^{+0.086}_{-0.075}$ to being statistically consistent with zero (the fourth row in Table 2). Therefore, as originally noted in GM&M19, the detection of the evolution of the slope of the linear SMHM relation requires the use of $m_{\text{gap}}$ in the analysis.

It is interesting that the offset $\alpha$ still shows statistically significant evolution when we ignore $m_{\text{gap}}$ in the analysis. We note that in Equation (7), the offset parameter ($\alpha_0$) is not a direct measure of the amplitude of the SMHM. Even when incorporating $n_\alpha$, the first term $\alpha(t)$ does not quantify the amplitude of the SMHM (because of the inclusion of $m_{\text{gap}}$). However, without M14, $\alpha(t)$ is simply the overall amplitude of the SMHM relation as a function of lookback time, which is characterized by $\alpha_0$ (at the redshift limit of the data) and $n_\alpha$. Thus, when we ignore $m_{\text{gap}}$ in our analysis, it appears that we are detecting significant evolution in the amplitude of the SMHM to $z = 0.6$. However, assuming that $\alpha(t)$ traces the evolution in the amplitude, the sign on the evolution $n_\alpha$ would imply that BCGs are getting more massive as we look back in time. This of course cannot be the case, and we explain this and how to best interpret the observed evolution of $\alpha(t)$ in the next subsection.

## 4.1. Comparison to Binned Results

By using Equation (7), we can characterize the evolution of the SMHM relation through the parameters $n_\alpha$, $n_\beta$, $n_\gamma$, and $n_{\sigma_{\text{int}}}$ simultaneously over the full-redshift parameter space. However, we can also apply a strong prior and set those parameters to zero, which reverts Equation (7) to Equation (5). We can then infer $\alpha(t = t_{\text{med}})$, $\beta(t = t_{\text{med}})$, $\gamma(t = t_{\text{med}})$ and $\sigma_{\text{int}}(t = t_{\text{med}})$, on data separated into discrete redshift bins, as shown in Figure 1, where $t = t_{\text{med}}$ is the median lookback time of the BCGs in each predefined bin. This allows us to make a direct comparison between the fit and a timed step evolution of the SMHM relation. In the binned analysis, we assume no evolution within the upper and lower limits on the redshift, which is likely to be true if the time intervals are small enough. We note that the full analysis of Equation (7) is the correct statistical analysis, since it does not require that assumption and because it does not require a somewhat arbitrary choice of binning. However, the binned analysis provides a good cross-check.
<figure><img src="image_12.png" /><figcaption>**Figure 9.** The effective offset, slope, stretch, and intrinsic scatter from Equation (7) as a function of lookback time. The green lines represent the result of the fit to the full equation using all of the data. The green error bands are the total error in each parameter as a function of lookback time. The brown dashed line represents the median posterior when we fit our model without the lowest redshift data bin ($z < 0.09$). The points represent the redshift-binned data when the evolution parameters (e.g., $n_\alpha$) are fixed to zero. The error bars contain the middle 67% of the 1D marginalized posterior. The SDSS and DES clusters are shown in purple and the simulation-based TNG 300-1 data are shown in blue.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 10.** The combined offset of $\alpha(t) + \gamma(t) \times \text{M14}$. This is for M14 = 2.0. This figure illustrates that the combined value of these parameters does not evolve. This indicates that we are not seeing any evolution in the overall amplitude of the SMHM relation, as characterized by the combination of the offset term ($\alpha(t)$) and the stretch term ($\gamma(t)$).</figcaption></figure>

In Figure 9, we compare the SMHM parameter values from Equation (7) with the evolution parameters fixed to zero (the purple dots) against the fully evolving parameters for the offset $\alpha(t)$, slope $\beta(t)$, stretch $\gamma(t)$, and intrinsic scatter $\sigma_{\text{int}}(t)$ (the green line and green band). The purple error bars are the $1\sigma$ error bars on the binned posteriors. The green error band incorporates the error on the parameter zero-points (e.g., $\alpha_0$, $\beta_0$, $\gamma_0$, and $\sigma_{\text{int0}}$), as well as the error on the corresponding evolution component ($n_\alpha$, $n_\beta$, $n_\gamma$, and $n_{\sigma_{\text{int}}}$). We find good agreement between the two separate analyses: binned and unbinned.

We first note the evolution in the slope parameter $\beta(t)$ of the SMHM, which is quantified via the Pareto function as $n_\beta = 0.263_{-0.075}^{+0.086}$ and is evident in Figure 9. The Pareto function does a good job of capturing the shape of the evolution, which is changing fast at low redshift, as originally reported in GM&M19. The slope of the SMHM relation becomes roughly constant beyond a lookback time of 3 Gyrs, corresponding to $z = 0.245$. However, in recent times, we clearly identify a steepening of the slope of the SMHM relation for massive clusters.

Figure 9 shows no evidence for evolution in the intrinsic scatter looking back 6 Gyrs ($z = 0.6$). The value of the scatter is $\sigma_{\text{int0}} = 0.085 \pm 0.005$. This is the same low value for intrinsic scatter found in GM&M18 and GM&M19, except extending to $z = 0.6$. We note that an outlier at a lookback time of ~4.3 Gyr exists in the binned analysis. We are unable to explain this feature of the data, which is in a bin containing exclusively DES data.

In contrast to the intrinsic scatter, Figure 9 shows statistically significant evidence for evolution in both the offset and stretch parameter over the last 6 Gyrs. As for the slope, the Pareto function captures the shape of the observed binned evolution, which shows a gradual increase (for $\alpha$) or decrease (for $\gamma$). Thus, we are clearly identifying an increase in the offset and a decrease in the stretch parameter at higher redshifts.

As discussed at the end of the previous section, despite detecting evolution in the offset, $\alpha$, we are not actually detecting evolution in the overall amplitude (or median stellar mass at fixed halo mass and $m_{\text{gap}}$) for the SMHM relation. As we show in Figure 10, we plot the combination of the offset
<figure><img src="image_14.png" /><figcaption>**Figure 11.** For the first, lowest redshift bin, we show the contours representative of the 2D posterior distribution for α and γ in purple. A strong covariance exists between these parameters. In green, we overlay the posterior distribution for these same parameters, as estimated from Figure 9, for the median lookback time of bin 1.</figcaption></figure>

and stretch terms, given mathematically as $α + γ × M_{14}$, as a function of lookback time, and detect no discernible evolution in the amplitude of the SMHM to $z = 0.6$ for a fixed M14 value of 2.0.

## 4.2. Statistical Correlations

In the following subsections, we address correlations in the inferred parameters, as well as in the observables.

### 4.2.1. Parameter Correlations

Figure 8 shows some interesting structures in the 2D posteriors. Besides the obvious correlation between the parameters and their corresponding evolution (i.e., α and $n_α$), there is also a weak correlation between the slope and the offset. We note that this correlation has been minimized by recentering the data using a pivot point selected to be the midway value of the extreme values of the observables.

More importantly, we note the correlation between the offset and the stretch parameters (α and γ), intertwined with their evolution parameters ($n_α$ and $n_γ$). This was also seen in GM&M19, but that analysis lacked the redshift depth to study the consequences of this correlation. In this work, we have enough data over a large enough lookback time to bin the data beyond where the evolution of the slope flattens.

In Figure 9, we notice that α(t) and γ(t) in the lowest redshift bin are ~1σ low (α(t)) and ~1σ high (γ(t)) when compared to unbinned fits. While some other bins have similar differences, this is the only bin where the binned values do not follow the general trend displayed by the green posterior distributions (even though the measured values are within 1σ). We explain this discrepancy via the covariance between α and γ, evident in Figure 8. As the parameter α scatters low in the MCMC sampling, γ scatters high. There is a clear degeneracy between these two parameters. We show this degeneracy just for clusters in the lowest redshift bin in Figure 11. We overplot the 2D posterior distribution between α(t) and γ(t) for the lowest redshift bin (in green) and the total posterior (in purple). Figure 11 highlights that we see the 2D posterior error ellipses associated with the single-bin analysis overlap with those measured based on the entire sample. However, we note that

there is likely a weak covariance between the two sets of 2D posteriors that may be responsible for part of this agreement.

Figure 11 exemplifies why one may want to avoid binning in this type of SMHM analysis, since unaccounted-for parameter covariances can lead to incorrect fits to the evolving SMHM relation. Our fully unbinned analysis and our hierarchical Bayesian formalism allow for these covariances to naturally be accounted for in the fitted parameters and their marginalized posteriors. However, we note that despite this reservation, the results from the binned analysis are largely consistent with our evolution analysis. For an additional test, we fit our model to the data after excluding those clusters in the lowest redshift bin ($z < 0.09$). We show the median posterior for α(t), β(t), γ(t), and $σ_{int}(t)$ as the brown dashed line in Figure 4. The entire posterior is given in line 5 of Table 1. While the median posterior values differ (likely due to the covariances between the parameters and their evolution), when plotted as a function of lookback time, there are no significant differences. We note that, as is evident from Figure 11, a higher stretch parameter at low redshift is still preferred.

### 4.2.2. Data Correlations

We make some assumptions in Equation (7) and in Section 2. Primarily, we assume that the observables (stellar mass, halos mass, and $m_{gap}$) are independent observables. If our data were strongly correlated to each other in some complicated way (or to some additional latent parameter), we would need to quantify those correlations and their impacts on the fitted SMHM parameters. Our main concern is that unlike GM&M18, we use richness as a proxy for halo mass, and a correlation could exist between richness and either $m_{gap}$ or the BCG stellar masses, which would affect our conclusions.

Hearin et al. (2013) reported a correlation between M12 and the cluster richness at fixed halo mass. While they do not quantify the correlation, they suggest that there is evidence that having a large $m_{gap}$ is correlated with being under-rich at a given halo mass. At fixed X-ray luminosity (as a proxy for halo mass), Erfanianfar et al. (2019) report a weak and positive correlation between the cluster richness and the stellar mass of BCGs (a Pearson correlation coefficient of ~0.2). Furnell et al. (2018) used dynamical masses to find a similar weakly positive correlation between richness and BCG stellar mass (a Spearman rank correlation coefficient of $r_s = 0.137$). On the other hand, Farahi et al. (2020) used the IllustrisTNG simulations to find a moderate anticorrelation between richness and BCG stellar mass at fixed halo mass (a Pearson correlation coefficient of ~−0.4). None of these reported correlations are strong or consistent. Compounding the issue is the fact that the richnesses, halo masses, stellar masses, and $m_{gap}$ measurements are not homogeneously measured in either the data or the simulations.

Given the above information, and the fact that we do not have other halo mass proxies like X-ray luminosity, weak lensing, or dynamical masses for our clusters, there is little we can do in terms of a precision exploration of the data correlations. However, we can calculate the Pearson correlation coefficient between our richness-inferred halo masses and $m_{gap}$ at fixed BCG stellar mass. If we fix the stellar mass to within $11.4 < \log_{10}(M_*/(M_\odot h^{-2})) < 11.6$, the stellar mass range that allows us to measure the correlation across the entire redshift range, we find a moderate anticorrelation of ~−0.4.
The statistical significance of correlation coefficients is not well defined. Most of the literature uses some form of jackknife sampling (e.g., Erfanianfar et al. 2019). Here, however, we use a Bayesian-like approach, where we apply Equation (5) to forward-model our data using the uncertainties given in Section 2.6. We can then apply a correlation between $m_{gap}$ and the halo mass, before the simulated observational uncertainties are incorporated. We then run 10,000 simulations with and without the correlation, and measure the standard deviation on the measured correlation coefficient as well as the probability of the correlation coefficient being observed in a purely noncorrelated data set (i.e., a null test). We find that the error on the correlation is ~0.04, and the probability of a purely randomized data set, showing the same level of correlation, we find to be $p=0.001$. We conclude that correlation between M14 and $M_{200m}$ (as inferred by richness) is significant. We note that we find a nearly identical anticorrelation between M14 and $M_{200m}$ in the IllustrisTNG sample (-0.36).

We can use this same forward-modeling technique to quantify the effect that this correlation could be having on the parameters we measure when assuming independence. We note that this is not the same as developing a new statistical model that incorporates correlations between the data, which we reserve for a future effort. Using this Bayesian-like approach, we do, however, estimate the impact of this correlation on the slope. We find that the correlation between M14 and $M_{200m}$ results in an increase in the slope by approximately 0.15. However, we note that because this correlation persists across the redshift space, we do not believe that it impacts our detected redshift evolution, but rather just the measured value of the slope. Thus, this analysis provides us with a good idea of the level of the impact of the correlation on the slope, offset, and stretch, without introducing one or more new free parameters to the model. We will explore these interesting correlations in a future analysis.

### 4.3. Comparison to IllustrisTNG 300-1

Figure 9 offers direct comparisons between the observed (SDSS and DES) results and the simulated TNG 300-1 measurements. Such a comparison allows us to understand whether the physical prescriptions built into the TNG 300-1 simulated universe yield observations that match those found in the observed universe. This analysis is designed to yield a fair comparison, since for both data sets we subtract off the same median in stellar and halo mass, which is based on the observed SDSS and DES values, allowing the posterior distributions to be directly compared. In such a comparison, $\alpha$ is related to the median stellar mass (at a halo mass of $\log_{10}(M_{halo}/(M_{\odot}/h)) = 14.65$) at a given $m_{gap}$ (since the $\gamma$ values agree). We note that for the simulated data, as for the observational data, $M_{halo}$ refers to $M_{200m}$.

The only similar result between the observed and simulated universes is the lack of evolution of $\sigma_{int}$; no evolution is detected in either. Interestingly, Pillepich et al. (2018a) detect modest evolution in $\sigma_{int}$, such that from $z=0.0$ to $z=0.5$ the value increases by ~0.04 dex. However, the results presented in Pillepich et al. (2018a) do not account for $m_{gap}$. In contrast, when we measure the evolution in the SMHM relation without incorporating $m_{gap}$, we do not find this evolution, though the size of our error bars may prevent us from detecting it.

One of the more significant results from using our approach to measure the 50 kpc mag for the TNG 300-1 data is the

absence of noticeable evolution in the slope of the SMHM relation for TNG 300-1. In our observed data set, late-time growth appears to occur primarily in the last 2 billion years; however, in the TNG 300-1 simulation, there is no detectable evolution over the entire time range analyzed. However, we note that the absence of redshift evolution in the slope with the TNG 300-1 data agrees with Pillepich et al. (2018a) and Engler et al. (2021), who claim no such evolution. Thus, unlike for observations, where GM&M19 found that the incorporation of $m_{gap}$ led to the detection of evolution, for the TNG 300-1 simulation this is not the case.

Another difference between our TNG 300-1 and the prior measurements from Pillepich et al. (2018a) and Engler et al. (2021) is the value of $\beta$. We measure a value of approximately 0.42 for the slope when the stellar mass is measured within 50 kpc when $m_{gap}$ is incorporated, and 0.48 when it is not. Our estimate is therefore in agreement with Pillepich et al. (2018a), who measure the stellar content within 30 kpc and find a slope of 0.49 (no error bars reported). We note that our slope is much shallower than that measured in Engler et al. (2021) and other slopes measured in Pillepich et al. (2018a), ≈0.70, which are measured using two times the stellar half-mass-radius, a radius far greater than the 50 kpc aperture we use. Therefore, we can conclude, as shown in GM&M19, that had we used a large aperture ($> = 100$ kpc) to measure the BCG stellar mass and magnitude, we would likely recover a steeper slope. One additional note is that both here and in GM&M19, we find that the slope of the SMHM relation is steeper when $m_{gap}$ is incorporated, which serves as evidence that incorporating information about the satellite galaxy population (via $m_{gap}$) yields a steeper slope than the traditional SMHM relation, in agreement with the general conclusions from Tinker et al. (2021). However, as shown in the Appendix, this trend is not shown in the TNG 300-1 data. Instead, we see that the slopes are within $1\sigma$, which may serve as the first bit of evidence that the BCGs and growth prescriptions in the TNG 300-1 simulation are overdominant.

The remaining two parameters $\alpha$ and $\gamma$ are also dramatically different. Unlike in our observational data, there is little evidence of any evolution in $\alpha$ or $\gamma$ out to high redshift. Given that $\alpha$ and $\gamma$ are covariant, it is unsurprising that if one of these two parameters shows no evolution, the other parameter also shows no evolution; and, as discussed in Section 5, this is likely related to the growth prescription used in TNG. Additionally, the values for $\alpha$ also differ significantly. At first glance, it appears as though the TNG 300-1 BCGs are undermassive. However, that is not the case. A more valid comparison would use the value of $\alpha + \gamma \times z_{med}$, which is representative of the median stellar mass of the BCGs at a given halo mass. This comparison suggests that the TNG systems are approximately 0.06 dex overmassive. It is important to note, when doing such a comparison, that the median M14 values for TNG are approximately 1–1.5 mag greater than the observed values; unlike in the observed universe, low $m_{gap}$ systems ($M_{14} < 1.0$) do not exist. While our measurements suggest that part of this difference is a result of the slightly overmassive BCGs, for such a scenario to occur, it is likely that the merging prescription used in TNG 300-1 also results in poorly populated red sequences, such that few intermediate-brightness galaxies exist, yielding substantially fainter fourth-brightest galaxies.
<figure><img src="image_15.png" /><figcaption>**Figure 12.** We display two sets of distributions for the low-z (unfilled) and high-z (filled) data. For both, we show two $m_{gap}$ regimes, the upper 80%-90% regime (orange and red) and the lower 10%-20% regime (purple and blue). The shaded regions represent the posterior distributions from our Bayesian model. For the low-redshift data, as shown, we see a steeper slope and more pronounced stratification, which results from a larger stretch.</figcaption></figure>

## 5. Discussion

In GM&M19, we introduced the novel observation of evolution in the slope of the SMHM relation, and used that observation to offer insight into the late-time hierarchical growth of BCGs. As shown here, by significantly expanding the parameter space out to higher redshifts/earlier lookback times using DES-redMaPPer data, we reach a much deeper understanding of how BCGs and the clusters that they reside within have grown and evolved over the last 6 billion years.

Currently, a clear consensus between observations, simulations, and models about how BCGs grow over this redshift range does not exist. Using semi-analytic models, researchers have found that at late times ($0.0 < z < 0.5$) BCG's grow by a factor of ≈1.5–2.0 (De Lucia & Blaizot 2007; Guo et al. 2011; Shankar et al. 2015). In contrast, observations suggest that over this redshift range much of the growth occurs in the BCG's outermost envelope, incorporating regimes that are often characterized as being part of the ICL (van Dokkum et al. 2010; Burke et al. 2015; Huang et al. 2018; Furnell et al. 2021), which highlights the necessity of looking at the BCG+ICL system jointly. However, in this work, we use the additional information provided via the inclusion of $m_{gap}$ into the SMHM relation to determine physically what growth is occurring in the BCG+ICL system over this redshift range.

In this work, we extend the redshift evolution of the cluster-scale SMHM presented in GM&M19 ($0.03 < z < 0.30$) out to $z_{red} = 0.6$. To briefly summarize our findings, we confirm all of the key results from GM&M18 and GM&M19: $m_{gap}$ is definitively a latent parameter within the SMHM relation; incorporating $\gamma$ and M14 into the SMHM relation reduces $\sigma_{int}$, and accounting for $m_{gap}$ yields significant evolution in the slope of the SMHM relation over late time. From this analysis, for the first time, we report evolution in both the $\alpha$ and $\gamma$ parameters, which represent the offset and stretch, respectively.

It is this observed evolution that drives our understanding of how BCGs evolve.

To understand how the stellar mass, halo mass, and $m_{gap}$ are changing as a function of lookback time (or redshift), in Figure 12 we plot the SMHM relation data for a low-redshift sample (the second and third bins) and a high-redshift sample (the seventh and eighth bins). We note that due to the lack of data and the larger difference in parameter values, we do not use the lowest redshift bin. For each sample, we plot the data in the 10th–20th and 80th–90th percentiles of the $m_{gap}$ distribution. This is shown by the filled (high-z) and unfilled (low-z) data points. We then overlay the results of the posterior distributions shown in Figures 8 and 9 as the shaded regions. Figure 12 highlights a few of our key findings. First, $\gamma$ significantly grows as one moves forward in lookback time, as evidenced by how much larger is the separation between the two shaded regions at low-z when compared to high-z. We note that if $\gamma$ were not evolving, the separations between the high and low M14 bins would not be growing, regardless of the change in the halo mass distribution of the data, since $\gamma$ does not vary with $M_{halo}$. Second, $\beta$ is growing as one moves forward in lookback time. Third, $\sigma_{int}$, the spread in the data at fixed $m_{gap}$, is unchanged between these two distributions, which supports our measurement that $\sigma_{int}$ is not evolving. Fourth, the most insightful observation shown here, as highlighted by the regions in the $M_{halo}$ distribution where these data sets overlap ($14.4 < \log_{10}(M_{halo}/(M_{\odot}/h)) < 14.7$), the BCG stellar mass distribution remains the same, and thus the BCG stellar mass within 50 kpc is not growing. This is also supported by the constant value of $\alpha + \gamma \times M_{14}$ given in Figure 10.

In GM&M19, given the absence of evolution in $\gamma$, we assumed that any growth observed was due entirely to growth in the BCG. However, as shown in Figure 12, the stellar mass within 50 kpc is not growing over this redshift range. This observation highlights that the driver behind all of the evolution we detect and have previously detected is instead $m_{gap}$. First, with respect to the slope, if the stellar mass is not changing, then the only way for the slope to increase would be for the stellar masses of the distributions of clusters that are linked by having similar $m_{gap}$ values to change, as a result of changing $m_{gap}$ values. For the slope to increase, this would likely be in such a manner that the most massive systems, with the more massive BCGs, have $m_{gap}$ values that are growing more efficiently and quickly, likely since these BCGs reside in richer clusters.

Recall that $m_{gap}$ is the difference in brightness between the BCG and the fourth-brightest cluster member (Dariush et al. 2010), and results from the hierarchical assembly of the BCG (GM&M18), such that we expect clusters characterized by larger $m_{gap}$ values to form earlier. Since the observed evolution results from changes in the $m_{gap}$ distribution, the greatest insight into what is physically happening can instead be gleaned from the evolution in $\gamma$. For $\gamma$ to evolve, the $m_{gap}$ distribution must be changing with time. This is not happening in a manner that changes the BCG stellar mass (within 50kpc). Instead, therefore, what is likely happening is that there are mergers between the bright satellite galaxies, and the BCG deposit the stellar material at radii beyond 50 kpc, what we interpret as the ICL. Therefore, the outer envelopes contain all of the recent BCG growth, and it is only through the incorporation of $m_{gap}$ that we are able to detect this evolution
without measuring the BCG + ICL profile, as done in Zhang et al. (2019a). As a result of this scenario, the separation of clusters with fixed $m_{gap}$ values, what we refer to as our stratification becomes larger, while the stellar mass distribution (within 50 kpc) remains fixed. Therefore, the incorporation of $m_{gap}$ has elucidated that BCGs continue to grow hierarchically in this redshift range, but all of the added stellar material goes directly into the growth of the ICL. This result is supported observationally by Furnell et al. (2021), who find evidence of ICL growth over $0.1 < z < 0.5$.

While the main takeaways of this paper are observational, we do want to comment on what the absence of evolution in TNG 300-1 means. Since we detect no evolution in either the slope, stretch, or offset parameters, clearly the same kind of hierarchical growth prescription is not occurring within the TNG 300-1 simulation. Additionally, the TNG 300-1 clusters are characterized by larger $m_{gap}$ measurements. Therefore, it is possible that the majority of the stellar mass within these BCGs is assembled at earlier times. Moreover, due to an overefficient merger process, there exists an absence of fainter satellite galaxies in the TNG 300-1 simulation, the same population that we observationally find must be responsible for the continued hierarchical assembly of the BCG + ICL systems.

In this work, we have focused on the late-time evolution of the SMHM relation out to $z \sim 0.6$. As shown here, $\beta$ shows significant late-time evolution, predominantly over the redshift range $0.0 < z < 0.15$, and for the first time we detect statistically significant evolution in $\alpha$ and $\gamma$, which clarifies that this evolution is driven by BCG hierarchical growth that is evident not in the stellar mass, but rather in the $m_{gap}$. We are left with a few paths forward. If we choose to further tighten the constraints on this late-time evolution, we must either incorporate more large, statistically complete samples of low-redshift clusters, $z < 0.1$ (there are fewer than 200 SDSS low-z clusters compared to $\sim 1300$ DES high-z clusters), which are difficult to obtain; or we can forge ahead to higher redshifts to determine whether these parameters continue to evolve out to $z = 1.0$, using a data set such as the DES-ACT overlap (Hilton et al. 2021) or the DES-South Pole Telescope overlap (Bleem et al. 2015, 2020; Huang et al. 2020), an approach that faces similar observational and modeling challenges as the results presented here, but presents the opportunity for us to further quantify and better constrain this evolution. Additionally, given that we have now statistically verified that the stellar mass-$m_{gap}$ trend exists both in observations and state-of-the-art hydrodynamic simulations, although we note that the evolution trends do not match, a key step forward may be to determine the physical meaning of this correlation between stellar mass and $m_{gap}$, what it may reveal to us about the formation history of the BCG and its host cluster dark matter halo, and quantify how the stellar mass, halo mass, and $m_{gap}$ parameter space map to a cluster's formation redshift. Lastly, as explored in GM&M19, we should continue to study the BCG light profiles out to large radii of 100 kpc and beyond. Another vital step forward as part of that effort will be to take advantage of the ICL measurements done by Zhang et al. (2019a) for the DES clusters, to determine whether we are able to detect significant growth in the ICL over this redshift range. Such a result would verify our conclusions that all of the recent growth is contained within the ICL, and that these recent mergers change the $m_{gap}$ distribution and yield our detected evolution.

We would like to thank the anonymous referee for their useful comments that helped to improve this paper. Additionally, this paper has gone through an internal review by the DES collaboration. J.G.M. would like to thank Emmet Golden-Marx for useful discussions about generating the photometric images used in this paper and Ying Zu for key discussions about the cosmological measurements that were used in this analysis. J.G.M. acknowledges support from the National Key Basic Research and Development Program of China (No. 2018YFA0404504) and the National Science Foundation of China (No. 11873038, 11890692).

Funding for the DES Projects has been provided by the U.S. Department of Energy, the U.S. National Science Foundation, the Ministry of Science and Education of Spain, the Science and Technology Facilities Council of the United Kingdom, the Higher Education Funding Council for England, the National Center for Supercomputing Applications at the University of Illinois at Urbana-Champaign, the Kavli Institute of Cosmological Physics at the University of Chicago, the Center for Cosmology and Astro-Particle Physics at the Ohio State University, the Mitchell Institute for Fundamental Physics and Astronomy at Texas A&M University, Financiadora de Estudos e Projetos, Fundação Carlos Chagas Filho de Amparo à Pesquisa do Estado do Rio de Janeiro, Conselho Nacional de Desenvolvimento Científico e Tecnológico and the Ministério da Ciência, Tecnologia e Inovação, the Deutsche Forschungsgemeinschaft, and the Collaborating Institutions in the Dark Energy Survey.

The Collaborating Institutions are Argonne National Laboratory, the University of California at Santa Cruz, the University of Cambridge, Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas-Madrid, the University of Chicago, University College London, the DES-Brazil Consortium, the University of Edinburgh, the Eidgenössische Technische Hochschule (ETH) Zürich, Fermi National Accelerator Laboratory, the University of Illinois at Urbana-Champaign, the Institut de Ciències de l'Espai (IEEC/CSIC), the Institut de Física d'Altes Energies, Lawrence Berkeley National Laboratory, the Ludwig-Maximilians Universität München and the associated Excellence Cluster Universe, the University of Michigan, NSF's NOIRLab, the University of Nottingham, The Ohio State University, the University of Pennsylvania, the University of Portsmouth, SLAC National Accelerator Laboratory, Stanford University, the University of Sussex, Texas A&M University, and the OzDES Membership Consortium.

Based in part on observations at the 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 National Science Foundation under Grant Numbers 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
<table><caption>**Table 3**

Posterior Distribution Results</caption><thead><tr><th>Data</th><th>z<sub>min</sub></th><th>z<sub>max</sub></th><th>z<sub>med</sub></th><th>log<sub>10</sub>(M<sub>halo<sub>min</sub></sub>)</th><th>n<sub>clusters</sub></th><th>α(z = z<sub>med</sub>)</th><th>β(z = z<sub>med</sub>)</th><th>γ(z = z<sub>med</sub>)</th><th>σ<sub>int</sub>(z = z<sub>med</sub>)</th></tr></thead><tbody><tr><td>GM&amp;M19 bin 4 (100 kpc)</td><td>0.208</td><td>0.300</td><td>0.247</td><td>14.39</td><td>210</td><td>-0.34 ± 0.03</td><td>0.32 ± 0.06</td><td>0.150 ± 0.013</td><td>0.082 ± 0.009</td></tr><tr><td>SDSS-Calibration</td><td>0.206</td><td>0.300</td><td>0.242</td><td>14.38</td><td>234</td><td>-0.26 ± 0.03</td><td>0.36 ± 0.05</td><td>0.153 ± 0.014</td><td>0.081 ± 0.008</td></tr><tr><td>DES-Calibration</td><td>0.206</td><td>0.300</td><td>0.243</td><td>14.24</td><td>351</td><td>-0.26</td><td>0.38 ± 0.04</td><td>0.159 ± 0.004</td><td>0.067 ± 0.006</td></tr><tr><td>Bin 1</td><td>0.030</td><td>0.090</td><td>0.075</td><td>14.03</td><td>112</td><td>-0.22 ± 0.05</td><td>0.52 ± 0.08</td><td>0.135 ± 0.020</td><td>0.070 ± 0.014</td></tr><tr><td>Bin 2</td><td>0.090</td><td>0.130</td><td>0.112</td><td>14.02</td><td>203</td><td>-0.25 ± 0.03</td><td>0.43 ± 0.06</td><td>0.168 ± 0.015</td><td>0.083 ± 0.009</td></tr><tr><td>Bin 3</td><td>0.130</td><td>0.170</td><td>0.151</td><td>14.03</td><td>289</td><td>-0.22 ± 0.03</td><td>0.45 ± 0.04</td><td>0.160 ± 0.012</td><td>0.075 ± 0.008</td></tr><tr><td>Bin 4</td><td>0.170</td><td>0.210</td><td>0.187</td><td>14.17</td><td>260</td><td>-0.25 ± 0.03</td><td>0.35 ± 0.05</td><td>0.159 ± 0.013</td><td>0.082 ± 0.008</td></tr><tr><td>Bin 5</td><td>0.210</td><td>0.270</td><td>0.236</td><td>14.28</td><td>404</td><td>-0.19 ± 0.02</td><td>0.36 ± 0.04</td><td>0.131 ± 0.011</td><td>0.079 ± 0.006</td></tr><tr><td>Bin 6</td><td>0.270</td><td>0.360</td><td>0.307</td><td>14.24</td><td>385</td><td>-0.18 ± 0.02</td><td>0.34 ± 0.04</td><td>0.120 ± 0.011</td><td>0.089 ± 0.006</td></tr><tr><td>Bin 7</td><td>0.360</td><td>0.470</td><td>0.407</td><td>14.35</td><td>317</td><td>-0.11 ± 0.03</td><td>0.33 ± 0.05</td><td>0.101 ± 0.011</td><td>0.053 ± 0.006</td></tr><tr><td>Bin 8</td><td>0.470</td><td>0.600</td><td>0.528</td><td>14.35</td><td>353</td><td>-0.08 ± 0.03</td><td>0.33 ± 0.06</td><td>0.097 ± 0.012</td><td>0.092 ± 0.006</td></tr><tr><td>TNG300-1</td><td>0.08</td><td>0.08</td><td>0.08</td><td>13.93</td><td>241</td><td>-0.37 ± 0.04</td><td>0.44 ± 0.03</td><td>0.132 ± 0.011</td><td>0.108 ± 0.005</td></tr><tr><td>TNG300-1 no γ</td><td>0.08</td><td>0.08</td><td>0.08</td><td>13.93</td><td>241</td><td>0.13 ± 0.02</td><td>0.51 ± 0.04</td><td>0.138 ± 0.006</td><td>0.138 ± 0.006</td></tr><tr><td>TNG300-1</td><td>0.11</td><td>0.11</td><td>0.11</td><td>13.92</td><td>238</td><td>-0.39 ± 0.05</td><td>0.43 ± 0.03</td><td>0.138 ± 0.012</td><td>0.100 ± 0.005</td></tr><tr><td>TNG300-1 no γ</td><td>0.11</td><td>0.11</td><td>0.11</td><td>13.92</td><td>238</td><td>0.12 ± 0.02</td><td>0.47 ± 0.04</td><td>0.126 ± 0.006</td><td>0.126 ± 0.006</td></tr><tr><td>TNG300-1</td><td>0.15</td><td>0.15</td><td>0.15</td><td>13.90</td><td>236</td><td>-0.34 ± 0.04</td><td>0.45 ± 0.03</td><td>0.127 ± 0.010</td><td>0.102 ± 0.005</td></tr><tr><td>TNG300-1 no γ</td><td>0.15</td><td>0.15</td><td>0.15</td><td>13.90</td><td>236</td><td>0.12 ± 0.02</td><td>0.47 ± 0.04</td><td>0.130 ± 0.005</td><td>0.130 ± 0.005</td></tr><tr><td>TNG300-1</td><td>0.18</td><td>0.18</td><td>0.18</td><td>13.91</td><td>235</td><td>-0.32 ± 0.04</td><td>0.44 ± 0.03</td><td>0.120 ± 0.011</td><td>0.102 ± 0.005</td></tr><tr><td>TNG300-1 no γ</td><td>0.18</td><td>0.18</td><td>0.18</td><td>13.91</td><td>235</td><td>0.12 ± 0.02</td><td>0.47 ± 0.04</td><td>0.126 ± 0.006</td><td>0.126 ± 0.006</td></tr><tr><td>TNG300-1</td><td>0.24</td><td>0.24</td><td>0.24</td><td>13.87</td><td>234</td><td>-0.37 ± 0.04</td><td>0.46 ± 0.03</td><td>0.141 ± 0.010</td><td>0.094 ± 0.004</td></tr><tr><td>TNG300-1 no γ</td><td>0.24</td><td>0.24</td><td>0.24</td><td>13.87</td><td>234</td><td>0.12 ± 0.02</td><td>0.46 ± 0.04</td><td>0.125 ± 0.006</td><td>0.125 ± 0.006</td></tr><tr><td>TNG300-1</td><td>0.31</td><td>0.31</td><td>0.31</td><td>13.84</td><td>233</td><td>-0.34 ± 0.04</td><td>0.44 ± 0.03</td><td>0.135 ± 0.011</td><td>0.106 ± 0.005</td></tr><tr><td>TNG300-1 no γ</td><td>0.31</td><td>0.31</td><td>0.31</td><td>13.84</td><td>233</td><td>0.12 ± 0.02</td><td>0.45 ± 0.04</td><td>0.137 ± 0.006</td><td>0.137 ± 0.006</td></tr><tr><td>TNG300-1</td><td>0.40</td><td>0.40</td><td>0.40</td><td>13.80</td><td>231</td><td>-0.29 ± 0.04</td><td>0.48 ± 0.03</td><td>0.130 ± 0.011</td><td>0.108 ± 0.005</td></tr><tr><td>TNG300-1 no γ</td><td>0.40</td><td>0.40</td><td>0.40</td><td>13.80</td><td>231</td><td>0.15 ± 0.03</td><td>0.47 ± 0.04</td><td>0.137 ± 0.006</td><td>0.137 ± 0.006</td></tr><tr><td>TNG300-1</td><td>0.52</td><td>0.52</td><td>0.52</td><td>13.71</td><td>236</td><td>-0.25 ± 0.04</td><td>0.49 ± 0.03</td><td>0.127 ± 0.009</td><td>0.111 ± 0.005</td></tr><tr><td>TNG300-1 no γ</td><td>0.52</td><td>0.52</td><td>0.52</td><td>13.71</td><td>236</td><td>0.18 ± 0.03</td><td>0.50 ± 0.05</td><td>0.149 ± 0.007</td><td>0.149 ± 0.007</td></tr></tbody></table>

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 U.S. Department of Energy, Office of Science, Office of High Energy Physics.

## Appendix

In Table 3, found in this appendix, we provide the fits to the parameters in the binned analysis of Figure 9 (including both the observational and simulated data), along with the fits for the data subsets that were used to calibrate the uncertainties on the observables.

### ORCID iDs

* Jesse B. Golden-Marx https://orcid.org/0000-0002-6394-045X

* Y. Zhang https://orcid.org/0000-0001-5969-4631

* R. L. C. Ogando https://orcid.org/0000-0003-2120-1154

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

* T. M. C. Abbott https://orcid.org/0000-0003-1587-3931

* M. Aguena https://orcid.org/0000-0001-5679-6747

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

* D. Bacon https://orcid.org/0000-0002-2562-8537

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

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

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

* M. Carrasco Kind https://orcid.org/0000-0002-4802-3194

* F. J. Castander https://orcid.org/0000-0001-7316-4573

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

* M. Crocce https://orcid.org/0000-0002-9745-6228

* L. N. da Costa https://orcid.org/0000-0002-7731-277X

* J. De Vicente https://orcid.org/0000-0001-8318-6813

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

* H. T. Diehl https://orcid.org/0000-0002-8357-7467

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

* A. E. Evrard https://orcid.org/0000-0002-4876-956X

* I. Ferrero https://orcid.org/0000-0002-1295-1132

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

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

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

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

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

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

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

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

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

* T. Jeltema https://orcid.org/0000-0001-6089-0365

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

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

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

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

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

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

* R. Morgan https://orcid.org/0000-0002-7016-5471

* F. Paz-Chinchón https://orcid.org/0000-0003-1339-2683

* A. Pieres https://orcid.org/0000-0001-9186-6042

* A. K. Romer https://orcid.org/0000-0002-9328-879X
E. Sanchez https://orcid.org/0000-0002-9646-8198

S. Serrano https://orcid.org/0000-0002-0211-2861

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

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

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

## References

Abbott, T. M. C., Abdalla, F. B., Allam, S., et al. 2018, ApJS, 239, 18

Aguerri, J. A. L., Girardi, M., Boschin, W., et al. 2011, A&A, 527, A143

Aihara, H., Armstrong, R., Bickerton, S., et al. 2018, PASJ, 70, S8

Alam, S., Albareti, F. D., Allende Prieto, C., et al. 2015, ApJS, 219, 12

Behroozi, P. S., Wechsler, R. H., & Conroy, C. 2013, ApJ, 770, 57

Bernardi, M., Hyde, J. B., Sheth, R. K., Miller, C. J., & Nichol, R. C. 2007, AJ, 133, 1741

Bleem, L. E., Stalder, B., de Haan, T., et al. 2015, ApJS, 216, 27

Bleem, L. E., Bocquet, S., Stalder, B., et al. 2020, ApJS, 247, 25

Brough, S., Tran, K.-V., Sharp, R. G., von der Linden, A., & Couch, W. J. 2011, MNRAS, 414, L80

Bruzual, G., & Charlot, S. 2003, MNRAS, 344, 1000

Burke, C., Hilton, M., & Collins, C. 2015, MNRAS, 449, 2353

Burke, D. L., Rykoff, E. S., Allam, S., et al. 2018, AJ, 155, 41

Colless, M. 1989, MNRAS, 237, 799

Croton, D. J., Springel, V., White, S. D. M., et al. 2006, MNRAS, 365, 11

Dariush, A. A., Raychaudhury, S., Ponman, T. J., et al. 2010, MNRAS, 405, 1873

De Lucia, G., & Blaizot, J. 2007, MNRAS, 375, 2

Emsellem, E., Cappellari, M., Krajnović, D., et al. 2007, MNRAS, 379, 401

Engler, C., Pillepich, A., Joshi, G. D., et al. 2021, MNRAS, 500, 3957

Erfanianfar, G., Finoguenov, A., Furnell, K., et al. 2019, A&A, 631, A175

Farahi, A., Ho, M., & Trac, H. 2020, MNRAS, 493, 1361

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

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

Furnell, K., Collins, C., Kelvin, L., et al. 2018, MNRAS, 478, 4952

Furnell, K. E., Collins, C. A., Kelvin, L. S., et al. 2021, MNRAS, 502, 2419

Garilli, B., Maccagni, D., & Andreon, S. 1999, A&A, 342, 408

Golden-Marx, J. B., & Miller, C. J. 2018, ApJ, 860, 2

Golden-Marx, J. B., & Miller, C. J. 2019, ApJ, 878, 14

Gozaliasl, G., Finoguenov, A., Khosroshahi, H. G., et al. 2016, MNRAS, 458, 2762

Gu, M., Conroy, C., & Behroozi, P. 2016, ApJ, 833, 2

Guo, Q., White, S., Boylan-Kolchin, M., et al. 2011, MNRAS, 413, 101

Harrison, C. D., Miller, C. J., Richards, J. W., et al. 2012, ApJ, 752, 12

Hearin, A. P., Zentner, A. R., Newman, J. A., & Berlind, A. A. 2013, MNRAS, 430, 1238

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

Hoshino, H., Leauthaud, A., Lackner, C., et al. 2015, MNRAS, 452, 998

Huang, N., Bleem, L. E., Stalder, B., et al. 2020, AJ, 159, 110

Huang, S., Leauthaud, A., Greene, J. E., et al. 2018, MNRAS, 475, 3348

Jones, C., & Forman, W. 1984, ApJ, 276, 38

Jones, L. R., Ponman, T. J., & Forbes, D. A. 2000, MNRAS, 312, 139

Krajnović, D., Ural, U., Kuntschner, H., et al. 1919, A&A, 635, A129

Kravtsov, A. V., Vikhlinin, A. A., & Meshcheryakov, A. V. 2018, AstL, 44, 8

Kron, R. G. 1980, ApJS, 43, 305

La Barbera, F., de Carvalho, R. R., de La Rosa, I. G., et al. 2010, MNRAS, 408, 1313

Lauer, T. R., Postman, M., Strauss, M. A., Graves, G. J., & Chisari, N. E. 2014, ApJ, 797, 82

Lauer, T. R., Faber, S. M., Richstone, D., et al. 2007, ApJ, 662, 808

Li, T. S., DePoy, D. L., Marshall, J. L., et al. 2016, AJ, 151, 157

Lin, Y.-T., Brodwin, M., Gonzalez, A. H., et al. 2013, ApJ, 771, 61

Lin, Y.-T., & Mohr, J. J. 2004, ApJ, 617, 879

Lin, Y.-T., Hsieh, B.-C., Lin, S.-C., et al. 2017, ApJ, 851, 139

Mancone, C. L., & Gonzalez, A. H. 2012, PASP, 124, 606

Matthee, J., Schaye, J., Crain, R. A., et al. 2017, MNRAS, 465, 2381

McClintock, T., Varga, T. N., Gruen, D., et al. 2019, MNRAS, 482, 1352

McDermid, R. M., Alatalo, K., Blitz, L., et al. 2015, MNRAS, 448, 3484

Miller, C. J., Nichol, R. C., Reichart, D., et al. 2005, AJ, 130, 968

Montes, M., & Trujillo, I. 2018, MNRAS, 474, 917

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

Moster, B. P., Naab, T., & White, S. D. M. 2013, MNRAS, 428, 3121

Moster, B. P., Naab, T., & White, S. D. M. 2018, MNRAS, 477, 1822

Murante, G., Giovalli, M., Gerhard, O., et al. 2007, MNRAS, 377, 2

Nelson, D., Springel, V., Pillepich, A., et al. 2019, ComAC, 6, 2

Oliva-Altamirano, P., Brough, S., Lidman, C., et al. 2014, MNRAS, 440, 762

Oser, L., Ostriker, J. P., Naab, T., Johansson, P. H., & Burkert, A. 2010, ApJ, 725, 2312

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

Pereira, M. E. S., Palmese, A., Varga, T. N., et al. 2020, MNRAS, 498, 5450

Petrosian, V. 1976, ApJL, 210, L53

Pillepich, A., Nelson, D., Hernquist, L., et al. 2018a, MNRAS, 475, 648

Pillepich, A., Springel, V., Nelson, D., et al. 2018b, MNRAS, 473, 4077

Proctor, R. N., de Oliveira, C. M., Dupke, R., et al. 2011, MNRAS, 418, 2054

Rhee, G. F. R. N., & Latour, H. J. 1991, A&A, 243, 38

Rodríguez-Puebla, A., Primack, J. R., Avila-Reese, V., & Faber, S. M. 2017, MNRAS, 470, 651

Rozo, E., Bartlett, J. G., Evrard, A. E., & Rykoff, E. S. 2014, MNRAS, 438, 78

Rozo, E., Rykoff, E. S., Bartlett, J. G., & Melin, J.-B. 2015, MNRAS, 450, 592

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

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

Salpeter, E. E. 1955, ApJ, 121, 161

Schombert, J. M. 1986, ApJS, 60, 603

Shankar, F., Buchan, S., Rettura, A., et al. 2015, ApJ, 802, 73

Simet, M., McClintock, T., Mandelbaum, R., et al. 2017, MNRAS, 466, 3103

Solanes, J. M., Perea, J. D., Darriba, L., et al. 2016, MNRAS, 461, 321

Springel, V., White, S. D. M., Tormen, G., & Kauffmann, G. 2001, MNRAS, 328, 726

Springel, V., Pakmor, R., Pillepich, A., et al. 2018, MNRAS, 475, 676

The Dark Energy Survey Collaboration 2005, arXiv:astro-ph/0510346

Tinker, J. L. 2017, MNRAS, 467, 3533

Tinker, J. L., Cao, J., Alpaslan, M., et al. 2021, MNRAS, 505, 5370

Tinker, J. L., Brownstein, J. R., Guo, H., et al. 2017, ApJ, 839, 121

Tonini, C., Bernyk, M., Croton, D., Maraston, C., & Thomas, D. 2012, ApJ, 759, 43

Trevisan, M., Mamon, G. A., & Khosroshahi, H. G. 2017, MNRAS, 464, 4593

van Dokkum, P. G., Whitaker, K. E., Brammer, G., et al. 2010, ApJ, 709, 1018

von der Linden, A., Best, P. N., Kauffmann, G., & White, S. D. M. 2007, MNRAS, 379, 867

Zhang, Y., Miller, C., McKay, T., et al. 2016, ApJ, 816, 98

Zhang, Y., Yanny, B., Palmese, A., et al. 2019a, ApJ, 874, 165

Zhang, Y., Jeltema, T., Hollowood, D. L., et al. 2019b, MNRAS, 487, 2578

Zu, Y., & Mandelbaum, R. 2015, MNRAS, 454, 1161

Zwicky, F. 1933, AcHPh, 6, 110

Zwicky, F. 1951, PASP, 63, 61