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

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σ), and detect negative evolution in the stretch parameter (4.0σ) and positive evolution in the offset parameter (5.8σ). 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.


INTRODUCTION
The stellar mass -halo mass (SMHM) relation is a 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 /h)) ≥ 14.0), which is 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 AGN feedback in central galaxies (Kravtsov et al. 2018).The intrinsic scatter in stellar mass at fixed halo mass (σ int ) can constrain 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 evolution of the slope and scatter provide insight into how BCGs grow and evolve over cosmic time (Gu et al. 2016;Golden-Marx & Miller 2019).
As a result of "inside-out" growth, all 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 100kpc centered on the BCG (Huang et al. 2018), and that the 100kpc boundary may represent a transitional regime between the BCG's outer envelope and the ICL (Zhang et al. 2019a).Therefore, when characterizing BCG evolution associated with the parameters of the SMHM relation, it is vital to measure BCG photometry within large radii, as opposed to the more commonly used 20-30kpc aperture radii (e.g., Lin et al. 2013;Zhang et al. 2016;Lin et al. 2017) as discussed in Golden-Marx & Miller (2019), referred to as GM&M19.More specifically, including stellar mass within large radii strengthens the correlation between BCG stellar mass and halo mass (Moster et al. 2018;Golden-Marx & Miller 2019).However, to yield a stronger correlation, one must also incorporate a third parameter related to BCG growth.
One observational measurement inherently tied to BCG hierarchical growth is m gap , the difference in r−band magnitude between the BCG and 4th 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 4th 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 that of non-central galaxies, as BCGs grow hierarchically, their stellar mass and magnitude increase, while those same parameters for the 4th brightest member galaxy 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 a statistical latent parameter within the cluster SMHM relation, as first presented by Golden-Marx & Miller (2018), which from here on is 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 2nd 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 2-parameter SMHM relation is larger than the intrinsic scatter in the SMHM relation after accounting for the 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.1dex, which is smaller than previous studies by as much as a factor of two (Gu et al. 2016;Tinker et al. 2017;Zu & Mandelbaum 2015;Kravtsov et al. 2018;Pillepich et al. 2018a, e.g.,).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 BCG's 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.0from 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 Illustris TNG300-1 cosmological hydrodynamic simulation, measure little change in the slope between z=1.0 and z=0.0.In addition to the slope, the expected redshift evolution of the intrinsic scatter, σ int , has also been investigated in models and simulations (Matthee et al. 2017;Pillepich et al. 2018a;Gu et al. 2016).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 DES clusters using the parameter µ * , the sum of the individual galaxy stellar masses weighted by their membership probability.Palmese et al. (2020) find no evolution in the correlation between µ * and richness (λ).In contrast, Pereira et al. (2020) explicitly accounts for redshift evolution in their relation and finds weak redshift evolution.However, Pereira et al. (2020) note that their evolution is within the accepted uncertainty of the total stellar mass (≈ 0.1 dex).Therefore, like for the SMHM relation for BCGs, it is currently unclear how the 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 ≈0.20dex or 40%.To expand upon those results, a primary goal of this paper is to characterize the evolution of the cluster SMHM out to z ∼ 0.6.
To constrain 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. 2015a;Abbott et al. 2018).We chose 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≈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 Hyper Suprime-Cam Subaru Strategic Program (HSC SSP 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 goal of this analysis is to create a homogeneous data set to study redshift evolution.Since the redMaPPer algorithm (Rykoff et al. 2014(Rykoff et al. , 2016) ) has been applied to DES data, and there exists a set of clusters observed by both SDSS and DES, this makes the process of creating 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 (Illustris TNG300-1) used to measure stellar masses, halo masses, and m gap values for our SMHM relation.In Section 3, we describe the hierarchical Bayesian 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.

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: the Sloan Digital Sky Survey Data Release 12 (SDSS DR12) (Alam et al. 2015) and the Dark Energy Survey Year 3 (DESY3) (Flaugher et al. 2015a;Abbott et al. 2018;Morganson et al. 2018), which are briefly summarized below.
DES is an optical-to-near-infrared photometric survey covering 5,000 deg 2 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 observation were taken over a 6-year period, beginning in 2013.The observations were taken at the Cerro Tololo Internation Observatory (CTIO) in Chile using the ≈3 deg 2 CCD Dark Energy Camera (DECam Flaugher et al. 2015b) on the Blanco 4-m telescope.The data used in this analysis were taken over the first three years of observations.
SDSS is an photometric survey with overlapping spectroscopic data collected by Baryon Oscillation Spectroscopic Survey (BOSS) that covers a footprint in the northern sky of 14,055 deg 2 in the SDSS griz bands (for the purpose of this analysis only the g−, r−, and i−bands were used) observed between 1998-2009.The observations were taken using the Sloan Foundation 2.5m 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 the prior studies is in the radii within which the SDSS BCG magnitudes are measured, as discussed in Section 2.2.
GM&M19 used clusters in the redshift range 0.08 ≤ z ≤ 0.12 to characterize any differences in the halo masses, richnesses, central galaxy magnitudes, stellar masses, and magnitudes gaps between the SDSS-C4 and SDSS redMaPPer data.By conducting a direct comparison on individual clusters in both data sets, they ruled out systematic differences in the mean observables (e.g., biases) between the two samples which 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 8 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.

BCG identification
The BCG in every cluster is identified using a combination of visual identification, magnitudes, 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., clusters that appear in two or more of the three catalogs), we visually confirm the BCG (a thorough discussion of miscentering in redMaPPer is provided in Hoshino et al. (2015); Zhang et al. (2019b)).We note that the BCG must lie within 0.5×R 200,crit of the cluster center (see Section below for details).The photometric algorithms and visual confirmations ensure 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-100kpc.Therefore, to homogeneously infer the slope, we use fixed physical aperture BCG magnitudes as opposed to alternatives such as the Petrosian or Kron magnitude (Petrosian 1976;Kron 1980).
For the SDSS BCGs, we use the SDSS pipeline processed radial light profiles to measure 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 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 magnitude in the i−band out to a radius of 2.5R 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 co-added 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 500kpc from the BCG.In Figure 3 we compare the r−band light profiles for a single BCG in both the 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 compared to 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 which are above the background to 100kpc.We use this subset of data to quantitatively characterize the radially dependent magnitude differences between DES and SDSS BCG photometry.Beyond 50kpc, the SDSS photometric measurements are consistently fainter than the DES photometry as shown in Figure 4.The differences begin to grow beyond 50kpc, with the DES magnitudes nearly 0.5 magnitudes brighter in the g−band compared to SDSS when measured at 100kpc.The differences are less pronounced in the r− and i−bands, but large enough to cause concern about using the 100kpc aperture magnitude for BCGs.Based on this analysis, we choose 50kpc as the BCG aperture magnitudes 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−band.
As a final test, we compare the colors of the 48 BCGs which exist in both 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 mean (error) is 0.035 magnitudes with a standard deviation of 0.137 magnitudes.Therefore, we find no bias between the two data sets for the BCG colors.

BCG stellar masses
We use the observed radial light profiles to measure the projected BCG luminosity/magnitude within a 50kpc aperture.The magnitudes in different bands allow us to estimate the stellar mass.We follow the same procedure as outlined in GM&M19, summarized here.For each cluster, redMaPPer assigns every galaxy a membership probability, P mem , which is dependent on the cluster's richness, density profile, and background density (Rykoff et al. 2014(Rykoff et al. , 2016)).The P mem > 0.7 members are then used to estimate the cluster's photometric redshift, which we use for the BCG.
Given the BCG apparent magnitudes, colors, and redshifts, we then fit a passively evolving spectral model using the SED modelling software package EzGal (Mancone & Gonzalez 2012).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 Funtion (IMF), a formation redshift of z=4.9, and a metallicity of either 0.008, 0.02, or 0.05.The choice of metallicity for each DES BCG is determined based on which model yielded the lowest chi-squared statistic between the measured and modelled 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.To determine the best fit SED, we use a Bayesian MCMC approach using emcee (Foreman-Mackey et al. 2013), where we treat the absolute magnitude (the Ez-Gal 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 α parameter.
Based on a comparison between 61 clusters between 0.20 ≤ z ≤ 0.35 in both SDSS and DES (with light pro-  2. For such a comparison, we scaled the SDSS photometry to match that of the DES photometry, since there are slight differences in the wavebands used for the analysis.The DES profile matches the SDSS within the inner 50kpc region, but then becomes brighter as more light above the background is observed in the deeper DES observations.files out to 50kpc), we find that the mean difference between our stellar mass measurements is 0.04 dex with a standard devation of 0.07.Thus, the difference between the 50kpc SED inferred BCG stellar masses from two independent imaging surveys is statistically consistent with zero.Therefore, we find excellent agreement between stellar masses estimated for BCGs with both DES and SDSS photometry.We also test to confirm that 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.-The difference between the SDSS and DES photometric measurements for 48 BCGs identified in both the SDSS and the 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 > 50kpc, particularly in the g−band where differences between the two surveys are larger than the average magnitude measurement error on the BCG magnitude.

DES Cluster Richnesses and Masses
the DES Year 1 redMaPPer data, given by Equation 2.
M halo /(h −1 M ) = 10 14.344 (λ/40) 1.356 ( 1 + z red 1 + 0.35 ) −0.30 (2) For both equations, M halo is M 200m , z red is the redMaP-Per photometric redshift and λ 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 DES redMaPPer McClintock et al. (2019) version.We note that we are actually using DES Y3 data and that a preliminary analysis of the DES-redMaPPer Y3 richness estimate is consistent when compared to Y1 analysis.We also compare the redMaPPer richnesses for the 61 clusters which are in both SDSS and DES and find excellent agreement between the two measurements, with a difference of 0.9 ± 10.1.This translates to an offset of 0.01± 0.11 dex in units of log 10 M in halo mass.

m gap
For the low redshift SDSS-C4 clusters, we use all available spectroscopic information to identify the 4 brightest cluster galaxies and measure m gap .We use a radius of 0.5× R 200 , where R 200 is the radius where the mass density reaches 200× the critical density.For the SDSS-C4 clusters, R 200 is determined from the clusters 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 m gap values.As discussed in GM&M19, we use all galaxies with a redMaPPer P mem ≥ 0.984 to define the 4th brightest galaxy within half the virial radius.This high P mem value was chosen because it yielded a match between 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 have 4 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): where λ is the redMaPPer richness, and R c is the redMaPPer cutoff radius, given by Equation 4: Using this estimate for R 200 , we define M14 as the difference in the r−band apparent model magnitude of 4th brightest cluster member with P mem ≥0.984 within 0.5R 200 and the BCG's 50kpc 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 50kpc magnitude ensures consistency between the DES and SDSS BCG magnitudes as previously shown in Figure 4, and m gap values.
We use the same sample of 61 clusters as before to compare M14 measurements between the DES and SDSS data and find ∆M14 = −0.08 ± 0.5, which is consistent with zero.We note that there is significant scatter in the data when comparing m gap measurements between DES and SDSS, which likely results from different galaxies being identified as the 4th brightest cluster member.

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 survey data.Recall that we use two surveys (SDSS and DES) and three cluster catalogs to make a combined sample which can be analyzed with and without redshift binning over the range 0.03 ≤ z ≤ 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 GM&M19, where we analyze a redshift bin which 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 ≤ z ≤ 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 by using strong priors on the four main parameters which describe the richness dependent SMHM relation, the offset α, the slope β, the stretch γ, and the intrinsic scatter σ int in the multivariate linear relation: Each BCG stellar mass is treated as a draw from a Normal distribution with mean defined by the above equation and an intrinsic scatter (standard deviation) defined as σ int .The priors used for α, β, γ and σ int were also de-fined 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 pivot = 14.65 and y pivot = 11.50).Doing this subtraction removes covariance between α and β.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 treated the uncertainties on the SDSS-redMaPPer observables as free parameters and regressed 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 which can be inferred from both the dynamics (caustic halo masses) and from a mass-richness relation which 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 which 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 as the initial rung, we ensure consistency between our measured values, in particular σ int and those of prior studies.In GM&M18, we found excellent agreement between our σ 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 follow the same procedure described above to calibrate the uncertainties for the DES BCG stellar masses and m gap values.The SDSS-redMaPPer and DES-redMaPPer samples overlap within ∼hundred 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 as SDSS-Calibration and DES-Calibration in the Appendix.
We note that we have remeasured the SDSS-redMaPPer BCG stellar masses using a 50kpc 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 stellar mass, the SDSS-redMaPPer BCG sample is larger since more BCGs have light profiles measured out to 50kpc than to 100kpc.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 at a fixed aperture in the light profiles (on average).Second, the deeper DES data should make it more difficult to identify the 4th brightest galaxy (on average), 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 mea-surements.
We conduct the simplified (no evolution) Bayesian analysis on the DES-Calibration sample and use the SDSS-Calibration posterior distributions for α, β, γ, and σ 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 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 magnitudes).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 since 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), which finds the value to be between 0.17-0.21.
We do one final test to ensure that our Bayesianinferred 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 β, γ, and σ int be free in the Bayesian regression.We note that we fix α to the value measured from the SDSS-Calibration, our remeasurement of the results from GM&M19.This is because of the strong degeneracy between α and γ.The results are presented in rows 2 and 3 of the Appendix.We also report the results from GM&M19 in row 1 of Table 3.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σ agreement between the GM&M19 measurements and our SDSS and DES-Calbration measurements for all parameters except α, which is impacted most significantly by the offset and is a 1.5σ 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 that the mean values of the observables are identical where we expect them to be (i.e., Sections 2.2, 2.3, 2.4 and 2.5) we avoid the possibility of inferring evolution where there is none.After quantifying those 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 SMHM parameters.

DES Final Sample
We analyze how the SMHM relation evolves with redshift using two approaches.First, we divide our data, in-corporating the SDSS-C4, SDSS-redMaPPer, and DES-redMaPPer data, into 8 bins sorted by redshift and measure our Bayesian MCMC posteriors for each bin with 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 over SDSS and DES clusters.We note that we include all clusters, considered in GM&M19, not just those in 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 multi-survey data set covers the redshift range 0.03 ≤ z ≤ 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 data from DES over SDSS since the photometry is deeper and we keep the 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 /h, 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 sample of data to have differing halo mass lower limits as we move to higher redshift ranges.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 a m gap completeness criterion based on the binning of the BCG and 4th 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;Golden-Marx & Miller 2018).Additionally, since the halo mass distribution can be approximated as Gaussian, the peak indicates the mass that 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 those 2323 clusters in our final sample, 1062 come from SDSS and 1261 come from DES.
Figure 5 visualizes our final sample and shows the 50kpc 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 show 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.For the simulated analysis of this study, we examine the evolution of the SMHM -m gap relation using the magneto-hydrodynamical cosmological galaxy formation simulations from the IllustrisTNG1 suite (Nelson et al. 2019;Pillepich et al. 2018b;Springel et al. 2018).Specifically, we use the Illustris TNG300-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 which best match the median of our 8 binned samples given in the Appendix.For this comparison, we identify the 260 clusters with log10(M 200m /(M /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 3D information provided directly from each snapshot of the TNG300-1 simulation, including the M 200m halo masses, measured within R 200 × ρ m ; the galaxy positions, x, y, z; R 200 ; 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 200 , and within 2σ from the red sequence, are identified as members.
For the observed data, we use a 50kpc aperture on the deblended BCG (see Section 2.2 and 2.3).In the simulation, the parent halo and sub-halo have been identified using Subfind (Springel et al. 2001).We therefore use a fixed 50kpc 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 3 dimensional positions of the satellites are known.In Figure 7, we show that all but 4 BCGs have a close companion within the 50kpc aperture.The peak of the distribution is around 3 or 4 satellites per BCG, with some -The distribution of the number of satellite galaxies within 50kpc of the TNG300-1 BCGs, which highlights that either masking or using the particle information directly is necessary to estimate stellar masses of TNG 300-1 BCGs.
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 which is 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.
We use the snapshot information for each BCG subhalo in each of the 8 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 4th brightest member and the r−band model magnitude of the TNG300-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.

BAYESIAN MCMC MODEL
We use a similar hierarchical Bayesian MCMC approach to what is described in GM&M19 to determine the values of α, β, γ, and σ int given in Equation 5. Any changes in the underlying model are designed to improve the efficiency of our analysis.The Bayesian formalism works by convolving 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 gap values at each step in our likelihood analysis, which are directly compared to our observed measurements.

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 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 (σ x0 , σ y0 , σ z0 ) as well as a stochastic component from a beta function, β(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 σ x , σ y , and σ z .

The Unobserved Quantities
Our aim is to constrain the parameters of the SMHM relation: the offset, slope, stretch, and intrinsic scatter (α, β, γ, and σ int ) as a function of redshift.In GM&M19, we modeled the evolution of these parameters as powerlaws (1 + z) ni , where n i defined the amount and shape of the redshift evolution for each of the four model parameters, α, β, γ, σ int .
In this work, we extend the data from the limit of z ∼ 0.3 to z ∼ 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 capture 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.
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 α, β, γ and σ 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.4Gyr, which is a redshift of ∼ 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 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: evolving offset, α(t) where x, y, z are the observed halo masses, BCG stellar masses, and m gap values and t is the lookback time, calculated using the photometric redshift.α 0 , β 0 , γ 0 are free parameter offsets which 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 σ int (t) that evolves with redshift: σ int0 + ((n σ ) * (0.4 nσ )/(t) 1+nσ ).n α , n β , n γ , and n σ measure the redshift evolution of α, β, γ, and σ int 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 vs. without) using only the posterior distribution.We also note that in Equation 7, the values of α 0 , β 0 , γ 0 , and σ int0 , represent the parameter values at large t, which when the shape is flat is represented by the maximum redshift of the sample.This is different than the α, β, γ, and σ int parameters in GM&M19 which represent the values at z = 0.0.Therefore, these two sets of parameters cannot be 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 gap 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.
This is a hierarchical Bayes model because the priors on the true halo masses (x i ) and M14 values (z i ) depend on models themselves (the observed halo mass and M14 distributions).

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 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, α 0 , β 0 , γ 0 , σ int0 , n α , n β , n γ , and n σ .
A negative (positive) value for the evolution parameters (the n's) indicates that the parameter itself (α(t), β(t), γ(t), σ int (t)) is growing (shrinking) with increasing lookback time.The 1D marginalized posteriors in Figure 8 and given in row 3 of Table 1 illustrate the evidence for evolution in the SMHM relation in the offset (n α ), slope (n β ), and stretch (n γ ).We find no evidence for evolution in the intrinsic scatter (n σint ), which is small, well within 2σ of 0.0, and consistent with 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 100kpc apertures for the BCG stellar masses.In addi- The high-mass power law slope Linear Regression Prior γ 0 The stretch parameter, which describes the stellar mass -M14 stratification Linear Regression Prior σ int0 The uncertainty in the intrinsic stellar mass at fixed halo mass U (0.0, 0.5) y i The underlying distribution in stellar mass Equation 7x i The underlying halo mass distribution N (14.23,0.18 2 ) z i The underlying mgap distribution N (2.51,0.62 2 ) nα The shape parameter associated with the evolution of α U (−5.0, 5.0) n β The shape parameter associated with the evolution of β U (−5.0, 5.0) nγ The shape parameter associated with the evolution of γ U (−5.0, 5.0) nσ The power law associated with the evolution of σ int U (−20.0, 20.0) σy 0i The uncertainty between the observed stellar mass and intrinsic stellar mass distribution 0.08 or 0.06 dex σx 0i The uncertainty associated with the mass-richness relation 0.087 dex σz 0i The uncertainty between the underlying and observed mgap distribution 0.15 or 0.31 tion, 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 50kpc aperture stellar masses, and higher redshift data out to z = 0.6.Therefore, we re-analyze the GM&M19 data using this new model (Equation 7) on that original data set, as well as with the new model and using 50kpc 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σ statistical level.Switching from their power law fitting function to the Pareto function, we find that the evolution in the slope is also significantly non-zero, n β = 0.573 +0.152  −0.141 or ∼ 3.8σ.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σ.Most of the error bars on the parameters have also decreased compared to the analysis on the SDSS data alone.While we have nearly doubled the sample size from GM&M19, the evolution in β(t) is in the late universe where DES does not provide new data.However, the DES data is useful to pin 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 α(t), β(t), γ(t) and σ int (t) are evolving between z = 0.3 and z = 0.6 , we would expect to see differences in the zero points α 0 , β 0 , γ 0 and σ int0 between the second and third rows because these parameters represent the value 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 β 0 or σ int0 after we extend the analysis to z = 0.6 using the DES data.However, we do find that α 0 is significantly higher and γ 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 α at 5.8σ and evolution in the stretch γ at 4.0σ.
We also consider what would happen if we excluded m gap from our model by dropping z i in Equation 7. We find that the significance of evolution of the slope drops from n β = 0.263 +0.086  −0.075 to being statistically consistent with zero (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 gap in the analysis.
It is interesting that the offset α still shows statistically significant evolution when we ignore m gap in the analysis.We note that in equation 7, the offset parameter (α 0 ) is not a direct measure of the amplitude of the SMHM.Even when incorporating n α , the first term α(t) does not quantify the amplitude of the SMHM (because of the inclusion of m gap ).However, without M14, α(t) is simply the overall amplitude of the SMHM relation as a function of lookback time, which is characterized by α 0 (at the redshift limit of the data) and n α .Thus, when we ignore m 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 α(t) traces the evolution in the amplitude, the sign on the evolution n α 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 α(t) in the next subsection.7. The posterior distribution for α 0 , β 0 , γ 0 , σ int0 , nα, n β , nγ , and nσ.As in GM&M18, we see that γ is significantly non-zero and σ 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 α 0 , β 0 , γ 0 , and σ 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α is 5.8σ from 0.0, n β is 3.5σ from 0.0, and nγ is 4.0σ from 0.0.
ter 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 α(t = t med ), β(t = t med ), γ(t = t med ) and σ int (t = t med ) on data separated into discrete redshift bins shown in Figure 1 and where t = t 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.
In Figure 9, we compare the SMHM parameter values from Equation 7with the evolution parameters fixed to zero (purple dots) against the fully evolving parameters for the offset α(t), slope β(t), stretch γ(t), and intrinsic Note.-Equation 7parameter fits for the SDSS data in GM&M19 reference (z < 0.3) compared to the fits in his paper which use DES data to extend the analysis to z = 0.6.a The same data from GM&M19 was re-analyzed using the model from this paper (e.g., equation 8) and the original 100kpc apertures.b The same data from GM&M19 was re-analyzed using the model from this paper (e.g., equation 8) and 50kpc 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σ agreement as shown in Figure 9.
scatter σ int (t) (green line and green band).The purple error bars are the 1σ error bars on the binned posteriors.The green error band incorporates the error on the parameter zero points (e.g., α 0 , β 0 , γ 0 , and σ int0 ) as well as the error on the corresponding evolution component (n α , n β , n γ and n σint ).
We find good agreement between the two separate analyses: binned and unbinned.We first note the evolution in the slope parameter β(t) of the SMHM, which is quantified via the Pareto function as n β = 0.263 +0.086 −0.075 and is evident in Figure 9.The Pareto function does a good job 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 3Gyrs, 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 σ int0 = 0.085 ± 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.3Gyr exists in the binned analysis.We are unable to explain this feature of the data, which is an 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.Like for the slope, the Pareto function captures the shape of the observed binned evolution, which shows a gradual increase (for α) or decrease (for γ).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, α, we are not actually detecting evolution in the overall amplitude (or median stellar mass at fixed halo mass and m gap ) for the SMHM relation.As we show in Figure 10, we plot the combination of the offset 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.

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

Parameter Correlations
Figure 8 shows some interesting structure 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 re-centering 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, however 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 in these two parameters.We show this degeneracy just for clusters in the lowest-redshift binned analysis in Figure 11.We overplot the 2D posterior distribution between α(t) and γ(t) for the low-redshift bin (in green) and the total posterior (in purple).Figure 11 highlights that we see that 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 2 dimensional posteriors that may be responsible for part of this agreement.
Figure 11 exemplifies why one would may want to avoid binning in this type of SMHM analysis, since unaccounted parameter covariances can lead to incorrect fits to the evolving SMHM relation.Our fully unbinned analysis and our hierarchical Bayesian formalism allow  -The combined offset of α(t) + γ(t)×M14.This is for when 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 (α(t)) and the stretch term (γ(t)).-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.
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 4 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 do note that as evident from Figure 11, a higher stretch parameter at low redshift is still preferred.

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 away (or to some latent parameter), we would need to quantify those correlations and their impact 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 that 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 Xray 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 (Pearson correlation coefficient ∼ 0.2).Furnell et al. (2018) used dynamical masses to find a similar weakly positive correlation between richness and BCG stellar mass (Spearman rank correlation coefficient r s = 0.137).On the other hand, Farahi et al. (2020) used the Illustris TNG simulations to find a moderate anti-correlation between richness and BCG stellar mass at fixed halo mass (Pearson correlation coefficient ∼ −0.4).None of these reported correlations are strong or consistent.Compounding the issue is 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 /h 2 )) < 11.6, the stellar mass range that allows us to measure the correlation across the entire redshift range, we find a moderate anti-correlation of ∼ −0.4.
The statistical significance of correlation coefficients is not well defined.Most in the literature use some form of jacknife sampling (e.g., Erfanianfar et al. 2019).However, here, we use a Bayesian-like approach where we apply Equation 5 to forward model our data using 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 non-correlated 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 anti-correlation between M14 and M 200 in the Illustris-TNG sample (−0.36).
We can use this same forward modeling technique to quantify the effect 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 which 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 200 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 effect 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 Illustris TNG300-1 Figure 9 offers direct comparisons between the observed (SDSS and DES) results and the simulated TNG300-1 measurements.Such a comparison allows us to understand whether the physical prescriptions built into the TNG300-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, α is related to the median stellar mass (at a halo mass of log 10 (M halo /(M /h)) = 14.65) at a given m gap (since the γ values agree).We note that for the simulated data, like 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 σ int ; we detect no evolution in either.Interestingly, Pillepich et al. (2018a) detect modest evolution in σ 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 50kpc magnitudes for the TNG300-1 data, is the absence of noticeable evolution in the slope of the SMHM relation for TNG300-1.In our observed data set, late time growth appears to occur primarily in the last 2 billion years; however, in the TNG300-1 simulation, there is no detectable evolution over the entire time range studied.However, we note that the absence of redshift evolution in the slope with the TNG300-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 TNG300-1 simulation this is not the case.
Another difference between our TNG300-1 and prior measurements from Pillepich et al. (2018a) and Engler et al. (2021) is the value of β.We measure a value of approximately 0.42 for the slope when the stellar mass is measured within 50kpc 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 30kpc and find a slope of 0.49 (no error bars reported).We note 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 the 2 times the stellar half mass radius, a radius far greater than the 50kpc aperature we use.Therefore, we can conclude, as shown in GM&M19 that had we used a large aperture (>= 100kpc) to measure the BCG stellar mass and magnitude, then 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, which agrees with the general conclusions from Tinker et al. (2019).However, as shown in Table 3, this trend is not shown in the TNG300-1 data.Instead, we see that the slopes are within 1σ, which may serve as the first bit of evidence that the BCGs and growth prescriptions in the TNG300-1 simulation are over-dominant.
The remaining two parameters α and γ are also dramatically different.Unlike in our observational data, there is little evidence of any evolution in α or γ out to high redshift.Given that α and γ are covariant, it is unsurprising that if one of these two parameters shows no evolution, the other parameter also shows not evolution, and as discussed in Section 5 likely related to the growth prescription used in TNG.Additionally, the values for α also significantly differ.At first glance, it appears as though the TNG300-1 BCGs are undermassive.However, that is not the case.A more valid comparison would be the value of α + γ × z med , which would be representative of the median stellar mass of the BCGs at a given halo mass.This comparison yields that the TNG systems are approximately 0.06 dex overmassive.Of note, when doing such a comparison, the median M14 values for TNG are approximately 1-1.5 magnitudes greater than the observed values; unlike in the observed universe, low m gap systems (M14<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 TNG300-1 also results in poorly populated red sequences, such that few intermediate brightness galaxies exist, thus yielding substantially fainter 4th brightest galaxies.

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, 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 grow and evolve over the last 6 billion years.
Currently, there is not a clear consensus between observations, simulations, and models about how BCGs grow over this redshift range.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 key results from GM&M18 and GM&M19: m gap is definitively a latent parameter within the SMHM relation; incorporating γ and M14 into the SMHM relation reduces σ int ; and accounting for m gap yields significant evolution in the slope of the SMHM relation over late time.From this analysis, we for the first time, report evolution in both the α and γ parameters, which represent the offset and stretch, respectively.It is this observed evolution that drives our understanding of how BCG's 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 (2nd and 3rd bin) and a high redshift sample (7th and 8th bin).We note that due to 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 in (high-z) and unfilled (low-z) data points.We then overlay the results of the posterior distributions shown in Figure 8 and Figure 9 as the shaded regions.Figure 12 highlights a few of our key findings.First, γ is significantly growing as one moves forward in lookback time, as evidenced by how much larger the separation between the two shaded regions are at low-z when compared to high-z.We note that if γ were not evolving, the separation between the two high and low M14 bins would not be growing, regardless of the change in halo mass distribution of the data, since γ does not vary with M halo .Second, β is growing as one moves forward in lookback time.Third, σ int , the spread in the data at fixed m gap is unchanged between these two distributions, which supports our measurement that σ 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 /h)) < 14.7), the BCG stellar mass distribution remains the same, and thus the BCG stellar mass within 50kpc is not growing.This is also supported by the constant value of α + γ × M 14 given in Figure 10.
In GM&M19, given the absence of evolution in γ, we assumed that any growth observed was due entirely to growth in the BCG.However, as shown in Figure 12 the stellar mass within 50kpc is not growing over this redshift range.This observation highlights that the driver behind all 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, the only way for the slope to increase would be for the stellar masses of the  -We display two sets of distributions for the the low (unfilled) and high-z (filled) data.For both, we show two mgap 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 models.As shown, for the low redshift data, we see a steeper slope and more pronounced stratification, which results from a larger stretch.
distribution 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 due to their residing in richer clusters.
Recall that m gap is the difference in brightness between the BCG and 4th 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 most insight into what is physically happening can be instead gleaned from the evolution in γ.For γ 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).Therefore, instead, what is likely happening is that mergers between the bright satellite galaxies and the BCG deposit the stellar material at radii beyond 50kpc, what we interpret as the ICL.Therefore, the outer envelopes contain all 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 50kpc) remains fixed.Therefore, the incorporation of m gap has elucidated that BCGs continue to grow hierarchically in this redshift range, but all of that added stellar material is going 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 TNG300-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 TNG300-1 simulation.Additionally, the TNG300-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 over efficient merger process, there exists an absence of fainter satellite galaxies in the TNG300-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 ∼ 0.6.As shown here, β shows significant late-time evolution, predominately over the redshift range 0.0 < z < 0.15 and we for the first time detect statistically significant evolution in α and γ, which has clarified that this evolution is driven by BCG hierarchical growth that is evident not in the stellar mass, but rather within the m gap .We are left with a few paths forward.If we choose to tighten the constraints further 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 ∼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-SPT overlap (Bleem et al. 2015(Bleem et al. , 2020;;Huang et al. 2020), an approach which faces similar observational and modelling 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 massm gap trend exists in both observations and state-of-theart 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 inform us about the formation history of the BCG and its host cluster dark matter halo, and quantify how the stellar mass, halo mass, m gap parameter space maps 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 100kpc and beyond.Another vital step forward as part of that effort is 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 conclusion that all the recent growth is contained within the ICL and that it's these recent mergers, which change the m gap distribution and yield our detected evolution.This paper has gone through internal review by the DES collaboration.JGM 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.JGM acknowledges the support by 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 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.
Fig. 1.-The redshift distribution of the combined sample of SDSS and DES data for this analysis.The dashed lines represent the edges of the 8 bins used in this analysis.

Fig. 2 .
Fig.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 inset in both panels is a representation of the BCG after masking.The black circle is the 50kpc 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.

For
Fig. 4.-The difference between the SDSS and DES photometric measurements for 48 BCGs identified in both the SDSS and the 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 > 50kpc, particularly in the g−band where differences between the two surveys are larger than the average magnitude measurement error on the BCG magnitude.
Fig.5.-SMHM-M14Relation for the SDSS-C4, SDSS-redMaPPer, and DES-redMaPPer Samples colored via M14.We see that a stellar mass -mgap 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.

Fig. 6 .
Fig.6.-Anexample of a TNG300-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 500kpc centered on the BCG.The insert shows the stellar particle information for just the inner 100kpc of the subhalo containing the BCG.The circle represents the 50kpc aperture.

4. 1 .
Fig. 8.-SMHM-M14 parameter posteriors from Equation7.The posterior distribution for α 0 , β 0 , γ 0 , σ int0 , nα, n β , nγ , and nσ.As in GM&M18, we see that γ is significantly non-zero and σ int0 is less than 0.1 dex.However, we note that as a result of the modified redshift evolution form given by Equation7, the values are not directly comparable to the results from GM&M19.Instead, the posteriors for α 0 , β 0 , γ 0 , and σ int0 represent the values that these parameters asymptote to.To see the values at the redshifts measured in our study, see Figure9.Using this model, nα is 5.8σ from 0.0, n β is 3.5σ from 0.0, and nγ is 4.0σ from 0.0.
Fig. 9.-The effective offset, slope, stretch and intrinsic scatter from Equation 7 as a function of lookback time.The green line 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α) 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 TNG300-1 data in blue.
Fig.10.-Thecombined offset of α(t) + γ(t)×M14.This is for when 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 (α(t)) and the stretch term (γ(t)).
Fig. 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 Figure9for the median lookback time of bin 1.

Fig
Fig.12.-We display two sets of distributions for the the low (unfilled) and high-z (filled) data.For both, we show two mgap 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 models.As shown, for the low redshift data, we see a steeper slope and more pronounced stratification, which results from a larger stretch.

TABLE 1 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 × log(1 + 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.

TABLE 2 Posterior
Distribution Results with Lookback Time Evolution The DES participants from Spanish institutions are partially supported by MICINN under grants ESP2017-89838, PGC2018-094773, PGC2018-102021, SEV-2016-0588, SEV-2016-0597, and MDM-2015-0509, some of which include ERDF funds from the European Union.IFAE is partially funded by the CERCA program of the Generalitat de Catalunya.Research leading to these results has received funding from the European Research Council under the European Union's Seventh Framework Program (FP7/2007-2013) including ERC grant agreements 240672, 291329, and 306478.We acknowledge support from the Brazilian Instituto Nacional de Ciênciae 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.