Inferring the dark matter splashback radius from cluster gas and observable profiles in the FLAMINGO simulations

The splashback radius, coinciding with the minimum in the dark matter radial density gradient, is thought to be a universal definition of the edge of a dark matter halo. Observational methods to detect it have traced the dark matter using weak gravitational lensing or galaxy number counts. Recent attempts have also claimed the detection of a similar feature in Sunyaev-Zel'dovich (SZ) observations of the hot intracluster gas. Here, we use the FLAMINGO simulations to investigate whether an extremum gradient in a similar position to the splashback radius is predicted to occur in the cluster gas profiles. We find that the minimum in the gradient of the stacked 3D gas density and pressure profiles, and the maximum in the gradient of the entropy profile, broadly align with the splashback feature though there are significant differences. While the dark matter splashback radius varies with specific mass accretion rate, in agreement with previous work, the radial position of the deepest minimum in the log-slope of the gas density is more sensitive to halo mass. In addition, we show that a similar minimum is also present in projected 2D pseudo-observable profiles: emission measure (X-ray); Compton-$y$ (SZ) and surface mass density (weak lensing). We find that the latter traces the dark matter results reasonably well albeit the minimum occurs at a slightly smaller radius. While results for the gas profiles are largely insensitive to accretion rate and various observable proxies for dynamical state, they do depend on the strength of the feedback processes.


INTRODUCTION
Galaxy clusters are the result of hierarchical structure formation, forming from the collapse of dark matter overdensities.The gravitational potential provided by these massive halos allows gas, galaxies and stars to reside within and form the galaxy clusters we observe.Defining an edge of these systems is not trivial.Currently, spherical overdensities are used, where the boundary radius is defined as the point within which the average density of a cluster reaches a certain value.For example,  200m and  200c correspond to the radii at which the average cluster density within reaches 200 times the mean and critical density of the universe respectively and  200m and  200c are the masses within these radii respectively.It has been proposed to instead define a cluster's boundary following the trajectories of dark matter particles.The splashback radius is a boundary between infalling and collapsed dark matter, defined as the radius of the apoc-★ E-mail: imogen.towler@manchester.ac.uk entre of the first orbit of dark matter particles (Diemer & Kravtsov 2014) and as such is a model-independent physical definition of the cluster edge.It has been shown using simulations that this radius can be identified using the local logarithmic slope of the dark matter density profiles, where the steepest slope corresponds to the dark matter particles piling up at the apocentre of their first orbit (Diemer et al. 2017).Diemer (2020) found that defining the mass function of a halo using the splashback radius/mass leads to a more universal mass function, i.e. more independent of redshift and cosmology, than using radii such as  200m .For ΛCDM models, they find the universality to be similar to using either the virial mass or spherical overdensities.However, in some alternative cosmologies, the splashback mass functions are significantly more universal.Therefore, the splashback radius is likely to be a more meaningful halo boundary definition.
The splashback radius has been detected observationally using a variety of methods.Often, this is done using a stacked sample of galaxy clusters by using a proxy to measure the mass content of the cluster, for example via weak lensing (e.g.Chang et al. 2018;Contigiani et al. 2019;Shin et al. 2021;Fong et al. 2022).In addition, the splashback radius has been detected using measurements of the galaxy number density (e.g.More et al. 2016;Baxter et al. 2017;Zürcher & More 2019;Shin et al. 2021;Kopylova & Kopylov 2022;Rana et al. 2023).However, past works have found that the method of galaxy cluster selection affects the obtained splashback radius (e.g.Busch & White 2017), with optically selected clusters resulting in smaller splashback radii than expected from simulations (More et al. 2016;Chang et al. 2018;Shin et al. 2019).Currently, there has only been one tentative measurement of the splashback radius of an individual cluster rather than a stacked sample, measured using intracluster light (Gonzalez et al. 2021).
Simulations, both -body and hydrodynamical, have been used to investigate the splashback radius.They allow both the measurement of the trajectories of the individual dark matter particles as well as of the slope of the density profiles to determine the splashback radius, which Diemer & Kravtsov (2014) showed are the same feature.Investigation has led to the well known negative correlation between specific mass accretion rate and the splashback radius of a cluster (normalised using the spherical overdensity radius,  200m ) (e.g.Diemer & Kravtsov 2014;Diemer et al. 2017;Mansfield et al. 2017;Deason et al. 2021;O'Neil et al. 2021;Diemer 2022).In addition, Diemer et al. (2017) and O'Neil et al. (2021) (the latter of which used IllustrisTNG, a full hydro simulation) have both shown further correlations between the splashback radius of a cluster and its mass and redshift.Furthermore, simulations allow the investigation of projection effects on the obtained splashback radius.For example, Deason et al. (2021), who studied both dark matter and stellar density profiles in simulated clusters, found that on average the location of the caustic in projected profiles is approximately ten per cent smaller than the caustic found in 3D density profiles.
It may also be possible to infer the position of the dark matter splashback radius from gas observables.Lau et al. (2015) and Aung et al. (2021) studied a set of 65 clusters from the Omega500 nonradiative cosmological simulation Nelson et al. (2014), finding the location of the minimum gas density gradient (and maximum entropy gradient) to occur around the same location as the splashback radius.O'Neil et al. (2021) studied the Illustris-TNG300 simulation (e.g.Nelson et al. 2018) and showed that the gas density gradient minima consistently occurred at smaller radii than for the dark matter profiles (around 20-30 per cent on cluster scales).Observationally, Anbajagane et al. (2023) found a similar feature when analysing SZ Compton- (i.e.projected gas pressure) profiles of around 10 5 clusters.
As the gas is collisional, it does not behave in the same way as the dark matter so any physical origin of the association between their density gradient minima is unclear.One possibility is that the gas feature is associated with the accretion shock.Shi (2016) studied selfsimilar spherical collapse models (Bertschinger 1985) and showed that this radius coincides with the splashback radius in clusters with  = 5/3 and moderate accretion rates.However, Aung et al. (2021), who locate the shock radius using the minimum entropy gradient, find cluster accretion shocks at around twice the splashback radius (a similar result was found by Baxter et al. 2021, based on SZ Compton-𝑦 profiles).An alternative possibility is that the gas profile shape is the result of the underlying (dark matter-dominated) gravitational potential.
In this work we use the FLAMINGO simulations (Schaye et al. 2023;Kugel et al. 2023) to investigate methods of identifying the splashback radius in galaxy clusters, focusing particularly on identifying a reflection of the splashback feature in baryonic profiles and whether it is possible to find a corresponding splashback feature in potentially observable projected profiles.The FLAMINGO simulation suite contains cosmological boxes run with full hydrodynamics up to a box of side length 2.8 Gpc.This results in hundreds of thousands of simulated galaxy clusters with  200m > 10 14 M ⊙ , giving an excellently sized sample to determine the best way to identify the splashback radius and potentially its reflection in cluster gas properties.
In Section 2, we summarise the FLAMINGO simulations and the different baryonic models we have analysed in this work.We also define the profiles we have extracted from the simulations, both 3D density profiles and 2D projected observable profiles.In Section 3, we discuss stacking the cluster profiles.Next, we discuss the minima obtained from the log-slope of various cluster profiles, how these depend on cluster properties and how well the minima in the dark matter and gas densities correspond.Furthermore, in Section 3.3, we discuss how the correspondence between the minima in the dark matter and gas densities varies between simulations with different baryonic physics and cosmological models.In Section 4, we present our results for projected observables, both gas and dark matter, and investigate which give a good estimate for the splashback radius.Finally in Section 5, we summarise our results.

FLAMINGO SIMULATIONS
The FLAMINGO simulations (Full-hydro Large-scale structure simulations with All-sky Mapping for the Interpretation of Next Generation Observations Schaye et al. 2023) are a suite of cosmological simulations for cluster physics and cosmology run using the Smoothed Particle Hydrodynamics (SPH) code SWIFT (Schaller et al. 2023).The suite contains a variety of cosmological simulation boxes with side lengths up to 5.6 Gpc when running dark matter only and 2.8 Gpc with full hydrodynamics.
FLAMINGO bases its subgrid prescriptions on those developed for the OWLS (Schaye et al. 2010) and EAGLE (Schaye et al. 2015) projects.This includes element-by-element radiative cooling and heating rates from Ploeckinger & Schaye (2020); stellar mass loss from stellar winds arising from core-collapse supernovae, type Ia supernovae, massive stars and asymptotic giant branch stars implemented as described in Wiersma et al. (2009) and Schaye et al. (2015); stellar feedback as in Dalla Vecchia & Schaye (2008) and Chaikin et al. (2023), which is implemented by kicking SPH neighbours of young star particles; placing black hole seeds in sufficiently massive regions following Di Matteo et al. (2008) and Booth & Schaye (2009); and thermally implemented AGN feedback following Booth & Schaye (2009).
A variety of models were run in addition to the fiducial model in a 1 Gpc box, this includes 8 alternative astrophysics variations and 4 additional cosmologies.The astrophysical variations were calibrated to different values of the low-redshift galaxy stellar mass function and galaxy cluster gas fractions (See Table 2 of Schaye et al. 2023, for details of the variations of the observable data).Varying four of the subgrid parameters allowed these observed quantities to be altered in the resulting simulation by fixed amounts (the variations of these four parameters are given in Table 1 of Schaye et al. 2023).In Section 3.3, in addition to the fiducial model, we look at results from the alternative astrophysics models that vary the cluster gas fraction (labelled as fgas+2, -2, -4 and -8).The subgrid parameters of the fiducial model were calibrated directly to observations whereas the different fgas models were calibrated to match observed errorbars of the cluster gas fraction using machine learning optimisation (Kugel et al. 2023).In addition, we also look at the two models that alter the AGN feedback mechanism from thermal injection to jet feedback (Jet and Jet_fgas-4) using the method of Huško et al. (2022).These were separately calibrated to observed data or their errorbars.This results in more energy output from AGN feedback being distributed to the outskirts of clusters, albeit non-isotropically.
In addition to the effect of the baryonic model, we also briefly look at the cosmology variations of the simulation (see Appendix B).The fiducial model uses the cosmological model given by the Dark Energy Survey Y3 (DES Collaboration et al. 2022) and the alternative models use results from Planck (Planck Collaboration et al. 2020) or the "lensing cosmology" from Amon et al. (2023) and include varying neutrino masses.
In this work we analyse the results of two simulation box sizes: 1.0 and 2.8 Gpc with a resolution giving a particle gas mass of  gas ≈ 10 9 M ⊙ , these are labelled as L1_m9 and L2p8_m9 respectively.The larger of these was only run with the fiducial model.We use data from the smaller box to compare the different cosmological models and baryonic physics runs as well as to investigate the effects of projection.In both cases, our sample is selected such that all clusters have a mass of  200m > 10 14 M ⊙ .

Profile definitions
In this work, profiles were extracted from galaxy clusters in the FLAMINGO simulation data set.3D density profiles for both gas and dark matter ( gas and  DM respectively) were obtained by centring a series of spherical shells on the VELOCIraptor (Elahi et al. 2019) determined halo centre within 0.1 -5  200m with 44 equally spaced logarithmic bins and measuring the total mass of each particle type in that shell.The density profiles were extracted for all clusters with a mass  200m > 10 14 M ⊙ (defined as the mass within  200m ), giving approximately 16,000 clusters for each model in the 1 Gpc box and 380,000 clusters in the 2.8 Gpc box.In addition to the 3D density profiles for gas and dark matter, we also measure the mass-weighted gas temperature, where we weight the temperature by the mass () of the  th particle.We use the gas density and temperature profiles to estimate the pressure, , and entropy, , profiles as where we assume homogeneous clusters with primordial abundance giving  = 0.59 as the mean molecular weight.
In addition to the 3D profiles, we extract potentially observable 2D profiles to investigate how well we can obtain the splashback radius from different observables.These include the total mass surface density profile (relevant to weak lensing), hot gas emission measure (a proxy for the soft X-ray band) and integrated Compton- profiles (SZ).For each of these, a series of cylindrical shells were placed around the cluster centre of potential with a total depth of 10  200m , this depth was found to be sufficiently deep to ensure the profiles were converged.Each cylindrical bin is split into  seg = 50 angular segments and a median value for each radial bin is calculated (more detail of this process is in Towler et al. 2022) to reduce the noise in the profile due to substructures (Mansfield et al. 2017;Deason et al. 2021).While this azimuthal median method has been used observationally in X-rays, the splashback radius is located in the very outskirts of clusters, where the signal to noise is limited.Therefore, it will only be possible to measure observable profiles using this technique in future surveys.We measure the surface density, Σ, by calculating the total mass density in each segment and taking a median over all segments (represented by ⟨⟩ seg ) in each radial bin, where  seg is the volume of the angular segment.The emission measure is computed as where  H = 0.76 represents the hydrogen mass fraction,  e = 1.14 the mean molecular weight per free electron and,  seg the cross-sectional area of the angular segment.Finally, the integrated Compton- is measured following We measure these profiles three times for each cluster, one for each perpendicular projection along each simulation axis and each of these are treated as an independent cluster profile.

Stacking
The splashback radius is large enough that observers normally need to stack profiles to be able to identify the splashback feature in the outskirts of clusters.In addition, clusters contain intrinsic scatter in their density profile due to the presence of substructure and so we can improve the noise levels of the extracted density profiles by stacking them and obtain an average profile for similar clusters.We split the cluster sample over bins of equal spacing in the chosen quantity for stacking.Once a set of profiles have been stacked, we smooth the profiles using the fourth order Savitzky-Golay smoothing algorithm (Savitzky & Golay 1964) with a window size of the 19 nearest bins to remove any remaining noise, but this effect is minimal due to the large sample used.From these profiles, we take the radial gradient of the log profiles and then extract the radius of the gradient minimum in the profile to get the splashback radius or the location of the gas minima.In this section, we discuss the criteria we use to stack, using criteria which are all indicators of cluster dynamical state, and the effect of the cluster selection and stacking on the profiles.

Theoretical criteria
It has been shown that the splashback radius is strongly correlated with the accretion rate of the cluster (Diemer & Kravtsov 2014).The gravitational potential in high accretion rate clusters deepens faster, decreasing the splashback radius because it reflects earlier infall.The correlation between these means that stacking in bins of accretion rate leads to profiles with similar splashback radii and so the stacked splashback feature is relatively narrow and will not be broadened due to stacking.Following Diemer et al. ( 2017), we define the specific accretion rate to be where  dyn is the dynamical time and corresponds to approximately ( −  dyn ) = 2 3 ( = 0.5) for halos at () = 1 ( = 0) (Diemer 2017;Deason et al. 2021).
The ratio of the kinetic and thermal energy ( kin and  therm respectively) within a cluster can also be used to measure the dynamical state of simulated clusters (e.g.Barnes et al. 2017).This is measured using only one snapshot of the simulation and so is a more instantaneous measure of the dynamical state of the cluster than the accretion rate, which is measured over a dynamical time.In addition, as it is measured using the energy of the gas, it can capture any dynamical differences that arise in the gas that might not exist in the dark matter, e.g.due to feedback processes.It is calculated via summing over the mass, cluster rest frame speed (  ) and temperature of the gas particles  within  200m .Similar to the accretion rate and the true mass of a cluster, this is a purely theoretical quantity and so cannot be derived from observational data.

Observable criteria
Due to the well known correlation between the splashback radius and the accretion rate (e.g.Diemer & Kravtsov 2014;Wetzel & Nagai 2015;Mansfield et al. 2017;O'Neil et al. 2021), we would ideally stack cluster profiles in bins of accretion rate.However, it is not directly observable, so instead we use methods of measuring the current dynamical state of clusters, as those that have recently accreted large amounts of mass, e.g. through mergers, are much more likely to be disturbed.We investigate a series of gas morphology criteria to investigate their correlations with the accretion rate and splashback radius.These were measured using emission measure maps of the clusters and therefore trace the gas distribution within a cluster.We use the following: • The concentration parameter, which compares the total emission measure within two apertures of 0.15 500c and  500c to find clusters with a brighter, cooler core, which tend to be more relaxed (Peterson & Fabian 2006) and have a larger .
• The symmetry statistic (Mantz et al. 2015), where a series of  el = 5 ellipses have been fit to an emission measure map (a proxy for surface brightness in this case) of a cluster at different brightness levels varying between 0.1 − 1.0 500c .The distances between the centres of the ellipses and the cluster centre,  , , are compared with the average of the minor and major axes of the th ellipse, ⟨ el ⟩  .This measures the symmetry of a cluster around its global centre, in this case the point of minimum potential.A higher value of  shows that a cluster is more symmetrical and therefore more relaxed.
• The alignment statistic, which is measured using the same fitted ellipses as the symmetry statistic.However, this parameter aims to measure how the amount of substructure shifts at different radii.Therefore, it instead compares the distances between the centres of adjacent ellipses,  , +1 , to the average of the ellipse axes of the same adjacent ellipses, ⟨ el ⟩.
Similarly to the symmetry statistic, a larger value of  shows a cluster is more relaxed.
• The centroid shift, (Maughan et al. 2012) measures the distance between the centroid of the surface brightness and the global centre of the cluster (Δ), averaged over  = 8 increasingly smaller apertures within 0.15 − 1.0 500 .Smaller values of ⟨⟩ correspond to a smaller shift and the clusters are therefore more relaxed.
In addition to these morphology criteria, we investigate whether the magnitude gap, the difference in magnitude between the brightest cluster galaxy (BCG) and the  th brightest galaxy, can be used as a proxy for the accretion rate.Over time, satellite galaxies get tidally disrupted and stripped of matter, and it has been shown that larger satellites are affected more by dynamical friction.Therefore, as a halo ages, the brightness gap between the BCG and the brightest satellites grows.Shin & Diemer (2022) proposed that the magnitude gap should negatively correlate with accretion rate as a large magnitude gap is an indicator of an old halo and hence low accretion rate.Following Farahi et al. (2020), we measure the magnitude gap between the BCG and fourth brightest galaxy and denote this as 14.We use the galaxy -band luminosities measured within a 50 pkpc 3D aperture around the galaxy centre provided by the Spherical Overdensity and Aperture Processor (SOAP1 ) catalogue to determine the galaxy luminosities and hence measure the magnitude gap of each cluster.
We compare both the theoretical and observational criteria taken from L1_m9 in Fig. 1, the relationships between quantities are shown in the off diagonals and the distribution of each quantity on the diagonal.In addition, in the top-right, we show a correlation matrix showing the Pearson correlation coefficients between the different criteria.The morphology criteria were calculated three times for each cluster, one for each perpendicular axis from emission measure maps.In Fig. 1, only one direction is chosen to keep the sample sizes the same between 2D and 3D criteria.We find that there is a weak correlation between the mass of a cluster and almost all other quantities.Therefore, when splitting clusters into bins of the other quantities, each bin will be dominated by low-mass clusters.Conversely, we find particularly strong correlations between the accretion rate and the energy ratio, symmetry statistic and centroid shift.In cases with higher mass accretion, we expect the gas within a cluster to be more disturbed and therefore the energy ratio increases with the accretion rate.In addition, both the symmetry statistic and centroid shift measure how visibly dynamically disturbed the cluster is and so we expect these to correlate well with the accretion rate.7), gas energy ratio (Equation 8), magnitude gap, concentration (Equation 9), symmetry (Equation 10), alignment (Equation 11) and centroid shift (Equation 12) for the L1_m9 cluster sample.The diagonal gives the distribution for each quantity.The upper right corner shows a correlation matrix with the Pearson coefficients for each of the different stacking criteria.The dashed box shows the area containing the correlations between the different 2D criteria.

RESULTS FROM 3D PROFILES
In this section, we investigate the effects of stacking the 3D profiles obtained from the  = 0 output of FLAMINGO's fiducial hydro run, the 2.8 Gpc box, L2p8_m9.We investigate both the stacked dark matter and gas profiles as well as the radius and depth of the minima of the slope found in each profile.2021), we find that the minimum (most negative) local gradient in the dark matter profile corresponding to the splashback radius depends on the accretion rate of the clusters (see top left panel).This is to be expected as a larger recent accretion rate results in a steeper potential which leads to a smaller splashback radius.We also find that the clusters with the lowest accretion rate have an additional feature in the dark matter profiles at a smaller radius than the splashback radius.Deason et al. (2021) suggest this is a "second caustic" feature (first discussed in Adhikari et al. 2014), corresponding to a build-up of dark matter particles at the apocentre of their second orbit.Clusters with low accretion rates tend to be older and more relaxed, and so the particles will have had enough time to enter their second orbit.

Effects of stacking DM profiles
We also find a weak mass dependence for the splashback radius   7), middle: mass, and right: gas kinetic-thermal energy ratio (Equation 8).
. mass dependence of the splashback radius may be due to a correlation between the mass and accretion rate (Diemer & Kravtsov 2014) as we expect less massive clusters to be more relaxed.However, from Fig. 1, we expect this correlation to be weak.We investigate this directly by splitting our cluster sample into both mass bins and bins of accretion rate within that.Fig. 3 shows the resulting dark matter and gas density gradient profiles when stacked in this manner.We find that the dark matter profiles are nearly independent of the cluster mass for a fixed range of accretion rate values.Therefore, it is likely that the small mass dependence in the dark matter profiles of Fig. 2 originates entirely from the accretion rate dependence.
In Fig. 2, we also investigate whether there is a correlation between the splashback radius and the gas kinetic-thermal energy ratio within clusters.In general, the density gradient of the more relaxed clusters, which have a lower fraction of kinetic energy ( E ), match that of the clusters with smaller accretion rates and have a larger splashback radius.In addition, the strong correlation between the mass accretion rate and energy ratio is clear from the fact that the second caustic feature in the least accreting clusters is also visible in the clusters with the lowest energy ratio.
Fig. 4 explicitly shows the parameter dependence of the minimum gradient radius (top row) and depth (bottom row) for the dark matter density, gas density and gas pressure profiles.The left panels shows how the radii and depths of the minima depend on the accretion rate.In the upper-left panel, we compare our FLAMINGO results with the More et al. (2015) model for dark matter density profiles, with the values for the free parameters, , ,  and  found in    (Diemer & Kravtsov 2014;Deason et al. 2021;O'Neil et al. 2021).In addition, Fig. 4 shows the dependence on mass (central panels) and energy ratio (rightmost panels).Typically, the mass dependence of the radius of the splashback feature in the dark matter is attributed to larger halos having a higher accretion rate (e.g.Diemer & Kravtsov 2014), which we also find to be true (see Fig. 3).

Effects on stacking gas property profiles
Beyond the accretion shock radius, the gas traces the dark matter component.Aung et al. (2021) found the accretion shock to be 20-100 per cent larger than the splashback radius, when the former is defined as minimum in the log-slope of the entropy.Farahi et al. (2022) also found that the dark matter and gas densities are tightly coupled beyond  200c and the correlation between the two quantities weakens within 0.3  200c of the centre of the cluster.Therefore, there is likely a coupling between the gas and dark matter density at scales close to the splashback radius, even if it is not particularly strong.
While minima in the log-slope of the dark matter density profile are dependent on the orbital dynamics within a halo, the cluster gas is strongly affected by shocks.Shi et al. (2016) finds that for an adiabatic index of  ≈ 5/3, the self-similar collapse model predicts that the splashback radius and the accretion shock radius align.However, other works investigating the accretion shock radius have found it to be much larger than the range in which we expect to find the splashback radius.As mentioned before, Aung et al. (2021) find that the shock radius is 1.89 times larger than the splashback radius.In addition, Anbajagane et al. (2022) measured the location of minima in stacked observed Compton- (projected thermal electron pressure) profiles.They find two minima, one at a large radius of 4.58  200m , which they find is consistent with accretion shocks seen in other works, and one at 1.08  200m , consistent with what Anbajagane et al. (2023) found in the SPT and ACT data.They attribute this latter depression to arise from the thermal non-equilibrium between electrons and ions in the intracluster medium.However, they also find that this feature does not appear when creating a comparable stacked sample from The Three Hundred simulations.When including only the most relaxed cluster sample from The Three Hundred simulations, they do find a reproduction of the same minimum even though the simulations do not model non-equilibrium effects between electrons and ions, implying that, in simulations, the existence of this minimum in the extracted Compton- profile depends on the dynamical state of the cluster sample.Both shocks and the splashback lead to a drop in the gas density profile (O'Neil et al. 2021), but the shocks lead to a wider and shallower minimum when profiles from multiple clusters are stacked, which we see in the right column of Fig. 2, showing the stacked gas density profiles.The radius of the minimum of the gas density slope often matches that of the dark matter, but this feature could be a reflection of the splashback in the gravitational potential or a result of shocks.
In Fig. 3, where we split the clusters into bins of both cluster mass and accretion rate, we see that the location of the gas minimum is dependent on both properties (particularly mass), whereas the splashback minimum is solely dependent on the accretion rate.This highlights the effect that the feedback and other baryonic processes are having on the gas density profiles (see also Fig. 4).
In addition to the gas density, we investigate the logarithmic slope profiles of the gas pressure (left) and entropy (right) in bottom two rows of Fig. 2. We find that the pressure profiles also have a minimum gradient at approximately  200m , but not a minimum corresponding to the accretion shock radius in the expected range, 2-3  200m .The location of the pressure minimum corresponds well with the minima at smaller radii found by Anbajagane et al. (2022) and Anbajagane et al. (2023) in observed Compton- profiles.However, they hypothesise that their pressure deficit arises from a thermal non-equilibrium between ions and electrons, which our simulations do not model.We find minima in the log-slopes of the entropy profile at a radius that corresponds to the expected location of a shock feature, and maxima at approximately the splashback radius.The shapes of the entropy profiles match that of Aung et al. ( 2021), but we find that our minima are much shallower and very strongly dependent on the mass of the cluster.However, we found that the way in which the entropy profiles are constructed affects the results in the outskirts.We combine our temperature and density profiles to obtain our entropy whereas, Aung et al. ( 2021) calculates the volume-weighted entropy.
Fig. 4 also shows the radii of the minima in the 3D pressure gradient profiles.The parameter dependence of the radius of the minimum of the pressure gradient profiles closely resembles that of the dark matter for the cluster mass and energy ratio.However, when looking at the radii of the minima when stacking according to the cluster mass accretion rate, we find almost no correspondence between the minima in the dark matter density gradient and gas pressure gradient.For example, in the most relaxed clusters (Γ < 1) the splashback radius is around 50 per cent larger than the gas pressure (and density) gradient minimum radius.
In summary, we find that the minimum of the log-slope in the gas density and pressure profiles often aligns with that of the dark matter.However, the radii of the minima in the gas properties are particularly dependent on the mass of the cluster, showing that there are other baryonic processes affecting these profiles.Furthermore, we have identified minima in the log-slope of the entropy profiles, corresponding to potential shock features which do not appear in the gas density or pressure profiles.

Alternative simulation models
To investigate the effect of the baryonic physics model on the splashback radius, we extract gas and dark matter density profiles from clusters from simulation runs with alternative physics, detailed in Section 2. This includes eight astrophysics variations, which are calibrated to vary the resulting galaxy stellar mass function and galaxy cluster gas fractions at low redshift.Furthermore, there are four alternative cosmology runs including varying neutrino masses(see Schaye et al. 2023, for further details), results from these runs are presented in Appendix B.

Dark-matter-only simulations
Before looking at the different baryonic physics models, we first investigate how the inclusion of baryons affects the dark matter density profiles and recovered splashback radii.We compare the dark matter only run (hereafter DMO), L1_m9_DMO, with the hydro L1_m9 run.The sample chosen from each simulation followed the previous mass cut with all halos with  200m > 10 14 M ⊙ included.However, halos tend to be slightly more massive in DMO simulations due to feedback blowing out baryonic matter in the hydrodynamical simulations.Consequently, there are approximately 600 more halos in the DMO sample than the hydro sample, sufficiently small to not have a large effect on the resulting profiles.We highlight the differences in the splashback radii obtained from each set of profiles in Fig. 5, when stacked according to the mass (left) and accretion rate (right).We find that the differences between the dark matter density profiles for the DMO and hydro runs are minimal.On average, the splashback radius taken from the DMO data is 1 per cent larger than for the same cluster bin in the hydro run.Therefore, the addition of baryons to the simulation has a minimal effect on the presence of the splashback feature in the dark matter, in agreement with O'Neil et al. (2021).The depth of the splashback feature is similarly unaffected, in the DMO clusters the gradient is on average 1 per cent lower than the clusters from the hydro simulation.

Alternative baryonic models
We now look at the splashback feature obtained from dark matter and gas profiles from simulations that vary the resulting cluster gas fraction (models fgas+2, -2, -4 and -8) and feedback mechanism (Jet and Jet_fgas-4, the latter also altering the cluster gas fraction).This allows us to probe whether the baryonic physics model used in a simulation affects the splashback feature and how sensitive it is to varying amounts of AGN feedback.
We compare the radius of the splashback minimum feature with different hydro runs relative to the fiducial run, L1_m9, in Fig. 6, with error bars showing an estimate on the amount of noise we expect due to sampling calculated by bootstrapping the clusters.The top row shows the differences between the splashback feature found in the dark matter density profiles and the second row shows the same for minima in the gas.We find, in agreement with O'Neil et al. (2021), that the location of the splashback radius itself (i.e. in the dark matter density) is mostly unaffected by the baryonic physics model used.Even in the most extreme cases shown, the splashback radius only varies by about 5-6 per cent from the fiducial run but the difference is not significantly larger than the sampling uncertainty.
However, we find that the baryonic physics model used has a much larger effect on the radius of the gas minimum density gradient.The stronger AGN feedback runs, e.g.fgas-4 and fgas-8, tend to have minima at larger radii and similarly the weaker AGN run, fgas+2, at smaller radii.The "Jet" runs result in the gas minima occurring at smaller radii.The location of the minimum is strongly dependent on the amount of feedback, i.e. energy given off by the cluster AGN, within a cluster and increased feedback effectively "blows" out the minimum to a larger radius.Lower mass clusters are more susceptible to the effects of feedback so this is most obvious in the middle panel where we can see models with higher levels of feedback have a stronger effect on the radius of the minima of lower mass clusters.This shows that the location of the minimum in the gas is not strictly defined by the location of the splashback (which is set by gravitational physics) and various other hydrodynamical effects within the gas can easily shape and move it away from the splashback radius.

RESULTS FROM PROJECTED PROFILES
Identifying the splashback radius in 3D profiles is useful to check how it is affected by cluster properties, baryonic physics and statistical effects.However, observers will only obtain projected images of clusters and so the resulting measurement of the splashback radius will be different.We discuss the effects of projection on the log-slope density profiles as well as pseudo-observable profiles to investigate the differences between the splashback radii obtained from 3D density profiles and the radii of the minima of observable profiles.

Projected splashback radius
We can calculate the radius where we expect to find splashback features in projected profiles by fitting a model to the stacked 3D dark matter density profile and projecting that model.Following Diemer & Kravtsov (2014), we fit the following model for the dark matter density: where is the Einasto model describing the inner halo, models the transition region and models the outer density profile.Once the ideal free parameters, { s ,  s ,  t , , , ,   ,   }, have been fitted to each dark matter density profile, we project the model using, and identify the minimum from the gradient of the resulting projected profile.We then use this as the expected splashback radius in 2D to compare with observable profiles.Observational works such as More et al. (2016) and Zürcher & More (2019) have shown that the splashback radius found from projected profiles is smaller than the radius found in 3D.In this work, we find that the minima in these projected gradient profiles are located on average at roughly 0.82±0.03times that of the 3D splashback radii.

Observable profiles
As described in Section 2.1, profiles have been extracted from the FLAMINGO halos to broadly represent what can be obtained from observations.Each halo was projected three times, once for each perpendicular axis of the simulation, and azimuthally averaged emission measure (X-ray), Compton- (SZ) and total surface density (weak lensing) profiles were extracted from each projection with a total thickness of 10  200m .In this section, we show how these profiles are affected by different stacking methods as well as if these observable profiles could be used to accurately measure the splashback radius.
We show maps of an example FLAMINGO cluster of mass  200m = 1.97 × 10 14 M ⊙ and accretion rate Γ = 1.99 in each of these three observables in Fig. 7.The maps include circles denoting the location of the 3D splashback radius for this cluster (1.27  200m , in white) calculated from the dark matter density gradient profiles and the location of the minimum,  min , in each observable's gradient profile (shown in black).The 3D splashback radius, as well as  min , in both emission measure and in surface density are approximately located where one would expect for a cluster of this mass (see Figs. 4  and 9).However,  min for the Compton- profile is smaller than we find in the stacked profiles.This is the result of looking at a singular cluster rather than a stacked cluster.Compton- profiles of individual clusters often contain multiple minima and so the radius of the deepest minimum is found at smaller radii than in stacked profiles where it is smoothed by the larger shock radius feature.
The observable gradient profiles stacked according to the accretion rate, mass and energy ratio are shown in Fig. 8.The log-slope of the emission measure profiles show a much deeper minimum than the log-slope of the gas density (Fig. 4), because of the density-squared dependence of the emission measure.The SZ profiles are similar to what we obtained from the 3D pressure profiles (Fig. 2), however we find much broader minima in the SZ due to projection effects.As the mass of the cluster is dominated by dark matter, we expect that the surface density gradient profiles are similar to the gradient profiles for the 3D dark matter density (Xhakaj et al. 2020).We find that this is true, and that the location of the minima moves to smaller radii and is broader than in the 3D dark matter profiles.
We study the radius of the minimum gradient more closely in Fig. 9.We also include the location of the minima we expect from projecting the same sample of 3D dark matter density profiles (black line).Overall, there is little correspondence between the radius of the minima in the three observables.In general, the surface density matches that of the projected model, again as expected as the surface density will be dominated by dark matter.We consistently identify minima in the emission measure and Compton-y profiles but these occur at a larger radius than the splashback radius.The differences between the locations of these minima and the splashback radius may have implications for the underlying physics of the components each of the observables probe.
Furthermore, in Fig. 9, we find that there is a mild anti-correlation between the accretion rate and the radius of the minimum in the total surface density gradient, similar to what we found for the dark matter density in Fig. 4 but with a slightly less steep curve at low accretion rates.We find that the minimum radius is smaller in the surface density than for the projected 3D dark matter density due to the gas contribution.The minimum of the gas density gradient tends to occur at smaller radii than the dark matter density (see Fig. 4) and so when these are combined in the total surface density, the radius of the gradient minimum is reduced with respect to the dark matter alone.
Fig. 9 also shows that there is no strong correlation between  min for the surface density and the mass or the energy ratio.In addition, we find that the radius of the minimum in the emission measure depends on the mass of the cluster.High-mass clusters have minima at smaller radii, at similar locations as the dark matter splashback radii.Fig. 4 showed that the minima found in 3D gas density gradient profiles do not show much of a correlation with either the accretion rate or the energy ratio.Similarly, we find that  min for the emission measure has a weak correlation with both quantities.Emission measure is expected to scale roughly with density squared and so we expect similar trends between the emission measure and 3D gas density.In addition, while there is not much of a trend between the minima in the Compton- gradient and the accretion rate or energy ratio, the radius of the minimum generally increases with mass.This is at odds with the minimum in the pressure profiles (Fig. 2) which decreases slightly for increasing mass, though this effect is not strong over the mass range used.In Fig. 8, we saw that the Compton- gradient profiles of the high mass clusters have a broader minimum than at lower masses.In the pressure, it is common to have additional minima in the gradient profiles of individual higher mass clusters, but as these profiles have been stacked, the extra minima have been smoothed out to create a single broader minimum.In the projected Compton- gradient profiles, this broadening of the minima tends to lead to  min being identified at larger radii but with an increased uncertainty.However, in the stacked 3D pressure profiles, Fig. 2, the minima at the larger radii are not as significant.The radii of minima are reduced by projection and so the outer minima are further out in 3D and have less of an effect on the pressure profiles around the splashback radius.

Stacking profiles using observables
While the splashback radius has a strong dependence on the accretion rate of the host halo, the accretion rate is not an observable property of clusters.Instead, we aim to find an observable proxy that would allow the cluster profiles to be stacked appropriately.In Fig. 10, we show the three projected gradient profiles (top: emission measure,  middle: Compton-, bottom: surface density) from clusters in the mass range 10 14.2 <  200m /M ⊙ < 10 14.4 .These were stacked in bins of the five properties introduced in Section 2.2.2 (left to right: concentration, symmetry, alignment, centroid shift and magnitude gap).When stacking according to the surface brightness concentration, we find that the profiles are all very similar, particularly towards the minimum.The differences appear in the core of the cluster, which is to be expected as the concentration parameter probes differences in the centre rather than the outskirts of the cluster.The profiles separate more when stacked according to the symmetry statistic, so much so that a small second caustic feature appears in the surface density gradient profiles of the most regularly shaped clusters ( > 1.4).Due to the strong correlation between the symmetry statistic and the accretion rate (see Fig. 1), we expect there to be a large number of low accretion rate clusters within the bins with high , increasing the likelihood for the second caustic to appear.However, this feature is not as clearly defined as seen previously for the accretion rate because the symmetry statistic and accretion rate do not correlate perfectly.
Motivated by the stronger correlations between the accretion rate and the symmetry statistic, centroid shift and magnitude gap, we plot  in Fig. 11 how the extracted minimum radii of the three projected profiles varies with each of these statistics.A smaller symmetry and a lower centroid shift both correspond to more dynamically disturbed clusters and we find that, for both statistics, the radius of the minimum extracted from the surface density and emission measure gradient profiles decreases for more dynamically disturbed clusters, matching the more dynamically disturbed clusters with higher accretion rates and smaller splashback radii.Shin & Diemer (2022) proposed that the accretion rate of halos negatively correlates with the magnitude gap between the BCG and the brightest satellite galaxy.Over time, satellites get stripped of matter and tidally disrupted, and larger galaxies sink to the centres of clusters due to dynamical friction.It has been shown that this effect is most prominent for the largest satellites, meaning that over time the brightness gap increases as a halo grows.We measured the magnitude gap following Farahi et al. (2020), see Section 2.2.2.We find a mild, negative correlation between the accretion rate and the magnitude gap (this was found to be slightly dependent on the resolution of the simulation, see Appendix A) and therefore expect a positive correlation between the magnitude gap and the splashback radius.The rightmost column of Fig. 11 compares  min obtained from the observable profiles with the magnitude gap.We find that there is a slight correlation between the two, most prominent in the surface density profiles.However, we find almost no correlation between the magnitude gap and the minima found in either the emission measure or Compton- gradient profiles.

CONCLUSIONS
In this work we have used clusters from the FLAMINGO cosmological hydrodynamical simulations to obtain dark matter density, gas density, gas pressure and gas entropy profiles.From these we have extracted the location of the minimum local gradient to identify the splashback radius from the dark matter or a potentially matching feature in the gas profiles.We investigated the effect of stacking the profiles in a variety of ways and how the location of the radius changes when using projected profiles.Our results can be summarised as follows: • The splashback radius identified using the dark matter density gradient profile has a strong anti-correlation with the accretion rate, in agreement with previous works (see Figs. 2 and 4).
• The splashback radius has a weak negative mass dependence.Previous works have suggested that this is due to the connection between cluster mass and accretion rate (Diemer & Kravtsov 2014).We find a weak correlation (see Fig. 1) and when stacking in bins of both mass and accretion rate (Fig. 3), while the accretion rate has obvious qualitative effects on the location of the splashback radius, there is not a significant mass dependence.
• We also identify a minimum in the cluster gas density gradient profiles, approximately corresponding to the dark matter splashback radius, but the radius of this minimum has a stronger mass dependence (Figs. 2 and 3).In addition, we find a similar feature in the gas pressure gradient, with no further minima indicating potential shock features at higher radii.However, the gas entropy gradient profiles have minima in the region we expect to find a shock feature in the gas, i.e. 2-3  200m (Aung et al. 2021).
• Comparison between hydrodynamical and dark matter only simulations finds minimal difference in the recovered splashback radii (see Fig. 5).
• Altering the astrophysical parameters used in the simulation results in an essentially unchanged splashback radius.However the radius of the gas minima is much more sensitive to the astrophysics models within the simulation, reducing the correspondence between the splashback and gas minimum (Fig. 6).
• In Section 4, we investigate the effects of projection on the recovered splashback radius.Fitting the dark matter density following Diemer & Kravtsov (2014) to 3D profiles out to 5  200m and projecting the fitted profiles results in the minima occurring at 0.8 times that of the minima in the 3D profiles (Section 4.1).
• We find that the minimum of the total mass surface density gradient profile has a similar radius to what we expect from the projected 3D dark matter profiles.However, although similar, it is found at a systematically smaller radius due to the contribution of the gas (see Fig. 9).
• We find that the minima of the observable profiles (emission measure, Compton- and total surface density) are dependent on the morphology measures used to stack the profiles: to varying degrees more regular clusters tend to have larger splashback radii (see Fig. 11).In addition, the radius of the surface density minima most closely resembles where we expect the splashback radius to be after projection.This is to be expected, as the surface density is dominated by dark matter.However, the minima of the emission measure and Compton- gradient profiles are located at substantially larger radii and are nearly independent of the dynamical state.
The gas density and pressure gradient profiles of the FLAMINGO clusters demonstrate that there exists a minimum at approximately the splashback radius.There are similarly also minima at similar radii in the gradient profiles of gas observables such as the emission measure and Compton-.However, due to the differing physics governing the gas and dark matter motions, it is unclear whether these minima are a true reflection of the splashback radius or coincidentally located at a similar radius.Observational works such as Anbajagane et al. (2023) have started to identify minima at radii similar to what we have found in Compton- gradient profiles of clusters.Future optical/IR, X-ray and SZ observations will be useful to compare with simulations and further our understanding of the behaviour of baryons and dark matter in cluster outskirts.

APPENDIX A: EFFECTS OF RESOLUTION ON MAGNITUDE GAP
The magnitude gap between the BCG and the fourth brightest galaxy within a cluster correlates with the cluster's accretion rate and therefore, its splashback radius.However, galaxies within the cosmological simulation are poorly resolved.In the case of smaller halos or halos with much more mass concentrated in the main halo, the fourth brightest galaxy can be difficult to resolve and therefore our estimates for the magnitude difference between the two galaxies may be more uncertain.The FLAMINGO suite contains multiple resolutions of the same simulation box, L1_m9 for the standard resolution used in this work (gas particle mass  gas ≈ 10 9 M ⊙ and L1_m8 a higher resolution where the particles are a factor of eight less massive (note the models are calibrated separately but to the same observable).This allows us to directly test whether the resolution of the simulation has an effect on the retrieved magnitude gap and its correlation with the accretion rate or splashback radius.
We find that the Pearson correlation coefficient between the accretion rate and the magnitude gap increases in the higher resolution simulation, from -0.35 in L1_m9 to -0.48 in L1_m8.We plot the splashback radius obtained from the dark matter density gradient profile stacked according to the magnitude gap in Fig. A1.We find that there is a slightly steeper relationship between the magnitude gap and the splashback radius in the higher resolution simulation, but the two results are qualitatively similar.

APPENDIX B: ALTERNATIVE COSMOLOGICAL MODELS
In addition to the varying baryonic physics models in Section 3.3, we also investigate the effects of varying the cosmological model on the splashback radius.Fig. A2 shows the variation of the minima in the dark matter (top panels) and gas (bottom panels) density gradient profiles relative to the fiducial model when using different cosmological models.We find that the splashback radius extracted from the dark matter density profiles is essentially unaffected by the change in cosmological model of the simulation, agreeing with the cosmological independence found in Diemer et al. (2017).In addition, we find that the location of the gas minimum is less affected by the change in cosmological model as expected as the cosmological model has less of an effect on the gas physics.Additionally, for both the gas and the dark matter, the effect of the variation of the cosmological model on the depth of the minima is only of the order of a few percent.Thus the depth of the minimum is also much more sensitive to the baryonic physics than to the cosmological model of the simulation.

Figure 1 .
Figure 1.Corner plot comparing the total cluster mass, specific accretion rate (Equation7), gas energy ratio (Equation8), magnitude gap, concentration (Equation9), symmetry (Equation10), alignment (Equation11) and centroid shift (Equation12) for the L1_m9 cluster sample.The diagonal gives the distribution for each quantity.The upper right corner shows a correlation matrix with the Pearson coefficients for each of the different stacking criteria.The dashed box shows the area containing the correlations between the different 2D criteria.

Fig. 2
Fig. 2 shows the stacked dark matter density (top row), gas density (second row), pressure (third row) and entropy (bottom row) gradient profiles in bins of accretion rate (left column), mass (middle column) and energy ratio (right column).In agreement with Diemer & Kravtsov (2014); Diemer et al. (2017); O'Neil et al. (2021), we find that the minimum (most negative) local gradient in the dark matter (see middle column of top row) over the range  200m > 10 14 M ⊙ (where the maximum mass is set by the most massive cluster in our sample,  200m = 10 15.58 M ⊙ ).This is in agreement with what has been found by Diemer et al. (2017); O'Neil et al. (2021).However, the

Figure 2 .
Figure2.Comparison of the effect of using different stacking criteria on the stacked (from top to bottom) dark matter density, gas density, pressure and entropy gradient profiles.Left column: accretion rate, (Equation7), middle: mass, and right: gas kinetic-thermal energy ratio (Equation8).

Figure 3 .
Figure 3. Stacked dark matter (left) and gas (right) density gradient profiles in bins of both accretion rate (Equation 7; different rows) and mass (different colours).

Figure 4 .
Figure 4.The dependence of the splashback radius on the three properties used to bin the halos, accretion rate (Equation 7, left), mass (centre) and ratio of kinetic and thermal energy of the gas (Equation 8, right).The error bars show the three-sigma error found from bootstrapping the cluster sample in each bin used for stacking the profiles.

Figure 5 .
Figure 5.Comparison of the splashback radius extracted from the L1_m9 and L1_m9_DMO simulations when stacking the dark matter density profiles into bins of mass (left) and accretion rate (right).The error bars show the uncertainty by propagating the bootstrap error of the splashback radius from both the DMO and hydro simulations.

Figure 6 .
Figure 6.The variation between the splashback radius obtained for different stacking bins (from left to right, halos are binned based on their accretion rate, mass and energy ratio respectively) for different baryonic physics runs in comparison to the fiducial run, L1_m9.The top row shows the splashback radius obtained from the dark matter density gradient profiles and the bottom row shows the radius of the feature reflected in the gas.The error bars show the uncertainty obtained by bootstrap resampling the clusters from the fiducial model.

Figure 7 .
Figure 7. Maps of an example FLAMINGO cluster with mass  200m = 1.97 × 10 14 M ⊙ and accretion rate Γ = 1.99 at  = 0. From left to right, these maps show the emission measure, Compton- and total mass surface density of the cluster.Each map shows the location of the splashback radius ( SP = 1.27 200m) obtained from the cluster's 3D dark matter density gradient profile in white and the location of the minimum in the gradient profile of the respective 2D observable in black.We find the location of the minima to be  min = 1.04, 0.82 and 0.94 200m for the emission measure, Compton-y and total mass surface density respectively.

Figure 8 .
Figure 8. Stacked projected gradient profiles for the three observables tested in this work.Top: X-ray emission measure, middle: Compton-, bottom: total mass surface density.These profiles have been stacked according to the mass accretion rate (Equation 7, left column), mass (central column) and gas energy ratio (Equation 8, right column) of the clusters.

Figure 9 .
Figure 9.Comparison of the location of the minima of the three different observable gradient profiles used in this work (emission measure, gold; Compton-, blue; and surface density, purple) as well as the expected 2D splashback feature obtained from projecting the 3D dark matter density profiles.Each panel compares different bins used to stack the profiles (left: accretion rate, centre: mass and right:energy ratio).Error bars show the 1  uncertainty from bootstrap resampling.

Figure 10 .
Figure10.Stacked projected gradient profiles of the three observables tested in this work, top: X-ray emission measure, middle: Compton-, bottom: total mass surface density obtained from weak lensing.These profiles have been stacked according to the four morphology criteria and the magnitude gap, see Section 2.2.2.From left to right, the plot shows bins in concentration, symmetry, alignment, centroid shift and magnitude gap.The clusters shown are restricted to the range 10 14.2 <  200m < 10 14.4 M ⊙ .

Figure 11 .
Figure 11.Same as Fig. 9 but instead compares the parameter dependence of the radius of the minimum in the gradient profile with the centroid shift (Equation12, left), symmetry statistic (Equation10, middle) and magnitude gap (right).

Figure A1 .
Figure A1.Comparison between the splashback obtained from dark matter density profiles when stacked according to the magnitude gap for two different resolution simulations.The clusters all have a mass of  200m > 10 14 M ⊙ .The error bars show the expected 1 sigma errors due to resampling.

Figure A2 .
Figure A2.The variation between the splashback radius obtained for different stacking bins (from left to right, halos are binned based on their accretion rate, mass and energy ratio) for runs with different cosmological models in comparison to the fiducial run, L1_m9.The top row shows the splashback radius obtained from the dark matter density gradient profiles and the bottom row shows the radius of the feature in the gas.