© The Authors 2025

1. Introduction

Density profiles of halos are fundamental to cosmological research. Their systematic analysis, first through analytical methods (Gunn et al. 1972; Fillmore & Goldreich 1984; Bertschinger 1985; Hoffman & Shaham 1985) and later through N-body simulations (Davis et al. 1985; Crone et al. 1994; Cole & Lacey 1996; Tormen et al. 1997) and fitting functions (Einasto 1965; Hernquist 1990; Navarro et al. 1996; Diemer 2023), has been crucial in advancing our understanding of the internal dynamics of galaxy clusters and their connection to the cosmic web. A consistent finding across these studies is the emergence of universal characteristics, despite the observable variability among individual structures.

This consistency is accompanied by variations that are significant for understanding dynamic transitions within clusters, notably virialization and splashback. Virialization is theorized to take place at a radius where a structure reaches dynamic equilibrium, no longer undergoing significant mergers (Cuesta et al. 2008). Although this radius is challenging to define in realistic systems, owing to its proximity to galaxy cluster cores, it is generally associated with the zone where the intracluster medium is in hydrostatic equilibrium. The splashback radius, marking the point where infalling particles create a caustic at the first apocenter of their orbits, has recently gained attention in both simulations (Diemer & Kravtsov 2014; Adhikari et al. 2014; More et al. 2015; O’Neil et al. 2021; García et al. 2023) and observational studies (More et al. 2016; Baxter et al. 2017; Chang et al. 2018; Zürcher & More 2019; Anbajagane et al. 2024) as an important measure for delineating cluster boundaries and assessing their accretion rate.

The turnaround radius is another critical scale in the study of halo formation, marking the transition from the collapsing structure to the Hubble flow. Within the spherical collapse model (SCM) framework (Gunn et al. 1972), turnaround represents the outermost boundary of halos occurring at a fixed density with respect to the background with its value being sensitive to the dark energy content of the Universe (Busha et al. 2003; Pavlidou & Tomaras 2014).

Building on these insights, in our previous work (Korkidis et al. 2020, 2023) we showed that the density at the turnaround radius in realistic simulated clusters closely matches the predictions of the SCM. In Pavlidou et al. (2020), we showed that a measurement of the present-day average matter density and its evolution with time at the turnaround scale for a few thousand clusters would be sensitive to the overall matter density in the universe and to the value of the cosmological constant. Such measurements have the potential to provide constraints on cosmological parameters that complement those obtained from Baryon Acoustic Oscillations (BAO), Type Ia Supernovae (SNIa), and the Cosmic Microwave Background (CMB) (Eisenstein et al. 2005; Amanullah et al. 2010; Komatsu et al. 2011).

Despite its potential to introduce new constraints on the cosmological parameters, the study of turnaround in clusters remains under-explored. This is primarily due to the assumption that its observation requires access to the phase space information of galaxies surrounding clusters – a task complicated by redshift-space distortions and projection effects. Consequently, there is a prevailing notion that the turnaround scale is presently unobservable, which has, in turn, discouraged systematic searches for turnaround signals in the density profiles of cluster outskirts.

In contrast to the complexities involved in acquiring phase space information, density profiles are notably easier to observe in actual systems. On that front, the past decade has seen substantial progress in our ability to analyze the outer densities of galaxy clusters. Observations such as the profile of satellite galaxies (Rines et al. 2013; Tully 2015), Sunyaev-Zel’dovich effect and weak lensing (Mandelbaum et al. 2006; Umetsu & Diemer 2017; Shin et al. 2019, 2021; Anbajagane et al. 2024) have now become commonplace for examining the outer regions of galaxy clusters. Inspired by these advancements, this work seeks to re-examine density profiles with a particular focus on the turnaround scale.

Our understanding holds that the turnaround radius should reside beyond the splashback feature in density profile, provided there is a considerable infalling component surrounding the halos under consideration. At these scales, analytical models of density profiles, particularly those assuming spherical symmetry (Prada et al. 2006; Tavio et al. 2008; Betancort-Rijo et al. 2006; Cuesta et al. 2008), have been very successful in describing stacked cluster halos in simulations and have demonstrated their universality up to the onset of significant shell-crossing. In particular, Prada et al. (2006) showed that beyond approximately three Virial radii (which closely corresponds to the turnaround radius at low redshifts in simulations), the effect of shell crossing diminishes significantly.

Motivated by these findings, in our recent work (Korkidis & Pavlidou 2024) we achieved similar success in analytically modeling the outer density profiles of simulated clusters. In this study, we used excursion set theory for the statistical representation of initial profile distributions of halos, and the spherical collapse model (excluding shell-crossing) was employed to follow their non-linear evolution. This approach effectively modeled the outer regions of average cluster density profiles as a function of their collapsed mass. Through this model, we successfully derived scaling relations connecting the collapsed mass of clusters with their turnaround mass, achieving accuracy better than 10% across various cosmologies.

These findings suggest that the SCM without shell-crossing may offer insights into identifying signatures of the turnaround radius within stacked cluster density profiles. Therefore, in this work, we aim to model the outer profiles of clusters and compare them with simulations to understand how deviations correlate with the turnaround radius. To accomplish this, we extend the methodology introduced in Korkidis & Pavlidou (2024) to develop an analytical description of the density profile as a function of its radial distance.

This paper is organized as follows: Section 2 briefly outlines the simulation and halo sample used in this study and details the methodology for constructing density profiles. Section 3 describes how the radial density profile is derived from the mass profile of Korkidis & Pavlidou (2024). In Sect. 4, we compare the profiles and their logarithmic slopes of our analytical model against those of the simulated profiles, demonstrating that the Hubble flow leaves a recognizable mark on the profiles. Finally, section Sect. 5 summarizes our findings.

2. Simulated data

2.1. Simulation overview and halo selection

To analyze the stacked density profiles of cluster sized halos we used one of the dark-matter-only boxes from the second run of the MulitiDark Planck simulations (MDPL2). On the scales of interest in this work baryonic effects have a negligible impact, thus a dark-matter-only box was sufficient for our needs. The simulated box that we employed used a Planck ΛCDM cosmology Planck Collaboration XIII (2016) with a box of 1000 cMpc/h on a side and particle mass of 1.51 × 10⁹ h⁻¹ M_⊙. The halo catalog was produced with the Rockstar algorithm (Behroozi et al. 2013).

For the halo selection we considered the same two samples as in Korkidis et al. (2023): one sample with the 3000 most massive halos in the range 10¹⁴ M_⊙ ≤ M₂₀₀ ≤ 6 × 10¹⁵ M_⊙ and another sample for which we randomly selected 780 halos from logarithmically spaced mass bins, with M₂₀₀ ≥ M_min. We opted for a random halo selection for our second sample to keep the computational costs of the analysis tractable. The value of M_min at each redshift snapshot scaled proportionately to the mass of the largest halo at that time, and was equal to 6 × 10¹³ M_⊙ at a z = 0. The resulting median M₂₀₀ were 7.4, 3.8 and 1.8 ×10¹⁴ M_⊙ for z = 0, 0.5 and 1 respectively. Our choice to analyze the largest halos at each redshift, instead of making a selection based on the halos’ merger history, or accretion rate was mainly related to the fact that these parameters are not expected to affect the profiles at the scales of interest in this work.

2.2. Building the profiles

To generate the profiles, we divided the area around each galaxy cluster’s center into 500 concentric shells, extending to a radius of 25 Mpc, ensuring that each halo was contained within the simulation’s boundaries. The size of our bins was large enough to contain at least 6 × R₂₀₀ in radius for both group and cluster-sized halos. Within the shells, we calculated the dark matter density and average radial velocities of dark matter particles. Following the methodology proposed by Korkidis et al. (2020), we identified the kinetically driven turnaround scale as the radius R_ta of the outermost shell that is not expanding, measured from the cluster’s center. Then we made sure to discard halos that lay within the turnaround radius of a bigger structure. Finally, we normalized the radius of each individual halo’s density profile to its turnaround radius R_ta.

3. Analytic profile

In Section 2 of Korkidis & Pavlidou (2024) we employed excursion set theory for a statistical description of dark matter halos density and environment, and complemented it with the SCM to follow the non-linear growth of these structures. This process allowed us to derive a simple expression for the mode density profile of cluster sized halos as a function of their collapsed mass m. Specifically, we established that the average matter density of a sphere enclosing mass βm (so β > 1 always) around a cluster-sized-halo, is universal, and equal to

$$\begin{array}{c}\hfill {\rho }_{\mathrm{avg}}(\beta )={\rho }_m{(1-{\beta }^{-\gamma })}^{-{\stackrel{\sim }{\delta }}_c},\end{array}$$(1)

where ρ_m is the average matter-density of the Universe at the desired redshift; −γ = dlnS/dlnm is the logarithmic slope of the density field variance S(m) when smoothed on a scale that encloses mass m; and ${\stackrel{\sim }{\delta }}_{\mathrm{c}}$ is the overdensity of a collapsing structure linearly extrapolated to its time of collapse. For detailed recipes to calculate γ and ${\stackrel{\sim }{\delta }}_{\mathrm{c}}$ refer to Appendices A and B of Korkidis & Pavlidou (2024).

From Eq. (1), it is straightforward to derive the local density within a narrow spherical shell of radius¹ $\stackrel{\sim }{R}\equiv R/R_{\mathrm{ta}}$. A detailed derivation can be found in Appendix A. The functional form of this radial density profile is given parametrically by the following two equations:

$$\begin{array}{cc}\hfill \rho (\beta )& ={\rho }_m{[1-{\beta }^{-\gamma }]}^{-{\stackrel{\sim }{\delta }}_c+1},\hfill \end{array}$$(2)

$$\begin{array}{cc}\hfill \stackrel{\sim }{R^3}(\beta )& =(1+{\delta }_{\mathrm{ta}}){[1-{(1+{\delta }_{\mathrm{ta}})}^{-1/{\stackrel{\sim }{\delta }}_c}]}^{1/\gamma }\beta {[1-{\beta }^{-\gamma }]}^{{\stackrel{\sim }{\delta }}_c}.\hfill \end{array}$$(3)

Equation (2) gives the local density on the surface of a sphere enclosing mass βm and Eq. (3) connects the mass enclosed by the sphere with its normalized radius. Since Eq. (3) cannot be solved for β in closed form, it has to be solved numerically or be approximated by a simpler function.

Importantly, equations (2) and (3) describe the local density profile as a function of the turnaround radius and the cosmological parameters on which ${\stackrel{\sim }{\delta }}_c,\gamma \mathrm{and}{\delta }_{\mathrm{ta}}$ depend. From Pavlidou et al. (2020) and Korkidis et al. (2023) we know that the turnaround density and its evolution with cosmic time are sensitive to the value of the cosmological constant. This is true both in the context of spherical collapse and in simulations. In our formalism, the dependence on redshift when the profile is expressed in terms of R_ta enters through δ_ta. This introduces a dependency of the profiles not only on the present-day matter density (Ω_m) but also on the value of the dark energy density Ω_Λ.

4. Comparison between simulated and analytic profiles

It is widely acknowledged that density profiles exhibit universal characteristics when expressed in units of a constant-overdensity radius. Thus, since the kinematically defined turnaround radius corresponds to a constant overdensity as described by the SCM (Korkidis et al. 2020, 2023), it is reasonable to anticipate that our simulated profiles, when normalized with respect to this radius, will exhibit similar behavior.

In Fig. 1 we plot the normalized density as a function of the normalized radius of our halo sample at present time. The density and radii have been normalized with respect to the mean matter density of the universe and the kinematically-defined turnaround radius, respectively. From the density of lines (represented by color in the plot) it is evident that the majority of individual profiles aggregate in a narrow band. Around ∼0.4R_ta there is a steep break in the slope of this band. Further out, the slope increases smoothly (from negative to zero) as the profiles will eventually asymptote to the mean matter density of the Universe. It is this behavior – the gradual change in slope – that we seek to capture with the analytic description of the outskirts of clusters discussed in the previous section.

Fig. 1. Dark matter density (in units of the background-Universe average matter density) plotted against the normalized radius of concentric spherical shells around MDPL2, z = 0 halos. The color scale indicates the number of lines per pixel on the plot. The radii of the shells have been normalised with respect to the kinetically defined turnaround radius.

To make meaningful comparisons between simulations and our analytical predictions, we grouped these profiles into 30 linear radial bins, from 0 to 3 × R_ta, and within each bin we calculated the median value of the density, thus constructing the “median profile”. We then computed the logarithmic slope of this median profile using finite differences and smoothing them with a Savitzky-Golay filter with a 5-point window and a 3rd order polynomial².

In the upper panel of Fig. 2, we overlay our analytical model with the median density profiles at various redshifts of the simulation. In the middle panel we display their logarithmic slope. Since our primary focus is the profile slopes, we opted to compare our analytic profile (which models the mode density) with the median because the different representations mainly affect the profile normalizations. Additionally, calculating the mode would necessitate smoothing the profiles, which can introduce difficult-to-track correlations with the various cluster masses in each radial bin. Finally, the median has the benefit of being more readily calculable and having easily defined errors. For a comparison between different profile statistics, we refer the reader to Appendix B.

Fig. 2. Comparison between simulated and analytic density profiles for different redshifts. Upper panel: solid lines depict the median density profile at different redshifts of our MDPL2 sample in spherical shells of radius R/Rta; shaded regions represent the 1σ spread of the density within each shell; dashed lines represent the analytical prediction for the outer portions of the mode profile discussed in Sect. 3. Middle panel: Logarithmic slope of the density profiles, with dashed lines illustrating the logarithmic slope of the analytical prediction. Lower panel: ratio of logarithmic slope in median simulated profile over the one of the analytically predicted profile. The analytic slope deviates from the one in simulations for R/Rta ≲ 1.

From the figure, three main observations emerge:

(A) Starting at the turnaround radius and for larger radii (R/R_ta ≳ 1), both the profile shapes and their normalization remain largely unchanged, with the slope increasing smoothly as it reaches for the constant background density. From z = 0 to z = 1, we anticipate significant variations in collapsed masses, accretion rates, and concentrations, all of which impact the central cluster regions differently. However, despite these changes in the dynamical state of our halo sample, the outer profiles remains universal and the turnaround is reached close to a slope value of −1.

(B) The shapes of the median profiles at large scales, where the Hubble flow dominates, are well described by the physics encoded in early Universe statistics and spherical collapse without the assumption of shell crossing. This is evident from the excellent capture of profile slopes by our analytical model, as depicted in the middle and lower panels of the figure. In contrast, these assumptions fail profoundly within the infall region (R/R_ta < 1) of the clusters. In these regions, as we approach the cluster cores, the profiles are shaped by a complex interplay of different particle trajectories surrounding the central cluster regions. As expected, this complexity exceeds what the simple assumption of simple spherical collapse can account for. A more involved approach to model this situation would at least require introducing the concept of shell crossing (see for example Prada et al. 2006).

(C) Redshift predominantly influences the inner regions of the density profiles R/R_ta ≤ 0.4. This is evident in the alteration of the profile normalization at small radii, as well as in the position of their splashback feature. The latter is represented by the steepening of the slope profile seen in the middle panel. Despite our simple criterion for sampling the largest structures in the simulated box, the variation in profile normalization can not be attributed to the decrease in the average M₂₀₀ of our halo sample. Rather, this change is mostly driven by the fact that the average R_ta of clusters decreases with redshift and is associated with a higher overdensity with respect to the background Universe.

The divergence between analytical and simulated profiles near the turnaround radius, coupled with the universality of profiles beyond this radius, has a significant implication for the observability of the turnaround radius: the turnaround scale is imprinted in 3-dimensional profiles as the scale beyond which the average/median/mean profile slope³ matches the one expected from early-Universe statistics and spherical collapse. We expect this feature to continue to be true in observable projected (column density) profiles. Projected density profiles observed in the plane of the sky are contaminated at all scales by material unrelated to orbiting parts of clusters. Beyond the turnaround radius, however, even projected material should be well-described by our analytic profile. Consequently, deviations between analytical modeling of the outer profiles and actual projected density profiles, similar to this work, could be closely connected to the location of the turnaround radius.

Another implication is related to the significant impact of the turnaround radius redshift evolution on the variation of profile normalizations. As this evolution is influenced by the values of the cosmological parameters (Pavlidou et al. 2020; Korkidis et al. 2023), we can expect that this dependence will also be reflected in the profiles. In Fig. 3 we investigate this hypothesis by plotting the analytic profile and its logarithmic slope for different values of the cosmological parameters. To construct the profiles we assumed a Hubble constant h₀ = 0.7, a power spectrum slope n_s = 0.96 and median virial mass for the cluster of M₂₀₀ = 10¹⁴ M_⊙. Then, we used the scaling relation of Eq. (11), introduced in Korkidis & Pavlidou (2024) and the known connection between overdensity masses and radii to derive the SCM turnaround radius for clusters with the quoted virial mass (M₂₀₀). The small error associated with this relation makes it a good proxy for our lack of access to the true turnaround radius of the individual clusters that would make these profiles (see Fig. D.1, Appendix D).

Fig. 3. Dependence of the analytical profile and its logarithmic slope on the cosmological parameters. Upper panels illustrate matter density across shells, normalized by the turnaround radius, at redshifts indicated by the titles. With different line colors we show how the profiles change for different values of Ωm (assuming ΩΛ = 0.7) and ΩΛ (assuming Ωm = 0.3). Bottom panels explore the change in the logarithmic slope of the respective density profiles.

Looking at the first two columns of the figure, depicting the model at z = 0, we observe a pronounced dependence on the Universe’s overall matter content. This effect is particularly prominent at low redshifts but slightly decreases as we approach z = 1. While the profile exhibits no discernible change with varying values of the cosmological constant at z = 0, this effect becomes more pronounced with increasing redshift (third and fourth columns of the figure). Hence, the change in the profile sensitivity to Ω_Λ with redshift should be attributed to the turnaround radius. This is because, as the cosmological parameters and redshift change, so does R_ta as it corresponds to a different overdensity as predicted by the SCM. Had we kept the overdensity at a fixed value for all redshifts and cosmologies, we would have mitigated this effect. This is clearly demonstrated in Fig. C.1 of Appendix C, where we illustrate the dependence of the profiles on the cosmological parameters when normalizing the analytic profile with respect to R₂₀₀.

This observation opens a significant avenue to impose constraints on the cosmological parameters through the turnaround density without directly observing the turnaround radius of individual clusters. Just as we used the scaling relation D.1 to find the turnaround radius from R₂₀₀, which were then used to normalize the profiles of Fig. 3, we could similarly apply these scaling relations on observed stacked clusters at different redshifts. Fitting our model to the outer portions of these profiles could, in principle, allow us to derive constraints on the cosmological parameters, while the scale where their slopes start to deviate from universality is expected to correlate strongly with the average turnaround radius of the clusters.

This is, of course, a statement of principle. Applying this technique to data requires a careful assessment of the systematics involved in simulated surveys and mock observations.

5. Conclusions

Despite the extensive research on cluster density profiles, a gap remains in our modeling of their transition to the Hubble flow. To address this gap, we normalized stacked cluster-density profiles from simulations to their turnaround radius which represents the scale at which radial velocities around clusters, on average, become zero. Subsequently, we compared these profiles with an analytical model for the most probable outer radial density profile. This model, a continuation of the work by Korkidis & Pavlidou (2024), is based on Gaussian early Universe statistics and the spherical collapse model without shell crossing and is found to be sensitive both to the overall matter content of the Universe and the cosmological constant (Ω_m and Ω_Λ).

From our analysis, we find the simulated profiles to be universal in shape and well-described by our analytical model beyond turnaround, where the Hubble flow dominates. However, notable deviations are observed within the turnaround radius, signaling the emergence of shell-crossing effects. At smaller scales, after a sudden steepening of the profile slopes at the splashback radius, this universality breaks down. The profiles exhibit variations correlated with changes in their dynamical state, as demonstrated by profile changes at various redshifts within our halo sample. These variations impact both the shapes and normalization of the profiles, challenging the validity of the assumptions made by our model.

Our findings suggest that analytical, spherically symmetric models of halo formation without shell-crossing could be instrumental in identifying the turnaround radius of stacked clusters in observational data. Deviations from such models in stacked density profiles are anticipated to closely correlate with the position of their average turnaround radius. Furthermore, the ability of these models to accurately describe the outer regions of density profiles is expected to be sensitive to the values of cosmological parameters.

Acknowledgments

We acknowledge support by the Hellenic Foundation for Research and Innovation under the “First Call for H.F.R.I. Research Projects to support Faculty members and Researchers and the procurement of high-cost research equipment grant”, Project 1552 CIRCE (GK, VP); by the European Research Council under the European Union’s Horizon 2020 research and innovation programme, grant agreement No. 771282 (GK); and by the Foundation of Research and Technology – Hellas Synergy Grants Program (project MagMASim, VP). The research leading to these results has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie RISE action, Grant Agreement n. 873089 (ASTROSTAT-II). The CosmoSim database used in this paper is a service by the Leibniz-Institute for Astrophysics Potsdam (AIP). The MultiDark database was developed in cooperation with the Spanish MultiDark Consolider Project CSD2009-00064. The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V. (http://www.gauss-centre.eu) and the Partnership for Advanced Supercomputing in Europe (PRACE, https://prace-ri.eu/) for funding the MultiDark simulation project by providing computing time on the GCS Supercomputer SuperMUC at Leibniz Supercomputing Centre (LRZ, http://www.lrz.de). The Bolshoi simulations have been performed within the Bolshoi project of the University of California High-Performance AstroComputing Center (UC-HiPACC) and were run at the NASA Ames Research Center. Throughout this work we relied extensively on the PYTHON packages Numpy (Harris et al. 2020), Scipy (Virtanen et al. 2020) and Matplotlib (Hunter 2007).

References

Appendix A: Derivation of the analytical outer density profile

In this Appendix we use the profile for the average density ρ_avg in a sphere around a collapsed structure of mass m, enclosing mass βm, to derive the local density within a thin spherical shell. A shell of mass mdβ and of volume 4πR²dR (where an increment of radius dR corresponds to an increment in mass mdβ) will have local density given by the ratio:

$$\begin{array}{c}\hfill \rho =md\beta /4\pi R^{2}dR,\end{array}$$(A.1)

whereas the average density within the entire sphere of mass βm will be connected to this mass and radius through

$$\begin{array}{c}\hfill {\rho }_{\mathrm{avg}}=3\beta m/4\pi R^3.\end{array}$$(A.2)

Differentiating this relation with respect to β and substituting in Eq. A.1 we obtain

$$\begin{array}{c}\hfill \rho (\beta )=\frac{m}{4\pi R^2}\frac{3\beta }{R(1-\frac{\beta }{{\rho }_{\mathrm{avg}}}\frac{d{\rho }_{\mathrm{avg}}}{d\beta })}=\frac{{\rho }_{\mathrm{avg}}}{1-\frac{\beta }{{\rho }_{\mathrm{avg}}}\frac{d{\rho }_{\mathrm{avg}}}{d\beta }}.\end{array}$$(A.3)

To calculate (β/ρ_avg)(dρ_avg/dβ), we differentiate Eq. (1) with respect to β to obtain:

$$\begin{array}{c}\hfill \frac{\beta }{{\rho }_{\mathrm{avg}}}\frac{d{\rho }_{\mathrm{avg}}}{d\beta }=\frac{-{\stackrel{\sim }{\delta }}_c\gamma {\beta }^{-\gamma }}{[1-{\beta }^{-\gamma }]}.\end{array}$$(A.4)

Substituting this equation in Eq. (A.3) we obtain

$$\begin{array}{c}\hfill \rho (\beta )=\frac{{\rho }_{m,0}{a}^{-3}{[1-{\beta }^{-\gamma }]}^{-{\stackrel{\sim }{\delta }}_c+1}}{1+({\stackrel{\sim }{\delta }}_c\gamma -1){\beta }^{-\gamma }}.\end{array}$$(A.5)

The product ${\stackrel{\sim }{\delta }}_c\gamma$ for a concordance ΛCDM cosmology and for cluster masses M₂₀₀ ≥ 6 × 10¹³ M_⊙ has a value very close to 1.

Thus, for β > 1, which concerns us in this work, the entire second term in the denominator will be suppressed, and we can write the much simpler expression

$$\begin{array}{c}\hfill \rho (\beta )={\rho }_{m,0}{a}^{-3}{[1-{\beta }^{-\gamma }]}^{-{\stackrel{\sim }{\delta }}_c+1}.\end{array}$$(A.6)

This relation connects the density in a thin shell at the surface of a sphere with its enclosed mass βm. We next turn to expressing this profile in terms of R/R_ta. We normalize Eq. (A.2) to the turnaround scale and use again Eq. (1) for ρ_avg to obtain

$$\begin{array}{c}\hfill {\left(\frac{R}{R_{\mathit{ta}}}\right)}^3=\frac{\beta {[1-{\beta }^{-\gamma }]}^{{\stackrel{\sim }{\delta }}_c}{\rho }_{\mathrm{ta}}}{{\beta }_{\mathrm{ta}}{\rho }_{m,0}{a}^{-3}}=\frac{\beta }{{\beta }_{\mathrm{ta}}}{[1-{\beta }^{-\gamma }]}^{{\stackrel{\sim }{\delta }}_c}(1+{\delta }_{\mathrm{ta}}),\end{array}$$(A.7)

where β_ta is the turnaround mass of a structure in units of its collapsed mass, ρ_ta is the turnaround density predicted by the SCM, and 1 + δ_ta = ρ_ta/ρ_m(a) is the density contrast of a structure at turnaround. For β_ta, Eq. (1) gives

$$\begin{array}{c}\hfill {\beta }_{\mathit{ta}}={[1-{(1+{\delta }_{\mathrm{ta}})}^{-1/{\stackrel{\sim }{\delta }}_c}]}^{-1/\gamma }.\end{array}$$(A.8)

Substituting Eq. A.8 in Eq. A.7 we obtain

$$\begin{array}{c}\hfill {\stackrel{\sim }{R}}^3=(1+{\delta }_{\mathit{ta}}){[1-{(1+{\delta }_{\mathit{ta}})}^{-1/{\stackrel{\sim }{\delta }}_c}]}^{1/\gamma }\beta {[1-{\beta }^{-\gamma }]}^{{\stackrel{\sim }{\delta }}_c}.\end{array}$$(A.9)

Equation A.7 cannot be solved for β in closed form. Thus our options are to either find $\beta =\beta (\stackrel{\sim }{R})$ numerically or approximate it by a polynomial. An example of such an approximation would be to take a Taylor’s expansion around $\stackrel{\sim }{R}=1$ (the turnaround radius) at a desired accuracy; we find that beyond second order of the expansion mostly affect the profile at radii smaller than 0.3R_ta where the profile fails by design. Throughout this work we take the latter approach.

Appendix B: Radial density profile through different statistics

In the upper panel of figure B.1, we plot the mode, the median, and the mean profiles (local density in radial shells, normalized to the mean matter density of the Universe, as a function of shell radius, normalized to the turnaround radius) for our z = 0 sample. To calculate the different the profiles from different statistics we followed the procedure described in section 4 for the median. That is, we grouped the profiles into 30 linear radial bins, from 0 to 3 × R_ta and within each bin we calculated the mode, median and mean value of the density. For the mode we employed a Gaussian kernel density estimator (Virtanen et al. 2020) to model the probability density function (PDF) of densities in the bin. The mode corresponded to the density for which the PDF has a maximum. The shaded region corresponds to the 1σ spread of profiles. In the lower panel we plot the logarithmic slope of each profile. Further out of the splashback feature, different profiles exhibit different normalization (highest for the mean, lowest for the mode). However, the slopes of the profile for radii beyond the turnaround radius appear mostly unaffected by the type of statistic used to form the profile. As expected, the mode profile is more noisy than the median or the mean.

Fig. B.1. Effect of choosing different statistics to represent the density profile. The two panels present the mode, median and mean density profiles along with their logarithmic slopes of our halo sample at z = 0.

Fig. B.2. Analytical profile dependencies on the cosmological parameters at different redshifts. As in Fig. 3, different panels show the change of the density profile and its logarithmic slope with the difference that in this figure the radii have been normalized with respect to R200.

Appendix C: Cosmological dependence of outer profile when normalized to R₂₀₀

Figure 3 illustrates that the most probable outer density profile of Section 3, when normalized with respect to R_ta exhibits a dependence on Ω_Λ that evolves with redshift. However, of the three model parameters ${\stackrel{\sim }{\delta }}_c,\gamma \mathrm{and}{\mathrm{R}}_{\mathrm{ta}}$only the turnaround radius has a redshift dependency, as it corresponds to the overdensity at which the SCM predicts turnaround to occur. It is thus reasonable to expect this sensitivity to Ω_Λ to stem from R_ta.

To test this hypothesis, in C.1 we plot our model and its logarithmic slope as a function of R₂₀₀, for different values of the cosmological parameters and redshift. From the second column of this figure, it is evident that by keeping the radius at a constant overdensity (is this case, 200), we eliminate any dependence the profile would otherwise have on Λ. Simultaneously, we observe that the dependence on Ω_m does not evolve with redshift, in contrast to Fig. 3, where the sensitivity to Λ come at the expense of the sensitivity to Ω_m.

Appendix D: Scaling relation between R_ta and R₂₀₀

Fig. D.1. Correlation between the true turnaround radius Rta,actual and R200 for our cluster sample at z = 0. The red solid line depicts the theoretical scaling relation of Eq. D.1. The blue solid line represents the median value of the blue points in bins of R200. The blue shaded region corresponds to the 16th and 84th percentiles.

In Korkidis & Pavlidou (2024) we demonstrated that our model for the most probable outer density profile of galaxy clusters allows for the derivation of analytical mass scaling relations that depend on cosmology and redshift. These relations connect cluster masses at different overdensities, including the turnaround, and were shown to perform well, with the peak of the error distribution remaining below 15% across all cosmologies, mass samples, and overdensities tested.

Since the kinetically defined turnaround radius corresponds to a constant overdensity as described by the SCM, we can use the known connection between overdensity masses and radii (corresponding to overdensities δ = 200 and δ = δ_ta, respectively) along with the scaling relation of to show

$$\begin{array}{c}\hfill {\mathrm{R}}_{\mathit{ta}}={({\rho }_{200}{\beta }_{\mathit{ta}}/{\rho }_{\mathit{ta}})}^{1/3}\times R_{200}.\end{array}$$(D.1)

Here β_ta is the scaling factor between the "collapsed mass" M₂₀₀ and M_ta (Eq. 11, Korkidis & Pavlidou (2024) for X ≡ X_ta) while ρ₂₀₀ is the density at δ = 200.

In Fig. D.1, we test that this scaling relation exists in simulations by plotting the actual turnaround radius R_ta,actual as a function of R₂₀₀ for our MDPL2 cluster sample at z = 0. The median along with the 16th and 84th percentiles of this scaling in bins of R₂₀₀, is shown with the solid blue line and the shaded region, respectively. When comparing the median with the scaling relation in Eq. D.1, depicted in red,we observe a strong correlation between the two, with a median offset close to 5%.

All Figures

	Fig. 1. Dark matter density (in units of the background-Universe average matter density) plotted against the normalized radius of concentric spherical shells around MDPL2, z = 0 halos. The color scale indicates the number of lines per pixel on the plot. The radii of the shells have been normalised with respect to the kinetically defined turnaround radius.
In the text

Fig. 2. Comparison between simulated and analytic density profiles for different redshifts. Upper panel: solid lines depict the median density profile at different redshifts of our MDPL2 sample in spherical shells of radius R/Rta; shaded regions represent the 1σ spread of the density within each shell; dashed lines represent the analytical prediction for the outer portions of the mode profile discussed in Sect. 3. Middle panel: Logarithmic slope of the density profiles, with dashed lines illustrating the logarithmic slope of the analytical prediction. Lower panel: ratio of logarithmic slope in median simulated profile over the one of the analytically predicted profile. The analytic slope deviates from the one in simulations for R/Rta ≲ 1.

In the text

Fig. 3. Dependence of the analytical profile and its logarithmic slope on the cosmological parameters. Upper panels illustrate matter density across shells, normalized by the turnaround radius, at redshifts indicated by the titles. With different line colors we show how the profiles change for different values of Ωm (assuming ΩΛ = 0.7) and ΩΛ (assuming Ωm = 0.3). Bottom panels explore the change in the logarithmic slope of the respective density profiles.

In the text

	Fig. B.1. Effect of choosing different statistics to represent the density profile. The two panels present the mode, median and mean density profiles along with their logarithmic slopes of our halo sample at z = 0.
In the text

	Fig. B.2. Analytical profile dependencies on the cosmological parameters at different redshifts. As in Fig. 3, different panels show the change of the density profile and its logarithmic slope with the difference that in this figure the radii have been normalized with respect to R200.
In the text

	Fig. D.1. Correlation between the true turnaround radius Rta,actual and R200 for our cluster sample at z = 0. The red solid line depicts the theoretical scaling relation of Eq. D.1. The blue solid line represents the median value of the blue points in bins of R200. The blue shaded region corresponds to the 16th and 84th percentiles.
In the text