The Effect of Modified Gravity on the Odds of the Bound Violations of the Turn-around Radii

Consider a massive group for which prior information on the virial mass and radius (M_vir and r_vir, respectively) is available. Its gravity will influence the peculiar motions not only of its satellites located in the infall zone, but also of the neighbor galaxies located in the bound zone that corresponds to the distance range of $(3\mbox{--}8){r}_{\mathrm{vir}}$. For the case of the GR+ΛCDM cosmology, the interplay between the gravity and the expanding space molds the bound-zone velocity profile around a massive group to have the following power-law shape (Falco et al. 2014):

$$\begin{eqnarray}&&\displaystyle \frac{v(r)}{{V}_{c}}=-A{\left(\displaystyle \frac{{r}_{\mathrm{vir}}}{r}\right)}^{-n},\end{eqnarray} \tag{ 2 }$$

where A and n are the amplitude and the slope parameters of the profile, respectively, and ${V}_{c}\equiv {({{GM}}_{\mathrm{vir}}/{r}_{\mathrm{vir}})}^{1/2}$. After Falco et al. (2014) reported that Equation (2) with $A\approx 0.8$ and $n\approx 0.42$ fitted well the average profile v(r) numerically obtained from the bound-zone DM particles around the group-size halos with ${M}_{\mathrm{vir}}\sim {10}^{14}\,{h}^{-1}\,{M}_{\odot }$ at z = 0. It was proven that with the same best-fit parameters, Equation (2) still describe well the average bound-zone velocity profiles obtained not from the DM particles but from the bound DM halos (Lee 2016).

In the DGP model, the presence of MG is more eminent in the bound zone than in the infall zone since the former has lower densities. In other words, in the DGP model the bound-zone objects around a massive group react more sensitively to a fifth force produced by the unscreened MG than the infall-zone satellites. We speculate that the fifth force would decrease the slope of the average bound-zone velocity profile since the effective gravity of the massive group enhanced by the fifth force is capable of resisting the expanding space at farther distances. Before quantitatively investigating if the decrement of the slope of the bound-zone velocity profile will be substantial, however, it is first necessary to confirm that the bound-zone peculiar velocity profiles for the nDGP models can be also described by the same formula, Equation (2), whose validity was only tested for the case of the ΛCDM cosmology.

Putting a mass threshold cut ${M}_{\mathrm{vir},\mathrm{th}}={10}^{13}\,{h}^{-1}\,{M}_{\odot }$ on the Rockstar halo catalog described in Section 2, we make a sample of the central groups with virial masses ${M}_{\mathrm{vir}}\geqslant {M}_{\mathrm{vir},\mathrm{th}}$. For each central group in the sample, we look for the neighbor Rockstar halos that satisfy two conditions. First, they should belong to the bound zone around the central group with their separation distances r lie in the range of $3\leqslant r/{r}_{\mathrm{vir}}\leqslant 8$. Second, the numbers of the particles, ${N}_{{\rm{p}}}$, that comprise a halo are equal to or larger than 20. The bound-zone halos composed of less than 20 dark matter particles are excluded to avoid possible contamination caused by incomplete condensation.

Let ${{\boldsymbol{v}}}_{G}$ and ${{\boldsymbol{v}}}_{b}$ denote the comoving velocities of a central group and a bound-zone halo, respectively, and let ${\boldsymbol{r}}$ be the separation vector from the central group to the bound-zone halo. We first subtract ${{\boldsymbol{v}}}_{G}$ from ${{\boldsymbol{v}}}_{b}$ to obtain the relative peculiar velocity of the bound-zone halo in the rest frame of the central group. Then, we perform the dot product between ${{\boldsymbol{v}}}_{b}-{{\boldsymbol{v}}}_{G}$ and $\hat{{\boldsymbol{r}}}\equiv {\boldsymbol{r}}/| {\boldsymbol{r}}|$ to project the relative peculiar velocity of the bound-zone halo onto the radial direction, $\hat{{\boldsymbol{r}}}$. Let v denote the magnitude of the projected relative peculiar velocity, $| ({{\boldsymbol{v}}}_{b}-{{\boldsymbol{v}}}_{G})\cdot {\boldsymbol{h}}|$, and call it the bound-zone velocity at the separation distance $r\equiv | {\boldsymbol{r}}|$.

Dividing v and r of each bound-zone halo by the circular velocity of its central group V_c and the virial radius r_vir, respectively, we express the rescaled bound-zone velocity profile $\tilde{v}(\tilde{r})\equiv v/{V}_{{\rm{c}}}$ as a function of the rescaled separation distance $\tilde{r}\equiv r/{r}_{\mathrm{vir}}$. Note that both of $\tilde{v}$ and $\tilde{r}$ are dimensionless. We also divide the range of $\tilde{r}$ into several intervals, $[\tilde{r},\tilde{r}+{\rm{\Delta }}\tilde{r}]$, each of which has the same length ${\rm{\Delta }}\tilde{r}$. We record the numbers of the bound-zone halos, N_b, whose values of the rescaled distances $\tilde{r}$ belong to each interval. The mean bound-zone velocity at each $\tilde{r}$-interval is computed by taking the average over those N_b halos.

Let ${N}_{b,k}$ and ${N}_{T,k}$ denote the numbers of the bound-zone halos around each central group and the number of the central groups, respectively, in the kth realization of each model. Let ${\tilde{v}}_{{ij},k}(\tilde{r})$ also denote the bound-zone velocity of the ith halo whose separation distance lies in the range of $[\tilde{r},\tilde{r}+{\rm{\Delta }}\tilde{r}]$ around the jth central group in the kth realization of each model. The bound-zone velocity profile in the kth realization, ${\tilde{v}}_{k}(\tilde{r})$, can be computed by taking the average first over the ${N}_{{\rm{b}},{\rm{k}}}$ bound-zone halos and then over the ${N}_{{\rm{T}},{\rm{k}}}$ central groups as

$$\begin{eqnarray}&&{\tilde{v}}_{k}(\tilde{r})=\displaystyle \frac{1}{{N}_{{\rm{T}},{\rm{k}}}}\displaystyle \sum _{j=1}^{{N}_{{\rm{T}}}}\left(\displaystyle \frac{1}{{N}_{{\rm{b}},{\rm{k}}}}\displaystyle \sum _{i=1}^{{N}_{{\rm{b}}}}{\tilde{v}}_{{ij},k}(\tilde{r})\right).\end{eqnarray} \tag{ 3 }$$

Finally, the bound-zone velocity profile $\tilde{v}(\tilde{r})$ for each model can be obtained by taking the average over the five realizations. The errors can be also computed using the bootstrap technique. Let ${\{{\tilde{v}}_{k}(\tilde{r})\}}_{k=1}^{5}$ denote the original sample of the bound-zone velocity profiles from the five realizations. From this sample, we draw five bound-zone velocity profiles with repetition allowed to create a boostrap resample. We create 1000 Boostrap resamples and calculate the errors associated with the measurement of the average bound-zone velocity profile as the one standard deviation scatter around the mean value:

$$\begin{eqnarray}&&{\sigma }_{v}^{2}(\tilde{r})=\displaystyle \frac{1}{1000}\displaystyle \sum _{\alpha =1}^{1000}{\left[\left(\displaystyle \frac{1}{5}\displaystyle \sum _{k=1}^{5}{\tilde{v}}_{k}^{\alpha }(\tilde{r})\right)-\left(\displaystyle \frac{1}{5}\displaystyle \sum _{k=1}^{5}{\tilde{v}}_{k}(\tilde{r})\right)\right]}^{2},\end{eqnarray} \tag{ 4 }$$

where ${\tilde{v}}_{k}(\tilde{r})$ represents the average bound-zone velocity profile from the original sample of the kth realization and ${\tilde{v}}_{k}^{\alpha }(\tilde{r})$ is from the αth Bootstrap resample of the kth realization.

Figure 1 plots the bound-zone peculiar velocity profile averaged over five realizations as filled circles with the Boostrap errors. There is a substantial difference in $\tilde{v}$ between the GR and the DF5 models. The latter has a lower slope than the former, as speculated. Furthermore, the magnitude of $\tilde{v}$ in the DF5 model is larger in the whole bound-zone range than in the GR model. Fitting the numerically obtained profile $\tilde{v}(\tilde{r})$ of each model to Equation (2) by employing the maximum-likelihood method, we search for the values of n and A which minimize the following ${\chi }^{2}$:

$$\begin{eqnarray}&&{\chi }^{2}=\displaystyle \sum _{i=1}^{{N}_{{\rm{r}}}}\displaystyle \frac{{[\tilde{v}({\tilde{r}}_{i})-{\tilde{v}}^{\mathrm{the}}({\tilde{r}}_{i}| {n}_{{\rm{v}}},\beta )]}^{2}}{{\sigma }^{2}({\tilde{r}}_{i})},\end{eqnarray} \tag{ 5 }$$

where ${\tilde{r}}_{i}$ denotes the ith interval of $\tilde{r}$ and ${\tilde{v}}^{\mathrm{the}}$ represents the theoretical prediction of Equation (2). To determine the uncertainties associated with the determination of n and A, we first determine the joint probability $p(A,n)\,=p[-{\chi }^{2}(a,b)/2]$. Using the probability density functions p(n) and p(A) calculated as $p(A)=\int {dn}\,p(A,n)$ and $p(n)\,=\int {dA}\,p(A,n)$, respectively, we also determine the marginalized errors, ${\sigma }_{n}$ and ${\sigma }_{A}$.

Zoom In Zoom Out Reset image size

Figure 2 shows the boundaries of three different regions in the space spanned by A and n that enclose those points, over which the integration of $P(A,n)$ becomes 0.68, 0.95, and 0.99 as the thickest, thick, and thinnest lines, respectively, for each model. The third and fourth columns of Table 1 display the best-fit values of A and n with the associated errors (${\sigma }_{A}$ and ${\sigma }_{n})$ for each model. The solid lines in Figure 1 correspond to the analytic model whose characteristic parameters are set at the best-fit values listed in Table 1. As shown, the power-law slope n has the largest (smallest) value for the GR (DF5) case and the difference in the value of n between the two models is statistically significant. Due to the non-vanishing fifth force, the bound-zone velocity profile in the DF5 model decreases less rapidly with the distance than in the GR model. For the case of the DF6 model, the fifth force is not strong enough to produce any statistically significant difference in $\tilde{v}(\tilde{r})$ from the GR case, as expected.

Zoom In Zoom Out Reset image size

The results shown in Figure 1 have one important implication. The peculiar velocity profiles for the nDGP models are still well fitted by Equation (2), which was empirically derived by Falco et al. (2014) from N-body simulations for the GR+ΛCDM cosmology. Having no analytic framework within which the bound-zone velocity profile can be derived from the first principles for nDGP models, we have assumed that no matter what background cosmology is used, the peculiar velocity of a bound-zone galaxy may be proportional to some power of the separation distance and thus that the same functional form of Equation (2) can still describe the average peculiar velocity profile even for the nDGP models. This assumption is justified by the very fact that the current work has found a good agreement in Equation (2) with the best-fit parameters and the numerical results obtained from the N-body simulations for the two nDGP models.

Before estimating the turn-around radii by using $\tilde{v}(\tilde{r})$, it may be worth addressing one crucial issue. The Rockstar algorithm counts only the bound particles to calculate the virial mass of a halo, which adopts the same definition of the boundedness regardless of the background cosmology. However, any departure of the gravitational law from GR may change the concept of being gravitationally bound, the mass of a halo computed by the Rockstar finder may not be the true virial mass, which in turn may change the shape of the peculiar velocity profile in the nDGP models. To address this issue, we include the unbound particles within the virial radius of each central group to compute its virial mass. Then, repeating the same procedure, we redetermine $\tilde{v}(\tilde{r})$ for the three models, which are plotted in Figure 3. As shown, there is almost no change between the results displayed in Figures 1 and 3, which implies that the difference in the definition of the boundness between GR and nDGP models is unlikely to have a significant effect on our estimates of the turn-around radii of the central groups from the average bound-zone peculiar velocity profiles.

Zoom In Zoom Out Reset image size

The GR+ΛCDM cosmology puts an upper bound, ${r}_{t,u}$, on the average turn-around radius of a galaxy group with mass M_vir (Pavlidou & Tomaras 2014):

$$\begin{eqnarray}&&{r}_{t,u}={\left(\displaystyle \frac{3{M}_{\mathrm{vir}}G}{{\rm{\Lambda }}{C}^{2}}\right)}^{1/3}.\end{eqnarray} \tag{ 6 }$$

This upper bound, however, limits the expectation value of the turn-around radius but not the individual values of r_t because the event of a turning around is a generically random process (Pavlidou & Tomaras 2014). In other words, the turn-around radii of individual galaxy groups can have values larger than ${r}_{t,u}$, although the occurrences of such bound violations are quite rare (Lee & Yepes 2016).

Strictly speaking, the turn-around radius r_t of a galaxy group is an attribute that it acquires at the end of its proto-group stage. The optimal way to estimate the turn-around radius of a central group in a N-body experiment is to track the trajectories of the component DM particles back to the proto-group regime and then to find the location at which the mean peculiar velocity equals the Hubble speed. This optimal routine, however, is simply impractical and thus not applicable to real data. Recently, Lee et al. (2015) formulated a less optimal but much more practical routine that makes it possible to estimate the turn-around radius of a galaxy group from the direct observables. This routine counts on Equation (2) to find r_t at which the following equation holds true:

$$\begin{eqnarray}&&\displaystyle \frac{{H}_{0}{r}_{t}({M}_{\mathrm{vir}})}{{V}_{c}}=A{\left[\displaystyle \frac{{r}_{t}({M}_{\mathrm{vir}})}{{r}_{\mathrm{vir}}}\right]}^{-n}.\end{eqnarray} \tag{ 7 }$$

In Section 3.1 we determine the best-fit values of A and n for the mean bound-zone velocity profile averaged over the central groups. To use Equation (7) to estimate r_t of each central group, however, it is necessary to determine the values of A and n by separately fitting the individual bound-zone velocity profile around each central group to Equation (2). Let ${\tilde{v}}_{j,k}(\tilde{r})$ be the bound-zone velocity profile of the jth central group in the kth realization of each model. Replacing ${\tilde{v}}_{k}$ in Equation (5) by ${\tilde{v}}_{j,k}(\tilde{r})$ and minimizing ${\chi }^{2}$, we determine the best-fit values of A and n and put them into Equation (7) to estimate the turn-around radius of the jth central group in the kth realization, say ${r}_{t,j,k}({M}_{\mathrm{vir}})$.

Dividing the range of the logarithmic masses of the central groups, ${m}_{\mathrm{vir}}\equiv \mathrm{log}({M}_{\mathrm{vir}}h\,{M}_{\odot }^{-1})$, into several intervals each of which has the same length ${\rm{\Delta }}{m}_{\mathrm{vir}}$, we calculate the mean turn-around radius, ${r}_{t,k}({M}_{\mathrm{vir}})$, by taking the average of ${r}_{t,j,k}({M}_{\mathrm{vir}})$ over the central groups whose logarithmic masses lie in a given interval of $[{m}_{\mathrm{vir}},{m}_{\mathrm{vir}}+{\rm{\Delta }}{m}_{\mathrm{vir}}]$. The mean turn-around radius, r_t, of a central group at each m_vir-interval is now evaluated as taking the average of ${r}_{t,k}$ over the five realizations for each model. Figure 4 plots r_t versus m_vir for the GR, DF6, and DF5 models as red, blue, and green solid lines, respectively. As shown, there is a notable difference in ${r}_{t}({m}_{\mathrm{vir}})$ between the GR and the DF5 models. To see whether or not this difference is statistically significant, we calculate the Bootstrap errors ${\sigma }_{r}$ in the estimation of r_t by generating 1000 Bootstrap resamples as done in Section 3.1.

Zoom In Zoom Out Reset image size

Figure 5 shows ${r}_{t}({m}_{\mathrm{vir}})$ as filled circles with the Bootstrap errors for the GR case. The red solid line represents the upper bound limit ${r}_{t,u}$ given in Equation (6). As shown, the mean turn-around radius is lower than the upper limit in the whole range of m_vir. The blue solid line is obtained by putting the global average values of A and n listed in Table 1 into Equation (7) and solving it for r_t. Note that the blue solid line is in good agreement with the filled circles in the entire range of m_vir, which is consistent with the claim of Lee & Yepes (2016) that the average turn-around radii can be computed by using the average bound-zone velocity profile with the two parameters set at the universal best-fit values. Figure 6 plot the same as Figure 5 but for the case of DF5 model. As shown, similar to the GR case, the DF5 model yields the average turn-around radii lower than the upper limit ${r}_{t,u}$ in the entire range of m_vir. Note, however, that the gap between the average turn-around radii and the upper bound limit is narrower in the DF5 model than in the GR case.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Now, we are ready to calculate the odds of the bound violations, for which we exclude those central groups whose logarithmic masses exceed 14.5. Given that the DF5 model produces more massive groups in the highest mass section than the other two models, it might cause a bias in the calculation of the odds of the bound violations if those central groups with ${m}_{\mathrm{vir}}\geqslant 14.5$ were not excluded. Let ${r}_{t,k}$ denote the turn-around radii of the central objects with mass m_vir in the kth realization of each model. Define η as the ratio of ${r}_{t,k}$ to the upper bound ${r}_{t,u}$ as $\eta \equiv {r}_{t,k}/{r}_{t,u}$. Dividing the whole range of η into several small bins, each of which has the same length ${\rm{\Delta }}\eta$, and counting the numbers of the central groups whose ratios belong to each bin, $[\eta ,\eta +{\rm{\Delta }}\eta ]$, we first compute the probability density function, ${p}_{k}(\eta )$, and then integrate ${p}_{k}(\eta )$ over η, to derive the cumulative probability distribution, ${P}_{k}(\geqslant \eta )$ from the kth realization of each model. The first five columns of Table 2 list the values of ${P}_{k}(\geqslant \eta )$ at $\eta =1$ for the three models. These values equals the ratio of the bound violating central groups to the total number of the central groups in each realization.

Table 2.  Odds of the Bound Violations for Three Models

models	${P}_{1}(\geqslant 1)$	${P}_{2}(\geqslant 1)$	${P}_{3}(\geqslant 1)$	${P}_{4}(\geqslant 1)$	${P}_{5}(\geqslant 1)$	$\bar{P}(\geqslant 1)$
GR	109/537	108/539	89/554	102/527	95/535	0.188 ± 0.007
DF6	110/534	105/541	88/554	116/548	91/509	0.190 ± 0.008
DF5	130/532	128/539	119/544	116/508	108/519	0.225 ± 0.006

Download table as:  ASCIITypeset image

Finally, we take the average of the cumulative probability functions over the five realizations for each model as

$$\begin{eqnarray}&&\bar{P}(\geqslant \eta )=\displaystyle \frac{1}{5}\displaystyle \sum _{k=1}^{5}{P}_{k}(\geqslant \eta ).\end{eqnarray} \tag{ 8 }$$

Figure 7 plots $\bar{P}(\geqslant \eta )$ versus η with the Bootstrap errors for the three models. As can be seen, the DF5 (GR) model yields the highest (lowest) values of $\bar{P}(\geqslant \eta )$ in the range of $\eta \geqslant 0.8$. In other words, in the DF5 model the bound violation occurs relatively more frequently than in the GR model. Meanwhile, no significant difference is found between the GR and the DF6 models, as expected.

Zoom In Zoom Out Reset image size

The sixth column of Table 2 lists the average odds of the bound violations, $\bar{P}(\geqslant 1)$, for the three models. In the GR model the odds of the bound violations is 0.188 ± 0.007, while in the DF5 model it is 0.225 ± 0.006. As speculated, in the DF5 model, due to the fifth force the bound, violations occur relatively more frequently. The signal-to-noise ratio for the difference in the odds of the bound violations between the two models is found to be as high as 4.2. This result indicates that the odds of the bound violations of the turn-around radii of the central groups can be powerful indicators of the presence of MG.

To examine if the odds of the bound violations depend on the mass threshold of the central groups, we increase the value of ${M}_{\mathrm{vir},\mathrm{th}}$ and rederived $P(\geqslant \eta )$ by repeating the whole process described in the above. Figures 8 and 9 plot the same as Figure 7 but for the cases of ${M}_{\mathrm{vir},\mathrm{th}}=3\times {10}^{13}\,{h}^{-1}\,{M}_{\odot }$ and ${M}_{\mathrm{vir},\mathrm{th}}=5\times {10}^{13}\,{h}^{-1}\,{M}_{\odot }$, respectively. As shown, no substantial change is made by increasing the values of ${M}_{\mathrm{vir},\mathrm{th}}$. Due to the lower numbers of the central groups for these two cases, however, the cumulative probabilities have larger errors and the odds of the bound violations have lower signal-to-noise ratios of ∼3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size