© The Authors 2024

1 Introduction

Galaxy clusters, the largest gravitationally bound structures in the universe, have long been recognized as valuable cosmological laboratories. During the past three decades, significant advances have been made in understanding their composition, dynamics, and integral role within the cosmic web. The well-studied interiors of these clusters, particularly their relaxed regions, have laid the foundations for our concordance model of hierarchical structure formation. This progress has now extended beyond the traditionally studied virialization regime to include the kinetically driven splashback regions (Diemer & Kravtsov 2014; Adhikari et al. 2014; More et al. 2015, 2016; O’Neil et al. 2021). However, the potential of the outer regions of galaxy clusters, extending outside the splashback radius, into the infalling region and beyond, remains largely untapped. These outer reaches hold substantial promise for deepening our understanding of structure formation and cosmology.

Recently, the turnaround scale has received significant attention in cosmological studies (Pavlidou & Tomaras 2014; Lee & Li 2017; Fong et al. 2022; Lopes et al. 2019; Capozziello et al. 2019). This scale represents the point at which galaxies transition from infall towards the central cluster to expansion with the background Universe. In previous work (Pavlidou et al. 2020), we identified the turnaround scale as a novel cosmological probe. We showed that the present average matter density on this scale (the turnaround density ρ_ta) probes the overall matter content of the universe, and that its evolution with time is influenced by the existence of a cosmological constant. The turnaround scale has thus the potential to provide new constraints on cosmological parameters, complementary to the ones obtained from the Cosmic Microwave Background (CMB) (e.g., Komatsu et al. 2011), Baryon Acoustic Oscillations (BAO)(e.g., Eisenstein et al. 2005), and Type Ia Supernovae (e.g., Amanullah et al. 2010).

Our subsequent studies (Korkidis et al. 2020, 2023) confirmed the utility of ρ_ta for testing cosmological models using N-body cosmological simulations. However, measuring the turnaround density in actual observations using current astronomical surveys presents significant challenges. The radius where galaxies join the Hubble flow has thus far been measured for only a few nearby superclusters (Karachentsev & Nasonova 2010; Nasonova et al. 2011; Lee 2018). Accurately measuring the total mass of galaxies at these scales presents an even more formidable challenge. Several issues complicate such observations on a cluster-by-cluster basis, making the associated errors unsuitable for precise cosmological analysis: (a) Accurate measurement of the mass at the turnaround scale would necessitate extensive spectroscopic and gravitational lensing surveys, mapping the galaxy distribution on very large scales around clusters. (b) In the context of the cold dark matter (CDM) model of structure formation, it is well-acknowledged that most of the mass surrounding galaxies is dark and, for such small objects as galactic halos, largely inaccessible; (c) Foreground and background galaxy contamination at these scales is non-negligible and can skew observations.

In this context, developing scaling relations that would allow us to infer the mass at turnaround scales from well-established observable masses becomes paramount. Current observational techniques, such as X-ray, Sunyaev-Zel’dovich, and weak-lensing surveys (Andrade-Santos et al. 2021; Bahar et al. 2022; Planck Collaboration Int. III 2013; Sifón et al. 2016; Hilton et al. 2021, 2018; Ruppin et al. 2018; Gruen et al. 2014), in conjunction with machine learning methods (see, e.g., Ho et al. 2019; Gupta & Reichardt 2020, 2021; Armitage et al. 2019a,b), routinely identify clusters and measure fixed-overdensity masses (masses on a scale where the average density of the cluster is a set multiple of the the mean matter density in the universe). These overdensity masses scale well with each other, leading to the plausible hypothesis that a straightforward scaling relation might exist among all overdensity masses. The turnaround mass is also a fixed-overdensity mass (for instance, in a ΛCDM cosmology at z = 0, the turnaround mass corresponds to an overdensity of δ = 11, see for example Korkidis et al. 2020). Identifying the origin of these scalings could thus provide a straightforward path towards establishing scalings for the turnaround mass as well.

One path to developing a predictive model for these scaling relations may pass through the the mass profile of cluster halos. Much of the effort in this field focuses on the cluster interior, resulting in the widespread use of common profiles like NFW or Einasto (Navarro et al. 1996; Einasto 1965). These profiles were developed to describe the innermost parts of such structures through extensive studies in cosmological simulations. More recently, Diemer (2023) has proposed a model that fits robustly the entire stacked density profile of simulated clusters, both interiors and exteriors. In practice, such comprehensive profiles could be leveraged to construct scaling relations.

In this study, we do exactly this, making use of a model for the mass profile of cluster exteriors derived analytically from first principles. The advantages of using an analytic model for the profile are twofold. First, the derivation process elucidates the physical processes that shape the profile. Second, the profile parameters are directly and causally relatable to cosmological and structure-formation parameters, without the need for expensive runs of large suits of cosmological simulations.

To derive an analytic relation for the (most probable) outer density profile of galaxy clusters, we make use of the model for the density distribution around dark matter halos proposed by Pavlidou & Fields (2005). Our aim is to arrive to an explicit relation between the overdensity mass in the inner, relaxed portions of cluster-sized halos and the overdensity mass in their outer regions, including the turnaround scale.

This paper is organized as follows: In Section 2, we develop an analytic formulation for the most-probable outer mass profile of galaxy clusters, and we illustrate how mass scaling relations can be directly derived from this formulation. In Section 3, we detail the cosmological simulations employed in our study to test our analytic profile, and describe the methodology used to evaluate the performance of our analytical models. In Section 4 we compare our derived profile with those obtained from simulations and test the validity of the associated mass scaling relations. Finally, in Section 5, we summarize our findings.

2 An analytic closed-form profile for cluster exteriors

To build a model of the density distribution around cosmic structures of varying masses, Pavlidou & Fields (2005) used excursion set theory to derive a joint distribution of structures with respect to mass and surrounding overdensity, thus providing a general, two-parameter statistical description. This double distribution extends the Press-Schechter formalism by introducing a clustering scale parameter, β, which quantitatively defines the “environment” around a structure of mass m to be a scale that includes mass βm (including the structure itself, so β > 1). By providing the statistics of the overdensity field as a function of enclosed mass, the double distribution offers a path to deriving the most probable density profile at large distances away from cosmic structures – in other words, to understanding how interior cluster profiles merge into the background density of the expanding Universe.

Models of the outer average density profile of clusters from first principles were pioneered by Barkana (2004); Prada et al. (2006); Betancort-Rijo et al. (2006); Tavio et al. (2008). Our methodology, while not fundamentally divergent from theirs, differs in four key choices. Firstly, in contrast to their use of a “Lagrangian variable q” – essentially the enclosed mass in units of radius – as an independent variable, we use the enclosed mass itself. This choice not only allows for a more intuitive understanding of our derivations, but, most importantly, straightforwardly leads to usable mass scaling relations. Secondly, we employ a much simpler, while still accurate, approximation for the transformation of linear-theory overdensities to spherical-collapse ones, which was proposed by Pavlidou & Fields (2005). This choice simplifies the calculations and results in a more transparent closed-form expression. Thirdly, our focus is primarily on the largest structures. This specificity allows us to omit, with minimal error, corrections for structures engulfed by larger ones, given the immense size and rarity of our targeted structures. Finally, we approximate the variance of the density field, S(m), with a power law of the smoothing mass scale, m. This approximation holds when concentrating on a limited range of masses, as we do in our study, and it is critical in demonstrating the (quasi-)universality of our derived profile. The cumulative effect of these modifications is a significant simplification of the mathematical framework without a significant compromise in predictive power.

Building upon the work by Pavlidou & Fields (2005), we extend their analysis by deriving the most probable density profile as a function of collapsed mass from the double distribution. In what follows, we summarize the mathematical formalism and the assumptions employed in this derivation.

2.1 The double distribution

In Pavlidou & Fields (2005), a joint distribution is constructed to describe the frequency, within a given volume, of collapsed structures in specific intervals of mass and local overdensity. The overdensity is defined as the density contrast with respect to the average matter-density of the background Universe, calculated on a smoothing scale enclosing mass βm (including the central structure itself). It is labelled by the corresponding linearly extrapolated to the present time overdensity, denoted by . Mathematically, then, the distribution provides the comoving number density of collapsed, relaxed halos with mass between [m, m + dm], residing within local overdensity in the range . The linearly extrapolated overdensity field corresponds to the overdensity field that would result if all structures continued to grow according to linear theory until the present time. The mathematical representation for the double distribution is the following: (1)

In this expression, (a) is the linearly extrapolated overdensity of a structure collapsing at redshift z (scale factor a = 1/(1 + z)), extrapolated to today; is the linearly extrapolated overdensity of a sphere enclosing the central collapsed structure, extrapolated to the present time; and S(m) is the variance of the density field when smoothed on a scale that encloses mass m. One implicit assumption made in Eq. (1) is that for steps in mass Δm, the density variance ΔS(m) (and hence, the density contrasts ) are uncorrelated (see Appendix A, for a more involved discussion).

The use of linearly extrapolating overdensities in the derivation of Eq. (1) is necessary for the treatment of halo statistics to remain analytically tractable. However, before Eq. (1) can become predictive for structures residing in the real, non-linear density field, we have to relate these linearly extrapolated over-densities to their non-linear counterparts. To this end, we use a relation introduced in Pavlidou & Fields (2005) which connects linearly extrapolated overdensities to overdensities δ determined using the spherical collapse model: (2)

In Eq. (2), is an overdensity linearly extrapolated to time a, δ_a is the real overdensity at that same time a, and is the overdensity of a collapsing structure linearly extrapolated to the time of collapse. The central assumption of the spherical collapse model in that the matter surrounding each overdensity is distributed isotropically and follows radial motion during its collapse. This relation is universal for all structures and all times of collapse. Eq. (2) has the correct asymptotic behavior both at low and high values of δ_a and is accurate to few percent throughout its domain.

2.2 From the double distribution to the most probable profile

Another way to state the content of Eq. (1) is that, for a given central collapsed structure of mass m at a cosmic epoch a, the double distribution describes the probability distribution of over-densities around structures of that mass m, in spheres of increasing mass βm. We can then derive a profile (a function of βm) of maximum-probability (linearly extrapolated) overdensity – the maximum probability profile. Mathematically, this profile is simply derived by differentiating the double distribution with respect to overdensity, and setting that derivative to zero.

In what follows, we describe a series of simplifications, both to the double distribution itself, and to the profile derivation process, that result in a simple, intuitive, closed-form profile, without significant loss of accuracy.

Simplification 1: dropping the structure-in-structure correction. In Eq. (1), the factor:

represents the fraction of points in space that, when smoothed on a scale βm, have a overdensity . The role of the second exponential in the fraction is to correct for points in space that are in fact part of a collapsed structure on a smoothing scale larger than βm. In random walk theory, this is referred to as the existence of an absorbing boundary at the collapse overdensity . For scales much larger than galaxy clusters themselves, it can be shown that the second exponential, already for β > 1.2 (for reference, turnaround is predicted to be around β > 1.6) is less than 20% of the first term. Dropping it can significantly simplify the final expression, and we do so for the remainder of this paper.

Simplification 2: first differentiate, then convert to nonlinear overdensities. As discussed earlier, Eq. (1) depends on the linearly extrapolated overdensity. Strictly speaking, in order to derive the most probable profile, we must first perform a change of variables using Eq. (2) to recast the double distribution in terms of δ_ℓ and then differentiate to obtain the value of δ_ℓ as a function of βm where the distribution reaches its maximum. However that would significantly complicate both the algebra involved, and the final resulting expression. Instead, we will adopt the approximate path of first calculating a most-probable profile of linearly extrapolated over-densities by differentiating Eq. (1), and then use the conversion relation (Eq. (2)) to convert the profile of linearly extrapolated overdensity to profile of non-linear overdensity. Physically, this would be equivalent to determining what the most probable density profile around a structure that has collapsed today looked like in the early Universe, when the density field was still in the linear regime; and then, evolving that early most-probable profile forward in time, using the spherical collapse model. We expect some inaccuracy – which we will evaluate through comparison with simulations – to stem from this choice, since the “most-probable profile” is a statistical quantity and does not represent, as a whole, the density profile around a single real structure (see discussion of the “typical profile” in Betancort-Rijo et al. 2006, and also in Barkana 2004 for a discussion of similar strategies).

Differentiating Eq. (1) with respect to , using the simplifications discussed above, and setting the result equal to zero returns the following simple expression for the most probable profile of overdensity linearly extrapolated to the present cosmic time: (3)

We note that is the overdensity of a structure collapsing at time a, linerarly extrapolated to a = 1 (today). It is related to (the overdensity of a structure collapsing at a extrapolated to time a, a universal value for all epochs a) through , where D(a) is the linear growth factor at redshift a. Using this relation, we can rewrite Eq. (3) in terms of quantities extrapolated to the time of collapse of the central structure a: (4)

Finally, we use Eq. (2) to convert the to δ_ℓ(a). This is now possible, because both the linearly extrapolated and the nonlinear (spherical collapse) overdensities are calculated at the same time a. This yields (5)

Solving for δ_ℓ we get (6)

which we can confirm has the correct asymptotic behavior, going to ∞ when β = 1, and going to 1 (i.e., δ_ℓ = 0, average matter-density of the universe) when β → ∞. We can also express this equation in terms of average densities within a sphere including mass βm (instead of overdensities), (7)

where ρ_m is the average matter-density of the Universe at the desired redshift.

We can simplify Eq. (7) even further if the variance of the density field can be approximated as a power law in mass, S(m) = m^−γ, with −γ = d ln S/d ln m (in Appendix A we explore the legitimacy of this assumption). With this substitution, we obtain (8)

Put in words, up to the validity of all the approximations we have discussed so far, the average matter density of a sphere enclosing mass β times central galaxy cluster mass (so β > 1 always), normalized to the mean matter-density of the universe at the time of observation, is universal, and equal to . For very large β, this profile asymptotes to 1 (i.e., eventually the matter density of increasingly large spheres encompassing the cluster eventually tends to the average matter-density of the Universe at the time of observation). For β → 1, the profile asymptotes to ∞, as is the expectation in the spherical collapse model, where a structure collapses to a singularity. In practice, this means that the validity of this (outer-region) profile will break down at some β > 1, to be determined by comparison with simulations.

2.3 Mass scalings

The very existence of a universal average density profile dependent only on β (and not on m) implies that the virial mass will scale with any other “constant overdensity” mass that is larger than the virial mass itself. If, then, we take the virial (collapsed, relaxed) mass of the central cluster to be the usually-assumed m₂₀₀ (i.e., a mass enclosed within a sphere which is, on average, 200 times denser than the average background-matter density of the Universe), and m_X to be a different fixed-overdensity mass (a mass enclosed in a sphere X times overdense with respect to the background Universe, with X < 200), it is straightforward to derive a scaling between the two. In order to see this we solve Eq. (8) for β to get: (9)

So now any constant overdensity criterion ρX = Xρ_m(a) will yield a specific value for β (a value independent of m, and so the same for clusters of all masses): (10)

And thus the sphere of average overdensity ρX = Xρ_m(a) will most probably contain a mass β_Xm (“most probably” because our formalism yields the most probable average density profile). Then, the most probable scaling of m_X with m (i.e., most commonly, with M₂₀₀), will be (11)

This will also hold for the turnaround mass, since the turnaround density is a constant for a given cosmology and a given redshift (for example, for the concordance ΛCDM cosmology today, X = 11 for the turnaround mass). (12)

This means that we have a theoretically predictable scaling between M₂₀₀ and M_ta, for any cosmology. Because γ has a slight dependence on mass, the logarithmic scaling between masses is expected to have some tilt compared to the the M_ta ∝ M₂₀₀ line.

Similarly, we can derive a scaling between any two overdensity masses, say m_X and m_Y, defined as the masses of spheres with average matter density Xρ_m(a) and Yρ_m(a), respectively, provided that both X and Y are lower than 200 (or whatever the threshold for the virialized part of the cluster is taken to be): (13)

Detailed recipes for deriving the values of γ and used in this work are given in Appendices A and B, respectively.

3 Simulations and methods

In this section, we discuss the simulations used for testing the validity and accuracy of the analytic profile and the predicted mass scalings. We present the criteria we adopted to select our halo sample, and the methodology we used to construct the most-probable average mass profiles for these halos.

3.1 Simulation overview

To validate our analytical models against realistic, simulated average mass profiles of dark matter halos in cluster-sized systems, we used the same cosmological simulations as in our prior work (Korkidis et al. 2023; a more in depth presentation of the simulation can be found therein). This includes the large-box concordance ΛCDM MDPL2 simulation (1000 h⁻¹ comoving Mpc on a side, particle mass 1.51 × 10⁹h⁻¹M_⊙), and the Virgo consortium simulations, with identical 239.5 h⁻¹ Mpc – sided boxes for three different cosmologies [a concordance ΛCDM, a no-Λ Ω_m = 0.3 (OCDM), and a “standard” Ω_m = 1 flat CDM (SCDM)], and particle masses equal to 6.86 × 10¹⁰ h⁻¹ M_⊙ for the ΛCDM/OCDM, and 22.7 × 10¹⁰ h⁻¹ M_⊙ for the SCDM respectively.

Fig. 1 Average matter density (in units of the background-Universe average matter density) within concentric spheres around MDPL2 z = 0 cluster-sized halos, as a function of mass that is enclosed in each sphere. The colorscale encodes the number of lines per pixel on the plot. The enclosed mass is normalized to M200, so all profiles converge at an average density 200 per the definition of M200. The red solid line depicts the profile mode. The independent variable is equivalent to β in the analytic profile of Eq. (8).

3.2 Building the profiles

The profile that we derived in Section 2, described the most probable profile of the average matter density within concentric spheres of increasing enclosed mass, as a function of that enclosed mass. To compare the analytical predictions against the simulated data, we analyzed the matter distribution around 1000/900 randomly selected galaxy-cluster-sized (M₂₀₀ ≥ 8 × 10¹³ M_⊙) dark matter structures from each of the MDPL2/Virgo simulations, respectively. Given that in hierarchical structure formation cosmologies the halo distribution follows a mass function in which larger structures are more rare, during our sampling we divided the halo catalog in 30 mass bins, and from each bin, we selected 40 random clusters. We segmented the region surrounding the center of each galaxy cluster into 500 concentric spheres extending up to 10 × R_200,m in radius, ensuring that each halo was contained within the simulation’s boundaries. Within each sphere, we computed the enclosed mass, and, dividing by the volume of the sphere, the corresponding dark matter density. The chosen number and size of bins were sufficiently large to accurately represent the profiles at all enclosed overdensity masses for both group and cluster-sized halos. Finally, we normalized the enclosed mass for each individual halo profile to its virial mass, specifically M_200.

4 Results

4.1 Mass density profile of simulated halos

We first confirmed that the type of profiles examined here do cluster around a most probable behavior which our analytical profile aims to model. Figure 1 presents the compiled profiles from the MDPL2 simulation at z = 0. The colorbar indicates the density of lines on the plot. This visualization clearly shows that most profiles tend to converge around a common profile. In Fig. 1, the independent variable is equivalent to β in the analytic profile of Eq. (8), assuming that M₂₀₀ is a reasonable approximation of the true relaxed mass m of the central cluster.

To formally capture and represent this pattern, we grouped these normalized profiles into 40 linear mass bins, from 0 to 4 × M_200. For each bin, we determined the most probable value (mode) of the density. This was achieved by employing a Gaussian kernel density estimator¹ (Virtanen et al. 2020) to model the probability density function (PDF) of densities in the bin, and subsequently identifying the density for which the PDF is maximum. The resulting mode profile is depicted by the red solid line in Fig. A.1, closely tracking regions with high profile density, as anticipated.

We next tested the extent to which our analytic profile matches the most probable profile seen in simulations. In Fig. 2, we overlay the analytical profile with the mode of the profiles derived from the MDPL2 simulation. Different curve colors stand for different redshifts of our halo samples. For the analytical profile, we incorporate the cosmological parameters from the simulation (Ω_m, Ω_Λ, h₀, n_s) and the median value of M₂₀₀ for the halo sample (6.3 × 10¹⁴ M_⊙ for z = 0 halos and 1.2 × 10¹⁴ M_⊙ for z = 1). The lower panel compares the simulated and analytical profiles by examining their ratio.

The analytical profiles, developed using the spherical collapse model to evolve linearly extrapolated overdensities into the non-linear regime, tend to infinity as the enclosed mass approaches M₂₀₀². For , the analytical profile shows an excellent alignment with the mode profile; fluctuations are within 15%. The profile exhibits little variation with redshift, beyond what is encoded in the overall increase of average background matter density with increasing redshift.

Simulations consistently show that the density profile around clusters exhibits near-universal characteristics when normalized against a constant overdensity radius or mass relative to the background. This behavior is replicated in Fig. 3, where, with solid lines, we plot the mode density profiles of MDPL2 clusters at z = 0, against their enclosed mass in concentric spheres of increasing radius. Different colors correspond to different average masses of the halo sample. The consistency of the profiles between different mass bins is even stronger in simulations than in our analytical profiles (shown with dashed lines), with deviations between simulations and analytic profile for the smaller and larger mass bins reaching 25%. This variation of the analytic profile with cluster mass originates in the mild dependence of γ on mass (each profile in the figure has a different γ value as shown in the legend; also, see Appendix A and Fig. A.1). We have confirmed that the expected mild mass-dependence of the analytic mode profile is real, and not a result of our approximations: it persists when the full expression for S(m) is used rather than the power-law approximation: the deviation between the two profiles is negligible for the range of considered here, when γ is taken to be equal to (as we do in all cases, using the median M₂₀₀ of our sample each time). For much higher , of course, when the transfer function T(k) → 1, S(m) does asymptote to a power-law (see Eq. (A.3)), and the full universality of the mode profile is recovered; however this becomes relevant for mass scales much larger than those considered in this paper. The dependence of the mode profile also persists if we replace the approximation of Eq. (2) with the more accurate but also more cumbersome Eq. (8) of Sheth & Tormen (2002), also used by Prada et al. (2006).

The lack of such dependence in simulations is plausibly an effect of the somewhat fuzzy correspondence between M₂₀₀ and the true “collapsed mass” entering the analytical mode profile. Even under the assumption of spherical symmetry, the SCM is only applicable for shells that have never undergone shell crossing. In this context, the appropriate collapsed mass is somewhat ambiguous. Here we have taken it to be M₂₀₀. However, if one were to require that all mass that has undergone shell crossing is considered to be part of the collapsed mass, then the splashback mass (Diemer & Kravtsov 2014; Adhikari et al. 2014; Diemer et al. 2017) should also be considered. This is likely one of the effects that partly contribute to the offset observed in Fig. 3 between the solid and the dashed blue curves.

This effect would also be exaggerated by the underlying assumption in Eq. (1) that steps in enclosed mass result in uncorrelated corresponding steps in enclosed underdensity (implemented formally by the adoption of a sharp-in-k filter for S(m); see also discussion in Appendix A). This assumption is required for the excursion-set formalism used to derive Eq. (1) to be strictly applicable; however, in the analysis of simulations, successive steps in overdensity corresponding to successive steps enclosed mass are correlated (they are calculated using a top-hat, rather than a sharp-in-k, filter). The extent of this mismatch would be dependent on the mass accretion rate.

Finally, another factor to consider is the finite accuracy with which we determine the mode density in every sphere of a given : the mode profile remains reasonably within uncertainties of the simulated mode across mass bins, simulations, and redshifts.

Fig. 2 Performance of our analytic density profile for different redshifts. Upper panel: most probable (mode) average matter density within spheres of enclosed mass M as a function of M (normalized to M200; is equivalent to β in Eq. (8)), from our analytic result (Eq. (8), dashed lines), and from the MDPL2 simulation halos (solid lines; the shaded region includes 68% of the densities PDF), for z = 0 (blue) and z = 1 (green). Lower panel: ratio of MDPL2 over analytic mode profile as a function of enclosed mass. In the bottom panel, the ratio of the mode to the analytic profiles is displayed.

Fig. 3 Most probable average matter-density profile from our analytical model (dashed lines) and MDP12 (solid lines), for different ranges of cluster masses M200. The independent variable is equivalent to β in Eq. (8).

Fig. 4 Correlation between turnaround mass Mta,actual and collapsed mass M200 for the MDPL2 cluster sample at z = 0. The red solid line depicts the theoretical scaling relation of Eq. (11). The blue solid line represents the mode value of the blue points in bins of M200. The blue shaded region corresponds to the 16th and 84th percentiles.

4.2 Mass scalings

We now turn to evaluate the performance of our scaling relation between different overdensity masses (Eq. 11) for our halo samples. As detailed in our prior works (Korkidis et al. 2020, 2023), the kinematically-defined turnaround mass M_ta,kinematic (the mass enclosed within the kinematically-identified turnaround radius) corresponds well to an overdensity mass, consistent with the prediction of the spherical collapse model. For concordance ΛCDM and z = 0, spherical collapse predicts that X_ta = 11. We also set the collapsed mass m = M_200. For each halo, employing Eq. (11), we can predict the turnaround mass M_ta,predicted ≡ M_11,predicted from its M₂₀₀ and also calculate its actual mass within that overdensity (X = 11, M_11,actual ≈ M_ta,actual).

We first test that such a scaling relation does exist in simulations, for the case that is of most interest to us (scaling of M₂₀₀ with the turnaround mass). In Fig. 4 we plot the turnaround mass of each cluster in our MDPL2 sample, M_ta,actual, as a function of its M₂₀₀. The mode along with the 16th and 84th percentiles of this scaling in bins of M₂₀₀ is shown with the solid blue line and the shaded region, respectively, and they do confirm that a tight scaling between the two masses is indeed present in the MDPL2 sample, although outliers do exist. For comparison, we overplot, with the red line, the mode scaling predicted by Eq. (11), and it is already obvious that the analytic mode scaling is in excellent agreement with the mode scaling seen in simulations.

To further quantify the performance of our scaling we plot, with the histogram in the left panel of Fig. 5, the fractional difference, as a percentage, between the turnaround mass M_11,actual measured from the simulation data, and the turnaround mass M_11,predicted predicted by our scaling relation using M₂₀₀ as input. The histogram has a pronounced peak near zero, suggesting that the scaling relation predicts the turnaround mass well. However, the distribution of residuals does feature a long tail, so the mean of this distribution has a non-negligible bias (up to 20%).

Another factor affecting the performance of the scaling relation is the mass range of the sample. This is explored in the right panel of Figure 5. Here, different colors correspond to subsamples of different mass ranges, with the median mass in each subsample shown in the legend. Although the qualitative features of the histograms remain consistent, we observe a systematic trend from negative to positive offsets of the peak of the distribution as the median sample mass increases. The absolute value of this offset is up to 15% for the range of masses considered here. This behavior is anticipated based on the performance of the analytic profile for different sample mass ranges seen in Fig. 3.

In Fig. 6 we evaluate again the performance of our mass scaling for different values of overdensity δ at which the mass is evaluated (corresponding to X in Eq. (11)). Each panel corresponds to a different value of δ, which decreases from left to right, and which is displayed in the legend of each panel. In all scenarios, the distribution peaks at < 10% error. The most probable error increases with increasing overdensity, a behavior we anticipate due to the analytic profile diverging at the assumed collapse overdensity δ = 200. However, the high-error outlier tails decrease with increasing overdensity, a behavior also anticipated from Fig. 1: closer to M₂₀₀ individual profiles are more clustered around the mode profile. For this reason, deviations of individual outliers from the analytic mode profile, and the associated deviations from the mass scaling derived from the mode profile, are more pronounced for higher mass scales (lower overdensity values).

The same analysis for X ≡ X_ta (turnaround mass) is then conducted in Fig. 7 for halo samples from three Virgo simulations, each corresponding to a different set of cosmological parameters: Ω_m = 0.3, Ω_Λ = 0.7 (VirgoLCDM); Ω_m = 0.3, Ω_Λ = 0 (VirgoOCDM); and Ω_m = 1, Ω_Λ = 0 (VirgoSCDM). In all cases, we analyzed z = 0 snapshots. The behavior in the first two cosmologies, featuring the same dark matter content at the present cosmic epoch, is very similar. This is consistent with our findings both analytically in the context of the spherical collapse model (Pavlidou et al. 2020) and in simulations (Korkidis et al. 2023) that the behavior on turnaround scales at some specific redshift depends primarily on the average matter density at that same redshift. The mild difference of the most-probable error compared to the MDPL2 sample (Fig. 5) likely stems from the smaller box size in the Virgo simulations, which results to a smaller median mass of the halo sample. As also seen in Figure 5, in this range of halo masses the scaling tends to overpredict the turnaround mass (negative mode error). In the VirgoSCDM box, where the median mass of the sample is larger, the error mode is almost zero, and a more pronounced high-error tail is present, also consistent with our findings for different cluster mass ranges in MDPL2. It would appear that sample mass drives the performance of the scaling relation of Eq. (11), more than the variation in cosmological parameters, which appear to be accounted for sufficiently by the model.

Finally, we also test the performance of Eq. (13) – a scaling relationship between any two overdensity masses m_X and m_Y. We employ this scaling relation to estimate the actual value of the turnaround mass M_ta,actual (X = 11) using an overdensity mass M₁₀₀ (Y = 100), and plot the distribution of its percentage difference from M_ta,actual in Fig. 8. In the same figure, we also overplot the same distribution for the prediction derived from Eq. (11). The comparison indicates an almost identical position for the peak of the distributions; however, when the turnaround mass is derived from M₁₀₀ rather than M₂₀₀ the error distribution is tighter around zero – a result of the analytic profile performing better away from its divergence point at the collapse mass.

Fig. 5 Probability density of the fractional (percentage) difference between the turnaround mass (here M11) predicted from the scaling relation of Eq. (11) and its actual value, for the full MDPL2 simulated cluster sample at z = 0 (left panel), and for halo subsamples of different median M200 (right panel). In each case, the solid line shows a Gaussian kernel estimator of the probability density function, and the vertical line marks the mode of the error distribution. The scaling relation tends to overpredict M11 for lower M200 and underpredict it for higher M200.

Fig. 6 Performance of our scaling relation for masses corresponding to different overdensity values δ (equivalent to X in Eq. (11)) in the MDPL2 sample. The solid curves and vertical lines are the same as in Fig. 5. As δ decreases (left to right), the most probable offset decreases, while high-error tails become more pronounced. Both behaviors can be traced to the performance of the analytic profile (see text for discussion).

Fig. 7 Performance of our scaling relation for three different Virgo simulations cosmologies. All panels correspond to z = 0 snapshots. The median M200 of the halo samples are 1.2 × 1014 M⊙ for VirgoLCDM/VirgoOCDM and 1.7 × 1014 M⊙ for VirgoSCDM. Dashed red lines correspond to ±15% error.

Fig. 8 Comparison between the scaling relations delineated in Eqs. (13) and (11). The blue histogram depicted herein corresponds to that of the left panel in Fig. 5. Superimposed on this histogram is the fractional difference between the actual value of Mta (X = 11) and its predicted counterpart assuming a prior knowledge of M100 (that is, if we set Y = 100).

5 Conclusions

In this study, we derived a theoretical model for the most-probable outer average density profile of large galaxy clusters as a function of enclosed mass, using excursion set theory. This model, based on Gaussian early universe statistics and the spherical collapse framework, has two parameters, γ and , both of which can be analytically derived from the cosmological parameters (Ω_m, Ω_Λ, and power spectrum slope n), and the median virial mass of the sample under consideration. We compared our profile with profiles around galaxy-cluster–sized halos in ΛCDM simulations, and found good agreement across mass ranges and redshifts.

From the profile we derived a scaling relation that links different overdensity masses, and found that it performs well with the the peak of the error distribution remaining below 15% for all cosmologies, sample masses, and overdensities we tested. We have traced residual inaccuracies of the scaling relation to two factors.

The offset from zero of the maximum of the error is dependent on the exact cluster mass range considered. In particular, the mode profile in simulations appears to be more universal (similar across mass ranges) than the analytic model predicts. This is plausibly an effect of an imperfect identification between the collapsed mass m of the analytic profile and the overdensity-200 mass M₂₀₀ in simulated clusters, leading to a cross-contamination of mass ranges.

The tail of the error distribution tends to grow as the mass scale we consider grows (or, equivalently, as the overdensity where the mass is calculated decreases). This is due the increasing spread of the profile distribution with increasing scale. As a result, masses at overdensities closer to 200 (say, 150 or 100) scale more tightly with M₂₀₀ than the turnaround mass (M₁₁ for z = 0 in concordance Λ CDM).

Our analytical model demonstrates the existence, in principle, of a scaling between any two overdensity masses, dependent on cosmology and redshift, and it provides a reasonable approximation to this scaling directly from these parameters. A fit of the parameter γ to a particular simulated sample of clusters with masses in the range of interest can provide an even more accurate scaling should one be required. We have confirmed that such fits return values of γ that always fall within those predicted theoretically (see Appendix A and Fig. A.1) for masses in the galaxy-cluster range.

Acknowledgements

We would like to thank the anonymous referee for his constructive report. 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); 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 no. 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. (www.gauss-centre.eu) and the Partnership for Advanced Supercomputing in Europe (PRACE, www.prace-ri.eu) for funding the MultiDark simulation project by providing computing time on the GCS Supercomputer SuperMUC at Leibniz Supercomputing Centre (LRZ, 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. The Virgo Consortium simulations used in this paper were carried out by the Virgo Supercomputing Consortium using computers based at Computing Centre of the Max-Planck Society in Garching and at the Edinburgh Parallel Computing Centre. The data are publicly available at www.mpa-garching.mpg.de/galform/virgo/int_sims. Throughout this work we relied extensively on the PYTHON packages Numpy (Harris et al. 2020), Scipy (Virtanen et al. 2020) and Matplotlib (Hunter 2019).

Appendix A Density field variance

In deriving Eq. (8) from Eq. (7) we assumed that, for a limited range of masses, S(m) can be approximated by a power law, that is (A.1)

where (A.2)

Here, we evaluate the validity of this assumption, and we provide a recipe for the evaluation of γ. The variance of the density field is given by (A.3)

In this equation, we have made two assumptions. First, that right after matter-radiation equality, ⟨|δ_k|²⟩ can be simply described in terms of a power-law in k modified by a transfer function, ⟨|δ_k|²⟩ ∝ T²(k)kⁿ. Second, that the window function for the calculation of the variance is sharp in k–space. This assumption is necessary for random-walk formalism (used to derive Eq. 1) to be strictly applicable: each “step” in the time-like variable (here S(m)) should produce a step in the space-like variable (here, the overdensity corresponding to the mass scale m) that is completely independent from the previous step. This requires that the k-modes producing a change in the spacelike variable not appear in any of the previous steps. A sharp-in-k window function enforces this condition. For a more extensive discussion on this assumption and its consequences see Bond et al. (1991); Lacey & Cole (1993); Pavlidou & Fields (2005).

The variance of the field has been normalized to the present time at scale of 8 comoving Mpc (σ₈). For T(k) we adopt the fitting formulae of Bardeen et al. (1986) for the adiabatic cold dark matter transfer function (Eq. G3). We have verified that using other approximations (e.g., per Eisenstein & Hu 1998) does not alter substantially any of the results presented in this work.

To obtain k(m), we integrate the sharp-in-k window function over all space and multiplying by pρ_m,0 we obtain (A.4)

Then, from Eq. (A.3), we can calculate S(m) and its logarithmic slope, which we plot in Fig. A.1. Clearly, this slope is not constant over the range of masses of interest (as it should be for a power law), but it does not vary wildly either. As a result, for a small range of masses, S(m) can be approximated reasonably well by a power law. Of course as m → ∞ the transfer function T → 1 and S(m) asymptotes to ∝ m⁻⁽ⁿ+^3)/3, but this occurs on scales much larger than those of interest in this work.

Appendix B Calculating

For a given cosmology, that is, a set of Ω_m, Ω_Λ, we first calculate the parameters (B.1)

Fig. A.1 Variation in the logarithmic slope of the density field dispersion S for diverse values of the collapsed mass m. The function’s minimal curvature indicates near-consistency with a straight line across each order of magnitude in mass. This pattern suggests that approximating the variance of the density profile with a power-law model could be a plausible approach.

and (B.2)

Note that for a flat cosmology, ϕ = 0, while for a cosmology without Λ, ω = 0.

Now let a_p be the (evolving) radius of an initially overdense region (density perturbation) normalized so that, had the specific region begun its evolution at the same density as the bacground Universe, a_p at the present cosmic epoch would have been equal to 1. The size of this perturbation at turnaround, a_pta, that is reaching collapse today (i.e. at scale factor of the background Universe equal to 1; note that this structure must have reached turnaround at sometime in the past), can be found through (e.g., Pavlidou & Fields 2005; Pavlidou et al. 2020)

Then, the extrapolated overdensity to the time of collapse³ will be (Pavlidou & Fields 2005): (B.3)

where D(a₀) is the linear growth factor for the present cosmic epoch a₀ = 1 (Carroll et al. 1992) (B.4)

Overall, the variations in in different cosmologies are small. For concordance ΛCDM (Ω;_m = 0.3, Ω_Λ = 0.7), = 1.6757. For SCDM (Ω_m = 1, Ω_Λ = 0), = 1.6865.

References

All Figures

Fig. 1 Average matter density (in units of the background-Universe average matter density) within concentric spheres around MDPL2 z = 0 cluster-sized halos, as a function of mass that is enclosed in each sphere. The colorscale encodes the number of lines per pixel on the plot. The enclosed mass is normalized to M200, so all profiles converge at an average density 200 per the definition of M200. The red solid line depicts the profile mode. The independent variable is equivalent to β in the analytic profile of Eq. (8).

In the text

Fig. 2 Performance of our analytic density profile for different redshifts. Upper panel: most probable (mode) average matter density within spheres of enclosed mass M as a function of M (normalized to M200; is equivalent to β in Eq. (8)), from our analytic result (Eq. (8), dashed lines), and from the MDPL2 simulation halos (solid lines; the shaded region includes 68% of the densities PDF), for z = 0 (blue) and z = 1 (green). Lower panel: ratio of MDPL2 over analytic mode profile as a function of enclosed mass. In the bottom panel, the ratio of the mode to the analytic profiles is displayed.

In the text

	Fig. 3 Most probable average matter-density profile from our analytical model (dashed lines) and MDP12 (solid lines), for different ranges of cluster masses M200. The independent variable is equivalent to β in Eq. (8).
In the text

	Fig. 4 Correlation between turnaround mass Mta,actual and collapsed mass M200 for the MDPL2 cluster sample at z = 0. The red solid line depicts the theoretical scaling relation of Eq. (11). The blue solid line represents the mode value of the blue points in bins of M200. The blue shaded region corresponds to the 16th and 84th percentiles.
In the text

Fig. 5 Probability density of the fractional (percentage) difference between the turnaround mass (here M11) predicted from the scaling relation of Eq. (11) and its actual value, for the full MDPL2 simulated cluster sample at z = 0 (left panel), and for halo subsamples of different median M200 (right panel). In each case, the solid line shows a Gaussian kernel estimator of the probability density function, and the vertical line marks the mode of the error distribution. The scaling relation tends to overpredict M11 for lower M200 and underpredict it for higher M200.

In the text

Fig. 6 Performance of our scaling relation for masses corresponding to different overdensity values δ (equivalent to X in Eq. (11)) in the MDPL2 sample. The solid curves and vertical lines are the same as in Fig. 5. As δ decreases (left to right), the most probable offset decreases, while high-error tails become more pronounced. Both behaviors can be traced to the performance of the analytic profile (see text for discussion).

In the text

	Fig. 7 Performance of our scaling relation for three different Virgo simulations cosmologies. All panels correspond to z = 0 snapshots. The median M200 of the halo samples are 1.2 × 1014 M⊙ for VirgoLCDM/VirgoOCDM and 1.7 × 1014 M⊙ for VirgoSCDM. Dashed red lines correspond to ±15% error.
In the text

	Fig. 8 Comparison between the scaling relations delineated in Eqs. (13) and (11). The blue histogram depicted herein corresponds to that of the left panel in Fig. 5. Superimposed on this histogram is the fractional difference between the actual value of Mta (X = 11) and its predicted counterpart assuming a prior knowledge of M100 (that is, if we set Y = 100).
In the text

	Fig. A.1 Variation in the logarithmic slope of the density field dispersion S for diverse values of the collapsed mass m. The function’s minimal curvature indicates near-consistency with a straight line across each order of magnitude in mass. This pattern suggests that approximating the variance of the density profile with a power-law model could be a plausible approach.
In the text