# Can the splashback radius be an observable boundary of galaxy clusters?

Théo Lebeau<fnref n="1" />,*, Stefano Ettori<fnref n="2" />,<fnref n="3" />,<fnref n="1" /> ID, Nabila Aghanim<fnref n="1" />, and Jenny G. Sorce<fnref n="4" />,<fnref n="1" />,<fnref n="5" /> ID

Received 27 March 2024 / Accepted 24 May 2024

## ABSTRACT

The splashback radius was proposed as a physically motivated boundary of clusters as it sets the limit between the infalling and the orbitally dominated regions. However, galaxy clusters are complex objects connected to filaments of the cosmic web from which they accrete matter that disturbs them and modifies their morphology. In this context, estimating the splashback radius and the cluster boundary becomes challenging. In this work, we use a constrained hydrodynamical simulation replicating the Virgo cluster embedded in its large-scale structure to investigate the impact of its local environment on the splashback radius estimate. We identify the splashback radius from 3D radial profiles of dark matter density, gas density, and pressure in three regions representative of different dynamical states: accretion from spherical collapse, filaments, and matter outflow. We also identify the splashback radius from 2D-projected radial profiles of observation-like quantities: mass surface density, emission measure, and Compton-y. We show that the splashback radius mainly depends on the dynamics in each region and the physical processes traced by the different probes. We find multiple values for the splashback radius ranging from 3.3 ± 0.2 to 5.5 ± 0.3 Mpc. In particular, in the regions of collapsing and outflowing materials, the splashback radii estimated from gas density and pressure radial profiles overestimate that of the dark matter density profiles, which is considered the reference value given that the splashback radius was originally defined from dark matter simulations in pioneering works. Consequently, caution is required when using the splashback radius as a boundary of clusters, particularly in the case of highly disturbed clusters like Virgo. We conclude with a discussion of the detection of the splashback radius from pressure radial profiles, which could be more related to an accretion shock, and its detection from stacked radial profiles.

**Key words.** methods: numerical – galaxies: clusters: intracluster medium – galaxies: clusters: individual: Virgo

## 1. Introduction

Clusters of galaxies are tracers of the cosmic matter density and distribution as they are the most massive gravitationally bound structures in the Universe. They are thus used to put constraints on cosmological parameters via their number count (e.g. Planck Collaboration XX 2014; Salvati et al. 2018; Aymerich et al. 2024) or the baryon fraction (e.g. Wicker et al. 2023), which requires a precise and unbiased mass estimation and calibration (e.g. Salvati et al. 2019; Gianfagna et al. 2021; Lebeau et al. 2024; Aymerich et al. 2024 and references therein). This, in turn, necessitates the definition of a cluster boundary in which the mass is estimated. Assuming that galaxy clusters formed via the spherical collapse of initial overdensities (Gunn & Gott 1972), we can estimate their mass in a sphere of a characteristic radius. Considering the self-similar evolution of clusters (Kaiser 1986), characteristic radii are also useful for comparing clusters through normalised quantities.

In the past decade, a new physically motivated radius setting the transition from free-fall to orbital motion of dark matter (DM) particles in N-body simulations has been proposed (Diemer & Kravtsov 2014; Adhikari et al. 2014; More et al.

2015): the splashback radius, $R_{sp}$. It has been studied in various N-body and hydrodynamical cosmological simulations using DM, baryons, and galaxy densities (e.g. Diemer & Kravtsov 2014; Diemer et al. 2017; Mansfield et al. 2017; O'Neil et al. 2021; Aung et al. 2021; Shin & Diemer 2023; Zhang et al. 2023; Towler et al. 2024) or intracluster light (Deason et al. 2021). Diemer & Kravtsov (2014) and Pizzardo et al. (2024) also used galaxy radial velocities. In particular, Pizzardo et al. (2024) proposed to define $R_{sp}$ as the transition point from a concave to a convex shape – that is, the second derivative equal to zero – of the radial profile of galaxy radial velocities.

Diemer (2020) shows that using the splashback radius to define the halo mass provides mass functions that are more universal – in other words, more independent of redshift and cosmology as expected by theoretical models (e.g. Press & Schechter 1974) – than mass functions based on the virial mass or spherical overdensities' masses. According to them, this confirms that $R_{sp}$ is a physically motivated definition of a DM halo boundary. However, $R_{sp}$, and thus $M_{sp}$, the mass enclosed in a sphere of this radius, was determined by following DM particles' trajectories in Diemer (2020) (see Mansfield et al. 2017 for details). In contrast, in observations, we can only try to derive $R_{sp}$ from radial profiles of either the galaxy number density or radial velocities, weak lensing, pressure, and intracluster light. Given that galaxy clusters are connected to filaments of the

<footnote n="1">¹ Université Paris-Saclay, CNRS, Institut d'Astrophysique Spatiale, 91405 Orsay, France</footnote>

<footnote n="2">² INAF, Osservatorio di Astrofisica e Scienza dello Spazio, via Piero Gobetti 93/3, 40129 Bologna, Italy</footnote>

<footnote n="3">³ INFN, Sezione di Bologna, viale Berti Pichat 6/2, 40127 Bologna, Italy</footnote>

<footnote n="4">⁴ Université de Lille, CNRS, Centrale Lille, UMR 9189 CRISTAL, 59000 Lille, France</footnote>

<footnote n="5">⁵ Leibniz-Institut für Astrophysik (AIP), An der Sternwarte 16, 14482 Potsdam, Germany</footnote>

<footnote>* Corresponding author; theo.lebeau@universite-paris-saclay.fr</footnote>
cosmic web (e.g. Cautun et al. 2014; Peebles 2020; Gouin et al. 2021, 2022, 2023; Galárraga-Espinosa et al. 2020, 2021, 2022), defining their boundary in this complex environment is complicated.

Moreover, $R_{sp}$ is located in clusters’ outskirts, which can be very challenging to observe due to the signal’s faintness. Still, numerous tentative detections in various wavelengths were performed by stacking the observed signal of a sample of clusters to reach a sufficiently high signal-to-noise ratio. Then, $R_{sp}$ was identified as the minimum gradient of the radial profile (see more details in Sect. 2.2) in every observable. In the optical, galaxy number density (e.g. More et al. 2016; Chang et al. 2018; Shin et al. 2019; Adhikari et al. 2021; Bianconi et al. 2021; Baxter et al. 2021; Rana et al. 2023), weak lensing (Contigiani et al. 2019; Shin et al. 2021; Fong et al. 2022; Giocoli et al. 2024), and intracluster light (Gonzalez et al. 2021) were used. However, in our case study we cannot stack over a cluster sample because we only consider one. Moreover, due to the small number of galaxies in a single cluster, the detection of $R_{sp}$ from the galaxy number density or radial velocity profile is hardly reliable. Finally, at sub-millimeter wavelengths, through the thermal Sunyaev-Zel’dovich (tSZ, Sunyaev & Zeldovich 1972) effect, a pressure decrease in the intracluster medium (ICM) was detected and associated with the splashback radius (e.g. Anbajagane et al. 2022, 2024; Rana et al. 2023).

In the present study, we investigate the sensitivity of $R_{sp}$ to the dynamics in cluster outskirts and its observability. To this end, we conduct a detailed study of a state-of-the-art constrained hydrodynamical simulation reproducing the Virgo cluster in its local environment (Sorce et al. 2021). This Virgo replica, like the real observed one (see Sorce et al. 2021; Lebeau et al. 2024 and references therein), is a highly disturbed cluster connected to multiple filaments; it is thus a particularly interesting test case. We first present the methodology in Sect. 2, including selecting regions of accreting, outflowing, and filament material to investigate the sensitivity of the splashback radius estimate to specific dynamics. We identify $R_{sp}$ either on DM or gas density 3D radial profiles as well as on pressure 3D radial profiles in Sect. 3. Since we use a constrained replica of the Virgo cluster in its local environment (see Sorce et al. 2019 for the constraining method), we can produce projected maps along the line of sight from the Milky Way to Virgo in our simulation, which is consistent with Virgo’s observations. We thus identify $R_{sp}$ on 2D-projected radial profiles extracted from maps of mass surface density (SD), X-rays’ emission measure (EM), and the Compton-y signal in Sect. 4. In Sect. 5, we summarise our results, draw perspectives about the identification of $R_{sp}$ in observations of Virgo or other clusters, discuss the use of pressure, electron density, and entropy as tracers of $R_{sp}$, and compare the impact of using median or mean profiles to identify $R_{sp}$. Finally, we present our conclusions in Sect. 6.

## 2. Methodology

In this section, we first introduce the simulated replica of Virgo, define characteristic radii of galaxy clusters, present the method adopted to compute 3D and 2D-projected radial profiles, and finally identify the regions our study focuses on.

### 2.1. The Virgo replica simulation

Recently, Sorce et al. (2019) produced replicas of Virgo within its cosmic environment by using constrained initial conditions

of the local Universe (see Sorce et al. 2016 for details). This set of DM simulations is in remarkable agreement with observations. The most representative simulated halo was then re-used to produce a high-resolution hydrodynamical simulation of Virgo (Sorce et al. 2021), which also agrees well with observations. More details about the simulation and its agreement with observations can be found in Sorce et al. (2021) and Lebeau et al. (2024).

Like all the CLONE (Constrained LOcal and Nesting Environment) simulations (e.g. Ocvirk et al. 2020; Libeskind et al. 2020; Dolag et al. 2023) using initial conditions of the local Universe from Sorce et al. (2016) and Sorce (2018), the Virgo replica zoom-in simulation relies on the Planck Collaboration XVI (2014) cosmological parameters; namely, the total matter density, $\Omega_m = 0.307$, dark energy density, $\Omega_\Lambda = 0.693$, baryonic density, $\Omega_b = 0.048$, amplitude of the matter power spectrum at $8\ \text{Mpc h}^{-1}$, $\sigma_8 = 0.829$, Hubble constant, $H_0 = 67.77\ \text{km s}^{-1}\ \text{Mpc}^{-1}$, and spectral index, $n_s = 0.961$. This particular CLONE was produced using the RAMSES (Teyssier 2002) adaptive mesh refinement (AMR) code. The zoom region is contained in a $500\ \text{Mpc h}^{-1}$ local Universe box. It is a $30\ \text{Mpc}$ diameter sphere with a resolution of $8192^3$ effective DM particles of mass $m_{\text{DM}} = 3 \times 10^7\ \text{M}_\odot$. The AMR grid has a finest cell size of $0.35\ \text{kpc}$. The simulation also contains sub-grid models for star formation, radiative gas cooling and heating, and kinetic feedback from the active galactic nucleus (AGN) and type II supernova (SN) similar to the Horizon-AGN implementation of Dubois et al. (2014, 2016). Furthermore, the AGN feedback model has been enhanced by orientating the jet according to the black hole (BH) spin (see Dubois et al. 2021 for details).

The properties of DM and gas cells were extracted using the rdramses code<fnref n="1" /> first used in Renaud et al. (2013). The Tweed et al. (2009) halo finder has been used to identify the Virgo DM halo and its galaxies using the stars. In this work, we define the gas as the baryonic component in the simulation cells, in which hydrodynamics equations are solved following an Eulérian approach, and thus it does not include the stars. In the first approach, we did not apply any density or temperature selection; the gas thus includes its cold component from galaxies. All the gas cells and DM particles were used to compute the radial density profiles. On the contrary, to compute 3D pressure radial profile and produce Compton-y maps, we selected cells of ionised gas – namely, cells with a temperature above $10^5\ \text{K}$, which is the ionisation threshold – and removed the cells associated with the galaxies (see Lebeau et al. 2024 for more details). To produce projected EM maps from the ICM gas, we only selected the cells with a temperature above $10^7\ \text{K}$, given that this is approximately the minimum observable temperature of X-rays telescopes.

### 2.2. Characteristic radii of galaxy clusters

Assuming a spherical collapse model (Gunn & Gott 1972), we estimate the mass of a galaxy cluster within a sphere of a given radius defining its boundary. A first physically motivated boundary can be estimated using the virial theorem (Gunn & Gott 1972; Lacey & Cole 1993; Peebles 2020). A cluster at dynamic equilibrium considered as virialised has a potential energy ($E_p$) equal to the opposite of twice its kinetic energy ($E_k$), $E_p + 2E_k = 0$. The distance at which this equation is valid is the virial radius ($R_{vir}$) and defines a boundary inside which the galaxy cluster matter is at dynamical equilibrium. Following this theorem, we can estimate cluster masses with galaxy velocities in optical

<footnote n="1">¹ https://github.com/florentrenaud/rdramses.</footnote>
<table><caption>**Table 1.** Virgo's characteristic radii (top) and splashback radii (bottom).</caption><thead><tr><th colspan="4">Characteristic radius (Mpc)</th></tr><tr><th>R<sub>500c</sub> = 1.1 Mpc</th><th>R<sub>200c</sub> = 1.7 Mpc</th><th>R<sub>vir</sub> = 2.1 Mpc</th><th>R<sub>200m</sub> = 2.9 Mpc</th></tr><tr><th colspan="4">Splashback radius (Mpc)</th></tr><tr><th>3D</th><th>DM (R<sub>sp,DM</sub>)</th><th>Gas (R<sub>sp,gas</sub>)</th><th>Pressure (R<sub>sp,press</sub>)</th></tr></thead><tbody><tr><td>Collapsing material</td><td>3.4 ± 0.2</td><td>3.9 ± 0.2</td><td>4.3 ± 0.3</td></tr><tr><td>Outflowing material</td><td>4.3 ± 0.3</td><td>4.9 ± 0.3</td><td>5.5 ± 0.3</td></tr><tr><td>Filament material</td><td>–</td><td>–</td><td>–</td></tr><tr><td>Full material</td><td>–</td><td>–</td><td>–</td></tr><tr><th>2D</th><th>Surface density (R<sub>sp,SD</sub>)</th><th>Emission measure (R<sub>sp,EM</sub>)</th><th>Compton-y (R<sub>sp,y</sub>)</th></tr><tr><td></td><td>3.3 ± 0.2</td><td>3.9 ± 0.2</td><td>4.9 ± 0.3</td></tr></tbody></table>

**Notes.** The splashback radii are identified using 3D radial profiles in the regions or using the full material and 2D-projected quantities. The value gives the centre of the bin of minimum gradient; the error bars are plus and minus one bin width at this radius.

galaxy surveys (e.g. de Vaucouleurs 1960 for the Virgo cluster) or in numerical simulations (see Vogelsberger et al. 2020 for a review). For this Virgo replica simulation, we found $R_{vir}$ = 2.1 Mpc using the Tweed et al. (2009) halo finder.

However, $R_{vir}$ is not directly accessible when observing the ICM through the thermal Sunyaev-Zel'dovich effect or in the X-rays, since these observables trace the clusters' thermodynamic properties. Consequently, we usually estimate cluster masses using scaling relations (Pratt et al. 2019, and references therein) or the hydrostatic equilibrium (e.g. Ettori et al. 2013) within radii enclosing a given overdensity. More precisely, we define a radius, $R_{Δc}$, at which the density is $\Delta$ times the critical density of the Universe:

$$R_{\Delta c} = R(\rho = \Delta\rho_c) \quad \text{with} \quad \rho_c = 3H^2/8\pi G. \qquad (1)$$

The mass enclosed in a sphere of this radius is thus

$$M_{\Delta c} = \Delta\rho_c \frac{4}{3} \pi R_{\Delta c}^3. \qquad (2)$$

The most widely used characteristic radii are $R_{500c}$ and $R_{200c}$. For our Virgo replica, $R_{500c}$ = 1.1 Mpc and $R_{200c}$ = 1.7 Mpc.

Instead of using the critical density of the Universe, we can also use the mean matter density, $\rho_m = \rho_c \Omega_m$, and define $R_{Δm}$. Usually, we use $R_{200m}$ because it has been shown (e.g. Kravtsov & Borgani 2012) that the ratio of the density at $R_{vir}$, $\rho_{vir}$, over $\rho_m$ is ~200 (~178 in an Einstein-de Sitter Universe, see Bryan & Norman 1998 for other cosmologies). To compare with other works, we normalise every radial profile by $R_{200m}$. For our Virgo replica, $R_{200m}$ = 2.9 Mpc. The values of these characteristic radii are summarised in Table 1 for comparison with the $R_{sp}$ identified in this work.

The focus of this study is the splashback radius, $R_{sp}$, introduced by works on N-body simulations (Adhikari et al. 2014; Diemer & Kravtsov 2014; More et al. 2015). It was initially defined as the apocentre of the first orbit of DM particles free-falling onto a massive DM halo. It is a dynamic definition of a galaxy cluster since it sets the boundary between the region of infalling motion on the cluster and the region in which orbital motion dominates. Infalling particles will pile up at this radius as they inverse their radial velocity, so the particle density will strongly decrease beyond this characteristic radius. It can thus be identified as the steepest slope in the radial density profile. In other words, it is the minimum of the gradient of the radial density profile:

$$R_{sp} = R \left[ \text{Min} \left( \frac{d\log(\rho)}{d\log(r)} \right) \right], \qquad (3)$$

with $\text{Min}$ for the minimum and $\rho$ the density of either DM, gas, or their sum.

## 2.3. 3D and 2D-projected radial profiles

In this work, we identify $R_{sp}$ as the radius of minimum gradient on radial profiles of DM, gas, and their sum, but also pressure, as other works (e.g. Shi 2016; Towler et al. 2024) relate it to $R_{sp}$. We also produce projected maps of observation-like quantities from which we extract 2D-projected radial profiles and once again identify $R_{sp}$ as the radius of minimum gradient, similar to Towler et al. (2024).

We adapted the procedure presented in Lebeau et al. (2024) to compute the radial profiles used in our analysis. We defined a binning of 0.025 on a log scale ranging from 200 kpc to the limit of the zoom-in region (~14.75 Mpc) minus the size of the biggest (~350 kpc) cells out of the zoom-in region. This avoided contamination from low-resolution cells at the borders of the zoom-in region, leading to non-physical profile shape. For 3D density profiles (DM, gas, and their sum), we summed the masses in each spherical shell, and divided the summed mass by the volume of the spherical shell. For 3D gas properties (pressure, temperature, electron density, and entropy), we computed their mass-weighted mean. The error bars are given by the mass-weighted standard deviation, and thus show the dispersion around the mean value.

We also computed observation-like quantities, from which we derived 2D-projected profiles and tried to identify $R_{sp}$. We built maps of SD akin to lensing signal, EM in X-rays, and integrated Compton-y (y). The projected maps were built using the method detailed in Lebeau et al. (2024). They contain 15 728² pixels of 1.4065 kpc length and are 22.122 Mpc wide. We summed the contribution of all the cells along the line of sight in each pixel. In equations (4), (5) and (6), $m_{cell,i}$, $l_{cell,i}$, $P_{e,cell,i}$, $n_{e,cell,i}$, and $n_{p,cell,i}$ are the mass, length, electron pressure, electron density, and proton density of the i-th gas cell and $m_{part,j}$ is the mass of the j-th DM particle.

We built the SD map by summing the DM particles and gas cell masses in each pixel. For the gas cells contributing to more than one pixel, the mass attributed to each pixel is the cell's total mass divided by the number of occupied pixels. We then divided each pixel by its size; the map is in M$_{⊙}$ kpc$^{-2}$. The method is summarised in the following equation giving the SD in a pixel, $\Sigma_{px}$:

$$\Sigma_{px} = \left( \sum_{i}^{N_{\text{baryon cell in px}}} \frac{m_{\text{cell},i}}{S_{\text{cell},i}/S_{\text{px}}} + \sum_{j}^{N_{\text{DM particle in px}}} m_{\text{part},j} \right) / S_{\text{px}}, \qquad (4)$$
<figure><img src="image_1.png" /><figcaption>Fig. 1. Projected maps of SD (left), EM (centre), and Compton-y (y, right). The projection follows the direction between the Milky Way’s assumed position and the Virgo cluster’s centre. The maps are 22.122 Mpc wide and contain 15 728² pixels. The circles represent R<sub>200m</sub> (solid white) and R<sub>sp</sub> (dashed white) identified in each 2D-projected profile extracted from the maps.</figcaption></figure>

where $S_{px}$ is the pixel surface and $S_{cell,i}$ is the surface of the i-th cell.

The pressure in the ICM of galaxy clusters can be studied via the thermal Sunyaev-Zel'dovich effect (see Birkinshaw 1999 for a review) since it is proportional to the electron pressure integrated along the line of sight. To build the Compton-y map, we summed the pressure in cells multiplied by their length in each pixel. The Compton-y intensity in a pixel, $y_{px}$, is

$$y_{px} = \frac{\sigma_T k_B}{m_e c^2} \sum_{i}^{N_{\text{baryon cell in px}}} P_{e,cell,i} l_{cell,i}, \qquad (5)$$

with $\sigma_T$ the Thomson cross section, $k_B$ the Boltzmann constant, $m_e$ the electron mass, and c the speed of light. Compton-y is unitless.

Finally, the density in the ICM is investigated through X-rays observations (Sarazin 1986). In the X-rays soft band, the EM is the proton density times the electron density integrated along the line of sight. To build the map from the simulation box, we summed the product of the proton density and the electron density in the cells with a temperature higher than $10^7$ K multiplied by their length in each pixel. The EM in a pixel, $EM_{px}$, in units of cm$^{-6}$Mpc, is

$$EM_{px} = \sum_{i}^{N_{\text{baryon cell in px}}} n_{p,cell,i} n_{e,cell,i} l_{cell,i}. \qquad (6)$$

In Fig. 1, we present the maps of SD (left), EM (centre), and Compton-y (right). The color bar is on the log scale below each map. The solid white circle is R<sub>200m</sub>; the dashed white circle is R<sub>sp</sub>, identified on the 2D-projected radial profile extracted from the map (see Sect. 4). All the maps show the filamentary structure connected to Virgo: on the SD map we observe the DM halos’ distribution, whereas on the EM and Compton-y maps we see the diffuse gas, and more particularly tracers of its density and pressure, respectively. We also see an almost empty area on the left part of the Compton-y map.

We then extracted 2D radial profiles from the maps by computing the mean (similarly to pressure profiles, we discuss using the median instead in Sect. 5) among the pixels in annuli, using

the same radial binning as for the 3D spherical shells. The error bars are the standard deviation over the pixels in the circular annuli. We chose a projection following the direction between the centre of the large-scale simulation box, almost centred on the Milky Way, and the centre of Virgo; it is thus comparable to a real Virgo observation. For every radial profile presented in this paper, both axes are on a log scale and the radius in the abscissa was normalised by R<sub>200m</sub> = 2.9 Mpc introduced in Sect. 2.2 (shown as a vertical dashed grey line) to compare with other studies.

## 2.4. Selection of studied regions in Virgo

Clusters originally accreted matter via spherical collapse before relaxing and virialising. At present, they can still accrete matter from spherical collapse but also from cosmic filaments connected to them. These accretion regimes can take place at the same time in different regions of a single cluster. Hence, this work takes advantage of the zoom-in Virgo simulation embedded in its local environment to study the impact of the local environment, and so the dynamical state, on the estimation of the splashback radius in different regions. We selected three regions that are representative of distinct accretion regimes: accretion due to spherical collapse, accretion from filament, and matter outflow. They have the same volume and radial range, enabling direct one-to-one comparisons between regions. We defined these regions using both visual inspection and a radial distance-radial velocity phase-space diagram.

First, the 3D visualisation<fnref n="2" /> of the pressure in the zoom-in box used for the visual inspection is presented in Fig. 2. The top panel shows the full box, and the bottom panel shows a slice in the cluster core along the simulation’s z axis. Then, the radial distance-radial velocity phase-space diagram of the ionised gas, defined as gas cells not belonging to galaxies and with $T > 10^5$ K, is presented in Fig. 3. The radial velocity axis is in the range of [−2000, 2000] km s$^{-1}$ divided into 50 bins, and

<footnote n="2">² This clipping and slicing visualisation was made using the PyVista python library (Sullivan & Kaszynski 2019).</footnote>
<figure><img src="image_2.png" /><figcaption>**Fig. 2.** 3D visualisation of the pressure in the Virgo replica simulation. The pressure is in keV.cm⁻³ and is shown on a log scale in the range of [-8, 0]. The top panel shows the full box and the bottom panel shows a slice in the cluster's core along the z axis. The orange sub-box is the region of outflowing material, the blue sub-box is the region of collapsing material, and the pink sub-box is the region of filament material.</figcaption></figure>

the radial distance is in the range of [0, 5] Mpc divided into 100 bins. All the ionised gas in the zoom-in region – Virgo and its neighbourhood – is shown in the top left panel. Throughout this paper, we refer to the full zoom-in region as ‘full material’.

In the pink box in Fig. 2, we see a large filament connected to Virgo that was already highlighted in Lebeau et al. (2024). When observing the phase-space diagram of this region, in the top right panel of Fig. 3, we observe that the majority of the ionised gas, in the range of [1, 5] Mpc, has a high negative radial velocity, indicating that it is falling onto the cluster without being slowed down while entering the ICM. It confirms the presence of a filament shown in other works (e.g. Gouin et al. 2022; Vurm et al.

2023). In the following, we call the matter in this region – that is, gas and DM – ‘filament material’.

Then, in the orange box in Fig. 2, we observe a minimal amount of pressure, and so matter, outside the cluster and an almost spherical cluster shape. In fact, its phase-space diagram, in the bottom left panel of Fig. 3, shows that most of the ionised gas in the range of [1.5, 4] Mpc, has a positive radial velocity, indicating that it is moving away from the cluster. There is no accretion in this region as there is almost no ionised gas with a negative radial velocity outside the cluster; we thus call the matter in this region ‘outflowing material’.

Finally, in the blue box in Fig. 2, Virgo is also quite spherical in this region, but there is more infalling matter outside the cluster compared to the region of outflowing material. Its phase-space diagram, in the bottom right panel of Fig. 3, shows that the ionised gas is falling on the cluster but with a much lower absolute radial velocity than for the filament material, particularly in the range of [1, 2] Mpc. It shows that it has been slowed down while entering the ICM. We consider that we are in the case of a spherical collapse; we thus call the matter in this region ‘collapsing material’.

## 3. Splashback radius from 3D radial profiles

In this section, we identify $R_{sp}$ on 3D radial profiles. We first show the density profiles of DM, gas, and their sum in each region and compare the values found for these quantities. We then identify $R_{sp}$ on the pressure radial profile. The values of the identified splashback radii are summarised in Table 1.

### 3.1. Splashback radius from 3D density profiles

We present the density profiles of DM, gas, and their sum in the regions of collapsing material (Fig. 4, left), outflowing material (Fig. 4, right), filament material (Fig. 5, left), and full material (Fig. 5, right). In each figure, the radial profiles are in the top sub-panel and their gradient is in the bottom sub-panel; the profile displayed in black is the DM, and the profile displayed in light colours (respectively blue, orange, pink, and grey) is the gas. Its dark-coloured counterpart is the sum of DM and gas.

On the one hand, we clearly identify $R_{sp}$ either from DM and gas profiles in the regions of collapsing material (Fig. 4, left) and outflowing material (Fig. 4, right). For the former, we find $R_{sp,DM} = 3.4 \pm 0.2$ Mpc (dashed vertical dark blue line) and $R_{sp,gas} = 3.9 \pm 0.2$ Mpc (dashed vertical light blue line), and for the latter we find $R_{sp,DM} = 4.3 \pm 0.3$ Mpc (dashed vertical dark orange line) and $R_{sp,gas} = 4.9 \pm 0.3$ Mpc (dashed vertical light orange line). Moreover, we observe that the total matter density profiles are very similar to the DM profiles. It is expected given that DM largely dominates the mass budget in clusters; consequently, we find that $R_{sp,tot} = R_{sp,DM}$.

These identified $R_{sp}$ all have different values. First, when comparing $R_{sp,DM}$ and $R_{sp,gas}$ found in the region of collapsing (or outflowing) material, we observe that $R_{sp,DM}$ is about 0.5 Mpc (0.6 Mpc) smaller than $R_{sp,gas}$ – that is, a bin width of 3.5 (2) at this radius, which is quite significant. This can be explained by the fact that DM and gas are not subject to the same physical processes; DM particles are collisionless and will thus not be slowed down while entering the cluster, whereas the gas has viscosity and will be slowed down and heated while encountering the ICM.

Then, when comparing $R_{sp,DM}$ (or $R_{sp,gas}$) between the regions of collapsing and outflowing materials, we observe a 0.9 Mpc (1 Mpc) difference, which is very significant given that
<figure><img src="image_3.png" /></figure>

**Fig. 3.** Radial distance-radial velocity phase-space diagram of the ionised gas using the full material in the zoom-in box (top left), the region of filament material (top right), the region of outflowing material (bottom left) and the region of collapsing material (bottom right). In the region of filament material, most of the gas has a negative radial velocity of around −1000 km s⁻¹, indicating that the gas is not slowed down while entering the ICM, which is typical of a filament. Then, in the region of collapsing material, most of the gas also has a negative radial velocity but this time of around −500 km s⁻¹, meaning that in this region the gas is slowed down while being heated, which is the case for a spherical collapse. Finally, in the region of outflowing material, most of the gas has a positive radial velocity, and we thus have an outflow.

<figure><img src="image_4.png" /></figure>

**Fig. 4.** Comparison of radial density profiles for the region of collapsing material (left) and the region of outflowing material (right). For both sides, the top panel is the density profile and the bottom panel is the gradient. The solid black profile is the DM, the dashed coloured profile (blue for collapsing material and orange for outflowing material) is the gas, and its dotted dark counterpart is the sum of DM and gas. Vertical dashed lines are $R_{200m}$ (grey), $R_{sp}$ from gas (coloured), and DM (dark coloured). In each region, we observe that $R_{sp,DM} < R_{sp,gas}$. We also note that $R_{sp}$ is smaller for both DM and gas in the region of collapsing material compared to the region of outflowing material.

it is three to four (three to five) bin widths at this radius. This reflects the dynamical state in each region; in the region of collapsing material, there is still free-fall matter accretion, so $R_{sp}$ is located at a smaller radius than in the region of outflowing mate-

rial where the accretion rate is much lower. This was expected and is coherent with other studies (e.g. Diemer & Kravtsov 2014; O’Neil et al. 2021 and Towler et al. 2024) given that a high accretion rate induces a steeper potential well, leading to
<figure><img src="image_5.png" /></figure>

**Fig. 5.** Comparison of radial density profiles for the region of filament material (left) and for the full material (right). For both sides, the top panel is the density profile and the bottom panel is the gradient. The solid black profile is the DM, the dashed coloured profile (pink for the region of filament material and grey for the full material) is the gas, and its dotted dark counterpart is the sum of DM and gas. The vertical dashed line is $R_{200m}$ (grey). $R_{sp}$ is not identified clearly in the region of filament material and when using the full material.

a smaller $R_{sp}$. According to Towler et al. (2024), it can also be due to a high kinetic over thermal energy ratio; in our case, it might be a combination of both.

On the other hand, no clear feature is found in the gradient of either DM and gas density profiles in the region of filament material (Fig. 5, left) or with the full material (Fig. 5, right). There is an important minimum on one bin at ~2.4 $R_{200m}$ for the DM in the region of filament material, but this is certainly due to the presence of a massive galaxy or group of galaxies in the filament like all the other narrow peaks in the density profiles inducing a local maximum followed by a local minimum. Two other local minima are located close to $R_{200m}$ both for DM and gas gradients, but they are quite weak compared to the $R_{sp}$ gradient in other regions, as is seen in Fig. 4. Given that, in this region, the matter is funnelled in the cluster by the filament, we are not in the case of a spherical collapse; identifying $R_{sp}$ was thus not expected. The density profiles of the full material are quite similar to their counterparts in the region of filament material. This was expected, given that this region dominates the matter budget. We can observe minima for both the gas and DM gradient in the range of [2.6, 4] Mpc for DM and [2.9, 4.3] Mpc for gas. However, the gradients in this range are not as low as the minimum gradients in the regions of outflowing and collapsing material. The contributions from all the regions are mixed, the identification of $R_{sp}$ is thus challenging. It only yields a mean value that does not encompass the full complexity of the cluster's dynamics.

## 3.2. Splashback radius from 3D pressure profiles

In the left panel of Fig. 6, we present the pressure radial profile (top sub-panel) and its gradient (bottom sub-panel). Similarly to Figs. 4 and 5, the profiles displayed in orange, blue, and pink are the regions of outflowing, collapsing, and filament material, respectively, and the profile displayed in black is the full mate-

rial. For the full material, we observe a peak in the pressure profile at ~3.5 $R_{200m}$ associated with a massive group of galaxies already identified in Lebeau et al. (2024). Apart from this feature, we find the same results as for the DM and gas density profiles (see Figs. 4 and 5): $R_{sp,press}$ is clearly identified in the regions of outflowing and collapsing material with more than one megaparsec difference whereas it is not in the region of filament material and for the full material. $R_{sp,press} = 5.5 \pm 0.3$ Mpc in the region of outflowing material, and $R_{sp,press} = 4.3 \pm 0.3$ Mpc in the region of collapsing material. They are highlighted by dashed vertical blue and orange lines.

These radii are larger than those extracted from the gas density profiles in those regions. However, $R_{sp,press}$ and $R_{sp,gas}$ in a given region are marginally compatible with their uncertainties. The pressure being the product of the electron density and the temperature means that this difference is due to the temperature. In fact, in the top right panel of Fig. 6, we observe an extended temperature increase of up to a factor of 1.6 in the region of collapsing material and a factor of 6.7 in the region of outflowing material compared to the full material, which is compatible with both $R_{sp,gas}$ and $R_{sp,press}$. Multiple reasons could explain this temperature increase. First, we could have heating due to AGN feedback, coherently with Towler et al. (2024), but it might not be powerful enough at these distances from the cluster's core to heat the ICM gas, and thus induce $R_{sp,press}$ to be hundreds of kpc further than $R_{sp,gas}$. $R_{sp}$ was originally defined from density profiles in N-body simulations. Thus, the most probable explanation is that we are not tracing $R_{sp}$ with pressure. We might rather be tracing accretion shocks, which is coherent with the strong temperature increase at these radii. We discuss this in Sect. 5. Moreover, in the region of outflowing material, the shock in this region could be located at a larger radius in part due to the matter inflow.

In addition to the aforementioned shock fronts in the regions of outflowing and collapsing materials, we notice a temperature
<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Left: comparison of radial pressure profiles (top) and their gradient (bottom) in the selected regions and for the full material. Right: comparison of radial temperature (top), electron density (middle), and entropy (bottom) profiles also in the selected regions and for the full material. The colours stand for the same regions as for other figures: orange for the region of outflowing material, blue for the region of collapsing material, pink for the region of filament material, and black for the full material. $R_{sp,press}$ identified with the pressure profile with the outflowing material and the collapsing material are shown by dashed vertical orange and blue lines. $R_{200m}$ is shown in grey. We observe that $R_{sp,gas} < R_{sp,press}$ and that the minimum in the gradient corresponds to temperature and entropy peaks. Therefore, it might be an accretion shock rather than $R_{sp}$.</figcaption></figure>

increase of up to a factor of 1.2 in the region of filament material but at a much smaller radius, identified in Lebeau et al. (2024) at ~850 kpc, which also indicates a shock front. This is due to the matter inflow in the filament penetrating deeply into the cluster without being slowed down and heated when encountering the ICM, as is shown in recent works (e.g. Gouin et al. 2022; Vurm et al. 2023). None of these temperature increases are visible in the full material temperature profile because they are averaged over the overall mean temperature in each bin.

## 4. Splashback radius from observation-like quantities

In this section, we use observation-like quantities projected along the Virgo-Milky-Way line of sight to identify $R_{sp}$. The maps and the method used to build them are presented in Sect. 2. The 2D-projected radial profiles extracted from the maps (see again Sect. 2.2 for details) are presented in Figs. 7, 8, and 9, showing respectively the SD, EM, and Compton-y. In each panel, $R_{200m}$ is represented by a dashed vertical grey line. We manage to identify $R_{sp}$ in each profile; the minimum in the gradient is clearly visible before a flattening of the profile induces an increase in the gradient. Their values can be found in Table 1, and their positions are highlighted on the profiles by a dashed vertical line, in black, orange, and turquoise, respectively.

When using the SD 2D-projected profile shown in Fig. 7, we find that $R_{sp,SD} \sim 3.3 \pm 0.2$ Mpc, in agreement with the value derived from the DM density profile in the region of collapsing material. This value is also in the range of the minimum gradient of the DM density profile of the full material (see right panel of Fig. 5). However, it differs from that derived from the gas density profile of the region of collapsing material. This was expected, given that DM dominates the mass budget in galaxy clusters, and so does the slope of the total mass density profile, as is shown in Figs. 4 and 5. Moreover, $R_{sp,SD}$ is not compatible with that found in the region of outflowing material, either using the DM, $R_{sp,DM} \sim 4.3 \pm 0.3$ Mpc, or the gas, $R_{sp,DM} \sim 4.9 \pm 0.3$ Mpc, density profile. The slope steepening of the radial profile that we associate with $R_{sp}$ extends over at least one order of magnitude

<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Projected radial profiles of SD. The profile is shown in the top sub-panel, and its gradient is shown in the bottom sub-panel. The dashed vertical lines represent $R_{200m}$ (grey) and $R_{sp}$ (black) identified in the 2D-projected profiles extracted from the SD map.</figcaption></figure>

in the DM profiles. Therefore, the steepening at the lowest radius among the regions determines that of the SD radial profile; it explains the agreement with the value found in the region of collapsing material. We can conclude that, in our case and in agreement with other works (e.g. Towler et al. 2024), when using SD map to identify $R_{sp}$, we trace the DM dynamics in the region of collapsing material.

The EM traces the gas density as it is the integral of the electron density times the proton density. We find that $R_{sp,EM} \sim 3.9 \pm 0.2$ Mpc (see Fig. 8), which is, this time, in agreement with the value derived from the gas density profile in the region of collapsing material. This was also expected, given that EM traces the gas density. Once again, the steepening at the lowest radius among the regions, the region of collapsing material in our case, determines that of the EM radial profile. The EM is thus a tracer of the gas dynamics in the region of collapsing material.

Finally, the thermal Sunyaev-Zel'dovich effect traces the pressure distribution; we find that $R_{sp,yc} \sim 4.9 \pm 0.3$ Mpc in
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Projected radial profiles of EM. The profile is shown in the top sub-panel, and its gradient is shown in the bottom sub-panel. The dashed vertical lines represent $R_{200m}$ (grey) and $R_{sp}$ (orange) identified in the 2D-projected profiles extracted from the EM maps.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Projected radial profiles of Compton-y. The profile is shown in the top sub-panel, and its gradient is shown in the bottom sub-panel. The dashed vertical lines represent $R_{200m}$ (grey) and $R_{sp}$ (turquoise) identified in the 2D-projected profiles extracted from the Compton-y maps.</figcaption></figure>

the Compton-y 2D-projected radial profile, shown in Fig. 9, which is much more distant from the centre than its SD and EM counterparts. This is in between the values extracted from the pressure profiles in the regions of outflowing and collapsing material. We would expect that $R_{sp}$ extracted from this 2D-projected profile agrees with that extracted from the pressure profile in the region of collapsing material, similarly to EM and SD. The low-pressure area visible on the left part of the map, inducing a substantial Compton-y intensity decrease, is located in the region of outflowing material and undoubtedly impacts $R_{sp}$ identification. The contribution of the regions of outflowing and collapsing material seems then to be more mixed in this 2D-projected Compton-y profile than in the SD and EM profiles due to this almost-empty region; it is another projection effect. The Compton-y map traces the projected pressure distribution and gives much larger $R_{sp}$ than the SD and EM maps. However, it is still compatible with $R_{sp,gas}$ and $R_{sp,press}$ identified in the region of outflowing material.

<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Comparison of the splashback radii obtained in this work. The values can be found in Table 1. In addition, the minimum and maximum gradients of entropy are shown in pink. In each region we observe that $R_{sp,DM} < R_{sp,gas} < R_{sp,press}$, and that $R_{sp,SD}$ is in good agreement with $R_{sp,DM,Collapsing}$ just like $R_{sp,EM}$ and $R_{sp,gas,Collapsing}$. We note that the minimum entropy gradient is in good agreement with $R_{sp,press}$, where the maximum entropy gradient corresponds well to $R_{sp,gas}$.</figcaption></figure>

# 5. Discussion

In this work, we used DM, gas, and total density as well as pressure radial profiles to identify $R_{sp}$. We also used projected observation-like quantities comparable to the aforementioned 3D profiles. We have found different values ranging from 3.3 Mpc to 5.5 Mpc; we compare them in Fig. 10. It shows that $R_{sp,DM}$, considered as the reference value since $R_{sp}$ was defined from DM simulations, has the smallest value in each region. $R_{sp,gas}$, and $R_{sp,press}$ even more, overestimates $R_{sp}$ compared to the $R_{sp,DM}$ reference value. However, as was discussed above, this is coherent since DM is collisionless whereas gas has viscosity, which leads to a larger gravitational potential well for the gas. In contrast, $R_{sp,press}$ might rather be an accretion shock, as we discuss in more detail below. Moreover, the radius found by each tracer is always at smaller radii in the region of collapsing material than in the region of outflowing material, which is consistent with other works (e.g. Towler et al. 2024) since the former region has a higher accretion rate, leading to a steeper gravitational potential well. Finally, $R_{sp}$ identified from 2D-projected radial profiles of SD, EM, and Compton-y are in quite good agreement with their 3D counterpart (i.e. DM for SD, gas for EM, and pressure for Compton-y) in the region of collapsing material. Figure 10 thus highlights how much the identified $R_{sp}$ depends on the accretion regime in a given region and the probe used to identify it.

## 5.1. Identification of $R_{sp}$ in observations

On the observational side, given that Virgo is located in our direct vicinity, at approximately 16 Mpc (Mei et al. 2007),
comprehensive studies of its galaxy population (e.g. the Next Generation Virgo Cluster Survey Ferrarese et al. 2016 or the Virgo Environmental Survey Tracing Ionised Gas Emission Boselli et al. 2018 surveys) have been carried out. We could thus derive its projected galaxy density or the galaxy radial velocities, but similarly to the present work, the limited number of detected galaxies in this mid-mass cluster, $M_{vir} \sim 6 \times 10^{14} M_{\odot}$ (Lebeau et al. 2024 and references therein), could prevent the identification of $R_{sp}$ from galaxies.

Observing Virgo in the optical wavelength would be very challenging, then, given that it is very extended ($R_{vir}$ angular size of 4°; see e.g. Planck Collaboration Int. XL 2016). Deriving its mass SD would therefore necessitate a dedicated observation program and enough luminous background sources to build a faithful lens model. To our knowledge, this has not been done so far; at best, galaxy-galaxy lensing could be performed in Virgo (Ferrarese et al. 2012). For other individual clusters, it might be possible to build a cluster lens model up to $R_{sp}$ provided there are sufficiently luminous and extended background sources.

Observing Virgo's $R_{sp}$, or any other cluster, in the X-rays would be arduous too because the surface brightness is proportional to the squared density, leading to a very faint signal in the outskirts. As a matter of fact, Urban et al. (2011), with XMM-Newton, and Simionescu et al. (2015), with Suzaku, observed Virgo beyond $R_{200c}$ but only in particular arms using mosaics. The recently launched XRISM telescope (XRISM Science Team 2020; Tashiro 2022) will not be able to study Virgo's outskirts due to its small collecting area; it will focus instead on the bright region Virgo where $M_{87}$ is located. Recently, McCall et al. (2024) took advantage of the SRG/eRosita all-sky survey to study the Virgo cluster up to 3 Mpc. It is unprecedented, but still does not reach the radii at which we identified $R_{sp}$ in this work, given that they are all beyond 3 Mpc (see Table 1). However, they conducted an azimuthal study similar to that performed both in observations (e.g. Ettori et al. 1998; Vikhlinin et al. 1999; Gouin et al. 2020) and simulations Gouin et al. (2022), and they probed an extended ICM emission and a cold front. With dedicated programs in the future, we might be able to conduct region studies at large radii like ours with observation data.

Finally, when considering the thermal Sunyaev-Zel'dovich effect, Virgo was only observed by the Planck telescope (Planck Collaboration Int. XL 2016); the SPT and ACT ground telescopes cannot observe Virgo as it is close to the North Pole. Due to its proximity, Virgo has the largest integrated thermal Sunyaev-Zel'dovich flux (Taylor et al. 2003), so it would be the best target for a first detection of $R_{sp}$ in a single cluster. Moreover, the thermal Sunyaev-Zel'dovich effect is the integrated pressure along the line of sight so it is proportional to electron density, whereas the surface brightness in the X-rays is proportional to the squared density, allowing one to study the ICM in the outskirts. However, Planck Collaboration Int. XL (2016) could only study Virgo up to ~1 Mpc, which would not be enough to identify $R_{sp,press}$. Future sub-millimetre experiments such as LiteBIRD (LiteBIRD Collaboration 2023) might be able to observe clusters' outskirts up to $R_{sp}$.

## 5.2. Pressure as a tracer of $R_{sp}$

Similarly to what has been done with DM and gas density profiles, we identify $R_{sp,press}$ from the gas pressure radial profile in the regions of collapsing and outflowing materials, but at much larger radii. However, $R_{sp}$ is defined as a dynamical boundary; it is thus supposed to be identified from radial density profiles as it traces the underlying gravitational potential well-driving

particles' velocity. The gas in the ICM experiences much more complex physical processes than DM, such as heating and turbulence due to AGN and SN feedback or shocks, particularly in the cluster's outskirts, where cosmic matter is accreted and shocked while entering the ICM. Therefore, when identifying the radius of the steepest slope on pressure radial profiles, we might identify an accretion shock rather than $R_{sp}$.

Some works based on simulations (e.g. Shi 2016; Towler et al. 2024) find a good agreement between the accretion shock radius, $R_{sh}$, and $R_{sp}$ whereas Aung et al. (2021) find a ratio of $R_{sh}/R_{sp} \sim 1.89$, where $R_{sp}$ is defined from DM particles' trajectories and $R_{sh}$ is defined as the radius of the minimum in the slope of the entropy profile. Moreover, Anbajagane et al. (2022) identify two minima in the stacked pressure profile of clusters observed with the South Pole Telescope (SPT; Plagge et al. 2010). They associate the first one with a non-thermal equilibrium between ions and electrons occurring around shocks, as kinetic energy is better converted to thermal energy by ions than electrons, given that they are more massive.

However, this minimum cannot be observed and confirmed in our simulation since we assume local thermal equilibrium between ions and electrons. Anbajagane et al. (2022) identify the second minima as being related to an accretion shock front, for which they define a lower limit of $R_{sh}/R_{sp} \sim 2.16 \pm 0.59$. For simplicity, in the present work, we noted the minimum in the pressure gradient, $R_{sp,press}$, but in fact the ratio of $R_{sp,press}/R_{sp,DM}$ (similar to $R_{sh}/R_{sp}$ in Aung et al. 2021) is roughly 1.28 in the region of outflowing material and 1.26 in the region of collapsing material. It is a bit higher than one but much lower than the ratios found by Aung et al. (2021) and Anbajagane et al. (2022). Consequently, we are not able to state if $R_{sh}$ is distinguishable from $R_{sp}$ using the $R_{sp,press}/R_{sp,DM}$ ratio in our case study, although the temperature increase at these radii tends to motivate the presence of an accretion shock.

## 5.3. Identification of $R_{sp}$ from electron density and entropy profiles

Other thermodynamic quantities could be used as tracers of $R_{sp}$; for instance, it might in principle be possible to identify it from electron density profiles. In our case, we could not identify $R_{sp}$ in any of the electron density profiles (see middle right panel of Fig. 6). Still, for the regions of outflowing and collapsing material, these radii seem to agree with their pressure counterparts since the steepening before the flattening of the profile occurs in the same range. Finally, the entropy, defined as $K = T/n_e^{2/3}$, is a tracer of the thermalisation state of the gas after being heated via shocks or adiabatic compression (see Tozzi et al. 2000 for detection of entropy and Markevitch & Vikhlinin 2007 for a review on shocks). Both Aung et al. (2021) and Towler et al. (2024) searched for the minimum gradient of the entropy profile to identify the accretion shock radius, as was discussed earlier. Towler et al. (2024) also find the maximum entropy gradient at approximately the splashback radius when defined from the gas density profile. This was more or less expected since when there is a minimum in the gas density gradient, there should be a corresponding maximum in the entropy gradient. In our work, in the regions of outflowing and collapsing material, the minimum gradient is at the same radius as $R_{sp,press}$ and the maximum is at approximately the same radius as $R_{sp,gas}$ (see the maximum and minimum entropy gradient lines displayed in pink in Fig. 10). This would tend to show that $R_{sp,press}$ is an accretion shock, even if, as was already discussed, we cannot confirm it
from the $R_{sp,press}/R_{sp,DM}$ ratio. Anyway, the entropy seems to be a tracer of both $R_{sp}$ and $R_{sh}$.

## 5.4. Comparison between the median and mean profiles

As a baseline for our analysis, we chose the mean pressure and 2D-projected radial profiles. Computing the median profiles is a better tracer of an average quantity for log-normal distributions, which can be the case for the ICM properties in observations (see e.g. Zhuravleva et al. 2013; Eckert et al. 2015). In our case, given that we selected the ionised gas cells (see Sect. 2), that we have a very high number of cells per bin, and that the pressure is a quite smoothly distributed quantity, using the median instead of the mean gives quite similar results for 3D pressure profiles (see Fig. A.1 in Appendix A) and for Compton-y 2D-projected profiles (see Fig. B.3 in Appendix B). However, the SD intensity in pixels is not Gaussian-distributed. Indeed, there are important density contrasts due to sub-structures. Therefore, the median 2D-projected profile (see Fig. B.1 in Appendix B) is quite different from its mean counterpart. In particular, there is an important dip within $R_{200m}$ causing a very important gradient. Though we doubt it can be associated with $R_{sp}$ at such a low radius, it emphasises the transition from the halo-dominated region to the background-dominated one. For EM, the median profile (see Fig. B.2 in Appendix B) is consistent with the mean profile, including its gradient, up to $R_{sp,EM}$ but then gets steeper, so its minimum gradient is located a few bins further away. In this case, we are not able to conclude if this is due to the temperature cut, a projection effect, or if the median really identifies $R_{sp,EM}$ at a larger radius, moreover because the EM pixel distribution is not log-normal. In conclusion, even more caution is required when identifying $R_{sp}$ from EM and SD since the method used to compute the radial profile has a significant impact.

## 5.5. Complementarity of this work with statistical studies

We have shown that $R_{sp}$, more precisely the minimum in the gradient of a radial profile, seems to be more a tracer of the dynamical state of a cluster component (DM or gas) in a given region than an observable natural cluster boundary. However, our work suffers one major limitation, since it is a case study of a particularly disturbed cluster, while other cited works conducted statistical studies; our conclusions are thus not to be taken as universal. Yet this first case study of $R_{sp}$ in regions of a highly refined zoom-in simulation highlights how complex a single cluster can be. It also shows that when generalising over a cluster sample we mainly get a mean value averaging the signals in all the cluster environments. Specifically, when stacking over a large cluster sample, the resulting signal is averaged out, which smooths the outskirts, preventing one from identifying $R_{sp}$ (e.g. see the stacking of 100 random projections of Virgo in Figs. B.1, B.2 and B.3 in Appendix B). Our case study is thus complementary to statistical studies that have tested the dependence of $R_{sp}$ on other parameters such as cluster mass, accretion rate, or redshift (Adhikari et al. 2014; Diemer & Kravtsov 2014; More et al. 2015; Diemer et al. 2017; Mansfield et al. 2017; O’Neil et al. 2021).

# 6. Conclusion

In this work, we identified the splashback radius, $R_{sp}$ – the minimum in the gradient of the radial profile – using the radial density profiles of DM, gas, or their sum as well as pressure in

regions of the simulated Virgo replica. These different regions, which contain outflowing, collapsing, and fragmented material, represent dynamical states occurring simultaneously in our Virgo case study. We then built projected maps of SD, EM, and Compton-y along the Virgo-Milky Way line of sight in our constrained simulation of the local Universe. In addition, we discussed the use of different tracers to investigate $R_{sp}$; in particular, the pressure could be a tracer of an accretion shock rather than $R_{sp}$. We also discussed their potential detectability for a single cluster like Virgo, using median profiles instead of mean profiles to identify $R_{sp}$. We pointed out that our conclusions about this particularly disturbed cluster are not to be taken as universal, although this case study is complementary to statistical studies.

We have found that $R_{sp}$, when its identification is possible, varies by ~1 Mpc between regions for a given quantity, which is about three to four bin widths at this radial distance from Virgo. We have also found that in the same region there is an almost 0.5 Mpc difference between the splashback radius extracted from the gas density profile and that extracted from the DM density profile, so two to three bin widths, which was expected, given that these components of galaxy clusters (i.e. gas or DM) do not follow the same physical processes. Moreover, we notice that the contributions from all the regions are mixed when using the matter in the full zoom-in box to identify $R_{sp}$, which prevents us from identifying it clearly.

When extracting the splashback radius from 2D-projected profiles of SD, EM, and Compton-y, we have also found that each traces the physics processes of their underlying component. Moreover, the identified splashback radii are more in agreement with those found in the region of collapsing material than that of outflowing material. This is certainly because the slope steepening in this region determines that of the projected profile, given that it is located at a smaller radius than in the other regions, except for the Compton-y profile, for which there is a consequent projection effect.

This study of regions in this simulated Virgo cluster shows the complex and various dynamics of clusters’ outskirts and their impact on the splashback radius estimate. Caution is thus required when using the splashback radius as a natural boundary of clusters, especially for a complex cluster, as we have demonstrated for Virgo, and when stacking.

**Acknowledgements.** The authors acknowledge the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) for providing computing time on the GCS Supercomputers SuperMUC at LRZ Munich. This work was supported by the grant agreements ANR-21-CE31-0019/490702358 from the French Agence Nationale de la Recherche/DFG for the LOCALIZATION project. This work has been supported as part of France 2030 program ANR-11-IDEX-0003. SE acknowledges the financial contribution from the contracts Prin-MUR 2022 supported by Next Generation EU (n.20227RNLY3 *The concordance cosmological model: stress-tests with galaxy clusters*), ASI-INAF Athena 2019-27-HH.0, “Attività di Studio per la comunità scientifica di Astrofisica delle Alte Energie e Fisica Astroparticellare” (Accordo Attuativo ASI-INAF n. 2017-14-H.0), from the European Union’s Horizon 2020 Programme under the AHEAD2020 project (grant agreement n. 871158), and the support by the Jean D’Alembert fellowship program. The authors thank the IAS Cosmology team members for the discussions and comments. TL particularly thanks Stefano Gallo and Gaspard Aymerich for helpful discussions. The authors thank Florent Renaud for sharing the *rdramses* RAMSES data reduction code.

# References

Adhikari, S., Dalal, N., & Chamberlain, R. T. 2014, JCAP, 2014, 019

Adhikari, S., Shin, T.-H., Jain, B., et al. 2021, ApJ, 923, 37

Anbajagane, D., Chang, C., Jain, B., et al. 2022, MNRAS, 514, 1645

Anbajagane, D., Chang, C., Baxter, E. J., et al. 2024, MNRAS, 527, 9378
Aung, H., Nagai, D., & Lau, E. T. 2021, MNRAS, 508, 2071

Aymerich, G., Douspis, M., Pratt, G. W., et al. 2024, A&A, submitted,
[arXiv:2402.04006]

Baxter, E. J., Adhikari, S., Vega-Ferrero, J., et al. 2021, MNRAS, 508, 1777

Bianconi, M., Buscicchio, R., Smith, G. P., et al. 2021, ApJ, 911, 136

Birkinshaw, M. 1999, Phys. Rep., 310, 97

Boselli, A., Fossati, M., Ferrarese, L., et al. 2018, A&A, 614, A56

Bryan, G. L., & Norman, M. L. 1998, ApJ, 495, 80

Cautun, M., Van De Weygaert, R., Jones, B. J., & Frenk, C. S. 2014, MNRAS, 441, 2923

Chang, C., Baxter, E., Jain, B., et al. 2018, ApJ, 864, 83

Contigiani, O., Hoekstra, H., & Bahé, Y. M. 2019, MNRAS, 485, 408

Deason, A. J., Oman, K. A., Fattahi, A., et al. 2021, MNRAS, 500, 4181

de Vaucouleurs, G. 1960, ApJ, 131, 585

Diemer, B. 2020, ApJ, 903, 87

Diemer, B., & Kravtsov, A. V. 2014, ApJ, 789, 1

Diemer, B., Mansfield, P., Kravtsov, A. V., & More, S. 2017, ApJ, 843, 140

Dolag, K., Sorce, J. G., Pilipenko, S., et al. 2023, A&A, 677, A169

Dubois, Y., Pichon, C., Welker, C., et al. 2014, MNRAS, 444, 1453

Dubois, Y., Peirani, S., Pichon, C., et al. 2016, MNRAS, 463, 3948

Dubois, Y., Beckmann, R., Bournaud, F., et al. 2021, A&A, 651, A109

Eckert, D., Roncarelli, M., Ettori, S., et al. 2015, MNRAS, 447, 2198

Ettori, S., Fabian, A., & White, D. 1998, MNRAS, 300, 837

Ettori, S., Donnarumma, A., Pointecouteau, E., et al. 2013, Space Sci. Rev., 177, 119

Ferrarese, L., Cote, P., Cuillandre, J.-C., et al. 2012, ApJS, 200, 4

Ferrarese, L., Côté, P., Sánchez-Janssen, R., et al. 2016, ApJ, 824, 10

Fong, M., Han, J., Zhang, J., et al. 2022, MNRAS, 513, 4754

Galárraga-Espinosa, D., Aghanim, N., Langer, M., Gouin, C., & Malavasi, N. 2020, A&A, 641, A173

Galárraga-Espinosa, D., Aghanim, N., Langer, M., & Tanimura, H. 2021, A&A, 649, A117

Galárraga-Espinosa, D., Langer, M., & Aghanim, N. 2022, A&A, 661, A115

Gianfagna, G., De Petris, M., Yepes, G., et al. 2021, MNRAS, 502, 5115

Giocoli, C., Palmucci, L., Lesci, G. F., et al. 2024, A&A, 687, A79

Gonzalez, A. H., George, T., Connor, T., et al. 2021, MNRAS, 507, 963

Gouin, C., Aghanim, N., Bonjean, V., & Douspis, M. 2020, A&A, 635, A195

Gouin, C., Bonnaire, T., & Aghanim, N. 2021, A&A, 651, A56

Gouin, C., Gallo, S., & Aghanim, N. 2022, A&A, 664, A198

Gouin, C., Bonamente, M., Galárraga-Espinosa, D., Walker, S., & Mirakhor, M. 2023, A&A, 680, A94

Gunn, J. E., & Gott, J. R., III 1972, ApJ, 176, 1

Kaiser, N. 1986, MNRAS, 222, 323

Kravtsov, A. V., & Borgani, S. 2012, ARA&A, 50, 353

Lacey, C., & Cole, S. 1993, MNRAS, 262, 627

Lebeau, T., Sorce, J. G., Aghanim, N., Hernández-Martínez, E., & Dolag, K. 2024, A&A, 682, A157

Libeskind, N. I., Carlesi, E., Grand, R. J., et al. 2020, MNRAS, 498, 2968

LiteBIRD Collaboration 2023, Prog. Theor. Exp. Phys., 2023, 042F01

Mansfield, P., Kravtsov, A. V., & Diemer, B. 2017, ApJ, 841, 34

Markevitch, M., & Vikhlinin, A. 2007, Phys. Rep., 443, 1

McCall, H., Reiprich, T. H., Veronica, A., et al. 2024, A&A, in press, https://doi.org/10.1051/0004-6361/202449391

Mei, S., Blakeslee, J. P., Côté, P., et al. 2007, ApJ, 655, 144

More, S., Diemer, B., & Kravtsov, A. V. 2015, ApJ, 810, 36

More, S., Miyatake, H., Takada, M., et al. 2016, ApJ, 825, 39

Ocvirk, P., Aubert, D., Sorce, J. G., et al. 2020, MNRAS, 496, 4087

O'Neil, S., Barnes, D. J., Vogelsberger, M., & Diemer, B. 2021, MNRAS, 504, 4649

Peebles, P. J. E. 2020, The Large-scale Structure of the Universe, 98 (Princeton University Press)

Pizzardo, M., Geller, M. J., Kenyon, S. J., & Damjanov, I. 2024, A&A, 683, A82

Plagge, T., Benson, B., Ade, P. A., et al. 2010, ApJ, 716, 1118

Planck Collaboration XVI. 2014, A&A, 571, A16

Planck Collaboration XX. 2014, A&A, 571, A20

Planck Collaboration Int. XL. 2016, A&A, 596, A101

Pratt, G., Arnaud, M., Biviano, A., et al. 2019, Space Sci. Rev., 215, 1

Press, W. H., & Schechter, P. 1974, ApJ, 187, 425

Rana, D., More, S., Miyatake, H., et al. 2023, MNRAS, 522, 4181

Renaud, F., Bournaud, F., Emsellem, E., et al. 2013, MNRAS, 436, 1836

Salvati, L., Douspis, M., & Aghanim, N. 2018, A&A, 614, A13

Salvati, L., Douspis, M., Ritz, A., Aghanim, N., & Babul, A. 2019, A&A, 626, A27

Sarazin, C. L. 1986, Rev. Mod. Phys., 58, 1

Shi, X. 2016, MNRAS, 461, 1804

Shin, T.-H., & Diemer, B. 2023, MNRAS, 521, 5570

Shin, T.-H., Adhikari, S., Baxter, E., et al. 2019, MNRAS, 487, 2900

Shin, T.-H., Jain, B., Adhikari, S., et al. 2021, MNRAS, 507, 5758

Simionescu, A., Werner, N., Urban, O., et al. 2015, ApJ, 811, L25

Sorce, J. G. 2018, MNRAS, 478, 5199

Sorce, J. G., Gottlöber, S., Yepes, G., et al. 2016, MNRAS, 455, 2078

Sorce, J. G., Blaizot, J., & Dubois, Y. 2019, MNRAS, 486, 3951

Sorce, J. G., Dubois, Y., Blaizot, J., et al. 2021, MNRAS, 504, 2998

Sullivan, B., & Kaszynski, A. 2019, J. Open Source Software, 4, 1450

Sunaev, R., & Zeldovich, Y. B. 1972, Comm. Astrophys. Space Phys., 4, 173

Tashiro, M. S. 2022, Int. J. Mod. Phys. D, 31, 2230001

Taylor, J. E., Moodley, K., & Diego, J. 2003, MNRAS, 345, 1127

Teyssier, R. 2002, A&A, 385, 337

Towler, I., Kay, S. T., Schaye, J., et al. 2024, MNRAS, 529, 2017

Tozzi, P., Scharf, C., & Norman, C. 2000, ApJ, 542, 106

Tweed, D., Devriendt, J., Blaizot, J., Colombi, S., & Slyz, A. 2009, A&A, 506, 647

Urban, O., Werner, N., Simionescu, A., Allen, S., & Böhringer, H. 2011, MNRAS, 414, 2101

Vikhlinin, A., Forman, W., & Jones, C. 1999, ApJ, 525, 47

Vogelsberger, M., Marinacci, F., Torrey, P., & Puchwein, E. 2020, Nat. Rev. Phys., 2, 42

Vurm, I., Nevalainen, J., Hong, S., et al. 2023, A&A, 673, A62

Wicker, R., Douspis, M., Salvati, L., & Aghanim, N. 2023, A&A, 674, A48

XRISM Science Team 2020, ArXiv e-prints, [arXiv:2003.04962]

Zhang, Y., Adhikari, S., Costanzi, M., et al. 2023, Open J. Astrophys., 6, 46

Zhuravleva, I., Churazov, E., Kravtsov, A., et al. 2013, MNRAS, 428, 3274
**Appendix A: Comparison of median and mean pressure radial profiles in the regions of collapsing and outflowing materials**

<figure><img src="image_12.png" /></figure>

**Fig. A.1.** Pressure radial profiles in the regions of collapsing (solid blue for the mean and dashed dark blue for the median) and outflowing materials (solid orange for the mean and dashed dark orange for the median). The vertical dashed lines are $R_{sp, Collapsing Material}$ (blue), $R_{sp, Outflowing Material}$ (orange), and $R_{200m}$ (grey).

Fig. A.1 compares the mean (light profiles) and median (dark profiles) pressure radial profiles in the regions of collapsing (solid blue and dashed dark blue) and outflowing materials (solid orange and dashed dark orange). It shows that using the median instead of the mean to identify $R_{sp}$ gives the same result in the region of outflowing material. In the region of collapsing material, $R_{sp}$ is identified at a smaller radius in the median profile but still in agreement with that identified in the mean profile given the uncertainties. However, in the median profiles a steepening of the profile occurs at a smaller radius, although the median and mean profiles are less than $2\sigma$ from each other at worst.

**Appendix B: 2D-projected radial profiles of 100 stacked random projections compared to the mean and median profile of one projection**

We present the SD (solid light red profile in Fig. B.1), EM (solid dark red profile in Fig. B.2), and Compton-$y$ (solid dark blue profile in Fig. B.3) 2D-projected radial profiles of 100 stacked random projections. The simulation box was rotated using randomly chosen angles to create 100 projections, from which we computed the maps and derived the 2D-projected mean profiles following the method presented in Sect. 2. We then calculated the average profile over the hundred realisations. For comparison, we also show 2D-projected mean (dotted black, red, and cyan, respectively) and median (dashed grey, orange, and green, respectively) radial profiles. It shows that the stacked profile tends to have a higher signal in the outskirts than the mean or median profiles from the projection (named Cen in the figures) used in this work. This was expected because the filaments have various orientations in the random projections, leading to a smoothed stacked profile. We thus cannot identify $R_{sp}$ in any of these stacked 2D-projected profiles. We discuss the differences between the mean and median profiles in Sect. 5.

<figure><img src="image_13.png" /></figure>

**Fig. B.1.** SD 2D-projected radial profile of 100 stacked random projections, displayed in solid light red. For comparison, we also show 2D-projected mean (dotted black) and median (dashed grey) radial profiles. The vertical dashed lines are $R_{200m}$ (grey) and $R_{sp}$ (black).

<figure><img src="image_14.png" /></figure>

**Fig. B.2.** EM 2D-projected radial profile of 100 stacked random projections, displayed in solid dark red. For comparison, we also show 2D-projected mean (dotted red) and median (dashed orange) radial profiles. The vertical dashed lines are $R_{200m}$ (grey) and $R_{sp}$ (red).

<figure><img src="image_15.png" /></figure>

**Fig. B.3.** Compton-$y$ 2D-projected radial profiles of 100 stacked random projections, displayed in solid dark blue. For comparison, we also show 2D-projected mean (dotted cyan) and median (dashed light blue) radial profiles. The vertical dashed lines are $R_{200m}$ (grey) and $R_{sp}$ (cyan).