Density profiles and voids in modified gravity models

TL;DR

This paper investigates how voids and halo properties in the universe can reveal differences in gravity strength on large scales, offering new ways to test modified gravity models beyond traditional methods.

Contribution

It demonstrates that void and halo density profiles, along with their abundances, can serve as complementary probes to constrain modified gravity models, highlighting effects on large-scale structure formation.

Findings

Weaker large-scale gravity models have fewer large voids and massive halos.

Stronger gravity leads to more concentrated halos and thinner void walls.

Void and halo profiles provide additional constraints on gravity models.

Abstract

We study the formation of voids in a modified gravity model in which gravity is generically stronger or weaker on large scales. We show that void abundances provide complementary information to halo abundances: if normalized to the CMB, models with weaker large-scale gravity have smaller large scale power, fewer massive halos and fewer large voids, although the scalings are not completely degenerate with $\sigma_8$. Our results suggest that, in addition to their abundances, halo and void density profiles may also provide interesting constraints on such models: stronger large scale gravity produces more concentrated halos, and thinner void walls. This potentially affects the scaling relations commonly assumed to translate cluster observables to halo masses, potentially making these too, useful probes of gravity.

Figures

Equations

$$\varphi (r)=Gm\frac{1+\alpha (1-e^{-r/r_s})-\alpha (1-e^{-r/r_c})}{r}$$

$$\frac{d^{2}r}{d{t}^2}=F=-\nabla \Phi (r).$$

$$\frac{1}{2}{\left(\frac{dr}{dt}\right)}^2+\Phi (r)=C,$$

$$S_i\equiv \int \frac{dk}{k}\frac{k^3{P}_i(k)}{2{\pi }^2}{W}^2(k{r}_i),$$

$$\frac{dn}{d\ln M}d\ln M=\frac{\bar{\rho }}{M}f(S_i)d{S}_i.$$

$$\Phi (r)$$

$$\mathcal{F}(x)=\frac{x\cosh (x)-\sinh (x)}{x^3/3}.$$

$${\nabla }_r\Phi =\alpha \frac{G\delta M_i}{r_i^2}\frac{r_i}{r}\frac{F_i(r_s)-F_i(r_c)}{2}$$

$$F_i(r_x)$$

$$F_{\alpha }(r)$$

Full text

††preprint: APS/123-QED

I Introduction

A wealth of observations, from WMAP, supernovae Ia, galaxy clustering on large scales, and cross-correlations between galaxies and the CMB, suggest that we live in a spatially flat Friedmann-Robertson-Walker universe currently dominated by either a cosmological constant or repulsive dark energy. The best fit for the dimensionless energy density parameters are $\Omega_{m}=0.28$ and $\Omega_{\Lambda}=1-\Omega_{m}$, within the concordance $\Lambda$CDM Model Dunkley et al. (2009). However, because this requires the vast majority of the energy density in the universe to be in two unknown substances, dark matter and dark energy, there is considerable interest in alternative interpretations of what the data imply.

The key equation of General Relativity, Einstein’s equation, relates the curvature and the expansion rate of the Universe to its matter and energy content. The current paradigm is to modify the matter content of the universe, by including dark matter and dark energy, to account for observations. Instead, however, we might modify how the universe curves in response to matter, which would mean modifying our theory of gravity. There are many ways that one could modify gravity B-D ; milgrom ; bekenstein ; CarrollST ; fr1 ; fr2 ; fr3 ; Cham ; dgp ; dgpSims ; dgpPT ; dgpSim2 ; lss ; mannheim ; diaferio ; huReview , but the one that we will focus on in this paper is purely phenomenological. The assumption is that it is merely a modification to the potential, changing the form from that of standard gravity to:

[TABLE]

This potential (with $r_{c}\to\infty$ has been studied before ngp ; shirata1 ; sealfon ; fritz ; shirata2 and in some sense, is an interesting limiting case of some $f(R)$ models fR ; fRsims . In what follows, we will use $r_{c}\gg r_{s}$ but finite, following mss09 . The precise choice of $r_{c}$ does not matter for any of our conclusions.

In mss09 we studied the formation of nonlinear objects in this model – the first study of virialized dark matter halo abundances in any modified gravity model. In Section II, we study the evolution of the structures that are in some ways the opposites of halos, namely voids. Along the way, we show that the density profiles of objects in these models exhibit interesting departures from standard gravity models. Section III shows how we use insights from the evolution model to estimate void abundances. Section IV shows that our model captures the essence of how void abundances depend on $\alpha$ and $r_{s}$ in numerical simulations of structure formation with this modified potential. A final section summarizes our results. An Appendix shows the corresponding trends for halo profiles; these potentially allow X-ray observations to provide powerful probes of modifications to gravity.

II Evolution of underdensities

This section describes the evolution of initially underdense spherically symmetric regions for the modified gravity model described in the previous section.

II.1 The spherical tophat model

Following gg72 , we consider the evolution of a spherical patch of density different from the background in an expanding universe. The spherical collapse calculation begins with the statement that the force driving the acceleration is related to the gravitational potential by

[TABLE]

This can be integrated once to get

[TABLE]

where $\Phi(r)$ is the integral of the potential over the mass distribution, and $C$ is the total energy of the patch, which is constant. In standard gravity, the potential of a shell of mass $M$ is the same as that of a point mass at the center of the sphere, so $\Phi(r)$ reduces to $GM(<r)/r$. The constant $C$ can be related to the initial overdensity and/or expansion rate of the patch: the initial expansion rate is given by the Hubble expansion rate of the background in which the patch is embedded, namely in comoving coordinates $\dot{x}_{i}=0$, so $\dot{r}_{i}=\dot{a}_{i}x_{i}=\dot{a}_{i}/a_{i}(a_{i}x_{i})$, so $(dr/dt)_{i}=H_{i}r_{i}$. Including a cosmological constant presents no conceptual difference.

In standard gravity one can directly solve this equation. The solution is a cycloid, which is usually written in parametric form involving sines and cosines (but see ls08 for simple but accurate approximation). Only objects which are sufficiently overdense initially will collapse (meaning the size of the patch has become vanishingly small) at the present time. The critical value of the initial overdensity required for collapse today, $\delta_{c}$, does not depend on the initial size of the patch gg72 ; lc93 . This scale independence of $\delta_{c}$ is a result of Birkhoff’s Theorem: the evolution of a tophat sphere is the same as that given by Friedmann’s equations, so the actual size of the patch drops out. (Indeed, for triaxial perturbations, $\delta_{c}$ depends on the size and shape of the initial patch smt01 .)

When gravity is modified, things are no longer so simple. For example, when the potential is given by equation (1), Birkhoff’s theorem no longer applies: A particle offcenter in a uniform spherical shell will feel a force from the shell because the force no longer varies as $1/r^{2}$. This has two consequences. First, equation (2) can still be integrated once to get equation (3), only now $\Phi(r)$ has contributions from both the internal and external mass distributions. We can get $\Phi$ by integrating equation (1) of the mass distribution, about which more shortly, leaving $C$ and the initial value for $dr/dt$ to be determined. As before, $C$ is the total energy (constant in time), and we set $(dr/dt)_{i}=H_{i}r_{i}$.

Second, whereas evenly spaced concentric shells remain evenly spaced in the standard tophat model, this is no longer the case when the potential is modified. This is most easily seen by writing $\nabla\Phi$ as the sum of three terms: that due to the Newtonian $1/r$ part, and those due to the Yukawa term from shells internal and external to the point (or shell) of interest (see Appendix A for details). If $\delta$ denotes the overdensity of the source shell, then the internal and external contributions are both proportional to $\alpha\delta$, but they have opposite signs. As a result, the initial tophat perturbation develops a nontrivial density profile (see Appendix A for examples), meaning interior shells may cross one another as the initial ‘boundary’ around them collapses. When $\alpha>0$, initially overdense perturbations become more centrally concentrated than when $\alpha=0$, as the contribution from the modified part of the potential pulls mass away from the boundary in both directions (provided the boundary is comparable to $r_{s}$). When $\alpha<0$, then the perturbation develops a ridge at its boundary, as mass is pulled towards the boundary from both directions.

II.2 Density profiles

For voids, much of the preceding discussion remains valid. In standard gravity ($\alpha=0$), the solution is not a cycloid, but similar: the sines and cosines are replaced by their hyperbolic counterparts (however, the simple approximation of ls08 remains accurate even in this regime). When $\alpha\neq 0$, then the previous discussion continues to apply. However, because the sign of $\delta$ is reversed, the dependence on $\alpha$ is also reversed. As a result, mass moves towards the boundary of the profile from both sides when $\alpha>0$, but away from it when $\alpha<0$.

Figure 1 shows the evolution of a number of concentric shells within which the initial density was uniform, but slightly below that of the background. We used the methods of mss09 to calculate the evolution of the shells. (Note, however, that whereas few shells are needed to accurately estimate the evolution of the boundary of the initial tophat perturbation, many more are needed to study the evolution of the density profile.) When $\alpha=0$, then the expansion of the inner shells pushing against the exterior shells leads to the formation of a ridge (svdw04, , and references therein). This remains true when $\alpha\neq 0$, although, when $\alpha>0$ then the ridge is smaller in physical size, whereas the opposite is true for negative $\alpha$. In particular, for the $\alpha<0$ case, the ridge extends further and trails off much more slowly. In addition, when $\alpha\neq 0$, then the density in the interior region does not remain constant – this is clearly more evident in the negative $\alpha$ case, but also evident in the positive $\alpha$ case. Whereas the depedence of the interior density profile on $\alpha$ may be difficult to detect in observations, the effect on the ridge may be observable. Of course, to really exploit this effect, one must study more realistic initial density profiles – this is beyond the scope of our work.

II.3 Critical underdensity for void formation

One sense in which underdensities are more difficult to study than overdensities is that the concept of a critical time, the analogue of the collapse time, is not well-defined. While it is clear that a halo is collapsed at that time when the radius of its outermost shell reaches zero, what is the appropriate condition for undersities? For standard gravity, one uses the condition of “shell-crossing” which occurs when initially interior shells first cross initially exterior shells: the objects for which is true are about $0.2$ times denser than the background svdw04 . Because of the ambiguity associated with shell crossing in modified gravity models, we will use this value, an underdensity of $-0.8$, to define voids even when $\alpha\neq 0$. This critical nonlinear underdensity corresponds to some critical initial underdensity, $\delta_{v}$. When $\alpha=0$, $\delta_{v}\approx 2.86$ independent of the size of the underdensity, for the same reason that $\delta_{c}\approx 1.686$ is independent of mass. However, when $\alpha\neq 0$, we know $\delta_{c}$ depends on halo mass mss09 , and so we expect $\delta_{v}$ to depend on void volume.

To find this dependence, we study the evolution of underdense patches starting from a large grid of initial underdensities and sizes. Once we have the scale dependence of $\delta_{v}(R_{i})$, we can relate it to a mass scale by $M=(4\pi/3)R_{i}^{3}\bar{\rho}_{i}(1+\delta_{i})\approx(4\pi/3)R_{i}^{3}$, the last relation holding because $\delta_{i}$ is small. Thus, we can calculate $\delta_{v}(M)$. Intuitively, we expect that if, in the course of its evolution, the size of a patch never exceeds $r_{s}$, then we expect the large-scale modification to gravity is inconsequential. Hence, for sufficiently small voids, which form from smaller patches in the initial field, we expect $\delta_{v}$ to be the same as in standard gravity. For large voids, the effect of the modification should be stronger. For positive $\alpha$ voids should be easier to form, so we expect $|\delta_{v}|$ to be smaller, whereas for negative $\alpha$ the opposite should be true.

The solid curves in Figure 2 show that $\delta_{\rm v}(M)$ has the expected dependence on $M$ and $\alpha$. Note that we have chosen to express the scale dependence of $\delta_{v}$ in terms of $M$, for ease of comparison with the scale dependence of $\delta_{c}$ shown in Fig. 2 of mss09 , and reproduced here as the dashed curves. Of course, we could have expressed it in terms of the initial size $R_{i}$, or in terms of the final size $r_{v}$. This is because we have defined voids as being 0.2 times the density of the background, making $r_{v}=5^{1/3}\,R_{i}\approx 1.7\,R_{i}$: the comoving radius of a void is 1.7 times larger than it was initially. A close inspection of the figure will show that $\delta_{v}$ departs from its $\alpha=0$ value at slightly smaller mass scales than does $\delta_{c}$. This is because voids expand, so smaller mass scales can eventually exceed $r_{s}$ in size, and so notice that $\alpha\neq 0$.

III Statistics of void abundances

When combined with the spherical evolution model, excursion set methods bcek91 ; lc93 ; rks98 ; st02 ; ls09 are commonly used to study the abundance of virialized objects. They can easily be extended to the study of voids svdw04 . In short, this is done by relating the properties of such objects with the initial (rather than linearly evolved) density fluctuation field mss09 .

III.1 The excursion set method

Now, as shown, when the potential is modified, then $\delta_{v}$ is no longer scale-independent. Because it depends on mass, the relevant excursion set problem is one with a ‘moving’ rather than ‘constant’ barrier, so it is of the type first studied by st02 . Figure 3 shows the result of using this formalism to estimate the abundance of voids. Briefly, making this estimate requires that one generate an ensemble of random walks in the $(\delta_{i},S_{i})$ plane, where $S_{i}\equiv\sigma^{2}_{i}(M)$ is the variance in the initial fluctuation field when smoothed on scale $r_{i}$. Since $\sigma_{i}(M)$ is a monotonic function of $M$, the variables $S_{i}$, $M$ and $r_{i}$ are essentially equivalent to one another. In particular,

[TABLE]

where $W(x)=(3/x^{3})\,(\sin x-x\cos x)$. For halos, one then finds the first crossing distribution $f(S_{i})dS_{i}$ of the ‘barrier’ $\delta_{ci}(M)=\delta_{ci}(S_{i})$. The abundance of objects is given by

[TABLE]

For voids all this is slightly more complicated. The relationship between $f(S_{i})$ and $dn/d\ln M\,d\ln M=(\bar{\rho}/M)\,f(S_{i})dS_{i}$ is unchanged, but one must account for a subtlety in the agrument. To estimate halo abundances, one is careful to choose the first upcrossing of the collapse barrier, thus giving the largest smoothing radius which is sufficiently overdense that it will collapse today. The larger scale environment is irrelevant: for example, if the larger scale density is negative, the end result would just be a halo which formed inside a void. With voids, though, one must worry about underdense regions which are surrounded by large scale overdensities which would have collapsed by the present time – the void-in-cloud problem. Namely, if on some scale the random walk crosses below the critical density for void formation, but on a larger scale it crossed above the barrier for halo formation, then this would be a void which was crushed as the region around it collapsed. Therefore, such walks should not contribute to the estimated void abundance. As a result, rather than a one barrier problem, we instead have a two barrier problem svdw04 : $\delta_{v}$ must be crossed before $\delta_{c}$, is crossed. (A more conservative version of this argument sets the barrier associated with the collapsing object to equal that for turnaround rather than full collapse.)

Because there are two barriers, the first crossing distribution has a peaked shape, with the distribution being cut off at both high and low masses, on the high end due to the rarity of such extreme underdensities, and at the low end, because $\delta_{c}<|\delta_{v}|$. In large part, we will not be interested in the intricacies of which halo formation barrier we should use, as this largely affects only the low mass end of the mass function, where the modification to gravity is smaller.

Figure 3 shows the resulting void mass functions for standard and modified gravity. (We do not show the extremely low mass regime, which is most sensitive to the void-in-cloud cutoff.) Notice that when $\alpha$ is positive there are more large voids and fewer small voids, whereas the opposite is true when $\alpha$ is negative. This is from two effects: one is that any given walk tends to cross later when gravity is weaker ($\alpha<0$) as the barrier is lower (recall $\delta_{v}$ is negative) and so voids tend to be shifted towards lower mass. The second effect is that because it is harder to form high mass halos in weaker gravity, some walks that cross $\delta_{c}$ in stronger gravity may not cross $\delta_{c}$ for weaker gravity, however due to the fact that they are near $\delta_{c}$ at high mass, it is still difficult for such a walk to cross $\delta_{v}$ even at low mass. The first effect is likely more significant, but both contribute to the difference at low mass.

III.2 Initial conditions

One thing to note about the histograms in figure 3 is that the difference between the standard gravity and modified gravity void mass functions comes from two sources. Perhaps the more obvious is that modifying gravity modifies the behavior of individual structures, changing how quickly they form. This information is encapsulated in the barrier shape, and hence is accounted for using excursion set methods. The result of a shifted barrier is a shift in the function $f(S)$, the first crossing distribution. The other source is the modification to the power spectrum. In the study above, we assumed that the initial power spectra were the same for all the models – what is usually called CMB-normalization. However, because gravity is modified, the late time power spectra (in both linear and nonlinear theory) depend on $\alpha$ (and $r_{s}$). So now there are two other possible cases. We could consider modified gravity but with the initial conditions chosen so that the $z=0$ power spectrum is that of standard gravity, or we could consider standard gravity but with the initial conditions modified so that the $z=0$ power spectrum is that of modified gravity. In excursion set terms, either of these choices affects the relationship between $M$ and $S$.

For halos, whether one chooses to match the initial power spectrum or the final power spectrum matters mss09 , because it changes the mapping between $S$ and $M$ (or $R$), and this mapping affects two steps of the calculation: the transformation of $\delta(M)$ into $\delta(S)$, and the relation between $f(S)$ and $n(M)dM$. The same is true for voids so, in what follows, we will consider some of these other possiblities as well.

IV Comparison with simulations

We now compare our spherical evolution predictions for void abundances with measurements in the simulations of fritz . These simulations followed the evolution of $128^{3}$ particles in a periodic box of size $100h^{-1}$Mpc, for various choices of $\alpha$ and $r_{s}$. In all cases, the background cosmology was flat $\Lambda$CDM with $\Omega=0.3$, and the particle mass was $1.1\,\times\,10^{10}M_{\odot}$. The $\alpha=0$ simulation, with standard initial conditions has $\sigma_{8}=1.0$ at $z=0$. The corresponding runs for $\alpha=0.5$ and $\alpha=-1$ have $\sigma_{8}=1.10$ and $0.84$ respectively. Following our discussion of how the halo and void counts depend on the shape and normalization of the initial power spectrum, we also study results from $\alpha=0$ simulations in which the initial power spectrum was modified so that, at $z=0$, it has the same shape as the two $\alpha\neq 0$ cases ($\sigma_{8}=1.10$ and 0.84).

The simulations were analyzed using the void finder of joergVF . Figure 4 shows the halos and voids found in three runs, all with the same iitial conditions. Notice that the voids are largest when $\alpha=0.5$ and smallest when $\alpha=-1$, in qualitative agreement with expectations. Since the $\alpha=0.5$ run has more large scale power at $z=0$ this is not surprising. Figure 5 shows the results of a more quantitative comparison: although qualitatively similar, the predicted effects (solid lines) for large voids are larger than we find in the simulations, and the effect on small voids is smaller than is observed.

We explore the dependence on choice of initial conditions in two steps. The top panel in Figure 6 shows the ratio of void counts in two standard gravity runs, one with initial conditions modified to produce the same $z=0$ power spectrum as $\alpha\neq 0$, and the other with the standard $\sigma_{8}=1$ initial conditions (note that the modification to the initial conditions is different for the two values of $\alpha$ shown). This plot is qualitatively similar to the previous one, because stronger large-scale gravity ($\alpha>0$) is qualitatively like having more large scale power, but the differences between these two plots shows that the result of modifying gravity is not degenerate with modifying the shape and normalization of the initial power spectrum. The bottom panel shows this information slightly differently: here, the void counts in the $\alpha\neq 0$ run are ratioed to the counts in the $\alpha=0$ run in which the initial conditions were tuned to produce the same $z=0$ power. The fact that the ratios are not unity implies that there is more information in the void size distribution than in the power spectrum itself.

Finally, Figure 7 shows the ratio of void counts when it is the initial conditions in the the modified gravity runs which have been tuned to produce the same $z=0$ power spectrum (rather than tuning the $\alpha=0$ initial conditions). Unfortunately, we do not have simulations of this case, but we again see that the ratio is predicted to differ from unity, indicating that modifications to gravity are not degenerate with changes to the initial conditions.

V Conclusions

We studied a particular modification to large-scale gravity, which has two free parameters: a scale $r_{s}$, and an amplitude $\alpha$ (equation 1). Previous work has shown that halo abundances can be modeled well enough mss09 to use them to constrain such models. We present an example of this in Appendix B. We argued that since voids are so large, and almost fill space, they too can be useful probes. In particular, the study of the formation of large voids provides a useful counterpart to the the study of the formation of massive halos because both probe the nature of gravity on large scales. This is timely, given that large, statistically representative void samples are only just becoming available (e.g. voidp, ; voids, ).

Whether stronger large-scale gravity produces more or fewer large voids depends on how one normalizes the models. When normalized to have the same power spectrum at early times (the most common choice), stronger gravity produces more large voids (Figure 3). This case also produces more massive halos, so the net result is qualitatively like having standard gravity with more large scale power (Figure 4). Although there are quantitative differences which allow one to distinguish between a larger normalization and modified gravity (Figures 5–7), we note that this fact alone – a mismatch between the CMB- and later-time determinations of the amplitude of the power spectrum on 10 Mpc scales when standard gravity is assumed – are generic signatures of modifications to large-scale gravity. Note that this particular signature has the same sign for halos and for voids – stronger large scale gravity means more massive voids and more large voids. In contrast, if gravity is standard but the initial conditions were non-Gaussian, then halo and void abundances are modified in qualitatively different ways, at least for the local non-Gaussian initial conditions of current interest lsd09 .

However, if normalized to have the same power at $z=0$, the dependence of void (and halo) abundances on whether large-scale gravity is stronger or weaker depends on how one chooses to do this. If this is done by increasing (decreasing) the initial power spectrum in models with $\alpha<0$ ($\alpha>0$), then one still expects the stronger gravity models to produce larger voids (Figure 7), although the predicted abundances are quantitatively different. However, if one modifies the standard gravity initial conditions to match those of the $\alpha\neq 0$ models (increase/decrease initial large scale power to match $\alpha>0$/$\alpha<0$), then one predicts more (fewer) large voids when $\alpha<0$ ($\alpha>0$) relative to the $\alpha=0$ case. This dependence on how the models were normalized is qualitatively similar to the trends seen for massive halos mss09 .

We presented a method for estimating these effects which is in good qualitative greement with the simulations, but there are quantitative differences. However, the agreement is not as good as it was for describing halos. Larger simulations are required to determine if this is due to the relatively small size of the simulation boxes available to us, or to some more fundamental problem with our analysis (see smt01 ; ls09 for discussion of the drawbacks of the excursion set approach).

We also showed that the density profiles of voids and halos may also provide interesting probes of modified gravity. When the initial profile is a tophat, then the evolved halo profile in the $\alpha<0$ case generally has a cusp at the virial radius (mass flows towards the boundary), whereas when $\alpha>0$ the halos are more centrally concentrated (Figure 8). The structure of the void walls also depends slightly on $\alpha$ (Figure 1). We provided tentative evidence that the profiles of halos in simulations do depend on $\alpha$ (Figure 9), and hope that these initial measurement motivate further study of this interesting problem. It will be interesting indeed if these trends persist in simulations with better mass and force resolution, because it is conceivable that effects of this magnitude will soon be measured by weak lensing surveys. The more dramatic effects associated with initially tophat profiles potentially provide an interesting X-ray signature of modified gravity – it would be interesting to simulate such models using SPH codes.

In conclusion, we expect our results will prove useful in studies which use the large scale distribution of galaxies, and the structure of galaxy clusters, to constrain large scale modifications of gravity.

Acknowledgments

We would like to thank Fritz Stabenau for providing the outputs from his simulations in electronic form, and Joerg Colberg for running his void finder algorithm on these outputs. This work was supported in part by NSF-AST 0908241.

Appendix A Halo density profiles

A complementary study to the one in the main text regarding the evolution of the density profile of voids is a study of the evolution of halo density profiles. Figure 8 shows one such study. In this case, the three objects all started as tophats with the same initial overdensity, but different values of $\alpha$. The solid and dashed curves show the density profile for two different values of $r_{s}$ and at three different times, one corresponding to turnaround, one to when one shell reaches a density of 200 times the background, and the final to when the density enclosed in the outer shell reaches 200 times the background. Since gravity has a different strength in the three cases, and the patches all started with the same overdensity, the final collapse happeds at different times. For positive $\alpha$, objects collapse faster, whereas for negative objects collapse slower, with the difference decreasing as $r_{s}$ increases. This is why the effect is noticeable earlier for negative $\alpha$: the turnaround time is later than for standard gravity, so the deviation is greater.

The evolved profiles in the two modified gravity cases are strikingly different from standard gravity, with the formation of sharp spikes at the edge of the halo in the negative $\alpha$ case, and the formation of a large central peak in the positive $\alpha$ case. Whether these profiles would hold through virialization is unknown, but since the effect is so strong during the period before collapse it seems likely that the profiles of virialized haloes would be affected. If so, then halo profiles might provide new interesting constraints on $\alpha$. However, to place such constraints, one would have to extend our analysis to more realistic initial profiles which do not have sharp edges – although it is likely that the formation of the cusp at the boundary when $\alpha<0$ will remain (see Figure 9 and related discussion). Note that for larger $r_{s}$, the difference from standard gravity is smaller, because the perturbation spends essentially all its time at scales which are smaller than $r_{s}$.

These trends are perhaps most easily understood by writing the potential due to an extended source at scales that are beyond its boundary. The potential at $r$ due to a top perturbation which contains mass $M$ within radius $R$ can be written as the sum of the standard gravity piece plus a term which arises from the modification:

[TABLE]

where $\delta M\equiv M-\bar{\rho}4\pi R^{3}/3$ and

[TABLE]

Note that ${\cal F}\geq 1$ for all $x$. However, ${\cal F}\to 1$ when $x\ll 1$, so it is easy to see that when $r_{s}\gg R$ then $\Phi$ looks just like it does standard gravity. If the overdensity profile is a power-law of slope $-1$ within $R$, so $\delta M=2\pi R^{3}$, then the expression for $\Phi$ remains true, but ${\cal F}(x)\to[\cosh(x)-1]/(x^{2}/2)$: again, ${\cal F}\geq 1$ for all $x$, and ${\cal F}\to 1$ when $x\ll 1$. In practice, we are most interested in the regime in which $R$ is smaller than a few times $2r_{s}$, for which the difference between these form factors is small. This means that one can approximate the force outside the perturbation by assuming it is a tophat, and not worrying about its internal structure. In fact, this suggests that one can solve for the evolution of the boundary of the perturbation without worrying about the fact that the internal structure is changing; in essence this is why the analysis in mss09 was so straightforward.

The evolution of the profile within the perturbation is more easily studied using the force, which is $\nabla_{r}\Phi(r)$, and can be split into a Newtonian part, plus terms due to the modification. These in turn can be split into the contributions from shells internal and exterior to $r$. If we write the mass overdensity in the $i$th shell as $\delta M_{i}=\bar{\rho}(a)\,\delta_{i}\,4\pi r_{i}^{2}\,dr_{i}$ then

[TABLE]

where

[TABLE]

when $r_{x}\leq r$ (i.e., from the internal shells) and

[TABLE]

if $r_{x}>r$ (the external shells). Given these expressions, we can (and do) integrate to get the total force at a radius $r$.

To get a qualitative feel for the effect of $\alpha\neq 0$, it is instructive to study the net force associated with a tophat perturbation, for which $\delta_{i}$ is the same for all shells within the perturbation. (Note that all the analysis in the main text, and the figures in this Appendix, were made by solving for the exact evolution, i.e., not assuming that the initial tophat remains a tophat.) In this case the modification to the Newtonian force $-GM/r^{2}$ is given by

[TABLE]

In the present context, we are also interested in the profile within the perturbation. A similar analysis of $\Phi(r)$ when $r<R$ shows that the modification to the force is

[TABLE]

The term in square brackets is always positive, so the correction is proportional to the product of $\alpha$ and $\delta M$. Hence, when both are positive, then the resulting profile becomes steeper than when $\alpha=0$: one might generically expect halos to be more concentrated. However, for overdense perturbations when $\alpha<0$, the net force towards the center is smaller than when $\alpha=0$, so we expect the profile to evolve away from a tophat in the sense of becoming more concentrated, not at the center, but at the outer boundary.

We have looked for this effect in the simulations we study in the main text. Unfortunately, because the particle mass in these simulations is rather large (the runs were designed to study larger scale structures), we can only estimate halo density profiles for the rarest most massive objects. In the $\alpha=(-1,0,0.5)$ runs, the most massive object had $\log_{10}(M/h^{-1}M_{\odot})=(14.85,15.15,15.58)$ and virial radii, $R_{\rm vir}/h^{-1}{\rm Mpc}=(1.31,1.66,2.36)$. Figure 9 shows the profiles of these three objects. In each case, we show the expected enclosed density profile for an object of this mass if $\alpha=0$. In all cases, the measured profiles lie below this line on scales smaller than $R_{\rm vir}$, suggesting that our simulations are not able to properly resolve small scale structures. However, there is a hint that the $\alpha=-1$ object has slightly more than the expected amount of mass just beyond the virial radius, whereas the opposite is true for the $\alpha=0.5$ object. It will be interesting if these trends persist in simulations with better mass and force resolution, because it is conceivable that effects of this magnitude will soon be measured by weak lensing surveys.

Finally, we note that there is another sense in which the modified gravity model is more complicated. As the object pulls itself together from the expansion, it will eventually reach turnaround. In standard gravity, all its shells reach turnaround at the same time, so the total energy at turnaround is potential, and energy conservation means that this equals its initial energy. When $\alpha\neq 0$, however, the different shells turnaround at different times, so the total energy when the outermost shell turns around is no longer all potential.

Appendix B Constraining $\alpha$ and $r_{s}$

The main text studied a particular modification to large-scale gravity, which has two free parameters: a scale $r_{s}$, and an amplitude $\alpha$. Previous work has sought to constrain the values of these parameters by using the shape of the galaxy power spectrum shirata1 ; shirata2 . Although, on large scales, the amplitude of the power-spectrum depends on $\alpha$, potentially providing a powerful constraint, this is offset by the fact that galaxies are biased tracers of the dark matter field, and this bias is unknown. Even allowing for the simplest, linear bias relation removes most of the sensitivity of this method to the information contained in the amplitude of the late-time $P(k)$. Our succesful model for halo abundances, however, is sensitive to precisely this difference, so, we can use current high mass cluster counts to constrain the allowed range of $(\alpha,r_{s})$.

Figure 10 shows the result of comparing our predicted mass functions with the measured counts of vikhlinin when we normalize to the same (CMB) initial conditions for all $\alpha$ (blue contours), and when normalized to have the same linear theory $P(k)$ at $z=0$ (black). The inner most contours represent one, two, and three $\sigma$ confidence levels. The constraints for the CMB normalized curves are slightly weaker than the others. This can be traced to Figure 6 in mss09 , which shows that the solid line is closer to unity than the dotted line in the mass range of clusters, $10^{14}-10^{15}\,M_{\odot}/h$. This is because the mass function is determined by the shape of the barrier and the shape of the initial power spectrum, so that in the case of CMB normalization only one of these components changes (the barrier), whereas for cluster normalization both change, leading to a larger change in the mass function – and hence to stronger constraints. Note in particular, that the constraints that we derive here are roughly comparable to those derived by shirata1 ; shirata2 . However, if halo density profiles are indeed sensitive to $\alpha$, as our analysis suggests, then this will affect the scaling relations that are usually assumed to convert X-ray observables into masses. This, in turn, will alter the derived constraints. We look forward to the day when better simulations are available, so the question of how X-ray cluster scaling relations depend on $\alpha$ and $r_{s}$ has been settled.