Issue 
A&A
Volume 639, July 2020



Article Number  A122  
Number of page(s)  9  
Section  Cosmology (including clusters of galaxies)  
DOI  https://doi.org/10.1051/00046361/201937337  
Published online  20 July 2020 
Turnaround radius of galaxy clusters in Nbody simulations
^{1}
Department of Physics and Institute for Theoretical and Computational Physics, University of Crete, 70013 Heraklio, Greece
email: gkorkidis@physics.uoc.gr, pavlidou@physics.uoc.gr
^{2}
Institute of Astrophysics, Foundation for Research and Technology – Hellas, Vassilika Vouton, 70013 Heraklio, Greece
^{3}
Scuola Normale Superiore, Piazza dei Cavalieri 7, 56126 Pisa, Italy
Received:
17
December
2019
Accepted:
27
April
2020
Aims. We use Nbody simulations to examine whether a characteristic turnaround radius, as predicted from the spherical collapse model in a ΛCDM Universe, can be meaningfully identified for galaxy clusters in the presence of full threedimensional effects.
Methods. We use The Dark Sky Simulations and IllustrisTNG darkmatteronly cosmological runs to calculate radial velocity profiles around collapsed structures, extending out to many times the virial radius R_{200}. There, the turnaround radius can be unambiguously identified as the largest nonexpanding scale around a center of gravity.
Results. We find that: (a) a single turnaround scale can meaningfully describe strongly nonspherical structures. (b) For halos of masses M_{200} > 10^{13} M_{⊙}, the turnaround radius R_{ta} scales with the enclosed mass M_{ta} as M_{ta}^{1/3}, as predicted by the spherical collapse model. (c) The deviation of R_{ta} in simulated halos from the spherical collapse model prediction is relatively insensitive to halo asphericity. Rather, it is sensitive to the tidal forces due to massive neighbors when these are present. (d) Halos exhibit a characteristic average density within the turnaround scale. This characteristic density is dependent on cosmology and redshift. For the present cosmic epoch and for concordance cosmological parameters (Ω_{m} ∼ 0.3; Ω_{Λ} ∼ 0.7) turnaround structures exhibit a density contrast with the matter density of the background Universe of δ ∼ 11. Thus, R_{ta} is equivalent to R_{11} – in a way that is analogous to defining the “virial” radius as R_{200} – with the advantage that R_{11} is shown in this work to correspond to a kinematically relevant scale in Nbody simulations.
Key words: largescale structure of Universe / methods: numerical / galaxies: clusters: general
© 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 welldefined boundary for largescale 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^{2}/8πG. (d) As a consequence, the radius of any nonexpanding structure of mass M in a ΛCDM universe can never exceed (3GM/Λc^{2})^{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 Nbody 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 Nbody 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. Nbody 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 treebased Nbody code. The data extraction and analysis were carried out using yt (Turk et al. 2011; b) The IllustrisTNG simulations (Nelson et al. 2015, 2019), which is a suite of gravomagnetohydrodynamical simulations run with the AREPO moving mesh code (Springel 2010). For the data extraction of the IllustrisTNG we used the methods presented on the project website^{1} 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.
Simulation number of particles, box size, particle mass, and Plummer equivalent gravitational softening.
Cosmological parameters used in simulations.
For producing halo catalogs, TDSS uses the ROCKSTAR algorithm (Behroozi et al. 2013) and IllustrisTNG uses friendsoffriends (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 sixdimensional 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 (TNG3002 and TNG3003). 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 halofinding 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 halofinding 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 halofinding 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 TNG300. More specifically, from TDSS we randomly chose 578 halos with M_{200} ≥ 10^{15} M_{⊙}^{2} while from TNG3002 (TNG3003) we chose 430 (449) halos with masses in the interval 10^{14} M_{⊙} ≤ M_{200} ≤ 10^{15} 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^{3}. 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 M_{200} = 2.15 × 10^{15} M_{⊙}, and one from TNG3003 (right panel) with M_{200} = 6.63 × 10^{14} M_{⊙}. Error bars indicate the ⟨V_{r}⟩ 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 TNG3002 (TNG3003) 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_{200} is labeled as “structure” and is retained. All the others are labeled as “substructure” and are removed from further consideration. Our substructurecleaning algorithm did not find any substructures at turnaround included in the TDSS halos, and it identified and removed 70 (51) substructures from the TNG3002 (TNG3003) 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 welldefined 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 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 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 at turnaround normalized to 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 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 Nbody 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 uncertainty (standard error of the mean) in each shell. In both panels, the red vertical line indicates the turnaround radius derived from the ⟨V_{r}⟩ profiles of Fig. 1. 
Fig. 3. Upper panel: distribution of in the turnaround shell in units of for halos in TDSS (orange histogram) and TNG3002 (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 is minimum in that same structure. 
In order to visualize the origin of this result, Fig. 4 shows 2D projections of the M_{200} = 6.63 × 10^{14} M_{⊙} structure from TNG3003 depicted in the righthand panels of Figs. 1 and 2. The turnaround radius is shown with a red circle, while a black circle indicates the location of R_{200}. 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. Twodimensional projections of the M_{200} = 6.63 × 10^{14} M_{⊙} structure from TNG3003 depicted in the righthand panels of Figs. 1 and 2. Left column: number of particles (representing column density). Middle column: sum of particle radial velocities (representing radialvelocityweighted column density). Right column: average radial velocity (representing projections of the radial velocity field). Black circles indicate the location of R_{200} and red circles the location of R_{ta}. 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 radialvelocityweighted 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^{4}. 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 IllustrisTNG and TDSS cosmological parameters, respectively. There is an extremely tight correlation between turnaround radius and enclosed mass. A powerlaw fit (black line) yields a scaling very close to (bestfit slope 0.353 ± 0.001), as expected from structures of comparable average density. Also, this correlation is resolution independent as halos from TNG3003 and TNG3002 have the same density behavior.
Fig. 5. Turnaround radius as a function of enclosed mass. Black points: TNG3002 halos; yellow points: TNG3003 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 latetime) and presenttime values of turnaround radius, respectively. The black line is a powerlaw fit (slope: 0.353 ± 0.001). 
In the same figure we overplot with the blue line the scaling with mass of the maximum (asymptotic latetime) 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 structuretostructure 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 , 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_{200} > 10^{13} 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_{200}, 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_{200}, 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 halofinding 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; TNG3002 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 R_{ta} 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
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 ≡ logR_{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 lowermass TNG3002 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 lowmass (M_{200} = 10^{14}) 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 pvalue 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 (pvalue 8 × 10^{−5}). However the trend is very weak (Spearman correlation coefficient −0.16).
Fig. 7. Absolute fractional deviation of R_{ta} 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 lowmass halos with massive neighbors. 
These results do not confirm the suspicion raised by Giusti & Faraoni (2019) that asphericity might induce a significant error in sphericalcollapsebased 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 lowermass simulated structure (≲10^{15} 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 gravityopposing 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 . 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
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 sphericalcollapse 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 TNG3002, no halo of comparable mass was found within 10 × R_{ta}. Although we include these halos in the plot assigning to them a (lowerlimit) 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 pvalue of 10^{−5}. 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 R_{ta} 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 (pvalue 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 TNG3002 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 TNG3002 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_{200} ≥ 10^{14} M_{⊙} (galaxy clusters and higher in mass) for two reasons. The first is practical: these are the lowestmass 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^{2})^{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 gravityopposing 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 individualparticle 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 alternativegravity 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 Nbody 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_{200} > 10^{13} 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 largescale structures, as it is clean from the subtleties of baryonic effects and can be easily identifiable via the matterdensity 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 halofinders 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^{5} ∼11. Consequently, in this context, R_{ta} is equivalent to R_{11} – in a way that is analogous to defining the “virial” radius as R_{200}. However, R_{11} has the added benefit, shown in this work, of corresponding well to the outermost zerocrossing of the spherically averaged radial velocity profile for structures with M_{200} > 10^{14} 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.
See however Diemer (2017) for an algorithm for the kinematic identification of the splashback radius. The “static mass” boundary of Cuesta et al. (2008) is also kinematically identified. See Sect. 5 for a detailed discussion of the differences between these different kinematicallyrelevant halo boundaries.
The only exception would be the presence of a massive neighbor that could create anisotropies to the potential. This could in turn create anisotropies in the velocity field at the turnaround scale, as can be seen in the yz velocity projection where a structure is outflowing in contrast to the bulk of structures around it.
We note that this value depends on the cosmological parameters, and in general evolves with cosmic epoch. See Pavlidou et al. (2020) for a derivation of its evolution with time for any value of Ω_{Λ} and Ω_{m}.
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/20072013)/ERC grant agreement no 617001. This project has received funding from the European Union’s Horizon 2020 research and innovation program under the Marie SklodowskaCurie RISE action, grant agreement no 691164 (ASTROSTAT).
References
 Adhikari, S., Dalal, N., & Chamberlain, R. T. 2014, JCAP, 2014, 019 [Google Scholar]
 Behroozi, P. S., Wechsler, R. H., & Wu, H.Y. 2013, ApJ, 762, 109 [NASA ADS] [CrossRef] [Google Scholar]
 Bertschinger, E. 1985, ApJS, 58, 39 [NASA ADS] [CrossRef] [Google Scholar]
 Bhattacharya, S., & Tomaras, T. N. 2017, Eur. Phys. J. C, 77, 526 [CrossRef] [EDP Sciences] [Google Scholar]
 Bhattacharya, S., & Tomaras, T. N. 2019, ArXiv eprints [arXiv:1911.06228] [Google Scholar]
 Bhattacharya, S., Dialektopoulos, K. F., Enea Romano, A., Skordis, C., & Tomaras, T. N. 2017, JCAP, 2017, 018 [CrossRef] [Google Scholar]
 Capozziello, S., Dialektopoulos, K. F., & Luongo, O. 2019, Int. J. Mod. Phys. D, 28, 1950058 [CrossRef] [Google Scholar]
 Cuesta, A. J., Prada, F., Klypin, A., & Moles, M. 2008, MNRAS, 389, 385 [NASA ADS] [CrossRef] [Google Scholar]
 Davis, M., Efstathiou, G., Frenk, C. S., & White, S. D. M. 1985, ApJ, 292, 371 [NASA ADS] [CrossRef] [Google Scholar]
 Diemer, B. 2017, ApJS, 231, 5 [CrossRef] [Google Scholar]
 Diemer, B., & Kravtsov, A. V. 2014, ApJ, 789, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Falco, M., Hansen, S. H., Wojtak, R., et al. 2014, MNRAS, 442, 1887 [NASA ADS] [CrossRef] [Google Scholar]
 Fillmore, J. A., & Goldreich, P. 1984, ApJ, 281, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Giusti, A., & Faraoni, V. 2019, Phys. Dark Univ., 26, 100353 [CrossRef] [Google Scholar]
 Gunn, J. E., & Gott, J. R., III 1972, ApJ, 176, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Haas, M. R., Schaye, J., & JeesonDaniel, A. 2012, MNRAS, 419, 2133 [NASA ADS] [CrossRef] [Google Scholar]
 Hoffmann, K., Planelles, S., Gaztañaga, E., et al. 2014, MNRAS, 442, 1197 [CrossRef] [Google Scholar]
 Karachentsev, I. D., & Kashibadze, O. G. 2005, ArXiv eprints [arXiv:astroph/0509207] [Google Scholar]
 Karachentsev, I. D., & Nasonova, O. G. 2010, MNRAS, 405, 1075 [NASA ADS] [Google Scholar]
 Knebe, A., Knollmann, S. R., Muldrew, S. I., et al. 2011, MNRAS, 415, 2293 [NASA ADS] [CrossRef] [Google Scholar]
 Kravtsov, A. V., & Borgani, S. 2012, ARA&A, 50, 353 [NASA ADS] [CrossRef] [Google Scholar]
 Lee, J. 2016, ApJ, 832, 123 [CrossRef] [Google Scholar]
 Lee, J. 2018, ApJ, 856, 57 [CrossRef] [Google Scholar]
 Lee, J., & Yepes, G. 2016, ApJ, 832, 185 [CrossRef] [Google Scholar]
 Lopes, R. C. C., Voivodic, R., Abramo, L. R., & Sodré, L., Jr. 2018, JCAP, 2018, 010 [NASA ADS] [CrossRef] [Google Scholar]
 Lopes, R. C. C., Voivodic, R., Abramo, L. R., & Sodré, L., Jr. 2019, JCAP, 2019, 026 [CrossRef] [Google Scholar]
 Mansfield, P., Kravtsov, A. V., & Diemer, B. 2017, ApJ, 841, 34 [CrossRef] [Google Scholar]
 More, S., Diemer, B., & Kravtsov, A. V. 2015, ApJ, 810, 36 [NASA ADS] [CrossRef] [Google Scholar]
 Nasonova, O. G., de Freitas Pacheco, J. A., & Karachentsev, I. D. 2011, A&A, 532, A104 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Nelson, D., Pillepich, A., Genel, S., et al. 2015, Astron. Comput., 13, 12 [NASA ADS] [CrossRef] [Google Scholar]
 Nelson, D., Springel, V., Pillepich, A., et al. 2019, Comput. Astrophys. Cosmol., 6, 2 [Google Scholar]
 Nojiri, S., Odintsov, S. D., & Faraoni, V. 2018, Phys. Rev. D, 98, 024005 [CrossRef] [Google Scholar]
 Pavlidou, V., & Tomaras, T. N. 2014, JCAP, 9, 020 [NASA ADS] [CrossRef] [Google Scholar]
 Pavlidou, V., Korkidis, G., Tomaras, T. N., & Tanoglidis, D. 2020, A&A, 638, L8 [CrossRef] [EDP Sciences] [Google Scholar]
 Shi, X. 2016, MNRAS, 459, 3711 [NASA ADS] [CrossRef] [Google Scholar]
 Skillman, S. W., Warren, M. S., Turk, M. J., et al. 2014, ArXiv eprints [arXiv:1407.2600] [Google Scholar]
 Springel, V. 2010, MNRAS, 401, 791 [NASA ADS] [CrossRef] [Google Scholar]
 Tanoglidis, D., Pavlidou, V., & Tomaras, T. N. 2015, JCAP, 12, 060 [CrossRef] [Google Scholar]
 Tanoglidis, D., Pavlidou, V., & Tomaras, T. 2016, ArXiv eprints [arXiv:1601.03740] [Google Scholar]
 Turk, M. J., Smith, B. D., Oishi, J. S., et al. 2011, ApJS, 192, 9 [NASA ADS] [CrossRef] [Google Scholar]
 Warren, M. S. 2013, ArXiv eprints [arXiv:1310.4502] [Google Scholar]
All Tables
Simulation number of particles, box size, particle mass, and Plummer equivalent gravitational softening.
All Figures
Fig. 1. Spherically averaged radial velocity profiles of two dark matter halos, one from TDSS (left panel) with M_{200} = 2.15 × 10^{15} M_{⊙}, and one from TNG3003 (right panel) with M_{200} = 6.63 × 10^{14} M_{⊙}. Error bars indicate the ⟨V_{r}⟩ 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 uncertainty (standard error of the mean) in each shell. In both panels, the red vertical line indicates the turnaround radius derived from the ⟨V_{r}⟩ profiles of Fig. 1. 

In the text 
Fig. 3. Upper panel: distribution of in the turnaround shell in units of for halos in TDSS (orange histogram) and TNG3002 (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 is minimum in that same structure. 

In the text 
Fig. 4. Twodimensional projections of the M_{200} = 6.63 × 10^{14} M_{⊙} structure from TNG3003 depicted in the righthand panels of Figs. 1 and 2. Left column: number of particles (representing column density). Middle column: sum of particle radial velocities (representing radialvelocityweighted column density). Right column: average radial velocity (representing projections of the radial velocity field). Black circles indicate the location of R_{200} and red circles the location of R_{ta}. 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: TNG3002 halos; yellow points: TNG3003 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 latetime) and presenttime values of turnaround radius, respectively. The black line is a powerlaw fit (slope: 0.353 ± 0.001). 

In the text 
Fig. 6. Distribution of the logarithmic difference between average matter density within R_{ta} 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 R_{ta} 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 lowmass halos with massive neighbors. 

In the text 
Fig. 8. Fractional deviation of R_{ta} 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 (pvalue 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 TNG3002 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 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.