LoCuSS: The Splashback Radius of Massive Galaxy Clusters and Its Dependence on Cluster Merger History

The outskirts of galaxy clusters are the new frontier to improve our understanding of the multifaceted physical processes impacting baryons and dark matter during structure formation. They mark the transition between pristine field sparseness and evolved collapsed halos. Here merger shocks are produced and galaxy evolution is jolted, due to accretion onto clusters (Walker et al. 2019). This makes cluster outskirts the prime targets to observe the transformation of galaxy properties (Haines et al. 2015), to assess biases in clusters' mass estimates (Reiprich et al. 2013), and test predictions of cosmological models (Pratt et al. 2019).

The collapse and growth of dark matter halos follows the evolution of primordial density perturbations whose gravitational pull locally overcomes the expansion of the universe. In the simplest scenario involving spherically symmetric, continuous accretion, virialization during collapse redistributes the potential energy of the primordial overdensities, allowing for a final stable structure characterized by a specific density contrast Δ_vir (Gunn & Gott 1972, see also Voit 2005 for a review). It follows that the mass enclosed within the virialized sphere is ${M}_{\mathrm{vir}}=\tfrac{4}{3}\pi {r}_{\mathrm{vir}}^{3}{{\rm{\Delta }}}_{\mathrm{vir}}\bar{\delta }$ (Cole & Lacey 1996). This definition is commonly used to indicate gravitationally bound structures, while allowing for different choices of Δ ∈ [180, 200, 500] and reference density $\bar{\delta }$, namely the mean matter density or critical density of the universe. This formalism is particularly advantageous since it permits a bridge between simulations and observations, and it has been used to probe self-similarity for a wide range of halo properties, involving both baryonic and non-baryonic matter, in regimes that are dominated by gravity (e.g., Navarro et al. 1996; Bullock et al. 2001; Tinker et al. 2008; Schaller et al. 2015; Springel et al. 2018; Farahi et al. 2019).

Modifications have been proposed to the idealized collapse model to allow for more general conditions of structure growth, such as those in ΛCDM cosmology (Bryan & Norman 1998). However, it is a challenge to fully describe the structure on the outskirts of these halos, with models generally unable to capture the rich extent of density fields that connect halos that are gravitationally bound despite being separated by many virial radii (Prada et al. 2006). Indeed, as the density contrast is defined against mean cosmic densities that have a redshift dependence, an intrinsic pseudo-evolution of halo masses arises from the formalism (Babul et al. 2002; Diemand et al. 2007; Diemer et al. 2013). Naively, the spherical collapse model suggests the presence of a sharp transition in the radial density profile of the forming halo, corresponding to the region that physically separates gravitationally bound matter accreted previously from newly infalling regions (Fillmore & Goldreich 1984; Bertschinger 1985; Adhikari et al. 2014; Shi 2016). Infalling particles that get trapped in the growing halo potential populate a shell at the apocenter of their first orbit—the so-called splashback region. For this reason, the splashback radius has been proposed as a physical boundary of collapsed halos and is detected as a steepening of the radial density profile of dark matter halos from large-scale numerical simulations (Diemer & Kravtsov 2014; More et al. 2015, see also Tomooka et al. 2020).

Observationally, the challenge of detecting the splashback radius has been recently undertaken using galaxies and gas measured in large-scale surveys as tracers of the underlying dark matter halo potential in galaxy clusters. Unfortunately, it has become apparent that the measured splashback radius can be significantly affected by the cluster selection method. Optical selection methods based on cluster richness measured within an aperture result in an underestimated splashback location due to projection effects and interloper contamination (Baxter et al. 2017; Busch & White 2017; Zu et al. 2017). Additional biases on the measure of the splashback location have been confirmed for methods relying on clusters selected using redMaPPer (Rykoff et al. 2014) by comparison with cluster mass profiles obtained using weak-lensing shear (Chang et al. 2018). More recently, selection methods that are unaffected by the same biases, and thus minimize the risk of spurious correlations between the splashback radius and the cluster selection, such as Sunyaev–Zel'dovich (SZ, Sunyaev & Zeldovich 1972; Shin et al. 2019; Zürcher & More 2019; Adhikari et al. 2020) and X-ray (Umetsu & Diemer 2017), have been attempted. To date, these methods have yielded less precise measurements of the splashback radius due to their smaller samples (More et al. 2015; Baxter et al. 2017; Murata et al. 2020), and to the typically noisy measurements from weak-lensing at these radii.

Simulations indicate that cluster outskirts are continuously bombarded by smaller halos, galaxy- or group-like, which is the prevalent ingredient for the fast evolution of cluster mass, especially below redshift z < 0.5 (McGee et al. 2009; De Lucia et al. 2012). The overall mass accretion rate is shown to impact the measured cluster density profiles. In particular, high accretion rates translate into steep slopes of the density profiles, even in projection at radii r > 0.5 R_200,m (Diemer & Kravtsov 2014). In addition, fast accretion impacts clusters across their entire extent, inducing rapid growth in both r_s and r₅₀₀, with respect to quieter clusters whose inner regions remain undisturbed while still accreting mass in their outskirts (Mostoghiu et al. 2019). Confirmation of this ongoing accretion is given from the increasing efforts of observations to census infalling halos at their first encounter with the cluster potential, such as ram-pressure stripped galaxies (Poggianti et al. 2016) and groups (Eckert et al. 2014). Therefore, a systematic cluster campaign targeting infalling halos enables us simultaneously to address galaxy evolution, cluster physics, and physics of hierarchical assembly (e.g., Bianconi et al. 2018; Haines et al. 2018). An observational study of a well-defined cluster sample that would enable direct comparison with theoretical predictions of the impact of merger history on the splashback feature has not been attempted so far.

In this work, we present the detection of the splashback feature in a sample of massive X-ray selected clusters at redshifts of 0.15 < z < 0.3 selected from the Local Cluster Substructure Survey (LoCuSS). Specifically, we take advantage of the rich multiwavelength LoCuSS data set and highly complete spectroscopic follow-up observations from the Arizona Cluster Redshift Survey (ACReS) to detect the feature in stacked projected density profile of the clusters, constrain the de-projected splashback radius of the whole sample, and examine how the splashback radius depends on observables that probe the assembly history of the clusters. The use of spectroscopically confirmed cluster member galaxies, and dissection of the cluster sample based on their properties, are important new steps that can help to shape future investigations of cluster assembly with upcoming surveys including the Vera Rubin Observatory's Legacy Survey of Space and Time, ESA's Euclid satellite, ESO's 4MOST instrument, and the German–Russian eROSITA X-ray satellite.

In Section 2 we define the cluster sample, give an overview of the LoCuSS/ACReS data set, and discuss the structure of the clusters in our sample. In Section 3, we describe the details of how we approach the modeling of the data. Section 4 contains our results, beginning with the empirical detection of the splashback feature in the stacked projected density profile of the clusters, before applying the Bayesian modeling scheme described in Section 3 to infer the splashback radius in three-dimensions, and finally examining how the results depend on the structure of the clusters. We close by summarizing and discussing our results in Section 5. We assume cosmology values presented in Planck Collaboration et al. (2016), with h = 0.678, H₀ = 100 h km s⁻¹ Mpc⁻¹, Ω_m = 0.309, and Ω_Λ = 0.691.

The Local Cluster Substructure Survey (LoCuSS, PI: G. P. Smith) is a multiwavelength survey of X-ray selected massive galaxy clusters from the ROSAT All Sky Survey (Ebeling et al. 1998; Böhringer et al. 2004). We study 20 clusters that are the union of the complete high-L_X sample of 50 clusters (Okabe et al. 2013; Martino et al. 2014; Okabe & Smith 2016; Smith et al. 2016; Farahi et al. 2019; Mulroy et al. 2019), and the galaxy evolution sample of 30 clusters (Smith et al. 2010; Haines et al. 2012, 2013, 2015, 2018; Bianconi et al. 2018). Therefore, the sample of 20 clusters considered here are an unbiased subsample of the complete high-L_X sample, for which high quality wide-field multiwavelength data are available. For reference, the high-L_X sample was selected using the following criteria: 0.15 ≤ z ≤ 0.3, L_X(0.1 − 2.4 keV)/E(z) ≥ 4.2 × 10⁴⁴ erg s⁻¹, and −25 < δ[deg] < +65, n_H ≤ 7 × 10²⁰cm⁻². The 20 clusters are listed in Table 1.

Table 1. List of Clusters Selected for This Study, Presented Together with Central Coordinates, Average Redshift, Weak-lensing Mass M_200,m from Okabe & Smith (2016), and Membership of the Main Subsamples Discussed in the Text, and Based on Haines et al. (2018) and Sanderson et al. (2009)

Name	R.A.(J2000)	Decl.(J2000)	Redshift	M200,m	Infalling	Low Entropy
			〈z〉	1014 h−1 M⊙	Groups?	Core?
Abell 68	00:37:06.84	+09:09:24.28	0.251	${8.39}_{-1.64}^{+2.00}$	Y	N
ZwCl 0104.4 + 0048	01:06:49.50	+01:03:22.10	0.253	${2.98}_{-1.26}^{+2.21}$	N	Y
Abell 209	01:31:53.45	−13:36:47.84	0.209	${17.01}_{-2.93}^{+3.70}$	Y	N
Abell 383	02:48:03.42	−03:31:45.05	0.189	${6.90}_{-1.64}^{+2.18}$	N	Y
Abell 586	07:32:20.22	+31:37:55.88	0.171	${8.32}_{-2.32}^{+3.54}$	N	N
Abell 611	08:00:56.81	+36:03:23.40	0.286	${11.77}_{-2.35}^{+2.77}$	Y	N
Abell 697	08:42:57.58	+36:21:59.54	0.282	${14.22}_{-3.73}^{+6.14}$	Y	N
ZwCl 0857.9 + 2107	09:00:36.86	+20:53:39.84	0.234	${3.52}_{-1.39}^{+1.97}$	N	Y
Abell 963	10:17:03.65	+39:02:49.63	0.204	${9.46}_{-1.79}^{+2.20}$	Y	Y
Abell 1689	13:11:29.45	−01:20:28.32	0.185	${13.15}_{-1.97}^{+2.32}$	N	N
Abell 1758	13:32:33.50	+50:30:31.61	0.279	${7.22}_{-1.83}^{+2.42}$	Y	N
Abell 1763	13:35:18.07	+40:59:57.16	0.232	${22.89}_{-4.32}^{+5.94}$	Y	N
Abell 1835	14:00:52.50	+02:52:42.64	0.252	${12.27}_{-2.28}^{+2.75}$	Y	Y
Abell 1914	14:25:59.70	+37:49:41.63	0.167	${12.51}_{-2.65}^{+3.55}$	Y	N
ZwCl 1454.8 + 2233	14:57:15.11	+22:20:34.26	0.257	${6.28}_{-2.69}^{+6.10}$	Y	Y
Abell 2219	16:40:22.56	+46:42:21.60	0.226	${15.17}_{-3.16}^{+4.53}$	Y	N
RXJ 1720.1 + 2638	17:20:10.14	+26:37:30.90	0.160	${7.23}_{-2.26}^{+3.46}$	Y	Y
RXJ 2129.6 + 0005	21:29:39.88	+00:05:20.54	0.234	${7.35}_{-2.48}^{+4.11}$	N	Y
Abell 2390	21:53:36.85	+17:41:43.66	0.229	${13.75}_{-2.42}^{+2.91}$	Y	Y
Abell 2485	22:48:31.13	−16:06:25.60	0.247	${7.56}_{-1.74}^{+2.27}$	N	N

Download table as:  ASCIITypeset image

The full wide-field multiwavelength data set on the sample of 20 clusters spans X-ray to millimeter wavelengths, and includes data from Chandra, XMM-Newton, GALEX, the Subaru 8.2 m telescope, the Hectospec instrument on the Multiple Mirror Telescope (MMT), UKIRT/WFCAM, the Mayall 4 m telescope, WISE, Spitzer/MIPS, the PACS and SPIRE instruments on Herschel, the Sunyaev Zeldovich Array, and ESA's Planck satellite. We concentrate on the X-ray, optical, and near-infrared data in this article, and summarize the cluster properties derived from these data, and previously published LoCuSS articles that contain detailed information in Table 2.

Table 2. Guide to the LoCuSS Data Sets Used in This Article

Telescope:	XMM-Newton	Chandra	Subaru	MMT	UKIRT
Instrument:	EPIC	ACIS-I	SuprimeCAM	Hectospec	WFCAM
Energy / wavelength:	0.5–2.4 keV	0.3–7 keV	$V/i^{\prime}$-bands	400–900 nm	J/K-bands
Physical radius probed:	1.5r200,m	0.8r200,m	1.5r200,m	3r200,m	2.5r200,m
Sensitivity:	10−14 erg s−1 cm−2	10−13 erg s−1 cm−2	V(5σ) = 27.5	90% complete	J(5σ) = 23
			$i^{\prime} (5\sigma )=26$	at K ≤ K⋆ + 1.5	K(5σ) = 22

Measurements:	Data Sets Used					References
Density profile				✓	✓	This paper
Total cluster mass	✓	✓	✓			1, 2
Infalling groups	✓			✓	✓	3, 4
Central entropy	✓	✓				5
X-ray centroid shift	✓	✓				1
X-ray concentration	✓	✓				1
Luminosity gap			✓	✓	✓	6
X-ray/BCG offset	✓	✓		✓	✓	6

Note. The references listed in the final column are as follows: (1) Martino et al. (2014); (2) Okabe & Smith (2016); (3) Haines et al. (2018); (4) Bianconi et al. (2018); (5) Sanderson et al. (2009); (6) Mulroy et al. (2019).

Download table as:  ASCIITypeset image

Highly complete spectroscopic follow-up of stellar-mass selected candidate galaxy cluster member galaxies out to ≈3r_200,m from the Arizona Cluster Redshift Survey (ACReS ⁶ ) is central to our analysis and results. This is because we aim to use spectroscopically confirmed member galaxies as test particles to measure the density profile of the clusters on scales well outside the physical scale on which the splashback feature is expected to be found. The combined UKIRT/WFCAM and MMT/Hectospec data set allows us to do this to projected clustercentric radii of 6–7 Mpc. The cluster mass measurements from weak-lensing analysis of Subaru observations reach a typical radius of r_200,c (Okabe & Smith 2016). This, combined with the mass information available to us from the X-ray analysis, allows us to experiment with different approaches to stacking the cluster luminosity density profiles—i.e., in physical comoving units, or scaled to a common overdensity radius such as r_200,m. A visual impression of the quality of the data is given in Figure 1, and full details of density profile construction are described in Section 4.

Zoom In Zoom Out Reset image size

The combined LoCuSS/ACReS data set also allow us to characterize the structure of clusters. This is important, because the splashback feature in cluster density profiles is associated with infall of material into the cluster potential, and the structure of clusters seen at different wavelengths is sensitive to their infall history. In addition to stacking all of the 20 clusters, we therefore present results in Section 4 based on sub-dividing the sample in to subsamples selected based on their structural properties that are motivated by infall history. The main results are based on subsamples selected on the presence of X-ray detected infalling galaxy groups, which indicate active ongoing accretion of galaxies, based on the analysis of Haines et al. (2018). Infalling X-ray groups were identified as extended X-ray sources in our XMM images down to ≈2 × 10⁴² erg s⁻¹ (M₂₀₀ ≈ 2.5 × 10¹³ M^⊙) and confirmed to be infalling into the clusters by identifying group members within the X-ray emitting region of the group with spectroscopic redshifts consistent with that of the host cluster. The X-ray data typically permit us to identify infalling groups within the radial range 0.35–1.3 r₂₀₀. In addition we included our own analysis of three clusters not in their sample (A 586, A 2485, and ZwCl 0104.4+0048). We excluded Abell 1689 from the subsample of clusters without infalling groups because of the rich structure of this cluster along the line of sight, which is strongly indicative of ongoing merging activity (Miralda-Escude & Babul 1995; Andersson & Madejski 2004; Peng et al. 2009). In addition, we use proxies for cluster merger history including central entropy, X-ray/BCG offset, X-ray brightness concentration and centroid shift, and luminosity gap, as discussed in Appendix A).

Empirical profiles have proven successful in describing the radial variations of density profiles of dark matter halos, featuring single or multiple power-law behaviors. Among the latter, Jaffe (Jaffe 1983), Hernquist (Hernquist 1990), and Navarro–Frenk–White (Navarro et al. 1996) share the same functional form, differing only in the inner and outer power-law slopes (see Binney & Tremaine 2008 for a review). Here, we choose to model the inner slope of the density profile with a single power-law prescription (Einasto et al. 1974), corresponding to the virialized part of the halo, together with an additional power-law describing the infall region. These two components are linked by a transitional term that allows for the radial profile to steepen due to the presence of the splashback feature. Following Diemer & Kravtsov (2014) and using the Colossus toolkit (Diemer 2018), we adopt the radial density profile ρ(r) accordingly:

$$\begin{eqnarray}&&\rho (r)={\rho }_{\mathrm{Ein}}(r)\times {f}_{\mathrm{trans}}(r)+{\rho }_{\mathrm{infall}}(r),\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{\rho }_{\mathrm{Ein}}(r)={\rho }_{s}\,\exp \left\{-\displaystyle \frac{2}{\alpha }\left[{\left(\displaystyle \frac{r}{{r}_{s}}\right)}^{\alpha }-1\right]\right\},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{f}_{\mathrm{trans}}(r)={\left[1+{\left(\displaystyle \frac{r}{{r}_{t}}\right)}^{\beta }\right]}^{-\gamma /\beta },\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{\rho }_{\mathrm{infall}}(r)=\displaystyle \frac{{\rho }_{0}}{1/{\rm{\Delta }}+{\left(r/{r}_{\mathrm{pivot}}\right)}^{{s}_{e}}},\end{eqnarray} \tag{ 4 }$$

with r_pivot = 1.5 Mpc. The parameters ρ_s and ρ₀ are inner and infall scale density, respectively. The parameter r_s is the scale radius of the inner profile and r_t is the transition radius between the inner and infall regimes. The exponent α sets the slope of the inner profile, β and γ set the shape and depth of the transition term, and s_e sets the slope of the outer term (Diemer & Kravtsov 2014). The transition function helps in reconstructing the density profile around the splashback feature, which gets smoothed by the radial averaging (Mansfield et al. 2017). The parameter Δ acts as a maximum cutoff density of the outer term to avoid its spurious contribution at small radii. Overall, implementing such truncation has a negligible impact, and we leave this threshold unconstrained. Hence α, β, γ, r_s , r_t , ρ_s , ρ₀, and s_e are the free parameters of the model. We considered also models without transition by setting the parameter γ = 0. We relate the 2D projected density Σ(R) observed from LoCuSS to the 3D density ρ(r) via

$$\begin{eqnarray}&&{\rm{\Sigma }}(R)=2{\int }_{R}^{\infty }{\rm{d}}l\,\displaystyle \frac{\rho (r){rdr}}{\sqrt{{r}^{2}+{R}^{2}}},\end{eqnarray} \tag{ 5 }$$

where R is the 2D projected distance from the cluster center. The splashback radius, r_sp, is a derived parameter in this model, and represents the minimum of the logarithmic derivative $d\,\mathrm{Log}(\rho )/d\,\mathrm{Log}(r)$ of the density profile.

We perform the estimation of best-fit parameters in the context of Bayesian inference, by computing the posterior distribution

$$\begin{eqnarray}&&P({\boldsymbol{\theta }}| d,m)=\displaystyle \frac{{ \mathcal L }({\boldsymbol{\theta }})P({\boldsymbol{\theta }}| m)}{{ \mathcal Z }},\end{eqnarray} \tag{ 6 }$$

where ${ \mathcal L }$ is the likelihood, P( θ ∣m) is the prior, and ${ \mathcal Z }$ is the evidence, for a given model m(x_i ∣θ) evaluated at the data point x_i using parameter values specified by θ . We adopt a Gaussian likelihood ${ \mathcal L }$ for the data d = (x_i , y_i ), given θ :

$$\begin{eqnarray}&&\mathrm{log}{ \mathcal L }({\boldsymbol{d}}| \theta )=-\displaystyle \frac{1}{2}\displaystyle \sum _{i}^{{N}_{d}}{\left(\displaystyle \frac{{y}_{i}-m({x}_{i}| \theta )}{{\sigma }_{i}}\right)}^{2}+\mathrm{log}{ \mathcal N }({N}_{d},\{{\sigma }_{i}\}),\end{eqnarray} \tag{ 7 }$$

where σ is the uncertainty on the data, and i cycles over the number of data points. ${ \mathcal N }({N}_{d},\{{\sigma }_{i}\})$ acts as normalization constant to ensure that the integral of the posterior distributions equals to unity.

In this work, we utilize the nested sampling algorithm to compute the Bayesian evidence ${ \mathcal Z }$, following the formalism by Skilling (2004) implemented in CPNest (Veitch et al. 2017). Nested sampling advantages include the simultaneous estimates of both evidence and posterior samples, and in being easily parallelizable. By splitting the prior volume over intervals of equal likelihood, the evidence results from simple numerical integration, with the addition of reducing the dimensionality of the problem to 1D. The algorithm allows us to evaluate the posteriors on the model parameters

$$\begin{eqnarray}&&\mathrm{ln}{ \mathcal P }({\boldsymbol{\theta }}| d,m)=\mathrm{ln}\left[{ \mathcal L }({\boldsymbol{d}}| {\boldsymbol{m}}({\boldsymbol{\theta }}))\Pr ({\boldsymbol{\theta }}| m)\right].\end{eqnarray} \tag{ 8 }$$

Nested sampling allows us to perform model selection based on marginalized likelihood. Therefore we can directly quantify and compare the evidence in the data in favor of models with and without the transition feature. We consider eight free parameters (ρ₀, ρ_s , r_t, r_s, α, β, γ , and s_e) from the cluster profile model. While the number of free parameters is large relative to the number of data points, our primary intention here is not to extract robust constraints on the model parameters, but rather to use the model fits to smoothly interpolate the data to extract constraints on the density profile. Despite this, we are able to run to convergence all the fitting procedures using flat priors (see Table 3), without imposing additional a priori knowledge on the parameter distributions, and to extract data-driven information.

4.1. Constructing the Projected Density Profile

4.2. Empirical Detection of the Splashback Feature

Figure 2 (top left panel) shows the luminosity density profile of the total cluster sample. A dip in the density profile can be noted around 4 Mpc. This feature appears more clearly (bottom left panel) when plotting the logarithmic slope of the profile $\epsilon =d\,\mathrm{Log}(\rho )/{\text{}}d\,\mathrm{Log}(R)$, and extends between 3 and 5 Mpc, peaking at 4 Mpc, where it reaches = − 2.30 ± 0.06. This has a significance of 5.9σ with respect to the mean slope value of = − 1.55 ± 0.11 at neighboring radii. The location of this feature in the density profile, and the run of slope with radius that we obtain are typical of what has been predicted from numerical simulations (Diemer & Kravtsov 2014), even when considering the 2D surface density. In Figure 2, we show also the density profiles of the cluster sample, split according to the presence of infalling X-ray groups in their surroundings. Interestingly, we notice that the sharp splashback feature, which clearly appears in both subsamples occurs at smaller radii when considering clusters without infalling X-ray groups. In particular the splashback feature peaks around 2.5 Mpc with a slope of = − 2.23 ± 0.07 at 6.5σ for clusters with no infalling groups and peaks at 3.8 Mpc a slope of = − 2.37 ± 0.06 at 6.1σ for the systems with infalling groups. Furthermore, we notice that the splashback feature appears consistently at the same location, whether considering scaled or nonscaled radii, for clusters without infalling groups. This is consistent with clusters not actively accreting groups having more self-similar structure than clusters that are actively accreting.

Zoom In Zoom Out Reset image size

We note a similar picture when classifying the clusters using the different structural parameters (see Appendix A: Figures 6 and 7). Specifically, cluster subsamples that are discussed in the literature variously as relaxed, undisturbed, or dynamically quiet (e.g., based on low central entropies, small X-ray centroid shift, and large luminosity gap) have a more prominent splashback feature, appearing at smaller clustercentric radii than their so-called unrelaxed, disturbed, or dynamically more active cousins. The consistency of this picture is very striking because the structural parameters used to define the different subsamples span a wide range of scales, from central entropy on scales of 20 kpc through to the presence of X-ray emitting infalling groups at 1–3 Mpc.

4.3. Measurement of Splashback Radius in Three-dimensions

We performed the model fit to the entire cluster sample and to subsamples classified by the presence and absence of infalling groups, as shown in Figure 3. The best-fit parameters, together with salient properties of the best-fit model are summarized in Table 4. From the full cluster sample, we recover a ratio between the 3D splashback radius and r_200,m of r_sp/r_200,m = 1.74 ± 0.34 and cluster masses consistent with Okabe & Smith's weak-lensing analysis.

Zoom In Zoom Out Reset image size

Table 4. Parameters of the Best-fit Model Resulting from the Fit of the Full Cluster Sample, and Classified According to the Presence of Infalling Groups

Parameter	Sample
and Marginalized Posterior	Total	Without	With
		Infalling Groups	Infalling Groups
ρ0 [103 × M⊙ h2/kpc3]	${4.43}_{-1.49}^{+1.41}$	${2.72}_{-0.88}^{+0.83}$	${4.11}_{-2.05}^{+2.01}$
se	${1.75}_{-0.16}^{+0.15}$	${1.61}_{-0.16}^{+0.15}$	${1.46}_{-0.20}^{+0.19}$
rs [Mpc/h]	${0.48}_{-0.01}^{+0.01}$	${0.48}_{-0.05}^{+0.05}$	${0.58}_{-0.02}^{+0.02}$
ρs [105 × M⊙ h2/kpc3]	${246.5}_{-24.4}^{+22.1}$	${150.7}_{-33.6}^{+30.3}$	${186.3}_{-12.3}^{+12.3}$
rt [Mpc/h]	${4.18}_{-0.87}^{+1.00}$	${1.81}_{-0.15}^{+0.14}$	${2.78}_{-0.12}^{+0.11}$
α	${0.27}_{-0.02}^{+0.02}$	${0.19}_{-0.05}^{+0.05}$	${0.20}_{-0.02}^{+0.02}$
β	>6	>6	>6
γ	>4	>4	>4
r200,m [Mpc/h]	2.19 ± 0.01	1.74 ± 0.01	2.35 ± 0.01
M200,m [1014 × M⊙/h]	14.11 ± 0.12	6.84 ± 0.13	17.64 ± 0.18
rsp [Mpc/h]	3.83 ± 0.75	2.01 ± 0.12	3.04 ± 0.15
max(Log slope 3D)	$-{3.4}_{-0.6}^{+0.7}$	$-{4.2}_{-0.2}^{+0.7}$	$-{4.4}_{-0.2}^{+0.9}$

Note. The errors quoted in the individual parameters are estimated from their posterior distribution and encompass the 15th and 85th percentile. The table includes also marginalized posteriors, i.e., r_200,m, M_200,m, and 3D splashback radius r_sp and their ± 3σ intervals. In addition we quote the minimum values of the 3D logarithmic slope of the density profile.

Download table as:  ASCIITypeset image

As expected based on the results in Section 4.2, the three-dimensional splashback radius of clusters without infalling groups, r_sp/r_200,m = 1.158 ± 0.071, is smaller than for clusters with infalling groups r_sp/r_200,m = 1.291 ± 0.062. This difference between the subsamples is significant at 4.2σ, and persists if we remove the spectroscopically confirmed infalling group members from the cluster member sample. Masses M_200,m for both halos with and without infalling groups are in agreement with the ones recovered from the weak-lensing analysis by Okabe & Smith (2016). We stress that our data set is sensitive to the detection of massive infalling groups and retains completeness down to galaxy stellar masses of M_⋆ = 2 × 10¹⁰ M_⊙ (Haines et al. 2018). Therefore, sampling the full mass range of halo accretion on clusters in beyond the reach of this data set, hence a direct comparison with simulations remains challenging because predictions that match our observational sample are not yet available.

Among the parameters considered in the fit, we note looser constraints obtained for β and γ. These parameters are known to be related to the accretion rate of halos, as shown by numerical simulations (Diemer & Kravtsov 2014), and can jointly span the prior space even in case of halos with similar properties. This degeneracy is known and has been mitigated by imposing stringent log-normal priors in the literature (Shin et al. 2019; Murata et al. 2020). We have chosen to not adopt this restriction in our study, after verifying the low impact of the flat priors choice on the fit. In particular, we obtain density profile slopes that are below values of −3, which is the lower limit of NFW profiles, when considering the full extent of the posteriors obtained for β and γ. This result provides further evidence of the splashback feature. The corner plot presenting the distribution of the posteriors of each parameter are shown in Figure 8 in Appendix A.

The Bayesian framework that we used through the CPNest implementation allowed us to directly compare the models, with and without splashback features, and determine which one is preferred according to the data-driven information. In particular, we can compute directly the Bayes factor ${ \mathcal B }$ from the evidence of the two models extracted from the fitting procedure in CPNest, under the assumption of equal and uniform priors. We report a Bayes factor in excess of ${ \mathcal B }\gt 100$ in favor of the model with a splashback feature when considering the profiles of clusters without infalling groups. This is the data set showing the strongest signature of the splashback feature among the ones considered here. The fit to the model without feature outputs a skewed posterior of the ρ₀ parameter toward the lower limit of the prior. This limit corresponds to the critical matter density of the universe at the redshift considered, and bounds the density domain of the infalling part of the model (Diemer & Kravtsov 2014). This is further evidence of the necessity of a model including a density transition to describe the data profiles. Furthermore, by fitting the same models with the same priors to the data excised of the transition region, we obtain a Bayes factor ${ \mathcal B }\lt 0.6$, favoring the model not allowing for the splashback transition.

5.1. Summary of Results

5.2. Comparison with Simulations

In Figure 4, we can see the trend of r_sp/r_200,m with respect to redshift. For comparison, the background lines mark the theoretical trend of the ratio r_sp/r_200,m as a function of halo accretion rate, which is defined as the logarithmic variation of the cluster mass within one dynamical timescale (≈1 Gyr) from Diemer et al. (2017). Simulations have shown that actively accreting clusters present a contraction of the splashback radius, correlating with accretion rate. In this work, we observe the opposite behavior. Sorce et al. (2020) analyzed a set of cluster halos from the MultiDark simulation to find that cluster halos with massive neighbors (with masses above about 10% of the cluster halos) within 2–4 × r_vir had quieter cluster assembly histories recently than on average, and were more active beyond z ≈ 1. These halos indeed did not accrete their close-by halos and thus did not empty their neighborhood. On the contrary, a low number of neighbors in the same distance range is linked to the opposite scenario, namely recently active and quieter in the past. This helps to reconcile our findings with the results from simulations by Diemer & Kravtsov (2014). The clusters showing massive infalling groups are about to enter a phase of substantial accretion, which will result in a contraction of the splashback radius. We nevertheless stress the challenges of capturing the full extent of halo accretion in observed clusters, which is composed of a continuous stream of halos of widely different masses. As shown in Diemer & Kravtsov (2014), Diemer et al. (2017) and in Figure 4, the ratio r_sp/r_200,m shows the strongest dependence with mass accretion rate. Additionally, at fixed r_sp, the ratio r_sp/r_200,m depends nontrivially on halo mass, accretion rate and redshift (Diemer et al. 2017). We note that our results do not show strong dependence on the redshift and halo mass as a result of the narrow interval of both for the cluster sample considered here. Hence, we can further ascribe the observed trend of the ratio r_sp/r_200,m to the accretion rate of clusters. We note that the cluster sample with and without groups are characterized by discrepant mean masses, with the former larger by a factor of ≈2.6. A similar mass separation is found when using the different proxies for cluster dynamical state. However, this is not a major source of systematic uncertainty in our results because the differences in the assembly rate of clusters over the mass and redshift range under consideration are minimal (Pizzardo et al. 2021). Following the Press and Schechter formalism (Lacey & Cole 1993), the typical mass of infalling halos is approximately 10% of the mass of the main accreting halo, implying that the typical infalling group mass is correlated with the cluster mass. Therefore, we argue that the detection of infalling groups via X-ray emission in less massive clusters will be impeded by survey limits more than in more massive clusters. This is coupled with the luminosity boost of bright X-ray cluster cores, which impacts favorably the detection of less massive clusters. Crucially, our study confirms that massive clusters are undergoing continuous accretion of group-like halos, which does not affect the presence of the splashback feature as drastically as for clusters classified via different dynamical proxies (e.g., central entropy). This could be related to the different timescales considered by the different proxies, and with respect to simulations (Diemer & Kravtsov 2014; Diemer et al. 2017), and requires further investigation. Overall, the whole cluster sample provides the less stringent constraints on the recovered splashback feature resulting in a ratio between the 3D splashback radius and r_200,m of r_sp/r_200,m = 1.74 ± 0.34. This supports the scenario in which the dynamical state of individual clusters dilutes the stacked signal of the splashback feature.

Zoom In Zoom Out Reset image size

5.3. Comparison with Previous Observations

5.4. Future Outlook

M.B., R.B., G.P.S., and S.L.M. acknowledge support from the Science and Technology Facilities Council through grant number ST/N021702/1. A.B. acknowledges support from NSERC (Canada). M.B. acknowledges Benedikt Diemer for the helpful discussions and support using the Colossus toolkit. ⁷ M.B. thanks Arya Farahi, August Evrard, and Surhud More for the insightful discussions. We warmly thank and acknowledge Maria Pereira and Eiichi Egami for their work on the Arizona Cluster Redshift Survey (ACReS), and for sharing their data with us. Computational work was performed using the University of Birmingham's BlueBEAR HPC service. We acknowledge the use of the Python package ChainConsumer (Hinton 2016) for the parameter corner plots.

We have explored how classifying clusters using different proxies of their dynamical state tracing the properties of the central regions, reverberates at greater radii on the splashback feature. These proxies include central entropy, X-ray surface brightness morphology, i.e., concentration and centroid-shift, K-band luminosity gap between the two most luminous galaxies, the offset between X-ray peak emission and BCG location. We follow Mulroy et al. (2019) in using the entropy measurements of Sanderson et al. (2009) to divide the clusters into those with stronger cooling based on K₀( < 20 kpc) < 80 keV cm⁻² and those with less strong cooling, i.e., larger values of K₀. The X-ray concentration parameter is defined as the ratio of surface brightness at two characteristic radii, the first encompassing the typical size of cool cores and the second the majority of X-ray emission ${c}_{\mathrm{sb}}=\tfrac{{S}_{{\rm{X}}}(\lt 40\,\mathrm{kpc})}{{S}_{{\rm{X}}}(\lt 400\,\mathrm{kpc})}$, to maximize the dichotomy between the surface brightness distribution of cool-core and non-cool-core clusters (Cassano et al. 2010). Here, we utilize the concentration measures from Mulroy et al. (2019), who extracted it from the X-ray surface brightness maps from Chandra/ACIS-I and XMM-Newton/EPIC observations. Similarly, the measure of the position of the X-ray emission centroid in circular apertures of increasing radii has been used to label cluster merger activity. In particular, high standard deviation (〈w[0.01 r_500,x]〉 > 1) of the centroid peak has been associated with dynamically disturbed clusters (Poole et al. 2006; Maughan et al. 2008; Mahdavi et al. 2013). Here we use the measures from Martino et al. (2014), cautioning about projection effects, which could hide ongoing line-of-sight mergers.

The properties of the BCG can help identify the formation history of the halo in which it is hosted. In particular, the older the halo, the higher the magnitude gap between the BCG and the second most luminous galaxies (Ponman et al. 1994; Gozaliasl et al. 2014; Farahi et al. 2020). This is due to the action of dynamical friction, which facilitates the fall of massive galaxies toward the deep end of the cluster potential where the BCG resides. Subsequent mergers result in a dominant central object surrounded by smaller galaxies. We select a threshold value of the K-band magnitude gap Δm_k,12 = 0.5 (Smith et al. 2010), dividing the bimodal distribution of the full sample.

Bridging between the information provided by the ICM emission and the BCG luminosity, the projected offset between the peak of the X-ray emission and the bulk of the stellar light from the BCG can reveal ongoing cluster dynamical activity (Bildfell et al. 2008). In particular, in the case of old and dynamically quiet clusters, the X-ray emission and the BCG should coincide indicating the deep core of the cluster potential well. Following Sanderson et al. (2009), we select a threshold value of 0.03 r_500,wl. Figure 5 summarizes the distributions of the dynamical proxies considered of the parent high-L_X LoCuSS sample, and of the cluster subsample used in this work.

Zoom In Zoom Out Reset image size

Figures 6 and 7 show the density profiles and slopes of clusters classified according to dynamical proxies as presented in Appendix A. Figure 8 shows corner plots of the parameters' posterior distribution from the fit of the density profiles of clusters with and without infalling groups.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size