© ESO 2020

1. Introduction

The turnaround radius naturally appears within the context of the collapse of a single, spherically symmetric structure in an otherwise homogeneous and isotropic expanding universe as the boundary between the nonexpanding structure and the Hubble flow. In recent years, the turnaround radius has attracted considerable attention (e.g., Pavlidou & Tomaras 2014; Tanoglidis et al. 2015; Lee & Yepes 2016; Bhattacharya & Tomaras 2017; Bhattacharya et al. 2017; Lee 2018; Lopes et al. 2018, 2019; Nojiri et al. 2018; Capozziello et al. 2019) because of its potential as a possible new probe of cosmological parameters, as a constraint on alternative theories of gravity, and as a well-defined boundary for large-scale structures.

The attractiveness of the turnaround radius as a boundary for cosmic structures stems from two factors. The first is its unambiguous definition based on the radial velocity profile, which on the one hand traces the dynamics on the structure, and on the other hand is straightforward to both calculate and explain: the turnaround radius is the scale where the edge of the structure joins the Hubble flow. The second is that structures defined on turnaround scales have only mildly evolved into the nonlinear regime. Thus, their behavior is expected to be closer to the predictions of simple analytic models, such as the “spherical collapse” model (Gunn & Gott 1972; Fillmore & Goldreich 1984; Bertschinger 1985).

The potential of the turnaround radius as a cosmological observable has indeed been primarily explored under the assumption of a single, spherically symmetric structure evolving in an otherwise unperturbed Universe. The “spherical collapse” calculations that follow from this assumption make the following predictions (Pavlidou & Tomaras 2014; Tanoglidis et al. 2016; Pavlidou et al. 2020): (a) The matter enclosed by the turnaround radius has a characteristic average matter density (the “turnaround density”, ρ_ta), which is the same for structures of all masses at a given cosmic epoch. (b) The value of ρ_ta and its evolution with cosmic time depends on (and probes) the cosmological parameters Ω_m and Ω_Λ. (c) At late cosmic times, when Λ fully dominates the dynamics of the expansion, ρ_ta asymptotically approaches a value which is only dependent on the value of the cosmological constant and is equal to 2ρ_Λ, where ρ_Λ = Λc²/8πG. (d) As a consequence, the radius of any nonexpanding structure of mass M in a ΛCDM universe can never exceed (3GM/Λc²)^1/3.

However, before any meaningful comparisons between these predictions of the spherical collapse model and observational data can be made, it is necessary to test whether the predictions persist in simulated structures with realistic shapes, and to quantify any deviations due to departures from spherical symmetry and the presence of neighboring structures. This is the scope of the current paper. In particular, we address the following questions: (1) Can a single turnaround radius meaningfully characterize a realistic 3D structure? (2) If so, how does this turnaround radius compare to the predictions of the spherical collapse model for objects of the same mass? (3) Does the prediction of spherical collapse, that there exists a single characteristic average turnaround density ρ_ta for all structures at a given redshift, persist in N-body simulations? (4) Does the turnaround radius depend on the shape of structures? (5) Does the turnaround radius depend on the presence of massive neighbors?

This paper is structured as follows. In Sect. 2 we describe the set of N-body simulations used in this work and the characteristics of the dark matter halos under study. In Sect. 3 we describe our methodology for the calculation of the turnaround radius and for substructure removal when the turnaround radius is used as the boundary of a structure. In Sect. 4 we compare the turnaround radius in simulated structures to the predictions of the spherical collapse model and investigate the effect of halo shape and environment on the turnaround radius. We discuss these findings in Sect. 5 and conclude in Sect. 6.

2. N-body simulations

In this work we use data from: (a) The Dark Sky Simulations (TDSS; Skillman et al. 2014). The TDSS were run using 2HOT (Warren 2013), a purely adaptive tree-based N-body code. The data extraction and analysis were carried out using yt (Turk et al. 2011; b) The Illustris-TNG simulations (Nelson et al. 2015, 2019), which is a suite of gravo-magnetohydrodynamical simulations run with the AREPO moving mesh code (Springel 2010). For the data extraction of the Illustris-TNG we used the methods presented on the project website¹ along with yt. We list the basic properties of the simulations we use in Table 1. Table 2 shows the cosmological parameters adopted in each simulation.

For producing halo catalogs, TDSS uses the ROCKSTAR algorithm (Behroozi et al. 2013) and Illustris-TNG uses friends-of-friends (FOF; Davis et al. 1985). The FOF algorithm is historically one of the first used to identify halos. It considers particles to be members of the same group if their distance is smaller than a given linking length (Kravtsov & Borgani 2012; Davis et al. 1985). The value of the linking length is usually chosen so that the defined halo has a density contrast that approaches that of a virialized structure (Kravtsov & Borgani 2012; Behroozi et al. 2013; Knebe et al. 2011). ROCKSTAR makes use of the full six-dimensional phase space of particle positions and velocities, as well as of time (Behroozi et al. 2013).

The snapshot that we used from each simulation was z = 0 (the present cosmic epoch).

The resolution provided by the simulations is more than adequate in order to have a sufficient number of particles within turnaround scales, though to be sure that there were no issues, we analyzed halos from the same simulations but with two different mass resolutions (TNG300-2 and TNG300-3). Some special consideration regarding the removal of substructure is also necessary for our analysis. Although each halo finder does include a substructure identification algorithm (Hoffmann et al. 2014), these operate on the virial scale (or a scale that is purposefully selected to be close to the virial scale, in the case of FOF halo finders). However, we are interested in placing the boundary of structures at the much larger turnaround scale. As a consequence, there may be structures that are located outside the virial radii of all nearby larger structures (therefore not considered “substructure” by the halo-finding algorithms), but within the turnaround radius of some structures (and are therefore considered “substructure” for our purposes).

A somewhat similar problem arises for structures in the process of merging. If the centers of mass of two structures are approaching each other, then the two structures are a part of a single “turnaround” structure (they are not receding from each other as the universe expands, hence the system as a whole has detached from the Hubble flow). However, the distance between the two centers may still be large enough that the halo-finding algorithm lists these as two distinct structures. This may result in the same structure appearing twice in our analysis. The larger “turnaround” structure encompassing the system may also have a center of mass that is significantly displaced from the center of either individual structure provided by the halo-finding algorithm. We discuss our strategies for dealing with this issue in Sect. 3.

For computational efficiency, we only analyzed a (randomly chosen) fraction of the structures found in the simulated boxes of TDSS and TNG-300. More specifically, from TDSS we randomly chose 578 halos with M₂₀₀ ≥ 10¹⁵ M_⊙² while from TNG-300-2 (TNG-300-3) we chose 430 (449) halos with masses in the interval 10¹⁴ M_⊙ ≤ M₂₀₀ ≤ 10¹⁵ M_⊙.

3. Calculating the turnaround radius

An important difference of the turnaround radius R_ta from other halo boundary definitions, such as for example the virial radius or the splashback radius (Diemer & Kravtsov 2014; More et al. 2015), is that R_ta is kinematically identified³. For this, we use the radial profile of ⟨V_r⟩ (the average radial velocity around a halo center of particles residing within a spherical shell). In our analysis, spherical shells are logarithmically spaced, with the outer radius of each shell for TNG being 0.68% (0.78% for TDSS) larger than the inner radius of the shell. The obtained value of ⟨V_r⟩ in each shell is then assigned to the average distance of all particles in that shell from the halo center. The uncertainty of ⟨V_r⟩ in each bin is taken to be the standard error of the mean.

The radial profile of ⟨V_r⟩ exhibits a common pattern for every halo: beyond some distance from the halo center, matter follows the Hubble flow, making the average radial particle velocities positive thereafter. This typical behavior is shown in Fig. 1, where radial velocity profiles are plotted for one halo from each simulation. Different halos, even in the same mass range, can feature different behavior in the spherically averaged radial velocities in their inner and intermediate parts, depending on their relaxation state and merging history (Cuesta et al. 2008); however, all show a clear transition to positive ⟨V_r⟩ at their outskirts. The spatial location of this transition is the turnaround radius, R_ta, depicted by a vertical red line in ⟨V_r⟩−⟨r⟩ plots.

Fig. 1. Spherically averaged radial velocity profiles of two dark matter halos, one from TDSS (left panel) with M200 = 2.15 × 1015 M⊙, and one from TNG-300-3 (right panel) with M200 = 6.63 × 1014 M⊙. Error bars indicate the ⟨Vr⟩ uncertainty (standard error of the mean) in each shell. In both panels, the red vertical line indicates the turnaround radius.

From each velocity profile, a value of R_ta is obtained as follows: starting at very large radii and moving inwards, the first crossing of ⟨V_r⟩ = 0 is located. The two points in the profile straddling that first crossing are identified; the intersection of the line defined by these two points with the r- axis is R_ta. On average, the number of particles found in the Turnaround shells in the case of TDSS is 612 and in the case of TNG-300-2 (TNG-300-3) is 2250 (275), and so the Turnaround shell is well sampled in all cases.

As discussed in the previous section, the value of R_ta does exhibit some sensitivity to the selected centre of the structure. For this reason, once a value of R_ta is estimated, the location of the center of the structure is reevaluated: the center of mass of all dark matter particles within a distance R_ta from the (original) structure center is calculated; subsequently, the spherically averaged radial velocity profile is recalculated with respect to that center of mass, and a new value of R_ta is evaluated as above. This process is repeated until the value of R_ta converges within 500 kpc, or until five iterations are performed. In the latter case, the structure is flagged.

Once we have estimated the locations of the centers of mass and R_ta, we address the substructure issue. For each structure, all neighboring structures with centers within R_ta are identified. Of these (selected neighbors and original structure), the object with the largest M₂₀₀ is labeled as “structure” and is retained. All the others are labeled as “substructure” and are removed from further consideration. Our substructure-cleaning algorithm did not find any substructures at turnaround included in the TDSS halos, and it identified and removed 70 (51) substructures from the TNG-300-2 (TNG-300-3) halo sample.

4. Results

4.1. Can a single turnaround radius meaningfully characterize a realistic 3D structure?

In an expanding universe, that is homogeneous and isotropic on large scales, the turnaround radius is always a well-defined scale, even for a structure that exhibits strong deviations from spherical symmetry in both its mass distribution and its kinematics. The collapsed center of the structure will not be expanding; eventually, sufficiently large scales will join the Hubble flow and will expand. Thus, the spherically averaged radial velocity profile will pass through zero, and the distance from the center where this occurs can always be defined as the “turnaround scale”. Whether this scale meaningfully characterizes a strongly nonspherical collapsed structure is however a different and nontrivial question. One could imagine that for a strongly nonspherical structure, particles in different directions turn around at different distances from the center. Individual particles within the “turnaround” shell could thus exhibit significant inward or outward motions in different directions that cancel out on average, a behavior similar to what is usually encountered at the centers of structures.

In order to test for such a scenario, we constructed the spherically averaged profiles of the square of the radial velocity. In a spherically symmetric structure, the turnaround radius is where all particles stand still, and the square of the radial velocity will become zero at the turnaround scale. If a single turnaround scale indeed characterizes realistic structures in our simulations, it should be similarly imprinted in the spherically averaged profiles of $V_r^2$ as a pronounced global minimum, very close to zero compared to typical squared velocities in the halo. This is indeed the case as we can see in Fig. 2. The two profiles correspond to the same structures as the ⟨V_r⟩ profiles of Fig. 1. This behavior of low or minimum $\langle V_r^2\rangle$ in the turnaround shell is typical for halos in the two simulations. We demonstrate this in Fig. 3. The upper panel shows the distribution of $\langle V_r^2\rangle$ at turnaround normalized to $V_{200}^2$ of each halo, and demonstrates that the radial velocity dispersion in the turnaround shell is small for the vast majority of halos in our sample. The lower panel shows the distribution of the ratio between the turnaround radius (measured from the linear radial velocity profiles; see Sect. 3) and the radius where $\langle V_r^2\rangle$ is minimum, demonstrating that the two are very close in all halos. We conclude that the turnaround radius obtained from the spherically averaged radial velocity profile is indeed a kinematically and dynamically meaningful scale that can describe the boundary of halos in N-body simulations despite the very strong deviations of the mass distribution of these halos from spherical symmetry.

Fig. 2. Spherically averaged profiles of the square of radial velocity, for the same objects as the corresponding panels of Fig. 1. Error bars indicate the $\langle V_r^2\rangle$ uncertainty (standard error of the mean) in each shell. In both panels, the red vertical line indicates the turnaround radius derived from the ⟨Vr⟩ profiles of Fig. 1.

Fig. 3. Upper panel: distribution of $\langle V_r^2\rangle$ in the turnaround shell in units of $V_{200}^2$ for halos in TDSS (orange histogram) and TNG-300-2 (blue histogram): the turnaround shell is indeed characterized by very small individual radial velocities of all its particles. Lower panel: distribution of the logarithm of the ratio of the turnaround radius of a structure over the radius of the shell where $\langle V_r^2\rangle$ is minimum in that same structure.

In order to visualize the origin of this result, Fig. 4 shows 2D projections of the M₂₀₀ = 6.63 × 10¹⁴ M_⊙ structure from TNG-300-3 depicted in the right-hand panels of Figs. 1 and 2. The turnaround radius is shown with a red circle, while a black circle indicates the location of R₂₀₀. The left column shows the halo column density along the three different coordinate axes, by plotting the number of dark matter particles that are projected in each bin on the plane. The mass distribution of the halo is very anisotropic on the turnaround scale. The turnaround scale itself is not clearly identifiable in the column density projections. Despite the presence of substructure, there is a very pronounced mass concentration at the center of the halo that represents a significant fraction of the total mass enclosed by the turnaround radius.

Fig. 4. Two-dimensional projections of the M200 = 6.63 × 1014 M⊙ structure from TNG-300-3 depicted in the right-hand panels of Figs. 1 and 2. Left column: number of particles (representing column density). Middle column: sum of particle radial velocities (representing radial-velocity-weighted column density). Right column: average radial velocity (representing projections of the radial velocity field). Black circles indicate the location of R200 and red circles the location of Rta. Middle and right columns: red hues correspond to outflow and blue hues to infall.

The middle and right columns on Fig. 4 represent different renderings of radial velocities. The middle column shows the radial-velocity-weighted dark matter column density by plotting the sum of particle radial velocities in each bin on the plane. The radial velocity field itself is shown in the right column, where the color in each bin on the plane corresponds to the average radial velocity of particles projected in that bin. Here the turnaround scale is immediately identifiable as a kinematic boundary of the structure, even without the guidance of the red circle. Interestingly, at the turnaround scale, the collapse or expansion is much less anisotropic than the mass distribution, and much less anisotropic than the velocity structure near the center. At the center, the collapse or bounce back of the dark matter is homogeneous only to the extent that the infalling satellite halos have been integrated to the central halo. In contrast, at the turnaround scale the dynamics appear to be dominated by the central mass concentration⁴. The gravitational potential at these large distances is then much closer to a central potential, and this is reflected in the considerably reduced anisotropy in the collapse or expansion around the kinematic boundary of the object.

4.2. How does the size of the turnaround radius of 3D structures compare to the predictions of the spherical collapse model for objects of the same enclosed mass?

We address this question in Fig. 5, where we plot the turnaround radii for the halos we have analyzed as a function of the mass enclosed within those radii. All halos were cleaned from substructure with the method described in Sect. 3. Because the simulations do not have exactly the same cosmological parameters, we rescaled all R_ta in TDSS by (ρ_TNG/ρ_TDSS)^−1/3 where ρ_TNG and ρ_TDSS are the turnaround densities predicted by spherical collapse for z = 0 for the Illustris-TNG and TDSS cosmological parameters, respectively. There is an extremely tight correlation between turnaround radius and enclosed mass. A power-law fit (black line) yields a scaling very close to $R_{\mathrm{ta}}\sim M_{\mathrm{ta}}^{1/3}$ (best-fit slope 0.353 ± 0.001), as expected from structures of comparable average density. Also, this correlation is resolution independent as halos from TNG-300-3 and TNG-300-2 have the same density behavior.

Fig. 5. Turnaround radius as a function of enclosed mass. Black points: TNG-300-2 halos; yellow points: TNG-300-3 halos; magenta points: TDSS halos. The two simulations use different cosmological parameters; TDSS data points were renormalized to scale with cosmological parameters of TNG. Blue and red lines represent theoretical predictions of the spherical collapse model for maximum (asymptotic late-time) and present-time values of turnaround radius, respectively. The black line is a power-law fit (slope: 0.353 ± 0.001).

In the same figure we overplot with the blue line the scaling with mass of the maximum (asymptotic late-time) turnaround radius (Pavlidou & Tomaras 2014), the normalization of which only depends on Λ. The red line represents the scaling predicted by the spherical collapse model for the turnaround radius of structures of all masses at the present cosmic time (see e.g., Tanoglidis et al. 2015; Pavlidou et al. 2020 for a derivation), and it is in striking agreement with the turnaround radii of simulated halos, despite the strong deviations of the latter from spherical symmetry. Halos indicated with triangles in this plot correspond to cases where the location of the center of mass of the halo failed to converge within 500 kpc after five iterations. Such cases of ambiguous halo center account for a small fraction of the outliers from the average scaling and the predictions of the spherical collapse model.

Despite the scatter that is pronounced for halos of lower virial masses, the majority of the structures lie well below the upper bound indicated by the blue line. However, there are several structures that are close to it, and it is entirely conceivable that additional inaccuracies in the observational determination of the turnaround radius and the enclosed mass of structures may result in individual measurements that are in violation of the bound. Our results are therefore consistent with the findings of Lee & Yepes (2016) that any single observation in violation of the bound is not directly at odds with ΛCDM, given observational uncertainties and realistic structure-to-structure scatter. However, unless a systematic bias is present in observations, the average behavior of the structures should accurately track the predictions of spherical collapse; if anything, simulated halos exhibit a slight bias towards smaller values of the turnaround radius compared to the expectations from spherical collapse for structures of the same enclosed mass (i.e., a bias towards a higher value of the average turnaround density; see also Sect. 4.3).

4.3. Is the average turnaround density of realistic 3D halos independent of halo mass at a given cosmic epoch?

Given that R_ta was shown in Fig. 5 to scale with enclosed mass as $M_{\mathrm{ta}}^{1/3}$, it is clear that a characteristic average density for structures defined on turnaround scales (a turnaround density, ρ_ta) does exist. The good agreement between the normalization of the R_ta − M_ta scaling and the predictions of the spherical collapse model implies that this characteristic density is indeed close to the predicted density of a single spherical structure in an otherwise unperturbed Universe turning around today: ρ_{ta, sph.coll., z = 0}.

In this section, we investigate the behavior of the distribution of ρ_ta in simulated halos with M₂₀₀ >  10¹³ M_⊙. This is important for two reasons. First, in the spherical collapse model, it is the evolution of ρ_ta with redshift that probes cosmological parameters on turnaround scales (Tanoglidis et al. 2016; Pavlidou et al., in prep.). Second, if indeed the distribution of ρ_ta is strongly peaked around the characteristic value, this could provide a straightforward way to estimate the turnaround radius of a structure based on its density profile alone, just as the virial radius, R₂₀₀, is identified in simulations and observations as the radius of a given density contrast with the mean matter density of the Universe. We stress however that, unlike R₂₀₀, which is obtained assuming that the enclosed structure has an average density of 200 times the average matter density of the Universe at the same cosmic epoch, R_ta is derived from the velocity profile of each structure, with no a priori constraint on the enclosed matter density. That such a characteristic density does appear to exist is a physical result of the dynamics of the problem rather than of our halo-finding algorithm.

To this end, in Fig. 6 we plot the distribution of the logarithm of the ratio between ρ_ta, sim measured in simulations and the value ρ_{ta, sph.coll., z = 0} prediced by spherical collapse. A value of log(ρ_ta, sim/ρ_{ta, sph.coll., z = 0}) = 0 corresponds to perfect agreement with spherical collapse predictions. Different colors correspond to different simulations and, due to the difference in the size of the simulation boxes, different halo masses (as can be also seen in Fig. 5). Both distributions are strongly peaked close to zero: TDSS halos have a median log(ρ_ta, sim/ρ_{ta, sph.coll., z = 0}) of 0.03, and a standard deviation of 0.04; TNG-300-2 halos have a median of 0.08 and a standard deviation of 0.12. The difference between halos in the two simulations is caused by the difference in the masses of the halos we have analyzed in each simulation. The larger TDSS halos are in better agreement with spherical collapse predictions. For both simulations, ρ_ta in simulations is consistent with spherical collapse predictions within one standard deviation. However, there is a systematic bias towards higher ρ_ta (lower R_ta for a given enclosed mass) in simulations; the bias decreases with increasing mass. The same trend can be seen in Fig. 5.

Fig. 6. Distribution of the logarithmic difference between average matter density within Rta and the characteristic value predicted by the spherical collapse model for all structures at z = 0 independently of their mass. The orange histogram corresponds to (higher mass) TDSS halos, while the blue histogram to (lower mass) TNG halos.

4.4. Is there a dependence of the turnaround radius on the shape of structures?

Although the correlation between R_ta and M_ta is clearly significant, has the correct slope, and has a normalization very close to the prediction of the spherical collapse model, it does feature appreciable scatter, especially for structures of lower masses, and a slight normalization offset.

We have already identified the ambiguity in the definition of the center of a turnaround structure as one of the sources of the scatter. Intuitively, one would argue that deviations of the matter distribution in a halo from spherical symmetry would also be a prime suspect as a cause for deviations of the turnaround radius from the predictions of the spherical collapse model. To quantify and evaluate this hypothesis, we construct a measure α of the asphericity of a halo on turnaround scales. To this end, we calculate the principal moments of inertia I_k with k = 1, 2, 3, and define α as

$$\begin{array}{c}\hfill \alpha =\frac{I_{k,\min }}{I_{k,\max }},\end{array}$$(1)

so that a value of 1 corresponds to a sphere, whereas a value of 0 corresponds to a prolate or oblate object of infinitesimal thickness.

In Fig. 7 we plot the logarithm of the absolute fractional deviation F ≡ log|R_ta, sim−R_{ta, sph.coll., z = 0}|/R_{ta, sph.coll., z = 0} of R_ta from the corresponding value predicted by spherical collapse for a structure of the same mass against the asphericity parameter α. For the lower-mass TNG-300-2 halos (blue dots) we observe a positive correlation, which is primarily driven by a few outlying points marked with a red triangle. These points correspond to low-mass (M₂₀₀ = 10¹⁴) halos found in environments with massive neighbors, the presence of which perturbs the potential of the central structure leading to a false identification of the turnaround radius (F >  −0.55). By performing a correlation analysis after neglecting these points, no discernible trend is seen, despite the presence of a large range of shapes and appreciable deviations from spherical collapse: that is, the dominant cause of these deviations in smaller halos is not their shape. This is confirmed by a formal correlation analysis, which yields a Spearman correlation coefficient of 0.08 and a p-value of 8.2%: the correlation is extremely weak, in the opposite direction than intuitively expected, and is not statistically significant. For the larger TDSS halos, the deviations from spherical collapse predictions are much smaller to begin with. Here some trend of decreasing deviations with decreasing asphericity (increasing α) can be seen, and a correlation analysis confirms this (p-value 8 × 10⁻⁵). However the trend is very weak (Spearman correlation coefficient −0.16).

Fig. 7. Absolute fractional deviation of Rta from the value predicted by spherical collapse for a structure of identical enclosed mass at z = 0, plotted against the asphericity parameter α. Higher α correspond to halos closer to spherical symmetry. Red triangles indicate low-mass halos with massive neighbors.

These results do not confirm the suspicion raised by Giusti & Faraoni (2019) that asphericity might induce a significant error in spherical-collapse-based estimates of the turnaround radius; rather, they are in agreement with the findings of Bhattacharya & Tomaras (2019), namely that asphericity has a very minor effect on the turnaround radius.

4.5. Is there a dependence of the turnaround radius on the presence of massive neighbors?

We have shown that asphericity is not the dominant factor driving deviations of the turnaround radius of a lower-mass simulated structure (≲10¹⁵ M_⊙ within R_ta) from the predictions of spherical collapse for a halo of identical mass. Nevertheless, we have not yet identified the culprit of the deviations in this mass range.

A hint comes from the average direction of the deviations. As seen in Fig. 6, the overall trend is towards higher values of ρ_ta. In the analytical treatment of Bhattacharya & Tomaras (2019), an increased value of the spherically averaged turnaround density (decreased value of the spherically averaged turnaround radius) is found in aspherical halos when the asphericity is driven by some effect opposing the gravitational attraction of the central potential, such as rotation. In contrast, a deviation from spherical symmetry itself, not driven by a gravity-opposing mechanism, produces to first order a zero net effect on the spherically averaged turnaround radius. It is therefore reasonable to conclude that we are looking for an effect which: (a) opposes the gravitational attraction of the central mass, and (b) is more likely to affect smaller halos.

The tidal effect of massive neighbors could play that role. It has a net effect opposing the gravity of the structure, and it is more likely to operate on small halos, as larger halos are too rare to be found near a neighbor of comparable or higher mass. To test this hypothesis, we select the neighbor with the dominant tidal effect using an environmental parameter similar to D_N, f of Haas et al. (2012); namely distance to the Nth nearest neighbor of a halo with M_ta at least f times the M_ta of the halo, normalized to R_ta of the neighbor. The tidal force due to the Nth nearest neighbor scales as $D_{N,f}^{-3}$. For simplicity, we select the neighbor within 10 × R_ta, halo from the center of the halo with the minimum value of D_N, 1 as the one that will have the dominant tidal effect on the halo, and use the D value of that neighbor as a proxy for tidal effects encountered by the central halo. Our “tidal parameter” D is therefore defined as

$$\begin{array}{c}\hfill D=min{\left\{\frac{r_{\mathrm{neighbor}-\mathrm{halo}}}{R_{\mathrm{ta},\mathrm{neighbor}}}\right|}_{M_{\mathrm{ta},\mathrm{neighbor}}\ge M_{\mathrm{ta},\mathrm{halo}}\mathrm{and}{r}_{\mathrm{neighbor}-\mathrm{halo}}\le 10\times R_{\mathrm{ta},\mathrm{halo}}}\}.\end{array}$$(2)

An advantage of D as a tidal parameter is that it is very weakly correlated with halo mass Haas et al. (2012).

In Fig. 8 we plot the logarithm of the absolute fractional deviation of R_ta from the spherical-collapse prediction against D as a scatter plot for all halos and in bins of log D. The pileup of 40 halos at D = 1 is because for the ten largest (and rarest) halos in TNG300-2, no halo of comparable mass was found within 10 × R_ta. Although we include these halos in the plot assigning to them a (lower-limit) value of D = 1, we do not include them in the correlation analysis. A significant, moderately strong correlation is indeed found, with a Spearman correlation coefficient of −0.22 (more tidally affected halos deviate more from spherical symmetry predictions) and a p-value of 10⁻⁵. We note that for log D >  0.6 the effect flattens off, likely because the tidal effects of the dominant neighbor are already very weak at that point.

Fig. 8. Fractional deviation of Rta from the value predicted by spherical collapse for a structure of identical enclosed mass at z = 0, plotted against the tidal parameter D (see text). Lower values of D correspond to a stronger tidal force. A moderately strong (Spearman coefficient −0.22) but very significant (p-value 10−5) correlation can be seen. The red data points correspond to the mean and standard error of the mean for the depicted bins in log D.

To verify that tidal effects of neighbors are indeed the dominant factor that drives the median ρ_ta to lower values in smaller halos, we show in Fig. 9 a histogram of the logarithmic difference between the turnaround density measured in TNG-300-2 halos and the predictions of spherical collapse, splitting halos according to the value of their tidal parameter, D. Indeed, halos with log D <  0.5 (which experience stronger tidal forces; blue histogram) have a higher median (0.1) than halos with log D ≥ 0.5 (median of 0.08). The latter is closer to the behavior seen in the larger TDSS halos.

Fig. 9. Histogram of the logarithmic difference between ρta measured in TNG-300-2 halos and the predictions of spherical collapse. Blue histogram: halos with strong tidal effects from their environment (log D <  0.5). Orange histogram: halos with smaller tidal effects (log D ≥ 0.5). The medians of the two distributions are clearly shifted (blue: median of 0.1; orange: median of 0.08.

5. Discussion

We focus our discussion on halos with M₂₀₀ ≥ 10¹⁴ M_⊙ (galaxy clusters and higher in mass) for two reasons. The first is practical: these are the lowest-mass objects for which an observational determination of the turnaround radius could even in principle be attempted. Objects of significantly lower mass do not include enough galaxies to allow a sufficiently accurate mapping of the Hubble flow and peculiar motions around them to determine the turnaround scale (e.g., Karachentsev & Kashibadze 2005; Nasonova et al. 2011; Karachentsev & Nasonova 2010; Lee 2018). The second is physical: most frequently halos with smaller masses tend to be found in crowded environments, and neighbors larger than themselves in close proximity are very likely. This is problematic, not only because of the tidal effects of these neighbors, but because their presence significantly contaminates the radial velocity profiles of the small halos, complicating the determination of the turnaround scale in many cases.

We explored in detail the possible effect of realistic shapes and environments on the turnaround radius of simulated dark matter halos. Similar issues were recently discussed by Lee (2016), who, extending the methodology of Falco et al. (2014), introduced empirical fits to peculiar velocity profiles as a means to probe the exterior of virialized structures. Lee (2016) investigated the behavior of peculiar velocity profiles in cosmological simulations in the case where entire halos rather than dark matter particles are used to determine the peculiar velocity profile, and applied it specifically to the case of determining the turnaround radius. Lee & Yepes (2016) used a similar methodology to study environmental effects on peculiar velocity profiles and estimate the likelihood that such effects might be responsible for observational estimates of the turnaround radius that apparently violate the bound (3GM/Λc²)^1/3 for the maximum radius of any nonexpanding structure of Pavlidou & Tomaras (2014). Our findings here are consistent with these conclusions.

More recently, Giusti & Faraoni (2019) proposed a path towards investigating analytically the effect of deviations from sphericity on the turnaround radius, and predicted that asphericities would significantly affect the turnaround size of structures. Our analysis does not confirm this expectation. The effect of asphericity on the turnaround radius is only dominant in the largest, rarest objects, and even then the effect is very weak.

In contrast, Bhattacharya & Tomaras (2019) calculated analytically that the effect of asphericities on the turnaround should be very small. When taking a spherically averaged turnaround radius, they predicted that the effect of asphericities should be nonvanishing to first order only when the asphericity is due to some gravity-opposing effect, such as rotation. Our analysis is consistent with these findings.

The turnaround radius is not the only kinematically motivated boundary of halos that can be defined. Other such boundaries proposed include the static radius (corresponding to the static mass) of Cuesta et al. (2008), and the splashback radius of Diemer & Kravtsov (2014). These are distinctly different from the turnaround radius we have studied here.

The static mass boundary corresponds to the innermost radius in which mean radial velocity is equal to zero – in other words, the innermost edge of the accretion region, within which the velocities are approximately randomized. The static part of the halo as used by Cuesta et al. (2008; and the corresponding mass and radius) is thus the kinematically defined equivalent of the virialized part of a halo. In contrast, the turnaround radius is the outermost zero crossing of the radial velocity – the outermost edge of any accretion region (which for the larger mass halos on which we focus is always present).

The splashback radius is another boundary that is meant to separate virialized from infalling material, and thus is also located at the inner edge of the accretion region. The splashback radius is defined as the radius where particles reach the apocenter of their first orbit, and corresponds to the first caustic of spherically symmetric accretion. The location of the splashback radius in simulated halos is identified through analyzing either the density field around a halo (Mansfield et al. 2017), or the distribution of individual-particle orbit apocenters (Diemer 2017). In the context of spherical collapse, it has been studied by Adhikari et al. (2014) and Shi (2016).

The turnaround radius is located at much larger distances from the halo center than both the static and splashback boundaries, and always outside the infall region. The motivation for studying it is that at these outermost nonexpanding scales, the effect of the cosmological constant or any alternative-gravity effects on individual structures (rather than the Universe as a whole) would be most pronounced (e.g., Pavlidou & Tomaras 2014; Bhattacharya et al. 2017; Lee 2018; Lopes et al. 2018; Nojiri et al. 2018; Capozziello et al. 2019; Pavlidou et al., in prep.).

6. Conclusions

Here, we use cosmological N-body simulations to probe the turnaround radius of dark matter halos with realistic shape, and compare its properties with those predicted from the spherical collapse model.

We find that a single turnaround radius can indeed accurately describe simulated halos, even in the presence of significant asphericities in their mass distribution. The value of the turnaround radius as measured from radial velocity profiles is in good agreement with spherical collapse for structures corresponding to large galaxy groups and galaxy clusters (M₂₀₀ >  10¹³ M_⊙).

Moreover, we show that deviations of the turnaround radius from the spherical collapse model are primarily driven by the tidal forces of large nearby neighbors. The effect of the deviation of halo shapes from spherical symmetry is much weaker, and only detectable in larger halos where the presence of nearby neighbors of comparable or larger mass is highly unlikely.

Our results indicate that the turnaround radius could indeed be used as an alternative boundary for studying the abundance of massive large-scale structures, as it is clean from the subtleties of baryonic effects and can be easily identifiable via the matter-density profile of structures. Indeed, although in our analysis the boundary of “turnaround structures” was identified using radial velocity profiles alone, these structures were all found to feature a characteristic average density.

This property suggests that halo-finders specifically geared towards the analysis of “turnaround” structures can be developed. For concordance cosmology, the predicted (by spherical collapse) density contrast of a region enclosed by the turnaround radius of a structure with the matter density of the background Universe at the present cosmic epoch is⁵ ∼11. Consequently, in this context, R_ta is equivalent to R₁₁ – in a way that is analogous to defining the “virial” radius as R₂₀₀. However, R₁₁ has the added benefit, shown in this work, of corresponding well to the outermost zero-crossing of the spherically averaged radial velocity profile for structures with M₂₀₀ >  10¹⁴ M_⊙.

Conversely, if the turnaround radius and mass of the halo on turnaround scales can be independently determined through observations, then the value of the turnaround density as a function of redshift can be used to probe cosmology.

Acknowledgments

G. K. would like to thank A. Zezas for his constructive comments throughout this work and G. Mouloudakis and E. Garaldi for fruitful discussions. We thank A. Cuesta, B. Diemer, A. Giusti, A. Kravtsov, and F. Prada for feedback and on our manuscript and helpful suggestions, and an anonymous referee for constructive comments that have improved this work. G. K. and K. T. acknowledge support from the European Research Council under the European Union’s Horizon 2020 research and innovation program, under grant agreement no 771282. E. N. as a MC Fellow is supported by a Marie Curie Action of the European Union (grant agreement number 749073). K. K. acknowledges funding from the European Research Council under the European Union’s Seventh Framework Program (FP/2007-2013)/ERC grant agreement no 617001. This project has received funding from the European Union’s Horizon 2020 research and innovation program under the Marie Sklodowska-Curie RISE action, grant agreement no 691164 (ASTROSTAT).

References

All Tables

All Figures

	Fig. 1. Spherically averaged radial velocity profiles of two dark matter halos, one from TDSS (left panel) with M200 = 2.15 × 1015 M⊙, and one from TNG-300-3 (right panel) with M200 = 6.63 × 1014 M⊙. Error bars indicate the ⟨Vr⟩ uncertainty (standard error of the mean) in each shell. In both panels, the red vertical line indicates the turnaround radius.
In the text

	Fig. 2. Spherically averaged profiles of the square of radial velocity, for the same objects as the corresponding panels of Fig. 1. Error bars indicate the $\langle V_r^2\rangle$ uncertainty (standard error of the mean) in each shell. In both panels, the red vertical line indicates the turnaround radius derived from the ⟨Vr⟩ profiles of Fig. 1.
In the text

Fig. 3. Upper panel: distribution of $\langle V_r^2\rangle$ in the turnaround shell in units of $V_{200}^2$ for halos in TDSS (orange histogram) and TNG-300-2 (blue histogram): the turnaround shell is indeed characterized by very small individual radial velocities of all its particles. Lower panel: distribution of the logarithm of the ratio of the turnaround radius of a structure over the radius of the shell where $\langle V_r^2\rangle$ is minimum in that same structure.

In the text

Fig. 4. Two-dimensional projections of the M200 = 6.63 × 1014 M⊙ structure from TNG-300-3 depicted in the right-hand panels of Figs. 1 and 2. Left column: number of particles (representing column density). Middle column: sum of particle radial velocities (representing radial-velocity-weighted column density). Right column: average radial velocity (representing projections of the radial velocity field). Black circles indicate the location of R200 and red circles the location of Rta. Middle and right columns: red hues correspond to outflow and blue hues to infall.

In the text

Fig. 5. Turnaround radius as a function of enclosed mass. Black points: TNG-300-2 halos; yellow points: TNG-300-3 halos; magenta points: TDSS halos. The two simulations use different cosmological parameters; TDSS data points were renormalized to scale with cosmological parameters of TNG. Blue and red lines represent theoretical predictions of the spherical collapse model for maximum (asymptotic late-time) and present-time values of turnaround radius, respectively. The black line is a power-law fit (slope: 0.353 ± 0.001).

In the text

	Fig. 6. Distribution of the logarithmic difference between average matter density within Rta and the characteristic value predicted by the spherical collapse model for all structures at z = 0 independently of their mass. The orange histogram corresponds to (higher mass) TDSS halos, while the blue histogram to (lower mass) TNG halos.
In the text

	Fig. 7. Absolute fractional deviation of Rta from the value predicted by spherical collapse for a structure of identical enclosed mass at z = 0, plotted against the asphericity parameter α. Higher α correspond to halos closer to spherical symmetry. Red triangles indicate low-mass halos with massive neighbors.
In the text

Fig. 8. Fractional deviation of Rta from the value predicted by spherical collapse for a structure of identical enclosed mass at z = 0, plotted against the tidal parameter D (see text). Lower values of D correspond to a stronger tidal force. A moderately strong (Spearman coefficient −0.22) but very significant (p-value 10−5) correlation can be seen. The red data points correspond to the mean and standard error of the mean for the depicted bins in log D.

In the text

	Fig. 9. Histogram of the logarithmic difference between ρta measured in TNG-300-2 halos and the predictions of spherical collapse. Blue histogram: halos with strong tidal effects from their environment (log D <  0.5). Orange histogram: halos with smaller tidal effects (log D ≥ 0.5). The medians of the two distributions are clearly shifted (blue: median of 0.1; orange: median of 0.08.
In the text