The Bispectrum of f(R) Cosmologies
Héctor Gil-Marín, Fabian Schmidt, Wayne Hu, Raul Jimenez, Licia Verde
PPrepared for submission to JCAP
The Bispectrum of f ( R ) Cosmologies
H´ector Gil-Mar´ın, a,b
Fabian Schmidt, c Wayne Hu, d Raul Jimenez, e,b and Licia Verde e,b a Institute of Space Sciences (IEEC-CSIC), Faculty of Science, Campus UAB, Bellaterra 08193, Spain b Institute of Sciences of the Cosmos (ICC-IEEC), University of Barcelona, Barcelona 08024, Spain c Theoretical Astrophysics, California Institute of Technology, Mail Code 350-17, Pasadena, California91125 d Kavli Institute for Cosmological Physics, Department of Astronomy & Astrophysics, University ofChicago, Chicago, IL 60637 e ICREA Instituci´o Catalana de Recerca i Estudis Avan¸catsE-mail: [email protected], [email protected], [email protected],[email protected], [email protected]
Abstract.
In this paper we analyze a suite of cosmological simulations of modified gravitationalaction f ( R ) models, where cosmic acceleration is induced by a scalar field that acts as a fifth force onall forms of matter. In particular, we focus on the bispectrum of the dark matter density field on mildlynon-linear scales. For models with the same initial power spectrum, the dark matter bispectrum showssignificant differences for cases where the final dark matter power spectrum also differs. Given thedifferent dependence on bias of the galaxy power spectrum and bispectrum, bispectrum measurementscan close the loophole of galaxy bias hiding differences in the power spectrum. Alternatively, changesin the initial power spectrum can also hide differences. By constructing ΛCDM models with verysimilar final non-linear power spectra, we show that the differences in the bispectrum are reduced( (cid:46) f ( R ) model. This weak dependence of the matter bispectrum ongravity makes it useful for breaking degeneracies associated with galaxy bias, even for models beyondgeneral relativity.— a r X i v : . [ a s t r o - ph . C O ] N ov ontents f ( R ) Gravity 23 Simulations 34 Power Spectrum and Bispectrum 35 Results 4
Observations of Type Ia supernovae suggest that the Universe has been accelerating since redshift z ∼ . f ( R ) class of models (see [19]and references therein). These models produce accelerated expansion through a modification of theEinstein-Hilbert action by an arbitrary function of the Ricci scalar R . As a consequence, an extrapropagating scalar field appears that mediates a fifth force on all forms of matter. The range of thisforce depends on the functional form of f ( R ). In order to satisfy solar system tests, f ( R ) models areoften chosen to present a chameleon behavior. The chameleon mechanism makes the extra scalar fieldbecome increasingly massive in higher-curvature regions, suppressing the range of the fifth force indense environments.In previous works, cosmological simulations [9] have been used to study the power spectrum[10] and halo statistics [13] of these kinds of models. More recent studies with higher resolutionhave confirmed these previous results [21] and extended the investigation to smaller scales. In thepresent work we focus on how the dark matter bispectrum is modified in this class of models. Whilethese models also predict a non-linear matter power spectrum different from the ΛCDM one, it isnevertheless interesting to look at the bispectrum for at least two reasons: a) except for gravitationallensing, measurements of clustering yield the galaxy or the baryon power spectrum, not the darkmatter one: as baryonic physics and galaxy formation are complicated phenomena, the observedpower spectrum may be biased, i.e. may differ significantly from the dark matter one; the bispectrumis well known for helping disentangle effects of gravity from effect of biasing e.g., [4, 20]. b) once weallow ourselves to consider non-standard models, the initial (linear) matter power spectrum does nothave to be the power-law ΛCDM one to reproduce the observations. The form of the bispectrumkernel is a possible “signature” of gravity as it gets modified by any modifications from GR behaviore.g., [15]. – 1 –ere we pay special attention to see whether the bispectrum can be used to break degeneraciesbetween models with the same observed power spectrum and the same cosmology, but different gravity.We begin in § f ( R ) models, in § § § § f ( R ) Gravity
The f ( R ) class of models generalizes the Einstein-Hilbert action to include a function f ( R ) of theRicci scalar R , S = (cid:90) d x √− g (cid:20) R + f ( R )16 πG + L m (cid:21) . (2.1)Here L m is the Lagrangian of matter and we have assumed c = (cid:126) = 1. For standard GR witha cosmological constant, f ( R ) = − πGρ Λ , whereas for modified gravity, the force modification isassociated with an additional scalar degree of freedom f R ≡ df /dR . In particular, in this paper weuse the model for f ( R ) proposed by [6], f ( R ) ∝ RAR + 1 , (2.2)where A is a constant with dimensions of length squared. We can write this equation as a functionof its derivative evaluated at ¯ R (the background curvature today), namely f R . We adjust theproportionality constant to match some effective cosmological constant ρ Λ in the limit where f R → AR (cid:29) f ( R ) can then be approximated as, f ( R ) = − πGρ Λ − f R ¯ R R . (2.3)The modified Einstein equations can be computed by varying the Einstein-Hilbert action (Eq. 2.1)with respect to the metric. We work in the quasistatic limit where the time derivatives are negligiblecompared to the spatial derivatives. In this regime, valid on scales much smaller than the horizon1 /H , the trace of the modified Einstein equations yields the f R field equation, ∇ δf R = a δR ( f R ) − πGδρ m ] , (2.4)where a is the scale factor, δf R = f R ( R ) − f R ( ¯ R ), δR = R − ¯ R and δρ m = ρ m − ¯ ρ m . Here ¯ R is thebackground curvature that can be approximated by a ΛCDM universe for | f R | (cid:28) ρ m (¯ ρ m ) isthe (background) matter density.On the other hand, the time-time component of the Einstein equations yields the modified Poissonequation, ∇ Ψ = 16 πG a δρ m − a δR ( f R ) (2.5)where Ψ = δg / (2 g ) is the Newtonian potential.For small fluctuations of the field, we can approximate δR (cid:39) ( dR/df R ) | ¯ R δf R . We will refer tothis linearization as the non-chameleon limit. Conversely if the field fluctuations are large enoughsuch that δR ( f R ) cannot be linearized, the chameleon mechanism operates. We will refer to use ofthe exact, as opposed to linearized, equations as full f ( R ) or just chameleon models.The linearized field equations formed by Eqs. 2.4 - 2.5 can be solved for the Newtonian potentialas a function of the density field. In the linear approximation for δR , these two equations in Fourierspace yield, k Ψ( k ) = − πG (cid:18) −
13 ¯ µ a k + ¯ µ a (cid:19) a δρ m ( k ) . (2.6)– 2 –his equation is identical to the one in GR but with a modification of the gravitational constant, G eff ( k , t ) ≡ G (cid:18) −
13 ¯ µ ( t ) a ( t ) k + ¯ µ ( t ) a ( t ) (cid:19) . (2.7)Here µ ( R ) ≡ (3 df R /dR ) − / is the effective mass of the scalar field f R and ¯ µ just stands for µ ( ¯ R ). Thedependence on time is introduced though ¯ R ( t ). Note that when f R → G eff → G and we recover theΛCDM limit, as expected. It is interesting to see that for a given value of f R there are two differentregimes for G eff , depending on whether the physical scale we are studying is larger or smaller thanthe inverse mass of the field. On large scales k (cid:28) µ ( t ) a ( t ), G eff → G and gravity behaves as GR,whereas on small scales k (cid:29) µ ( t ) a ( t ) G eff → G/ /
3. In other words, in Eq. 2.6, one assumes that the mass of the scalar field µ only depends ontime and is the same in all regions of the Universe at a given epoch. However, for cosmologicallyinteresting values of µ the field is then essentially massless within the Solar System. The presenceof such a scalar field (fifth force) is ruled out by light deflection and time delay measurements inthe Solar System, which are all consistent with GR. In the full non-linear f ( R ) theory, R ∝ f − / R can become very large in dense environments, suppressing the field and restoring the GR relation δR = 8 πGδρ m (Eq. 2.4). Thus, gravity is not modified in the same way everywhere, but depends onenvironment. In regions with large potential wells (inside halos) the mass of the scalar field becomeslarge and therefore the effective range of interaction of this field shrinks recovering GR. We call thisthe chameleon mechanism. The simulations used in this paper are described in previous works [9, 10, 13]. Briefly, the fieldequation for f R (Eq. 2.4) is solved on a regular grid using relaxation techniques and multigrid iter-ation. The potential Ψ is computed from the density and f R fields following Eq. 2.5 using the fastFourier transform method. The dark matter particles are then moved according to the gradient ofthe computed potential, −∇ Ψ, using a second order accurate leap-frog integrator.The simulations were run using the values of | f R | = 10 − , 10 − , 10 − (both for chameleonand non-chameleon cases) and 0, which is equivalent to ΛCDM . The background expansion historyfor all cases differ from ΛCDM only at O ( f R ) and are hence practically indistinguishable. Thecosmology used is Ω Λ = 0 .
76, Ω m = 0 .
24, Ω b = 0 . H = 73 km/s/Mpc and initial power incurvature fluctuations A s = (4 . × − ) at k = 0 .
05 Mpc − with a tilt of n s = 0 . a = 0 .
02 and are integrated in time in steps of ∆ a = 0 . L = 256 and 400 Mpc/ h with 512 grid cells and 256 particles. For each box size, we have 6 runs for each value of f R , with differentrealizations of the initial conditions. The simplest statistic of interest of the matter density field is the power spectrum P ( k ), defined bythe second moment of the Fourier amplitude of the density contrast, (cid:104) δ ( k ) δ ( k (cid:48) ) (cid:105) ≡ (2 π ) δ D ( k + k (cid:48) ) P ( k ) , (4.1) In this paper we are always using negative values for f R , so when we talk about f R we refer to its absolute value. – 3 –here (cid:104) . . . (cid:105) denotes the ensemble average over different realizations of the Universe. By statisticalisotropy, the power spectrum does not depend on the direction of the k -vector. In practice we onlyhave one observable Universe, so the average (cid:104) . . . (cid:105) cannot be computed. However, using the isotropyof the power spectrum we can compute the average over all different directions for each k -vector.Note also that P ( k ) is defined to be real. Since k = − k (cid:48) , δ ( k ) δ ( k (cid:48) ) ∼ | δ ( k ) | , which is a real number.The second statistic of interest is the bispectrum B , defined by, (cid:104) δ ( k ) δ ( k ) δ ( k ) (cid:105) ≡ (2 π ) δ D ( k + k + k ) B ( k , k , k ) . (4.2)The Dirac delta function δ D , ensures that the bispectrum is defined only for k -vector configurationsthat form closed triangles: (cid:80) i k i = 0. Note that once the average is taken, the imaginary part of δ ( k ) δ ( k ) δ ( k ) goes to zero.It is convenient to define the reduced bispectrum Q ≡ Q ( k , k , k ) as, Q ≡ B ( k , k , k ) P ( k ) P ( k ) + P ( k ) P ( k ) + P ( k ) P ( k ) , (4.3)which takes away most of the dependence on scale and cosmology. The reduced bispectrum is usefulwhen comparing different models, since it has a weak dependence on cosmology and one can thusbreak degeneracies between cosmological parameters to isolate the effects of gravity. Hereafter, whenwe speak of the bispectrum we are always referring to the reduced bispectrum. In this paper, we present two ways of comparing the f ( R ) and ΛCDM reduced bispectra we obtainfrom N-body simulations. The differences depend on whether the models are matched in their initialor final power spectra. In method A we compare the output bispectra from N-body simulations withthe same initial power spectra. Thus some of the difference in the bispectra can be attributed tothe different amounts of final nonlinear power in the two sets. Method B tries to separate thesecontributions by generating modified initial power spectrum ΛCDM simulations whose power spectraat z = 0 match those of the f ( R ) simulations.For both methods we compute the bispectrum randomly drawing k -vectors from a specified bin,namely ∆ k and randomly orientating the triangle in space. We make the number of random trianglesto depend on the number of fundamental triangle per bin, that scales as k k k ∆ k [14]. In thispaper we always choose ∆ k = 3 k min . For the equilateral case, at scales of k ∼ . h /Mpc we aregenerating ∼ × triangles. We have verified that increasing the number of triangles beyond thisvalue does not have any effect on the measurement. Our first test of f ( R ) vs ΛCDM bispectra utilize the same initial power spectrum. In our f ( R ) models,modifications to gravity go to zero rapidly with redshift and the expansion history differs negligiblyfrom ΛCDM. Thus models with the same initial power spectra as ΛCDM fit observations at highredshift, such as primary CMB anisotropy, equally well.Since all N-body simulations start from the same initial power spectrum, the f ( R ) modificationsto gravity during the acceleration epoch lead to differences in the dark matter power spectra at lowredshift that increase with | f R | as was noted in Fig. 2 of [10]: these differences reach up to ∼ | f R | = 10 − and ∼
10% for | f R | = 10 − for k (cid:39) . h/ Mpc with respect to the ΛCDM model.Bispectra for matched initial power spectra but differing final power spectra is what is usuallycomputed analytically [1, 2]: one predicts (using the modified Euler and continuity equations inperturbation theory) the reduced bispectra for different gravity models starting from a given initial δ k field. One might expect that the reduced bispectra differences are independent of the power spectra,but this is only strictly true for equilateral configuration and only in the tree-level regime (for k < . h /Mpc at z = 0). This is the main caveat of method A: the differences seen in the reduced bispectrum– 4 – . 2 0 . 4 0 . 6 0 . 8- 0 . 10 . 00 . 1 k = 2 k = 0 . 4 h / M p c Qf(R)/QLCDM-1 | f
R 0 | = 1 0 - 4 f u l l f ( R ) n o n - c h a m e l e o n
E q u i l a t e r a l
Qf(R)/QLCDM-1 | f
R 0 | = 1 0 - 4
Qf(R)/QLCDM-1 | f
R 0 | = 1 0 - 5
Qf(R)/QLCDM-1 | f
R 0 | = 1 0 - 5 q / p Qf(R)/QLCDM-1 | f
R 0 | = 1 0 - 6 k [ h / M p c ]
Qf(R)/QLCDM-1 | f
R 0 | = 1 0 - 6
Figure 1 . Relative dark matter reduced bispectrum deviations (following method A) between ΛCDM and f ( R ) models for k = 2 k = 0 . h /Mpc (left panels) and equilateral configuration (right panels) at z = 0as a function of the angle between k and k , namely θ (left panel) and as a function of k (right panel)for | f R | = 10 − , 10 − , 10 − (top to bottom). Blue points (squares) correspond to chameleon simulationsand red points (circles) to non-chameleon. Both ΛCDM and f ( R ) bispectra have been computed from N-body simulations with the same initial conditions. As a consequence, the corresponding final ( z = 0) powerspectra of the compared models are different. Error bars are the 1- σ standard deviation of the ratio of Q values amongst the 6 independent runs. Because of that, the errors due to cosmic variance cancel out. Only L = 400 Mpc /h side-box runs are used. could be due to differences in the final matter power spectrum and not unique signatures of f ( R )gravity.In Fig. 1 we show the dark matter reduced bispectra deviation between ΛCDM and f ( R ) for k = 2 k = 0 . h/ Mpc (left panels) and for equilateral configurations (right panels) for z = 0 accordingto method A. The top panels correspond to f ( R ) theories with | f R | = 10 − ; | f R | = 10 − for middlepanels; and | f R | = 10 − for bottom panels. The blue points correspond to full f ( R ) theories whereasthe red points to non-chameleon ones. Deviations of f ( R ) bispectra with | f R | = 10 − with respectto ΛCDM present a characteristic shape dependence, where the difference is maximal for θ ∼ π , and minimal for θ ∼ . π , for both chameleon and non-chameleon and it increases as the scaleis reduced. A similar trend is present for non-chameleon theories with | f R | = 10 − . On the otherhand, chameleon theories with | f R | = 10 − present a constant deviation from ΛCDM of ∼ | f R | = 10 − , both chameleon and non-chameleon present a constant deviation from ΛCDM of (cid:46)
5% and is consistent with 0. The errors of Fig. 1 are suppressed compared to the individual cosmic– 5 –ariance errors because we are taking the ratio of N-body simulations with the same initial powerspectrum and phases. Therefore, we can conclude that having the same initial conditions, the darkmatter bispectra of ΛCDM and f ( R ) theories is significantly different for | f R | (cid:38) − , especially forelongated triangles ( θ (cid:39) , π ) and can reach deviations in the reduced bispectra up to ∼
10% for | f R | = 10 − and up to ∼
12% for | f R | = 10 − . We see a similar dependence on triangle shape asshown in Fig. 5 of [1] (note that β as defined there is 1 / √ f ( R )).Although differences in the final dark matter power spectra between ΛCDM and f ( R ) theoriesof the same initial power are large and potentially easier to test than those in the bispectra, it ispossible that the galaxy power spectra for ΛCDM and f ( R ) models could still be similar for someparticular galaxy bias model [18]. Since the galaxy bias acts differently on the power spectrum andon the bispectrum, it would be very unlikely that the same galaxy bias could make P ΛCDMgal = P f ( R )gal and Q ΛCDMgal = Q f ( R )gal simultaneously.Conversely, changes in the initial power spectra between the models might conspire to make an f ( R ) model look like a ΛCDM model for the power spectrum at z = 0. These can be hidden fromthe CMB at high redshift if they only occur at high k . Because of that, in the next section we assessthe differences between the ΛCDM and f ( R ) dark matter reduced bispectra in the case where bothmodels have the same final power spectrum. The final power spectra of the f ( R ) models deviate significantly from that of the ΛCDM model withthe same initial conditions (see Fig. 2 of [10]). These deviations reach ∼
50% for | f R | = 10 − ; ∼ | f R | = 10 − ; all at k (cid:39) h /Mpc and at z = 0. That begs the question of whether bispectrumdifferences seen in Method A are driven by these final power spectrum differences or by uniquelygravitational modifications.To address this question, we would like to adjust the initial conditions of the f ( R ) simulationsuntil the final power spectra match that of ΛCDM at z = 0. However, the f ( R ) simulations are com-putationally very expensive (a factor of ∼
20 increase over ordinary GR simulations). Instead we do theconverse: we adjust the initial conditions of the ΛCDM model until its final power spectrum matchesthe f ( R ) simulations. Matching ΛCDM to the f ( R ) simulations still tests whether the remaining bis-pectra difference between the two models reflects gravitational modifications, independently of powerspectrum differences. In order to match final power spectra we need a means of quickly predicting the impact of adjustinginitial conditions in ΛCDM. HaloFit [17] provides an approximate analytic mapping between theinitial and final power spectra. Using HaloFit we can determine the desired initial conditions, run thematching ΛCDM simulations and compare the bispectra with those of the f ( R ) simulations. We first test the accuracy of HaloFit in modeling the ΛCDM simulation results (see Fig. 2).In the HaloFit computation we use the same transfer function as employed in the ENZO code. Wesee that for the L = 400 Mpc/ h runs the data-points with k > k N / (cid:39) . h /Mpc underestimatethe value of the predicted power spectrum by HaloFit. Here, k N is the Nyquist mode defined as k N = πN / p / (2 L ). We can in principle solve this limitation by using smaller boxes. For the L = 256Mpc/ h runs, k N / (cid:39) . h /Mpc and we see that up to this scale the simulation agrees with thetheoretical prediction. However the errors increase considerably as we reduce the box-size. Errorbars in Fig. 2 correspond to 1- σ standard deviation amongst 6 independent runs. Likewise at low k ,the simulations carry large sampling errors even for the largest boxes. To evade these problems, we An alternative (much faster) approach would be, instead of re-running N-body simulations and computing theevolved bispectrum from there, to use fitting formulae from the literature to predict the ΛCDM bispectrum from thematched power spectrum which add a small running of the tilt. We have attempted this route; unfortunately we havefound that available fitting formulae are not sufficiently accurate for our purposes [5]; as we will see we need a relativeaccuracy of better than 4% on the reduced bispectrum which exceeds that of available fitting formulae in the mildlynon-linear regime of interest. – 6 – . 0 1 0 . 1 11 0 k [ h / M p c ] L = 4 0 0 M p c / h L = 2 5 6 M p c / h
P(k) [(Mpc/h)3]
PLCDM/PHALOFIT-1 k [ h / M p c ]
Figure 2 . HaloFit accuracy for the power law ΛCDM simulations with box size L = 256 Mpc/ h (bluesquares) and L = 400 Mpc/ h box size (red circles) Left panel : linear power spectrum (dashed line) andnon-linear prediction from HaloFit (dotted line) plotted with the simulation results.
Right panel : fractionaldifference of simulation results and the HaloFit prediction. In all cases the errors correspond to 1- σ standarddeviation amongst 6 independent runs. use HaloFit to model only relative differences between simulations of the same L = 400 Mpc/ h size,resolution and initial phases as we shall now describe.In order to match the excess small scale final power in the f ( R ) model, we add an extra runningof the spectral tilt parameter to the ΛCDM initial power spectrum. Specifically, we assume a 3free-parameter initial power spectrum model: P i ( k ) = P (cid:18) kk p (cid:19) n + α ln( k/k p ) , (5.1)where P is the amplitude of the power spectrum at k p = 0 . h Mpc − and z = 0 without the effect ofthe transfer function and α is the running of the tilt, d ln P i d ln k = n + α ln( k/k p ) . (5.2)We therefore have 3 parameters p = { P , n , α } which specify the initial conditions.To find the best-fitting three parameters for a given model, we take the simulation results forthe power spectrum ratio (following method A), R sim ( k ) ≡ (cid:42) P sim f ( R ) ( k ) P simΛCDM ( k ) (cid:43) . (5.3)Next we use the HaloFit prescription P HF ( k ; p match ) for the non-linear matter power spectrum at z = 0 to find the best parameter set p match , by minimizing the χ given by χ = (cid:88) j σ R sim ( k j ) (cid:20) P HaloFit ( k j ; p match ) P HaloFit ( k j ; p ) − R sim ( k j ) (cid:21) , (5.4)where the sum runs over bins in k . Here, p describes the initial power spectrum used for the f ( R )simulations (see first line in Tab. 1). Finally we simulate a matched ΛCDM simulation with the sameinitial phases as the original but with a rescaling of the initial power P IC ( k ; p match ) = P IC , orig . ( k ) P i ( k ; p match ) P i ( k ; p ) . (5.5)In order to avoid confusion, we designate these ΛCDM simulations as “matched”; whereas ΛCDMwithout this modifier denotes the standard, power law, initial conditions (the one used in § f R | P / [10 (Mpc h − ) ] n α σ p − − − − − − Table 1 . Best-fit linear ΛCDM initial power spectrum parameters P , n , α (see text) that match the N-bodypower spectrum at z = 0 for each f ( R ) simulation. These values have been computed minimizing the χ of theN-body non-linear power spectrum and the non-linear HaloFit power spectrum for k ≤ k N / (cid:39) . h /Mpc.The same transfer function and the same cosmology is used in all cases. Also σ is shown for clarity. Pf(R)/Pmatched-1 k [ h / M p c ]
Figure 3 . Relative power spectrum offset between f ( R ) simulations and corresponding matched-ΛCDMsimulations for different | f R | values: for the chameleon simulations with | f R | = 10 − (red circles), 10 − (blue triangles), 10 − (green stars); and for the non-chameleon simulations with | f R | = 10 − (pink diamonds),10 − (orange hexagons), 10 − (cyan pentagons). The matched simulations make use of initial power spectrumconditions shown in Table 1 found by fitting relative deviations with HaloFit. The advantage of this matching method is that we only model relative deviations with the HaloFitprescription. Thus the cosmic variance of the original simulations scale out as do absolute errors inHaloFit, initial condition generators, resolution, etc. In Table 1 we show the best-fit values of the 3initial power spectrum parameters. We have only used R sim ( k ≤ . h/ Mpc) for the minimization.In Fig. 3 we show P f ( R ) /P matched −
1, where P f ( R ) is the power spectrum of the f ( R ) simulationsas before and P matched is the power spectrum of the matched ΛCDM simulations. Fig. 3 indicatesthat HaloFit is an excellent tool to predict relative differences in non-linear power spectra even for(some) non-standard ΛCDM models. We see that the differences between the matched ΛCDM and f ( R ) power spectra are up to ∼
4% in the range 0 . h/ Mpc < k < h/ Mpc, although for most of thescales and the cases are about 2-3%.As an aside, we can also test the absolute accuracy of HaloFit’s prediction for the power spectraof the matched models. Examples for different matched models are shown in Fig. 4. HaloFit producesa good fit compared with the sample variance errors for all k < . h . As in the pure ΛCDMcase, the sample variance at low k in the simulations is quite large. Deviations up to 10% for k < . h likewise appear due to the limited simulation resolution. Our modeling of relative effectseliminates these small differences. – 8 – . 0 1 0 . 1 1- 0 . 4- 0 . 20 . 00 . 20 . 4 k [ h / M p c ] Pmatched/PHALOFIT-1
Figure 4 . HaloFit accuracy for a representative set of matched ΛCDM models (matched to chameleon | f R | = 10 − , black squares; 10 − red circles; 10 − , blue triangles). HaloFit is accurate within the errors forall k < . h for these models which contain running of the tilt. Deviations from sample variance andresolution in the simulation seen here are largely absent in the relative matching technique shown in Fig. 3 .In reality one does not observe at a single z but in a wide z-range. As mentioned above it is notpossible to match the power spectrum at widely separated redshifts simultaneously and this featurecan provide observational signatures independent from the bispectrum. We can quantify this furtherby estimating over what redshift interval the power spectrum matching is expected to hold. Changesin the P f ( R ) /P ΛCDM were studied in detail by [10] and the excess evolves on the Hubble time scale.Therefore we generically expect that the matching evolves across a redshift interval of ∆ z = 1, i.e.no faster than any other aspect of the modeling. With the simulations of the matched ΛCDM models, we can now compare the bispectra for ΛCDMand f ( R ) models whose final power spectra match to a few percent.In Fig. 5 we show Q f ( R ) ( k ) /Q matched − k = 2 k = 0 . h/ Mpc (left panel) and for equilateraltriangle configuration (right panel), where Q f ( R ) ( k ) is the reduced bispectrum for f ( R ) simulations,and Q matched is the reduced bispectrum for the matched ΛCDM simulations. Red points show theratio for non-chameleon simulations whereas blue points for chameleons ones. Top panels correspondto | f R | = 10 − , middle panels to | f R | = 10 − and bottom panels to | f R | = 10 − . In particular,we see that for the chameleon and non-chameleon cases with | f R | = 10 − and 10 − the deviation isvery close to 0 ( (cid:46) | f R | = 10 − some differences appear: for the non-chameleon casethere is an excess of ∼
4% and for the chameleon case there is a deficit of ∼
4% in Q f ( R ) respectto Q matched , both within 5-6 σ . The value | f R | ∼ − is special in that it marks the onset of thechameleon mechanism in the largest structures in the simulations. The chameleon effect may havea small but measurable impact on Q in this transition region where the chameleon effect is presentfor some but not all structures. Analogous transient enhancements appear in the mass function [7].One should bear in mind though that this difference is of order the difference in the matched powerspectra which varies between the full and no-chameleon cases.Thus, for all values of f R deviations are below ∼ θ (cid:39) , π ) present higher deviations between different gravity modelsas has been observed in method A (Fig. 1) and predicted from theoretical models that followed thesame assumptions as adopted in method A [1, 2].Finally, we found that it is better to analyze the deviation between reduced bispectra Q ratherthan between bispectra B . This is because the power spectrum dependence is partially canceled in– 9 – . 2 0 . 4 0 . 6 0 . 8- 0 . 0 8- 0 . 0 40 . 0 00 . 0 40 . 0 8 k = 2 k = 0 . 4 h / M p c Qf(R)/Qmatched-1 | f
R 0 | = 1 0 - 4 f u l l f ( R ) n o n - c h a m e l e o n
E q u i l a t e r a l
Qf(R)/Qmatched-1 | f
R 0 | = 1 0 - 4
Qf(R)/Qmatched-1 | f
R 0 | = 1 0 - 5
Qf(R)/Qmatched-1 | f
R 0 | = 1 0 - 5 q / p Qf(R)/Qmatched-1 | f
R 0 | = 1 0 - 6 k [ h / M p c ]
Qf(R)/Qmatched-1 | f
R 0 | = 1 0 - 6
Figure 5 . Relative reduced bispectrum deviation for matched final power spectra (method B) between f ( R ) and ΛCDM simulations for k = 2 k = 0 . h/ Mpc as a function of θ (left panels), and equilateralconfiguration as a function of k (right panels), all at z = 0. Upper panels correspond to | f R | = 10 − ,middle panels to | f R | = 10 − and bottom panels to | f R | = 10 − . Red points correspond to non-chameleonsimulations, whereas blue points to full f ( R ) simulations. The error-bars are the 1 σ standard deviationamongst the ratio of 6 independent runs. Since we are taking the ratio between runs with the same initialphases, the cosmic variance errors are not present. the reduced bispectra. In spite of having run ΛCDM simulations to match the f ( R ) power spectra,several percent differences between f ( R ) and matched ΛCDM power spectra are still present (Fig. 3).These lead to higher deviations in B between the models (up to ∼
8% in some cases) than in Q .Thus, using Q instead of B is much more robust if we want to compare models with similar powerspectra. Of course, one should keep in mind that not all the P ( k )-dependence is cancelled when using Q , as evidenced by comparing with the results of method A. Strictly speaking, this is only true forequilateral configurations and up to tree-level in Eulerian perturbation theory.Finally, one may want to make a connection between these results and some analytic model,namely perturbation theory (PT). Since at tree level in PT the reduced bispectrum is independent ofthe power spectrum (at least for equilateral configuration), the differences observed between Fig. 1and 5 should be due to higher order corrections in ΛCDM. At 1-loop, corrections to the bispectrumcan be found in e.g., [14, 16]. One can see that the leading terms depend on the linear and one-looppower spectrum and weakly on cosmology and gravity through the standard tree level bispectrumkernel (see [1] for a modification of this kernel for some f ( R ) theories). Thus, a small modification ofthis formula could be expected due to f ( R ) gravity. However, this interpretation should be considered– 10 –ore in a qualitative way than in a strictly quantitative way. In fact one should take into account thatthe precision of 1-loop PT for the bispectrum is not much better than the other (phenomenological)analytic formulae [5]. As we have already mentioned, currently there is no analytic model that predictsthe bispectrum at scales of interest here with an accuracy of few percent. If the remaining ∼
4% deviation for | f R | = 10 − reflects gravity and not the residual mismatchin power spectra, then it is in principle measurable with large-volume surveys. In this work, con-sidering only the 6 runs of 400 Mpc/ h box-size and provided that h = 0 .
73, the total volumeis 6 × (0 . /h ) (cid:39) . We expect that future surveys will cover larger volumes: BOSS V ∼ /h ) , DES V ∼
10 (Gpc/h) or EUCLID V ∼
100 (Gpc /h ) . As the 6 runs have dif-ferent initial conditions we can use them to estimate the expected error on Q in the limit that it isdominated by cosmic variance. We have measured that the error in Q for our simulations at scales of k ∼ . h/ Mpc is about 5%. We assume that the variance scales as the inverse of the number of modes,and thus the standard deviation approximately scales as V / . Therefore, for a 10 Gpc survey theerror bars could, in principle, be as much as √ ∼ >
10 Gpc volume (e.g., DES, EUCLID) would yield an error on Q ∼
2% at thesescales. Since the expected deviation may be of order 4%, having smaller errors would help us toconfirm or discard possible deviation of the bispectrum due to modifications of gravity.On the other hand, we have analyzed the dark matter bispectrum which is not directly observable.In practice, sources of error that we have neglected here may appear: i) galaxy-surveys provide a biasedinformation about dark matter, ii) additional effects such as redshift distortions change the observedbispectrum (in fact we expect modified gravity to affect redshift distortions more strongly than thedensity field itself). Also as we go to higher z , we expect less deviations at a given scale. Conversely,the matched power spectra at z = 0 would become mismatched and provide other observable effects.The results from Fig. 5 provide another important result. We have seen that two f ( R ) theories ofgravity with indistinguishable non-linear dark matter power spectrum, have very similar and possiblyindistinguishable dark matter bispectra. This opens up the possibility of using these two statisticsto break degeneracies in the galaxy bias in a way that is robust to the assumptions about the trueunderlying gravity model.In fact the f ( R ) effects on the power spectrum are at the 20-50% level. A modification of galaxybias achieving similar effects would likely affect the reduced bispectrum at least at the 10% level (forexample, a linear bias term affect the power spectrum ∝ b and the reduced bispectrum ∝ /b ),significantly larger than the f ( R ) effects on the reduced bispectrum. In this work we have analyzed the deviations in the reduced bispectrum produced by a modificationof gravity, specifically the f ( R ) class of models, both with and without the chameleon mechanism. Inorder to do that, we make use of a suite of f ( R ) and ΛCDM simulations. We have proceeded in twodifferent ways to analyze the bispectrum deviation from these simulations, methods A and B, whichdiffer in whether the initial or final power spectra of the two cosmologies are set equal.Method A compares the bispectrum output of f ( R ) and ΛCDM N-body simulations with thesame initial power spectrum. Fig. 1 shows the bispectrum deviation obtained using this method.We observe a considerable deviation (up to 10 − f ( R )models and the ΛCDM one. Such differences in the bispectrum could be easily detected by surveyscovering volumes > | f R | and for squeezed triangle configurations. In this method, both ΛCDM and f ( R ) gravity runs start from the same initial δ k values. Because of that, the different evolution ofthe gravity models naturally leads to different power spectra (as was observed in [10]) and also todifferent bispectra. This way of proceeding is equivalent to the theoretical works of Bernardeau & Baryon Oscillation Spectroscopic Survey Dark Energy Survey – 11 –rax [1], Borisov & Jain [2]. In order to explain discrepancies between the matter power spectrum in f ( R ) and the observed galaxy power spectrum, one could invoke a scale-dependent galaxy bias. Sincegalaxy bias enters into the reduced galaxy bispectrum in a different way than in the power spectrum,bispectrum measurements can in principle close this loophole.Alternatively, the large power spectrum differences can be eliminated by changing the shape ofthe initial power spectrum to instead match the final dark matter power spectrum at z = 0. This isat the base of method B. In this method, we compute the bispectrum deviation between a ΛCDMand a f ( R ) model, both with the same final power spectra. Thus, we simulate a ΛCDM modelwith certain initial power spectrum parameters (summarized in Table 1) that are adjusted to bestmatch the f ( R ) power spectrum at z = 0. From the simulations outputs we compute power spectraand bispectra. For the power spectra, residual differences are never higher than 4% in the range0 . h/ Mpc < k < h/ Mpc.Likewise the differences in the reduced bispectrum are also smaller in the matched comparison.For the | f R | = 10 − and 10 − cases, the Q deviation is consistent with 0 within 1 σ . For | f R | = 10 − deviations in Q at most reach the 4% level with 5 − σ significance. These deviations are potentially asignature of the onset of the chameleon mechanism in the largest structures in the Universe. Howevergiven that this is the same order as the power spectrum difference it is unclear whether these differencesindicate power-spectrum-independent modified gravity effects or that the two power spectra are notperfectly matched. In the former case, larger surveys like EUCLID will allow for a measurement ofthe bispectrum with enough precision to obtain a > σ significance, even when exactly matching thepower spectra.On the other hand, the effect of deviations from GR gravity on the reduced bispectrum are weakcompared to those on the power spectrum (at least for the cases considered here), opening up thepossibility of breaking the galaxy-bias degeneracy. In fact the effect of galaxy bias is expected to bedifferent in the power spectrum and in the bispectrum, which is why, in the context of GR gravity,the bispectrum is used to constrain galaxy bias. While the shape of the non-linear power spectrumseems to carry information about the underlying gravity model, one may always argue that a shapeof the evolved power spectrum not compatible with the GR predictions could be due to biasing. Forthe cases we have considered here, the dependence of the reduced bispectrum on deviations fromGR is weaker than the effects of bias modifications necessary to explain the deviations in the powerspectrum. While we have only studied f ( R ) models here, there is no apparent reason why this resultshould be specific to f ( R ). Hence, if our findings were to remain qualitatively true for other gravitymodifications, this would confirm the usefulness of employing the reduced bispectrum together withthe power spectrum to constrain bias parameters. We thank Christian Wagner for useful discussions and help with N-body simulations. HGM is sup-ported by a CSIC JAE grant, and thanks the Kavli Institute for Cosmological Physics (KICP) atUniversity of Chicago for hospitality. Part of this work stemmed from discussions at the Centrode Ciencias de Benasque Pedro Pascual. WH is supported by the KICP under NSF contract PHY-0114422, DOE contract DE-FG02-90ER-40560 and the Packard Foundation. LV and RJ are supportedby MICINN grant AYA2008-03531. LV acknowledges support from grant FP7 ERC- IDEAS Phys.LSS240117. FS is supported by the Gordon and Betty Moore Foundation at Caltech.
References [1] Bernardeau, F., & Brax, P. 2011, arXiv:1102.1907[2] Borisov, A., & Jain, B. 2009,
Phys. Rev. D , 79, 103506[3] Eisenstein, D. J., & Hu, W. 1998,
ApJ , 496, 605[4] Fry, J. N. 1994, Physical Review Letters, 73, 215[5] Gil-Mar´ın, H., Wagner, C., Fragkoudi, F., Jimenez, R., & Verde, L. 2011, arXiv:1111.4477 – 12 –
6] Hu, W., & Sawicki, I. 2007,
Phys. Rev. D , 76, 064004[7] Li, Y., Hu, W. arXiv:1107.5120[8] O’Shea, B. W., Bryan, G., Bordner, J., Norman, M. L., Abel, T., Harkness, R., & Kritsuk, A. 2004,arXiv:astro-ph/0403044[9] Oyaizu, H. 2008,
Phys. Rev. D , 78, 123523[10] Oyaizu, H., Lima, M., & Hu, W. 2008,
Phys. Rev. D , 78, 123524[11] Perlmutter, S., et al. 1999,
ApJ , 517, 565[12] Riess, A. G., et al. 1998, AJ , 116, 1009[13] Schmidt, F., Lima, M., Oyaizu, H., & Hu, W. 2009, Phys. Rev. D , 79, 083518[14] Scoccimarro, R. 1997,
ApJ , 487, 1[15] Sealfon, C., Verde, L., & Jimenez, R. 2005,
Phys. Rev. D , 71, 083004[16] Sefusatti, E. 2009,
Phys. Rev. D , 80, 123002[17] Smith, R. E., et al. 2003,
MNRAS , 341, 1311[18] Song, Y. S., Peiris, H., Hu, W. 2007,
Phys. Rev. D , 76, 063517[19] Sotiriou, T. P., & Faraoni, V. 2010, Reviews of Modern Physics, 82, 451[20] Verde, L., et al. 2002,
MNRAS , 335, 432[21] Zhao, G.-B., Li, B. and Koyama, K. 2010,
Phys. Rev. D , 83, 044007, 83, 044007