The $f(R)$ halo mass function in the cosmic web

TL;DR

This paper investigates how modified gravity, specifically $f(R)$ models, influences halo properties and abundance within the cosmic web, using an extended formalism that aligns well with simulations on large scales.

Contribution

It generalizes the excursion set formalism to produce a halo mass function conditioned on large-scale structure, revealing environment-independent effects of $f(R)$ gravity on halo properties.

Findings

$f(R)$ affects halo abundances, spin, and concentration-mass relation.

No significant deviation from $ ext{Lambda}$CDM in fixed environment density.

Greater deviations occur in underdense environments and for high-mass haloes.

Abstract

An important indicator of modified gravity is the effect of the local environment on halo properties. This paper examines the influence of the local tidal structure on the halo mass function, the halo orientation, spin and the concentration-mass relation. We generalise the excursion set formalism to produce a halo mass function conditional on large-scale structure. Our model agrees well with simulations on large scales at which the density field is linear or weakly non-linear. Beyond this, our principal result is that $f(R)$ does affect halo abundances, the halo spin parameter and the concentration-mass relationship in an environment-independent way, whereas we find no appreciable deviation from $\Lambda$CDM for the mass function with fixed environment density, nor the alignment of the orientation and spin vectors of the halo to the eigenvectors of the local cosmic web. There is a general trend for greater deviation from $\Lambda$CDM in underdense environments and for high-mass haloes, as expected from chameleon screening.

Figures

a](#S3.F2.sf1)), and the peak-backgound split $\delta_{c}-\delta_{\textrm{env}}$ is a montotonically-decreasing function of $\delta_{\textrm{env}}$([Figure 2b](#S3.F2.sf2)). Compared to $\Lambda$CDM, we expect haloes to form at higher masses (earlier in the random walk) and to have more in low-density regions (where the fifth-force potential is unscreened). For a moving barrier like this, and under the assumption that trajectories in $(\delta,S)$ take uncorrelated steps, the diffuson equation admits a solution [[40](#bib.bib40)]:

c](#S5.F12.sf3) shows the spin distribution for all halos in the $\Lambda$CDM and $f(\mathcal{R})$ simulations and [Figure 12d](#S5.F12.sf4) the ratio of the distributions between both cosmologies for the four different environment types. The spin-parameter is seen to be higher in modified gravity which is in agreement with what was previously found in [[57](#bib.bib57)]. We also split the halo catalog into two samples $M<3\cdot 10^{13}M_{\odot}$ and $M>3\cdot 10^{13}M_{\odot}$ and computed the PDF in each of the samples. For the high-mass sample the PDF is very close to $\Lambda$CDM while for

Tables

Table 1. Table 1 : The N -body simulations used in this paper. The cosmological parameters used to generate initial conditions was: Ω m ​ 0 = 0.27 subscript Ω 𝑚 0 0.27 \Omega_{m0}=0.27 , Ω Λ ​ 0 = 0.73 subscript Ω Λ 0 0.73 \Omega_{\Lambda 0}=0.73 , h = 0.704 ℎ 0.704 h=0.704 , n s = 0.966 subscript 𝑛 𝑠 0.966 n_{s}=0.966 and σ 8 = 0.8 subscript 𝜎 8 0.8 \sigma_{8}=0.8 .

Box size (Mpc/h)	Particles per box	Cosmology	MG algorithm
$1024$	${1024}^3$	$\Lambda$CDM, $\left|f_{ℛ0}\right|={10}^{-5}$	Approximate
$256$	${512}^3$	$\Lambda$CDM, $\left|f_{ℛ0}\right|={10}^{-5}$	Full

Equations

$$S_J=\frac{1}{2}\int d^{4}x\sqrt{-g}\left[\mathcal{R}+f(\mathcal{R})\right]+\int d^{4}x\sqrt{-g}{\mathcal{L}}_{\mathrm{m}}[{\Phi }_i,g_{\mu \nu }],$$

$$S_E=\frac{1}{2}\int d^{4}x\sqrt{-\tilde{g}}\tilde{\mathcal{R}}+\int d^{4}x\sqrt{-\tilde{g}}\left[-\frac{1}{2}{\tilde{g}}^{\mu \nu }{\tilde{\nabla }}_{\mu }\varphi {\tilde{\nabla }}_{\nu }\varphi -V(\varphi )\right]+S_{\text{matter}}[{\Phi }_i,e^{\beta \varphi }{\tilde{g}}_{\mu \nu }]$$

$${\tilde{g}}_{\mu \nu }=e^{2\omega }{g}_{\mu \nu }\text{requiring}{e}^{-2\omega }\left(1+f_{\mathcal{R}}\right)=1,\varphi \equiv \frac{\omega }{\beta }$$

$$V(\varphi )=\frac{1}{2}\frac{\mathcal{R}{f}_{\mathcal{R}}(\mathcal{R})-f(\mathcal{R})}{{\left(1+f_{\mathcal{R}}(\mathcal{R})\right)}^2}$$

$$f(\mathcal{R})=-H_0^2{\Omega }_m\frac{\left(c_1{\left(\frac{\mathcal{R}}{m^2}\right)}^{n_{f(\mathcal{R})}}\right)}{1+c_2{\left(\frac{\mathcal{R}}{m^2}\right)}^{n_{f(\mathcal{R})}}}\text{where c1=Ωm6ΩΛc2}$$

$$|f_{\mathcal{R}0}|=n_{f(\mathcal{R})}\frac{|c_1|}{c_2^2}{\left(\frac{{\Omega }_m}{3({\Omega }_m+4{\Omega }_{\Lambda })}\right)}^{1+n_{f(\mathcal{R})}}$$

$$F=\frac{G{M}_{\mathrm{TH}}m}{r^2}\left(1+\frac{1}{3}\frac{\Delta R}{R_{\mathrm{TH}}}\right)$$

$$\frac{\Delta R}{R_{\text{TH}}}=\min \left\{\frac{3|f_{\mathcal{R}}^{\mathrm{TH}}-f_{\mathcal{R}}^{\mathrm{env}}|}{2{\Phi }_N},1\right\}$$

$$f_{\mathcal{R}}({\rho }_m)=f_{\mathcal{R}0}{\left(\frac{1+\frac{4{\Omega }_{\Lambda 0}}{{\Omega }_{m0}}}{\frac{{\rho }_m}{{\rho }_{m0}}+\frac{4{\Omega }_{\Lambda 0}}{{\Omega }_{m0}}}\right)}^{n_{f(\mathcal{R})}+1}$$

$$F=\frac{G{M}_{\mathrm{TH}}m}{r^2}\left[1+\frac{1}{3}{F}_{\mathrm{eff}}\right]$$

$$F_{\text{eff}}(a,R_{\mathrm{TH}},{\rho }_{\mathrm{TH}},{\rho }_{\mathrm{env}})$$

$$\frac{\Delta R}{R_{\text{TH}}}$$

$$S(M)\equiv \frac{1}{2{\pi }^2}{\int }_0^{\infty }dk{k}^{2}P(k,z=0)|W(kR){|}^2.$$

$$F(>M)=2{\int }_{{\delta }_c}^{\infty }P(\delta (S))d\delta =\text{erfc}\left(\frac{{\nu }_h}{\sqrt{2}}\right),$$

$$F(>M)=\frac{\exp \left(-c{\nu }_h^2\right)}{1+a{\nu }_h^b}\text{where}\{\begin{array}{ll}a& =1.529\\ b& =0.704\\ c& =0.412\end{array}$$

$$f(M)\equiv \frac{dF(<M)}{d\ln M}=\text{f}(S)\left|\frac{dS}{d\ln M}\right|=-⟨\text{f}\left({\delta }_c(S),S|{\delta }_{\text{env}},S_{\text{env}}\right){⟩}_{\mathrm{env}}\frac{dS}{d\ln M}$$

$$0=\frac{d^{2}y}{d\ln a^2}+\left(2-\frac{3}{2}{\Omega }_m(a)\right)\frac{dy}{d\ln a}+\frac{1}{2}{\Omega }_m(a)\left(\frac{1}{y^3}-1\right)y$$

$$M=\frac{4}{3}\pi {\rho }_{m0}{a}_i^3{R}^3=\frac{4}{3}\pi \rho r_i^3⟹y_i=1-\frac{{\delta }_i}{3}\text{and}\frac{d{y}_i}{d\ln a}=-\frac{{\delta }_i}{3}.$$

$$0=\frac{d^{2}y}{d\ln a^2}+\left(2-\frac{3}{2}{\Omega }_m(a)\right)\frac{dy}{d\ln a}+\left(1+F_{\mathrm{eff}}\right)\frac{1}{2}{\Omega }_m(a)\left(\frac{1}{y^3}-1\right)y,$$

$$g\left(S|S_{\mathrm{env}},{\delta }_{\mathrm{env}}\right)$$

$$k(S,x)$$

$$h(S)$$

$${\mathsf{T}}_{ij}=\frac{{\partial }_i{\partial }_j\varphi }{4\pi G\bar{\rho }}.$$

$${\nu }_e$$

$$\theta$$

$$\rho$$

$${\nu }_e-{\nu }_{\mathrm{th}}$$

$$f_1\left(\rho ,\theta |\alpha \right)$$

$$f_2\left(\rho ,\theta |\alpha \right)$$

$$p\left(\rho ,\theta ,{\nu }_e\right)=225\frac{\sqrt{5}}{2\pi }\rho \left({\rho }^2-{\theta }^2\right)\exp \left[-\frac{1}{2}\left(15{\rho }^2+5{\theta }^2+{\nu }_e^2\right)\right]$$

$$F_{\mathrm{vol}}={\int }_0^{\infty }d\rho {\int }_{-\rho }^{\rho }d\theta {\int }_{{\nu }_{\mathrm{th}}+f_1\left(\rho ,\theta |\alpha \right)}^{{\nu }_{\mathrm{th}}+f_2\left(\rho ,\theta |\alpha \right)}d{\nu }_{e}p(\rho ,\theta ,{\nu }_e).$$

$$p\left({\delta }_h|{\delta }_{\text{env}}\right)d{\delta }_h=\frac{d{\nu }_h}{\sqrt{2\pi \left(1-{ϵ}^2\right)}}\exp \left[-\frac{{\left({\nu }_h-ϵ{\nu }_e\right)}^2}{2\left(1-{ϵ}^2\right)}\right]$$

$$g\left(S|\alpha \right)$$

$$g\left({\nu }_{\mathrm{eff}}\right)$$

Full text

1 Introduction

The discovery of the accelerated expansion of the Universe [1, 2] is one of the biggest puzzles in modern cosmology. Having accepted the Einstein equations as the correct description of the interaction between matter and geometry, one is left with no alternative but to modify either side of the equation in order to account for this acceleration. On the one hand, modifying the stress-energy tensor necessitates the presence of an unknown substance dubbed dark energy; on the other, modifying the Einstein tensor require modifications to relativity at a more fundamental level. This poses the question: should one of these choices be favoured over the other? Dark energy—at its simplest a vacuum energy which acts in the form of an effective cosmological constant—is in excellent agreement with most observations so far, but suffers from some theoretical problems including fine-tuning [3]. Extended gravity theories—at their simplest adding extra complications to the Einstein-Hilbert action—have a secure theoretical motivation, but need to be carefully designed to satisfy observational constraints.

General Relativity has been tested to great precision on Earth and in the Solar System [4, 5, 6]. Thus, if it is modified then some form of screening mechanism [7, 8, 9, 10] is required to hide the effects in these regimes. As we have tested general relativity only in high density regimes (relative to the cosmic mean), this naturally leads to an environmental dependence on the modifications of gravity. In turn this would translate into an environmental dependence on observables [e.g. 11, 12]. Thus the very screening which modified theories of gravity invoked to evade local, Solar System tests can be harnessed as a detection mechanism by examining different regions of the cosmic web. This can also be used to distinguish it from other dark energy scenarios.

This paper examines a variety of halo properties in the different geometric structures of the cosmic web, namely: the abundance of halos of different masses; the alignment of the halo spins and shapes with the local tidal structure; the halo spin parameter and the concentration-mass relation. We provide an excursion-set-based approximation to the N-body results for the halo mass function. This synthesises two, major existing results: the unconditional mass function in $f(\mathcal{R})$ calculated by [13, 14, 15] and the classification of halo statistics in the cosmic web by e.g. [16, 17, 18]. In the same vein as those papers, our aim is to build a semi-analytic model which combines an empirically-derived halo mass function in the unconditional case for $\Lambda$CDM with a physical model derived using the tools of excursion set theory. Such a model, if sufficiently accurate, can then be used to cheaply compute observables for quantities that otherwise would require expensive numerical simulations.

The structure of the paper is as follows: We summarise our chosen model of modified gravity in Section 2, before reviewing the excursion set theory in Section 3.1 (general relativity) and Section 3.2 (modified gravity) and by extending this to dependence on the cosmic web in Section 3.3. The N-body simulations against which we benchmark our semi-analytical model are described in Section 4. In Section 5 we compare the accuracy of our theoeretical model to N-body simulations (Sections 5.1 and 5.2) and we analyse deviations from $\Lambda$CDM for each of the halo properties (Sections 5.3 and 5.4). We suggest possible generalisations and applications of this method in and conclude in Section 6.

2 Review of modified gravity

In this section we outline the differences between the modified gravity theory used in this paper and general relativity. We motivate our choice of $f(\mathcal{R})$ theory and describe some necessary conditions for which it is observationally viable.

An $f(\mathcal{R})$ theory [19] adds a scalar function to the Einstein-Hilbert action for General Relativity. It can be defined in the so-called Jordan frame via the action

[TABLE]

where we have chosen units such that $8\pi G=1$, the function $f(\mathcal{R})$ is a general function of the Ricci scalar $\mathcal{R}$ and $\Phi_{i}$ denotes all matter fields.

For a modified gravity model of this kind to be observationally viable, it must exhibit the chameleon screening mechanism [20]. Via a field redefinition and a conformal transformation we can turn the Jordan frame action into an equivalent Einstein frame one [21]

[TABLE]

where a tilde denotes Einstein frame quantities and we have performed a conformal transformation

[TABLE]

where $\beta=\sqrt{1/6}$, $f_{\mathcal{R}}\equiv df(\mathcal{R})/d{\mathcal{R}}$ and $V(\phi)$ is the scalar field self-interaction potential given by

[TABLE]

The $f(\mathcal{R})$ modification in the Jordan frame translates to a scalar-tensor theory of gravity in the Einstein frame where the scalar field $\phi$ is coupled to matter. In this formulation we can see that this model behaves as standard general relativity with the inclusion of a fifth-force mediated by the scalar field with coupling strength $1/3$, i.e. in the unscreened limit, the gravitational force can be enhanced by a factor of up to $4/3$. In the fully screened limit the gravitational force is not enhanced, i.e. this modification does not manifest itself in dense environments or on small scales [22]. In order to exhibit this behaviour, the scalar field self-interation potential must satisfy a number of constraints (viz. [20]), which again translates into constraints on the functional form of $f(\mathcal{R})$ [23].

The explicit $f(\mathcal{R})$ model we are working with in this paper is the Hu-Sawicki model of [24]. This is a well studied model known to exhibit chameleon screening [25]. It is defined by

[TABLE]

Rather than using $c_{2}$ as the free parameter we will instead express it in terms of $f_{\mathcal{R}0}$, the value of $f_{\mathcal{R}}$ in the cosmological background evaluated at $z=0$, using

[TABLE]

In this paper we will only consider $n_{f(\mathcal{R})}=1$ and $|f_{\mathcal{R}0}|=10^{-5}$.

To see how screening works, let us consider a top-hat overdensity of radius $R_{\textrm{TH}}$ and mass $M_{\rm TH}$. The Newtonian potential of the overdensity is $\Phi_{N}=\frac{GM_{\rm TH}}{R_{\rm TH}}$. As shown in [21] the gravitational force on a test-mass of mass $m$ outside the top-hat is approximately given by

[TABLE]

where the screening factor $\frac{\Delta R}{R_{\rm TH}}$ is given by

[TABLE]

where $f_{\mathcal{R}}^{\rm TH}=f_{\mathcal{R}}(\rho_{\rm TH})$ and $f_{\mathcal{R}}^{\rm env}=f_{\mathcal{R}}(\rho_{\rm env})$ are the scalar field values inside and outside the body respectively. For the Hu-Sawicki model the scalar field value in a region density $\rho_{m}$ and time (scale-factor) $a$ can be estimated to be

[TABLE]

where $\rho_{m0}=3\Omega_{m0}H_{0}^{2}M_{\rm Pl}^{2}$ is the average matter density in our Universe at the present time. When the overdensity is massive ($\frac{1}{\Phi_{N}}$ is very small) or is located in a very dense environment ($|f_{\mathcal{R}}^{\rm TH}-f_{\mathcal{R}}^{\rm env}|$ is very small) then $\frac{\Delta R}{R_{\textrm{TH}}}\ll 1$ and the fifth-force is screened. In contrast, when the overdensity is not massive (so $\frac{1}{\Phi_{N}}$ is very large) then $\frac{\Delta R}{R_{\textrm{TH}}}\approx 1$ and the force is $4/3$ the value of the Newtonian prediction. Thus we see that the modification to gravity is sensitive to both the halo mass and the environment density.

A more accurate expression for the screening factor, which we will use in this paper, was derived in [13] from results in [26, 15]. The gravitational force on a test-particle at the surface of a top-hat overdensity of radius $R_{\textrm{TH}}$ collapsing in an expanding background is given by

[TABLE]

where

[TABLE]

The background cosmology of this model is always very close to $\Lambda$CDM as long as $|f_{\mathcal{R}0}|\ll 1$ (which is required to satisfy local gravity constraints). We have three regimes: on cosmological scales, the background solution mimics $\Lambda$CDM; on local scales in high-density regions, the modification is screened in order to evade tight, Solar System constraints; on local scales in low-density regions the modifications of gravity are in full effect.

3 Halo abundances in cosmology

In this section we outline the excursion set in general relativity (Section 3.1) and then discuss the modifications induced by modified gravity (Section 3.2). Finally we present the key integral of this paper, which incorporates the dependence on the cosmic web (Section 3.3). In particular, we choose the Peacock mass function Equation 3.3 as our unconditional model in $\Lambda$CDM and we discuss which aspects of the excursion set theory in modified gravity are compatible with our model for the cosmic web. In this way, we ensure that our base model is an excellent fit to simulations, while also retaining the semi-analytical nature of the model.

3.1 The excursion set theory in $\Lambda$CDM

The Ansatz of the excursion set approach is to relate the fraction of the density field $\delta$ above a certain critical density $\delta_{c}$ to the cumulative fraction $F(>M)$ of mass contained in haloes above mass $M$.

The overdensity field $\delta$ is obtained from the Gaussian fluctuations in the fractional, linearly-evolved overdensity smoothed over a scale $R$ defined by a given Fourier-space window function $W(kR)$. In what follows we will relate the smoothing scale to a mass scale $M$ defined as the mass encompassed by the kernel $W(kR)$ in a homogeneous Universe.

The fluctuations of the overdensity field have some variance $\sigma^{2}$, and we define the variable $\nu\equiv\delta/\sigma$. In what follows we will refer to the variance at redshift zero $S\equiv\sigma^{2}(M,z)/D(z)^{2}$ (where $D(z)$ is the linear growth factor normalised to $D(z=0)=1$) as the resolution, which is related to the linear matter power spectrum as:

[TABLE]

The foundation of the excursion set formalism is the relation of the linear overdensity field in real space collapsing to a halo of mass $M$ to trajectories in the density field being absorbed by a critical density at resolution $S$. As shown in [27, 28], where they assume top-hat window functions in $k$-space, this problem reduces to finding the fraction of random walks in the plane $(\delta,S)$ that are absorbed by the collapse overdensity $\delta_{c}$ on masses larger than $M$. For Gaussian random fields this has a simple solution given by the Press-Schechter mass function [29]:

[TABLE]

where $\nu_{h}\equiv\delta_{c}/\sqrt{S}$.

The qualitatively-useful Press-Schechter distribution of [29] has long been supplanted—at least in general relativity—by a variety of empirical fits to N-body simulations (e.g. see [30, 31, 32, 33]), often preserving the assumption of universality of the collapsed mass fraction. In this context, expressing $\mathrm{d}F(>M)/\mathrm{d}\ln\nu$ in units of $\nu_{h}$ makes the resulting expression independent of the choice of cosmological parameters [34]. In this paper, we adopt the Peacock mass function222This simple functional form matches the Warren et al. fit [31] to the mass function to sub-percent accuracy. [35]

[TABLE]

Instead of the collapsed mass fraction (Eq. 3.3), we will use its differential, known as the multiplicity function333Note that some authors define the multiplicity function as the derivative with respect to $\sigma$ rather than $\ln M$. Moreover, although $f(M)$ and $\raisebox{-0.86108pt}{\text{{\large{f}}}}{}(S)$ are related, they do not denote the same functions, despite the same letter being conventionally used for both.:

[TABLE]

where we have defined the first-crossing distribution $\raisebox{-0.86108pt}{\text{{\large{f}}}}(S)$, the probability that the random walk will be absorbed by the barrier at resolution $S$. We will see in the next section that this becomes much more important in modified gravity than general relativity.

We finish this section by noting that, in a strict sense, when relating the excursion-set predictions to simulated or real data, the scale of the environment $S_{\rm env}$ should correspond to its Lagrangian scale (i.e. in the initial conditions) instead of its Eulerian size, and as pointed out by [36, 37] these effects could be relevant for modified-gravity theories. However, for the environment scales studied in this work ($R_{\rm env}\geq 10\,{\rm Mpc}/h$), we expect these two quantities to be very similar, the evolution of the density field since on such large scales is close to self-similar. We can verify this with a quick test: the effect of the Eulerian evolution can be roughly captured by making a log-normal transformation on the linear (Lagrangian) density field [38]. Then, since the shape of density profile around any point in either the Lagrangian or Eulerian fields is, on average, determined by the correlation function of the field, we can explore the effects of this distinction by studying the correlation function in either case. Under a lognormal transformation, the correlation function of the log-normalized (i.e. Eulerian) field would be given by $\xi_{\rm Eu}(r)=\exp[\xi_{\rm Lag}(r)]-1$, where $\xi_{\rm Lag}$ is the correlation function of the linearized (i.e. Lagrangian) field. Figure 1 shows both correlation functions for the smallest smoothing scale explored in this work ($R=10\,{\rm Mpc}/h$), and shows the negligible effect of ignoring the true Eulerian scale of the environment. This completes the excursion set theory in general relativity.

3.2 The excursion set theory in $f(\mathcal{R})$

We now summarise the differences between excursion set theory in $f(\mathcal{R})$ gravity compared to $\Lambda$CDM. This section presents a formula for the overdensity required for collapse $\delta_{c}(M,z,f_{\mathcal{R}0},\delta_{\rm env})$. The main difference in modified gravity is that $\delta_{c}$ depends both on mass, environment density and redshift apposed to just depending on redshift as in $\Lambda$CDM. We then discuss the consequences for calculating the first-crossing distribution $f(S)$ and state which of these we will include in our model.

The collapse of the environment surrounding the halo is equivalent to collapse in general relativity. Let us assume that the initial overdensity is a spherical top-hat in Eulerian space. We can utilise the resulting axisymmetry to simplify the gravitational collapse equation to:

[TABLE]

where $y(a)$ is the ratio of the physical radius of the halo $R_{\rm TH}(a)$ to the physical radius of the filter $a(t)R$ [15]. The initial conditions are obtained by mass conservation, which forces the equality [13]:

[TABLE]

An important corollary of Equation 3.6 is that $\delta=y^{-3}-1$, which contributes the nonlinearity in the final term of Equation 3.5. Now we have an expression for the halo density in general relativity and the environment density in modified gravity.

It is possible to show that, in $f(\mathcal{R})$, Eq. 3.5 takes the form:

[TABLE]

where $F_{\mathrm{eff}}=F_{\mathrm{eff}}(a,R_{\rm TH},\rho_{\rm TH},\rho_{\rm env})$ was introduced in Equation 2.11, with $R_{\rm TH}=yaR$, $\rho_{\rm TH}=\rho_{m0}a^{-3}y^{-3}$. Since $F_{\mathrm{eff}}\ll 1$ at early times, we can use the same initial conditions as we used for $\Lambda$CDM. This is the general form for modified gravity collapse of a spherical top-hat. Substituting $F_{\mathrm{eff}}$ from Equation 2.11 determines the ode for collapse of the halo in modified gravity.

How does this affect the excursion set formalism? After solving Eqs. Equations 3.5 and 3.7 we can derive values for $\delta_{c}$ in general relativity and modified gravity, and for the starting point of the random walk $\delta_{\textrm{env}}$. So far we have neglected the additional complications that the resolution $S$ is sensitive to $f(\mathcal{R})$ via the power spectrum and that the collapse density is now a function of $S$. We address these now.

We have chosen to use the linear $P(k)$ computed in the $\Lambda$CDM model in Equation 3.1 when computing the $f(\mathcal{R})$ gravity predictions. The same approach was taken in [39]. The reason for this choice is that linear theory massively overestimates the clustering in modified gravity theories such as $f(\mathcal{R})$: Figure 3 shows that the ratio of the linear power-spectrum in $\left|{f_{\mathcal{R}0}}\right|=10^{-5}$ to the linear power-spectrum in $\Lambda$CDM shows a much larger deviation than the ratio of the corresponding non-linear power spectra obtained from simulations. The reason is that screening is not present in linear theory. One can easily construct modified gravity models of this type where non-linear results for the power spectrum and the mass function can be arbitrarely close to $\Lambda$CDM while at the same time having a linear power-spectrum that deviates significantly from $\Lambda$CDM. In such scenarios excursion set theory is not able to give results that agree with simulations if the modified gravity linear power spectrum is used in the calculations. We have checked that this gives result that are in better agreement with simulations by comparing the results using the full linear power-spectrum. An alternative would be to consider a compromise approach where we use the linear $\Lambda$CDM power-spectrum corrected with a boost-factor $P_{f(\mathcal{R})}(k)/P_{f(\Lambda\rm CDM}(k)$ computed from simulations, however this would require explicit simulations which goes against the appeal of using the excursion set approach—which is to extract observables without having to perform expensive numerical calculations.

The barrier density is no longer flat: $\delta_{c}(S,\delta_{\textrm{env}})$ is a montotonically-decreasing (-increasing) function of $S$ ($M$; Figure 2a), and the peak-backgound split $\delta_{c}-\delta_{\textrm{env}}$ is a montotonically-decreasing function of $\delta_{\textrm{env}}$(Figure 2b). Compared to $\Lambda$CDM, we expect haloes to form at higher masses (earlier in the random walk) and to have more in low-density regions (where the fifth-force potential is unscreened). For a moving barrier like this, and under the assumption that trajectories in $(\delta,S)$ take uncorrelated steps, the diffuson equation admits a solution [40]:

[TABLE]

It is straightforward to show that this integral equation has an analytical solution for a constant barrier [41, 40], including the constant case $\delta_{c}^{\Lambda}\approx 1.676$ which is the $\Lambda$CDM solution. Once the first crossing distribution is found after solving this equation, the procedure would be to marginalise over the possible values of the environmental density. This procedure is numerically complicated and, since its derivation is based on purely uncorrelated random walks, which do not apply to the environment definitions used in this work, we took a simpler approach, based on including the effects of $f(\mathcal{R})$ simply by substituting the expression for the mass- and environment-dependent collapse threshold $\delta_{c}(S,\delta_{\textrm{env}})$ in the $\Lambda$CDM result described in the next section.

3.3 Tidal classification of the cosmic web

Now we build on the unconditional Press-Schechter result in Section 3.1 to find the equivalent result dependent upon the local environment. As done in e.g. [17, 18], we characterise the environment in terms of the properties of the local tidal tensor $\mathsf{T}_{ij}$, defined as the the Hessian of the normalised, Newtonian potential $\phi$ smoothed with a given kernel of size $R_{\rm env}$

[TABLE]

The environment is then classified into one of four different classes, depending on the values of the three eigenvalues of ${\sf T}$, $\lambda_{1}\geq\lambda_{2}\geq\lambda_{3}$, in relation to a given eigenvalue threshold. We thus define voids, sheets, filaments and knots as regions in which $\alpha\in\left\{0,1,2,3\right\}$ eigenvalues are above the threshold respectively. Since the tidal tensor quantifies the direction and intensity of the local tidal forces, this classification thus informs about the number of dimensions in which extended objects are contracted or stretched. We will call these environment types “elements” of the cosmic web.

For convenience, let $S_{\rm env}$ be the variance of the overdensity field convolved with the kernel defining the environment, and let us introduce the three useful quantities:

[TABLE]

The classification scheme outlined above can be cast in terms of the values of these quantities as:

[TABLE]

where $\lambda_{\text{th}}$ is the eigenvalue threshold, and $\alpha$ is the number of eigenvalues above it.

It can be proven that the probability distribution for the environmental tidal field in terms of $\nu_{e},\,\rho$ and $\theta$ is given by

[TABLE]

after integrating over the irrelevant rotation angles that diagonalize the tidal tensor. Thus we can estimate the volume fractions taken up by each element of the cosmic web as:

[TABLE]

Now we use the Gaussian statistics to find the conditional Press-Schechter Ansatz. The simplification that we are only interested in spherical collapse has the corollary that the distribution of the overdensity field conditional on environment only depends explicitly on the value of the environmental overdensity $\delta_{\textrm{env}}=\nu_{e}\sqrt{S_{e}}$:

[TABLE]

where $\nu_{h}=\delta_{c}(S,\delta_{\textrm{env}})/\sqrt{S}$ relates the halo mass $M$ to the corresponding value of the overdensity field $\delta_{h}$, with $\nu_{e}=\delta_{\textrm{env}}/\sqrt{S_{\textrm{env}}}$ the analagous quantity for the environment and $\epsilon^{2}=S_{eh}/\left(S_{\textrm{env}}S\right)$ is the correlation coefficient between $\delta_{h}$ and $\delta_{\textrm{env}}$.

We are now in a position to present the key integral of this paper. First we must take into account that we require $S>S_{\textrm{env}}$ for the excursion set theory to work. This corresponds to the physical requirement that the halo be embedded in the environment and not vice-versa. Under the simplification that the moving barrier $\delta_{c}$ does not affect Equation 3.14, we find that the conditional collapsed mass fraction is that same as the result in Equation 3.3, subject to the change of variables from $\nu$ to $\nu_{\mathrm{eff}}$ (defined below). Thus we find the corollary that the first-crossing distribution becomes:

[TABLE]

This is the first-crossing distribution which we substitute into Equation 3.4 to obtain the $f(\mathcal{R})$ multiplicity function conditional to the tidal classification of the environment.

4 Simulations and algorithms

The simulations used in this paper were performed with the ISIS code [42], which is a modified gravity modification of the N-body code RAMSES [43].

We performed two large-box simulations: one for $\Lambda$CDM and one for $\left|{f_{\mathcal{R}0}}\right|=10^{-5}$ using an approximate method [44] to incorporate the scalar field. We also performed two small-box simulations: one for $\Lambda$CDM and one for $\left|{f_{\mathcal{R}0}}\right|=10^{-5}$ fully solving for the scalar field in the simulation box [45]. The approximate method has previously been shown to produce very good results (accuracy to a few percent) with respect to the matter power spectrum and the halo mass function [44]. The modified gravity simulations were performed using the same inital conditions as the corresponding $\Lambda$CDM simulations. This allows us to compare the two models without requiring multiple realisations of the intial conditions. The cosmological parameters are: $\Omega_{m0}=0.27$, $\Omega_{\Lambda 0}=0.73$, $h=0.704$, $n_{s}=0.966$ and $\sigma_{8}=0.8$ and the other simulation parameters are listed in Table 1.

We located haloes using a Friend-of-Friend algorithm444https://github.com/damonge/MatchMaker as well as the spherical overdensity halo finder AHF (Amiga Halo Finder)555http://popia.ft.uam.es/AHF/Download.html [46]. The N-body particles were binned to a $N_{\mathrm{grid}}=512^{3}$ grid using Cloud-in-Cell interpolation and then smoothed with a Gaussian kernel of widths $R_{\rm env}=10$ and $20$ Mpc$/h$ in this paper. The local tidal tensor and its eigenvalues were then computed using a public code666https://github.com/damonge/DensTools. The error bars in the plots were computed using jack-knife resampling of the simulation box with 8 samples. To generate error-bars for ratio plots coming from two different simulations, we used:

[TABLE]

where $r$ is the correlation coefficient. The quantities whose ratios we compute are highly correlated (since they are taken from simulations with the same initial conditions) so for simplicy we have simply used a fixed value of $r=0.9$ (or $r=0$ for cases where we wanted to be conservative) to estimate the error-bars.

5 Results

We compare the behaviour of the $\Lambda$CDM and $f(\mathcal{R})$ mass functions in various environments. Then we evaluate the accuracy of our semi-analytic model using N-body simulations and describe its limitations and possibles paths to improvement. Finally, we analyze the effect of modified gravity on a number of internal halo properties as a function of tidal environment.

5.1 Volume fractions

We can indirectly measure the accuracy of the Gaussian model for the eigenvalue probability distribution via the volume fractions in different cosmic web elements. In Figure 4 we compare the values from our model to those from the N-body simulations.

The behaviour with $\lambda_{\mathrm{th}}$ permits us to select a suitable value for $\nu_{\mathrm{th}}$ in Equation 3.11. We require a value for which there is a (relatively) even distribution between the different morphologies (amongst other requirements discussed in [18]). Thus we select $\lambda_{\mathrm{th}}=0$ for $R_{\textrm{env}}=20$ Mpc$/h$ and $\lambda_{\mathrm{th}}=0.1$ for $R_{\textrm{env}}=10$ Mpc$/h$. We utilise the same threshold for the smoothing in our semi-analytic model and the N-body simulations.

We compare the volume fractions calculated from Equation 3.13 to the N-body volume fraction in Figure 4. At scales of $R_{\rm env}=20\,{\rm Mpc}/h$ Equation 3.12 performs well, whereas at $R_{\rm env}=10\,{\rm Mpc}/h$ we find that $p(\rho,\theta,\phi)$ is a poor approximation of the real overdensity field. This should not be surprising given the non-linear evolution of the density field on all but the largest scales.

The Gaussian volume fraction in Figure 4 illustrates the large-scale behaviour of the scalar field modification to general relativity. Compared to $\Lambda$CDM (solid lines in Figure 4) there is a greater fraction of voids and a smaller fraction in collapsed structures. At larger scales, the change in volume fraction due to modified gravity is more pronounced. This is due to the scalar field leaving the chameleon regime at small scales and moving to the linear regime at large scales (although on cosmological scales it returns to mimic general relativity behaviour). More compact morphologies exhibit a greater deficiency, so that the overall fraction remains at unity. However, the difference between even $\left|{f_{\mathcal{R}0}}\right|=10^{-5}$ and $\Lambda$CDM is dwarfed by the difference between the theoretical volume fraction and the N-body result at all smoothing scales. This provides a quantitative basis for the assertion in [13] that there is insufficient evidence to warrant modelling the environemnt density using the $f(\mathcal{R})$ equations.

5.2 Multiplicity functions in each environment

We compare the behaviour of the $\Lambda$CDM mass functions to $f(\mathcal{R})$ with $\left|{f_{\mathcal{R}0}}\right|=10^{-5}$ in the different structures of the cosmic web and averaged over all environments. Our semi-analytic model performs best at successively larger scales and in lower-density environments, but reproduces the main trend found in the simulations: the $f(\mathcal{R})$ abundances are amplified in underdense environments and at low masses.

The smoothing scale determines the range of masses in which our results are reliable. Recall from Section 3.1 that the major shortcoming of the excursion set formalism is that whenever $S<S_{\textrm{env}}$ the process is not well defined, since the condition $S<S_{\textrm{env}}$ implies that the halo has already collapsed at resolution $S$, before the starting points of its trajectory at $S_{\textrm{env}}$. Therefore care must be taken when interpreting any results for masses with $S\sim S_{\textrm{env}}$. For the filter sizes used in this work, $R=10\,{\rm Mpc}/h$ and $R=20\,{\rm Mpc}/h$, the corresponding masses enclosed within them are $M_{10}=1.67\cdot 10^{15}\,M_{\odot}/h$ and $M_{20}=1.34\cdot 10^{16}\,M_{\odot}/h$, which act as upper bounds for the mass range where the excursion set predictions are meaningful.

A desirable property of any model of the mass function that incorporates environmental dependence is that, when marginalised over all possible environments, one should recover the unconditional mass function given by the same model. Currently our model is based on simply substituting the argument of the unconditional mass function from $\nu_{h}$ to $\nu_{\rm eff}$, inspired by the equivalent result found in the excursion-set formalism. However, this rescaling of the mass function is only mathematically consistent, in the sense described above, for the specific form of the Press-Schechter mass function [40]. Thus, most simple attempts at modelling the conditional mass function have only been able to provide a qualitative description of it, with a poor quantitative performance. This problem is also illustrated in Figure 5, which shows the mass function from our simulations together with the Peacock unconditional mass function and the excursion set prediction after rescaling the Peacock mass function and marginalising it over the environment. Our results agree with [40], who showed that the only solution of Equation 3.8a which matches the unconditional case is Press-Schechter with a $\delta_{c}$ linear in $S$. Thus we do not expect the cosmic-web-sensitve result to match the percent-level fits to N-body which motivated the unconditional form of the Peacock mass function.

The additional dependence of the collapse density in $f(\mathcal{R})$ amplifies this same problem: compare Figures 6a, 8a, 7a and 9a to Figures 6b, 8b, 7b and 9b, which show the multiplicity functions in different environments measured from the simulations together with the predictions of our method for a range of smoothing kernels (with $R_{\rm env}=10$ and 20 ${\rm Mpc}/h$ in Figures 8 and 9 and Figures 6 and 7 respectively, and using Gaussian and sharp-$k$ filters in Figures 8 and 6 and Figures 9 and 7 respectively). We concur with [13], who state that the environment-averaged result simplifies to the result with the random walk starting at $\left(\delta_{\textrm{env}}=0,S_{\textrm{env}}=0\right)$ (i.e. to the unconditional mass function) only in $\Lambda$CDM and not in $f(\mathcal{R})$. Thus we expect our fits to be less accurate (when compared to N-body results) in $f(\mathcal{R})$ than in $\Lambda$CDM, a hypothesis borne out by our results.

The trends with environment smoothing scale and density shown in Figures 6, 7, 8 and 9 can be understood in terms of deviations of the non-linear environment distribution with respect to the Gaussian prediction. Having seen in Figure 4 that the $R_{\rm env}=20\,{\rm Mpc}/h$ volume fraction is well-approximated by the integral of $p(\rho,\theta,\nu_{e})$, the density field on such large scales in still evolving linearly from the primordial Gaussian density field. Thus in Figures 6 and 7, we find good agreement in both $\Lambda$CDM and $f(\mathcal{R})$, regardless of halo mass or the location within the cosmic web. On environment scales of $R_{\rm env}=10\,{\rm Mpc}/h$, the results in Section 5.1 suggest that the density field evolution is weakly-linear. We expect deviations from linear collapse to occur more quickly in increasingly dense environments, where the local structures have already collapsed along one (sheets), two (filaments) or all three (knots) axes. In Figures 8 and 9 our predictions lie within the N-body error bars at progressively fewer points for voids, sheets, filaments and knots. These trends occur independent of cosmology, which futher reinforces the hypothesis that the Gaussian model for the density is responsible, as we model the large-scale structure according to general relativity in both $\Lambda$CDM and modified gravity.

Like [18], we find that the excursion-set prediction for the conditional mass function in $\Lambda$CDM becomes a worse model for the N-body data in higher-density environments, underpredicting the number of low-mass halos and over-predicting the abundance in the high-mass tail at $R_{\rm env}=10\,{\rm Mpc}/h$ (see Figures 8a and 8b), and this is true also for $f(\mathcal{R})$. For a smoothing of $20\,{\rm Mpc}/h$, Figures 6a and 6b show that the prediction becomes a much better fit, with a slight deficiency of haloes at high masses in knots.

Despite the shortcomings of the excursion set formalism when making accurate predictions for the conditional mass function, we can still use this framework to understand the differences between the $\Lambda$CDM and $f(\mathcal{R})$ predictions. The $f(\mathcal{R})$ mass functions are amplified compared to $\Lambda$CDM. The key trends visible in Figures 6c, 8c, 9c and 7c are that the amplification increases with lower environment density. The amplification tends to increase with halo mass in voids and sheets, while there is a maximum enhancement mass in filaments and knots. This behaviour can be explained by the peak-background split. For equivalent values of $S$, the barrier density is lower for $\left|{f_{\mathcal{R}0}}\right|=10^{-5}$ than $\Lambda$CDM, which makes it more probable for the first-crossing of the excursion set trajectory to occur. This effect is amplified by a larger value of $\delta_{c}(\delta_{\textrm{env}},S)-\delta_{\textrm{env}}$, i.e. in voids and sheets, but damped in dense regions where $\delta_{c}(\delta_{\textrm{env}},S)-\delta_{\textrm{env}}\rightarrow 0$ via screening. Thus, the modified gravity behaviour compared to $\Lambda$CDM is driven by screening at all scales.

In order to improve the semi-analytical model for the mass function used here, three problems need to be overcome:

• •

The Gaussian eigenvalue distribution $p(\rho,\theta,\nu_{e})$ is a poor description on small smoothing scales.

• •

At resolutions smaller than $S_{\textrm{env}}$, the excursion set theory which underpins the first-crossing distribution is undefined.

• •

The $\nu\rightarrow\nu_{\rm eff}$ rescaling is known to be a poor approximation to the conditional mass function even in $\Lambda$CDM.

The first issue can be resolved in three ways. Following [18], we could replace the functional form with the distribution of values from the N-body simulations. We find this solution unsatisfactory, due to the requirement to run many realisations to avoid cosmic variance and the dependence upon individual simulation parameters. However, as [18] note, this remains a useful technique to separate the effects of the Gaussian density field from the other approximations. Alternatively, it may be possible to add small-scale corrections to the probability density function, or to replace the Gaussian density field of [47] with the log-normal density field proposed by [38].

The second item appears to be insurmountable. This is because $M(S_{\textrm{env}})$ is completely controlled by our smoothing scale $R_{\textrm{env}}$ via the window function. Fortunately, on mildly non-linear scales and larger, there are so few haloes with $M>M(S_{\textrm{env}})$ that the problem is only an issue in mass bins which are already uncertain due to small-number statistics and cosmic variance. Resolving these problems will improve the quantitative accuracy of our model beyond the large-scale, low-density regime. [18] found a better fit to the data by using an “effective-universe” approach, where halo abundances in an environment with overdensity $\delta_{\textrm{env}}$ are predicted as the unconditional abundance in a universe with effective cosmological parameters governed by the value of $\delta_{\textrm{env}}$, and the $S_{\textrm{env}}$ mass cut is accounted for in an ad-hoc manner by limiting the Fourier modes available in that effective universe to those smaller than $\sim R_{\rm env}$. While there remain quantitative issues to overcome, our extremely simple model performs well qualitatively on all scales and in all density environments.

The last item is not readily solved, and is a common limitation of the conditional mass function. In order to produce an environment definition more akin to the excursion-set formalism we have also studied cases where the density field is filtered using a sharp-$k$ window function. As shown in Figures 9 and 7, this produces a better agreement with the N-body simulations. A possible improvement may be to use an alternative to excursion set, e.g. the survival probability approach of [48] or the Markovian Velocity excursion set of [49].

We have seen in the preceding two sections that the model illustrated in this paper is qualitatively useful but not always quantitatively accurate. At large smoothing scales, where we are in the linear regime, our results agree well with the full N-body calculations. As we move to smaller scales, we encounter the increasing non-Gaussianity of the environemnt overdensity and the chameleon regime of the $f(\mathcal{R})$ modification. While we can extend our model into this regime in low-density environments, it performs badly in knots, and by non-linear scales we only reproduce qualitative behaviour. Nevertheless, we can attribute the difficulties to certain approximations and assumptions in our model, which in turn suggests avenues for improving our results.

5.3 Multiplicity functions at fixed environment density

One of the key predictions of the excursion-set formalism in $\Lambda$CDM cosmologies with Gaussian initial conditions is that halo abundances only depend on the density of the environment, i.e. not on any other local environmental parameters. In spite of the approximate nature of this formalism, this result was found to be true in [18] in comparison with N-body data. Therefore, it is interesting to explore whether the same is true for modified gravity theories.

Figure 10 shows the conditional mass function for fixed environmental densities in different types of environments defined in terms of the eigenvalues of the tidal tensor. We do not find significant differences between environments in terms of halo abundances. This is also a prediction of our semi-analytic model: we treat the environment as if it were $\Lambda$CDM (rather than with a full $f(\mathcal{R})$ treatment) not only in calculating the environment density ODE, an approximation justified in [13]; but also (via the power spectrum) the environment variance. Furthermore, the only change from $\Lambda$CDM to $f(\mathcal{R})$ in our model is the barrier density, which itself is only a function of $\delta_{\textrm{env}}$ and not of $(\rho,\theta)$. This is an a priori simplification which is a corollary of only examining the modified gravity collapse assuming spherical symmetry, instead of the full ellipsoidal collapse. Nevertheless, as in the case of $\Lambda$CDM, we find that this phenomenological result is in agreement with the simulated data. It could be that $f(\mathcal{R})$ does modify the $\Lambda$CDM result in certain regimes, but $R_{\textrm{env}}=10/h\textrm{ Mpc}$ is a sufficiently large scale that the modifications are negligible, i.e. the $f(\mathcal{R})$ modification is in the cosmological regime where it must mimic the $\Lambda$CDM background. Alternatively, perhaps the influence of $f(\mathcal{R})$ is minimal on the structure of the cosmic web (e.g. by late times non-linearities in the density field or the merger history of haloes may overwhelm any contribution from modified gravity) so the same physics dominates in both cosmologies. Regardless, the fact that we encounter the same result in both the semi-analytical and the numerical solution in encouraging.

5.4 Internal halo properties

The relation of certain internal halo properties with the tidal structure of the environment could also potentially contain relevant signatures for modified gravity models. To this end, we have studied the relative orientation of the halo shapes and halo spins with respect to the eigenvectors of the local tidal tensor, as well as the distribution of the halo spin parameter and the concentration-mass relation as a function of environment type. In order to obtain reliable estimates of these internal properties, we selected a subset of the halo catalog containing only virialized structures (as done in [42]) with more than 100 particles per halo. While [50] suggest a cut-off of 300 particles per halo, on the basis of ensuring that the resolution bias is less than 10% compared to the ellipticity result from sampling a halo of 1000 particles, our smaller cut-off produces a maximum resolution bias of $\sim 60\%$, which declines to $\sim 30\%$ as the ellipticity increases. This affects the uncertainty on the shape of the haloes: thus, any difference between $\Lambda$CDM and $f(\mathcal{R})$ would have to exceed this uncertainty to be a non-null result. We select only virialised structures by combining the kinetic energy $T$, gravitational potential energy $W$ and surface pressure term $E_{S}$ of the halo to define a parameter $\eta$ such that:

[TABLE]

The first equality is the virial theorem, which defines a virialised structure by virtue of having reached a state of equilibrium [51]. The possibly non-zero value of $\eta$ allows for uncertainty in the measurement of the kinetic energy of the halo, so we set haloes with $\eta<0.2$ to be virialised in accordance with the literature [42]. The decision to remove non-virialised haloes, rather than keeping all gravitationally-bound structures, stems from the different distributions of halo properties within the two classes of haloes [51]. We compute the relative orientation of the halo shape by estimating the inertia tensor of the halo as [52, 42]:

[TABLE]

where the position of the $n$-th halo is $x_{n}$ along each axis $i,j\in[1,3]$ and the sum is over all particles in the halo. The principal axis of the halo is then defined as the eigenvector corresponding to the largest eigenvalue of $\hat{I}$, the (reduced) inertia tensor. For each halo we then compute $\mu_{11}$ and $\mu_{13}$, the cosine of the angle between its principal axis and ${\bf e}_{1}$, and ${\bf e}_{3}$, the eigenvectors of the local tidal field (defined with a smoothing scale $R_{\rm env}=10\,{\rm Mpc}/h$) corresponding to the largest and smallest eigenvalues respectively. Figure 11 shows the average $\mu_{11}$ and $\mu_{13}$ as a function of halo mass in the four different environment types. We find no significant difference between $\Lambda$CDM and $f(\mathcal{R})$ in the relative halo orientations with respect to the cosmic web. We also performed a similar analysis studying the relative orientation between the principal directions of the cosmic web and the halo spin. We found no significant difference between $f(\mathcal{R})$ and $\Lambda$CDM in this case either. Given the interpretation of [53] that the alignment between dark matter haloes and the cosmic web is driven by the dynamics of cosmic flows during the merger history of the haloes, it is unlikely that this behaviour should be altered by modifications to gravity. Specifically, the presence of the scalar field will affect the gravitational interaction during halo mergers, which are the key contributors to the final orientation and spin of the haloes. However, if haloes flow preferentially towards increasingly collapsed structures, ([53]; e.g. along filaments towards knots), then the merger history of these haloes will occur in increasingly screened areas. If the halo formed at early times (even in an unscreened void), each additional contribution to its orientation and/or spin vector will tend towards the $\Lambda$CDM equivalent, washing out any initial $f(\mathcal{R})$ modification to the large-mass haloes. In $\Lambda$CDM, [53] find that haloes above a threshold of $M>(8\pm 2)\cdot 10^{12}M_{\odot{}}$ exhibit this large-mass behaviour, which is the range we have tested here. Thus our null result is consistent with current interpretations of dark matter halo formation.

It is also customary to define the halo spin parameter [54]

[TABLE]

where ${\bf J}$ and $R$ are the halo angular momentum and virial radius respectively. The distribution of $\lambda$ is known to be only mildly dependent on mass777This dependence is described by [55] as “small but real” but the precise nature of it is so dependent upon the halo finder used that they do not quanitify the dependence for the FOF finder used in this paper. and redshift in $\Lambda$CDM [55, 56], and therefore it is interesting to explore departures from this independence in modified gravity theories. Figure 12c shows the spin distribution for all halos in the $\Lambda$CDM and $f(\mathcal{R})$ simulations and Figure 12d the ratio of the distributions between both cosmologies for the four different environment types. The spin-parameter is seen to be higher in modified gravity which is in agreement with what was previously found in [57]. We also split the halo catalog into two samples $M<3\cdot 10^{13}M_{\odot}$ and $M>3\cdot 10^{13}M_{\odot}$ and computed the PDF in each of the samples. For the high-mass sample the PDF is very close to $\Lambda$CDM while for the low-mass sample the PDF is enhanced in $f(\mathcal{R})$ as shown in Figure 12. In this figure we also show how the result changes if we only include halos with more than $300$ particles ($M>2.2\cdot M_{\odot}/h$).

This agrees with [57] where the conclusion (for the same $f(\mathcal{R})$ model, but with $\left|{f_{\mathcal{R}0}}\right|=10^{-6}$) was that modified gravity spins up galactic-sized halos. We can intuitively understand why spin has to be modified via Equation 5.3. Let us approximate the spin as $\left|{\mathbf{J}}\right|=MR\left|{\mathbf{v}}\right|$. Recalling Equation 2.10, we see that $G$ is enhanced by the factor $(1+F_{\mathrm{eff}})$ with $F_{\mathrm{eff}}\in\left[0,\frac{1}{3}\right]$. Thus we indirectly increase $\left|{\mathbf{v}}\right|$, with $v\propto G$ in the linear regime, leading to the approximation $\lambda\propto\sqrt{(1+F_{\mathrm{eff}})G}$. In the absence of screening we would expect $\lambda^{f(\mathcal{R})}/\lambda^{\Lambda\mathrm{CDM}}\approx\sqrt{1+\frac{4}{9}}$, a factor of $1.2$. In practice, this is an upper bound because we have neglected the non-linear regime (in which $v\propto G$ does not hold) and the distribution of partial screening. The shift in the PDF we find is roughly $\lambda^{f(\mathcal{R})}/\lambda^{\Lambda\mathrm{CDM}}\approx 1.04$ which is within our upper bound. We also see a clear dependence of the spin-parameter in the cosmic web. Figure 12d demonstrates that the spin parameter is boosted mostly in low density environments (in voids and sheets) whereas in high-density environments (in filaments and knots) the value is close to $\Lambda$CDM for all values of $\lambda$. This is as expected due to the environmental dependence of screening: a halo that is not screened if placed by itself might still be screened if placed in a high density environment. Given the density ranges in Equation 3.11, the trend of enhancement with tidal morphology is what we expect from theory.

In Figure 13 we see the dependence on the NFW concentration-mass parameter with mass and environment. The concentration-mass is enhanced in our $f(\mathcal{R})$ model for halos with mass in the range $M\lesssim 10^{14}M_{\odot}/h$, whereas for halos of larger mass the results are close to $\Lambda$CDM. This agrees with the results of [58], which only considered cluster-sized halos. There is some dependence on the cosmic web, but it is insignificant in comparison with the scatter in each mass-bin.

6 Discussion

We study the properties of dark matter haloes in the context of the environmental tidal classification of the cosmic web in $f(\mathcal{R})$ theories. This classification defines four different types of environments based on the directionality of the tidal forces: voids are defined as regions of space where tidal forces will expand an extended object in all directions, sheets and filaments have instead one or two compressing directions respectively, and in knots tidal forces will only compress structures. The relation between the gravitational potential and the density field also implies that these environment types will sample partially overlapping ranges of densities, and thus this classification can be a useful tool to study the effect of the chamaleon screening, present in viable $f(\mathcal{R}{})$ theories, on the properties of dark matter haloes.

We have described an approximate method to predict the abundance of dark matter haloes in each environment type. Recall that our purpose was to use excursion set theory to take into account as much of the structure of the cosmic web and the non-linear collapse in $f(R)$ as possible in a semi-analytic model, in order to obtain results that can fit numerical simulations. Thus, we incorporated both the moving barrier of the excursion-set formalism—inherent in screened models—and the analytical description of the mass function conditional upon tidal environment described in [18]. When comparing this prediction with the data from simulations we find that, even though the method is able to reproduce the correct behaviour of modified gravity theories qualitatively, the predictions become unreliable in high-density environment and for mildly non-linear filter scales. The method is, however, able to predict the enhanced abundance of high-mass objects found in low-density environments, a direct consequence of chamaleon screening, and the prediction is accurate in voids, where the departure from $\Lambda$CDM is largest. The shortcomings of this model are directly related to the inability of the traditional excursion-set formalism to accurately describe the conditional mass function, and we have outlined a number of ways forward to improve the predictions.

In spite of these shortcomings, we have gained further insight into the possible impact of modified gravity on halo properties conditional on environment by directly comparing the results of $f(\mathcal{R})$ and $\Lambda$CDM simulations:

• •

One of the key predictions in $\Lambda$CDM is that halo abundances should only depend on the enviromental density, and not on the directionality of the tidal forces. Our results show no significant departure from this behaviour in $f(\mathcal{R})$, at least on mildly non-linear filter scales.

• •

We have identified a number of internal halo properties that are not affected by modified gravity in this context. In particular, we find that the alignment of halo orientations and spins with respect to the principal tidal directions of the environment.

• •

We have also seen that, even though the mass-concentration relation is different in $f(\mathcal{R})$ and $\Lambda$CDM, its dependence on the tidal environment only shows marginally significant differences between both models.

• •

We observe a dependence in the distribution of the halo spin parameter upon tidal environment, with stronger deviations from the $\Lambda$CDM results in voids and sheets, as expected in the presence of chamaleon screening.

Incorporating morphological classification of the cosmic web into the existing analysis of dark matter halo properties increases their utility as a probe of modified gravity. By doing so we can distinguish between the collapsed mass fraction in $f(\mathcal{R})$ and $\Lambda$CDM far more than merely using the unconditional result, especially in low-density environments. Furthermore, we can draw a boundary between the halo properties that are affected by screening and those which are not by explicitly examining the results in screened and unscreened regions, rather than averaging the effect over the distribution of environmental properties. While we find three null results, we also have three avenues for detecting the presence of $f(\mathcal{R}{})$ theories which are designed to evade local tests.

Acknowledgements

We would like to thank Pedro Ferreira for many useful comments and discussions. HAW and DA are supported by the Beecroft Trust and the Oxford Martin School.