key: cord-0461159-ieooezme
authors: Deason, Alis J.; Oman, Kyle A.; Fattahi, Azadeh; Schaller, Matthieu; Jauzac, Mathilde; Zhang, Yuanyuan; Montes, Mireia; Bah'e, Yannick M.; Vecchia, Claudio Dalla; Kay, Scott T.; Evans, Tilly A.
title: Stellar splashback: the edge of the intracluster light
date: 2020-10-06
journal: nan
DOI: nan
sha: 4bf04033c02048d9ac6258e8dda640a0b62b11ab
doc_id: 461159
cord_uid: ieooezme

We examine the outskirts of galaxy clusters in the C-EAGLE simulations to quantify the"edges"of the stellar and dark matter distribution. The radius of the steepest slope in the dark matter, otherwise known as the splashback radius, is located at ~r_200m; the strength and location of this feature depends on the recent mass accretion rate, in good agreement with previous work. Interestingly, the stellar distribution (or intracluster light, ICL) also has a well-defined edge, which is directly related to the splashback radius of the halo. Thus, detecting the edge of the ICL can provide an independent measure of the physical boundary of the halo, and the recent mass accretion rate. We show that these caustics can also be seen in the projected density profiles, but care must be taken to account for the influence of substructures and other non-diffuse material, which can bias and/or weaken the signal of the steepest slope. This is particularly important for the stellar material, which has a higher fraction bound in subhaloes than the dark matter. Finally, we show that the"stellar splashback"feature is located beyond current observational constraints on the ICL, but these large projected distances (>>1 Mpc) and low surface brightnesses (mu>>32 mag/arcsec^2) can be reached with upcoming observational facilities such as the Vera C. Rubin Observatory, the Nancy Grace Roman Space Telescope, and Euclid.

The dark matter haloes that underpin our hierarchical structure formation paradigm do not have uniquely defined boundaries. Several common definitions are used in the literature: the "friends-offriends" distance (Davis et al. 1985) , the virial radius (Bryan & Norman 1998) , and a radius within which the mean density equals a fixed value times the critical density or the cosmic mean value (spherical overdensity halo boundaries 1 ). Regardless of the exact definition, the use of a halo boundary is essential in order to define halo masses, distinguish between field and satellite galaxies, and, importantly, contrast the predictions of simulations with observations. Often the choice of halo boundary depends on the mass-scale under consideration, and whether the study is observationally or theoretically motivated. A common definition across halo mass and redshift that is also observationally motivated (or even applicable) is crucial.

Recent work has argued that the most physical definition of the halo boundary is related to the transition between collapsed and infalling material, or the one-halo and two-halo regimes, and has been termed the "splashback" radius (e.g Adhikari et al. 2014; Diemer & Kravtsov 2014; More et al. 2015; Diemer 2020) . Diemer & Kravtsov (2014) showed that the splashback radius in cosmological simulations can be identified from the radius of steepest slope in the density profiles of the dark matter. This splashback radius does not only define a physical halo boundary, it also crucially depends on the mass accretion rate of the collapsing halo (e.g Adhikari et al. 2014; Diemer & Kravtsov 2014; Diemer et al. 2017) . Importantly, there is now considerable evidence that the splashback radius has been identified in galaxy clusters, either through stacked satellite galaxy surface density profiles (More et al. 2016; Baxter et al. 2017; Murata et al. 2020) , weak lensing (Chang et al. 2018; Contigiani et al. 2019; Tam et al. 2020) or the Sunyaev-Zel'dovich effect (Shin et al. 2019; Zürcher & More 2019) . The location of the observed splashback radius tends to be lower (by ∼20%) than the predictions of ΛCDM simulations, however, selection effects from the optical cluster finding algorithms (e.g. Busch & White 2017; Xhakaj et al. 2019; Murata et al. 2020) , and the effect of dynamical friction on massive satellites (Adhikari et al. 2016 ) could explain the apparent discrepancy. However, while the effects mentioned above appear more likely, these differences with the ΛCDM predictions could also signal more exotic solutions, such as self-interacting dark matter (More et al. 2016; Banerjee et al. 2020) .

The theoretical background to spherical overdensity halo boundaries, and the more recently promoted splashback radius, is based almost entirely on the dark matter distribution. This is perhaps unsurprising, as the outer reaches of galaxy and cluster haloes are dominated by dark matter, and the majority of the visible material is concentrated at the very centre. However, the use of an observationally motivated halo boundary, defined using the baryonic material, is, in some cases, more attractive. Satellite galaxies are an obvious way forward: they can be identified with photometric and spectroscopic surveys, and can reach out to large distances in galaxy haloes (e.g van den Yang et al. 2005; Robotham et al. 2011; Budzynski et al. 2012; McConnachie 2012) . The main drawback is that the number of visible satellite galaxies can be low for individual systems, and care must be taken to understand the selection effect of a stellar mass limited sample (see e.g. Adhikari et al. 2016 ) and the impact of satellite galaxy colour (e.g. Baxter et al. 2017; Shin et al. 2019; Adhikari et al. 2020 ). An obvious, complementary, probe of the halo is the remains of destroyed satellite galaxies (i.e. the stellar halo), which are commonly referred to as the intracluster light (ICL) at cluster scales (e.g. Mihos 2016; Montes 2019).

Recently, Deason et al. (2020) provided the first foray into defining the stellar edges of galactic haloes (with halo mass ∼ 10 12 M ). These authors used high resolution cosmological hydrodynamic simulations of Milky Way-mass haloes to explore the boundary of the halo stars. Curiously, they found that the stars have a well-defined edge, but this is not coincident with the splashback radius of the dark matter. Rather, the edge of the halo stars appears to be related to a secondary dark matter caustic (termed the "second caustic" in this work), which is likely related to the apocenter of the last significant major merger. However, extrapolating these findings to other mass-scales is non-trivial owing to the nonlinear stellar mass-halo mass relation (Moster et al. 2010; Behroozi et al. 2013 ) and the varying importance of "smooth" accretion onto galaxy haloes with halo mass (e.g Genel et al. 2010; Fakhouri & Ma 2010) .

The stellar haloes of Milky Way-mass galaxies are primarily built from the left-over debris from destroyed dwarf galaxies (e.g Bullock & Johnston 2005; Cooper et al. 2010; Deason et al. 2015 Deason et al. , 2016 . On cluster-mass scales, the ICL is built predominantly from the destroyed remnants of Milky Way-mass galaxies (e.g. Murante et al. 2004; Conroy et al. 2007; Purcell et al. 2007; Puchwein et al. 2010; Contini et al. 2014; Montes & Trujillo 2014 De-Maio et al. 2018 ). This self-similarity from dwarf galaxies to Milky Way-mass galaxies to clusters is a beautiful example of hierarchical structure formation in action. Nonetheless, the detailed properties of galactic stellar haloes and the ICL have important differences, most notably the significance of this component to the total stellar mass, and their radial distributions (e.g. Pillepich et al. 2014 Pillepich et al. , 2018 . Indeed, while both galactic stellar haloes and the ICL form via merg-ers, important galaxy formation physics underpins their differences and similarities. Thus, studying these diffuse halo components over a range of mass scales allows for a critical view on both structure formation and models of galaxy formation.

In this work, as a complement to the Deason et al. (2020) study, we focus on the stellar haloes of clusters, i.e. the ICL. Recent work has shown an intriguing similarity between the dark matter density profiles of clusters and their ICL (e.g. Pillepich et al. 2018; Montes & Trujillo 2019; Alonso Asensio et al. 2020 ). This finding particularly motivates an investigation into the stellar edges of cluster-mass haloes. Here, we use the Cluster-EAGLE (C-EAGLE) suite of simulations to study the outer density profiles of both stars and dark matter, and their relation to each other. The paper is arranged as follows. In Section 2, we describe the C-EAGLE simulations, and in Section 3, we probe the edges of these galaxy clusters using both stars and dark matter. We explore the observationally motivated projected density profiles in Section 4, and discuss the implications for current and future observational probes of the ICL. Finally, we summarise our main conclusions in Section 5.

In this work we use the Cluster-EAGLE (C-EAGLE) project (Bahé et al. 2017; Barnes et al. 2017) to study the outer density profiles of galaxy clusters. This suite is a set of = 30 zoom-in cosmological hydrodynamical simulations of galaxy clusters in the mass range 10 14.0 < 200c /M < 10 15.4 . The simulations are run with the EAGLE galaxy formation model (AGNdT9 calibration, Schaye et al. 2015) , with a gas particle mass of 1.8 × 10 6 M , a dark matter particle mass of 9.7 × 10 6 M , and a physical softening at < 2.8 of 0.7 kpc. The high resolution volumes are set up such that they are devoid of any low resolution particles within at least 5 200c , and the clusters were selected to have no massive neighbours within 10 200c . Here, 200c is the radius at which the average density drops to 200 times the critical density at = 0. A subset (24) of the C-EAGLE sample has been simulated at high resolution out to at least 10 200c ; these are called the Hydrangea simulations (Bahé et al. 2017) . The simulations assume a flat ΛCDM cosmology with parameters (Planck Collaboration et al. 2014) : Ω m = 0.307, Ω b = 0.04825, Ω Λ = 0.693, ℎ = 0.6777, 8 = 0.8288 and = 0.9611. The EAGLE model is described in detail in Schaye et al. (2015) and Crain et al. (2015) , and includes subgrid models for baryonic processes such as star formation, stellar winds, gas cooling, metal production and stellar and black hole feedback. These subgrid recipes were calibrated to reproduce the present day stellar mass function, the galaxy size-stellar mass relation and the black hole mass-host galaxy mass relation. Note that since the EAGLE model was calibrated on galaxy properties, and not specifically on clusters, the properties of the C-EAGLE cluster sample are predictions of a model that produces realistic galaxies in the field. Projected images for two example clusters are shown in Fig. 1 . Here, we show the dark matter (left panels) and stellar mass (right panels) distributions.

The low redshift global properties of the C-EAGLE sample are described in Barnes et al. (2017 , see also Bahé et al. 2017 . These works showed that the total stellar content, metal content (see also Pearce et al. 2020 ) and black hole masses are in good agreement with the observations. However, the clusters are too gas rich, their central temperatures are too high, and they have larger entropy cores than observed. These mismatches with observations are likely driven by shortcomings in the AGN feedback model. . The color scale is logarithmic, with projected density values ranging from 1 × 10 −1 − 1 × 10 4 M /pc 2 and 1 × 10 −3 − 3 × 10 2 M /pc 2 for the dark matter and stellar distributions, respectively. This image was produced using the open source project yt (Turk et al. 2011). to this work, Alonso Asensio et al. (2020) recently studied the ICL of the C-EAGLE clusters and found that the shape of the stellar mass distribution closely follows that of the total matter, in good agreement with observations (Montes & Trujillo 2019). Moreover, Bahé et al. (2020, in preparation) also find that the ICL surface density profiles agree with observations. In this work, we focus on the "edges" of these clusters and compare the stellar and dark matter halo boundaries.

In this section, we probe the dark matter and stellar density profiles of the C-EAGLE clusters. Following the work by Diemer & Kravtsov (2014) , we use the differential logarithmic density profiles to identify the steepest slope, which signifies a transition between the collapsed (one-halo) and infalling (two-halo) material. We consider all 2 dark matter and stellar particles in the simulations out to a 3D radius of 4 200m from the halo centre. Here, 200m is the radius at which the average density drops to 200 times the universal matter density at = 0. Throughout this work, we scale physical radii with this radius. Our outer boundary of 4 200m sometimes contains a small fraction of low resolution dark matter particles. However, this makes little difference to our results as we are mainly interested in the region within ∼2 200m , which is completely devoid of any low The logarithmic profile is computed using the fourth-order Savitzky-Golay smoothing algorithm over the 15 nearest bins (Savitzky & Golay 1964) . The clusters are ordered according to the recent mass accretion rate, Γ, increasing from the top left panel to the bottom right panel. The cluster ID is also indicated (see Barnes et al. 2017 Table A1 ). The solid vertical lines show the most prominent minimum, defined as Caustic , or the splashback radius. We also show with the dotted lines cases with clear second caustics. These are much weaker than the splashback radii, and are more common amongst haloes with low recent mass accretion rates. resolution particles. Note, for ease of comparison, 200m ∼ 1.7 200c at the cluster mass scale.

For both dark matter and stars we construct density profiles in 40 evenly spaced logarithmic bins between 0.1 and 4 / 200m . We follow a similar approach to , see their Section 4.3) in order to construct the angular median density pro-file. Namely, for each logarithmic radial shell, the density profile is computed in = 50 (equally spaced) solid angle segments. We construct the density profile by taking the median of these profiles in each radial shell. This procedure minimizes the influence of massive substructures and other non-diffuse structures on the density profile. As we will show in Section 4, the median angular profile is far more effective at isolating the steepest halo slope than the more commonly used mean; this is particularly important for the stellar distribution. Finally, we compute the logarithmic slope profiles using a fourth-order Savitzky-Golay smoothing algorithm (Savitzky & Golay 1964 ) over the 15 nearest bins. This bin size and smoothing was chosen to minimize noise, while allowing us to identify the strongest features in the profile.

The logarithmic slope profiles of the dark matter density profiles for the = 30 C-EAGLE clusters are shown in Fig. 2 . The clusters are ordered according to the recent mass accretion rate, Γ, increasing from the top left panel to the bottom right panel. Here, we define mass accretion rate as:

where 1 = 1( 1 = 0) and 2 = 0.667( 2 = 0.5) (Diemer & Kravtsov 2014) . Note that the mass accretion rate is calculated Γ > 2.0 (N=12) 1.0 < Γ < 2.0 (N=10) Figure 4 . The stacked logarithmic slope profiles of the dark matter (solid black) and stellar (dashed red) density profiles of the C-EAGLE clusters. Each individual profile is determined using the angular median method and we show the median over all clusters. In the left panel the median of all = 30 clusters is shown. In the remaining panels, we show the median profiles for relatively low (Γ < 1.0), medium (1.0 < Γ < 2.0) and high (Γ > 2.0) recent mass accretion rates, respectively. The stellar distribution shows a caustic feature coincident with the dark matter. This caustic is stronger and located at smaller radii for higher mass accretion rates.

using the virial mass, vir , which is defined using the Bryan & Norman (1998) formalism (with density contrast of ∼102 relative to the critical density at = 0). For the majority of systems, a prominent "dip" is seen in the slope profiles, which previous works have labeled as the splashback radius (e.g. Diemer & Kravtsov 2014; More et al. 2015) . In this work, we also define this radius of steepest slope as the splashback, and indicate these with the solid blue lines in Fig. 2 . In agreement with previous work, this feature tends to become more pronounced at higher mass accretion rates. We note that CE-05 is currently undergoing a major merger, and hence the splashback feature, particularly in the dark matter distribution, is washed out.

In some cases, we also identify a secondary caustic feature in the density profiles (shown with the dotted blue lines). These are located at smaller radii, and have shallower slopes than the splashback. Deason et al. (2020) labeled these features as "second caustics", and we adopt this terminology here. Note, however, that these features don't necessarily relate to the classical definition of second caustic from spherical (or ellipsoidal) collapse models (see e.g. Adhikari et al. 2014) , and could have multiple origins. As these features are much weaker than the splashback, we must caution against fitting to noise. To this end, we only consider second caustics that have slopes steeper than −2.5 and the difference between the local minimum and maximum is greater than 0.5 dex. Interestingly, the second caustics are more common amongst haloes with low mass accretion rates, which is what the Adhikari et al. (2014) models predict. These features certainly deserve further scrutiny, and this will be a topic of future work. In this work, we focus on the splashback radii, and now turn our attention to the stellar distribution.

In Fig. 3 we show the logarithmic slope profiles for the stellar material. Here, we consider all stars in the cluster, and do not try to distinguish between the brightest cluster galaxy, diffuse stellar material, or the stars bound in subhaloes. Any biases caused by massive substructures are mitigated by the angular median method used to calculate the density profiles. Note, however, that these profiles are not "pure" ICL, as we have not explicitly removed stars bound to subhaloes. There are a variety of different definitions of the ICL in the literature, which can lead to significant differences in the derived ICL properties (see e.g. Rudick et al. 2011; Montes & Trujillo 2019) . Here, we consider all distant halo stars, and use the angular median method to minimise the effects of massive substructures and other non-diffuse structures. In practice, this approach is appealing as it can, potentially, be used in both simulations and observations. Conversely, simply removing all stars that are bound to subhaloes is an approach that is not directly applicable to observations. Furthermore, the identification of bound subhaloes depends on the algorithm used, and the resolution of the simulation (see Section 4 for further discussion).

The stellar profiles in Fig. 3 look similar to the dark matter profiles: a prominent dip is seen in almost all cases, and in some cases a second caustic-type feature is also apparent. There are, however, some differences. Most notably the scales in Fig. 3 are different. While the outer caustics in the dark matter tend to have slopes of ∼ −4.5, the stellar caustics are much steeper, with steepest slopes around −6.7. Note that this is not simply due to the entire stellar distribution having a steeper density profile (i.e. a vertical shift in the logarithmic slope profiles). In fact, the stars have similar slopes to the dark matter at smaller radii, and are only steeper by ∼ 0.5 − 1 dex (see also Schaller et al. 2015; Montes & Trujillo 2018; Pillepich et al. 2018) . Thus, although the stellar profiles are generally a bit steeper, the caustics are also more prominent. We also note that the spread of steepest slopes is larger for the stars; the dark matter profiles typically have slopes of −4.5 ± 0.6, while the stars have slopes of −6.7 ± 1.5.

In Fig. 4 , we consider the stacked density profiles of both the dark matter (solid black lines) and stars (long-dashed red lines). Here, all of the systems are stacked in the left panel, and the remaining panels show subsets of low (middle-left, Γ < 1.0), medium (middle-right, 1 < Γ < 2), and high (right, Γ > 2) mass accretion rates. In each panels we give the estimated caustic radius and associated uncertainty. Here, we use a bootstrap method (without replacement) to estimate the uncertainty in the caustic for the stacked profiles. Two things are immediately obvious from this figure: (1) The location of the "splashback" in dark matter coincides with the steepest slope of the stars, i.e. a "stellar splashback", and (2) the location and strength of this splashback radius, in both dark matter and stars varies with mass accretion rate: the caustic is stronger and located at smaller radii for higher mass accretion rates. We investigate these two key points in the following subsection. The square symbols show the (small number) of cases where a second caustic is identified in both the dark matter and the stars. The symbols are colour coded according to the recent mass accretion rate (Γ, see Section 3.1). The dotted line shows the one-to-one relation: the stellar caustics are located at almost the same radius as the dark matter. As seen in previous work, the caustics tend to be located at smaller radii when the recent mass accretion rate is higher.

The apparent coincidence between the dark mater and stellar splashback radius seen in the stacked profiles is compelling. However, in order to determine whether or not these two radii are really related, we need to compare each individual halo. This is shown in Fig. 5 , where the stellar caustics are shown against the dark matter caustics.

Here, the filled circle symbols indicate the splashback radii, and, for completeness, the filled squares show the second caustics. Note that we only show second caustics when one is robustly identified in both the dark matter and stars; this occurs in = 8 haloes (27%). In contrast, splashback radii in both stars and dark matter are found for all but one halo (CE-05 being the exception, which is currently undergoing a major merger). Remarkably, the stellar and dark matter caustics follow a tight one-to-one relation (with RMS scatter Δ ( / 200m ) = 0.11). The points in Fig. 5 are colour coded according to the mass accretion rate. Here, we can see the trend alluded to in the previous figures: the splashback is located at smaller radii for higher mass accretion rates. We look at this more explicitly below. Apart from being a more physically meaningful halo boundary, one of the most compelling reasons to probe the splashback radius in galaxy haloes is due to its strong link with mass accretion rate. Indeed, measurements of this radius can be used to classify galaxies by mass accretion rate, and can thus be used to probe aspects of halo formation like assembly bias (see e.g. More et al. 2016; Busch & White 2017) . In Fig. 6 , we show the location of the stellar (left) and dark matter (right) splashback radii as a function of Γ. The points are coloured according to the steepest slope at the caustic. The solid black line shows the predicted relation from Diemer (2020), assuming the median halo mass of the C-EAGLE sample (log 10 200m /M = 14.8). Our results are in good agreement with the Diemer (2020) predictions, and the stellar and dark matter caus- Haloes with higher mass accretion rates tend to have smaller splashback radii and steeper slopes. The solid black line shows the relation between splashback radius and mass accretion rate given by Diemer (2020) for the median halo mass of the C-EAGLE sample.

tics have 0.26 and 0.14 dex scatter in Caustic / 200m about fixed Γ, respectively. Perhaps most remarkable, however, is that the stellar caustics follow the Diemer (2020) trend (albeit with slightly larger scatter than the dark matter). Indeed, these results suggest that not only can detection of an outer caustic in the stars be used to define the physical boundary of the halo, the stellar splashback radius can also be used to measure the mass accretion rate when 200m is known. So far, we have only considered 3D density distributions. In reality, these are measured in projection, and thus to make connections with current and future observations we explore the projected density profiles in Section 4.

Before turning to the observational consequences of these theoretical results, it is worth discussing why we see these stellar splashback features in the simulations. The dark matter density profile, and the associated splashback radius, have been studied extensively in previous works. However, the corresponding stellar profiles have received much less attention. This is perhaps unsurprising: the hydrodynamical simulations required to form stars are far more expensive than dark matter only simulations, and, perhaps more importantly, include uncertain subgrid galaxy formation prescriptions. Deason et al. (2020) studied the edges of stellar haloes using high resolution simulations of Milky Way-mass galaxies. However, the contrast with the results for cluster-mass scales is striking! In particular, Deason et al. (2020) found that the stars did not generally reach out to the splashback radius of the dark matter, and, in fact, the edge of the Galactic-sized stellar haloes more often coincide with the second caustic of the dark matter. We suggest that this difference is mainly owing to three mass-dependent effects: (1) the stellar mass-halo mass relation, (2) the importance of smooth accretion, and (3) the formation age or concentration of the host halo.

First, the stellar mass-halo mass relation is non-linear and varies as a function of halo mass (Moster et al. 2010; Behroozi et al. 2013) . Milky Way-mass haloes accrete most of their mass from small subhaloes, which themselves have high dark matter fractions, and, in some cases, no stars at all. On the other hand, the diffuse light on cluster-mass scales is dominated by the remains of massive galaxies (∼10 10 − 10 12 M ), which form stars efficiently (see e.g. As for the 3D case, these closely follow a one-to-one relation. Conroy et al. 2007; Purcell et al. 2007; Puchwein et al. 2010) . This leads to Galactic stellar haloes being dominated by a small number of progenitors (see e.g. Deason et al. 2016) , while a larger number of progenitors contribute to the ICL. Thus, the stellar mass-halo mass relation can partly explain why the stellar density profiles of clustermass haloes are more strongly related to the underlying dark matter distribution (e.g. Montes & Trujillo 2018; Pillepich et al. 2018) . Second, the importance of smooth accretion, in the form of dark matter particles not bound to any halo, varies as a function of halo mass; smooth accretion is dominant on Milky Way-mass scales, but mergers dominate the mass-growth in clusters (e.g. Fakhouri & Ma 2010; Genel et al. 2010; Wang et al. 2011) . Moreover, the fraction of mass in substructures is much larger in clusters than galactic haloes (e.g. Gao et al. 2011) . Thus, there are many more massive objects losing stars in a cluster than there would be in a typical Milky Way-mass halo. Finally, the cluster-mass haloes tend to form later, and have lower concentration than Milky Way-mass haloes. The luminous satellites that are accreted more recently tend to deposit stars at larger radii (see e.g. Cooper et al. 2015) , and thus the stripped material can reach out to the splashback radius. Note, here we have discussed the main factors that we believe determine the location of the stellar edges on different mass scales. However, there are many other mass-dependent effects that could be important. For example, the relevance of pre-processed satellite galaxies (e.g. Bahé et al. 2019) , the (stellar and dark) density profiles of the disrupting satellites (e.g. Peñarrubia et al. 2008; Watson et al. 2012) , and the survival times of satellite galaxies (e.g Bahé et al. 2019 ).

We end this interlude by noting that the results found here for cluster-mass haloes and the previous Deason et al. (2020) work on Milky Way-mass scales, span a significant mass range, but further work is warranted to fill in the remaining "mass-gap". For example, does the stellar edge smoothly change with radius between the second caustic of the dark matter and the splashback radius, or is there a sudden transition at a particular mass scale? This, and other related questions encourage a separate, more extensive, study of stellar halo edges across a range of halo masses.

In this Section, we investigate the projected density profiles, which are more relevant for observational studies of clusters. Projected density profiles are constructed in a similar way to the 3D profiles. We use the same radial bin size and smoothing, and, by default, compute the projected density in each radial shell using the angular median method. In the angular median method used above, each radial bin in the 3D density profile is split into 50 evenly spaced solid angles, and the median value is computed. For the 2D projected profiles, the same number of angular bins are used, but these bins are angles instead of solid angles. As we are considering particles within a 4 200m spherical aperture, there can be projection effects at larger radii (where we are artificially running out of particles). However, we find that these effects are minimal within 2 200m , within which the outer caustics are typically located.

In Fig. 7 , we show the dark matter (left panel) and stellar (right panel) caustics for individual haloes in 2D projection versus 3D. The filled circles show the splashback radii (or outer caustic) and the filled squares show, where applicable, the second caustics. Three different coloured symbols are used to indicate projections computed along different axes (in this case, along , , or in the simulation box). The results for the different projections of a halo are connected with a solid vertical line. In some cases, these differences can be substantial. The dotted lines in the panels show the one-toone relation, but we also show a solid line which best describes the relation between the projected and 3D quantities where, caustic ∼ 0.9 caustic . Finally, in the right hand panel we show the projected caustics for stars versus the dark matter. Like the 3D cases, these caustics line up on the one-to-one line and are directly related.

Up to now, we have computed density profiles for individual clusters using an angular median method ). In Fig. 8 we show the stacked profiles when two different methods are used to compute the individual profiles: (1) the mean in each (projected) radial shell, (2) the angular median density in each (projected) radial shell; our fiducial method. Unsurprisingly, the caustics for both the dark matter and the stars are weaker when the mean density profile is used. Indeed, this is one of the reasons the angular median method was proposed by Mansfield ) are excluded. The caustics for both dark matter and stars are weaker when the mean profile is used. In addition, the influence of substructures is more pronounced in the stellar profiles, which significantly contribute to the star light at large radii.

as the mean values in radial shells are more affected by outliers. However, it is also apparent that the difference between the mean and angular median profiles is much more relevant for the stellar material. The stacked caustic is hardly identifiable with the mean profile, but is very prominent when the angular median profile is used. The reason for this becomes clear when we consider profiles with bright galaxies explicitly removed (shown with the red dashed lines). Here, we have excluded all star particles bound to a subhalo with Star > 10 9 . This is approximately the stellar mass limit for which cluster member galaxies can be masked in observations (e.g. Zibetti et al. 2005; Zhang et al. 2019) . When the bright galaxies are excluded the splashback feature is discernible, even when the mean method is used. However, the angular median profile is largely unchanged, and is still more prominent than the mean profile. This is because there are additional structures, such as streams, plumes and clouds, that can effect the derived density profile. Thus, using a technique such as the angular median method is essential in order to detect the edge of the stellar material. Note that stacking a significant amount of systems could help alleviate this problem (we only have 30 haloes to stack with C-EAGLE), but the method used to measure the density profiles will still affect how strong the derived signal, if any, is.

The dotted lines in Fig. 8 show the stacked profiles when all subhaloes identified by the algorithm (Springel et al. 2001 ) are omitted. The outer profiles are steeper when subhaloes are removed (see e.g. Fielder et al. 2020) , however, this makes a much larger difference to the stellar density profiles. This is because the fraction of stars bound in subhaloes at large radii is much greater than the fraction of dark matter (e.g. Gao et al. 2011; Pillepich et al. 2018 ). As we discussed earlier, the process of removing stars bound to satellite galaxies cannot be directly replicated in observations (and, importantly, the identification of subhaloes in the simulations is not perfect, e.g. Cañas et al. 2019) . Relatively bright galaxies can be identified, but the contribution of fainter satellites is either ignored, or roughly estimated by making assumptions about the satellite galaxy luminosity function and their radial distribution (e.g. Zibetti et al. 2005; Zhang et al. 2019; Sampaio-Santos et al. 2020) . Fortunately, when the angular median method is used, the stellar splashback feature is present with or without the contribution of satellite galaxies. However, in order to robustly compare the location of this splashback feature with observations (see following section), it will be desirable to perform mock observations where the cluster light profiles are computed in the same way as the data.

The results of the previous sections predict a well-defined "edge" at the outskirts of the stellar distribution of cluster-mass haloes. Recent work has shown that the ICL closely follows the total matter distribution (Montes & Trujillo 2019; Alonso Asensio et al. 2020) . In agreement with these results, we have now shown that the stars at the very outskirts of the cluster have a splashback radius, or radius of steepest slope, at the same location as the dark matter. This prediction begs the question: can this edge be observed in the ICL?

In Fig. 9 we show the surface brightness profile for the stacked sample of = 0.25 clusters from Zhang et al. (2019) . This profile was derived from = 300 galaxy clusters from the Dark Energy Survey (DES) Year 1 data, with median halo mass of 200m = 2.5 × 10 14 M . The shaded pink region shows their fiducial profile, while the gray line indicates the profile when the flux is defined to be zero at 1 Mpc (cf. zero flux at 1.8 Mpc for the fiducial profile). . The projected stellar mass density profiles of the C-EAGLE clusters have been converted into surface brightness assuming a constant mass-to-light ratio of / = 5, and then scaled to the same (average) halo mass as the DES sample. We have also included the (1 + ) 4 dimming factor applicable for = 0.25. We show the mean and the angular median stacked profiles with the dark blue and orange lines, respectively. The corresponding profiles when bright galaxies or all subhaloes are removed are shown with the dashed and dotted lines, respectively. Right: The logarithmic slope profiles of the projected stellar density profiles for the C-EAGLE clusters and the data. Identifying the caustic feature in the light distribution requires probing to larger projected distances, and thus fainter surface brightness limits ( ∼ 34 mag arcsec −2 at = 0.25). In addition, the caustic will only be detected if stacking methods take into account the influence of outliers, for example by using an angular median method.

This latter profile was made to directly compare with the Zibetti et al. (2005) results from a stack of SDSS clusters (shown with the purple line-filled region). The dark blue and orange lines show the (median) surface brightness profile for the C-EAGLE clusters using the mean, or angular median method to compute the density profiles. The dashed (dotted) lines show the profiles when bright galaxies (all subhaloes) have been explicitly removed. The projected stellar mass density profiles of the C-EAGLE clusters have been converted into surface brightness assuming a constant mass-to-light ratio of / = 5M /L (as in Alonso Asensio et al. 2020 ). In addition, we have scaled the derived intensity profiles to the median halo mass of the DES sample surface mass density, assuming the stellar mass of the ICL scales with halo mass (see e.g. Pillepich et al. 2018 ). Finally, we have included the (1 + ) 4 dimming factor to take into account the = 0.25 redshift of the observed clusters. Note that the absolute surface brightness level of the C-EAGLE clusters should be taken with a pinch of salt, as this depends on the agreement with the observed stellar masses and the exact way in which the ICL is defined. Nonetheless our simple conversion of projected surface density to surface brightness shows reasonable agreement with the observed data, and gives a good indication of the surface brightness levels needed to probe to the cluster edges.

In the right-hand panel of Fig. 9 we show the logarithmic slope profiles of the projected stellar density profiles. The lines from C-EAGLE are taken from Fig. 8 . We show the two profiles from Zhang et al. (2019) using different flux offsets with the filled pink and line-filled gray regions, respectively. For the observed profile, we scale relative to the median 200m of the sample. The errors in the observed fluxes are taken into account by computing the slope profiles many ( = 10 4 ) times and scattering the values according to the error distributions. Here, we ignore bins where the errors are particularly high (greater than 10%). In addition, we rebin the observed data to have = 21 logarithmic radial bins between −1.0 < log ( / 200m ) < −0.3. Fig. 9 that, regardless of the flux offset, the DES data is consistent with a constant slope. However, even though the data reaches to an impressive 1 Mpc, in order to probe the predicted stellar splashback, the data would need to go out to at least 2 Mpc -or surface brightness levels of ∼ 32 − 36 mag arcsec −2 . Although this appears to be unfeasible with current observations, this is certainly achievable with future observations. The Vera C. Rubin Observatory Legacy Survey of Space and Time (LSST, Ivezić et al. 2019) will survey an area that is > 10 times larger than the DES Y1 footprint, and is predicted to reach a depth of at least 2 magnitudes deeper. With the increase in both area and depth, and hence an increase in ICL flux collection by a factor of ∼ 100, LSST is capable of reaching these extreme projected distances and surface brightness levels. In the shorter term, the increased area and depth from the final DES data release may also provide some constraints on the stellar splashback; i.e. the DES data may be able to rule out such a feature if the profile does not show the predicted sharp drop. It is also worth considering the methods used to create the stacked ICL profiles. For example, while care is often taken to remove bright galaxies from the profiles, further improvement could be made by applying an angular median method on each individual cluster, as used in this work.

In addition to stacked profiles, there is also scope to probe these extreme outskirts for individual clusters. The Beyond Ultradeep Frontier Fields and Legacy Observations (BUFFALO) Hubble Space Telescope treasury program (Steinhardt et al. 2020) can probe the ICL of individual clusters out to ∼1 Mpc, and an extension of this program could push to even larger distances. For individual clusters, care must be taken to avoid confusion between the stellar splashback with the second caustic (this secondary feature is normally washed out in stacked profiles). However, these secondary features tend to be less prominent amongst haloes with very high (recent) mass accretion rates, and the BUFFALO clusters are particularly active systems. Finally, the depth and field-of-view of upcoming spacebased missions such as Euclid (Laureijs et al. 2011 ) and the Nancy Grace Roman Space Telescope (Spergel et al. 2015) will almost certainly push the boundaries of low-surface brightness science to reach these extreme stellar edges.

We have used the C-EAGLE suite of simulations to explore the outskirts of dark matter and stars on cluster-mass scales. Density profiles of each C-EAGLE system are constructed using an angular median method, which limits the influence of substructure and other non-diffuse components. The outer caustics, or splashback radii, of both the dark matter and stellar components are compared, both for individual clusters and with stacks based on mass accretion rate. Our main conclusions are summarised as follows:

• The stellar density profiles of clusters have a well-defined edge, defined by the radius of steepest log-density slope, which coincides with the outer dark matter caustic, or the splashback radius. This radius is typically located at 200m , in good agreement with previous work using dark-matter-only simulations.

• The location of the stellar and dark matter splashback radius depends on the mass accretion rate: slowly accreting haloes tend to have an edge at a larger radius, and a shallower steepest slope. The stellar profiles have more prominent outer caustics than the dark matter. In some cases (∼27%), a secondary caustic can be identified in the stellar and dark matter profiles: these are likely the result of earlier, massive mergers, but the features are much weaker than the radius of steepest slope, and hence more difficult to detect.

• The radius of steepest slope can also be identified in projection, where the 2D and 3D radii are related by Caustic ∼ 0.9 caustic . The method used to identify the caustic is crucial, as massive substructures can significantly dilute the signal of the steepest slope. This is especially true for the stellar material: there is a higher fraction of stars than dark matter bound in subhaloes (see e.g. Pillepich et al. 2018 ).

• Current observations of the intracluster light can reach out to ∼ 1 Mpc, either for individual systems or from stacking many systems. Detecting the stellar splashback will require probing out a further 1 Mpc, to surface brightnesses of ∼32 − 36 mag arcsec −2 . However, this challenging feat will be achievable with upcoming facilities ideally suited to low surface brightness studies, such as the Vera C. Rubin Observatory (see e.g. Brough et al. 2020) , the Nancy Grace Roman Space Telescope, and Euclid.

The relation between the visible and dark matter is complex, mass-dependent and depends on galaxy formation physics, hierarchical structure formation and cosmology. The stellar haloes of galaxies and clusters offer a unique way to probe the dark matter, as these are mainly built from mergers. This work shows that measuring the "edge" of the intracluster light offers an alternative way to define the halo boundary and quantify the mass accretion rate of clusters. Moreover, by comparing the stellar splashback with independent measures of the splashback radius, e.g. from satellite galaxies or weak lensing, we can quantify the link between the stellar and dark material, and thus test the predictions of our galaxy formation models. Ultimately, learning how the outskirts of dark and stellar haloes change with mass and time will provide an invaluable way to critically examine, and inform, our state-of-the-art cosmological models of structure formation.

The data presented in the figures are available upon request from the corresponding author. The raw simulation data can be requested from the C-EAGLE team (Bahé et al. 2017; Barnes et al. 2017) .

The General Assembly of Galaxy Halos: Structure, Origin and Evolution

AD thanks Stefano Zibetti for kindly sharing his surface brightness data. The authors also thank Anthony Gonzalez for providing useful comments that have improved this manuscript. Our gratitude is extended to all of the essential workers that support our livelihood, especially during the Coronavirus pandemic. Last, and certainly not least, AD thanks the staff at the Durham University Day Nursery who play a key role in enabling research like this to happen.AD is supported by a Royal Society University Research Fellowship, and by the Science and Technology Facilities Council