Marked clustering statistics in $f(R)$ gravity cosmologies
César Hernández-Aguayo Carlton M. Baugh Baojiu Li
astro-ph.CO


firstpage–lastpage
Precision mass measurements on neutron-rich rare-earth isotopes at JYFLTRAP - reduced neutron pairing and implications for the r-process calculations
    J. Äystö
    December 30, 2023
=====================================================================================================================================================





We analyse the two-point and marked correlation functions of haloes and galaxies in three variants of the chameleon f(R) gravity model using N-body simulations, and compare to a fiducial ΛCDM model based on general relativity (GR). Using a halo occupation distribution prescription (HOD) we populate dark matter haloes with galaxies, where the HOD parameters have been tuned such that the galaxy number densities and the real-space galaxy two-point correlation functions in the modified gravity models match those in GR to within 1∼3%. We test the idea that since the behaviour of gravity is dependent on environment, marked correlation functions may display a measurable difference between the models. For this we test marks based on the density field and the Newtonian gravitational potential. We find that the galaxy marked correlation function shows significant differences measured in different models on scales smaller than r≲ 20 h^-1 Mpc. Guided by simulations to identify a suitable mark, this approach could be used as a new probe of the accelerated expansion of the Universe. 





gravitation – cosmology: theory  – large-scale structure of Universe – methods: statistical – methods: data analysis












§ INTRODUCTION



The Λ cold dark matter (ΛCDM) model is currently the most widely accepted description of the Universe (e.g. see ). In this model small ripples in the density of the Universe at early times, seeded during a period of rapid expansion called inflation, were boosted by gravity to form the cosmic web of voids, galaxies and clusters of galaxies that we see today. 
The ΛCDM model works remarkably well on large scales, with the highlight being the prediction of the temperature fluctuations in the cosmic microwave background (CMB) radiation (e.g. ). However, the model arguably runs into difficulties on small scales, which could be related to the nature of the dark matter particle or could be solved by appealing to the physics of galaxy formation (for a review see ). 

Several large international survey projects are underway which aim to determine what is behind the accelerating expansion of the Universe, such as the Dark Energy Survey (DES) <cit.>, the Dark Energy Spectrographic Instrument (DESI) survey <cit.> and the European Space Agency's Euclid mission <cit.>. These surveys aim to make ambitious measurements which, for the first time, will be dominated by systematics rather than by the volume of the Universe surveyed. Accurate theoretical predictions are urgently required to compare with these expected measurements and to rule out competing scenarios for the accelerating cosmic expansion.

Theoretically, the ΛCDM model is somewhat unappealing due to the presence of the cosmological constant, the agent behind the accelerating cosmic expansion. The magnitude of the cosmological constant is hard to motivate from a particle physics perspective. There is also a “Why now?” problem: a strong coincidence seems to be required for us to be at the right point in cosmic history to experience comparable energy densities in matter and the cosmological constant, with the latter dominating the current expansion. As a result, alternatives to the cosmological constant have been studied extensively in recent years: one possibility is adding more matter species to the energy-momentum tensor (the so-called dark energy models; see e.g., ); on the other hand there are models that change the left-hand side of Einstein's equations (these models are called modified gravity models, for reviews see ). 

Here, we focus our attention on a particular class of modified gravity models – <cit.> chameleon f(R) gravity. This model is obtained by adding a general function of the Ricci scalar, f(R), to the Einstein-Hilbert action. This modification gives rise to a new scalar degree of freedom in gravity <cit.>. In order to recover GR, the new scalar field becomes massive in high-density regions (i.e., the Solar System) and its interactions are suppressed by the so-called chameleon screening mechanism <cit.>. 

The standard tool to model the growth of large scale structure into the non-linear regime is N-body simulation. In order to robustly test gravity on cosmological scales, reliable N-body simulations of modified gravity models are essential. The non-linear nature of the scalar field equation requires the implementation of novel numerical techniques, which is what makes N-body simulations of modified gravity challenging <cit.>. Once such simulations are ready, we can use measurements of clustering statistics from surveys to test and constrain cosmological models (see e.g., ).

Several recent works have studied clustering in f(R) gravity. For example, <cit.> predicted the matter and velocity divergence power spectra and their time evolution measured from several large-volume N-body simulations with varying box sizes and resolution. <cit.> predicted the clustering of dark matter in redshift space, finding significant deviations from the clustering signal in standard gravity, with an enhanced boost in  power on large scales and stronger damping on small scales in the f(R) models compared to GR at redshifts z<1. More recently, <cit.> compared the time evolution of the two-point correlation function of dark matter haloes in real and redshift space in modified gravity and GR. 

An approach related to the two-point correlation function has been proposed to test modified gravity models, called the marked  correlation function <cit.>.  Marked statistics offer the possibility of testing how galaxy properties correlate with environment. Previous applications of the marked correlation function range from the analysis of the environmental dependence of bars and bulges in disc galaxies (see e.g., ) to breaking degeneracies in halo occupation distribution modelling <cit.>.

<cit.> proposed that marked statistics might provide a means to distinguish between modified gravity models and GR, due to the environmental dependence of the strength of gravity (for earlier work  on environmental dependence in modified gravity, see, e.g., ). Marks can be designed which down-weight high-density regions, for which modifications to gravity are screened, and up-weight low-density, unscreened regions to maximise the differences in the clustering signal. <cit.> tested this idea by using the dark matter particle distribution from N-body simulations of symmetron and f(R) modified gravity models. In the case of f(R) gravity with |f_R0| = 10^-4, they found a maximum difference of 37% with respect to GR. <cit.> studied the galaxy marked correlation function by up-weighting low and high-density regions using marks in function of the galaxy density field and the host halo mass of galaxies, they found significant differences between the f(R) and GR models.  Here, we focus our attention on the clustering of dark matter haloes and HOD galaxies for f(R) gravity models to make a more direct connection with observations.  

This paper is organized as follows. In Section <ref> we give a brief review of the theoretical description of f(R) gravity and the physical motivation behind this model. Section <ref> explains the numerical set-up of the simulations and the generation of halo and galaxy catalogues. The main results are presented in Section <ref>. Finally, in Section <ref> we present a brief discussion and our general conclusions.




§ F(R) GRAVITY THEORY


A popular family of modified gravity models is obtained by replacing the Ricci scalar R in the usual Einstein-Hilbert Lagrangian density by some function f(R) (see 
for recent reviews). These models are often considered as an alternative solution for the accelerating cosmic expansion. However, it should be noted that they are able to accelerate the expansion only because a cosmological constant is added “through the back door” <cit.>, and so they do not offer a real solution to the cosmological constant problem. Nevertheless, they constitute a  representative model with which to study cosmological constraints on possible deviations from GR.




 §.§ Theoretical framework


The action of f(R) theories is given by <cit.>:

    S = ∫d^4x√(-g)1/16π G (R + f(R)) + S_ m(g_μν,ψ_i)  ,

where G is Newton's constant, g is the determinant of the metric g_μν and S_ m is the action of the matter fields ψ_i (including the contributions from cold dark matter, baryons, radiation and neutrinos). 

Varying the action given in Eqn. (<ref>), with respect to the metric, g_μν, one obtains the modified Einstein equations

    G_μν + f_R R_μν - g_μν(1/2f(R) -  f_R ) - ∇_μ∇_ν f_R = 8π G T^ m_μν ,

where G_μν = R_μν - 1/2g_μν R is the Einstein tensor, ∇_μ is the covariant
derivative, = ∇^μ∇_μ the d'Alambertian, and T^ m_μν is the energy-momentum tensor for matter.
Here

    f_R ≡df(R)/dR ,

is the extra degree of freedom of this model, known as the scalaron field.

Taking the trace of Eqn. (<ref>), we obtain the equation of motion for the scalaron field,

    f_R = 1/3 (R - f_R R + 2f(R) + 8π Gρ_ m) ,

where ρ_ m is the matter density of the Universe.

Since we are interested in the cosmological evolution of the model, we derive the perturbation equations using the flat Friedmann-Robertson-Walker (FRW) metric in the Newtonian gauge

    ds^2 = (1+2Ψ) dt^2 - a^2(t)(1-2Φ)γ_ijdx^i dx^j ,

where Ψ and Φ are the gravitational potentials, t is the cosmic time, x^i represent the comoving coordinates, γ_ij is the 3D metric, and a is the scale factor, with a(t_0) = a_0 = 1 at the present time.

In the quasi-static and weak field limits <cit.>, structure formation in this model is determined by the following equations:

    ∇^2 Φ = 16π G/3 a^2 δρ_ m + 1/6 a^2 δ R ,


    ∇^2 f_R = -a^2/3 [δ R + 8π Gδρ_ m] ,

in which ∇^2 is the three-dimensional Laplacian operator, 
δρ_ m = ρ_ m - ρ̅_ m and δ R = R(f_R) - R̅ are, 
respectively, the density and Ricci scalar perturbations (overbars denote 
background cosmological quantities). 



 §.§ The chameleon mechanism


In f(R) gravity the modifications to Newtonian gravity can be considered as a fifth force mediated by the scalaron field, f_R. The chameleon mechanism (see e.g., ) was introduced to give scalar fields an environment-dependent effective mass allowing the scalar mediated force to be suppressed under certain environmental conditions. 

Because the scalaron field itself is massive, this fifth force is of Yukawa type, i.e. decaying exponentially as exp (-m_eff r), where m_eff is the scalaron mass

    m^2_eff≡d^2V_eff/d f^2_R ,

and the effective potential is related to the trace of the modified Einstein Eqn. (<ref>)

    dV_eff/d f_R = 1/3 (R - f_R R + 2f(R) + 8π Gρ_ m) .

In f(R) gravity this fifth force enhances gravity in regions of weak gravitational potential and on scales below the Compton wavelength of the scalar field, λ = m_eff^-1.
While in high density regions, m_eff is heavy and the fifth force is more strongly suppressed and gravity reverts to GR.



 §.§ The Hu-Sawicki model


We consider the  model:

    f(R) = -m^2 c_1/c_2(-R/m^2)^n/(-R/m^2)^n + 1 ,

where m is a characteristic mass scale, defined by m^2 = 8π G ρ̅_ m0/3 = H^2_0 Ω_ m, Ω_ m is the cosmological matter density parameter, H_0 is the present-day value of the Hubble parameter, and n, c_1 and c_2 are dimensionless parameters of the model. The scalaron field, Eqn. (<ref>), takes the form:

    f_R = -c_1/c^2_2n(-R/m^2)^n-1/[(-R/m^2)^n + 1]^2 .


To match the background expansion to that in the ΛCDM model, we set

    c_1/c_2 = 6Ω_Λ/Ω_ m ,

where Ω_Λ≡ 1 - Ω_ m. 

The expression of the scalaron field, Eqn. (<ref>), simplifies when the background value of the Ricci scalar satisfies |R̅| ≫ m^2. From

    -R̅≈ 8π Gρ̅_ m - 2f(R̅) ≈ 3 m [ a^-3  + 2/3c_1/c_2] ,

when Ω_Λ≈ 0.7 and Ω_ m≈ 0.3 we find |R̅| ≈ 40 m^2 ≫ m^2.
Thus

    f_R ≈ -nc_1/c^2_2(m^2/-R)^n+1 .


The model then has two remaining free parameters, n and c_1/c^2_2. The latter is related to the present-day value of the background scalaron, f_R0,

    c_1/c^2_2 = -1/n[ 3(1 + 4Ω_Λ/Ω_ m) ]^n+1 f_R0 .

Hence, the choice of f_R0 and n fully specifies the model. Here we focus on the case of n=1, which is the most well-studied case in the literature.



§ SIMULATIONS AND HALO/GALAXY CATALOGUES





Here we present a description of the simulations used, the construction of halo catalogues, and the HOD prescription used to populate dark matter haloes with galaxies. 




 §.§ Numerical simulations


As we are interested in the effects of f(R) gravity on large scales, we choose three Hu-Sawicki models with n=1 and |f_R0| = 10^-6,10^-5,10^-4 (which we hereafter refer to as F6, F5 and F4, respectively) and the ΛCDM model which assumes GR. Despite the observational tensions faced by f(R) models with |f_R0| > 10^-5 (see e.g., ) it is interesting to consider a wide range of f(R) models to study their impact on the halo/galaxy clustering.

We use the elephant (Extended LEnsing PHysics using ANalaytic ray Tracing) simulations executed using the code ecosmog <cit.>, which is based on the adaptive mesh refinement (amr) N-body code ramses <cit.>. Table <ref> lists the properties of the simulations used in our analysis. The cosmological parameters were adopted from the best-fitting values to the WMAP 9 year CMB measurements <cit.>. All simulations use N_ p = 1024^3 particles with a mass of m_ p=7.798× 10^10 h^-1 M_⊙ to follow the evolution of the dark matter distribution in a volume of V_box = (1024 h^-1Mpc)^3. The initial conditions were generated at z_ini=49 using the MPgrafic code <cit.>. All simulations were run using the same initial conditions up to the present time, z_fi = 0, generating 37+1 snapshots. Here, we analyse the outputs at z=0.5. 




 §.§ Halo catalogues and mass function



Dark matter haloes are the building blocks of large-scale structure and the hosts of galaxies. Therefore, the study of their statistical properties, such as their abundance and clustering, is of great importance in understanding the nature of gravity. 
The halo catalogues were produced using the rockstar halo finder code <cit.>. rockstar calculates halo masses using the spherical overdensity (SO) approach <cit.>, including all  particles and substructures in the halo. We keep only the `parent' halo, omitting other substructures from our analysis.

We define the mass of a halo as M_200c, the mass within a sphere of radius r_200c,  which is the radius within which the mean overdensity is 200 times the critical density of the universe ρ_c,

    M_200c = 4π/3 200ρ_c r^3_200c .
 
 
The dark matter halo mass function (HMF) quantifies the number density of dark matter haloes as a function of their mass. The HMF is sensitive to the cosmological parameters, Ω_ m, Ω_Λ, and σ_8, and to modifications to gravity.
The ΛCDM model predicts an HMF in which the number of haloes increases with decreasing halo mass. f(R) models predict more haloes than the ΛCDM model at almost all masses due to the enhancement of gravity. 
Theoretically, the halo mass function is given by <cit.>

    d n/ d M_200c = f(σ) ρ̅_ m/M^2_200c|  dlnσ^-1/ dln M_200c| ,

where σ is the linear theory variance in the matter perturbation, ρ̅_ m is the mean density of the Universe and f(σ) is an analytical fitting formula. The cumulative number density of haloes above the mass M_200c is:

    n(> M_200c) = ∫^∞_M_200c d n/ dlog_10 M_200c dlog_10 M_200c .

We compare the fitting formula of <cit.> (hereafter Tinker10) to the simulation results. Tinker10 calibrated their fitting formula using a SO algorithm to identify dark matter haloes in numerical simulations which is consistent with the approach used in rockstar. The analytical predictions were computed by using the online tool HMFcalc[<http://hmf.icrar.org/>] <cit.>.

Fig. <ref> shows the cumulative halo mass function (cHMF) measured from the simulations and the relative difference between the f(R) models and GR at z=0.5.  As expected, the largest deviation from GR is displayed by the F4 model (red line; ).  The lower panel of Fig. <ref> shows that the cumulative halo mass function in the F4 model reaches a difference with respect to GR of >50 percent for haloes of mass M_200c > 10^14.3 h^-1 M_⊙. The maximum difference found between F5 (green line) and GR reaches 25 percent for haloes with mass M_200c≈ 10^13.2 h^-1 M_⊙. On the other hand, for the F6 model (blue line) we see that for very massive haloes, the halo mass function is the same as that in GR. This is because the chameleon mechanism works efficiently for such haloes to suppress the effects of the enhancement to gravity. These differences are purely the result of the modified gravitational force in f(R) models. The stronger deviation of F4 from GR is due to the inefficient screening mechanism in this model compared to the screening in the F6 model.

To make a direct comparison between the halo and galaxy clustering we select a halo population from the simulations by fixing the halo number density (in this case we take the number density of the BOSS-CMASS-DR9 sample at z=0.5, n_ h = n_ g = 3.2 × 10^-4 h^3 Mpc^-3; ) and selecting haloes above the mass threshold corresponding to that number density. The horizontal dashed line in Fig. <ref> corresponds to the halo number density used to define our halo sample. The minimum mass that defines the halo sample for each model is: 7.643 × 10^12 h^-1M_⊙ (GR), 7.798 × 10^12 h^-1 M_⊙ (F6), 9.124 × 10^12 h^-1 M_⊙ (F5) and  8.734 × 10^12 h^-1 M_⊙ (F4). The fact that F5 has a higher minimum mass is because this model produces more haloes with mass M_200c∼ 10^13 h^-1 M_⊙ (as we can see from the lower panel of Fig. <ref>) than F6 and F4. For F4, many of the medium-mass haloes have merged to form more massive haloes, hence this model predicts fewer smaller haloes than F5.






 §.§ HOD prescription and galaxy catalogues



 

To compare the simulations with observations one has to populate dark matter haloes with galaxies. This can be done using one of a number of empirical techniques depending on the physical application we are interested in, such as subhalo abundance matching <cit.>, the conditional luminosity function <cit.> or the halo occupation distribution <cit.>. These empirical descriptions of the galaxy-halo connection have the flexibility to give accurate reproductions of observational estimates of galaxy clustering.   
A second, more expensive but physically motivated method is hydrodynamical simulation <cit.>. A third possibility, which retains the physical basis of hydrodynamical simulation at a fraction of the computational cost is semi-analytical modelling of galaxy formation <cit.> in which an N-body dark matter-only simulation is populated with galaxies after solving a set of coupled differential equations. To date, little work has been done to study galaxy formation and clustering in modified gravity models (see   for an example), so here we will resort to the empirical approach of HOD modelling. 

We populate haloes using a functional form for the halo occupation distribution (HOD; ) with five parameters, as used by <cit.>. 

In this form, the mean number of galaxies in a halo of mass M_ h (in our case M_ h = M_200c) is the sum of the mean number of central galaxies plus the mean number of satellite galaxies,

    ⟨ N(M_ h) ⟩   =   ⟨ N_ c(M_ h) ⟩ + ⟨ N_ s(M_ h) ⟩ ,
    ⟨ N_ c(M_ h) ⟩   =   1/2[ 1 + erf( log_10 M_ h - log_10 M_min/σ_log M) ] ,
    ⟨ N_ s(M_ h) ⟩   =   ⟨ N_ c(M_ h) ⟩( M_ h - M_0/M_1)^α ,

and ⟨ N_ s(M_ h) ⟩ = 0 if M_ h < M_0. ⟨ N_ c/s(M_ h) ⟩ is the average number of central or satellite galaxies, respectively, in a halo of mass M_ h. The model depends on five parameters: M_ min, M_0, M_1, σ_log M and α. From Eqns. (<ref>) and (<ref>) we can see that M_min and M_0 represent the halo mass threshold to host one central or one satellite galaxy, respectively. Also, we assume that central galaxies are placed at the centre of their host haloes and satellite galaxies are orbiting inside haloes with M_ h≥ M_0. The satellite galaxies are radially distributed, between r=[0,r_200c], following the Navarro-Frenk-White (NFW) profiles of their host halo <cit.>.

We generate five galaxy catalogues (one for each independent realisation of the density field) for every gravity model following the prescription described above. The galaxy catalogues match the galaxy number density of the BOSS-CMASS-DR9 sample at z=0.5 (n_ g = 3.2 × 10^-4 h^3  Mpc^-3 ; ) and the galaxy two-point correlation function across all gravity models (more details are presented below). The BOSS-CMASS sample is dominated by LRGs which are massive galaxies typically residing in haloes with M_h ∼ 10^13 h^-1M_⊙ <cit.>. Hence, given the mass resolution of the elephant simulations, these runs are suitable to study the impact of f(R) gravity on galaxy clustering. 

The values of the HOD parameters used to populate the GR simulations are those inferred from the abundance and clustering measured for the BOSS-CMASS-DR9 galaxy sample <cit.>:

    log_10(M_ min/[h^-1M_⊙])    =    13.09 ,
    log_10(M_1/[h^-1M_⊙])    =    14.00 ,
    log_10(M_0/[h^-1M_⊙])    =    13.077 ,
    σ_log M   =    0.596 ,
    α   =    1.0127.


To find the f(R) HOD parameters, we use the simplex algorithm of <cit.> to search through the 5D parameter space. We start the algorithm with an initial guess at the values of  the HOD parameters, then the code walks through the 5-dimensional HOD parameter space looking for the values that minimize the root-mean-square difference (rms) of the two-point correlation function between f(R) and GR models. We measure the correlation function using 40 logarithmically spaced radial bins between 1 - 80 h^-1 Mpc. The fractional difference of the galaxy number density is also included (with a weight of 8) in the rms value calculated in order to ensure similar numbers of galaxies in all catalogues. We stop the search when rms < 0.02 (this means that the overall agreement is better than 2%). For the F4 model, the minimum value of the rms we could obtain in practice was ∼ 0.03.



Here, we are interested in the marked correlation function (mCF), which was proposed to highlight the environmental dependence in modified gravity models <cit.>. Hence, the most natural choice is to make the unmarked two-point correlation functions (2PCFs) of the different models as close to each other as possible; otherwise when there is a difference in their mCFs we can not be sure how much of this is due to the different 2PCFs. 

The f(R) HOD parameters were tuned for each model and realisation to match the clustering displayed in the counterpart simulation from the GR suite. The best-fitting values of the HOD parameters for the different realisations are listed in Table <ref>. The variation in the best-fitting parameter values is larger as the modification to gravity increases. We note that, despite the differences between the values of the best-fitting parameters between different models and realisations, the resulting correlation functions and galaxy number densities agree with the GR results to within our target accuracy. 

Note that the HOD parameters are degenerate to some extent, so that a comparison of the values of any single parameter across realisations or models should not be over interpreted. For instance, consider the parameter α that governs the number of satellite galaxies in haloes of a given mass: in the case of the most extreme model, F4, the variation between realisations is ∼ 3%. This difference is small compared with the 1σ scatter of HOD parameter fittings, e.g., <cit.>.   

The left panel of Fig. <ref> shows the HOD for CMASS galaxies at z=0.5. The gradual transition from zero to one galaxy per halo is determined by the values adopted for log M_ min and σ_log M for central galaxies (dashed lines). The appearance of satellites in haloes (dotted lines) is dictated by the values of M_1 and M_0, and the rapid increase in the satellite content of haloes with increasing halo mass is governed by α. The HOD parameters are adjusted in the f(R) models to approximately reproduce the abundance and clustering of CMASS galaxies realised in GR. We note that the resulting HODs are very similar between f(R) gravity and GR.

The right panel of Fig. <ref> shows the distribution of the number of galaxies as a function of the host halo mass (M_200c). We see that most galaxies are found in haloes with mass 10^13 < M_200c/[h^-1 M_⊙] < 10^14. We also note that the F5 and F4 models produce more galaxies than GR and F6 in this mass range. This is because the abundance of haloes in this mass range is boosted in F5 and F4, as we can see from the relative differences of the cHMFs presented in the lower panel of Fig. <ref>. Analysing the distribution of galaxies, we find good agreement between the five realisations. 

From the distribution of galaxies as a function of host halo mass plotted in Fig. <ref> we note that ≈ 0.1% galaxies reside in poorly resolved haloes (M_200c < 10^12h^-1 M_⊙). The inclusion of these galaxies in the final catalogues does not affect the clustering results.





§ RESULTS


Upcoming galaxy surveys will allow us to measure the clustering of galaxies to an unprecedented level of accuracy with the aim of developing a better understanding of the nature of dark matter, dark energy and the evolution of galaxies through cosmic time. 
In this section we present the statistical tools which can be used to characterise the halo and galaxy distributions in different gravity models. 
This is the first time that the halo and galaxy clustering has been studied to this level of detail in f(R) gravity models. 




 §.§ 2-point correlation function


To characterise the clustering of dark matter haloes and galaxies, we use the two-point correlation function, ξ(r). This is defined as the excess probability, compared with that expected for a random distribution, of finding two haloes (or galaxies) contained in volume elements dV_1 and dV_2 at a separation r 
<cit.>:


    dP_12(r) ≡n̅^2 [1 + ξ(r)]dV_1 dV_2 ,

where n̅ is the mean halo (galaxy) number density. The 2-point correlation functions and therefore the marked correlation functions are measured within the range 1-80 h^-1Mpc (for details see Sec. <ref>).

First we study the clustering of dark matter haloes, ξ_ h, (left panel of Fig. <ref>), for our halo samples with n_ h = 3.2× 10^-4 h^3Mpc^-3. 
Although this statistic is not directly observable, it is instructive to study the properties of ξ_ h, since this is a first step towards understanding differences in the clustering of galaxies.  

The first thing we notice is that the deviation from GR does not show a monotonic dependence on |f_R0|. More explicitly, F6 and F4 models have weaker clustering than GR, while F5 haloes are more clustered than GR. 

These perhaps counterintuitive results can be explained by considering the following two effects of the enhanced gravity.

Firstly,  stronger gravity means a faster growth of initial density peaks, and therefore more massive structures at late times. This generally leads to a higher mean halo number density above a fixed halo mass threshold. The enhancement of halo formation is not uniform: when screening is efficient, it is stronger in low-density regions than it is in high-density regions; when the screening is less efficient, then the growth of haloes is boosted in all environments, and those in dense regions can be boosted more because they have more matter around them to accrete.

Secondly, enhanced gravity generally leads to a stronger clustering of the structures that are formed from these initial density peaks. However, stronger gravity also means that we can expect more mergers in dense regions, reducing the number of haloes there. The latter effect can be seen by comparing the cHMFs of F5 and F4 in Fig. <ref>.

As we choose the halo mass cut to ensure that we consider the same number of haloes in each model any  differences in ξ_ h come from the different spatial distributions of haloes in the  models. For F6, the deviation from GR is weak and the fifth force is suppressed in high-density regions. As a result small density peaks in low-density regions grow faster than similar density peaks in high-density regions, and more of them make it into the fixed number density halo catalogue than in GR. This makes the haloes less clustered and ξ_ h(r) smaller.

For F5, the enhancement of gravity is stronger and the screening is weaker, so that haloes in all regions experience faster growth; those in high-density environments have a larger supply of raw materials for accretion and growth, so that they are more likely ending up in the final halo catalogue, leading to a stronger clustering and ξ_ h. For F4, the even stronger enhancement of gravity causes more mergers of haloes in dense regions to form even larger haloes, and to maintain the same n̅_ h more haloes in low-density regions have to be included into the halo catalogue, leading to less clustering and smaller ξ_ h.

In the case of the galaxy correlation function (right panel  of Fig. <ref>), as we said before, the HOD catalogues for the f(R) models were created by tuning the parameters (<ref>; see Table <ref>) to approximately match the two-point correlation function in GR (to within 1-3%).






 §.§ Marked correlation function


In this subsection we consider the marked correlations, in which one weights galaxies[For simplicity, we talk about galaxies here, but the same calculation can (and will) be applied to haloes.] by some property or `mark' when estimating clustering statistics. Marked correlations are particularly well-suited to quantifying how the properties of galaxies correlate with environment <cit.>. 
Here, we test the idea proposed by <cit.> that marked correlation functions may show a  clearer signature of modified gravity in the large-scale clustering of galaxies, by up-weighting low density regions, where screening is weak and deviations from GR are strong.   

The marked correlation function is defined as <cit.>:

    M(r) ≡1/n(r)m̅^2∑_ij m_i m_j = 1 + W(r)/1 + ξ(r) ,

where the sum is over all pairs with a given separation, r, n(r) is the number of such pairs and m̅ is the mean mark for the entire sample. In the second equality ξ(r) is the two-point correlation function in which all galaxies (or haloes) are weighted equally. W(r) is derived from a similar sum over galaxy (halo) pairs separated by r, as used to estimate ξ(r), but now each member of the pair is weighted by the ratio of its mark to the mean mark of the full sample. 
The marked correlation function M(r) can be estimated approximately using the simple pair count ratio WW/ DD (where DD is the count of data–data pairs and WW represents the corresponding weighted counts). Hence, no random catalogue is needed for its computation. 
 
The choice of the mark is flexible and depends on the application. Since we are interested in isolating the effects of the chameleon screening mechanism on structure formation, we study the clustering of HOD galaxies using two different definitions of environment: a) the number density field and b) the Newtonian gravitational potential <cit.>.



  §.§.§ Density


For this environment definition we use three marks with an adjustable dependence on density:

    m    =   ( ρ_* + 1/ρ_* + ρ_R)^p ,
    
    m    =   log_10 (ρ_R + ρ_*) ,
    
    m    =   1/√(2π)σ_ρ_Rexp(- (ρ_R - ρ_*)^2/2σ_ρ_R^2) ,

where ρ_R is the galaxy number density in units of the mean galaxy number density, ρ̅, and p, ρ_* and σ_ρ_R are adjustable parameters. 

A crucial step in the estimation of the marked correlation function is the definition of the density. 
We measure the galaxy number density using counts-in-cells (see, for example, ). We divide the simulation box into cells (or cubical boxes) of the same size, and then count the number of galaxies inside each cell. Hence, we can compute the overdensity, δ, as:

    1 + δ = N/N̅≡ρ_R ,

where N is the number of galaxies in each cell and N̅ is the mean number of galaxies in cells of a given size over the simulation volume. To compute the density we have used 60^3 cells of size ∼ 17 h^-1 Mpc. Given the mean galaxy number density of CMASS galaxies, n_ g = 3.2 × 10^-4 h^3 Mpc^-3, we have a mean number of galaxies in the cells of N̅ = 1.59. We checked that reducing the number of cells to 30^3 - 40^3 does not affect our results significantly, while further reducing the number of cells makes the signal weaker; in the limit of 1^3 cell, W(r) becomes identical to ξ(r), as expected. 

The first mark, Eqn. (<ref>), was proposed by <cit.> (hereafter the White–mark), with the motivation being that by up-weighting low density regions (i.e. by choosing p > 0), one might be able to find a signature of modified gravity, since previous studies have shown that the properties of voids are different in modified gravity theories than in GR <cit.>. 
The log mark, Eqn. (<ref>), allows us to up-weight regions with ρ_R > 1, i.e., intermediate and high-density regions. Finally, using the Gaussian-ρ_R mark, Eqn. (<ref>), we are able to control the regions we want to up-weight. Previously, <cit.> found that by using a Gaussian transformation of the density field is it possible to up-weight intermediate density regions and find bigger differences between the clustering of objects in modified gravity and GR models. Keeping this in mind, we use the Gaussian-ρ_R mark to up-weight only intermediate density regions.

It is evident that by using Eqn. (<ref>) one can control the up-weighting by varying the power p and the parameter ρ_*. For simplicity we chose p=1 and ρ_* = 10^-6. With the log-mark, a natural choice of the parameter which controls the enhancement is ρ_* = 1, given m=0 for voids (ρ_R = 0). The parameters we chose for the Gaussian-ρ_R mark are: ρ_* = 1.5 and σ_ρ_R = 0.2, which ensures that we up-weight  intermediate-density regions of interest. The functional form of the marks, Eqs. (<ref>) – (<ref>), is shown in the left-hand panels of Fig. <ref>. We have tried using different values of p, ρ_R and ρ_* but  found that our results do not show significant differences on varying these parameters. We refer to low-, intermediate- and high-density regions as those for which the cells contain N = 1, 2-3 and >4 objects or, equivalently, to cells with ρ_R = 0.62, 1.25-1.88 and >2.51, respectively (see Eqn. (<ref>)).

Fig. <ref> shows the marked correlation functions (mCFs) at z=0.5 measured from the halo (middle panels) and the HOD (left panels) catalogues in the f(R) and GR models. In all cases the marked correlation function goes to unity on large scales as expected (see right-hand expression of Eqn. (<ref>)). The first row of plots in Fig. <ref> shows the mCF using the mark defined by Eqn. (<ref>), the White-mark, with p=1 and ρ_* = 10^-6, the second row shows the log-mark, Eqn. (<ref>) with ρ_*=1, and the third row shows the Gaussian-ρ_R mark with ρ_* = 1.5 and σ_ρ_R = 0.2.  We observe different behaviours: for the White-mark,  Eqn. (<ref>), the marked correlation function is M(r) ≤ 1 at small separations, for the log mark, Eqn. (<ref>), we have M(r) ≥ 1, while for the Gaussian-ρ_R we notice a transition from M(r) ≤ 1 to M(r) > 1 at intermediate scales. 

Analysing the behaviour of the halo marked correlation functions (see middle panel of Fig. <ref>) we find the following features: 


  * The clustering of F6 is almost indistinguishable from that of GR for all three marks, because of the efficient screening.


  * For F5, the stronger growth (see Sec. <ref>) means more clustering of haloes on small scales, which is why W(r) and therefore the marked correlation function is more affected at smaller r.


  * In the case of F4, the higher production rate of massive haloes, driven by the more frequent mergers of lower mass haloes (see Sec. <ref>), leads to the incorporation of haloes into the fixed number density sample which correspond to low density peaks and which are more likely to come from low-density regions. Hence, the probability of finding a pair of tracers (haloes or galaxies) increases at intermediate separation r due to presence of these low mass haloes.  


The right columns of Fig. <ref> show that galaxies qualitatively mimic the marked clustering of haloes (at least for the White and log marks). Hence, the behaviour of the galaxy marked correlation functions can be understood following the same explanation as presented above for haloes. It is interesting to notice that even with the added complexity of populating haloes with HOD galaxies, the qualitative behaviour of the marked correlation functions preserves, suggesting that a true physical feature is being observed here.

For the Gaussian-ρ_R mark, Eqn. (<ref>), 
which enhances intermediate-density regions (cells with 2 or 3 haloes/galaxies), we found that the F4 galaxy marked correlation function reaches a maximum of 20% for the lowest separation bin, while F6 predicts a difference of 5% and F5 keeps closer to GR with a difference of ∼ 3%.






  §.§.§ Gravitational potential



Our second definition of environment is based on the Newtonian gravitational potential produced by dark matter haloes. The dark matter haloes in our simulations are reasonably well described by a NFW density profile <cit.>:

    ρ_NFW = ρ_ s/(r/r_ s)(1+r/r_ s)^2  ,

where r_ s is the scale radius where the profile has a slope of -2 and ρ_ s is the characteristic density.
The Newtonian gravitational potential is obtained by solving the Poisson equation, ∇^2 Φ_ N = 4π G ρ_NFW, for the NFW density profile Eqn. (<ref>) 
<cit.>:

    Φ_ N = - G M_200c/r_200cln(1 + c)/ln(1 + c) - c/(1 +c) ,

where G is Newton's gravitational constant, M_200c was defined in Eqn. (<ref>) and c is the concentration parameter defined as c≡ r_200c/r_ s. Previous studies have used the Newtonian gravitational potential in modified gravity to characterise local variations in the strength of gravity (see e.g., ).

For this environment definition we define a Gaussian mark which allows us to up-weight galaxies in some regions of interest,

    m = 1/√(2π)σ_Φexp[- (log_10(|Φ_ N|) - Φ_*)^2/2σ_Φ^2] ,

where Φ_* and σ_Φ are free parameters of the mark which control the amplitude and width of the regions highlighted. As we can see from the distribution of galaxies as a function of host halo mass (right panel of Fig. <ref>), most galaxies live in haloes with masses between 10^13 < M_200c/[h^-1 M_⊙] < 10^14 (which correspond to the mass range of groups of galaxies). Hence, we use the Gaussian-Φ_ N mark to up-weight galaxies contained in these haloes. In principle we should be able to find a bigger difference in the clustering between the GR and f(R) models using this mark, as suggested by the cumulative halo mass function (lower panel of Fig. <ref>).

The value of the centre of the Gaussian is Φ_* = -5.295. This value was found by computing the Newtonian gravitational potential for each galaxy in Box 1 for GR, then we pick the maximum value found for haloes with M_200c = 10^14 h^-1 M_⊙ (Φ_ max) and the minimum value for haloes with M_200c = 10^13 h^-1 M_⊙ (Φ_ min), finally we take Φ_* = (Φ_ max + Φ_ min)/2. We tried different values of the width, finding that σ_Φ = 0.1 best ensures that we only up-weight galaxies in the haloes of interest. The functional form of the Gaussian-Φ_ N mark, Eqn. (<ref>), is shown in the left-hand panel of Fig. <ref>. 

The halo and galaxy marked correlation function is presented in the middle and right panels of Fig. <ref>, respectively. The results can be summarised as follows:


  * In the case of the halo/galaxy marked correlation function (middle and right panels of Fig. <ref>, respectively), the two-point correlation function (used as the denominator of Eqn. (<ref>)) is lower than the weighted correlation function, leading to M(r) ≥ 1 for all gravity models, due to the stronger clustering of the up-weighted haloes in the mass range M_200c/[h^-1 M_⊙] = [10^13,10^14].

  * F6 predicts almost an identical halo/galaxy marked clustering to that in GR, which is consistent to our understanding that the screening mechanism in this model works efficiently in haloes of the mass range up-weighted.

  * For F4 haloes, the mCF is higher than the 2PCF for the reason given in the first bullet point above. However, in this model a larger fraction of haloes in the mass range M_200c/[h^-1 M_⊙] ∈ [10^13,10^14] are formed from low initial density peaks (due to stronger gravity) which are not very strongly clustered, such that the up-weighting of them – while making M(r)≥1 – does not lead to a M(r) as large as in GR. This leads to Δ M(r)/M_ GR(r)<0 for F4. For F5 haloes, the fifth force is strong enough to enhance their clustering, but not too strong to produce excessive merging, and so the up-weighting using the Gaussian mark increases the mCF as significantly as in GR.

  * For galaxies, a key difference from haloes is that a halo can host several galaxies while some haloes do not host galaxies at all. In F4 and F5, more relatively low initial density peaks have been promoted to the halo mass range M_200c/[h^-1 M_⊙] ∈ [10^13,10^14] due to the enhanced gravity, and at the same time some high density peaks have grown out of this mass range. This means that if we up-weight galaxies whose host haloes are in this mass range, we end up with more central and fewer satellite galaxies, and more of them are hosted by haloes from lower initial density peaks. By the same reasoning as above, while we still have M(r)>1 for these models, it is smaller than in GR and F6. In particular, we have noticed that Δ M(r)/M_ GR(r) reaches 5∼10% for F5 and 20∼30% for F4 in r=2∼5h^-1Mpc. These results are very stable, and change very little across the different simulation realisations.


Also we note that the differences between the f(R) and GR models are boosted when we use additional information to the density field. This can be seen by comparing the right panels of Fig. <ref> with the right panel of Fig. <ref>. The differences get larger in such cases because the galaxy density field and galaxy distribution have been tuned to match between the different models (see Sec. <ref>). In all cases we observe that signals above 20 h^-1 Mpc become identical between models. This is because the marked correlation function is the ratio of two correlation functions (see right hand expression of Eqn. (<ref>)) and we have ξ(r) ∼ W(r) for r>20 h^-1 Mpc.

From the observational point of view, we can measure the Newtonian gravitational potential from the X-ray temperature of galaxy clusters (see e.g., ), the gas mass fractions of clusters and the escape velocity profile, v_esc(r) <cit.>. Hence, if we reconstruct the gravitational potential from the observations mentioned above and use a mark that is a  function of the potential, similar to Eqn. (<ref>), then we can test this approach and potentially find a measurable signature of modified gravity. One caveat is that the gravitational potential constructed in this way is the dynamical potential, while in this study we have used lensing potential of haloes <cit.>.




§ DISCUSSION AND CONCLUSIONS



We study the clustering of haloes and galaxies in four different cosmologies: a ΛCDM 
model which is based on general relativity and three  f(R) chameleon models with fixed n=1 and |f_R0| = 10^-6, 10^-5,10^-4 (denoted F6, F5 and F4). We analyse the output of dark matter-only N-body simulations related to these models at z=0.5. 

First, we study the cumulative halo mass function, finding that the F4 model predicts more haloes than GR at all masses probed by our simulations, with the maximum difference reaching an excess of more than 50 percent for haloes with mass M_200c > 10^14.3 h^-1 M_⊙. These differences occur due to the enhancement of gravity in f(R) gravity, which results in the production of more massive haloes in F4 than GR through faster accretion and more frequent merging of small haloes. The differences found in F5 reach 25% for haloes with masses 10^13 h^-1 M_⊙ < M_200c < 10^14 h^-1 M_⊙ where the screening mechanism at this mass scale is inefficient for this model. F6  shows the smallest difference from GR because in this model the chameleon screening is strong in haloes with mass M_200c > 10^13 h^-1 M_⊙, thereby suppressing the effects  of the fifth force. 

We populate dark matter haloes with galaxies using a halo occupation distribution, using a five-parameter model which treats separately central and satellite galaxies, with the values of the parameters as used in  <cit.> to reproduce the clustering of CMASS galaxies with a density number, n_ g = 3.2 × 10^-4 h^3  Mpc^-3 <cit.> for our GR simulations. We tuned the parameters to match the galaxy number density and two-point correlation function of GR to within 1-3% for the f(R) models. 
The galaxy two-point correlation functions for the f(R) and GR models are presented in the right plot of Fig. <ref>.

Then we study the two-point clustering of dark matter haloes. We choose samples of haloes with fixed halo number density, n_ h = n_ g, resulting to different mass cutoffs in our halo catalogues for all gravity models: 7.643 × 10^12 h^-1M_⊙ (GR), 7.798 × 10^12 h^-1M_⊙ (F6), 9.124 × 10^12 h^-1M_⊙ (F5) and  8.734 × 10^12 h^-1M_⊙ (F4). We find significant differences in the clustering of dark matter haloes for f(R) models with respect to the GR predictions. The maximum difference between F4 and GR is ∼20%, while for F5 and F6 it is less than 5%. Also we note that haloes in F5 are more clustered than those haloes in the ΛCDM model, whereas for F6 and F4 haloes are less clustered than their GR counterparts. These results are the effects of the enhancement of gravity which means a stronger growth of density peaks and therefore more massive structures at late times which gives a stronger clustering of the structures that are formed from these density peaks.  

To investigate whether or not these differences could be boosted by using an alternative approach to measure galaxy clustering, we used the marked correlation function <cit.>. For this purpose we use two definitions for the environment of galaxies/haloes: a) the number density field and b) the Newtonian gravitational potential of the host halo. For the former we analyse three marks: i) an inverse power-law which enhances low-density regions (see Eqn. (<ref>)), ii) a log-transform mark which up-weighs intermediate- and high-density regions (see Eqn. (<ref>)) and iii) a Gaussian mark given by Eqn. (<ref>) which allows us to up-weight only intermediate-density regions, ρ_R = 1.25 - 1.88. For the latter we use a Gaussian mark which allows us to up-weight haloes (and galaxies within those haloes) with mass 10^13 < M_200c/[h^-1 M_⊙] < 10^14. 

We found that the halo and galaxy marked correlation functions for F6 is indistinguishable from GR using all marks, except for the galaxy Gaussian (ρ_R and Φ_ N) marked correlation functions which predict differences at most ∼5% from GR. 
For the F5 (F4) model, we notice that galaxies mimic the marked clustering at least for the White and log density-marks finding differences of 5% (2.5%) and 2.5% (2.5%), respectively. 
On the other hand, we observe that with the Gaussian marks (density field and gravitational potential) the difference in the galaxy marked correlation function is boosted, especially for F4, producing a difference of 20% (using density field) and 30% (using gravitational potential) with respect GR.

The galaxy marked correlation functions show smaller differences between the f(R) and GR models for the density marks, Eqns. (<ref>)–(<ref>), than in the case when using the gravitational potential mark, Eqn. (<ref>), this is because the galaxy density field has been tuned to match between the different models. One caveat for our results is that there will be systematics when estimating the mark for observational samples. 

Another important feature we observe from marked correlation functions is that the signal above 20   h^-1 Mpc does not distinguish between models (see corresponding plots of Figs. <ref> and <ref>). Instead, the measurable differences are on small scales. To improve our predictions on sub-Mpc scales we need to perform higher resolution  simulations, but we leave this for future work.

<cit.> recently found that using the dark matter distribution and the White-mark, Eqn. (<ref>) with ρ_* = 4 and p=10, the difference between F4 and GR marked correlation functions can reach a maximum of 37% at r=1.81 h^-1 Mpc. These results can not readily be compared with ours, since we consider dark matter haloes and galaxies rather than the dark matter itself. Furthermore, we employ a different definition of density (counts-in-cells versus the cloud-in-cell smoothing used by Valogiannis & Bean). Although their simulations are similar resolution to the ones we use, the volume of our boxes is ∼ 60 times larger, which allows more robust clustering measurements. 

If we consider the statistical errors presented in this paper, then the differences between the f(R) signals with respect GR are significant. It should be feasible to test the differentiating power of marked clustering statistics using real data from current galaxy surveys. In future work we need to improve the resolution of the simulations and make more realistic mock galaxy catalogues to allow a  fairer comparison with upcoming  observations.

Here we have demonstrated the potential of the marked correlation function to differentiate between gravity models. The next step is to extend these calculations, which were presented for massive galaxies, to the emission line galaxies that will be selected by the DESI and Euclid redshift surveys. 



§ ACKNOWLEDGEMENTS

The authors wish to thank Lee Stothert for providing the code to compute the 2-point correlation functions and the weighted number counts.
Also we would like to thank the referee for their positive and helpful comments about our paper.
CH-A is supported by the Mexican National Council of Science and Technology (CONACyT) through grant No. 286513/438352. BL is supported by the European Research Council (ERC-StG-716532-PUNCA). This project has received funding from the European Union's Horizon 2020 Research and Innovation Programme under the Marie Skłodowska-Curie grant agreement No 734374, LACEGAL. We acknowledge support from STFC Consolidated Grants ST/P000541/1, ST/L00075X/1. 
This work used the DiRAC Data Centric system at Durham University, operated by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (<www.dirac.ac.uk>). This equipment was funded by BIS National E-infrastructure capital grant ST/K00042X/1, STFC capital grants ST/H008519/1 and ST/K00087X/1, STFC DiRAC Operations grant ST/K003267/1 and Durham University. DiRAC is part of the National E-Infrastructure.







mnras
 























