Anna Niemiec
email: anna.niemiec@lam.fr

Eric Jullo

Marceau Limousin

Carlo Giocoli

Thomas Erben

Hendrik Hildebrant

Jean-Paul Kneib

Alexie Leauthaud

Martin Makler

Bruno Moraes

Maria E S Pereira

Huanyuan Shan

Eduardo Rozo

Eli Rykoff

Ludovic Van Waerbeke

Stellar-to-halo mass relation of cluster galaxies

Keywords: Cosmology, Lensing, Galaxy Clusters 1 INTRODUCTION

In the formation of galaxy groups and clusters, the dark matter haloes containing satellite galaxies are expected to be tidally stripped in gravitational interactions with the host. We use galaxy-galaxy weak lensing to measure the average mass of dark matter haloes of satellite galaxies as a function of projected distance to the centre of the host, since stripping is expected to be greater for satellites closer to the centre of the cluster. We further classify the satellites according to their stellar mass: assuming that the stellar component of the galaxy is less disrupted by tidal stripping, stellar mass can be used as a proxy of the infall mass. We study the stellar to halo mass relation of satellites as a function of the cluster-centric distance to measure tidal stripping. We use the shear catalogues of the DES science verification archive, the CFHTLenS and the CFHT Stripe 82 surveys, and we select satellites from the redMaPPer catalogue of clusters. For galaxies located in the outskirts of clusters, we find a stellar to halo mass relation in good agreement with the theoretical expectations from Moster, Naab & White (2013) for central galaxies. In the centre of the cluster, we find that this relation is shifted to smaller halo mass for a given stellar mass. We interpret this finding as further evidence for tidal stripping of dark matter haloes in high density environments.

part of their dark matter (Hayashi & White 2008;Giocoli, Tormen & van den Bosch 2008). Subhaloes which have been accreted earlier have experienced tidal stripping from the host for a longer time, and have thus lost a higher fraction of their initial mass. See also Gao et al. (2004); van den Bosch, Tormen & Giocoli (2005); Limousin et al. (2009); Giocoli et al. (2010).

In addition, there is a correlation between the distance of the subhalo to the cluster centre and the time since accretion, with subhaloes accreted earlier residing on more tightly bound orbits. This effect is due on one hand to the insideout assembly of dark matter haloes (i.e. at higher accretion redshift the host halo was smaller, leading to satellites being accreted at smaller cluster-centric distance compared to later redshifts (Tormen 1998)), and on the other hand to dynamical friction that slows subhaloes down and make them sink into the centre of the host as time passes (Gao et al. 2004). So observationally, to study the evolution of subhaloes during infall, we can use the distance from the satellite to the centre of the cluster as an indicator of the time since the accretion of the satellite. In observations, the distances that are measured are projected along the line of sight, but van den Bosch et al. (2016) shows the correlation between this projected distance and the accretion redshift is still very strong, although weaker than for the 3-dimensional distance.

Note that van den Bosch et al. (2016) shows that the segregation is much stronger if the present subhalo mass is normalized by its mass at accretion. Indeed looking at the global satellite population, tidal stripping and dynamical friction have opposite effects on the radial distribution of subhalo masses in the cluster. While stripping tends to reduce more strongly the mass of satellites close to the centre, the satellites which are bigger at accretion are subjected to stronger dynamical friction and thus more driven to the centre.

Since in observations the mass at accretion is not an observable, we use as a proxy the stellar mass. As the host tidal forces strip preferentially the outer part of the subhaloes, the baryonic part at their centre is not significantly disturbed. While it is also possible that stellar mass is created during infall, most galaxies have their star formation quenched at accretion (Zu & Mandelbaum 2015, 2016). The stellar mass is thus well correlated to the subhalo mass at accretion (Vale & Ostriker 2006;Cooray 2006;Behroozi, Conroy & Wechsler 2010;Smith et al. 2016).

In summary, we measure M DM sat /M star sat as a function of the projected distance from the satellite to the cluster centre, to see how the dark matter halo of the satellites is affected by the tidal stripping of the host throughout infall.

A well-established tool to measure the total projected mass of a galaxy, including the dark matter halo, is gravitational lensing. The light coming from background sources is deviated when it passes by a massive (lens) object, and the amplitude of the deviation depends on the total mass of the lens object. Lensing measurements of the mass of an object are independent of its nature (baryons, dark matter, etc.) and its state (equilibrium, etc.).

The first measure of the mass of subhaloes in a cluster was performed by Natarajan et al. (1998) in the cluster AC114, measuring the perturbation by subhaloes of the shear distribution of background sources. In a later work on the cluster Cl0024+16 using HST images, Natarajan et al. (2009) present the first indication for tidal stripping on the dark matter haloes of satellite galaxies. Later, Okabe et al. (2014) performed a similar analysis for the Coma cluster using data from the Subaru telescope.

Alternatively, the dark matter mass of satellites can be measured in a statistical way using galaxy-galaxy weak lensing over a large sample of galaxies. It consists of the measurement of the average tangential shear in the shape of background galaxies in circular bins centred on a lens object, in order to measure the mass of the lens. In galaxygalaxy lensing the lens is the dark matter halo of a galaxy, and so has a relatively low mass: the shear induced by a single galactic halo is very low, so to obtain a measurable lensing signal the shear induced by numerous lenses needs to be stacked together (Brainerd, Blandford & Smail 1996). Yang et al. (2006) first suggested to apply the galaxygalaxy lensing method to measure the mass of subhaloes in clusters, and Li et al. (2014) first measured the mass of satellite galaxies in groups using data from SDSS and the CFHT Stripe 82 (CS82) survey. Using this method, Gillis et al. (2013) measured the total mass of galaxies in low and high density environments, arguing that the galaxies in high density environment were subject to tidal stripping, compared to the ones in low density environment. In the 100 deg 2 overlap between the DR2 of the KiDS lensing survey and the GAMA spectroscopic survey Sifón et al. (2015) repeat the analysis but the statistical errors prevent them from measuring tidal stripping of dark matter for galaxies in clusters with a mass M host > 10 13 h -1 M . Finally, in the 170 deg 2 of the CFHT Stripe 82 survey, Li et al. (2016) measured the mass of the dark matter haloes and the evolution of the mass-to-light ratio for satellite galaxies as a function of their distances to the centres of the redMaPPer clusters, and obtain a significant signal consistent with tidal stripping. Finally, van Uitert et al. (2016) measured the stellarto-halo mass relation for galaxies from the GAMA survey, comparing central to satellite galaxies, and found no significant difference.

We extend these works by considering the stellar mass distribution of the satellites, and further splitting our galaxy samples according to it. As mentioned above, it appears from simulations that for satellites, the quantity which is more segregated with respect to the distance to the centre of clusters is the present subhalo mass normalized by the mass at accretion. We use the stellar mass as an observational proxy for the mass at accretion. This paper is organized as follows: in Section 2, we present the lensing method and the halo model we use to fit the data. In Section 3, we present our source and lens catalogues, then in Section 4 we show the results of the analysis, and we discuss them in Section 5. We conclude in Section 6.

We assume a WMAP7 cosmology (Komatsu et al. 2011), with (ΩM, ΩΛ, h, σ8, w) = (0.27, 0.73, 0.70, 0.81, -1). When relevant, the dependence on h is clearly stated.

METHODOLOGY

Lensing

In the weak lensing regime, the distortion induced on the image of source galaxies by a single lens galaxy is so weak that it is too tiny to be detected: the observed shape of a background source is the sum of its intrinsic ellipticity and of the lensing shear, and the shear represents only a few per cent of the total ellipticity. However, by stacking the measurements of many individual lensing signals together, the average tangential shear generated by a sample of lenses can be detected.

We measure the galaxy-galaxy lensing observable, the excess surface mass density profile ∆Σ(R) in comoving units as a function of the distance to the (stacked) lenses, and this quantity is related to the tangential shear γt by:

∆Σ(R) = Σcritγt(R) = Σ(< R) -Σ(R), (1)
where Σ(R) is the mean surface density at a projected distance R from the lens centre, Σ(< R) is the mean surface density in a disk of radius R centred on the lens, and Σcrit is the critical surface density in comoving units. It is defined as

Σcrit = c 2 4πG Ds D l D ls 1 (1 + z l ) 2 (2)
where D l and Ds are the angular diameter distances respectively to the lens and to the source, and D ls is the angular diameter distance between the lens and the source. z l is the redshift of the lens.

Halo Model

The lensing signal around the satellite galaxies is described by a four-term halo model (e.g. Cooray & Sheth 2002;Sheth & Jain 2003;Giocoli et al. 2010;Gillis et al. 2013):

∆Σ = ∆Σstar + ∆Σsat + fsat∆Σ host + ∆Σ 2h . ( 3
)
Here ∆Σstar is the baryonic component of the satellites, ∆Σsat corresponds to the contribution of the satellite dark matter (ie the subhalo), ∆Σ host is the contribution of the host (cluster) haloes of the satellites, and ∆Σ 2h is the twohalo term, produced by the neighbouring haloes. The factor fsat represents the fraction of satellites in the lens sample, and we fix fsat = 1 as we study the lensing signal produced by a theoretically pure satellite sample. We note that in the case fsat = 1 the term ∆Σsat should actually be written fsat∆Σsat + (1 -fsat)∆Σcent, with a contribution from both satellite (in subhaloes) and central galaxies (in haloes). We write only ∆Σsat for clarity. We do not take into account a central baryonic component in the host clusters.

The model is expressed in comoving units. We describe the different terms in this section.

The stellar component

We consider the baryonic component of the satellite galaxy mass to be a point source, with mass M * equal to the median stellar mass of the sample. Thus:

∆Σstar(R) = M * R 2 . (4)

The satellite term

The galaxies are assumed to live in dark matter haloes, characterized by a NFW density profile (Navarro, Frenk & White 1996). In the case of satellite galaxies, we use a smoothly truncated NFW profile, whose spatial density distribution is defined by (Baltz, Marshall & Oguri 2009):

ρ(r) = ρcritδc (r/rs)(1 + r/rs) 2 (1 + (r/rt) 2 ) , ( 5
)
where rt is the truncation parameter which ensures that the total NFW mass does not diverge. We fix it at rt = 2r200 (Hilbert & White 2010; Oguri & Hamana 2011), a value for which the profile hardly deviates from a standard NFW inside the virial radius. This profile has two free parameters: the halo mass M200 and the concentration c = r200/rs. The virial radius r200 defines a sphere with a density 200 times higher than the critical density of the Universe, which gives

ρcrit = 3 800π M200 r 3 200 , (6)
and δc can be expressed as

δc = 200 3 c 3 ln(1 + c) -c 1+c . ( 7
)
To reduce the number of free parameters, we use a massconcentration relation, defined at z = 0 in Neto et al. (2007):

c0 = 4.67 × M200 10 14 h -1 M -0.11 (8)
and add the redshift dependence c = c0/(1 + z). To verify the influence of this choice of a mass-concentration relation, we compare ∆Σsat computed using different massconcentration relations (Neto et al. 2007;Duffy et al. 2008;Dutton & Macciò 2014;Shan et al. 2017) at Msat = 10 12 h -1 M , and find a maximum deviation in the amplitude of the signal of ∼ 1%, which is negligible compared to the error bars on our measurements.

The host halo term

As we measure the lensing signal of satellite galaxies in clusters, an important component of the total signal is the one induced by the host clusters. We stack the lensing signal for satellites located at different projected distances to the centre of their cluster. The mean host halo signal is expressed as

∆Σ host (R, M host ) = r 2 r 1 ∆Σ 1host (R, Rs, M host )P (Rs) dRs r 2 r 1 P (Rs) dRs , (9)
where ∆Σ 1host (R, Rs, M host ) is the mean contribution of one host with mass M host located at a projected distance Rs from the satellite, and P (Rs) is the probability for a satellite to be located at a distance Rs from the centre of its host cluster. Therefore the term ∆Σ host (R, M host ) is the mean contribution of a host located at a distance Rs ranging between r1 and r2, and weighted by P (Rs) which is described later. Figure 1 shows a sample of lensing signals produced by a host located at different distances to the satellite, both with and without the weighting. We assume that all the host haloes have the same mass, and we measure this average mass M host .

The contribution of one host halo located at a distance Rs from its satellite is defined as

∆Σ 1host (R, Rs) = Σ1host (< R, Rs) -Σ1host (R, Rs) (10) 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 R[h 1 Mpc] 5 0 5 10 15 DS[M hpc 2 ] P(R s ) = N sat (R s ) P(R s ) flat Figure 1.
Lensing signal produced by host haloes at different distances from the satellite Rs ∈ [0.12; 0.52]h -1 Mpc. The dashed red curves represent hosts which all have the same weight, whereas for the blue lines hosts are weighted by P (Rs) = Nsat(Rs). The total host term is a continuous sum of weighted single halo terms, as in equation 9. For this plot we fixed M host at 10 14 h -1 M and use Nsat(Rs) from the satellites with 10.5 < log(M * /M ) < 11 and 0.1 < Rs/h -1 Mpc < 0.55.

∆Σ 1host (R, Rs) = 1 πR 2 R 0 R 2π 0 Σ 1host (R g ) dθ dR - 1 2π 2π 0 Σ 1host (Rg) dθ (11)
where Σ 1host (Rg) is the projected surface density of the host halo, a NFW halo of mass M host , measured at a distance

Rg = | R -Rs| = √ R 2 + R 2
s -2RRs cos θ from the host centre, where θ is the angle between the vector joining the centre of the satellite to the centre of the host, and the vector joining the centre of the satellite to the point of measurement. As we average the signal in circles or disks centred on the satellite, we integrate over θ.

To obtain the probability function P (Rs), we use the distribution of satellites in the data, stacked for all the clusters, in each bin of cluster-centric distance and stellar mass, ie. P (Rs) = Nsat(Rs). In order to smooth the obtained distribution, for each stellar mass bin we fit the distribution with two second order polynomials, one for each clustercentric bin. The Nsat distributions and the best-fit curves are shown in figure 2.

Satellite concentration. Alternatively, we fit the surface number density profile of the satellites with a NFW profile to measure their concentration. We plot the surface number density profiles Σsat, for each of the three stellar mass bins on figure 3. Each galaxy is weighted by the probability that it is a true member of its host Pmem. The left panel of the figure presents the measurement for the whole satellite sample, while the right panel plot is for our fiducial sample which contains only the galaxies with Pmem > 0.8 (see section 3.2). For each profile we fit a projected NFW profile. We normalize the measurement and the model to have the same amplitude around the scale radius. 10 -1 10 0

R s [h -1 Mpc] 10 -1 10 0 Σ sat [h 2 Mpc -2 ] All P mem 10 -1 10 0 R s [h -1 Mpc] P mem > 0. 8 Figure 3
. Surface number density of satellites in the three stellar mass bins. The continuous lines represent the stacked measurements, where each satellite is weighted by its probability of membership to the host (see 3.2). The left plot presents the measurements for all the satellites, while on the right plot the measurement is made for our fiducial sample, which contains only the satellites with Pmem > 0.8. The dashed lines are the corresponding best-fit NFW profiles.

Once normalized, the profiles for galaxies in the different stellar mass bin are similar to each other, and the best-fit concentrations we obtain have very close values between the different bins. From low to high stellar mass bins, we obtain concentrations of 3, 3.4 and 3.7 for the complete satellite sample, and 4.3, 4.4 and 4.5 for the sample with Pmem > 0.8.

While the full satellite sample distribution is well described by a NFW density profile, the sample with Pmem > 0.8 is steeper, which is why we fit the satellite distribution with polynomials in the model.

We note that if we do not weight the galaxies by their membership probability, the full sample gives a flatter distribution than the NFW profile, and the Pmem selected sample remains steeper.

The two-halo term

On large scales the lensing signal is dominated by neighbouring mass concentrations, e.g., nearby haloes and filaments. The contribution of the neighbours is accounted for by the two-halo term, defined as:

∆Σ 2h (R) = Σ2h (< R) -Σ2h (R) (12) with Σ2h (R) = 2ρc,0Ωm,0 ∞ 0 ξ 2h gm ( R 2 + χ 2 ) dχ (13) and Σ2h (< R) = 2 R 2 R 0 R Σ2h (R ) dR . (14)
In practice, we integrate up to 50Mpc in Σ2h (R), and verify that the function has converged at this point. ρc,0 and Ωm,0 are respectively the critical density and the matter density at current time, and ξ 2h gm is defined as

ξ 2h gm (r) = b h (M )ζ(r)ξm(r) (15)
where b h is the halo bias from Seljak & Warren (2004) and ζ(r) is the scale dependency of the bias as defined in equation B7 in Tinker et al. (2005):

ζ(r) = [1 + 1.17ξm(r)] 1.49 [1 + 0.69ξm(r)] 2.09 . ( 16
)
ξm(r) is the non-linear matter correlation function, computed as in Takahashi et al. (2012), and using the linear matter correlation function from Eisenstein & Hu (1998).

The two-halo halo term becomes predominant at scales of a few Mpc (see Fig A1), and has little influence as we measure the lensing signal only up to 1.8h -1 Mpc. We still include it to keep the model as generally applicable as possible.

The verification of the model predictions using numerical simulations are presented in appendix A1.

DATA

The source catalogues

We use three shear catalogues to compute the weak lensing signal: CFHTLenS, CS82 and DES-SV, covering a total area of 393 deg 2 . The total effective weighted source density is

n eff = 1 Ω ( w i ) 2
w 2 i = 7.8 galaxies/arcmin 2 , with Ω the total effective area (Heymans et al. 2012).

In the three catalogues, the estimation of photometric redshift, galaxy shape and stellar mass varies. To ensure that no catalogue gives biased results, we compare the results using each catalogue alone in appendix B1.

CFHT Stripe 82

The CFHT Stripe 82 Survey (CS82, Moraes et al. 2014) is a i-band imaging survey covering the 173 deg 2 of the SDSS Stripe 82 region, with a limiting magnitude of iAB ∼ 24.0 and seeing between 0.4 and 0.8 arcsec.

Each object has been attributed a mask flag indicating the quality of the photometry. Following Erben et al. (2013) we use all objects with mask 1. The remaining unmasked area is 129 deg 2 . We use photometric redshifts computed with BPZ (Benítez 2000) from ugriz SDSS photometry as described in Bundy et al. (2015). The initial catalogue contains all the sources with a photometric redshift z > 0, but we test alternative cuts: we examine the balance between having more sources or having more confident redshifts by either keeping only sources with the parameter odds > 0.5, or sources with 0.2 < z < 1.3. The odds parameter is a measure of the peakiness of the redshift probability distribution around the best redshift, which means that objects with odds close to unity have a peaked and unimodal p(z), and thus have a more reliable redshift. After testing the different cuts we choose to apply 0.2 < z < 1.3 for the rest of the work. We present the effect on our results of the different cuts in appendix B2 (see also Leauthaud et al. 2017).

The galaxy shapes were measured by the CS82 collaboration as described in Shan et al. (2017), using the lensfit method (Miller et al. 2013) following the procedure developed by the CFHTLenS collaboration (Erben et al. 2009(Erben et al. , 2013)). The shear systematics and calibration verifications also follow CFHTLenS as described in Heymans et al. (2012). Two calibration corrections need to be applied to the shear: one additive c2 and one multiplicative m (see section 3.4). The pipeline also provides an inverse variance weight for each source, and objects with weight = 0 are not used in the analysis. We select galaxies using the object classification obtained with lensfit, which separates galaxies, stars and badly fit objects (fitclass parameter, see details in Miller et al. 2013).

In summary, we keep objects with mask 1, weight > 0, 0.2 < z < 1.3 and fitclass = 0.The final source catalogue contains 3,791,129 sources with shear and photometric redshift measurement, that is an effective weighted source density of 6.7 galaxies/arcmin 2 (while n eff = 12.3 galaxies/arcmin 2 without any cuts). The completeness magnitude after the cuts is iAB ∼ 23.5.

CFHTLenS

The Canada-France-Hawaii Lensing Survey (CFHTLenS) is a 154 deg 2 (146.5 after masking) multi-band u*g'r'i'z' survey based on the Wide component of the Canada-France-Hawaii Telescope Legacy Survey (Heymans et al. 2012).

The photometric redshifts were measured over the whole survey using the BPZ code with PSF-matched photometry (Hildebrandt et al. 2012). We verify the effect on the lensing signal of the same source cuts as for the CS82 catalogue in appendix B2. Similarly to the CS82 catalogue, we finally use the cut 0.2 < z < 1.3 for the rest of the work.

The weak lensing data processing was performed with THELI (Erben et al. 2013), and the shape measurement with lensfit (Miller et al. 2013). The same calibration parameters as for the CS82 catalogue need to be applied to the shear.

In addition to the fitclass parameter, we cut our catalogue according to the star flag parameter to separate galaxies from stars. It classifies objects as galaxies or stars depending on their size, magnitude and colour information.

Similarly to the CS82 catalogue, we keep objects with mask 1, star flag = 0, fitclass = 0, 0.2 < z < 1.3 and weight > 0. The final catalogue contains 5,615,617 sources on an effective area after masking of 125 deg 2 , which gives an effective weighted source density of 10.7 galaxies/arcmin 2 (and n eff = 14.2 galaxies/arcmin 2 without cuts).

DES Science Verification Data

The Dark Energy Survey (DES, Flaugher 2005;Flaugher et al. 2015;Dark Energy Survey Collaboration et al. 2016) is an ongoing wide-field multi-band imaging survey, which will cover around 5000 deg 2 in the grizY filters, to a depth of iAB < 24. In this analysis we use data taken during the Science Verification period (SVA1), which covers an area of 139 deg 2 after masking, with a depth comparable to the full survey 1 . We use the shear catalogue obtained using the NGMix 2 algorithm (Sheldon 2014) and described in Jarvis et al. (2016), and photometric redshifts from Bonnett et al. (2016). After applying the cuts sva1 flag = 0 (good galaxy) and ngmix flag = 0 (good shape measurement), the catalogue contains 3,446,533 galaxies, or a weighted density of 6.8 galaxies/arcmin 2 (n eff = 17.2 galaxies/arcmin 2 without cuts). The completeness magnitude after the cuts is iAB ∼ 23.

The multiplicative calibration factor, which needs to be applied to the shear as described in section 3.4, is given in the NGMix catalogue as the sensitivity estimate sens avg and is related to the usual factor by sens avg = 1 + m.

The DES-SV data release contains a second catalogue with galaxy shape measurement, using an alternative algorithm Im3shape 3 . We made some lensing measurements using both catalogues and found no significant difference, and thus chose to use the NGMix catalogue.

The redshift distributions of the sources in the three catalogues are presented in figure 4.

1 https://des.ncsa.illinois.edu/releases/sva1 2 https://github.com/esheldon/ngmix 3 https://bitbucket.org/joezuntz/im3shape-git

Lens selection

Our sample consists of the overlap of the satellite galaxies from the redMaPPer catalogue with the CFHTLenS, CS82 and DES-SV fields.

The redMaPPer cluster catalogue is derived from photometric galaxy samples using the red-sequence Matched-filter Probabilistic Percolation cluster finder (Rykoff et al. 2014(Rykoff et al. , 2016)). For each galaxy the algorithm estimates the membership probability according to its color, position to the centre of the cluster and luminosity. To reduce contamination by line-of-sight galaxies, we select satellites with membership probability Pmem > 0.8; as shown in Zu et al. (2016), cutting at this value should eliminate most of the contamination. In our analysis we use the catalogues extracted from the 5-band (ugriz) photometry of the SDSS Data Release 8 (DR8, Aihara et al. 2011), described in Rykoff et al. (2014) and covering the CS82 and CFHTLenS fields, and from the griz photometry of the Dark Energy Survey Science Verification Data described in Rykoff et al. (2016).

For the lenses we use the redMaPPer redshifts, as they are more robust than photometric redshifts estimated using other methods (see Bundy et al. 2015, section 5). In order to have the same redshift range in the three samples, we keep only galaxies with 0.2 < z < 0.55. The redshift distribution of the satellite galaxies in the two radial bins is shown in the left panel of figure 5. In this redshift range, the redMaPPer cluster catalogue contains 289 clusters in the CFHTLenS footprint, 491 in CS82 and 349 in DES-SV. We note that the whole CFHTLenS field does not overlap with the redMaPPer cluster catalogue: there are around 30 deg 2 not covered, which gives an effective unmasked area of ∼ 115 deg 2 for CFHTLenS. The cluster density is 2.5 deg -2 in the CFHTLenS and DES-SV fields, and 3.8 deg -2 in the CS82 field. We suspect that this difference is due to the presence of a large structure in the CS82 field, as the cluster distribution appears to have a significant excess at z ∼ 0.4 compared to the other fields.

We select satellite galaxies according to their projected cluster-centric distance, in the two following bins: the satellites in the inner part of the clusters, with Rs ∈ [0.1; 0.55]h -1 Mpc, and the satellites in the outer part of the clusters, with Rs ∈ [0.55; 1]h -1 Mpc. Indeed, the redMaP-Per clusters used in the analysis have a richness between 20 and 180, which, according to the mass-richness relation from Rykoff et al. (2012) corresponds to a virial mass between 10 14 and 10 15 , which gives a virial radius in the range 0.6 -1.4h -1 Mpc.

A possible alternative binning would be in fractions of Rvir, ensuring that all the satellites are taken for each cluster. This choice would in addition avoid the selection effect discussed in section 4.3. However, we do not have precise estimations of the virial radius or mass for each cluster, and we should therefore use some mass-richness relation. These relations always show an important scatter (see for example Rozo & Rykoff 2014;Rozo et al. 2015), which would add noise to our bining.

We further split our lens sample in three stellar mass mean stellar mass and mean redshift are shown in table 1. We describe in the next section how the stellar masses were obtained.

Stellar mass

The CFHTLenS catalogue provides stellar masses for the galaxies, computed as described in Velander et al. (2014), by fitting spectral energy distribution (SED) templates, using the software le Phare (Arnouts et al. 1999;Ilbert et al. 2006) following the method of Ilbert et al. (2010). For CS82 and DES-SV we estimate the stellar masses using the software le Phare through the gazpar4 web service, and the available photometry for each survey:

• for the CS82 catalogue we use the ugriz SDSS Stripe 82 Coadd photometry and YJHK photometry from the UK Infrared Deep Sky Survey (UKIDSS) Large Area Survey;

• for DES-SV we use the griz photometry.

For both catalogues the fitted SED templates are built from the stellar population synthesis (SPS) package from Bruzual & Charlot (2003), with a Chabrier (2003) initial mass function (IMF), and the star formation history is described as a decreasing exponential function e -t/τ with 0.1 τ 30Gyr. Dust extinction was applied following two laws, a starbust (Calzetti et al. 2000) and an intermediate slope (λ 0.9 ) law, and reddening excess E(B -V ) ranging from 0 to 0.7.

The stellar mass distribution of the satellite galaxies in the two radial bins is shown in the right panel of figure 5. The distributions for the satellites in the inner part of clusters and in the outer part are quite similar. Moreover, assuming that the stellar masses of galaxies do not vary much on average during infall (due to the star formation quenching and the relative insensitivity of the stellar component to stripping, see Smith et al. 2016), we consider that any measured mass difference between the satellites in the two radial bins is due to the evolution of the subhaloes during accretion. We study the validity of the assumption that the stellar mass remains constant during accretion in section 5.1.

Computation of the lensing signal

We measure ∆Σ by stacking lens-source pairs in 13 radial bins from 0 to 1.8h -1 Mpc. For each sample of lenses, ∆Σ(R) is estimated using

∆Σ(R) = ls w ls γ ls t Σcrit(z l , zs) ls w ls , (17)
where γ ls t is the tangential shear, w ls = ws/ (Σcrit (z l , zs)) 2 , and ws is an inverse variance weight factor associated to each source galaxy and introduced to account for shape measurement error and intrinsic scatter in galaxy ellipticity (Heymans et al. 2012). The sum is calculated over all the lenssource pairs.

The multiplicative calibration factor m needs to be taken into account in a statistical way (Miller et al. 2013), we apply it to the mean shear measurement using the correction proposed in Velander et al. (2014) and Hudson et al. (2015):

1 + K(z l ) = ls w ls (1 + m) ls w ls , (18)
which gives a calibrated lensing signal:

∆Σ cal (R) = ∆Σ(R) 1 + K(z l ) = ls w ls γ ls t Σcrit ls w ls (1 + m) . ( 19
)
We only use lens-source pairs with zsource > z lens + z err lens + z err source to ensure that no sources are at lower redshift than the lens. The error bars on the lensing signal are obtained with a block bootstrap on the data: the field is divided in blocks, and the lensing signal is measured on the resampled blocks to estimate the variance of the measurement.

We compute the lensing signal for each lens sample, using a modified version of the athena5 software, a 2d-tree code estimating second-order correlation functions from input galaxy catalogues.

RESULTS

Fitting procedure

We fit the lensing signal with the model described in section 2.2, which has two free parameters: the mean satellite mass Msat and the mean host mass M host . We obtain the best-fit parameters by maximizing the likelihood L, and to obtain the intervals of confidence we use a Markov chain Monte Carlo (MCMC) method using emcee (Foreman-Mackey et al. 2013) which is a Python implementation of an affine invariant MCMC ensemble sampler.

The likelihood L is expressed as:

L = 1 2π|C| exp(- 1 2 13 i=1 (x obs i -x mod i ) 2 σ 2 i ) (20)
where x obs i are the measurements and x mod i the model predictions in the 13 radial bins. The σi are the diagonal terms of the covariance matrix C computed with the block bootstrap. We use only the diagonal as the full matrix is quite noisy. The average amplitude of a non-diagonal term of the covariance matrix is around 10% of a diagonal term, which is consistent with the noise in matrices computed with bootstrap as shown in Viola et al. (2015). The use of only diagonal terms might result in overestimating the quality of the fit, which could underestimate the size of the errorbars on the mass measurements. We assume flat and broad priors: for the satellite mass we choose log(Msat/h -1 M ) ∈ [9.5; 13.5] and for the cluster mass log(M host /h -1 M ) ∈ [13.5; 16].

Dependance of the SHMR on the projected cluster-centric distance

We plot the lensing signal for each of the lens samples in figure 6, with columns from left to right showing the bins in log stellar mass [10;10.5], [10.5;11], [11;11.5]. The top line shows the signal for satellites in the inner parts of clusters (Rs ∈ [0.1; 0.55]h -1 Mpc), and the bottom line for satellites in the outer part (Rs ∈ [0.55; 1]h -1 Mpc). In addition, for each sample we plot the best-fit model: the blue continuous line is the full model, and the blue shaded area is the 68% credible interval.

The maximum likelihood results with the 68% credible intervals are presented in table 1, with the average stellar mass and redshift for each lens sample. In addition, we present the joint 2-dimensional and marginalized 1-dimensional posterior probability distributions for the two parameters Msat and M host for each of the six samples in figure 7.

For the low mass satellites, the signal-to-noise is too low to obtain strong constraints. For the intermediate and high mass samples, we find that at a given stellar mass, the dark matter halo mass is shifted toward lower mass for satellites in the inner part of clusters, compared to the satellites in the outer part. A summary of the measured satellite dark matter mass as a function of the mean stellar mass in each sample is shown in figure 8.

We expect the satellites in the outer parts of the clusters to be similar to field/central galaxies. As the distance to the centre of the cluster correlates with the accretion redshift, most of them have only started recently their accretion process, and have not yet been submitted to strong tidal forces from the cluster for a long time. On the contrary, the satellites close to centre of the cluster have undergone the influence of the cluster much more strongly and for a longer time, and have thus been stripped of a significant part of their dark matter. This hypothesis is in agreement with our measurements. We plot on figure 8 Assuming that our measurements do reflect the dark matter stripping scenario, we quantify the effect for the two higher mass bins by computing stripping factors which represent the amount of stripped dark matter, defined as

τstrip = 1 - M inner sat M outer sat . ( 21
)
For the intermediate mass bin, we obtain τ inter strip = 0.83 +0.15 -0.69 , and for the high mass τ high strip = 0.69 +0.26 -1.14 . Our error bars do not allow us to have strong constrains, but we still compute a theoretical value for the stripping factor as a comparison. Using the equation from Gao et al. (2004) for the infall mass, the stripping rate can be expressed as:

τstrip(Rs) = 1 - M sub (Rs) M infall = 1 -0.65 Rs Rvir 2/3 . ( 22
)
We take for the distance cluster centre-satellite the median value in the inner bin sample Rs = 0.28h -1 Mpc and for the host virial radius the value corresponding to the median between our best-fit host masses Rvir = 1.16h -1 Mpc, and obtain a stripping rate τstrip = 0.75. We compute the relative likelihood of obtaining a value of τstrip equal to 0 (no stripping) or 0.75 (theoretical value for stripping) using our posterior probability distributions as

L(τ ) = L M inner sat = (1 -τ )M outer sat = L Minner = (1 -τ )M × L Mouter = M dM (23)
assuming the independence of our samples. We obtain for the high mass sample P (τ = 0) = 0.43 and P (τ = 0.75) = 0.84, and for the intermediate mass sample P (τ = 0) = 0.28 and P (τ = 0.75) = 0.59, which shows that our results clearly favour the stripping scenario over the no stripping one. In addition we perform a Kolmogorov-Smirnov twosample test to verify if the posterior probability distributions for M inner sat and M outer sat are statistically different, and reject the null hypothesis (no stripping) since the p-value is below 1%.

Finally, we compute the ratio subhalo mass over stellar mass for the subhaloes in the intermediate and high stellar mass bins, and plot it in figure 9. We also plot the results from Sifón et al. (2015) and Li et al. (2016) and find them to be broadly consistent with our results.

Mass of the redMaPPer clusters

We now look at our measurement of the mass of the host redMaPPer clusters. In each of the cluster-centric bins, we find that the best-fit host halo mass has a consistent value across the different stellar mass bins. However, the average host mass is lower for satellites in the inner part of clusters (log M host ∼ 14.20) than for satellites in the outskirts (log M host ∼ 14.30). This is due to a selection effect: larger host clusters have more satellites at a distance ranging between 0.55 and 1h -1 Mpc, and are thus more represented in the outer part satellite bin, pushing the mean host halo mass towards higher values.

M * range [M ] Rs[h -1 Mpc] N lenses < log M * /M > < z l > log(Msat/h -1 M ) log(M host /h -1 M ) 10-
R[h -1 Mpc] 0.2 0.6 1.0 1.4 R[h -1 Mpc] Rsat = [0.1 : 0.55]h -1 Mpc Rsat = [0.55 : 1]h -1 Mpc
We verify this assumption by computing the mean host richness in each satellite bin. For galaxies in the inner part of clusters, we find that the mean richness of their host varies between 43 and 45 (depending on the stellar mass bin), while for galaxies in the outer part the mean richness is between 58 and 63.

In addition, we compute from the richness the halo mass using the mass-richness relation from Rykoff et al. (2012): we find that for the satellites in the centre the mean host halo mass in each of the different samples computed from the mean richness is in the range log(M host /h -1 M ) ∈ [14.21; 14.22], and for the satellites in the outskirts the mean host mass is in log(M host /h -1 M ) ∈ [14.31; 14.34], which agrees very well with our mass measurements given in table 1.

DISCUSSION

Stellar mass evolution during infall

In this study, our observations suggest that for a given stellar mass, satellite galaxies have a more massive dark matter halo if they are in the outskirts of their host cluster than if they are close to its centre. We explain these observations with the following scenario: during its accretion to a cluster, a galaxy has its star formation quenched while its dark matter halo is stripped by the tidal forces of the host. In this section, we consider whether some alternative infall

1 1 1 2 1 3 log(M sat ) 1 4 . 2 1 4 . 3 1 4 . 4 log(M host ) 1 4 . 2 1 4 . 3 1 4 . 4 log(M host )
(a) 10 < log M * < 10.5, 0.1 < Rs < 0.55 11 -11.5]. In the 1d distributions, the blue lines are the maximum likelihood parameters, and the black dashed lines are the 68% credible intervals. In the 2d distributions, the blue crosses are the maximum likelihood solutions, and the black contours are the 68% and 95% joint credible regions. Masses are in units of h -1 M and distances in h -1 Mpc.

1 1 . 2 1 2 . 0 1 2 . 8 log(M sat ) 1 4 . 1 1 4 . 2 1 4 . 3 1 4 . 4 log(M host ) 1 4 . 1 1 4 . 2 1 4 . 3 1 4 . 4 log(M host ) (b) 10.5 < log M * < 11, 0.1 < Rs < 0.55 1 1 . 2 1 2 . 0 1 2 . 8 log(M sat ) 1 4 . 1 0 1 4 . 2 5 1 4 . 4 0 log(M host ) 1 4 . 1 0 1 4 . 2 5 1 4 . 4 0 log(M host ) (c) 11 < log M * < 11.5, 0.1 < Rs < 0.55 1 1 1 2 1 3 log(M sat ) 1 4 . 2 1 4 . 3 1 4 . 4 log(M host ) 1 4 . 2 1 4 . 3 1 4 . 4 log(M host ) (d) 10 < log M * < 10.5, 0.55 < Rs < 1 1 1 . 2 1 2 . 0 1 2 . 8 log(M sat ) 1 4 . 1 1 4 . 2 1 4 . 3 1 4 . 4 log(M host ) 1 4 . 1 1 4 . 2 1 4 . 3 1 4 . 4 log(M host ) (e) 10.5 < log M * < 11, 0.55 < Rs < 1 1 1 . 2 1 2 . 0 1 2 . 8 log(M sat )
processes could be considered, such as evolution of the satellite stellar mass.

A first possibility is that if the stellar component of the galaxy gets stripped along the dark matter halo, the stellar mass and the dark matter mass would decrease during infall. However, as tidal forces are expected to remove matter from a galaxy starting from the outside towards the centre, stellar matter is expected to remain undisturbed much longer than dark matter. Using hydrodynamical simulations, Smith et al. (2016) studied how the stellar component of a galaxy is affected by stripping compared to the dark matter halo. Using their equation 2, we compute the fraction of remaining stellar matter after stripping, assuming that dark matter has been stripped by a factor of τstrip. Using the best-fit values obtained in section 4.2, we find that for the high mass sample 99% of the stellar component remains, and 91% for the intermediate mass sample. The possibly removed amount of stellar matter is thus very small and would not change our conclusions. Even if the equation from Smith et al. (2016) underestimates the amount of stellar stripping, considering that stars are partly stripped during accretion only strengthens our conclusions: if a satellite galaxy has been stripped of dM * stars during infall, its progenitor was a galaxy of stellar mass M * + dM * which would have a more massive dark matter halo than a progenitor of stellar mass dM * (assuming that the relation between stellar mass and dark matter mass is monotonic, as in Moster, Naab & White 2013).

We now focus on the more interesting topic of star formation during accretion. Indeed, what we explain as a stripping of the dark matter at constant stellar mass might also be an increase of the stellar mass at fixed (or less shifted) dark matter mass. We verify here how this scenario would fit with our observations, by estimating how much stars could have formed during the infall process for our sample of galaxies in the inner part of clusters. First, using the equation from Giocoli, Tormen & van den Bosch (2008) we estimate the time since accretion t inf and the redshift at accretion Figure 9. Subhalo to stellar mass ratio as a function of the cluster-centric distance, calculated for the satellites in the intermediate (green) and high (blue) stellar mass bins. We do not plot the low stellar mass bin as it is not very constrained. We plot as a comparison the results from Li et al. (2016) and Sifón et al. (2015).

z inf :

t inf = ln 0.65 Rs Rvir 2/3 t dyn (z), (24)
where (van den Bosch, Tormen & Giocoli 2005) . Same as figure 8 zoomed on the two higher stellar mass bins. We show two possible scenario to explain our observations: the black arrows would represent a stripping of the dark matter at fixed stellar mass, while the green arrows represent a coevolution of the stellar and dark matter, with star formation rates from Buat et al. (2008).

t dyn (z) = 2 ∆vir(z) ∆vir(0) -0.5 H(0) H(z) , (25)
where ∆vir(z) is the virial overdensity at redshift z. Using the same values for Rs and Rvir as for the stripping factor in section 4.2, we find t inf = 2.03Gyr and z inf = 0.62. We can then obtain the amount of star formation that occurred during the infall time, assuming that there is no quenching due to the cluster environment at all. We use star formation rates computed in Buat et al. (2008), that is, for the three stellar mass bins, from low to high: SF R(M 1 * , z inf ) ∼ 5M /yr, SF R(M 2 * , z inf ) ∼ 11M /yr and SF R(M 3 * , z inf ) ∼ 23M /yr. We can then compute the stellar mass at infall for the three stellar mass bins:

M * ,inf = M * -t inf * SF R(z inf ), (26)
which gives log(M 1 * ,inf /M ) = 10.00, log(M 2 * ,inf /M ) = 10.50 and log(M 3 * ,inf /M ) = 11.00. Using the relation from Moster, Naab & White (2013), we infer the corresponding dark matter halo mass before accretion, and show in figure 10 the corresponding evolution scenario compared to the scenario with constant stellar mass. In the case where we consider an evolution in the stellar mass (green arrows), the shift in dark matter mass is less important than for the case with fixed stellar mass (black arrows), but still consistent with some stripping of the dark matter halo.

It is important to note that we use here star formation rates that correspond to star-forming galaxies (Buat et al. 2008) and that many studies show that cluster galaxies are at least partly quenched in this dense environment. The two scenarios we present are therefore only the two extreme limits (with no star-formation and with no quenching), and a more detailed study of the coevolution of stellar and dark matter should be carried out.

Systematic error

We perform different tests to verify the robustness of our results. The analyses are described in appendix A and B, and we summarize here the conclusions.

To verify if no systematic bias is introduced in the lensing profiles by the differences in the source catalogues (shape measurements, redshift estimations, etc.) we perform the same measurement using each catalogue individually. The subhalo masses agree within one sigma between the different lensing catalogues.

We then examine the influence of the redshift quality in the CS82 and CFHTLenS catalogues, by testing different quality cuts and comparing the results. We chose for our fiducial source catalogue all galaxies with 0.2 < z < 1.3.

Finally, using the BigMultiDark simulations (Giocoli et al. 2016;Klypin et al. 2016, described in appendix A1), we test the influence of line-of-sight projections on our measurements. In the redMaPPer cluster catalogue we cut satellites according to their cluster membership probability to improve the purity of the satellite catalogue, but some contamination by non-member galaxies is still possible. We quantify the effect of this on the lensing signal using simulations, and find that contamination can decrease the amplitude of the lensing profile by ∼ 1%.

CONCLUSION

Many numerical simulations and observations suggest that dark matter haloes of satellite galaxies are subjected to the tidal stripping induced by the gravitational potential of their host cluster. In this work, we study this effect by measuring the mass of satellite galaxies using galaxy-galaxy lensing: we compare the mass of satellites located in the outer part of clusters (which have just started their accretion process and are therefore less perturbed by the cluster), with satellites in the inner parts of clusters (which, in average, have been subject to the influence of the cluster for a longer time). We further divide our satellites in stellar mass bins, as the stellar mass of a galaxy appears to be a good tracer of the mass of the galaxy before infall.

We find good agreement with Moster, Naab & White (2013) for the mass of satellite galaxies in the outer radial bin and a suggested shift to smaller subhalo masses in the inner radial bin, in agreement with the dark matter stripping scenario. For the intermediate and high stellar mass bins, we find stripping factors of τ inter strip = 83 +15 -69 % and τ high strip = 69 +26 -114 %. Using the posterior probability distributions, we find that the theoretical stripping scenario τstrip = 0.75 is favoured with respect to the no stripping scenario τstrip = 0. We find for the intermediate mass sample P (τstrip = 0.75) = 0.59 and P (τstrip = 0) = 0.28, and for the high mass sample P (τstrip = 0.75) = 0.84 and P (τstrip = 0) = 0.43.

While our estimations of the subhalo masses appear to be consistent with a stripping of the dark matter by the gravitational potential of the host clusters, some other effects can have an influence on the results. We consider two possible evolutions in stellar mass, stripping or star-formation, and find that while stellar stripping can be neglected in the process, stellar formation could partly explain our observations. To disentangle the contributions of dark matter stripping and star formation, a study of the coevolution of dark and stellar matter during accretion should be carried out.

To improve the measurements of the subhalo mass, different solutions can be considered. Future lensing surveys will cover thousands of square degrees instead of hundreds, decreasing statistical errors. In addition, spectroscopically confirmed clusters can be studied, to remove contamination and have more precise redshifts. The improvement of the lensing measurement of subhaloes will allow not only to constrain their mass evolution more precisely, but also to study their density profiles and measure the tidal truncation radius. (2.20 ± 0.20)10 14 h -1 M . The lack of signal at small scales (R < 0.3h -1 Mpc) is due to the resolution of the simulation.

10 -1 10 0 10 1 R[h -1 Mpc] 10 0 10 1 ∆Σ[M hpc -2 ]
particle) is 2.35×10 10 h -1 M . The effect of this limited mass resolution has an impact on the lensing measurements, and induce a damping of the signal at small scales.

For our measurements we use 15 lightcone realizations. Each realization corresponds to a rectangular field covering 55 deg 2 .

A2 Model verification

To verify the predictions of the model, we use the BigMul-tiDark simulated data. We compute the lensing signal for all the haloes with 1 × 10 12 < M h /h -1 M < 2 × 10 12 for the 15 lightcones. We fit the model to the mean lensing signal over the 15 computed, to obtain a best-fit halo mass of M 1h = (1.42 ± 0.07)10 12 h -1 M and a best-fit host mass of M host = (2.20 ± 0.20)10 14 h -1 M , which is consistent with the input. In this case, as our lens sample is a mixture of central and satellite galaxies, the satellite fraction is not equal to one, and according to the data we fix fsat = 0.09. The lensing signal and best-fit model are shown in figure A1.

The model gives a good description of the simulated signal, except at small scales (R < 0.3h -1 Mpc), where the limited resolution of the simulation damps the signal.

A3 Contamination by non-satellite galaxies

Even if we select only the satellite galaxies in the redMaPPer catalogue with a probability of membership higher than 0.8, our samples can still contain line-of-sight galaxies which are not member of the clusters. We test their influence on the lensing signal using the BigMultiDark simulations.

We use as hosts the haloes with M h > 10 14 h -1 M in the redshift range 0.2 < z host < 0.55, and select the satellites in the same cluster-centric distance bins as in our measurement: 0.1h -1 Mpc < Rs < 0.55h -1 Mpc and 0.55h -1 Mpc < Rs < 1h -1 Mpc. We make the following red- A1. Best-fit masses for a sample of pure satellite galaxies and a sample contaminated by line-of-sight effects. The distances are expressed in h -1 Mpc and the masses in h -1 M .

10 -1 10 0 R[h -1 Mpc] 0 10 20 30 40 50 60 ∆Σ[M hpc -2 ] R s = [0.1 : 0.55]h -1 Mpc 10 -1 10 0 10 1 R[h -1 Mpc] R s = [0.55 : 1]h -
shift selections to obtain one sample with only real satellites, and one sample contaminated by galaxies in the line-of-sight:

• the pure satellite sample: we select galaxies with |z galz host | 0.0007. In the redshift range 0.2 -0.55, this corresponds to |d gal -d host | 1.5h -1 Mpc ∼ R host vir . • the contaminated sample: here we have |z gal -z host | 0.05. This gives a contamination rate of 37% in the inner radial bin and 65% in the outer radial bin.

We show in figure A2 the lensing signal computed in both radial bins, for the pure and the contaminated samples. In both radial bins, when contaminants are added, the amplitude of the signal decreases at small scales, inducing a drop in the subhalo mass: for the satellites in the inner part (respectively the outer part) of clusters, the best-fit satellite mass goes from 13.02 ± 0.03 (respectively 13.00 ± 0.04) for the pure sample to 12.88 ± 0.03 (respectively 12.86 ± 0.04) with contaminants. All the best-fit masses are summarized in table A1. To take into account this drop, we should increase our measurements by 1%, which does not change our conclusions.

In addition to line-of-sight projections, the miscentring of clusters could be another source of contaminations (George et al. 2012): if the galaxy we use as centre of the cluster is not its true centre, we can identify outer part satellites as inner ones and vice-versa. To verify the importance of this effect on the signal, we make the measurements again, selecting this time only clusters with a probability of having their true centre identified higher than 0.5. This selection has no significant impact on the signal, other than increasing the error bars due to the decrease in the number of lenses. We measured the haloes of satellite galaxies in three bins in stellar mass, and in two ranges in distance to the centre of the cluster, using shear catalogues from the CS82, CFHTLenS and DES-SV surveys. Each of these surveys have different photometric measurements, as well as different estimators of the photo-z, stellar mass and shape. To ensure that no bias is introduced by either of the catalogues, we repeat the same measurements using each of the catalogues separately. We compare the obtained masses in figure B1, and although the uncertainties are high, the different measurements of the satellite mass appear to be consistent.

B2 Redshift accuracy

Different cuts are possible in the source catalogues to have photometric redshifts as clean as possible. As shown in Heymans et al. ( 2012), keeping only galaxies with 0.2 < z < 1.3 ensures a relatively accurate photometric redshift, with a typical uncertainty of σz ∼ 0.04(1 + z). In addition, cutting objects with ODDS < 0.5 secures even further the accuracy of redshifts. We test both cuts on the CS82 and CFHTLenS . Top panel: best-fit results for the satellite masses computed using only the CFHTLenS source catalogue. The blue symbols are the masses for satellites from the inner part of the clusters, and the red symbols for satellites from the outer part of the clusters. We compute the masses using three different source catalogues: all the sources, sources selected with ODDS > 0.5 and sources with 0.2 < z < 1.3. Bottom panel: same but for the CS82 source catalogue.

catalogues: for CS82 the initial catalogue contains 4,381,917 objects, and keeps 87% with 0.2 < z < 1.3, and 64% with ODDS > 0.5. For CFHTLenS, the uncut catalogue contains 7,511,368 objects, going down to 75% with 0.2 < z < 1.3 and 89% with ODDS > 0.5. We show in figure B2 the best-fit subhalo masses computed using the different cuts. We find that the different results are very similar, and we choose to use the cut 0.2 < z < 1.3 which appears to give the best contraints.

Figure 2 .

 2 Figure 2. Satellite spatial distribution in three bins in stellar mass: 10 < log(M * /M ) < 10.5, 10.5 < log(M * /M ) < 11 and 11 < log(M * /M ) < 11.5. The dashed and dotted curves are the best-fit second order polynomials in each cluster-centric bin.

Figure 4 .

 4 Figure 4. Redshift distribution of the source galaxies in the CS82, CFHTLenS and DES-SV surveys.

Figure 5 .

 5 Figure5. Redshift (left panel) and stellar mass (right panel) distribution of the satellite galaxies in the combined catalogue, for the galaxies in the inner radial bin (Rs ∈ [0.1; 0.55]h -1 Mpc) in blue and the outer radial bin (Rs ∈ [0.55; 1]h -1 Mpc) in red. The distributions are normalized so their integrals sum to one, and the vertical lines in the right panel represent the limits of the stellar mass bins.

  the stellar-tohalo mass relation for field/central galaxies measured from N-body simulations in Moster, Naab & White (2013). Their relation is indeed consistent with the one we measure for satellites in the outskirts of clusters. In addition, the shift of the SHMR towards lower subhalo mass is consistent with the results obtained in Rodríguez-Puebla, Drory & Avila-Reese (2012); Rodríguez-Puebla, Avila-Reese & Drory (2013) using extended abundance matching technique.

  (M * /M ) < 10.5 10.5 < log(M * /M ) < 11 11 < log(M * /M ) < 11

Figure 6 .

 6 Figure 6. Lensing by the redMaPPer satellite galaxies. The top line shows the lensing signal for the satellites in the inner part of the cluster, with Rs ∈ [0.1; 0.55]h -1 Mpc, and the bottom one for the satellites in the outer part, with Rs ∈ [0.55; 1]h -1 Mpc. The left column is for satellites with log(M * /M ) ∈ [10; 10.5], the middle column for log(M * /M ) ∈ [10.5; 11] and the right one for log(M * /M ) ∈ [11; 11.5]. The black dots are the data points with bootstrap errors, the blue curve is the best-fit model, and the blue surface shows the 68% confidence interval.

Figure 7 .

 7 Figure 7. Joint 2-dimensional and marginalized 1-dimensional posterior probability distributions for our two parameters Msat and M host for each of the six samples: the top line presents the bins with Rs ∈ [0.1; 0.55]h -1 Mpc, and the bottom line Rs ∈ [0.55; 1]h -1 Mpc. The columns, from left to right represent the stellar mass bins log(M * /M ) ∈ [10 -10.5], [10.5 -11], [11 -11.5]. In the 1d distributions, the blue lines are the maximum likelihood parameters, and the black dashed lines are the 68% credible intervals. In the 2d distributions, the blue crosses are the maximum likelihood solutions, and the black contours are the 68% and 95% joint credible regions. Masses are in units of h -1 M and distances in h -1 Mpc.

Figure 8 .

 8 Figure8. Stellar-to-halo mass relation measured for the satellite galaxies in our sample, in blue for the galaxies in the inner part of clusters (Rs ∈ [0.1; 0.55]h -1 Mpc) and in red for the galaxies in the outer part of clusters (Rs ∈ [0.55; 1.]h -1 Mpc). The black line is the SHMR for field/central galaxies at z = 0.35, computed from simulations inMoster, Naab & White (2013).

  Figure10. Same as figure 8 zoomed on the two higher stellar mass bins. We show two possible scenario to explain our observations: the black arrows would represent a stripping of the dark matter at fixed stellar mass, while the green arrows represent a coevolution of the stellar and dark matter, with star formation rates fromBuat et al. (2008).

Figure B1 .

 B1 Figure B1. Comparison of the mass measurements using the CS82 and CFHTLenS and DES-SV catalogues. The left column shows the satellite mass and the right panel the host mass. We plot the x = y line to guide the eye.

  Figure B2. Top panel: best-fit results for the satellite masses computed using only the CFHTLenS source catalogue. The blue symbols are the masses for satellites from the inner part of the clusters, and the red symbols for satellites from the outer part of the clusters. We compute the masses using three different source catalogues: all the sources, sources selected with ODDS > 0.5 and sources with 0.2 < z < 1.3. Bottom panel: same but for the CS82 source catalogue.

Table 1 .

 1 Description of the lens samples and best-fit parameters.

	10.5 10.5-11 11-11.5	0.1-0.55 0.55-1 0.1-0.55 0.55-1 0.1-0.55 0.55-1	48891 1292 6935 1836 2126 677	10.31 10.32 10.73 10.73 11.17 11.17	0.35 0.37 0.37 0.36 0.38 0.36	11.55 +0.26 -1.05 11.18 +0.63 -1.26 11.76 +0.23 -0.75 12.54 +0.16 -0.48 12.44 +0.19 -0.64 12.95 +0.13 -0.48	14.17 +0.02 -0.02 14.31 +0.04 -0.05 14.17 +0.02 -0.02 14.29 +0.05 -0.05 14.20 +0.03 -0.03 14.34 +0.08 -0.07

  Lensing signal produced by haloes with M 1h ∈ [1 × 10 12 ; 2 × 10 12 ]h -1 M , from the BigMultiDark simulation. The best-fit model gives a best-fit halo mass of M 1h = (1.42 ± 0.07)10 12 h -1 M and a best-fit host mass of M host =

	Dsigma 2h Dsigma host Dsigma 1h Dsigma total
	Figure A1.

http://gazpar.lam.fr/index

http://www.cosmostat.org/software/athena/

ACKNOWLEDGEMENTS

Based

on observations obtained with MegaPrime/MegaCam, a joint project of CFHT and CEA/DAPNIA, at the Canada-France-Hawaii Telescope (CFHT), which is operated by the National Research Council (NRC) of Canada, the Institut National des Science de l'Univers of the Centre National de la Recherche Scientifique (CNRS) of France, and the University of Hawaii. The Brazilian partnership on CFHT is managed by the Laboratório Nacional de Astrofísica (LNA). We thank the support of the Laboratório Interinstitucional de e-Astronomia (LIneA). We thank the CFHTLenS team. This work was granted access to the HPC resources of Aix-Marseille Université financed by the project Equip@Meso (ANR-10-EQPX-29-01) of the program "Investissements d'Avenir" supervised by the Agence Nationale pour la Recherche. This work is based in part on data products produced by GAZPAR located at the Laboratoire d'Astrophysique de Marseille. M. L thanks CNRS for financial support. MM is partially supported by CNPq and FAPERJ. Fora Temer. This work is supported by the Deutsche Forschungsgemeinschaft in the framework of the TR33 'The Dark Universe'. Carlo Giocoli thanks CNES and acknowledges support from the Italian Ministry for Education, University and Research (MIUR) through the SIR individual grant SIMCODE, project number RBSI14P4IH. HH was supported by an Emmy Noether grant of the Deutsche Forschungsgemeinschaft (Grant no.: Hi 1495/2-1).

APPENDIX A: NUMERICAL SIMULATIONS A1 Description of the simulations

We use the BigMultiDark Planck cosmological simulation (BigMDPL, Klypin et al. 2016, www.multidark.org) which contains 3840 3 particles in a box of comoving side 2.5h -1 Gpc. Haloes and subhaloes are identified from density peaks using the Bound Density Maximum (BDM) halo finder (Klypin & Holtzman 1997;Riebe et al. 2013). The construction of the light-cones is described in Giocoli et al. (2016): the particles from 24 snapshots of the simulation were projected onto 24 lens planes distributed along the line-of-sight up to redshift z = 2.3. To simulate lensing, they trace light-rays through the successive lens planes of the light-cones using the ray-tracing code Glamer (Metcalf & Petkova 2014;Petkova, Metcalf & Giocoli 2014).

The simulation mass resolution (ie mass of the smallest


Dalia S Youssef

Abdel Massih

Abdul-Hamid Soubra

Jacques Harb

Mohamed Rouainia

DYNAMIC SLOPE STABILITY ANALYSIS BY A RELIABILITY-BASED APPROACH

The analysis of earth slopes situated in seismic areas has been extensively investigated in literature using deterministic approaches where average values of the input parameters are used and a global safety factor or a permanent displacement at the toe of the slope is calculated. In these approaches, a decoupled analysis is usually performed in which a constant critical seismic acceleration is first calculated based on a pseudo-static representation of earthquake effects. Then, the permanent soil displacement is computed by integration of the real acceleration record above the critical acceleration. A reliability-based approach for the seismic slope analysis is introduced in this paper. This approach is more rational than the traditional deterministic methods since it enables to take into account the inherent uncertainty of the input variables. Furthermore, the deterministic model used in the reliability analysis is based on a rigorous coupled analysis that simultaneously captures the fully non linear response of the soil and history of the real acceleration record. This model is based on numerical simulations using the dynamic option of the finite difference FLAC 3D software. The acceleration time history records used in this analysis is the one registered at the Lebanese Geophysical Center of the National Council for Scientific Research. The performance function used in the reliability analysis is defined with respect to the horizontal permanent displacement at the toe of the slope. The random variables considered in the analysis are the cohesion c and the shear modulus G of the soil since it was shown that these parameters have the most influence on slope displacements. The Stochastic Response Surface Methodology (SRSM) is utilized for the assessment of the probability distribution of the horizontal permanent displacement at the toe of the slope. The results have shown that the mean value of the permanent displacement is highly influenced by the coefficient of variation of the cohesion, the greater the scatter in c the higher the horizontal permanent displacement. Based on the probability distribution of the permanent displacement, the failure probability with respect to the exceedance of an allowable maximal permanent displacement was evaluated and discussed. At the end of the paper, a case study of a typical Lebanese slope is presented and analyzed to illustrate this approach.

INTRODUCTION

Recent earthquake events in South Lebanon in addition to the accidental topography and heterogeneous Lebanese geology have trigerred this research on the analysis of slopes subjected to seismic loads using a probabilistic approach.

The stability analysis of earth slopes located in seismic areas has been extensively investigated in literature using deterministic approaches. In these methods, average values of the input parameters (angle of internal friction, cohesion, seismic coefficient, etc.) are used and a global safety factor is calculated (Refer for instance [START_REF] Chen | Earthquake-induced slope failure in nonhomogeneous, anisotropic soils[END_REF][START_REF] Leshchinsky | Pseudostatic seismic stability of slopes: Design charts[END_REF][START_REF] You | Displacement charts for slopes subjected to seismic loads[END_REF][START_REF] Michalowski | Stability charts for uniform slopes[END_REF]Loukidis et al., 2003 among others). A reliability-based approach for the slope stability analysis in seismic zones is more rational than the traditional deterministic methods since it takes into account the inherent uncertainty of each input variable.

The reliability theory was introduced by several authors in the analysis of slope stability but without taking into consideration the seismic loading (Refer for instance [START_REF] Vanmarcke | Reliability of earth slopes[END_REF][START_REF] Chowdhury | Rational polynomial technique in slope-reliability analysis[END_REF][START_REF] Christian | Reliability applied to slope stability analysis[END_REF][START_REF] Low | Slope reliability analysis using generalized method of slices[END_REF][START_REF] Malkawi | Uncertainty and reliability analysis applied to slope stability[END_REF]Auvinet & Gonzales, 2000;[START_REF] El-Ramly | Probabilistic slope stability analysis for practice[END_REF]and Sivakumar Babu & Mukesh, 2004 among others). Most of these studies are based on approximate deterministic models as the limit equilibrium method which is founded on a priori assumptions concerning (i) the form of the slip surface and (ii) the normal stress distribution along this surface or the inter-slice forces. The application of the reliability theory in slope stability problems subjected to seismic loadings is much less investigated (Refer for instance [START_REF] Christian | Probabilistic evaluation of earthquake-induced slope failure[END_REF][START_REF] Al-Hamoud | Reliability analysis of three-dimensional dynamic slope stability and earthquakeinduced permanent displacement[END_REF][START_REF] Michalowski | Stability charts for uniform slopes[END_REF][START_REF] Abdel Massih | Reliability analysis applied on lebanese slopes subjected to seismic loads[END_REF].

In this paper, a reliability-based analysis of an earth slope situated in a seismic area is performed. The deterministic model used in the reliability analysis is based on a rigorous coupled analysis that simultaneously captures the fully non linear response of the soil and the acceleration time history record. This model is based on numerical simulations using the dynamic option of the finite difference FLAC 3D software. The time history of acceleration records used in the analysis is the one digitized at the Lebanese Geophysical Center of the National Council for Scientific Research. The random variables considered in the analysis are the cohesion c and the shear modulus G of the soil since it was shown that these parameters have the most influence on the slope displacements (Dhakal, 2004). The performance function used in the reliability analysis is defined with respect to the horizontal permanent displacement at the toe of the slope. The Stochastic Response Surface Methodology (SRSM) is utilized for the assessment of the Γfailure probability is evaluated with respect to the exceedance of an allowable maximun permanent displacement. After a brief description of the concepts underlying the stochastic response surface methodology, the seismic deterministic numerical modeling of a slope using FLAC 3D is described. Then, the probabilistic numerical results are presented and discussed. A case study of a typical Lebanese slope illustrates this approach at the end of the paper.

STOCHASTIC RESPONSE SURFACE METHODOLOGY

The stochastic response surface method (SRSM) is a technique for the reliability analysis of complex systems with low failure probabilities, for which Monte Carlo simulation (MCS) is too computationally expensive and for which other approximate methods (such as FORM) may be inaccurate. It allows one not only to determine the failure probability of a system due to a prescribed load as is the case of the traditional response surface methodology (RSM) (e.g. Youssef Abdel [START_REF] Abdel Massih | Reliabilitybased analysis of strip footings using response surface methodology[END_REF], but also to find the full probability distribution of the system response. Typically, SRSM approximates an analytically-unknown limit state function with a multidimensional chaos polynomial by fitting the polynomial to a number of sampling points. Once the coefficients of the polynomial chaos are determined, the probability distribution of the response can be calculated using the classical and wellknown Monte Carlo simulation based on the obtained analytical form of the response. Then, the statistical moments of the response (mean, variance, skewness, kurtosis, etc.) may be computed. Furthermore, the analytical equation of the response may be used to calculate the reliability index and the failure probability using the classical methods such as FORM, Monte Carlo, Importance Sampling, etc. The steps used in SRSM can be summarized as follows:

Step I: Representation of Stochastic Inputs The first step in the application of the SRSM is the representation of all the model inputs in terms of a set of "standardized random variables". Normal random variables are selected from a set of independent standard distributed normal random variables,   1 ,........., M   , where M is the number of independent inputs, and each i  has zero mean and unit variance. One of two approaches may be taken to represent the uncertain inputs in terms of the selected "standardized random variables": (a) direct transformation of inputs in terms of the "standardized random variables", and (b) series approximations in "standardized random variables". [START_REF] Devroye | Non-Uniform Random Variate Generation[END_REF] presents transformation techniques and approximations for a wide variety of probability distributions of random variables (cf. Table 1) Step II: Analytical approximation of outputs An output of a model (e.g. the permanent displacement at the toe of the slope or the amplification of the slope in the case of the present seismic study) may be influenced by any number of model inputs. Hence, any general functional representation of uncertainty in model outputs, should take into account uncertainties in all inputs. For a deterministic model with random inputs, if the inputs are represented in terms of the set

  1 M i i
  (as descrribed in step I), the output metrics can also be represented in terms of the same set, as the uncertainty in the outputs is obviously due to the uncertainty of the inputs. The present SRSM method addresses one specific form of representation: the series expansion of normal random variables in terms of Hermite polynomials. The expansion is called "polynomial chaos expansion". Using this series expansion, an output can be approximated by a polynomial chaos expansion on the set   1 M i i   , given by:

    1 1 0 ( ) ( ) P M j j k k j y a          (1)
where y is any output metric (or random output) of the model. The a j represent deterministic constants to be estimated, j  are multi-dimensional Hermite polynomials of order less than or equal to p. They are given by :

            T T e e ip i p p ip i p 2 1 1 2 1 1 ... 1 ,...,       
Γ p are random variables since they are functions of the random variables   1 n i i   . In equation (1), P denotes the size of the polynomial chaos basis which is given by:

 ! ! ! M p P M p   (2)
For a fourth order polynomial and for two random variables, the size of the polynomial chaos basis (which is equal to the number of unknown coefficients) is found to be equal to 15.

Step III: Choice of the chaos order and collocation points In general, the accuracy of the approximation increases as the order of the polynomial chaos expansion increases. The order of the expansion can be selected to reflect accuracy needs and computational constraints. After choosing the order of the chaos and the number of the random variables, one can obtain the number of available collocation points and their positions in the standard space of random variables. These points will be used for the evaluation of the system response. In the case of two random variables (as in this paper), the choice of these points is made so that each standard variable (ξ 1 or ξ 2 in this case) can take the values of the roots of the univariate Hermite polynomial of order n+1. Hence, the available collocation points are obtained by assigning to each root of the first variable the different roots of the second variable.

Table 2 provides the roots and the general shape of the Polynomial Chaos Expansion that will approximate the response of the model (i.e. the toe horizontal permanent displacement or the amplification), for the 4 th order of the polynomial chaos (i.e., p = 4). Figure 1 provides the positions of the collocation points (25 points) in the standard space of the random variables.

Table 2 :

Roots to be used for the 4 th order poynomial chaos and the expression of the corresponding output.

Type of distribution

Transformation*

Uniform (a,b)     1 1 / 2 2 2 a b a erf           Normal (μ,σ) μ+σξ Lognormal (μ,σ) exp(μ+σξ) Gamma (a,b) 3 a 9 1 1 a 9 1 ab            Exponentiel (λ)            2 / erf 2 1 2 1 log 1 Weibull (a) a 1 y Extreme Value -log(y)
Roots to be used for the polynomial chaos of order 4

Expression of the output

  0 ; 5 10                         4 0,0 0,0 1,0 1,0 1 0,1 0,1 2 2,0 2,0 1 1,1 1,1 1 2 0,2 0,2 2 3,0 3,0 1 2 , 1 2 , 1 1 2 1 , 2 1 , 2 1 2 0 , 3 0 , 3 2 4,0 4,0 1 3,1 3,1 1 2 2,2 ... ... , ... ... , , ... ... , U a a a a a a a a a a a a a                                                  2,2 1 2 1,3 1,3 1 2 0,4 0,4 2 , ... ... , a a           
1 , 1, 2 , 2, C m m C m m H                   (3)
where:

  H Chol   (4) and , , 1 1 c G c G           (5) Notice that   m C m C , 2 , 1 ,  in equation (5) is the standard correlated collocation point corresponding to   1, 2, , m m  
. Once the correlation between random variables has been performed, the standard correlated variables have to be transformed into the physical correlated variables in conformity with the marginal distribution of each variable. For the well-known statistical distributions (such as lognormal or beta), a direct transformation is possible using a transformation function which depends on the distribution. The transformation from standard to physical variables can also be done in Matlab using the following expression:

    1 X F     (6)
In this equation, X is the physical random variable, ξ is the standard normal variable and

    and    F
are the corresponding CDFs. In the present case:

        1 1 , 1 2 , m G Cm m c Cm G F c F              (7)
For the collocation points in the standard space, it is then possible to find the corresponding points in the physical space, and to call the deterministic model for each of these points and determine the M corresponding system responses.

Step IV: Estimation of parameters in the analytical functional approximation of the output Once the system responses are determined, several methods exist for the estimation of the unknown coefficients (a ik … 's) (cf. equation 1) in the polynomial chaos expansion. The one used in this paper is the Regression-Based Method. This method is presented in Sudret ( 2007). It is based on the traditional least square method which minimizes the residual between the exact solution and the approximated one. For the choice of the collocation points among those available, we have used here the total number of the available collocation points.

Step V: Estimation of the statistics of the output Once the coefficients (a ik … 's) used in the series expansion of model outputs are estimated, the statistical properties of the outputs, such as the density functions, moments, joint densities, joint moments, correlation between two outputs, or between an output and an input, etc…, can be calculated. To accomplish this task, a generation of a large number of realizations of the "standard normal variables" is performed. This results in a large number of samples of inputs and outputs. The samples can then be statistically analyzed using standard methods. It should be noticed that the calculation of model inputs and outputs involves evaluation of simple analytical expressions, and does not involve numerical model simulations runs which means that a significantly lower computational time is required.

As an example, if the inputs i

x 's are represented as

  i i i x F   ,
and if the output i y 's are estimated as

  1 2 , ,... i i n y G     , then
the following steps are involved in the estimation of the statistics of the inputs and outputs.

 generation of a large number of samples of  

1 2 , , ,... n i    ,
 calculation of the values of input and output random variables from the samples,  calculation of the moments using equations 8, 9, 10 and 11 given below,  calculation of density functions and joint density functions using the sample values.

From a set of N samples, the moments of the distribution of an output i y can be calculated as follows:

1

1 n i i y n     (8)   2 2 1 1 1 n i i y n        (9)      3 1 3 1 1 2 n i i n g y n n         (10)              2 4 2 4 1 1 3 1 1 2 3 2 3 n i i n n n g y n n n n n               ( 11
)
where g 1 and g 2 are respectively the Skewness and Kurtosis of the output random variables y i . Similarly, higher moments of outputs, or the correlation between two outputs, or between an input and an output can be directly calculated.

Step V: Evaluation of Convergence of Approximation

Once the coefficients of the polynomial chaos expansion are obtained, the convergence of the approximation is determined through comparison with the results from a higher order approximation. The next order polynomial chaos expansion is used, and the process for the estimation of unknown coefficients is repeated. If the estimates of probability density functions (pdfs) of output metrics agree closely, the expansion is assumed to have converged, and the higher order approximation is used to calculate the pdfs of output metrics. If the estimates differ significantly, yet another series approximation, of the next order, is used, and the entire process is repeated until convergence is reached. In the present work, it was shown that a chaos polynomial of the fourth order is sufficient to achieve a good convergence of the method.

PROBABILISTIC NUMERICAL MODELING OF A SLOPE

FLAC 3D code FLAC 3D (Fast Lagrangian Analysis of Continua) is a commercially available three-dimensional finite difference code in which a Lagrangian explicit calculation scheme and a mixed discretization zoning technique are used. In this code, although a static (i.e. non-dynamic) mechanical analysis is required, the equations of motion are used. The solution to a static problem is obtained through the damping of a dynamic process by including damping terms that gradually remove the kinetic energy from the system. The application of velocities, stresses or forces on a system creates unbalanced forces in this system. Damping is introduced to remove these forces or to reduce them to very small values compared to the initial ones. The stresses and deformations are calculated at several small time-steps or cycles until a steady state of static equilibrium or plastic flow is achieved. The convergence to this state may be controlled by a prescribed value of the unbalanced force.

Slope Model

This section focuses on the determination of the probability distribution of the horizontal permanent displacement at the toe of a slope in the presence of earthquake loading. The Stochastic Response Surface Methodology based on FLAC 3D numerical simulations is used.

A slope with height equal to 7 m and an angle equal to 70 o is considered in the analysis. The soil domain and mesh used are shown in fig. 2. Since this is a 2D case, all displacements in the Z direction are fixed. For the displacement boundary conditions in the (X,Y) plane, the bottom horizontal boundary was subjected to an earthquake acceleration signal and absorbent boundaries were applied to the right and left vertical boundaries.

A conventional elastic-perfectly plastic model based on the Mohr-Coulomb failure criterion is used to represent the soil. The cohesion c and the shear modulus G of the soil are random variables since it was shown that these parameters have the highest influence on the slope displacements (not shown here).

For the statistical parameters of the random variables, the illustrative values used for the mean of the cohesion and the shear modulus are given as follows: 15 kPa . A damping of 5% is used for the soil through the Rayleigh damping which is proportional to the stiffness and mass of the soil system.

(

) 12
where C is the damping matrix, M is the mass matrix and K is the stiffness matrix. The Rayleigh coefficients α and β can be estimated from the modal damping coefficient as follows:

(13)

Where w is the natural frequency of the soil.

Fig. 2. Soil domain and mesh used in FLAC 3D

A time-history earthquake acceleration record registered at the Lebanese station "BHL" at Bhannes city was used (cf. Figure 3). This record is registered 60 km far from the epicenter of a 5.1 magnitude earthquake that occurred on February 15, 2008. The epicenter of this earthquake was located at: 33.35 N ; 35.36 E (i.e. 16 km NE of the city of Sour). Only, the 100 time scaled North-South corrected horizontal component of the acceleration record is used. The scaled maximal acceleration obtained from this record is found equal to 2.68 m/s 2 .

Concerning the cost of the calculations, the computational time for one deterministic simulation was found to be approximatly equal 45 minutes using an Intel 2.40 GHz PC. In order to evaluate the coefficients of a 4 th order polynomial chaos expansion of the slope response using SRSM, 25 deterministic calculations (as shown previously in step III of the SRSM) were necessary for given soil variability and slope geometry.

RESULTS

Deterministic permanent horizontal displacement at the slope toe

For the slope geometry and the mean values of the soil characteristics described in the previous section, figure 4 presents the toe relative displacement as a function of the dynamic time. This displacement is computed as the difference between the horizontal toe displacement and the horizontal base displacement during the earthquake shaking. The maximum value of this displacement is called the toe permanent horizontal displacement. It is found to be equal to 21 cm as shown in fig. 4.

Fig. 4. Relative displacement at the toe of the slope versus the dynamic time

The contour of shear strain increment that gives an idea about the probable failure surface is presented in figure 5.

Fig. 5. Contour of shear strain increment

System response, limit state surface and probability distribution of the toe permanent displacement - For different values of an allowable maximal permanent displacement u max , one can plot the limit state surfaces defined by G=u-u max =0 that divides the combinations of   , c G that would lead to failure (i.e. the exceedance of u max ) from the combinations that would not. These limit state surfaces are plotted on figure 7 in the standard space of the random variables for different values of u max varying from 30 cm to 90 cm. One can notice that the limit state surfaces are not linear. Consequently, the classical FORM method that linearly approximates the response surface for the evaluation of the failure probability may lead in this case to inaccurate results. Thus, the determination of the entire probability distribution of the response using the SRSM method is necessary.

It was shown earlier in the present paper that the SRSM method allows one to determine the full probability distribution of the system response. The cumulative distribution function of the permanent horizontal displacement is presented in figure 8 From these functions, the failure probability with respect to the exceedance of an allowable maximum permanent displacement can be easily evaluated. If we consider an allowable maximal permanent displacement u max of 70 cm, the failure probability can be evaluated as follows:

    max max 1 1 0.9633 0.0367 (3.67%) f P u u u u           .
Effect of the coefficients of variation of the random variables on the mean value of the slope permanent displacement  is plotted versus the coefficient of variation of G and c . For each curve, the coefficient of variation of a parameter is held to a constant value (20% for the cohesion and 30% for the shear modulus) and the coefficient of variation of the second parameter is varied. The results have shown that u  is higher than the deterministic permanent displacement calculated using the mean values of the random variables (i.e. ). It can also be noted that the mean value of the permanent displacement is highly influenced by the coefficient of variation of the cohesion, the greater the scatter in c the higher is the permanent displacement. This means that the accurate determination of the distribution of this parameter is very important in obtaining reliable probabilistic results. . In order to understand and interpret this observation, the mean value of the amplification A  at the crest of the slope is plotted versus the coefficient of variation of G in figure 11. The amplification in this study refers to the ratio of the maximum output horizontal acceleration at the crest of the slope to the input maximal acceleration at the base of the slope. The amplification is considered here as an output random variable determined by the SRSM. It was noticed that the maximal value of A  corresponds to the maximal value of the mean permanent horizontal displacement (i.e. for 50% Cov G 

). One can conclude that for these values of the soil variability, the amplification of the wave and consequently the mean permanent horizontal displacement of the toe are maximal. The geologic formation of the slope is sandstone, composed by cretaceous sands, mixed with irons, and some intercalation of clay (cf. Figure 12). Since the developed probabilistic models are 2D, the study is done in a 2D cross section. The most critical cross section (AA') at the middle of the slope (i.e. at the greatest height) was considered. The height of this section is equal to 10 m and the average slope angle is found to be approximately equal to 35 o . . Since it was shown earlier in this paper that the cohesion and the shear modulus of the soil have the greatest influence on the slope permanent displacement when performing a dynamic analysis, the uncertainties on the friction angle is neglected. It is considered as a deterministic parameter with value equal to its mean value    . The range of the Young modulus found in the literature for a sandstone formation mixed with clay is between 25 and 200MPa. The value of Poisson ratio ranges between 0.25 and 0.33. The corresponding range of the shear modulus varies between 10 and 80 MPa. Consequently, a mean value of 50MPa is used and a high coefficient of variation of 50% is considered.

Probability distribution of the permanent horizontal displacement of Baabdat slope Figure 13 presents the PDF of the permanent displacement of the Baabdat Slope for the same scaled earthquake record presented in figure 3 (a max =2.64cm/s 2 ). By considering an allowable maximal permanent displacement u max of 70 cm, the failure probability is found equal to 85%. This slope is a typical one of the studied region. It has to be noticed also that the mean value of the permanent displacement is found equal to 1.29 m higher than the deterministic permanent displacement equal to 1.16 m calculated for the mean values of the soil characteristics.

CONCLUSION

This paper examines a dynamic analysis of slope stability using the reliability approach. A numerical simulation is carried through finite difference software. For the acceleration timehistory, a record of a local seismic event scaled to 2.68 m/s 2 has been used and applied on a local case study slope stability problem. The performance function used in the relability analysis is defined with respect to the permanent displacement at the toe of the slope. It was found that the considered variables cohesion c and shear modulus G of the soil have the most significant effect on displacement.

As for the estimation of the probability distribution of the permanent displacement at the toe, a Stochastic Response Surface Methodology has been adopted. From this approach, it was shown that the mean value of the permanent displacement is very sensitive to the coefficient of variation of the cohesion, the greater the scatter of c the higher is the permanent displacement. Besides, the failure probability is evaluated with respect to the exceedence of an allowable maximal permanent displacement. It was found that the maximal value of the amplification mean A  corresponds to the maximal value of the mean permanent displacement (i.e. 50% Cov G 

). One can conclude that for adopted values of the soil variability, the amplification of the wave and consequently the mean permanent displacement at the toe are maximal. Finally, the analysis is applied on a case study of a slope in the Baabdat area of the Metn Mountains North of Beirut. This case study represents a typical slope of the studied region. The results showed that when the permanent displacement u max is 70 cm, the failure probability is found equal to 85%. The probabilistic mean value of the permanent displacement found equal to 1.29 m is higher than the deterministic permanent displacement equal to 1.16 m calculated for the mean values of the soil characteristics.

Fig. 1 :

 1 Fig. 1: Position of the collocation points in the standard space of the random variables Translation of the collocation points from the standard to the physical space For each collocation point   1, 2, , m m   , one has to determine the corresponding physical point   , m m c G (i.e the cohesion and the shear modulus) to be used in the deterministic model. If there is a correlation , c G  between the random variables m c and m G , one should correlate the standard variables by multiplying the vector of standard uncorrelated variables by Cholesky transform H of the standard covariance matrix Σ (obtained from the covariance matrix C) as follows:

  probability distributions, the two random variables are considered as lognormally distributed. The soil unit weight  , the angle of internal friction  , the dilation angle  and the bulk modulus K are considered as deterministic parameters. Their values are as follows:

Figure 6

 6 Figure 6 presents the shape of the response surface of the permanent horizontal displacement approximated by the

Fig. 7 .

 7 Fig. 7. Limit state surfaces in the standard space of random variables for different values of u max

Fig. 9 .

 9 Fig. 8. Cumulative probability function of the permanent horizontal displacement at the toe slope

Figure 10

 10 Figure10presents the variation of the mean value of the permanent horizontal displacement u  at the toe of the slope as a function of the coefficient of variation of the random variables. In this figure, u

Fig. 10 .

 10 Fig. 10. Mean of the slope permanent displacement as a function of coefficient of variation of c and G From the same figure, one can notice that u  presents a maximum value for 50% Cov G . In order to understand and interpret this observation, the mean value of the amplification

Fig. 12 .

 12 Fig. 12. View of a slope cut Several soil samples were taken in order to perform "Laboratory Direct Shear tests" for the quantification of the variability of the soil shear strength parameters c and . The number of tests carried out at this slope area is 10. The statistical parameters of the soil shear strength were found equal to : 10  c  kPa,

Fig

  Fig. 13. PDF of the permanent displacement when a max =2.64cm/s 2

Table 1 :

 1 Representation

	of common univariate distributions as
	function of normal random variables
	* y is the exponential (1) distribution