Lensing Constraints on the Mass Profile Shape and the Splashback Radius of Galaxy Clusters

The lensing signal around galaxy clusters can, in principle, be used to test detailed predictions for their average mass profile from numerical simulations. However, the intrinsic shape of the profiles can be smeared out when a sample that spans a wide range of cluster masses is averaged in physical length units. This effect especially conceals rapid changes in gradient such as the steep drop associated with the splashback radius, a sharp edge corresponding to the outermost caustic in accreting halos. We optimize the extraction of such local features by scaling individual halo profiles to a number of spherical overdensity radii, and apply this method to 16 X-ray-selected high-mass clusters targeted in the Cluster Lensing And Supernova survey with Hubble. By forward-modeling the weak- and strong-lensing data presented by Umetsu et al., we show that, regardless of the scaling overdensity, the projected ensemble density profile is remarkably well described by an NFW or Einasto profile out to $R \sim 2.5h^{-1}$Mpc, beyond which the profiles flatten. We constrain the NFW concentration to $c_{200c} = 3.66 \pm 0.11$ at $M_{200c} \simeq 1.0 \times 10^{15}h^{-1}M_\odot$, consistent with and improved from previous work that used conventionally stacked lensing profiles, and in excellent agreement with theoretical expectations. Assuming the profile form of Diemer&Kravtsov and generic priors calibrated from numerical simulations, we place a lower limit on the splashback radius of the cluster halos, if it exists, of $R_{sp}/r_{200m}>0.89$ ($R_{sp}>1.83h^{-1}$Mpc) at 68% confidence. The corresponding density feature is most pronounced when the cluster profiles are scaled by $r_{200m}$, and smeared out when scaled to higher overdensities.


INTRODUCTION
Galaxy clusters, as the largest gravitationally bound objects formed in the universe, play a fundamental role in our understanding of cosmology and structure formation. A key ingredient for cluster-based cosmology is the distribution of dark matter (DM) in and around cluster halos. In this context, the standard Λ cold dark matter (ΛCDM) model and its variants, such as self-interacting DM (SIDM, Spergel & Steinhardt 2000) and wave DM (ψDM, Schive et al. 2014), provide distinct, observationally testable predictions.
For the case of collisionless DM, high-resolution N -body simulations exhibit an approximately "universal" form for the spherically averaged density profile of halos in gravitational quasi-equilibrium (Navarro et al. 1997, hereafter NFW), ρ(r) ∝ (r/r s ) −1 (1 + r/r s ) −2 , where r s is the characteristic scale radius at which the logarithmic density slope d ln ρ/d ln r equals −2. In this context, the halo concentration, c vir ≡ r vir /r s , is a key quantity that characterizes the structure of a halo (all relevant symbols are defined in detail at the end of this section). In the hierarchical ΛCDM picture of structure formation, concentration is predicted to correlate with halo mass, M vir , because r s stays nearly constant after an early phase of rapid accretion, whereas r vir continues to grow through a mixture of physical accretion and pseudoevolution (Bullock et al. 2001;Wechsler et al. 2002;Cuesta et al. 2008;Zhao et al. 2009;Diemer et al. 2013). As clus-keiichi@asiaa.sinica.edu.tw * Based in part on data collected at the Subaru Telescope, which is operated by the National Astronomical Society of Japan. 1 Institute of Astronomy and Astrophysics, Academia Sinica, P.O. Box 23-141, Taipei 10617, Taiwan 2 Institute for Theory and Computation, Harvard-Smithsonian Center for Astrophysics, 60 Garden St., Cambridge, MA 02138, USA ter halos are, on average, still actively growing today, they are expected to have relatively low concentrations, c vir ∼ 4-5 (Bhattacharya et al. 2013;Dutton & Macciò 2014;. These general trends are complicated by the large scatter in halo growth histories, which translates into a significant diversity in their density profiles (Ludlow et al. 2013).
Recently, closer examination of the outer halo density profiles in collisionless ΛCDM simulations has revealed systematic deviations from the universal NFW or Einasto (1965) form (Diemer & Kravtsov 2014, hereafter DK14). In particular, the profiles exhibit a sharp drop in density, a feature associated with the last shell that has reached the apocenter of its first orbit after accreting onto a halo in spherical collapse models (Gunn & Gott 1972;Fillmore & Goldreich 1984;Bertschinger 1985). The location of this "splashback radius," R 3D sp , is within a factor of two of r 200m and depends on the mass accretion rate of halos, with a secondary dependence on redshift (Diemer & Kravtsov 2014;Adhikari et al. 2014;More et al. 2015;Shi 2016;More et al. 2016;Adhikari et al. 2016;Mansfield et al. 2016). The splashback radius constitutes a physically motivated halo boundary because (at least in the spherical case) material outside R 3D sp is on its first infall into the halo, whereas material inside of it is orbiting in the halo potential. More et al. (2016) first observed the splashback feature in stacked galaxy surface density profiles around clusters (see also Tully 2015;Patej & Loeb 2016). Their measured R 3D sp is somewhat smaller than expected from the numerical calibration of More et al. (2015). This intriguing disagreement could be due to subtle effects in the analysis, errors in the numerical calculation, baryonic physics affecting cluster member galax-ies, or hitherto undetected properties of the DM itself, such as self-interaction. Given this wide range of possible reasons for the disagreement, other observational probes are of great interest. In particular, we are looking for a test that is subject to different systematic uncertainties than cluster member density profiles, but is still applicable to high-mass galaxy clusters where the splashback signal is expected to be strongest (Diemer & Kravtsov 2014;Adhikari et al. 2014). One potential probe that fulfills these requirements is gravitational lensing, because it measures the total mass profile rather than the distribution of subhalos, while the signal is strongest in galaxy clusters.
Cluster gravitational lensing offers a well-established method for testing halo structure, through observations of weak shear lensing (e.g., von der Applegate et al. 2014;Gruen et al. 2014;Umetsu et al. 2014;Hoekstra et al. 2015;Melchior et al. 2016), weak magnification lensing (e.g., Hildebrandt et al. 2011;Umetsu et al. 2011b;Coupon et al. 2013;Ford et al. 2014;Jimeno et al. 2015;Chiu et al. 2016;Ziparo et al. 2016), strong gravitational lensing (e.g., Broadhurst et al. 2005b;Zitrin et al. 2009;Coe et al. 2010;Jauzac et al. 2014;Diego et al. 2015), and the combination of all these effects (e.g., Broadhurst et al. 2005a;Umetsu & Broadhurst 2008;Umetsu et al. 2011aUmetsu et al. , 2012Umetsu et al. , 2015Umetsu et al. , 2016Coe et al. 2012;Medezinski et al. 2013). Over the last decade, cluster lensing observations (Broadhurst et al. 2005a;Okabe et al. 2010Okabe et al. , 2013Umetsu et al. 2011aUmetsu et al. , 2016Newman et al. 2013) have established that the projected total mass distribution within individual and stacked clusters is well described by a family of density profiles predicted for cuspy DM-dominated halos, such as the NFW (Navarro et al. 1997), Einasto (Einasto 1965), and DARKexp  models. Subsequent systematic studies targeting lensing-unbiased cluster samples (e.g., Okabe et al. 2013;Umetsu et al. 2014Umetsu et al. , 2016Merten et al. 2015;Du et al. 2015) show that the cluster lensing measurements are also in agreement with the theoretical c-M relation that is calibrated for recent ΛCDM cosmologies with a relatively high normalization (Bhattacharya et al. 2013;Dutton & Macciò 2014;Meneghetti et al. 2014;. In principle, any feature in the density profiles predicted by numerical simulations is directly accessible by lensing observations of a large sample of galaxy clusters (e.g., Okabe et al. 2010Okabe et al. , 2013Umetsu et al. 2014Umetsu et al. , 2016Miyatake et al. 2016). In reality, however, two effects make such measurements difficult. First, projecting the density profile into two dimensions smooths out features because any given sightline crosses a range of three-dimensional cluster radii. Second, in order to average lensing observations of individual clusters, we need to stack their density profiles using some radial scale. Conventionally, physical length units are chosen for this rescaling (e.g., Okabe et al. 2010;Sereno et al. 2015). If the cluster sample spans a wide range of masses, such a stacking procedure is likely to smooth out the intrinsic density profiles of the clusters, and sharp features such as the splashback radius in particular. Instead, we wish to rescale the profiles by a halo radius in units of which the features we are interested in are universal, i.e. they appear at the same rescaled radius. Numerical simulations show that this choice of scaling radius is far from trivial: while the inner profiles (r < ∼ r vir ) are most universal with halo radii that scale with the critical density of the universe (such as r 200c ), the outer profiles are most univer-sal when expressed in units of radii that scale with the mean cosmic density, such as r 200m (Diemer & Kravtsov 2014;Lau et al. 2015). These predictions have not hitherto been tested observationally.
In this paper, we develop new methods for scaling and modeling stacked cluster lensing profiles, and undertake the first investigation of the splashback radius based on lensing observations. We use the data presented in Umetsu et al. (2016, hereafter U16), who performed a joint analysis of stronglensing, weak-lensing shear and magnification data sets for 20 high-mass clusters targeted in the Cluster Lensing And Supernova survey with Hubble (CLASH, Postman et al. 2012;Umetsu et al. 2014Umetsu et al. , 2016Zitrin et al. 2015;Merten et al. 2015). Their analysis combines constraints from 16-band Hubble Space Telescope (HST) observations and wide-field multicolor imaging taken primarily with Suprime-Cam on the Subaru telescope. Such a joint analysis of multiple lensing probes allows us not only to improve the precision of mass reconstructions, but also to calibrate systematic errors inherent in each probe (Rozo & Schmidt 2010;Umetsu 2013). The large radial range covered by the combination of weak-and strong-lensing data allows us to explore a range of scaling overdensities, and to investigate their impact on the stacked ensemble fit. Thanks to the improved stacking procedure, we derive tighter constraints on halo concentration than in U16, and put a lower limit on the splashback radius of the stacked CLASH density profile.
The paper is organized as follows. In Section 2 we summarize the characteristics of the CLASH sample and describe the data used in this study. We then outline our procedure for modeling the cluster lensing profiles, and test its robustness using synthetic CLASH weak-lensing data. In Section 3 we apply our methodology to the CLASH lensing data and fit them with NFW, Einasto, and DK14 profiles. We discuss the results in Section 4. Finally, a summary is given in Section 5.
Throughout this paper, we adopt a spatially flat ΛCDM cosmology with Ω m = 0.27, Ω Λ = 0.73, and a Hubble constant H 0 = 100h km s −1 Mpc −1 with h = 0.7. We denote the mean matter density of the universe as ρ m and the critical density as ρ c . We use the standard notation M ∆c or M ∆m to denote the mass enclosed within a sphere of radius r ∆c or r ∆m , within which the mean overdensity equals ∆ c × ρ c (z) or ∆ m × ρ m (z) at a particular redshift z, such that M ∆c = (4π/3)∆ c ρ c (z)r 3 ∆c and M ∆m = (4π/3)∆ m ρ m (z)r 3 ∆m . We generally denote three-dimensional cluster radii as r, and reserve the symbol R for projected clustercentric distances. We define the splashback radius of a three-dimensional density profile, R 3D sp , as the radius where the logarithmic slope of the profile is steepest. Similarly, we use R 2D sp to denote the splashback radius derived from the steepest slope of the projected profile. We compute the virial mass and radius, M vir and r vir , using an expression for ∆ vir (z) based on the spherical collapse model (Appendix A of Kitayama et al. 1998). For a given overdensity ∆, the concentration parameter is defined as c ∆ = r ∆ /r s . All quoted errors are 1σ confidence limits (CL) unless otherwise stated.

CLASH Sample and Data
The CLASH survey (Postman et al. 2012) is a 524-orbit HST Multi-Cycle Treasury program targeting 25 high-mass galaxy clusters. Of these, 20 CLASH clusters were selected to be X-ray hot (T X > 5 keV) and to have a high degree of regularity in their X-ray morphology, with no lensing information used a priori. Another subset of five clusters were selected by their high-magnification properties. These highmagnification clusters often turn out to be complex, massive merging systems (e.g., Zitrin et al. 2013;Medezinski et al. 2013). A complete definition of the CLASH sample is given in Postman et al. (2012).
In this work, we shall focus on the analysis of the Xray-selected subsample to simplify the interpretation of our results. Numerical simulations suggest that this subsample is mostly composed of relaxed clusters (∼ 70%) and largely free of orientation bias (Meneghetti et al. 2014). Specifically, we use a lensing-unbiased subset of 16 CLASH X-ray-selected clusters taken from U16, who performed a comprehensive analysis of the strong-lensing, weak-lensing shear and magnification data. Our cluster sample lies in the redshift range 0.19 < ∼ z l < ∼ 0.69 and over a mass range 5 < ∼ M vir /(10 14 h −1 M ) < ∼ 30 (Table 1), spanning a factor of ∼ 6 in halo mass M ∆ , or a factor of ∼ 1.8 in r ∆ ∝ M 1/3 ∆ . A full description of the data used in this work is given by U16. Here, we provide only a brief summary of the most relevant aspects of the lensing reconstructions.
The U16 analysis uses the cluster lensing mass inversion (CLUMI) code developed by Umetsu et al. (2011b) and Umetsu (2013), in which lensing constraints are combined a posteriori in the form of azimuthally averaged radial profiles. U16 used constraints spanning the radial range 10 -960 obtained from 16-band HST observations (Zitrin et al. 2015) and wide-field multicolor imaging (Umetsu et al. 2014) taken primarily with Subaru/Suprime-Cam. The position of the brightest cluster galaxy (BCG) is adopted as the cluster center ( Table 1 in U16). Umetsu et al. (2014) obtained weaklensing shear and magnification measurements in 10 logarithmically spaced radial bins (N WL = 10) over the range 0.9 -16 for all clusters observed with Subaru, and 0.9 -14 for RX J2248.7−4431 observed with ESO/WFI. Zitrin et al. (2015) constructed detailed mass models for each cluster core from a joint analysis of HST strong and weak-shear lensing data. U16 constructed enclosed projected mass constraints for a set of four equally-spaced integration radii (10 -40 , N SL = 4) from the HST lensing analysis of Zitrin et al. (2015). U16 combined these full lensing constraints for individual clusters in their joint likelihood analysis to reconstruct binned surface mass density profiles with N bin = N SL + N WL + 1. We have N bin = 15 bins for all clusters, except N bin = 11 for RX J1532.9+3021 with N SL = 0, for which no secure identification of multiple images was made (Zitrin et al. 2015). The Σ profiles used in this work are shown in Figure 11 of U16.
U16 accounted for various sources of errors. Their analysis includes four terms in the total covariance matrix C ij of the Σ profile, where C stat represents statistical observational errors, C sys contains systematic uncertainties due primarily to the residual mass-sheet uncertainty (Umetsu et al. 2014), C lss is the cosmic-noise covariance matrix due to projected uncorrelated large-scale structures (Hoekstra 2003;Umetsu et al. 2011a), and C int accounts for the intrinsic variations of the projected cluster lensing signal at fixed mass due to variations in halo concentration, cluster asphericity, and the presence of correlated halos (Gruen et al. 2015). Overall, the reconstruction uncertainty is dominated by the C stat term (Figure 1 of U16). The relative contribution from the C int term becomes increasingly dominant at small cluster radii, especially at θ < ∼ 2 . The impact of the C lss term is most important at large cluster radii, where the cluster signal is small. Table 1 summarizes CLASH lensing determinations of the mass and concentration parameters (M ∆ , c ∆ ) for our 16 clusters based on the full lensing analysis of U16. These values were obtained from spherical NFW fits to individual cluster Σ profiles using the total covariance matrix C (Equation (1)), restricting the fitting range to R ≤ 2h −1 Mpc (∼ 2r 500c ∼ r 200m for the CLASH sample) to avoid systematic effects (Becker & Kravtsov 2011;Meneghetti et al. 2014). In Table 1, we also list effective overdensity masses M eff ∆ of the sample, which were obtained by U16 from a spherical NFW fit to the stacked surface mass density profile of the 16 CLASH clusters. The stacked ensemble has an effective halo mass M eff vir = (11.99 ± 0.93) × 10 14 h −1 M and an effective halo concentration c eff vir = 4.69 ± 0.35 (c eff 200c = 3.76 ± 0.28; see Table 1), and lies at a sensitivity-weighted average redshift of z eff l = 0.337, close to the median redshift of z l = 0.352. U16 quantified potential sources of systematic uncertainty in their mass calibration (see their Section 7.1), such as the effect of dilution of the weak-lensing signal by cluster members (2.4%), photometric-redshift bias (0.27%), shear calibration uncertainty (5%), and projection effects of prolate halos (3%). Combining them in quadrature, the total systematic uncertainty in the absolute mass calibration was estimated to be 6%. This is in close agreement with the value ∼ 8% empirically estimated from the shear-magnification consistency test of Umetsu et al. (2014).
Another potential source of systematic errors is smoothing of the central lensing signal from miscentering effects (Johnston et al. 2007;Umetsu et al. 2011a;Du & Fan 2014). On average, the sample exhibits a small positional offset between the BCG and the X-ray peak, characterized by an rms offset of σ off 11h −1 kpc (Umetsu et al. 2014(Umetsu et al. , 2016, which is much smaller than the typical resolution limit of our HST lensing data (θ min = 10 , corresponding to 35h −1 kpc at z l = 0.35). Hence, the miscentering effects are not expected to significantly affect our ensemble lensing analysis.

Radial Scaling of the Profiles
One of our main goals in this work is to investigate how the radial scaling of stacked profiles influences the fit results. Ideally, one would like to stack profiles as a function of the exact halo radius where a particular feature is expected, e.g. the NFW scale radius or the splashback radius. However, since their locations are a priori unknown, we need to resort to an alternative halo radius that can be measured for individual clusters, generally a spherical overdensity radius r ∆ . Now, the goal is to choose a definition in which the location of the feature in question is universal, i.e. independent of halo mass, and possibly of redshift.
Unfortunately, there is no guarantee that any one definition will be ideal for multiple features. In fact, DK14 discovered that halo density profiles in N -body simulations prefer different scaling radii in different regions of the profile: while the inner profiles (and thus concentrations) are most universal when expressed in units of halo radii that scale with the critical density ρ c (z) of the universe, such as r 200c , the outer profiles (and thus the splashback radius) are most universal in units of halo radii that scale with the mean cosmic density ρ m (z), such as r 200m . These scalings were confirmed in hydrodynamical simulations of galaxy clusters (Lau et al. 2015, see also Shi 2016). As we wish to investigate both the inner and outer profiles, we repeat our analysis using a number of different scalings covering a range between 2500ρ c and 200ρ m 94ρ c at z eff l = 0.337. To construct the scaled surface mass density profiles for the CLASH sample, we use our full lensing constraints on the NFW parameters M ∆ and c ∆ of each individual cluster as obtained by U16 (Section 2.1; see Table 1), and normalize their observed Σ(R) profiles to a given overdensity ∆ of interest. We stress that we do not rely on scaling relations (e.g., the c-M relation), but use observational lensing constraints on r ∆ and Σ(r ∆ ) for each cluster. The high-quality, multi-scale weak-and strong-lensing constraints from the CLASH survey enable us to explore the wide range of overdensities listed above.
We choose the NFW model for the scaling because recent cluster lensing observations show that the "projected total" matter distribution in the intracluster region (R < ∼ r 200m ) is in excellent agreement with the NFW form (Umetsu et al. 2011a(Umetsu et al. , 2014(Umetsu et al. , 2016Beraldo e Silva et al. 2013;Newman et al. 2013;Okabe et al. 2013;Niikura et al. 2015;, as predicted for collisionless DM-dominated halos in N -body cosmological simulations Meneghetti et al. 2014). We demonstrate the consistency of this choice with our results in Section 3.

Fitting Functions
The second goal of this work (besides investigating radial scalings) is to constrain the splashback radius of the CLASH cluster sample. While NFW and Einasto fitting functions were sufficient for the radial rescaling, they describe only the 1-halo term and do not take the steepening due to the splashback radius into account. Thus, we also fit the cluster lensing profiles with the more flexible fitting function of DK14 which was calibrated to a suite of ΛCDM N -body simulations. We emphasize that the quality of the data is insufficient to distinguish among the different profile models, which all describe the data very well (see Section 3.1). Nevertheless, in order to constrain the location of the splashback radius, we assume that the ΛCDM simulations of DK14 describe real cluster halos and use the DK14 profile as a fitting function in conjunction with generic priors.
Furthermore, we note that the "true" location of the splashback radius is not strictly equivalent to a particular location in the spherically averaged density profile. However, we follow More et al. (2015) in defining R 3D sp as the radius where the logarithmic slope of the three-dimensional density profile is steepest. According to this definition, R 3D sp is expected to lie within a factor of two of r 200m (More et al. 2015). Furthermore, the steepest slope would need to be steeper than that expected from the sum of the Einasto profile and the 2-halo term at large scales if a detection were to be claimed .
The DK14 fitting formula is described by eight parameters, and is sufficiently flexible to reproduce a range of fitting functions for the DM density profile, such as the halo model Hikage et al. 2013). We use the publicly available code COLOSSUS (Diemer 2015) for many of calculations relating to density profiles. The DK14 model is given by with r piv = 5r 200m and ∆ max = 10 3 . Here ∆ max has been introduced as a maximum cutoff density of the outer term to avoid a spurious contribution at small halo radii (Diemer 2015). The Einasto profile ρ Einasto describes the intracluster mass distribution, with the shape parameter α describing the degree of profile curvature and r s the scale radius at which the logarithmic slope is -2. The transition term f trans characterizes the steepening around a truncation radius, r t . The outer term ρ outer , given by a softened power law, is responsible for the correlated matter distribution around clusters, also known as the 2-halo term, ρ 2h . DK14 found that this fitting function provides a precise description ( < ∼ 5%) of their simulated DM density profiles at r < ∼ 9r vir . At larger radii, the outer term is expected to follow a shape proportional to the matter correlation function (e.g., Oguri & Takada 2011).
To scale out the mass dependence of the density profile, we use an arbitrary overdensity radius r ∆ as a pivot radius, and define a scaled version of ∆ρ(r) as a function of r = r ∆ x, namely This scaled DK14 model is described by a set of seven dimensionless parameters, namely the halo concentration c ∆ = r ∆ /r s , the Einasto shape parameter α, the dimensionless truncation radius τ ∆ = r t /r ∆ , the relative normalization of the outer term B ∆ ∝ b e , and three additional shape parameters (s e , β, γ) for the transition and outer terms. The model reduces to the Einasto model (specified by c ∆ and α) when f trans = 1 (τ ∆ → ∞) and f outer = 0 (B ∆ = 0). We refer the reader to Appendix A for a more detailed description of the scaled DK14 profile function. Similarly, we define a scaled version of the NFW profile, By projecting ∆ρ(r) along the line of sight, we derive the scaled surface mass density, which is a lensing observable, normalized as y ∆ (x) = 1 at x = 1, where The projected y ∆ (x) profile is modeled in terms of the scaled density profile f ∆ as (7) We note that it is straightforward to generalize our approach to shear-only weak-lensing observations where the differential surface mass density ∆Σ(R) = Σ(< R) − Σ(R) is a direct observable in the weak-lensing limit (e.g., Bartelmann & Schneider 2001).

Parameter Inference
We use a Bayesian approach to infer the shape and structural parameters of the mass distribution of the CLASH sample. We restrict the optimization to (c ∆ , α, τ ∆ , B ∆ ), the primary model parameters that describe the scaled DK14 model, and marginalize over the remaining parameters, (s e , β, γ), using priors based on the N -body simulations of DK14. In particular, following More et al. (2016), we adopt Gaussian priors of log 10 β = log 10 6±0.2 and log 10 γ = log 10 4±0.2 5 , allowing a wide, representative range of values. We assume a Gaussian prior on s e of 1.5 ± 0.1, centered around the value found by DK14. We use unconstraining flat priors on (c ∆ , α, B ∆ ). For τ ∆ , we assume τ 200m ∈ [0, 5], where the upper bound corresponds approximately to r t = 5r 200m 10h −1 Mpc for the CLASH sample (Table 1), which is larger than the maximum data radius of ∼ 5h −1 Mpc. We translate this prior to a given overdensity ∆ using the effective r ∆ radius of the sample (Section 2.1 and Table 1 We use a Gaussian log-likelihood − ln L(p) = χ 2 (p)/2 + const. with a χ 2 function given by where , C ∆ is the total covariance matrix of the scaled data y ∆ , and y ∆ (x i |p) represents the model for y ∆,i with parameters p. Combining all 16 clusters in our sample, we have a total of 236 data points (Section 2.1).
We use a Markov chain Monte Carlo (MCMC) approach with Metropolis-Hastings sampling to sample from the posterior distribution of the parameters, (c ∆ , α, τ ∆ , B ∆ , s e , log 10 β, log 10 γ), given the data and the priors stated above. We largely follow the sampling procedure of Dunkley et al. (2005) but employ the Gelman-Rubin statistic (Gelman & Rubin 1992) as a convergence criterion of the generated chains. Once convergence to a stationary distribution is achieved, we run a long, final chain of 10 5 sampled points, which is used for our parameter estimation and error analysis. We also use the MINUIT package from the CERN program libraries to find the global maximum a posteriori estimate of the joint probability distribution. This procedure allows a further refinement of the best-fit solution with respect to the one obtained from the MCMC sampling (see discussions in Planck Collaboration et al. 2014; Penna-Lima et al. 2016).

Tests with Synthetic Data
Given the simulation results of DK14, we expect the signature of the splashback radius in lensing data to be weak. Thus, one concern is that fitting with the DK14 profile function might introduce a spurious "detection" of a splashback feature due to systematics or overfitting in the presence of noise. In order to address this potential issue, we have tested our procedure (including data analysis, mass reconstruction, stacking, and the fitting process) on simulated lensing data. We focus  . Test of our forward-modeling method. We create 50 realizations of CLASH-like weak-lensing data by computing synthetic shear and magnification catalogs from analytically modeled cluster lenses. We simulate two configurations of the outer density profile, one without a splashback feature (top panels) and one with a splashback-like feature (bottom panels, see Appendix B for details on the lens models). The left panels show the scaled surface mass density Σ(R)/Σ(r 200m ) of the synthetic observations, where the black solid line represents the noise-free, sensitivity-weighted profile of the 16 clusters in the synthetic sample. The blue shaded region indicates the 1σ bounds on the ensemble of DK14 fits to the synthetic data. The gray lines show scaled Σ(R) profiles of individual clusters reconstructed from each particular source realization. The right panels show the logarithmic slope of the three-dimensional DK14 profiles, d ln ∆ρ(r)/d ln r, as a function of r/r 200m . For the lens model with a splashback-like feature, the range from the 16th to the 84th percentile of R 3D sp (gray vertical shaded area) inferred from the synthetic data is consistent with the sensitivity-weighted expectation value of the sample (red vertical line), R 3D sp /r 200m 0.87. In each panel, the blue dotted lines indicate the errors scaled to match the overall CLASH weak-lensing sensitivity. This test demonstrates the robustness of our fitting method, and in particular that fitting the DK14 profile does not introduce a spurious splashback feature. See Section 2.5 and Appendix B for details. on the recovery of the lensing signal in the noisy outer regions around r 200m where the splashback feature is expected. Hence, we consider only the wide-field weak-lensing observables, namely the shear and magnification effects in the subcritical regime. To this end, we create 50 source realizations of synthetic shear and magnification catalogs for our 16 CLASH clusters, each modeled as an NFW halo specified by its redshift z l (Table 1) and (M 200c , c 200c ) parameters fixed to the observed central values ( Table 2 of U16). For each NFW cluster, we consider two configurations of the outer density profile, one with and one without a splashback-like feature. For technical reasons, we do not use a DK14 profile for the synthetic cluster lenses and instead substitute a profile that in-troduces a similar density drop (see Appendix B for details).
The results of applying our methods to the synthetic weaklensing data are summarized in Figure 1. The lower and upper panels correspond to the simulations with and without a splashback-like feature, respectively. The blue shaded region in the left panels shows the mean and standard deviation of the best-fit DK14 profiles inferred from 50 realizations of the synthetic data for each configuration. On average, this fit is in excellent agreement with the noisefree, sensitivity-weighted averaged input profile (black solid curve). In the right panels, we show the corresponding logarithmic density slope d ln ∆ρ(r)/d ln r. For the model with a splashback-like feature, the range from the 16th to the 84th Table 2 Best-fit NFW, Einasto, and DK14 parameters for the CLASH Sample Note.
-The global best-fit parameters and their 1σ fractional errors (in per cent) for each model, with five different values of the scaling overdensity ∆. For the DK14 parameter B∆, we give the best-fit parameter and its 1σ uncertainty. For the DK14 parameter τ∆, we provide lower limits (68% CL) and best-fit parameter values in parentheses.
percentile of R 3D sp (0.74 ≤ R 3D sp /r 200m ≤ 1.27; vertical shaded area) inferred from the synthetic data is consistent with the sensitivity-weighted expectation value of the sample, R 3D sp /r 200m 0.87. In Figure 1, we also indicate the errors scaled to the CLASH weak-lensing sensitivity (dotted lines). The errors are computed by matching the synthetic (NFW) to the observed (CLASH) total signal-to-noise ratio (S/N) of the Σ profiles, where (S/N) 2 = N halo n=1 Σ t C −1 Σ n . This comparison suggests that the overall uncertainty in our synthetic observations is underestimated by 36%. 6 Nevertheless, the 1σ uncertainty after this correction is small compared to the absolute value of the slope at r < ∼ 2r 200m . Even though the determination becomes noisy at radii around and beyond r 200m , the steepening relative to the NFW profile is marginally identified at the 1.7σ level at the expected location, r/r 200m ∼ 0.9.
This test demonstrates that our analysis methods are able to accurately reproduce the input sensitivity-weighted density profile and its logarithmic gradient even at realistic noise levels. The results also show that the priors adopted (Section 2.4) are generic and flexible enough to reproduce the NFW-like shape of the profile, as well as a splashback feature. Importantly, we note that our analysis pipeline does not introduce spurious gradients that mimic the characteristic splashback feature, namely a steepening followed by an upturn due to the contribution from the 2-halo term.

RESULTS
Our main results are the best-fit NFW, Einasto, and DK14 profiles resulting from a simultaneous fit to the scaled Σ profiles of 16 CLASH clusters (Table 2 and Figures 2, 3, and 4). Table 2 lists the best-fit parameters, their uncertainties, and the χ 2 values for each of five pivot overdensities with which the analysis was performed (∆ = 2500c, 500c, 200c, virial, and 200m).
In Figure 2, we show the projected density profile y ∆ (x) = Σ(r ∆ x)/Σ(r ∆ ) of the individual CLASH clusters (gray lines), as well as the best-fit DK14 (blue), NFW (dashed black), and Einasto (solid black) fits. The results are shown for two different choices of the pivot overdensity, namely 200m (left panel) and 2500c (right panel). In the upper panels of Figure 3, we show the corresponding threedimensional logarithmic slope d ln ∆ρ(r)/d ln r of the DK14 6 The underestimation is partly due to the assumed ellipticity dispersion of source galaxies, σg = 0.3 (Appendix B). This is 29% lower than the observed value, σg 042, which includes contributions fro both intrinsic shape and measurement noise. The rest (∼ 20%) can be accounted for by the intrinsic clustering and other error contributions in the weak-lensing magnification measurements (Umetsu et al. 2014(Umetsu et al. , 2016 as well as by the cosmic noise contribution, C lss . fit with 1σ errors (shaded area), as well as the NFW and Einasto slopes. Similarly, the lower panels show the logarithmic slope d ln Σ(R)/d ln R of the best-fit surface mass density profiles. Finally, Figure 4 shows the one-and twodimensional marginalized posterior distributions for the complete set of scaled DK14 parameters with ∆ = 200m, p = {c 200m , α, τ 200m , B 200m , s e , β, γ}. For all parameters, the global best-fit values coincide well with their respective peak values of the one-dimensional marginalized posterior distributions.
In the following sections, we discuss the fit quality as well as the inferred parameters for the inner (r < ∼ r vir ) and outer regions of the profile.

Fit Quality
The ensemble mass profile in projection is remarkably well described by an NFW or Einasto profile out to R ∼ 1.2r 200m or R ∼ 4.5r 2500c (Figures 2 and 3), beyond which the data exhibit a flattening that is not modeled by those fitting functions. Note that to calculate r ∆ and Σ(r ∆ ) for each individual cluster (Section 2.1), we employed the spherical NFW fits of U16 obtained with a restricted fitting range of R ≤ 2h −1 Mpc ∼ r 200m . The results shown here thus ensure the self-consistency of our analysis.
As the outer profiles are expected to be most universal with respect to ∆ = 200m (DK14), that definition is of particular relevance for the splashback radius. We thus use the DK14 model with ∆ = 200m as a baseline model. This model has the best-fit χ 2 of 180 for 232 degrees of freedom, corresponding to a probability of 99.5% to exceed the observed χ 2 value, assuming the standard χ 2 probability distribution function. The model is therefore in good agreement with the data. However, as we will discuss in Section 4.2, the improvement in the fit is not significant compared to the NFW or Einasto fit, implying that the parameters that describe the transition region and outer terms are not well constrained by the data. This is not surprising, because Figure 2 shows that the CLASH lensing data do not resolve the profile curvature in the transition region particularly well. Hence, the shape of the gradient feature at r ∼ r 200m , which locally deviates from the three-dimensional NFW and Einasto profiles (Figure 3), is specific to the assumed DK14 profile form.
We note that the χ 2 values in Table 2 decrease with increasing overdensity ∆, independent of the fitting function. The reason for this trend is that the inner Σ profiles are more tightly constrained by the data, especially from the HST lensing analysis, so that scaling the Σ profiles to higher overdensities reduces the overall scatter, which is dominated by the inner regions. We also note that the reduced χ 2 values in Table  2 are systematically smaller than unity, which may indicate that the number of degrees of freedom is overestimated ow-     Figure 2, but for the gradient of the profiles. Upper panels: logarithmic gradient of the inferred three-dimensional density profile as a function of the scaled cluster radius r/r ∆ . As in Figure 2, the results are shown for ∆ = 200m (left) and 2500c (right). The blue solid line and the blue shaded area represent the best-fit DK14 model and its 1σ uncertainty, and are compared to the NFW (black dashed) and Einasto (black solid) fits. The gray vertical shaded area indicates the range from the 16th to the 84th percentile of the marginalized posterior distribution of the splashback radius, R 3D sp /r ∆ (see Figure 5). Lower panels: same as the upper panels, but showing the logarithmic slope of the surface mass density profiles. The best-fit DK14 profile is shown as blue dots at the locations of the data points. For comparison, the slope of the conventionally stacked Σ profile (U16) is shown in rescaled units (green circles with error bars).
ing to the effects of nonlinear modeling (Andrae et al. 2010) and/or that the errors are conservatively overestimated. Since U16 found the reduced χ 2 for their fits to the same input data to be > ∼ 1 (their Table 4), it is unlikely that the errors are sig-nificantly overestimated. The best-fit values of concentration shown in Table 2 agree well between the different fitting models. This similarity is also apparent in Figures 2 and 3, as the models have very similar shapes for both scaling overdensities shown. Furthermore, the best-fit values for the shape parameter α agree between the Einasto and DK14 fits, and lie in the range 0.18 < ∼ α < ∼ 0.22 with a typical 1σ uncertainty of 0.06 for the DK14 model, and 0.17 < ∼ α < ∼ 0.21 for the Einasto model. An Einasto density profile with α ∼ 0.2 closely resembles an NFW profile (e.g., Ludlow et al. 2013).
The uncertainties on c and α allow us to assess the impact of the scaling overdensity. If the profiles are most universal as a function of a particular radius definition, we expect the fractional uncertainty on the fit to be smallest in that definition. DK14 investigated the universality of halo density profiles, and found that the inner profiles are most universal in units of r 200c . However, this statement refers primarily to the redshift scaling of different definitions, which we cannot test here owing to the limited redshift range of the CLASH sample. At fixed redshift and fixed mass, we expect any overdensity within a range around r 200c to lead to reasonably uni-versal inner profiles, whereas for very extreme definitions the scatter in the profiles might increase.
These expectations are borne out in our results. The fractional uncertainties on NFW-c and Einasto-α (the primary parameter that determines the shape of the profile for each model) are smallest for the 200m and virial scalings, and increase toward the highest overdensities (despite a lower χ 2 ). However, the uncertainties on Einasto-c are slightly lower at somewhat higher overdensities such as ∆ = 500c. Overall, it appears that a rescaling with densities around ∆ = 200c leads to relatively low uncertainties. Regardless of the profile model, scaling with ∆ = 2500c results in significantly higher uncertainties. Since the determination of spherical overdensity radii (or masses) is not less certain at high overdensities (U16, their Section 4.2), we conclude that the increased uncertainty arises because the inner profiles are less universal at very high overdensities. Furthermore, we note that the relative insensitivity of the inner profiles to ∆ is likely in part due to the sample selection based on X-ray regularity (Section 2.1), which is understood to significantly reduce the scatter in concentration (Meneghetti et al. 2014). The results have been confirmed by the CLASH full lensing analysis of U16 (see their Section 6.2, as well as Merten et al. 2015).
Most importantly, we find that scaling the profiles by any overdensity radius (except for ∆ = 2500c) improves the constraints on the best-fit parameters. As a result, the limits on concentration are tighter than those of U16 who scaled the profiles in physical units and found uncertainties of 8% on NFW-c, 11% on Einasto-c, and 18% on Einasto-α (see their  Table 4).

The Outer Profile and Splashback Radius
We now turn toward the outer profiles and particularly the inferred gradient profiles shown in Figure 3. A comparison of the three-and two-dimensional slopes highlights why detecting the splashback radius in surface density profiles is challenging in practice: even though there is a noticeable steepening in the three-dimensional slope, the two-dimensional slope drops very little, owing to projection effects (see Section 4.4).
Nevertheless, we can derive a splashback radius and mass from each of the MCMC samples of the DK14 parameters shown in Figure 4 (Equations A6 and A7). In Figure 5, we show the corresponding one-dimensional marginalized posterior distributions for R 3D sp and M sp ≡ M (< R 3D sp ), both in scaled units. The results are shown for two different values of the chosen pivot overdensity, ∆ = 200m (left panels) and 2500c (right panels). In each panel, the vertical dashed lines indicate the 16th, 50th, and 84th percentiles of the marginalized posterior distribution. The locations of R 3D sp /r ∆ and M sp /M ∆ derived from the best-fit model ( Table 2) are indicated by vertical solid lines in Figure 5, and are in close Note. -Lower limits (68% CL) and best-fit model values (in parentheses) for the three-dimensional splashback radius and mass. The splashback radius and mass in physical length units were converted using the effective overdensity radius r eff ∆ and mass M eff ∆ of the sample, respectively. agreement with the respective median values of the distributions. The posterior distributions show a tail extending toward large positive values of R 3D sp /r ∆ and M sp /M ∆ , associated with large values of the truncation parameter, τ ∆ = r t /r ∆ (Figure 4). A large τ ∆ indicates a profile without a welldefined steepening feature. We thus place uninformative upper bounds on R 3D sp /r ∆ and M sp /M ∆ . On the other hand, we obtain tighter lower bounds on these parameters because the inner Σ profiles of the clusters are well constrained by the combination of strong-lensing, weak-lensing shear and magnification data (U16). Table 3 summarizes the 68% confidence lower limits and best-fit model values (in parentheses) on the splashback radius and mass. We also list R 3D sp and M sp in physical length units converted with the effective overdensity radius (r eff ∆ ) and mass (M eff ∆ ) of the sample (Section 2.1). Using the fiducial scaling overdensity of ∆ = 200m, these lower limits are R 3D sp /r 200m > 0.89 and M sp /M 200m > 0.93, corresponding to R 3D sp > 1.83h −1 Mpc and M sp > 1.21 × 10 15 h −1 M . We note that the location of the steepest slope in projection, R 2D sp /r ∆ , is smaller than R 3D sp /r ∆ (Figure 3) because of projection effects ). Using our best-fit base DK14 model (∆ = 200m in Table 2), we find R 2D sp /R 3D sp 0.8. We discuss the difference between the two-and threedimensional splashback radii further in Section 4.4. From a comparison of the lower bounds on the splashback radius and mass in Table 3, we see that the steepening feature is most pronounced when the cluster profiles are scaled by r 200m , and is smeared out when scaled to higher overdensities, resulting in less strict lower limits (see also the right panels of Figure 3). This trend is consistent with the prediction of DK14 that halos reveal self-similar behavior in their outskirts when their profiles are expressed in units of spherical overdensity radii defined with respect to the mean density of the universe, especially r 200m .

DISCUSSION
In this section, we compare our results for the shape of the profile, concentration, and splashback radius with predictions from N -body simulations. We discuss observational and simulation effects that could potentially complicate our analysis.

The Impact of Priors
We imposed a number of priors on the parameters of the DK14 profile, namely on the slopes β and γ, on the steepening radius τ ∆ , and on the outer profile slope s e (see Figure 4). These priors were based on the results of DK14 for the median profiles of halo samples spanning a wide range of masses and mass accretion rates, and chosen conservatively, i.e. allowing a much larger range of parameter values than found for high-mass cluster halos in DK14. We find that imposing constraining theoretical priors leads to an inflated sensitivity to the splashback feature. As stated in Section 2.3, we use the DK14 profile as a flexible fitting function to determine the location of the steepest slope, R 3D sp , which is an observable quantity. In this context, p = {c ∆ , α, τ ∆ , B ∆ , s e , β, γ} are considered to be merely fitting parameters, and we allow them to take on values not expected from simulations, such as very large α.
Nevertheless, one might worry that our inferences regarding R 3D sp are informed by the priors because they clearly inform the posterior of some parameters (Figure 4). In particular, the asymptotic slope of the 1-halo term, γ, is essentially unconstrained by the fit and relatively steep due to the prior (for example, γ = 0 is effectively excluded). However, it is important to note that γ (as well as β) has no impact on the DK14 profile if τ ∆ is large, because the steepening then moves out of the observed region of the profile. Thus, the most critical prior is that on τ ∆ .
Our prior allows values of τ ∆ up to 5, placing the steepening at r t < ∼ 10h −1 Mpc which lies far outside the maximum radius of our data (R 5h −1 Mpc). Thus, profiles with τ 200m > ∼ 2.5 effectively reduce to an Einasto profile with a 2halo term (f trans ≈ 1) in the observed radial range. We have already demonstrated that a fit with these priors can reproduce a profile without steepening in the analysis of synthetic weaklensing data (Figure 1). However, we note that, in combination with a higher value of α and a positive value of B ∆ , even profiles with τ 200m > 2.5 can reproduce a splashback feature in the observed radial range (with an inner profile steeper than the best-fit NFW/Einasto profile). We have confirmed that the resulting profiles are, in fact, a reasonable description of the data. Our priors also allow a profile with a negative 2halo normalization, B ∆ ≤ 0 (i.e., underdense regions), which can produce a steepening gradient without an upturn feature. Here, the ensemble fits yield positive ∆ρ (i.e., ρ > ρ m ) in the observed range because Σ is constrained to be positive (U16).
In order to confirm that the τ ∆ prior does not significantly affect our results, we have performed a fit with a flat prior of τ 200m < 20, corresponding to r t < ∼ 40h −1 Mpc for the CLASH sample. With this relaxed τ ∆ prior, we find the same best-fit solution and a 68% CL lower bound of R 3D sp /r 200m > 0.90, compared to > 0.89 obtained with τ 200m < 5. This difference represents a 1% change, which is not significant given the current sensitivity. Moreover, we find no noticeable changes in the constraints on the density and gradient profiles shown in Figures 2 and 3. Therefore, we conclude that our results are sufficiently robust against the choice of the τ ∆ prior.

Which Density Profile Does the CLASH Sample Prefer?
A robust outcome of our analysis is that, regardless of the pivot overdensity ∆ chosen, the ensemble CLASH mass profile in projection is in full agreement with the NFW or Einasto profile out to ∼ r 200m , consistent with previous lensing results (e.g., Umetsu et al. 2011aUmetsu et al. , 2014Umetsu et al. , 2016Newman et al. 2013;Okabe et al. 2013;Niikura et al. 2015;). This also ensures the self-consistency of our analysis, in which we used the NFW fits of U16 to calculate r ∆ and Σ(r ∆ ) for individual clusters, where the fitting range was restricted to R ≤ 2h −1 Mpc ∼ r 200m . Our base DK14 model with ∆ = 200m gives a slightly better fit than the corresponding NFW or Einasto model in terms of the best-fit χ 2 values ( Table 2). The relative improvements are ∆χ 2 = χ 2 NFW − χ 2 DK14 1 and ∆χ 2 = χ 2 Einasto −χ 2 DK14 1 for three and two additional free parameters, respectively. Hence, the improvement is not statistically significant, implying that our inference of the outskirt feature depends on the choice of the fitting function and priors.
As discussed in Section 2.3 and by More et al. (2016), we would apply two requirements to claim a detection of the splashback radius with a DK14 model, namely (1) that the location of the steepest slope in three dimensions with respect to r ∆ can be identified at high significance and (2) that this steepening is greater than that expected from a DK14 model with f trans = 1 (i.e., τ ∆ 1, reducing to an Einasto profile). The second criterion is important to ensure that the steepening is actually associated with a density caustic rather than the transition to the 2-halo term.
Given these criteria, we do not find sufficient evidence for the existence of a splashback feature in the CLASH lensing data because the data do not have the sensitivity necessary to resolve the profile curvature in the transition region. This result is not surprising, as demonstrated in our simulated experiment (Section 2.5). On the other hand, assuming the DK14 profile form and generic priors calibrated with numerical simulations, we have placed lower limits on the splashback radius of the CLASH clusters (Table 3), if it exists. Since we cannot rule out models with f trans = 1, it is possible that the observed gradient feature in the outskirts is a statistical fluctuation. In Section 2.5, we showed that our analysis pipeline produces unbiased results in terms of both Σ profiles and ensemble DK14 fits, and does not create spurious gradient features. Hence, the inferred outskirt feature is unlikely to arise from systematic errors.
An additional source of uncertainties in our scaling analysis is the statistical errors on the NFW parameters of individual clusters, which propagate into uncertainties in their scaling radius R ∆ and scaling density Σ(R ∆ ). In this work, these scaling parameters of individual clusters are fixed to their best-fit values from the NFW fits of U16. Hence, although our scaling analysis has led to significant improvements in the constraints on the inner profiles relative to the conventional stacking (see Section 4.3), these errors are likely to have smeared out local features to some level (Niikura et al. 2015). The degree of smearing can be assessed by marginalizing over uncertainties in the NFW parameters for all clusters, and such effects need to be accounted for in future studies with a larger statistical sample of clusters and with a higher statistical sensitivity.

Comparison with Simulation Results
We begin by making sure that our results for the standard profile parameters, namely the concentration and Einasto shape parameters, are congruent with expectations from simulations. The best-fit values for concentration are relatively independent of the fitting function chosen (Table 2) Since the DK14 profile assumes an Einasto profile for the 1-halo term, the Einasto parameters are of particular interest. The Einasto profile varies in steepness with halo radius, at a rate given by a shape parameter α. Gao et al. (2008) showed that this parameter is, to a good approximation, a function of only the peak height, ν (see also Dutton & Macciò 2014). The CLASH mass corresponds to a peak height of ν 200c = 3.76 and thus α = 0.29, significantly higher than the values of about 0.2 found in our fits. While the results of neither Gao et al. (2008) nor Dutton & Macciò (2014) are well constrained at such extreme peak heights, it is clear that α in their simulations exceeds 0.2 significantly at high ν. This tension already emerged as a ≈ 1σ difference in the analysis of U16. With the improved radial rescaling applied in this paper, the significance of the difference increases to 3.5σ.
We have tested whether the steepening due to the splashback radius can bias the fitted α high compared to the value preferred by the inner profile. We find that such a bias can indeed appear depending on the fitted radial range (a larger range leads to larger bias) and the weights given to each radial bin. However, Dutton & Macciò (2014) fitted out to 1.2r vir and weighted the bins by the number of particles, and for these parameters the bias is smaller than 1%. Another possible explanation is that the X-ray-selected CLASH clusters are preferentially relaxed systems compared to the average population (Meneghetti et al. 2014). However, the relation of Gao et al. (2008) and Dutton & Macciò (2014) for α is based on halo samples that exclude unrelaxed halos, although this selection may not capture the entire effect present in the CLASH sample. On the other hand, taking into account baryonic effects in nonradiative hydrodynamical N -body simulations, Meneghetti et al. (2014) find that the Einasto shape parameter for cluster-size halos lies in the range α = 0.21 ± 0.07, indicating that our measurement is only moderately in tension with simulations.
Finally, we compare our inferences regarding the splashback radius R 3D sp with the simulation results of More et al. (2015) who predicted that the ratio of R 3D sp and r 200m depends primarily on the mass accretion rate of halos, with a less important dependence on redshift. In the left panel of Figure 6, we compare the ratio R 3D sp /r 200m = 1.23 +2.33 −0.34 (blue shaded band), inferred for our sample at z eff l = 0.337, with the mean relation in simulations (gray band), where Γ ≡ d ln M vir /d ln a is the mass accretion rate (More et al. 2015). This comparison shows that the CLASH lensing constraints on the splashback radius overlap with a broad range of mass accretion rates Γ. Our lower 1σ (16th percentile) limit of R 3D sp /r 200m > ∼ 0.89 translates into an upper limit on the accretion rate of Γ < ∼ 4.0. This limit is not particularly informative, since only a very small fraction of halos experience such rapid accretion (DK14).
Since we cannot directly measure the mass accretion rate of the CLASH cluster sample, we cannot independently verify whether this prediction is congruent with observations. For this reason, More et al. (2015) also provide an approximate relation for the mean splashback radius as a function of peak height and redshift, The right panel of Figure 6 shows a comparison of this fitting function with our results for the splashback radius and the peak height of the CLASH cluster sample. Our inferred range of possible splashback radii covers the entire range of values suggested by the simulation results, showing that our results are generally compatible with simulations. While the results of More et al. (2015) were based on DM-only simulations, Lau et al. (2015) broadly confirmed the results of DK14 in hydrodynamical simulations of individual clusters. Thus, we have currently no reason to assume that baryons affect the location of the splashback radius significantly.

Projection Effects in Halo Profiles
In the previous sections, we defined R 3D sp as the halo radius where the logarithmic slope of the three-dimensional density profile is steepest. However, the quantity actually measured is the two-dimensional density profile, both when observing cluster member density profiles as in More et al. (2016) or the lensing signal of clusters as in this work. Thus, it is important to understand the relation between the location of the radii of the steepest slope in thee dimensions and in projection, R 3D sp and R 2D sp . We have investigated this relation using the simulated halo sample of More et al. (2015). In all cases, the steepening of the two-dimensional mass profile is less pronounced than that of the three-dimensional profile, as highlighted in Figure 14 of DK14. At high mass accretion rates (Γ ≥ 3), R 2D sp /R 3D sp approaches a fixed ratio of about 0.8, in agreement with our measurements (Section 3.3). At lower mass accretion rates, however, the R 2D sp derived from the profiles exhibits huge scatter and a seemingly random pattern.
The difficulty in deriving a valid R 2D sp from projected measurements can be understood by considering a few realizations of the DK14 profile with a power-law outer profile representing the 2-halo term (see DK14 for details). The location of the steepest slope in three dimensions is a trade-off between the steepening 1-halo term and the 2-halo term. The steepening is less pronounced for halos with lower mass accretion rates. Thus, in projection, the 2-halo term has a substantial impact on the apparent location of R 2D sp , which emerges at a much smaller radius that is unrelated to the steepening term ( Figure 3), a problem that becomes more serious at small peak heights.
These results highlight the importance of forward-modeling the effects of the steepening based on the underlying threedimensional density profile, rather than attempting to derive R 3D sp from R 2D sp directly (e.g., using Gaussian process modeling).

Compatibility with Measurements from Cluster Member
Density Profiles The splashback radius was first unambiguously detected in observations by More et al. (2016) who stacked surface number density profiles of cluster member galaxies for a large number of clusters. Their clusters were split into two subsamples with high and low concentrations of member galaxies (c gal ) at fixed richness and redshift ).
These high-and low-c gal samples were expected to represent populations of high and low Γ, respectively. However, Zu et al. (2016) have shown that the parameter c gal used in More et al. (2016) is strongly contaminated by projection effects and is likely sensitive to the large-scale environment of the clusters rather than their internal structure. We thus limit our comparison to the full sample of More et al. (2016) whose inferred splashback radius is R 3D sp /r 200m = 0.837 ± 0.031, with a weak-lensing mass of M 200m 1.87 × 10 14 h −1 M at z = 0.24 (S. More 2016, private communication). Figure 6 shows that our lower limit on R 3D sp is higher than this value, but overlaps with the full-sample measurement at the ∼ 1σ level. We note that we expect the average mass accretion rate of the CLASH sample to be low due to a high fraction of relaxed objects (Meneghetti et al. 2014). Hence, owing to the selection effects, there could be a bias toward higher values of R 3D sp /r 200m in our sample.

SUMMARY
We have developed methods for modeling averaged cluster lensing profiles scaled to a chosen halo overdensity ∆, which can be optimized for the extraction of features that are local in radius, in particular the steepening due to the splashback radius in the outskirts of collisionless DM halos. We have examined the ensemble mass distribution of 16 CLASH X-ray-selected clusters by forward-modeling the gravitational lensing data obtained by Umetsu et al. (2016). Our main conclusions are as follows.
• Regardless of the overdensity chosen, the CLASH ensemble mass profile in projection is remarkably well described by a scaled NFW or Einasto density profile out to R ∼ 1.2r 200m ∼ 2.5h −1 Mpc, beyond which the data exhibit a flattening with respect to the NFW or Einasto profile.
• We constrain the NFW halo concentration to c 200c = 3.66 ± 0.11 at M 200c 1.0 × 10 15 h −1 M , consistent with previous work based on the same input data (Umetsu et al. 2016). Our new analysis using scaled profiles provides tighter constraints on the halo shape and structural parameters (c and α) than the conventional stacking.
• We do not find statistically significant evidence for the existence of the splashback radius in the CLASH lensing data. At the current sensitivity, this result is in line with expectations from simulated, synthetic observations. Assuming the DK14 profile form and generic priors calibrated with simulations, we have placed a lower limit on the splashback radius of the clusters, if it exists, of R 3D sp /r 200m > 0.89 or R 3D sp > 1.83h −1 Mpc at 68% confidence. This constraint is in agreement with ΛCDM predictions.
• The gradient feature in the outskirts is most pronounced for a scaling with r 200m , consistent with simulation results of Diemer & Kravtsov (2014) and Lau et al. (2015).
The results obtained here are generally favorable in terms of the standard explanation for DM as effectively collisionless and nonrelativistic on sub-megaparsec scales and beyond, with an excellent match between lensing data and ΛCDM predictions for high-mass clusters. This study represents a first step toward using cluster gravitational lensing to examine detailed predictions from collisionless ΛCDM simulations regarding the shape and universality of the outer density profiles. Such predictions can, in principle, be unambiguously tested across a wide range of halo masses, redshifts, and accretion rates, with large statistical samples of clusters from ongoing and planned lensing surveys such as the Subaru Hyper Suprime-Cam survey (Miyazaki et al. 2015), the Dark Energy Survey, and the WFIRST and Euclid missions.
This work was made possible in part by the availability of high-quality lensing data produced by the CLASH team. We are grateful to the CLASH team who enabled us to carry out the work. We thank all authors of Umetsu et al. (2014Umetsu et al. ( , 2016 and Zitrin et al. (2015) for their contributions to the lensing analyses used here. We thank our referee for his valuable suggestions that improved the paper significantly. We thank Andrey Kravtsov and Surhud More for important suggestions and detailed comments on a draft of this paper. We acknowledge very fruitful discussions with Nobuhiro Okabe, Ho Seong Hwang, Arman Shafieloo, Zuhui Fan, and Congyao Zhang. This work is partially supported by the Ministry of Science and Technology of Taiwan under the grants MOST 103-2112-M-001-030-MY3 and MOST 103-2112-M-001-003-MY3.