The splashback radius of optically selected clusters with Subaru HSC Second Public Data Release

Abstract

1 Introduction

Galaxy clusters are the most massive gravitationally bound structure in the Universe, which form in dark matter halos around rare high peaks in the initial density field after interactions between gravitational dynamics and baryoni process related to galaxy formation (see Allen et al. 2011; Kravtsov & Borgani 2012; Weinberg et al. 2013; Wechsler & Tinker 2018; Pratt et al. 2019; Walker et al. 2019; Vogelsberger et al. 2020, for recent reviews). The splashback radius has been proposed as a physical boundary of dark matter halos at the outskirts which separates orbiting from accreting materials as a sharp density edge, which is seen even after stacking of halos through numerical simulations (Diemer & Kravtsov 2014; More et al. 2015) and can be explained with a simple analytical model (Adhikari et al. 2014). At the splashback radius, the logarithmic derivative of density profiles is predicted to drop significantly over a narrow range of radius in the outskirts due to the piling up of materials with small radial velocities at their first orbital apocenter after infall on to the halo. The splashback radius primarily depends on the accretion rate, redshift, and halo mass (e.g., Diemer et al. 2017). The secondary infall models with the spherical collapse model (Gunn & Gott 1972) predict such a sharp density edge analytically (e.g., Fillmore & Goldreich 1984; Bertschinger 1985; Adhikari et al. 2014; Shi 2016). Furthermore, the splashback features have been investigated from various aspects with simulations (Oguri & Hamana 2011; Mansfield et al. 2017; Diemer 2017; Okumura et al. 2018; Sugiura et al. 2020; Xhakaj et al. 2019), including dark energy (Adhikari et al. 2018), modified gravity (Adhikari et al. 2018; Contigiani et al. 2019b), and self-interacting dark matter (Banerjee et al. 2020).

Observationally, the splashback feature has been detected and constrained statistically for galaxy clusters selected by the optical cluster-finding algorithm redMaPPer (Rykoff et al. 2012, 2014, 2016; Rozo & Rykoff 2014; Rozo et al. 2015a, 2015b) through cross-correlation measurements between clusters and galaxies (More et al. 2016; Baxter et al. 2017; Chang et al. 2018; Shin et al. 2019) or weak lensing measurements (Chang et al. 2018) with high precision thanks to a large number of optically selected clusters. Less-precise measurements have also been done for clusters selected by the Sunyaev–Zel’dovich effect (Zürcher & More 2019; Shin et al. 2019) or X-ray flux (Umetsu & Diemer 2017; Contigiani et al. 2019a). We note that the signal-to-noise ratios for weak lensing measurements are smaller than those for cluster–galaxy cross-correlation measurements to constrain the splashback feature. It is important to validate theoretical predictions against observed splashback features with the help of weak-lensing mass calibrations, where weak-lensing measurements constrain the mass-observable relation mainly from their amplitudes, instead of the location of the splashback radius imprinted in lensing measurements. Interestingly, the high-precision measurements of the splashback radius for the redMaPPer clusters (More et al. 2016; Baxter et al. 2017; Chang et al. 2018) from the Sloan Digital Sky Survey (SDSS) or the Dark Energy Survey (DES) data show that their observed splashback radius from cluster–galaxy cross-correlation measurements is smaller than expected from N-body simulations under the standard Λ cold dark matter (ΛCDM) cosmology at the level of |$20\% \pm 5\%$|⁠. Shin et al. (2019) also show a similar trend at the ∼2σ level for the redMaPPer clusters at more massive mass scales.

One of the possible origins of the inconsistency is some systematics in the optical cluster-finding algorithm, such as projection effects in optical cluster-finding algorithms for which non-member galaxies along the line-of-sight as member galaxies in optical richness estimation (e.g., Cohn et al. 2007; Zu et al. 2017; Busch & White 2017; Costanzi et al. 2019; Sunayama & More 2019). In particular, Busch and White (2017) and Sunayama and More (2019) investigated projection effects on the splashback features with a simplified redMaPPer-like optical cluster-finding algorithm in the Millennium Simulation (Springel et al. 2005). However, their cluster abundance densities in the simulations is ∼3 times larger than observations as a function of richness, mainly due to the use of a mock galaxy population from a semi-analytical model, which could lead to the overestimation of the projection effects. Also, Busch and White (2017) and Sunayama and More (2019) did not analyze their mock data by fitting the model profile proposed in Diemer and Kravtsov (2014) to the observed projected cross-correlation function, thus not entirely mimicking the procedure adopted in observational studies (e.g., More et al. 2016; Chang et al. 2018). Therefore, we can improve this aspect in mock simulation analyses to investigate possible bias effects from optical cluster-finding algorithms by more closely resembling the comparison procedure in the observations. Thus, further independent investigation would be useful to confirm whether the inconsistency is caused by artifacts due to optical cluster-finding algorithms or not.

The distribution of galaxies after magnitude or color cuts in observations is expected to have the same splashback radius as dark matter distribution to the first-order approximation, unless such cuts select particular orbits of recently accreted galaxies in clusters due to fast quenching effects, for instance. The difference between dark matter and subhalos, where galaxies should reside, can be caused by dynamical friction (Chandrasekhar 1943) that affects subhalos with a ratio of mass to host halo mass M_sub/M_host ≳ 0.01 (Adhikari et al. 2016) to decrease splashback radii significantly (⁠|${\gtrsim}5\%$|⁠) compared to the dark matter. In observations, fainter galaxies are expected to live in less-massive subhalos, and thus dynamical friction effects are smaller. Hence, fainter galaxies are expected to have splashback radii closer to those for the dark matter. However, relations between subhalo masses and galaxy luminosities (including scatter between them) have not been fully understood. Dynamical friction effects on splashback radii are investigated in the literature (More et al. 2016; Chang et al. 2018; Zürcher & More 2019) with different absolute magnitude thresholds. However, a magnitude dependence of splashback radii for massive clusters has not been detected even with smaller error bars in More et al. (2016) at a |${\lesssim}5\%$| precision level, even though a |${\gtrsim}5\%$| shift of splashback radii is expected in simulations especially for more massive subhalos when selecting subhalos with maximum circular velocities throughout their entire history to reproduce galaxy abundances in observations (More et al. 2016; Chang et al. 2018). The scatter between the maximum circular velocities (or subhalo masses) and galaxy luminosities could wash out dynamical friction effects (More et al. 2016). It would be interesting to explore dynamical friction effects with an independent optical cluster-finding algorithm and fainter galaxy samples.

One of the most striking features within galaxy clusters is the presence of a larger number of red and elliptical galaxies with little ongoing star formation (Oemler 1974; Dressler 1980; Dressler & Gunn 1983; Balogh et al. 1997; Poggianti et al. 1999), exhibiting a tight relation in the color–magnitude diagram (e.g., Stanford et al. 1998). These quenching effects are expected from intra-cluster astrophysical processes including tidal disruption, harassment, strangulation (Larson et al. 1980), and ram-pressure stripping (Gunn & Gott 1972), or from the age-matching model (e.g., Hearin et al. 2014) which employs expectations that galaxies in larger overdensities form earlier. Quenching of satellite galaxies at the outskirts of galaxy clusters has been investigated (e.g., Wetzel et al. 2013; Fang et al. 2016; Zinger et al. 2018; Adhikari et al. 2019). The splashback radius separates orbiting from accreting galaxies as the physical boundary of halos, and thus investigating relations between the splashback and quenching features around clusters over a wide range of redshift would be informative to constrain galaxy formation and evolution physics.

The Hyper Suprime-Cam Subaru Strategic Program (HSC-SSP) is a wide-field optical imaging survey (Miyazaki et al. 2012, 2015, 2018; Furusawa et al. 2018; Kawanomoto et al. 2018; Komiyama et al. 2018) with a 1.77 deg² field-of-view camera on the 8.2 m Subaru telescope. The HSC survey is unique in its combination of depth and high-resolution image quality, which allows us to detect galaxy clusters over a wide redshift range up to z_cl ∼1 and to investigate their splashback features with photometric galaxies down to fainter apparent magnitudes compared to other ongoing wide-field surveys, although the survey area is smaller than those of other imaging surveys in the literature (More et al. 2016; Chang et al. 2018).

In this paper, we present constraints on the splashback features of the HSC optically-selected clusters in the cluster redshift range of 0.1 ≤ z_cl ≤ 1.0 and the richness range of 15 ≤ N ≤ 200 detected by the independent cluster-finding algorithm CAMIRA (Oguri 2014) from the HSC-SSP second public data release (Aihara et al. 2019). The main difference between the redMaPPer and CAMIRA cluster-finding algorithms is that CAMIRA subtracts the background galaxy levels to estimate the richness in cluster selection locally, but redMaPPer does this globally. Moreover, CAMIRA employs a larger aperture size (≃1 h⁻¹ Mpc in physical coordinates) independent of richness values in cluster selections, whereas redMaPPer uses a smaller richness-dependent aperture (<1 h⁻¹ Mpc in physical coordinates for almost all clusters). To constrain the splashback features, we measure the projected cross-correlation function between the HSC CAMIRA clusters and the HSC photometric galaxies. We investigate the dependence of splashback features for the CAMIRA clusters on cluster redshift, richness, galaxy magnitude limit, and galaxy color. We compare results from a model fitting to the projected cross-correlation measurements with predictions from the halo–matter cross-correlation function in Dark Emulator (Nishimichi et al. 2019), which is constructed from a suite of high-resolution N-body simulations, and the mass–richness relation of the HSC CAMIRA clusters (Murata et al. 2019) constrained by stacked weak gravitational lensing and abundance measurements with a pipeline developed in Murata et al. (2018). We also investigate how projection effects could bias splashback feature measurements using a simplified estimation of the impact of the optical cluster-finding algorithms partly following approaches presented in Busch and White (2017) and Sunayama and More (2019). For this purpose, we construct mock galaxy catalogs by populating galaxies in N-body simulations with a halo occupation distribution (HOD; Zheng et al. 2005) model to match cluster abundance densities and lensing profiles to observations approximately.

The structure of this paper is as follows. In section 2, we describe the Subaru HSC data and catalogs used in our splashback feature analysis and the projected cross-correlation function measurements between the HSC clusters and galaxies. In section 3, we summarize a model fitting method with a parametric density profile from Diemer and Kravtsov (2014) to constrain splashback features from the projected cross-correlation measurements. In section 4, we show the resulting constraints on the splashback features and compare them with model predictions. We discuss the robustness of our results in section 5. Finally, we conclude in section 6.

Throughout this paper, we use natural units where the speed of light is set equal to unity, c = 1. We use |$M \equiv M_{\rm 200m}=4 \pi (R_{\rm 200m})^3 \bar{\rho }_{\rm m0} \times 200/3$| for the halo mass definition, where R_200m is the three-dimensional comoving radius within which the mean mass density is 200 times the present-day mean mass density. In this paper, we use a radius and density in comoving coordinates rather than in physical coordinates, unless otherwise stated. We adopt the standard flat ΛCDM model as the fiducial cosmological model with cosmological parameters from the Planck15 result (Planck Collaboration 2016): Ω_b0 h² = 0.02225 and Ω_c0 h² = 0.1198 for the density parameters of baryon and cold dark matter, respectively, Ω_Λ = 0.6844 for the cosmological constant, σ₈ = 0.831 for the normalization of the matter fluctuation, and n_s = 0.9645 for the spectral index. We employ the convention used in More et al. (2015, 2016) for the splashback radius, where |$r_{\rm sp}^{\rm 3D}$| and |$R_{\rm sp}^{\rm 2D}$| correspond to the locations of the steepest logarithmic derivative of the radially averaged density profiles after stacking of halos in the three-dimensional or projected space, respectively, due to its accessibility in observations.

2 Data and measurement

2.1 HSC-SSP survey

HSC is a wide-field prime focus camera with a |${1{_{.}^{\circ}}5}$| diameter field-of-view mounted on the 8.2 m Subaru telescope (Miyazaki et al. 2012, 2015, 2018; Furusawa et al. 2018; Kawanomoto et al. 2018; Komiyama et al. 2018). With its unique combination of a wide field-of-view, a large aperture of the primary mirror, and excellent image quality, HSC can detect relatively higher redshift clusters and fainter galaxies compared with other ongoing wide-field surveys. To fully exploit its survey power, the HSC-SSP (Aihara et al. 2018a, 2018b, 2019) conducted a multi-band wide-field imaging survey over six years with 300 nights. The HSC-SSP survey consists of three layers: Wide, Deep, and UltraDeep. The Wide layer aims at covering 1400 deg² of the sky with five broad-bands (grizy) and a 5σ point-source depth of r ∼ 26. Bosch et al. (2018) describes the software pipelines used to reduce the data.

In this paper, we use HSC cluster and galaxy catalogs based on the photometric data in the Wide layer with the five broad-bands in the HSC-SSP second Public Data Release (PDR2; Aihara et al. 2019). PDR2 includes the datasets from 2014 March through 2018 January, about 174 nights in total. The median seeing for i-band images is |$\sim {0_{.}^{\prime \prime }6}$| in terms of the point-spread function (PSF) full width at half-maximum. The typical seeing, 5σ depth for extended sources, and saturation magnitude for z-band images in the Wide layer are |${0_{.}^{\prime \prime }68}$|⁠, 25.0 mag, and 17.5 mag, respectively, after averaging over the entire survey area included in PDR2 (Aihara et al. 2019).

2.2 HSC CAMIRA cluster catalog

We use the CAMIRA (Cluster-finding Algorithm based on Multi-band Identification of Red-sequence gAlaxies) cluster catalog from the HSC-SSP Wide dataset in PDR2, which was constructed using the CAMIRA algorithm (Oguri 2014). This catalog is an updated version of the one presented in Oguri et al. (2018) based on the first-year HSC-SSP dataset. The CAMIRA algorithm is a red-sequence cluster finder based on a stellar population synthesis model (Bruzual & Charlot 2003) to predict colors of red-sequence galaxies at a given redshift and to compute likelihoods of being red-sequence galaxies as a function of redshift. The stellar population synthesis model in the CAMIRA algorithm is calibrated with spectroscopic galaxies to improve its accuracy (Oguri et al. 2018). The richness in the CAMIRA algorithm corresponds to the number of red cluster member galaxies with stellar mass M_* ≳ 10^10.2 M_⨀ (roughly corresponding to the luminosity range of L ≳ 0.2 L_*) within a circular aperture with a radius R ≲ 1 h⁻¹Mpc in physical coordinates. Unlike the redMaPPer algorithm (Rykoff et al. 2012), the CAMIRA algorithm does not employ a richness-dependent aperture radius to define the richness. The HSC images are deep enough to detect cluster member galaxies down to M_* ∼ 10^10.2 M_⨀ even at redshift z_cl ∼ 1, which allows us to detect such high-redshift clusters reliably without the richness incompleteness correction. The algorithm employs a spatially-compensated filter that subtracts the background level locally in deriving the richness. CAMIRA defines cluster centers as locations of the identified brightest cluster galaxies (BCGs). The bias and scatter in photometric cluster redshifts of Δz_cl/(1 + z_cl) are shown to be better than 0.005 and 0.01, respectively, with 4σ clipping over most of the redshift range by using available spectroscopic redshifts of BCGs (Oguri et al. 2018). We use the cluster catalog after applying bright-star masks for grizy-bands presented in Aihara et al. (2019).

The catalog contains 3316 clusters at 0.1 ≤ z_cl ≤ 1.0 with the richness 15 ≤ N ≤ 200 with almost uniform completeness and purity over the sky region. We divide the whole sample into subsamples with different redshift and richness selections to investigate redshift and richness dependences of splashback features. The sample selections and characteristics are shown in table 1. The total area for this version of the HSC CAMIRA clusters is estimated to be 427 deg². In addition to the cluster catalog, we use a random cluster catalog with the same redshift and richness distributions as the data in the same footprint. There are 100 times as many random points as real clusters to measure projected cross-correlation functions with galaxies accurately.

We also employ the HSC CAMIRA member galaxy catalog to define color cuts to separate red and blue galaxies at each cluster redshift, following a red-sequence-based method in Nishizawa et al. (2018), which we will use to investigate galaxy color dependence of splashback features in subsection 4.3. These member galaxies are red-sequence galaxies identified by the CAMIRA algorithm as described in Oguri (2014). We use all the red-sequence member galaxies for clusters with richness N ≥ 10 to define separation criteria for red and blue galaxies (see Nishizawa et al. 2018). The total number of member galaxies is 331168 at 0.1 ≤ z_cl ≤ 1.0.

We use the richness–mass relation presented in Murata et al. (2019) to model splashback features for each redshift and richness selection, as described in appendix  1. Murata et al. (2019) constrained the richness–mass relation P(N|M, z) from the stacked weak gravitational lensing profiles from the HSC-SSP first-year shear catalog (Mandelbaum et al. 2018a, 2018b) and the photometric redshift catalog (Tanaka et al. 2018), and abundance measurements for the HSC CAMIRA clusters with 0.1 ≤ z_cl ≤ 1.0 and 15 ≤ N ≤ 200. We note that Murata et al. (2019) used the previous version of the HSC CAMIRA cluster catalog from the first-year dataset presented in Oguri et al. (2018). The richness differences from the updated version for matched clusters in both the versions are negligibly small compared to the scatter values in the mass–richness relation. Murata et al. (2019) also constrained the off-centering between real cluster centers and the identified BCGs for the HSC CAMIRA clusters from the stacked lensing profiles. In section 3, we use these results as inputs to our prior distributions of off-centering parameters to marginalize over its effect on the projected cross-correlation measurements between the HSC CAMIRA clusters and photometric galaxies.

2.3 HSC photometric galaxy catalog

We use the photometric galaxy catalog from PDR2 in Aihara et al. (2019) to measure projected cross-correlation functions with the clusters. We employ |$\tt {cmodel}$| magnitudes (see Bosch et al. 2018, for more details) for total z-band magnitudes in absolute-magnitude cuts (Nishizawa et al. 2018) because the z-band magnitude is more suitable for measurements of brightness of high-redshift galaxies at z_cl ∼ 1.0 than the i-band magnitude used in the literature for lower redshifts (e.g., More et al. 2016), given the wavelength of 4000 Å break at such high redshifts . We define magnitudes in the other broad-bands such that we derive colors of individual galaxies with afterburner aperture photometries (Bosch et al. 2018; Oguri et al. 2018; Nishizawa et al. 2018) to avoid inaccurate photometries in crowded regions, including cluster centers, due to a failure of the deblender in such regions. The HSC pipeline measures the afterburner aperture photometries after the pipeline matches the PSF sizes on the undeblended images between all five broad-bands to have accurate colors even in crowded regions. We use the target PSF size of |${1_{.}^{\prime \prime }3}$| and the aperture size of |${1_{.}^{\prime \prime }5}$| in diameter (⁠|$\tt {undeblended\_convolvedflux\_3\_15}$|⁠) for PDR2 (Aihara et al. 2019).

We apply the following quality cuts to select galaxies for the measurements. First, we select galaxies observed in all five broad-bands by imposing cuts in the number of exposure images used to create a co-add image for each galaxy, as described below. While the nominal definition of the full depth are four exposures for gr bands and six for izy bands in the Wide layer of PDR2, we adopt the following relaxed conditions to avoid gaps in the galaxy distribution due to CCD gaps.

• |$\tt {[gr]\_inputcount\_value}\ge 2$|

• |$\tt {[izy]\_inputcount\_value}\ge 4$|

Secondly, we apply the following basic flag cuts for all five broad-bands to remove galaxies that can be affected by bad pixels or poor photometric measurements. We also focus only on unique detections(see Aihara et al. 2018b, 2019).

• |$\tt {[grizy]\_pixelflags\_edge}=\tt {False}$|

• |$\tt {[grizy]\_pixelflags\_interpolatedcenter}=\tt {False}$|

• |$\tt {[grizy]\_pixelflags\_crcenter}=\tt {False}$|

• |$\tt {[grizy]\_cmodel\_flag}=\tt {False}$|

• |$\tt {[grizy]\_undeblended\_convolvedflux\_3\_15\_flag}=\tt {False}$|

• |$\tt {isprimary}=\tt {True}$|

Thirdly, we use galaxies with z-band cmodel magnitude brighter than 24.5 mag after dust extinction corrections (⁠|$\tt {a\_z}$|⁠), flux signal-to-noise ratio larger than 5, and with an i-band star-galaxy separation parameter indicating extended sources. We use the i-band star–galaxy separation because i-band images are taken in better seeing conditions on average. We also place weak constraints on the ri-band magnitudes as r < 28 mag and i < 28 mag to remove artifacts.

• |$\tt {z\_cmodel\_mag - a\_z}<24.5$|

• |$\tt {z\_cmodel\_flux} > 5\times \tt {z\_cmodel\_fluxsigma}$|

• |$\tt {i\_extendedness\_value}=1$|

We then apply bright-star mask cuts for all five broad-bands.

• |$\tt {[grizy]\_mask\_s18a\_bright\_objectcenter}$|=False

We use angular positions (⁠|$\tt {ra}$| and |$\tt {dec}$|⁠) of the galaxies for the measurements. The catalog after these quality cuts includes 48882039 galaxies for analyses in subsections 4.1 and 4.2. Also, we apply absolute-magnitude cuts for additional galaxy selections in each cluster redshift bin, as described in subsection 2.4.

We note that these quality cuts for the galaxy catalog are the same as those used for constructing the cluster catalog from PDR2 in subsection 2.2 except that |$\tt {z\_cmodel\_mag - a\_z}<24.0$| and |$\tt {z\_cmodel\_magsigma}<0.1$| are used for the galaxies instead of |$\tt {z\_cmodel\_mag - a\_z}<24.5$| and |$\tt {z\_cmodel\_flux} >5\times \tt {z\_cmodel\_fluxsigma}$|⁠. In the former case, the catalog includes 31719028 galaxies. These cuts are almost the same as those used for constructing the cluster catalog from the HSC first-year dataset presented in Oguri et al. (2018). Nishizawa et al. (2018) also used similar quality cuts for the HSC first-year datasets. We employ the catalog with these cuts for analyses with red and blue galaxies in subsection 4.3 because we define red and blue galaxies by the CAMIRA member galaxy catalog with the same galaxy quality cuts.

We also employ a random catalog for the photometric galaxies from PDR2 with the same quality cuts when related quantities are available for the random catalog in PDR2. Specifically, we do not apply the |$\tt {cmodel\_flag}$|⁠, |$\tt {undeblended\_convolvedflux\_3\_15\_flag}$|⁠, |$\tt {cmodel\_mag}$|⁠, |$\tt {cmodel\_fluxsigma}$|⁠, |$\tt {cmodel\_magsigma}$|⁠, or |$\tt {extendedness\_value}$| cuts described above for the random catalog. The input number density of random points is 100 points arcmin⁻² when constructing the random catalog (Aihara et al. 2019). Because we use the conservative magnitude and flux signal-to-noise ratio (or magnitude error) cuts for the galaxy catalogs compared with the detection limits for the z-band described in subsection 2.1, the footprint of the random catalog follows the data footprint. We use this random galaxy catalog for measurements of projected cross-correlation functions in subsection 2.4.

2.4 Projected cross-correlation measurement

We apply an absolute-magnitude cut M_z − 5log₁₀h < −18.8 for our fiducial galaxy selection in subsection 4.1. This absolute magnitude limit corresponds to an apparent magnitude limit of z-mag = 24.5 at z_cl = 1.0. We use the mean redshift value z_{cl, i} for each subsample to convert the absolute magnitude to each apparent magnitude limit. We employ absolute-magnitude cuts in order to make fair comparisons between results at different redshifts. Nishizawa et al. (2018) employed a similar absolute-magnitude cut on z-band magnitude with the first-year HSC-SSP datasets. In subsection 4.2, we also use several different absolute-magnitude cuts to investigate how constraints on the splashback feature change mainly to investigate possible dynamical friction effects. Furthermore, in subsection 4.3, we apply a red-sequence-based red and blue separation cut in Nishizawa et al. (2018) for the photometric galaxies with a slightly brighter absolute-magnitude cut, M_z − 5log₁₀h < −19.3, to check the dependence of splashback features on galaxy colors.

We then estimate the covariance describing the statistical uncertainties of the projected cross-correlation measurements from the jackknife method (e.g., Norberg et al. 2009) following the literature (More et al. 2016; Baxter et al. 2017; Chang et al. 2018; Nishizawa et al. 2018; Shin et al. 2019; Zürcher & More 2019). We divide the cluster footprint into 200 rectangular subregions with almost equal areas, and we estimate the jackknife covariance by repeating the measurements by excluding each subregion for 200 times. The side length for each jackknife subregion corresponds to ∼10 h⁻¹ Mpc for clusters at z_cl = 0.1, which is larger than the maximum radial distance of our cross-correlation measurements. We calculate the signal-to-noise ratio of the measurements for each cluster and galaxy selection and list these values in the tables below. For the full cluster sample with the fiducial galaxy selection, the signal-to-noise ratio is 59.2, whereas the signal-to-noise ratio in More et al. (2016) is ∼250 with measurements from the larger survey area.

We note that the projected cross-correlation function is related to surface galaxy density profiles Σ_g after subtracting mean surface galaxy densities 〈Σ_g〉 as Σ_g(R) = 〈Σ_g〉ξ_2D(R). We estimate the mean surface galaxy density for each cluster and galaxy selection, which we also list in tables below for reference. First, we estimate mean surface density values in cluster redshift subsamples by employing the number of galaxies used for measurements after galaxy selections and survey areas in comoving coordinates from the median cluster redshift values. We then measure the averaged surface galaxy density profile |$\widehat{\Sigma }_{\rm g}(R)$| in a manner similar to equation (2) to estimate the mean surface density values for each cluster and galaxy selection as |$\langle \Sigma _{\rm g} \rangle ={\rm Median}[\widehat{\Sigma }_{\rm g}(R)/\widehat{\xi }_{ {\rm 2D} }(R)]$|⁠, where we use a median value from the radial bins for simplicity. As the mean surface density values do not affect constraints on the splashback feature, we employ the projected cross-correlation function for our analyses following Zürcher and More (2019). We use the mean surface galaxy density values to estimate three-dimensional galaxy density profiles in subsection 3.1.

3 Model

3.1 Profile model for projected cross-correlation measurement

We employ the constraints on the off-centering parameters f_cen and R_off in Murata et al. (2019) from stacked lensing profiles around the HSC CAMIRA clusters as the prior distributions in our splashback analysis. Note that Murata et al. (2019) assumed the same offset distribution as given by equation (10) when expressed in the projected two-dimensional space. Murata et al. (2019) used an empirical parametric model for f_cen with mean richness and redshift dependences in each cluster sample, and one R_off parameter in each redshift bin with the same redshift binning as our Low-z, Mid-z, and High-z samples. For the prior distribution in our analyses, we employ the posterior distribution in Murata et al. (2019) to calculate the constraints of f_cen from the mean richness and redshift values in table 1, and those of R_off using the redshift bin in Murata et al. (2019) with the closest mean redshift value for each cluster selection. Gaussian distributions can approximate the constraints well, and we derive their mean and standard deviation values for each cluster sample. We show these values in table 2. We ignore a weak correlation between the constraints of f_cen and R_off in Murata et al. (2019) for simplicity.

To summarize, our model of cross-correlation signals contains 10 parameters, including two parameters for off-centering effects. The model parameters and their prior distributions are given in table 2. The literature used Gaussian distributions for the prior distributions of log₁₀α, log₁₀β, and log₁₀γ, which were determined conservatively compared to the simulation results in Gao et al. (2008) and Diemer and Kravtsov (2014). For the Gaussian priors of log₁₀α, log₁₀β, and log₁₀γ, we follow moderate assumptions made by More et al. (2016) and Shin et al. (2019).¹ We use flat priors for the other model parameters in the DK14 profile with conservatively large widths, as shown in table 2. We employ a constant value for the mean of the log₁₀α parameter in table 2 for simplicity, ignoring mass and redshift dependence in simulation calibrations for dark matter detailed in Gao et al. (2008). We have confirmed that our results do not change significantly (⁠|${\lesssim}1\%$|–|$2\%$| for splashback locations and derivatives at these locations) even when using mean values for log₁₀α from Gao et al. (2008) with mean mass and redshift values for each sample with a fixed standard deviation.

3.2 Model fitting

4 Results

In subsection 4.1, we show results for the cluster samples in table 1 with the fiducial absolute-magnitude cut M_z − 5log₁₀h < −18.8 for galaxy selections. We present results for the Low-z cluster sample with different absolute-magnitude cuts for galaxy selections in subsection 4.2 to investigate possible dynamical friction effects. In subsection 4.3, we also show results with separations for red and blue galaxies and the absolute-magnitude cut M_z − 5log₁₀h < −19.3 over 0.1 ≤ z_cl ≤ 1.0.

4.1 Splashback features around HSC CAMIRA clusters with cluster redshift and richness dependences

We show our projected cross-correlation function measurements and the model predictions from the MCMC chains in the first row of figures 1 and 2 for each cluster sample in table 1 with the fiducial absolute-magnitude cut M_z − 5log₁₀h < −18.8 for galaxy selections. The second row shows constraints on the logarithmic derivative of the projected cross-correlation function based on our model fitting results. The third and fourth rows present constraints on the three-dimensional galaxy density profile and the logarithmic derivative of the three-dimensional profile, respectively. In the second and fourth rows, we also show the model predictions from the halo–matter cross-correlation function (i.e., dark matter distributions around clusters) in Nishimichi et al. (2019) and the mass–richness relation (Murata et al. 2019) presented in appendix  1 for comparisons.

Fig. 1.

In the top row, we show observed projected cross-correlation measurements between the clusters and galaxies (points with error bars in red color) from left to right, for the Full, Low-z, Mid-z, and High-z cluster samples from table 3 with the fiducial absolute-magnitude cut M_z − 5log₁₀h < −18.8 for galaxy selections. Shaded colored regions show the 16th and 84th percentiles of the model predictions from the MCMC chains for each cluster sample. The dashed profile curves in the top row show contributions from off-centered clusters to projected cross-correlation measurements at the best-fitting model parameters. The second row shows constraints on the logarithmic derivative of the projected cross-correlation profiles without off-centering effects in equation (7). The third row presents constraints on the three-dimensional galaxy density profile calculated as ρ_g(r) = 〈Σ_g〉ξ_3D(r)/(2R_max) (see subsection 3.1 for more details). The last row shows the constraints on the logarithmic derivative of the three-dimensional profiles. The vertical black dashed lines denote mean values of R_200m shown in table 1 for the cluster samples. Black vertical shaded regions show model predictions from the halo–matter cross-correlation (i.e., dark matter distributions around clusters) and mass–richness relation in appendix  1 for |$R_{\rm sp}^{\rm 2D}$| in the second row and |$r_{\rm sp}^{\rm 3D}$| in the third and last rows from table 1, whereas colored vertical shaded regions show their constraints from the MCMC chains in each cluster sample. Black shaded profile regions in the second and last rows show derivative profiles from the model predictions for each cluster sample. (Color online)

Open in new tabDownload slide

Fig. 2.

Same as figure 1, but for the Low-N, Mid-N, and High-N cluster samples from table 3. (Color online)

Open in new tabDownload slide

In table 3, we show constraints on the model parameters and derived splashback features. The constraints on log₁₀α, log₁₀β, log₁₀γ, f_cen, and R_off parameters are determined mainly from the Gaussian priors, whereas those for the other parameters are well constrained by the measurements compared to their prior distributions. The χ² values of the best-fitting model parameters in table 1 indicate that our fitting results are acceptable. Following Baxter et al. (2017), we list the values of the logarithmic derivative of the inner profiles [i.e., ρ_in(r)f_trans(r)] at the location of |$r_{\rm sp}^{\rm 3D}$| in table 1. These values are smaller than −3 with more than ∼2σ significance for all the samples. This result suggests that classical fitting functions like the Navarro–Frenk–White profiles (NFW; Navarro et al. 1996) have difficulty in fitting the observed cross-correlation measurements and provides the evidence of the splashback features (Baxter et al. 2017).

The 1σ precisions of our constraints on |$r_{\rm sp}^{\rm 3D}$| compared to the median values are approximately |$12\%$|⁠, |$16\%$|⁠, |$17\%$|⁠, |$25\%$|⁠, |$28\%$|⁠, |$18\%$|⁠, and |$13\%$| for the Full, Low-z, Mid-z, High-z, Low-N, Mid-N, and High-N selections, respectively, after marginalizing over the off-centering effects with the fiducial absolute magnitude galaxy cut. These constraints are consistent with the model predictions from the halo–matter cross-correlation function in appendix  1 within 0.5σ, 1.9σ, 0.2σ, 0.2σ, 0.5σ, 1.6σ, and 1.1σ levels for Full, Low-z, Mid-z, High-z, Low-N, Mid-N, and High-N, respectively, whereas our constraints are consistent with the |$20\%$| smaller model predictions (i.e., |$0.8 \times r_{\rm sp, model}^{\rm 3D}$|⁠) within 1.2σ, 0.2σ, 1.0σ, 0.9σ, 0.3σ, 0.2σ, and 0.7σ levels. Given the larger error bars due to the smaller survey area, our constraints are consistent with both the model predictions and their |$20\%$| smaller values. However, we note that our constraints on |$r_{\rm sp}^{\rm 3D}$| for the Low-z sample at 0.1 < z_cl < 0.4 are on the lower side than the model prediction at the level of 1.9σ, which is in line with the constraints in the literature (More et al. 2016; Baxter et al. 2017; Chang et al. 2018) with clusters selected by another optical cluster-finding algorithm redMaPPer. Further investigation is warranted given our larger error bars. We note that we compare constraints from only red galaxies with the model predictions again for the Low-z, Mid-z, and High-z cluster samples in subsection 4.3, as these constraints are more precise than those without separations of red and blue galaxies. Our constraints with the fiducial absolute magnitude galaxy cut in figures 1 and 2 suggest that the three-dimensional profiles are slightly steeper (i.e., smaller derivative values) around the splashback radius than the model predictions from dark matter. We discuss a possible origin for these results in subsections 4.2 by comparing constraints from different absolute magnitude galaxy cuts.

4.2 Splashback features from different absolute magnitude galaxy cuts

We employ the Low-z cluster sample to derive constraints on splashback features with different absolute magnitude galaxy cuts, as we can use fainter absolute-magnitude cuts for the galaxies with lower redshift clusters given a fixed apparent magnitude limit, and the constraint on |$r_{\rm sp}^{\rm 3D}$| for the Low-z sample is better among the cluster samples as shown in table 3 with the fiducial absolute-magnitude cut. Specifically, we shift the absolute-magnitude cuts by one or two magnitudes to both fainter and brighter sides with the z-band magnitude compared to the fiducial absolute-magnitude cut M_z − 5log₁₀h < −18.8 in subsection 4.1. The faintest absolute-magnitude cut (M_z − 5log₁₀h < −16.8) corresponds to an apparent z-band magnitude limit of 24.1 at z_cl = 0.4, which is shallower than the apparent magnitude cut for the galaxy catalog.

We show these constraints on the model parameters and splashback features in table 4 after marginalizing over the off-centering effects. We present our projected cross-correlation function measurements for the Low-z sample with the different absolute magnitude galaxy cuts and the model predictions from the MCMC chains in the left-hand panel of figure 3. The measurements with fainter absolute-magnitude cuts have lower amplitudes, indicating that brighter galaxies are more strongly clustered around clusters than fainter galaxies. The right-hand panel in figure 3 shows constraints on the three-dimensional galaxy density profile from the MCMC chains. We also show histogram comparisons of the constraints on splashback features in figure 4. Figure 4 indicates that the constraints on |$R_{\rm sp}^{\rm 2D}$| and |$r_{\rm sp}^{\rm 3D}$| are robust against the different absolute-magnitude cuts within the error bars. Specifically, the 1σ precisions of the constraints on |$r_{\rm sp}^{\rm 3D}$| compared to the median values are |$13\%$|⁠, |$15\%$|⁠, |$16\%$|⁠, |$20\%$|⁠, and |$30\%$| from the brightest to the faintest absolute-magnitude cut including the fiducial cut, respectively. Within these error bars, we do not see any systematic trend consistent with dynamical friction effects (i.e., smaller values for brighter galaxies), which is not surprising given expectations (a |$\lesssim 20\%$| shift of |$r_{\rm sp}^{\rm 3D}$| for our brightest galaxy sample and |${\lesssim}5\%$| shifts for other galaxy samples) from subhalos in simulations (More et al. 2016; Chang et al. 2018) with simplistic abundance matching methods. Note that the constraints on |$r_{\rm sp}^{\rm 3D}$| are preferred to be smaller than that predicted by the model across the five absolute magnitude range within these error bars, although moderate correlations exist among these constraints. Indeed, we confirm that there are moderate positive correlations in the jackknife cross-covariance for the projected cross-correlation measurements among these different magnitude cuts. Further investigation with reduced error bars from a larger survey area is warranted to constrain dynamical friction effects at an |$\sim 5\%$| level for the CAMIRA clusters.

Fig. 3.

We show the results from table 4 with the Low-z cluster sample and the different absolute-magnitude cuts for galaxy selections after marginalizing over the off-centering effects. The left-hand panel shows the projected cross-correlation measurements and the 16th and 84th percentiles of the model predictions from the MCMC chains for each selection. The right-hand panel presents constraints on the three-dimensional galaxy density profile ρ_g(r) = 〈Σ_g〉ξ_3D(r)/(2R_max) (see subsection 3.1 for more details) from the MCMC chains for each selection. (Color online)

Open in new tabDownload slide

Fig. 4.

We show the probability distribution functions of splashback features from the MCMC chains for the analyses with the Low-z cluster sample and the different absolute-magnitude cuts for galaxy selections in table 4 and figure 3. The vertical black lines denote the model predictions from the halo–matter cross-correlation emulator and mass–richness relation in appendix  1 for comparisons. The vertical dashed lines show the mean value of R_200m. (Color online)

Open in new tabDownload slide

On the other hand, the logarithmic derivatives of the three-dimensional profiles at |$r_{\rm sp}^{\rm 3D}$| are larger for fainter absolute-magnitude cuts, as shown in the right-hand panel of figure 3 and figure 4. More et al. (2016) shows larger logarithmic derivative values at |$r_{\rm sp}^{\rm 3D}$| with decreasing maximum circular velocities of subhalos in simulations, which is consistent with our results. We expect that these differences in the derivative values are mainly caused by different ratios between the inner and outer material number densities (=ρ_inf_trans/ρ_out) around the splashback radii, as splashback radii mark a halo boundary between the multi-stream and single-stream regions. We thus expect that derivative values would be smaller (i.e., steeper profiles) when the ratios are large at |$r_{\rm sp}^{\rm 3D}$|⁠, whereas these values should be larger (i.e., shallower profiles) when the ratios are smaller. We confirm this expectation by checking that constraints on ρ_inf_trans/(ρ_inf_trans + ρ_out) = [1 + (ρ_inf_trans/ρ_out)⁻¹]⁻¹ at |$r_{\rm sp}^{\rm 3D}$| for the fainter galaxy samples are systematically smaller than those for the brighter galaxy samples as presented in tables 3 and 4 although the error bars are not very small. Note that this correlation between the density ratios and the logarithmic derivative values at |$r_{\rm sp}^{\rm 3D}$| also exists in table 3 for the fiducial absolute-magnitude cuts and in table 5 for the color cuts presented in subsection 4.3.

In particular, figure 4 shows that the constraint on the logarithmic derivative at |$r_{\rm sp}^{\rm 3D}$| from the faintest absolute-magnitude cut is consistent with the model prediction from the dark matter distribution, suggesting that the ratio between the inner and outer materials for the faintest galaxies should be on average similar to the one for the dark matter distribution around splashback radii for the Low-z cluster sample. Also, our results might indicate that a larger fraction of fainter galaxies accrete into cluster regions (⁠|$r \lesssim r_{\rm sp}^{\rm 3D}$|⁠) compared to brighter galaxies at 0.1 < z_cl < 0.4, or that fainter galaxies have been selectively destroyed or merged compared to brighter galaxies to decrease number densities after infalls into cluster regions. We note that the blue fraction is higher at fainter magnitudes, thus the interpretation is also related to the results of subsection 4.3, although the absolute magnitude limits are not the same. It might be informative to compare these results with hydrodynamical simulations to constrain galaxy formation and evolution, possibly including relations between subhalo masses and galaxy luminosities.

4.3 Splashback features from red and blue galaxies

We employ a method from Nishizawa et al. (2018) to define red (quenched) and blue (star-forming) galaxy populations for each cluster redshift bin with the help of red-sequence measurements from the CAMIRA red member galaxy catalog described in subsection 2.2, without resorting to photometric redshifts of individual galaxies. This method uses a color–magnitude diagram for the CAMIRA member galaxies to derive a linear red-sequence relation at each cluster redshift within the redshift interval Δz_cl = 0.05 (see subsection 2.4), where we employ g–r, r–i, and i–y color for clusters at 0.1 < z_cl < 0.4, 0.4 < z_cl < 0.7, and 0.7 < z_cl < 1.0, respectively, and we use the z-band magnitude. These color combinations are chosen to be sensitive to the 4000 Å break from red galaxies at each redshift range. This method defines red galaxies with colors above a linear color–magnitude relation, which is 2σ lower (bluer) than the derived linear red-sequence mean relation at each cluster redshift bin and as a function of the z-band magnitude. We calculate these σ error values by averaging over the absolute magnitude limit for each cluster redshift bin, and these errors include the photometric errors as well as the intrinsic errors in the red sequence. With this method, we can define red and blue galaxy populations even outside the aperture size of the CAMIRA finder (R ≃ 1 h⁻¹ Mpc in physical coordinates). Nishizawa et al. (2018) showed that number ratios between correlated red and blue galaxy populations around clusters marginally change by ∼0.1 at the filter transitions for the color combinations (i.e., z_cl = 0.4, 0.7) at R < 1h⁻¹ Mpc. We have confirmed that a difference of these selection criteria with the PDR2 is negligibly small from those in Nishizawa et al. (2018) with the first-year HSC data. In our analysis, we use the galaxy catalog after the same selection for the CAMIRA member catalog with |$\tt {z\_cmodel\_mag - a\_z}<24.0$| and |$\tt {z\_cmodel\_magsigma}<0.1$|⁠, and we employ the absolute-magnitude cut M_z − 5log₁₀h < −19.3 as described in subsection 2.3.

We note that this method is different from those in the literature. Specifically, Baxter et al. (2017) and Shin et al. (2019) employed photometric redshifts of individual galaxies to select galaxies for deriving color–magnitude diagrams as a function of redshift and to infer red/blue separation criteria based on percentile ranges in color values. On the other hand, Zürcher and More (2019) inferred g–r color cuts to separate red from blue galaxy population as a function of redshift with spectroscopic samples, although these spectroscopic samples are quite shallow and thus consist of galaxies that are brighter than their galaxies for cross-correlation measurements.

In figure 5 and table 5, we show the projected cross-correlation measurements and derived constraints. Consistently with the literature (Baxter et al. 2017; Nishizawa et al. 2018; Shin et al. 2019; Zürcher & More 2019), we find that red galaxies are more concentrated around cluster centers than blue galaxies in the projected measurements. We detect splashback features for both the red and blue galaxy populations defined by our selection criteria over 0.1 < z_cl < 1.0, although the significances from the derivative of the inner profiles (Baxter et al. 2017) are smaller for blue galaxy populations of the Low-z and High-z cluster samples (∼1.5σ). Within the error bars, the constraints on |$r_{\rm sp}^{\rm 3D}$| are consistent among those for the fiducial absolute-magnitude cuts without the red and blue separation in subsection 4.1 and those for the red and blue galaxy populations at each redshift range. Due to their sharper profiles around |$r_{\rm sp}^{\rm 3D}$| (i.e., smaller derivative values), constraints on |$r_{\rm sp}^{\rm 3D}$| for only the red galaxies are more precise than those without the red and blue separation in subsection 4.1. Specifically, the 1σ precisions of the constraints on |$r_{\rm sp}^{\rm 3D}$| are approximately |$11\%$|⁠, |$13\%$|⁠, and |$16\%$| for the Low-z, Mid-z, and High-z cluster samples, respectively, when employing only the red galaxies, which are more precise than those in subsection 4.1 by |$\sim\! 10\%$| for the High-z cluster sample. These more precise constraints are consistent with the model predictions in table 1 within 1.0σ, 0.4σ, and 1.0σ levels, whereas these constraints are consistent with the |$20\%$| smaller model predictions within 1.0σ, 1.9σ, and 2.0σ levels for the Low-z, Mid-z, and High-z samples, respectively. We note that the constraints for the Mid-z and High-z cluster samples are more consistent with the model predictions (≲1σ levels) than the |$20\%$| smaller values of the model predictions (∼2σ levels). We also note that these analyses with the red galaxy populations approximately correspond to simulation results in Adhikari, Dalal, and Chamberlain (2014), which selectively uses simulation particles orbiting within clusters. They find that using simulation particles with small radial motions, which corresponds to preferentially selecting simulation particles within clusters, does not significantly bias the splashback radius measurement within the level of |$\sim 5\%$|–|$10\%$|⁠. Also, based on the analysis of N-body simulations, Shin et al. (2019) argue that the splashback radius depends on the accretion time (that may correlate with the galaxy colors) of the subhalos used to measure it, although the dependence appears to be weak except for subhalos accreted very recently. Given our larger error bars from the limited area, further investigations are warranted to constrain |$r_{\rm sp}^{\rm 3D}$| values at the higher redshift ranges (0.4 < z_cl < 1.0) more precisely with upcoming datasets.

Fig. 5.

Same as figure 1, but for the analyses with separations of red (red color) and blue (blue color) galaxies and an absolute-magnitude cut M_z − 5log₁₀h < −19.3 for galaxy selections in subsection 4.3 and table 5 for the Low-z, Mid-z, and High-z cluster samples after marginalizing over the off-centering effects. (Color online)

Open in new tabDownload slide

We expect that the density profiles of blue galaxies would be consistent with a pure power-law without any splashback features when these blue galaxies are still on their first infall passage (Baxter et al. 2017; Shin et al. 2019; Zürcher & More 2019). We thus employ a pure power-law model that only marginalizes over the off-centering effects in analyses for blue galaxies, instead of the DK14 model for the three-dimensional profile, to check whether the power-law model can have reasonable |$\chi _{\rm min}^2$| values or not. We find that these analyses result in |$\chi ^{2}_{\rm min}/{\rm dof}=29.2/13, 21.7/13, 31.9/13$| for the Low-z, Mid-z, and High-z cluster samples, respectively. Their low p-values suggest that the DK14 model (table 5) can fit the measurements better than the power-law model given the difference in the minimum χ² and the degree-of-freedom values over 0.1 < z_cl < 1.0 (i.e., |$\Delta \chi ^{2}_{\rm min} \gtrsim 10$| by adding the six parameters for the DK14 model with the Gaussian priors on three of them). Note that, as shown in figure 5, the logarithmic derivative of the three-dimensional profile at |$r=r_{\rm sp}^{\rm 3D}$| for the blue galaxy population at 0.4 < z_cl < 0.7 is smaller than those at the other redshift bins at the level of ∼2σ, possibly due to the different combinations of filters used to define the blue galaxy populations in Nishizawa et al. (2018) as described above.

At the same time, we also expect red galaxies to scatter into the blue population, especially at fainter magnitudes and higher redshift due to larger photometric errors. We estimate the contamination fraction from the red to the blue galaxy population to be f_red ≃ 0.025 because we employ the 2σ lower limit for the separation criteria in Nishizawa et al. (2018). Compared to Zürcher and More (2019), we employ the same galaxy selection for the cross-correlation measurements as the one to define the red/blue separation criteria over all the redshift range, in order to apply this contamination fraction more consistently for fainter galaxies at higher redshift. This contamination could potentially alter the cross-correlation measurements for the blue populations, causing the observed splashback feature in the blue population. We follow a method found in appendix C of Zürcher and More (2019) to assess such contamination. First, we assume a pure power-law function for the blue galaxies with S_e = 1.5 and the best-fitting parameter value for log₁₀ρ_o in table 5 at each redshift bin. We then calculate a weighted sum of this power-law profile and the DK14 model profile for the red galaxy population with the best-fitting parameters in table 5 with a contamination fraction from the red to the blue population at each redshift bin (i.e., ρ_{g, blue} + f_redρ_{g, red}). We find that the |$2.5\%$| contamination fraction alters the logarithmic derivative of the three-dimensional profile around |$r_{\rm sp}^{\rm 3D}$| of the red population by 0.1–0.2, which is smaller than values needed to explain the observed profiles (≳ 0.5 differences in the logarithmic derivative of the three-dimensional profiles around |$r_{\rm sp}^{\rm 3D}$|⁠) in figure 5. We would require |$10\%$|–|$20\%$| contamination fractions from the red to the blue galaxy population to be consistent with the observed profiles for the blue populations when we assume a pure power-law function before such contaminations.

From these results, consistently with Zürcher and More (2019) from the different selection method on the red and blue galaxy populations, we conclude that some fraction of blue galaxies stay blue even after reaching their first apocenters of orbits within host clusters (i.e., splashback radii for these clusters) from the presence of the splashback feature in the blue galaxy populations when we define the blue galaxy populations using the method in Nishizawa et al. (2018). We note that Shin et al. (2019) found that their bluest galaxy sample defined by their criteria is consistent with a pure power-law, which is consistent with their bluest galaxy sample being on their first infall into their cluster samples. We expect that the difference from our results might come from the difference of the selection criteria of red and blue galaxy populations, as also discussed in Zürcher and More (2019).

Also, we investigate how the constraints on three-dimensional galaxy density profiles of the red and blue galaxy populations change around the locations of |$r_{\rm sp}^{\rm 3D}$| more directly in figure 5, whereas Baxter et al. (2017) and Shin et al. (2019) showed that their red/blue fractions abruptly change around |$r_{\rm sp}^{\rm 3D}$| in their projected profiles. Figure 5 shows that the three-dimensional galaxy density fractions between the red and blue galaxy populations change rapidly around the model predictions of |$r_{\rm sp}^{\rm 3D}$| in table 1 for all the redshift bins and that the blue galaxy populations in our separation criteria dominate at the outskirts of clusters over 0.1 < z_cl < 1.0 after the 15 < N < 200 selection. The constraints on three-dimensional galaxy densities in figure 5 do not evolve significantly within the error bars over 0.1 < z_cl < 1.0 for both the red and blue galaxy populations around clusters when selected by 15 < N < 200, including concentrations for red galaxy populations, possibly due to our off-centering marginalization compared to Nishizawa et al. (2018). In particular, we find that the three-dimensional galaxy densities for the red and blue populations become the same value at |$r=1.63^{+0.24}_{-0.15}\, h^{-1}\:\mbox{Mpc}$|⁠, |$1.50^{+0.29}_{-0.20}\, h^{-1}\:\mbox{Mpc}$|⁠, and |$1.76^{+0.28}_{-0.31}\, h^{-1}\:\mbox{Mpc}$| for the Low-z, Mid-z, and High-z cluster sample, respectively, where we calculate these values from the MCMC chains for each analysis. We note that we ignore positive cross-covariances between the projected cross-correlation measurements for the red and blue galaxy populations at each redshift bin for simplicity, and thus our error bars on these radial values are conservatively larger than those with these positive cross-covariances. These constraints are consistent with the model predictions of |$r_{\rm sp}^{\rm 3D}$| in table 1 within these error bars (≲1σ level) for all the redshift bins.