The mass–size relation of galaxy clusters

The outskirts of accreting dark matter haloes exhibit a sudden drop in density delimiting their multistream region. Due to the dynamics of accretion, the location of this physically motivated edge strongly correlates with the halo growth rate. Using hydrodynamical zoom-in simulations of high-mass clusters, we explore this deﬁnition in realistic simulations and ﬁnd an explicit connection between this feature in the dark matter and galaxy proﬁles. We also show that the depth of the splashback feature correlates well with the direction of ﬁlaments and, surprisingly, the orientation of the brightest cluster galaxy. Our ﬁndings suggest that galaxy proﬁles and weak-lensing masses can deﬁne an observationally viable mass-size scaling relation for galaxy clusters, which can be used to extract cosmological information.


I N T RO D U C T I O N
In the CDM paradigm, structure in the Universe arises from the initial density perturbations of an (almost) homogeneous dark matter distribution. Due to gravitational evolution, this leads to the appearance of collapsed structures, i.e. dark matter haloes. Some of the baryonic matter, following this process, cools down and settles at the centres of the gravitational potentials where it forms galaxies.
This mechanism has been studied through models of so-called spherical collapse (Gunn & Gott 1972;Bertschinger 1985), whose main prediction is the existence of a radius within which the material orbiting the halo is completely virialized. In general, this virial radius depends on cosmology and redshift, but both in numerical simulations and observations, fixed overdensity radii are widely used as proxies for this quantity. An example of this is r 200 m , defined as the radius within which the average density is 200 times the average matter density of the Universe, ρ m . The corresponding enclosed mass is known as M 200 m .
Halo mass functions constructed with these idealized definitions can capture the effects of cosmology (Press & Schechter 1974), the nature of dark matter (Angulo, Hahn & Abel 2013), and dark energy (Mead et al. 2016) on the growth of structure. In the real Universe, however, this picture is complicated by the triaxiality of haloes (Dubinski & Carlberg 1991;Monaco 1995) and the existence of clumpy (baryonic) substructure (Bocquet et al. 2015).
Because the process of structure formation is hierarchical, massive haloes contain subhaloes, some of which host galaxies themselves. The resulting clusters of galaxies are the focus of this work. What makes these objects particularly unique is the fact that they are not fully virialized yet. To this day, they are still accreting both ambient material and subhaloes through filamentary structures E-mail: contigiani@strw.leidenuniv.nl surrounding them (Bond, Kofman & Pogosyan 1996). Because of their definition, however, traditional overdensity definitions of mass are not only affected by halo growth, but also by a pseudo-evolution due to the redshift dependence of ρ m (Diemer, More & Kravtsov 2013). Diemer & Kravtsov (2014) and More, Diemer & Kravtsov (2015) were the first to note that this growth process leads to the formation of a sharp feature in the density profile that separates collapsed and infalling material. This feature therefore defines a natural boundary of the halo. The location of this edge, i.e. the splashback radius r sp , has an obvious primary dependence on halo mass, but also a secondary dependence on accretion rate. While this behaviour can be qualitatively explained using simple semi-analytical models of spherical collapse, none of the analytical models currently proposed (Adhikari, Dalal & Chamberlain 2014;Shi 2016) can fully describe its dependency on mass and accretion rate . Despite this, the corresponding definition of halo mass is particularly suited to define a universal mass function valid for a wide range of cosmologies (Diemer 2020a).
In this paper, we try to bridge the gap between the theoretical understanding of the splashback feature and observational results, both past and future. The outer edge of clusters has already been extensively measured through different tracers: the radial distribution of galaxies from wide surveys (More et al. 2016;Baxter et al. 2017;Chang et al. 2018), but also their velocity distribution (Tomooka et al. 2020;Fong & Han 2021), and in the weak-lensing signal of massive clusters (Umetsu & Diemer 2017;Chang et al. 2018;Contigiani, Hoekstra & Bahé 2019b). Furthermore, forecasts have already set expectations for what will be obtainable from nearfuture experiments (Fong et al. 2018;Xhakaj et al. 2019;Wagoner et al. 2020). Despite the wealth of data and studies, however, not many splashback observables have been proposed. The only robust application of this feature found in the literature is related to the study of quenching for newly accreted galaxies (Adhikari et al. 2020). To achieve our goal, we make use of hydrodynamical simulations of massive galaxy clusters, which we introduce in Section 2. We focus mainly on z = 0, but also make use of snapshots at redshifts z = 0.474 and z = 1.017. In Section 3, we start our discussion by introducing the physical interpretation of splashback and consider the connection between the galaxy and dark matter distributions. We then continue in Sections 4 and 5, where we explain how galaxy profiles and weak-lensing mass measurements can be combined to construct a mass-size relationship for galaxy clusters. Finally, we summarize our conclusions and suggest future developments in Section 6.

H Y D R A N G E A
The Hydrangea simulations are a suite of 24 zoom-in hydrodynamical simulations of massive galaxy clusters (log 10 M 200 m /M between 14 and 15.5 at redshift z = 0) designed to study the relationship between galaxy formation and cluster environment (Bahé et al. 2017). They are part of the Cluster-EAGLE project (Bahé et al. 2017;Barnes et al. 2017) and have been run using the EAGLE galaxy formation model (Schaye et al. 2015), which is known to reproduce galaxy observables such as colour distribution and star formation rates. To better reproduce the observed hot gas fractions in galaxy groups, the AGNdT9 variant of this model was used (Schaye et al. 2015).
The zoom-in regions stretch to between 10 and 30 Mpc from the cluster centre, roughly corresponding to ࣠10r 200 m . For the definition of the cluster centre, in this work, we choose the minimum of the gravitational potential. We note, however, that this choice will not impact our conclusions since we will focus on locations around r 200 m . The particle mass of m ∼ 10 6 M for baryons and m ∼ 10 7 M for dark matter allows us to resolve galaxy positions down to stellar masses M * ≥ 10 8 M and total masses M sub ≥ 10 9 M , respectively.
In Fig. 1, we show the log-derivative of the stacked subhalo density n s (r) at large scales. This is the result of a fit obtained using the model of Diemer & Kravtsov (2014), and we refer the reader to the aforementioned paper and Contigiani et al. (2019b) for a detailed explanation of the model and its components. The choice to employ this profile is based on its ability to capture the sharp feature visible Table 1. The Hydrangea clusters used in this paper and their z = 0 properties. 0.3 is the accretion rate measured between z = 0 and z = 0.297. The three splashback radii r sp , r g sp , r s sp refer to the splashback radius measured, respectively, in the dark matter, galaxy, and subhalo distributions (see Section 3). For two clusters, CE-28 and CE-18, the radius r sp is not used in this work because the dark matter distribution displays a featureless profile at large scales. All quantities are in physical units. around r 200 m , which is the focus of this work. We optimally sample its 8D parameter space using an ensemble sampler (Foreman-Mackey et al. 2013).
In the same plot, we also include the stacked subhalo profile of the accompanying dark matter only simulations, initialized with matching initial conditions. The two profiles match almost exactly, suggesting that baryonic effects do not alter this feature to a significant extent (see also O'Neil et al. 2020). While not shown, we report that the same conclusion can be reached by looking at the full matter distribution ρ(r) in the two sets of simulations. Similarly, this feature is also visible in the number density of galaxies, n g (r). Due to our focus on all three of these profiles, we choose not to work with background subtracted quantities.
For reference, we present a full list of the simulated clusters used in this paper and their relevant properties, some of them defined in the following sections, in Table 1.

Definition
For haloes that continuously amass matter, material close to its first apocentre piles up next to the edge of the multistream region, where collapsed and infalling material meets (Adhikari et al. 2014). A sudden drop in density, i.e. the feature visible in the profiles of Fig. 1, is associated with this process.
This intuitive picture leads to three characterizations of the splashback radius, depending on the approach used to measure or model it: While these definitions have all been previously hinted at in the introduction, in this section, we explicitly present them and discuss the connections existing between them. This also justifies our adopted definition, based on the density profile.
The first definition is clearly motivated in the spherical case but fails once it is applied to realistic haloes. The presence of angular momentum and tidal streams from disrupted subhaloes (see e.g. Vogelsberger & White 2011), smooth out this feature and make its description murky. The second definition was the first suggested in the literature. Introduced by Diemer & Kravtsov (2014), it is based on the study of dark matter profiles in N-body simulations and has been linked to the first, more dynamical, definition (Adhikari et al. 2014;Shi 2016). The third was first suggested by Diemer (2017), who showed that this location can be calibrated to the second one ) by choosing specific percentiles of the apocentre distribution.
To clarify the relationship between the outermost caustic and apocentre, it is educational to use a self-similar toy model based on Adhikari et al. (2014) to show the phase-space distribution of a constantly accreting halo with an NFW-like mass profile (Navarro, Frenk & White 1997).
In the absence of dark energy, we follow the radial motion of particles, between their first and second turnaround in the mass profile: We impose that the total mass evolves as M(R, t)∝t 2 /3 , R∝t 2(1 + /3)/3 , and the dimensionless NFW profile is defined as: In this set of equations, is the dimensionless accretion rate, R represents the turnaround radius, and the scale parameter r s is defined by the infall boundary condition This condition, combined with the turnaround dynamics, imposes that M(R, t)∝(1 + z) − (Fillmore & Goldreich 1984). We point out that the dependence on the time-sensitive turnaround properties M(R, t), R(t) can be factored out from the equations above, meaning that the entire phase-space at all times can be obtained with a single numerical integration.
In Fig. 2, we show the result of this calculation, denoting the location of the outermost caustic as r c sp . The caustic is formed by the outermost radius at which shells at different velocities meet (r/r c sp = 1 in the plot) and the location of shells at apocenter is defined by the intersection between the zero-velocity line and the phase-space distributions. From the figure, two things are noticeable: material at r c sp has not reached its apocentre yet, and the ratio between these two locations depends on the accretion rate. It is beyond the scope of this work to quantify this dependence since it depends heavily on the mass profile inside r c sp . Qualitatively, however, the difference between caustic and apocentre is easy to understand once the dynamical nature of this feature is considered: the halo is growing in size, and while some material is now reaching its apocentre, mass accreted more recently has the chance Figure 2. The phase-space structure of accreting dark matter haloes depends on the accretion rate . We employ a toy model of spherical collapse to describe the multistream region of NFW-like haloes. The figure shows that the material at the outermost caustic, r c sp , is not necessarily at apocentre (where v = 0) and that the ratio of these two radii is a function of accretion rate. For ease of readability, we have rescaled the coordinates by r c sp , and the velocity of collapsed material at this point.
to overshoot it and form the actual caustic. In a static picture, this would not be the case.
In realistic haloes, this dependence on accretion rate is only one of many factors that biases and adds scatter to the relationship between the halo boundary and apocentres. Other factors include, e.g. nonspherical orbits and the presence of multiple accretion streams. Despite this, Diemer (2017) has shown that there is a clear link between the apocentre distribution and splashback. The percentile definition introduced there is particularly suited to theoretical investigations, but its usefulness in the very low-regime is still uncertain Xhakaj et al. 2019), and it has not been explored in the presence of modifications of gravity (Adhikari et al. 2018;Contigiani, Vardanyan & Silvestri 2019a).
For this work, we define the splashback radius as the location of the steepest slope as defined by a profile fit. In Table 1 we report, for each cluster, this radius measured in the distribution of galaxies, subhaloes, and total matter (r g sp , r s sp , r sp ). The model is a modified Einasto profile (Einasto 1965) with the addition of a power law to take into account infalling material (Diemer & Kravtsov 2014). Regarding the goodness of fit, we find that up to and around r 200 m the standard deviation of the residuals is of order 10 per cent. On the other hand, the presence of substructure superimposed on a shallow density profile results in normalized residuals of order 50 per cent in the outer regions.
To further justify our approach, we show in Fig. 3 how this simple definition of splashback radius is able to capture the phase-space boundary of different haloes, even when a sudden drop in density is absent. The main benefit of this definition is that it avoids the arbitrariness of the apocentre definition, or the bias induced by multiple caustics in the minimum slope definition ). Its main caveats, however, are that (1) it is computationally expensive since it requires high-resolution simulations and a multiparameter fit procedure, and (2) it might not apply to low-mass clusters and galaxy groups. We leave this last question open for future investigations.
We wrap this subsection up by stressing that this definition of 'the' splashback radius is, like any other, useful only to study its correlation with other properties, or quantify the impact of different . Fitting simulated subhalo profiles with a smooth model. In the top panels, we show the radial subhalo distributions of two clusters (CE-16, left and CE-9, right), together with the best-fitting profiles used to reconstruct the log-derivative. In the bottom panels, we show how the inferred location of the log-derivative minimum (vertical line) identifies the phase-space edge of relaxed (left) and perturbed (right) galaxy clusters. In the phase-space plots, the cluster on the left is formed by collapsed particles, while the stream visible on the lower right is infalling material. The right-hand panels demonstrate how our approach is effective even in the presence of an on-going merger when the splashback feature is not visible as a sharp transition in the density profile. physical processes. While the flexibility of the chosen model is not surprising given the number of free parameters, the clear connection between the phase-space and the log-derivative in individual haloes is a powerful and seemingly general result. Ultimately, however, the observational results focus on stacked projected density profiles, and so should the predictions.

Accretion
It is well established (Diemer & Kravtsov 2014;More et al. 2015;Mansfield et al. 2017;Diemer 2020b) that the location of the halo boundary correlates with the accretion rate In this work, this ratio is calculated in the redshift range z = 0 to z = 0.293, since this time interval roughly corresponds to one crossing time for all clusters considered here, i.e. how long ago the material currently at splashback has been accreted (Diemer 2017). Although this choice is partially arbitrary, we have investigated the dependence of our results on the redshift upper limit and we have verified that our main conclusions are not affected. The archetypical relation demonstrating this idea is plotted in Fig. 4, where we have also included the relations found in More et al. (2015), Diemer et al. (2017), and Diemer (2020b), to provide additional context. We find good agreement, even though a perfect  . The distribution of subhaloes and galaxies outside the cluster edge as a function of accretion rate. Faster growing haloes display a more concentrated distribution of satellites outside of their boundary. This behaviour seen in individual clusters is not explained by simple models of spherical collapse (blue shaded area), but the average profile (marked by a star) matches the expectation. This suggests that non-isotropic processes shape this relation. match is not necessarily expected. The Hydrangea clusters represent a biased sample, selected to be mostly isolated at low redshift (Bahé et al. 2017). While the effect of this selection on the accretion rate distribution is not fully known, we show below that a connection between cluster environment and this quantity exists, and the presence of mergers might therefore influence it. This is not surprising since a connection between accretion and large-scale bias is already known (e.g. Fakhouri & Ma 2010).
We show this relationship explicitly in Fig. 5 by using one of the parameters of the profile model. As visible in the figure, the power-law index of material outside of splashback correlates with the accretion rate. We find that this is true for both subhaloes and galaxies and that the difference between the two is consistent with sample variance.
To try and explain this behaviour, we use a fully consistent model of spherical collapse introduced by Bertschinger (1985), which was also used in Contigiani et al. (2019a). The setup of this toy model is the same as what is shown in equation (1), but with a mass profile that also needs to be solved for. Starting from an initial guess for M(r, t) = M(r/R(t)), orbits are integrated and their mass distribution is calculated. Iterating this process multiple times returns a self-similar density profile and orbits consistent with each other.
The result of this calculation is also shown in Fig. 5. Because the mass-profile prediction is not a power law, we plot a filled line displaying the range of logarithmic slopes allowed between r c sp and 2r c sp . The fact that this prediction is not a function of accretion rate implies that the correlation between the slope and the accretion rate seen in the simulations is not purely dynamical, and suggests a connection between the cluster environment and accretion rate.
We stress here that previous splashback works have mostly focused on stacked halo profiles, for which the expectation of the spherically symmetric calculation shown above is roughly verified, even in the presence of dark energy (Shi 2016). We also recover this result for our sample (see the star symbol in Fig. 5), but we point out that this is a simple conclusion. Because Newtonian gravity is additive, stacking enough clusters should always recover the spherically symmetric result. Despite this, we also note that results from the literature do not always agree with this prediction. However, we do not linger on these discrepancies since (1) this was never the focus of previous articles, and (2) different methods to extract the power law have been employed.

Anisotropy
This departure from the spherical case implies that anisotropies play a role in shaping the accretion rate . To study the impact of accretion flows on the cluster boundary, we study 72 sky projections of the Hydrangea clusters (3 each, perpendicular to the x, y, and z axes of the simulation boxes) and rotate them to align the preferred accretion axes in these planes. For each projection, we define this direction θ ∈ (− π /2, π /2) in two ways: (1) to capture the filamentary structure around the cluster between r 200 m and 5r 200 m , we divide the subhalo distribution in 20 azimuthal bins and mark the direction of the most populated one, and (2) to capture the major axis of the BCG, we use unweighted quadrupole moments of the central galaxy's stellar profile within 10 kpc from its centre. The mean projected distributions according to these two methods are presented in the left and right top panels of Fig. 6, respectively.
Looking at the top-left panel of the figure, it is not surprising that filamentary structures of the cosmic web are visible around the central cluster -this is by construction. Because of the higher contrast between outside and inside regions, the subhalo distribution exhibits a sharper feature in the directions pointing towards voids (see central panel of Fig. 6). More surprisingly, however, these same traits are also noticeable in the mean distributions aligned according to the central galaxy's axis (see lower panel).
This result implies that the distribution of stellar mass within the central 10 kpc of the cluster contains information about the distribution of matter at radii which are a factor 10 2 larger. In fact, the connection between the shape of the dark matter halo and the ellipticity of the brightest cluster galaxy (BCG, which is also the central galaxy for massive galaxy clusters) is known (Okumura, Jing & Li 2009;Herbonnet et al. 2019;Ragone-Figueroa et al. 2020). And, similarly to other results (Conroy, Wechsler & Kravtsov 2007;De Lucia & Blaizot 2007), the Hydrangea simulations predict that the stellar-mass buildup of the BCG is driven by the stripping of a few massive satellites after their first few pericentre passages (Bahé et al., in preparation). Because these galaxies quickly sink to the centre, the material they leave behind is therefore a tracer of their infalling direction.

T H E M A S S -S I Z E R E L AT I O N
In our sample, we find that the splashback feature seen in the galaxy, subhalo, and total matter profiles are all at the same location. The mean fractional difference between any two of r g sp , r s sp , or r sp is consistent with zero, with a mean standard deviation of 3 per cent. We also verified that this statement is unaffected by cuts in subhalo mass or galaxy stellar mass. Due to the limited size of our sample, the effects of dynamical friction on the distribution of high-mass subhaloes are not visible (Adhikari, Dalal & Clampitt 2016).
We emphasize, however that this does not mean that galaxy selection effects have no impact on these quantities. For example, it is an established result, both in the Hydrangea simulations (Oman et al. 2020) and in observations (Adhikari et al. 2020) that the location of a galaxy in projected phase-space correlates with its colour and star-formation rate. This is because a red colour preferentially selects quenched galaxies that have been orbiting the halo for a longer time.
Until their first apocentre after turnaround, galaxies act as test particles orbiting the overdensity as the halo grows in mass. In the standard cold dark matter paradigm, based on a non-interacting particle, it is not surprising then that the edge formed in their distribution is identical to the one seen in the dark matter profile. It  (5) obtained for β = 0 and how modifications of gravity are expected to affect this relation (blue dashed line, see text for more details). In this relation, a secondary dependence on the accretion rate 0.3 is a source of scatter that can be captured if β = 0. As visible in the residuals in the bottom panel, a simple power-law form well reproduces this dependence. In the considered sample, we find that half of the total scatter (0.25 dex) is due to the mass accretion rate distribution.
should be noted, however, that this is not necessarily true in extended models in which dark matter does not act as a collisionless fluid. Due to their infalling trajectories, the distribution of galaxies will always display a splashback feature, even if the dark matter profile does not exhibit one.
In the cold dark matter scenario, our result implies that galaxies can be used to trace the edge of clusters. We note, in particular, that this measurement has already been performed several times using photometric surveys (Baxter et al. 2017;Nishizawa et al. 2017;Chang et al. 2018;Shin et al. 2019;Zürcher & More 2019). Furthermore, due to the large number of objects detected, galaxy distributions obtained through this method offer the most precise measurements of splashback. The accuracy of the results, however, depends heavily on the details of the cluster finding algorithm (Busch & White 2017;Shin et al. 2019).
With this in mind, we build an observational mass--size relation between the location of this feature in the galaxy distribution (r g sp ) and the mass enclosed within it (M g sp ). In Fig. 7, we present the correlation between the two for the Hydrangea clusters. Because the splashback radius is roughly located at r 200 m (see Fig. 4), this relationship can be understood as a generalization of the virial mass-radius relation, where we have introduced a dependence on accretion rate. Surprisingly, we find that the dependence on 0.3 is well captured by a simple form: While we do not constrain β precisely, we find that β = 1.5 provides an adequate fit by reducing the total scatter from 0.25 dex to about half of this value. This choice of exponent and functional form is supported by the model of self-similar collapse used for Fig. 5, where we find that a power-law β = 1.45 fits this relationship with the same precision as the exponential functions calibrated to numerical simulations shown in Fig. 4. For a more extensive comparison with these predictions, we refer the reader to Section 6. The virial relation is a trivial connection between the mass and size of haloes based on an overdensity factor, but its observational power is limited by the fact that these masses are usually extracted from parametric fits to weak-lensing profiles that do not extend to the respective overdensity radii. Because of this, the overdensity masses have a strong dependence on the assumed mass-concentration relation (see e.g. Umetsu et al. 2020). The splashback feature, on the other hand, naturally predicts a mass-size relation for galaxy clusters and does so without the need for external calibrations.
In Fig. 7, we also plot the expected change in this relation to due modifications of gravity. We use the symmetron gravity model of Contigiani et al. (2019a) with parameters f = 1 and z ssb = 1.5, and assume that the change affects only the splashback radius and not the mass contained within it. The exact result depends on the theory parameters, but the expected change in this relation is around 0.15 dex.
Experimentally, we argue that this relation can be probed using a combination of galaxy density profiles (to extract r g sp ) and weaklensing measurements. Aperture masses (Clowe et al. 2000), in particular, can be used to extract in a model-independent fashion the average projected mass within a large enough radius. If necessary, the aperture mass can also be deprojected to obtain a low-bias estimate (Herbonnet et al. 2020).

R E D S H I F T E VO L U T I O N
So far, we have only considered the simulation predictions at z = 0. In this section, we extend our analysis to higher redshifts by exploring two other snapshots of the Hydrangea simulations at z = 0.474 and z = 1.017.
At these higher redshifts, we find that the scatter in the splashback relation for individual haloes is large. This is visible in Fig. 8, where we plot the equivalent of Fig. 4 for these two snapshots. We recover the general result of Diemer (2020b) that the average values of r sp /r 200 m and should be higher at early times, but the correlation between the two is completely washed out by z = 1. We connect this to three causes: (1) The fixed time interval between the snapshots does not allow us to reliably estimate at higher redshift when the crossing times are smaller. (2) The lower number of resolved galaxies and subhaloes means that the residuals of the individual profile fits are larger around the virial radius. And finally, (3) the higher frequency of mergers at high redshift means that the numbers of haloes with profiles not displaying a clear splashback feature increases.
We find that equation (5) is still valid, even if our ability to constrain the scatter at high redshift is impeded by the sample variance. Furthermore, we report that the splashback overdensity M sp /r 3 sp has a redshift dependence. Or, in other words, that the logarithmic zero-point that was not specified in equation (5) is a function of redshift. Not accounting for the dependence, our best-fitting values for the logarithm of the average overdensity log 10 (M sp /M ) − 3log 10 (r sp /Mpc) are [13.3, 13.8, 14.1] ± 0.3 at redshifts [0, 0.5, 1].
Regarding the anisotropy in the splashback feature due to filamentary structures, we report that this phenomenon exists also at high redshift. In Fig. 9, we compare the sky-projected subhalo profiles s (R) towards different directions, similarly to what we have done for Fig. 6. In this case, however, we explicitly discuss the connection with Figure 8. The splashback radius and its correlation with the accretion rate as a function of redshift. This plot is an extension of Fig. 4 for redshifts z = 0.474 (orange crosses) and z = 1.017 (light blue plus symbols). The ratio between the splashback radius and the 200 m overdensity radius should correlate with the accretion rate , but for the Hydrangea snapshot at z = 1.017, the large sample variance washes out this correlation. Despite this, we still recover the expectation of previous results (plotted lines), a larger average r sp /r 200 m at higher redshift. Figure 9. The impact of filaments and accretion flows on the outer density profile of massive haloes as a function of redshift. We plot the mean value and variance of the ratio between the 2D subhalo distributions in quadrants perpendicular ( ⊥ ) and parallel ( ) to the accretion direction defined through two tracers. This ratio is closely related to what can be measured in observations. While the difference in profile towards and away from filamentary structures is visible at all redshifts, the orientation of the central galaxy is not a good tracer of the splashback anisotropy at high redshift. This is because the central BCG is still forming and its orientation is not yet finalized.
observations by plotting directly the ratio of the density profiles inside quadrants oriented towards and perpendicular to the two accretion directions, instead of focusing on the result of the profile fits. The mean and variance of these ratios are calculated assuming that the different projections are independent. We find that the orientation of the major axis of the brightest cluster galaxy does not correlate with a splashback anisotropy at z = 1. This is because, in most cases, the identification of a central, brightest galaxy is not straightforward at this redshift. At early times, the future central galaxy is still in the process of being created from the mergers of multiple bright satellites located close to the cluster's centre of potential.
To conclude this section, we point out that in the region around r 200 m , the difference between the profiles perpendicular and parallel to the central galaxy's major axis is about 10 per cent at redshift z ࣠ 0.5. This departure is well within the precision of galaxy profiles extracted from large surveys (e.g. Adhikari et al. 2020). Therefore, this measurement might already be possible using such catalogues.

D I S C U S S I O N A N D C O N C L U S I O N S
On its largest scales, the cosmic web of the Universe is not formed by isolated objects, but by continuously flowing matter distributed in sheets, filaments, and nodes. For accreting (and hence non-virialized) structures such as galaxy clusters, the splashback radius r sp represents a physical boundary motivated by their phase-space distributions. To exploit the information content of this feature, in this paper, we have introduced and studied two observable quantities related to it.
First, we have shown that the full galaxy profile can be used to define a cluster mass, i.e. the mass within r sp . This is an extension of the traditional approach of using richness as mass proxy (see e.g. Simet et al. 2016). Because of the dynamical nature of the equivalent feature in the dark matter profile, we conclude that, observationally, the splashback feature in the galaxy profile defines the physical halo mass. Moreover, we have shown here that the natural relation between the mass and size of haloes according to this definition (see Fig. 7) can be used to constrain new physics at cluster scales. Because this boundary is delimited by recently accreted material, we found that a majority of the scatter in this mass-size relation can be explained through a secondary dependence on accretion rate .
Secondly, we have explored how this connection to the accretion rate might be interpreted as a connection between the geometry of the cosmic web and how clusters are embedded in it. The relation between the two is made explicit in Figs 5, 6, and 9. In these figures, we have investigated how the cluster environment affects both the halo growth and the stellar distribution of the central galaxy. This information, combined with the scatter of the mass-size relation, can therefore be used as a consistency check for any property that claims to select for accretion rate.

The role of simulations
In the last few years, the study of the splashback feature has evolved into a mature field both observationally and theoretically. We use this section to discuss explicitly the connection between the two, in light of this work and its connection to previous endeavours.
In the context of splashback, simulations have guided the formulation of theoretical principles and hypotheses. However, as more measurements become viable, it becomes necessary to provide clear and powerful observables. Following this spirit, we used high-resolution hydrodynamical simulations to explore directly the connection between measurements based on sky-projected galaxy distributions and theoretical predictions.
Our conclusions regarding the mass-size relation and its redshift evolution are similar to the results of Diemer et al. (2017) and Diemer (2020b), which are based on more extensive N-body simulations. For the sake of completeness, it important to note that, in the same papers, it was also found that the splashback overdensity is not universal, but has both a mild dependency on M 200 m and a strong dependency on cosmology, especially at low redshift (z ≈ 0.2). Due to our limited sample, we are clearly unable to model these effects in this work. None the less, we point out that our goal here is to construct a pure splashback scaling-relation based on galaxy profiles and weaklensing mass measurements. Every other dependency, if present, should be captured either as additional scatter or through different parameter values.
We also point out that that these previous works are based on the apocentre definition of splashback (see Section 3). In contrast, we defined the splashback radius as the point of steepest slope according to a model fit to the density profiles of galaxy clusters. While we do not necessarily expect the two definitions to differ, our choice is based on its connection to observations, and the desire to highlight the fact that the splashback radius is not only some abstract halo property but can be defined as a characteristic of individual profiles, such as, e.g. the concentration parameter (Navarro et al. 1997).
An alternative method employed by other studies (Diemer & Kravtsov 2014;Mansfield et al. 2017;Xhakaj et al. 2019) to obtain a measure of r sp makes use of the minimum of the logarithmic slope in smoothed profiles. While this approach is much faster than profile fitting when only r sp is of interest, it does not describe the full shape visible in Fig. 1. In particular, a model that captures the width of this feature is necessary to define the slope of the outer region without an arbitrary choice of which radial scales to consider. Because the model used here contains an asymptotic outer slope, this definition is unique.
Our decision has, of course, its drawbacks. The versatility of the fitted model is necessary to capture the variance of the individual profiles, but the resulting intrinsic scatter is large and not the best suited to study tight splashback correlations (such as Fig. 4). At the same time, the large parameter space might also be seen by some as a chance to study a multitude of correlations between different model parameters. However, we resist this temptation, as inferences based on such correlations might say more about the particular model employed than provide any physical information.
A subtler difference between our method to characterize splashback and other ones present in the literature is related to the definition of spherical density profiles. Mansfield et al. (2017) and Deason et al. (2021) found that the most successful method to achieve a clear splashback feature for individual haloes is to measure the median profile along multiple angular directions. In light of the results of Fig. 6, we argue that the distribution of splashback as a function of direction is skewed by the presence of a few dense filaments and hence the difference between a median and mean splashback can be substantial. Therefore, we stress that future works should exercise caution when employing such methods. The use of median profiles smooths substructure by focusing on the halo boundary in the proximity of voids, but because this process is itself correlated with the halo growth rates (see Fig. 5), the connection with observations is not as simple as one might expect.

Next steps
Because in this work we have focused only on high-mass objects (M ∼ 10 14.5 M ), a natural future step is to investigate if the results apply also in other regimes. For example, a larger sample over a wide range of masses and redshifts is required to confirm the simple form of equation (5) and verify if it applies to lower mass groups (M ∼ 10 13.5 M ).
Exploring a wider range in mass, both in observations and simulations, can also be used to confirm a key prediction: because the median accretion rate is expected to be a function of mass and redshift (More et al. 2015), we expect the mass-size relation for an observed halo sample to not necessarily follow a simple form.
Finally, we point out that our results encourage a concentrated effort towards understanding the relationship between cluster environment and splashback. What is discussed in Figs 5, 6, and 9 suggests that the connection between accretion-flows, filaments, and cluster boundary is not a simple one. To better understand this process, it will become necessary to complement the usual insideout theoretical approaches to splashback, that look at haloes to define their boundaries, with outside-in approaches, that connect the cosmic web to its nodes. In this context, the amount of splashback data gathered by projects such as the Kilo Degree Survey (de Jong et al. 2013), Dark Energy Survey (The Dark Energy Survey Collaboration 2005), and, in the future, LSST (LSST Science Collaboration 2009) and Euclid (Laureijs et al. 2011) will provide a powerful probe for the study of structure formation.

AC K N OW L E D G E M E N T S
We thank Benedikt Diemer for providing valuable feedback on the manuscript. OC is supported by a de Sitter Fellowship of the Netherlands Organization for Scientific Research (NWO). HH acknowledges support from the VICI grant 639.043.512 from NWO. YMB acknowledges funding from the EU Horizon 2020 research and innovation programme under Marie Skłodowska-Curie grant agreement 747645 (ClusterGal) and the NWO through VENI grant 639.041.751. The Hydrangea simulations were in part performed on the German federal maximum performance computer "HazelHen" at the maximum performance computing centre Stuttgart (HLRS), under project GCS-HYDA/ID 44067 financed through the largescale project "Hydrangea" of the Gauss Center for Supercomputing. Further simulations were performed at the Max Planck Computing and Data Facility in Garching, Germany.

DATA AVA I L A B I L I T Y
The data underlying this article will be shared on reasonable request to the corresponding authors.