Weak lensing constraints on splashback around massive clusters

The splashback radius $r_\text{sp}$ separates the physical regimes of collapsed and infalling material around massive dark matter haloes. In cosmological simulations, this location is associated with a steepening of the spherically averaged density profile $\rho(r)$. In this work, we measure the splashback feature in the stacked weak gravitational lensing signal of $27$ massive clusters from the Cluster Canadian Comparison Project with careful control of residual systematics effects. We find that the shear introduced by the presence of additional structure along the line of sight significantly affects the noise at large clustercentric distances. Although we do not detect a significant steepening, the use of a simple parametric model enables us to measure both $r_\text{sp}=3.5^{+1.1}_{-0.7}$ comoving Mpc and the value of the logarithmic slope $\gamma = \log \rho / \log r$ at this point, $\gamma(r_\text{sp}) = -4.3^{+1.0}_{-1.5}$.


INTRODUCTION
In the concordance ΛCDM model, collisionless dark matter acts as the building block of cosmic structure, contributing about 25% of the total energy density in the Universe and the majority of the total mass (Planck 2015 results, Ade et al. 2016). In this framework, gravity is the primary force behind the growth of structure in the matter field and is able to form the present-day cosmic web from an almost homogeneous initial state. Fully collapsed structures, known as haloes, are thought to grow both through mergers of smaller ones (hierarchical clustering) and continuous infall of ambient dark matter (smooth accretion).
An intuitive understanding of this second mechanism is given by the study of spherical collapse in an expanding Universe (see Gunn & Gott 1972;Fillmore & Goldreich 1984, for some historic landmark results). Shells of material surrounding an overdensity eventually decouple from the Hubble flow and start collapsing towards it. As more shells orbit the halo, the wrapping in phase-space of different streams results in caustics visible in the density profile. Of particular interest is the region around the outermost caustic, where the physical regimes of accreting and collapsed material meet.
More recently, Diemer & Kravtsov (2014, DK14 from now on) studied the spherically averaged density profile ρ(r) of these regions in dark matter only simulations and have reported a change in slope compared to the collisionless E-mail: contigiani@strw.leidenuniv.nl equilibrium profile (Einasto or NFW, Einasto 1965;Navarro et al. 1997). More et al. (2015) argued that the splashback radius r sp , corresponding to the minimum logarithmic slope γ(r) = log ρ(r)/log r, could function as a physically motivated definition for the boundary of dark matter haloes. This role is usually assumed by proxies for the virial radius such as r 200m , defined as the radius inside which the average halo density is 200 times the average matter density of the Universe ρ m . While this radius has a clear definition based on analytical solutions of idealized virialization scenarios, the mass contained within it, known as M 200m , is an imperfect measure of the halo mass. This is because it is subject to a pseudo-evolution caused by the redshift dependence of ρ m (Diemer et al. 2013). In contrast, because the caustic associated with splashback is connected to the apocenter of recently accreted material, all the material within r sp is necessarily collapsed material and should rightfully contribute to the halo mass.
At larger distances, the presence of correlated structure surrounding the halo is expected to shape the density profile. Using the language of the halo model (see, e.g., Cooray & Sheth 2002, for a review), this is a transition region from the 1-halo term to the 2-halo term. DK14 have however reported that in the outermost regions (r 9r 200m ), the 2-halo term based on the matter correlation function provides a worse fit to simulations compared to a simple power-law.
Because the slope of the density profile at r sp is found to be, on average, a decreasing function of the halo mass, DK14 first pointed out that large overdensities are the ideal target for the detection of this feature -i.e., measuring a significant departure from the equilibrium profile. This makes galaxy clusters the ideal candidates since they correspond to the most massive haloes in the Universe. For this mass range, r sp is expected to be located around r 200m , at a cluster-centric distance of the order of a few Mpc.
The splashback feature should also be present in the radial distribution of galaxies. This was first detected by More et al. (2016) using the large sample of redMaPPer clusters from Rykoff et al. (2014), and studied further in Baxter et al. (2017). However, these studies find a discrepancy between the inferred splashback radius and the expected distribution of subhaloes based on dark matter only simulations. Known physical processes (e.g., tidal disruption and dynamical friction) are not expected to induce a mismatch between the galaxy and subhalo distributions at splashback scales and this deviation is still unexplained. In particular, while the results have been shown to depend on the details of the cluster finding algorithm (Zu et al. 2017;Busch & White 2017), it is still uncertain if this can fully explain the discrepancy (Chang et al. 2018). Chang et al. (2018) studied a sample of redMaPPer clusters selected in Dark Energy Survey year 1 data. For this large sample, they detected a splashback feature in the galaxy distribution and from weak lensing measurements. The latter has the advantage that the lensing signal probes the matter distribution directly (see e.g. Hoekstra et al. 2013, for a review). The first attempt to detect the splashback feature using weak gravitational lensing was presented in Umetsu & Diemer (2017), who used a sample of 20 highmass clusters in the Cluster Lensing and Supernova survey with Hubble (CLASH). Unfortunately, the limited field of view (foV) of Suprime-Cam prevented precise measurements in the outer regions, and as a result Umetsu & Diemer (2017) could only provide a lower limit on the splashback radius.
In this work we provide a measurement 1 of splashback using weak lensing observations for a sample of 27 massive clusters of galaxies that were observed as part of the Canadian Cluster Comparison Project (CCCP; Hoekstra et al. 2012). Hence our strategy is similar to that employed by Umetsu & Diemer (2017), but we take advantage of the fact that the CCCP observations were obtained using MegaCam, which has a foV of 1 deg 2 , and enables use to measure the lensing signal beyond the splashback radius. The paper is organized as follows: in Section 2 we present our dataset and describe our lensing analysis, in Section 3 we show the results of our fit and the implications for splashback, and in Section 4 we draw our conclusions. Throughout the paper we employ a flat ΛCDM cosmology with H 0 = 70 Mpc/km/s, Ω m = 0.3, Ω c = 0.25 and σ 8 = 0.80.

CLUSTER LENSING
In this section we discuss how the sheared images of distant galaxies can be used to constrain the matter distribution of clusters along the line of sight. After introducing our cluster 1 In the interest of reproducibility we make our splashback code publicly available at https://github.com/contigiani/splash/ sample, we present the weak lensing measurements and explain our methodology, with a particular focus on systematic effects and noise estimation.

Sample characterization
Our dataset is based on the Canadian Cluster Comparison Project (CCCP), a survey targeting X-ray selected massive clusters at z 0.5 introduced for the first time in Hoekstra et al. (2012) and re-analysed in Hoekstra et al. (2015, H15 from now on). The starting points of our analysis are the r-band images of 34 clusters captured by MegaCam at the Canada-France-Hawaii Telescope (CHFT). We exclude from this initial sample the entries corresponding to on-going mergers: Abell 115, Abell 222/3, Abell 1758, and MACS J0717.5+3745 because these systems display a visible double peaked matter distribution for which two splashback surfaces might intersect each other. After this cut, we are left with 27 massive clusters, which compose the sample used in this study.
The objects are characterized by masses 3.8 < M 200m /(10 14 M ) < 26.4 and cover a redshift range 0.15 < z < 0.55, with only 6 of them located at z > 0.3. Table 1 reviews the sample and presents the quantities relevant for the present work. For more details about the cluster sample we refer the reader to Hoekstra et al. (2012), H15 for a description of the weak lensing analysis, and the companion paper Mahdavi et al. (2013) for the analysis of X-ray observations.
In simulations, DK14 found a correlation between the splashback feature and the halo mass. We therefore define a high-mass subsample of our clusters, containing the 13 most massive objects. The average M 200m of the sample and the subsample, weighted by signal-to-noise ratio (SNR), equal 1.7 and 2.0 × 10 15 M , respectively. We choose to employ the gas mass M g within r 500c reported by Mahdavi et al. (2013) to define our high-mass threshold. This is because this value is found to be a robust estimator of the weak lensing mass and is statistically almost independent from it. A weak dependence between the two is left due to the lensing-based definition of r 500c .
Targeted observations such as the ones discussed in this work currently represent the most efficient approach to study clusters of virial mass around 10 15 M . In particular, such a sample cannot be obtained by present-day or nearfuture wide surveys, e.  (Miyatake et al. 2016) and 3684 clusters with M 200m = 3.6 × 10 14 M for DES Y1. In contrast, our dataset is much closer in nature to the CLASH sample used in Umetsu & Diemer (2017), also based on targeted observations. In particular, the mass of their stacked ensemble, M 200m = 1.9 × 10 15 M , matches ours. Nevertheless, we want to mention one feature unique to CCCP: the FoV of MegaCam (1 × 1 deg) is significantly larger than that of Suprime-Cam (34 × 27 arcmin), the instrument used for the CLASH profile reconstruction at large scales (Umetsu et al. 2016). This is particularly suited for  Table 1. The full cluster sample, "CCCP all", used in this paper. RA and DEC are the sky position of the cluster centre (brightest cluster galaxy, or X-ray peak for coordinates marked with §), z is the cluster redshift, β is the average value of D LS /D S (see Sec. 2.2), M g is the gas mass within r 500c , defined as the radius of the sphere inside which the mean halo density is 500 times the critical density of the Universe at redshift z and M 200m is the mass enclosed within r 200m . These values are recovered from the NFW fit performed in H15. The values for M g are taken from the X-ray analysis of M13 or, for values marked with †, they are defined using the scaling relations found in the same paper. Clusters listed below the horizontal line belong to the high mass subsample.
our purposes since it allows us to better cover cluster-centric distances where the splashback radius is located.

Tangential shear
In the weak lensing regime the shear field is found by averaging the PSF-corrected ellipticities of a sample of background sources. We follow H15 and use sources in the magnitude range 22 < m r < 25. The lower limit reduces the presence of foreground objects such as bright galaxies belonging to the clusters, which are abundant in the central regions and are not sheared by the cluster's mass distribution. Because this operation is unable to completely remove cluster members, we chose to model the residual contamination statistically, as explained in Sec. 2.3. Shapes are measured using an improved KSB method (Kaiser et al. 1995;Luppino & Kaiser 1997;Hoekstra et al. 1998). The quadrupole moments of the galaxy images are used to construct a polarization tensor e, which is then corrected for the point spread function (PSF) of the observing instrument. In Section 2.3 we address this step in more detail and mention the improvements we have implemented since H15. The shear polarizabilityP γ quantifies how the observed polarization of an individual galaxy is related to the gravitational shear. For an ensemble of sources the shear components are hence measured as a noise-weighted average, e i /P γ , where the individual weights are written as (Hoekstra et al. 2000) w = 1 2 + σ e /P γ 2 . (1) In this expression two sources of noise are included: the scatter introduced by the intrinsic variance of galaxy ellipticities 2 and the uncertainty in the measured polarization σ e due to noise in the imaging data. Following Hoekstra et al. (2000) we use 2 1/2 = 0.25.
For an isolated circular overdensity, the induced shear is purely tangential, i.e., the deformation is parallel to the radial direction. In general, this shear component is related to the projected mass surface density Σ(R) as a function of the radial coordinate R: In these expressions the profile ∆Σ(R) is called excess surface density (ESD) and the critical density Σ cr is a geometrical factor quantifying the lensing efficiency as function of the relative position of source and lens. The definition above  Figure 1. Lensing signal. The top panel shows the noise-weighted stacked signal of the 27 clusters in our sample as a function of comoving clustercentric distance, together with a best fit NFW profile to the first five data points (see Sec 2.2 for more information). The arrow points to the inferred location of r 200m ; in galaxy clusters the splashback feature is located around this position. The larger error bars are the full 1σ errors for the data points, while the inner error bars account only for statistical uncertainty. The difference between the two is apparent only in the last few data points. The bottom panel shows an estimate of the expected residual systematics left after the corrections discussed in Sec. 2.3 are applied, expressed as a fraction of the total uncertainty. These effects are found to be consistent with the error bars.
applies for a lens at distance D L shearing an ensemble of sources. β is the average of the quantity max [0, D LS /D S ] for each source, with D LS being the individual lens-to-source distance 2 and D S the distance to the source.
Because we work with single-band observations, we are unable to derive individual photometric redshifts. Fortunately, a representative photometric redshift distribution is sufficient to estimate β. This distribution is obtained for all clusters by magnitude-matching the most recent COSMOS photmetric catalogue (COSMOS2015, Laigle et al. 2016) to our source r-band magnitude range.
We point out that the measured average ellipticity is an estimator of the reduced shear however, because the splashback feature we focus on in this work is located in low-density regions, we simply assume 2 Note that D LS is negative for foreground sources.
that the average ellipticity is a measure of the shear γ i . From our source catalogs we extract the tangential component γ t (θ j ) in radial bins and estimate for each cluster the data covariance matrix as as the sum of two terms: The first is a diagonal matrix accounting for the statistical error on the weighted average of the measured ellipticities and the second quantifies the additional shear variance caused by the presence of cosmic structure between viewer and source (Hoekstra 2003;Umetsu et al. 2011) where P κ ( ) represents the projected convergence power spectrum for the multipole number . For an angular bin θ extending from θ − to θ + , g(l, θ) is defined using the Bessel functions of the first kind of order zero and one, J 0 and J 1 : For a given cosmology, P κ ( ) can be evaluated using the Limber projection starting from a source redshift distribution and a model for the non-linear matter power-spectrum (Kilbinger 2015). For this work this is done using CAMB 3 (Lewis 2013) and HALOFIT (Takahashi et al. 2012). A third term accounting for the intrinsic variance in a particular realization of galaxy clusters should be added to the matrix in Eq. 5. For massive clusters in the considered redshift range, this term is found to be dominated by Poissonian scatter in the number of haloes contained within the correlated neighbourhood (Gruen et al. 2015). We neglect this term because in similar lensing analyses (e.g., Umetsu et al. 2016;Miyatake et al. 2018) it is always found to be sub-dominant to statistical and large-scale structure noise, especially on the scales of interest for this work.
The top panel of Fig. 1 presents the average noiseweighted signal of the full cluster sample. The double error bars in the figure illustrate how the inclusion of C lss has an impact on the uncertainties at large scales. An NFW fit, obtained using the virial overdensity from Bryan & Norman (1998) at an assumed redshift z = 0.25, is also shown. The position of r 200m for the best-fit model is also indicated in the same figure.

Residual systematics
In this section we address the effects of the corrections we have implemented to tackle two systematic effects that are particularly important for our analysis: PSF anisotropy and cluster member contamination. In particular, we estimate the amplitude of any residual systematic effects as plotted in the bottom panel of Fig. 1. In the KSB method, the observed galaxy polarizations are corrected for PSF anisotropy using where the smear polarizability P sm quantifies how susceptible a source is to PSF distortions and p j is the PSF anisotropy measured using point-like sources (see, e.g., Hoekstra et al. 1998).
The observed polarizations and polarizabilities are, however, biased because of noise in the images. If unaccounted for, this leads to biased cluster masses. For the shear, these corrections can be expressed in terms of a multiplicative and additive bias, µ and b: To ensure accurate mass estimates, H15 focused on the impact of multiplicative bias. To do so, they used image simulations with a circular PSF to calibrate the bias as a function of source SNR and size. However, the actual PSF is not round and H15 therefore quantified the impact of an anisotropic PSF on the multiplicative bias correction. This study also revealed, however, that about 4 per cent of the PSF ellipticity was still visible in the average shear as additive bias (see their fig. A1). This is because because the smear polarizability itself is affected by noise.
For this study we have improved our calibration by empirically increasing P sm by a factor 1.065, such that no residual additive bias remains visible, see Fig. 2. We also verified that this latest correction does not introduce significant trends with source characteristics. We use the difference between the ensemble lensing signal measured before and after this improvement as a (conservative) estimate of any unknown systematics affecting the shape measurement method.
The second effect we account for is the presence of cluster members in our source catalogues. Note that in this case, we have not updated the methodology from H15, but we still report it here for completeness. If we assume that cluster members are randomly oriented, as found by Sifón et al. (2015), their inclusion among our sources has the effect of diluting the measured shear. To correct for this, we multiply γ t (R) by a boost factor B(R) defined as a function of the projected comoving distance R: The contamination term f cont accounts for the decrease of the ellipticity average due to the presence unsheared sources and, by comparison with blank fields, it is found to be 1) a decreasing function of distance from the cluster centre and 2) negligible beyond a distance r max . An extra factor f obs is also introduced to model the reduced background galaxy counts due to obscuration by cluster members. This factor is computed by stacking the cluster images with simulated blank fields and measuring how many simulated sources are obscured.
The functions appearing in the boost factor are written empirically as: 1 f obs (R) = 1 + 0.021 0.14 + (R/r 500 ) 2 and (11)  where n 0 and R c are fitted independently for each cluster and B(R) = 0 for R > r max ≡ 4(1 + z) Mpc.
To quantify the amplitude of residual systematics for this second correction, we refer to H15, where a residual scatter of about 2 per cent around the ensemble correction was reported.

SPLASHBACK
In this main part of the paper we fit the observed weak lensing signal using the spherical density profile presented in DK14. This profile is designed to reproduce the expected flattening of the density profile at large radii due to the presence of infalling material, as seen in numerical simulations.

Fitting procedure
The projected surface density profile Σ(R) for a spherical lens with matter density ρ(r) is: where we limit the integration range to [0, 40] Mpc for our numerical calculations. For cosmological overdensities, this profile can be connected to the lensing signal through eq. 2. In this section we use a model for ρ(r) with the following components: an Einasto profile ρ Ein (Einasto 1965) to model the inner dark matter halo, a transition term f trans (r) to capture a steepening effect at the halo edge and a powerlaw ρ out (r) to model the distribution of infalling material in the outer regions. The mathematical expressions are the following: In DK14 the infalling term includes an offset corresponding to the average matter density, but this is not present in our fitting function because the tangential shear in eq. 2 is completely insensitive to it.
In its general form, this model depends on a large number of parameters. In order to reduce its degrees of freedom we therefore choose to set strong priors on a few parameters. As done in Baxter et al. (2017) and Chang et al. (2018) we do not fit both ρ 0 and r 0 , but choose to fix one of them, as they are degenerate. We impose Gaussian priors log(0.2)±0.1, log(4) ± 0.2 and log(8) ± 0.2 on the logarithms of the exponents log α, log β and log γ respectively. The loose prior on the Einasto shape parameter α is motivated by dark matter only simulations and its 1σ interval covers the expected scatter due to the redshift and mass distribution of our sample (Gao et al. 2008;Dutton & Macciò 2014), while for the exponents in the transition term the stringent priors are centered on the values suggested by DK14. We also set a Gaussian prior on the truncation radius r t , 4 ± 2, based on the same results. The location of the median is based on the r 200m inferred from our NFW fit and the selected standard deviation covers the expected range due to the mass distribution of our sample. Finally, based on previous measurements, we also set a minimum value of 1 for the outer slope s e and a physically motivated minimum value of 0 for the density parameters ρ s and ρ 0 .
A rescaling of the radial coordinate with an overdensity radius (e.g. r 200m ) is often employed when fitting the profile described above. While the uncertainties on individual cluster profiles do not allow us to fully account for this effect when stacking, we still attempt to remove the redshift dependence of the average matter density of the Universe by using comoving coordinates.
We follow Umetsu & Diemer (2017) and do not include a miscentering term in our tangential shear model. In general, a shift in position of the cluster centres reported in Table 1 would cause a smoothing of the lensing profile in the central region. An estimate of the area affected by such an effect can be obtained by considering the difference between two independent estimators of the halo centre: the position of the brightest cluster galaxy or the X-ray luminosity peak. Our sample is found to be well centered (see M13) with the root mean square of the offset between the two σ off = 33 kpc. For the scales plotted in Fig. 3 we therefore do not expect our data to be affected by miscentering.
A fit the to input data γ t (R) with the covariance matrix defined in Sec. 2 is performed by sampling the posterior distribution of the parameters [ρ s , r s , log α, r t , log β, log γ, ρ 0 , s e ] using the Markov Chain Monte Carlo ensemble sampler emcee 4 (Foreman-Mackey et al. 2013, based on Goodman & Weare 2010). For both CCCP samples considered a minimum of the slope is identified. At larger distances, the results are the least interesting. In these regions, the power-law term becomes dominant and the value of the slope is set exclusively by the exponent s e . In particular, its lower limit is artificially imposed by our prior.

Interpretation
What is more relevant to our study is the minimum value of the slope γ(r) and its location, i.e., the splashback radius r sp . The 68 per cent credible intervals of both quantities are indicated as shaded sections of the vertical and horizontal histograms. Our measured 99.7 per cent confidence interval of γ(r sp ) for the full sample is [−12.1, −2.3], meaning that we are unable to measure a significant departure from the slope expected for an NFW profile (about −2.5). Despite this, we are still able to constrain the value of both the splashback radius and the logarithmic slope at this point, r sp = 3.6 +1.2 −0.7 Mpc and γ(r sp ) = −4.2 +1.0 −1.8 . We also highlight that the high-mass sample returns similar constraints with only half the sample size, r sp = 3.5 +1.2 −0.8 and γ(r sp ) = −3.6 +1.0 −1.8 . As a point of reference, we also show the expected profiles from a suite of zoom-in hydrodynamical simulations of massive clusters (Hydrangea, Bahé et al. 2017). From the full Hydrangea sample, we have selected the 8 most massive clusters for this comparison in order to obtain a sample with an average value of M 200m = 1.7 × 10 15 M , similar to our dataset, but evaluated at z = 0. Our slope measurements are found to be agreement with what is seen in simulations. Note that the amplitude of the signal plotted in Fig. 3 is lower than the observed sample due due to the different redshift.
As done in Umetsu & Diemer (2017), we study the impact of the the model parameters on the predictions for r sp and γ(r sp ) to verify that our dataset is informative and we 4 https://emcee.readthedocs.io/. Gaussian prior on r t Flat prior on r t Figure 4. Impact of the prior on the truncation radius r t on our results. The corner plot presents the two-dimensional and maginalized posterior distributions for the DK14 parameter r t , the inferred splashback position r sp , and logarithmic slope γ(r sp ). If, instead of a Gaussian prior (dashed red line), a flat prior is assumed (dashed black line), the parameter r t has no upper bound. This translates into weaker constraints on r sp .
are not simply sampling our model priors. Of crucial importance is the truncation radius r t , which, in the original definition of the DK14 profile, explicitly sets the position of the splashback feature. Similarly to Umetsu & Diemer (2017), we also find that we are unable to fully constrain this parameter. This can be seen in Fig. 4, where we plot the posteriors of three relevant parameters for two different choices of the r t prior: the Gaussian assumed in our main study and a flat prior in the range [0, 20] Mpc. While the posterior for γ(r sp ) (middle row) is mostly unaffected by this choice, we obtain a looser upper limit on the splashback radius (top panel) in the second case: r sp = 3.9 +2.5 −1.0 . As visible in the bottom-left panel, this is due to a clear correlation with r t .
We find no correlation between r sp and r t for r t 10 Mpc. In this regime the location of the minimum of γ(r) is controlled by the presence of the infalling term ρ in (r) ∝ r −s e . Because the slope s e is relatively gentle, if r t is large enough the truncation happens in a region dominated by the infalling material and cannot be constrained. Because the truncation is expected to be visible in the transition regime, our Gaussian prior on r t effectively forces it to a physically motivated position and, from the figure, we confirm that it does not introduce a biased posterior peak.

CONCLUSIONS
We have shown in this work that targeted weak lensing observations of massive clusters can be used to measure the splashback feature and that particular care is required when correcting for residual PSF contaminations, which should be well understood, and estimating the data covariance matrix, which should take into account the presence of additional structure along the line of sight. Using a stack of 27 massive clusters from CCCP we have fully constrained for the first time the splashback radius around massive clusters, r sp = 3.6 +1.2 −0.7 , and similar precision has also been achieved with as little as 13 objects. We stress that, because of the purely gravitational nature of weak lensing, minimal assumptions are required to interpret our signal.
In the last few years, the study of the physics of accretion at the outskirts of massive dark matter haloes has become observationally viable. Splashback offers a unique view into the phase-space configuration of haloes, which has not yet been explored in observations. In particular, the physics behind it appears to be remarkably uncomplicated and semianalytical models of spherical collapse for cold dark matter are able to reproduce the expectations from N-body simulations (e.g., Adhikari et al. 2014;Shi 2016). The fact that these results are based only on the dynamics of collapsing dark matter in an expanding Universe makes splashback a remarkable prediction of general relativity and dark matter. More generally, its connection to the growth of cosmological structures makes it a test for ΛCDM. As an example, it has also been shown recently that modifications of gravity have a significant impact on this feature (Adhikari et al. 2018). As the first results are starting to appear in the literature, we argue that splashback solicits further investigation exactly because it is a falsifiable prediction of the current paradigm.
We found that at the relevant scales a significant contribution to the lensing signal is cosmic noise. In the near future, this term can be reduced significantly with larger cluster samples. Looking further ahead, deep wide-area surveys such as Euclid (Laureijs et al. 2011) and LSST (LSST Science Collaboration et al. 2009) will provide unprecedented depth and survey area, and thus deliver the data required to study splashback over a wider mass and redshift range.