Astronomy&Astrophysics587,A132(2016) (cid:13)cESO2016 March3,2016 N-body simulations of γ gravity MarceloVargasdosSantos1,2,(cid:63),HansA.Winther3,(cid:63)(cid:63),DavidF.Mota2,(cid:63)(cid:63)(cid:63),andIoavWaga1,(cid:63)(cid:63)(cid:63)(cid:63) 1 InstitutodeFísica,UniversidadeFederaldoRiodeJaneiroC.P.68528,CEP21941-972,RiodeJaneiro,RJ,Brazil 2 InstituteofTheoreticalAstrophysics,UniversityofOslo,Postboks1029,0315Oslo,Norway 3 Astrophysics,UniversityofOxford,DWB,KebleRoad,Oxford,OX13RH,UK March3,2016 6 1 ABSTRACT 0 Wehaveinvestigatedstructureformationintheγ gravity f(R)modelwithN-bodysimulations.Theγ gravitymodelisaproposal 2 which,unlikeotherviable f(R)models,notonlychangesthegravitationaldynamics,butcaninprinciplealsohavesignaturesatthe r backgroundlevelthataredifferentfromthoseobtainedinΛCDM(Cosmologicalconstant,ColdDarkMatter).Theaimofthispaper a istostudythenonlinearregimeofthemodelinthecasewhere,atlatetimes,thebackgrounddiffersfromΛCDM.Wequantifythe M signaturesproducedonthepowerspectrum,thehalomassfunction,andthedensityandvelocityprofiles.Toappreciatethefeatures ofthemodel,wehavecomparedittoΛCDMandtheHu-Sawicki f(R)models.Fortheconsideredsetofparameterswefindthatthe 2 screeningmechanismisineffective,whichgivesrisetodeviationsinthehalomassfunctionthatdisagreewithobservations.Thisdoes notruleoutthemodelperse,butrequireschoicesofparameterssuchthat|fR0|ismuchsmaller,whichwouldimplythatitscosmic ] expansionhistorycannotbedistinguishedfromΛCDMatthebackgroundlevel. O C Keywords. Gravitation–Cosmology:darkenergy–Galaxies:clusters:general–Cosmology:large-scalestructureofUniverse– Galaxies:halos . h p -1. Introduction 1980),isdrivenbyatermofthetype f(R)=αR2 (α>0)andis o stillinexcellentaccordancewithobservations(Adeetal.2014). r Sincethediscoveryin1998(Riessetal.1998;Perlmutteretal. t More recently, the idea of an acceleration driven by late-time s1999)thattheUniverseisspeedingupinsteadofslowingdown a curvaturehasalsobeenexploredinCapozzielloetal.(2003)and (aswouldbeexpectedifgravityisalwaysattractive),consider- [ Carrolletal.(2004).Theseauthorsconsideredatheoryinwhich ableefforthasbeendevotedtounderstandingthephysicalmech- f(R)=−αR−n(n>0andα>0).However,thesemodelsdonot 2anism behind this cosmic acceleration. The two main theoreti- have a regular matter-dominated era and are incompatible with vcal approaches considered in the literature to explain this phe- structureformation(Adeetal.2014). 3 nomena are (1) to assume the existence of a new component 3with a sufficiently negative pressure (p < −ρ/3), generically Tobuildacosmologicallyviable f(R)theory,somestability 4 conditions have to be satisfied (Pogosian & Silvestri 2008): (a) 5denoted dark energy, and (2) to consider that general relativity f ≡d2f/dR2 >0(notachyons);(b)1+f ≡1+df/dR>0[the 01.0hcdouaafsrtsevttaaoittsuebrEeewin(mmstoeod≡idnifii’fisepedcdoa/sgtρmrlaaovrligotey=g)is.c−caTal1hl.ecesoH,nsoiosmwtra,penmvlteeos(rtΛr,edi)naawrcskciptuihetrneaaetenorlgfyey,iqtasuctaavntlieodorwiyn- ecaRffhnRadenclgtiiemvesRi→ggrn∞av(fniRtoat=giohno0aslt(sGc)]oe;nn(scet)raanaltf,tRe(rGeliaenfftfliva=ittiyoGRniN,sl/ir(me1cR+o→v∞efRref))(dRda)o/teRsea=nroly0t 6 DE DE DE times) ; (d) |f | is small at recent times, to satisfy solar system good accordance with current observations, Λ has some theo- R 1 and galactic scale constraints. In addition to these conditions, retical difficulties such as its tiny value as compared with the- : therearesomedesirablecharacteristicsthataviablecosmologi- voretical predictions of the vacuum energy density, the cosmic calmodelhastosatisfy(Amendolaetal.2007).Itshouldhavea Xicoincidenceproblem,andrelatedfine-tuning.Thissituationhas radiation-dominatederaatearlytimesandasaddle-pointmatter- motivatedthesearchforalternativeslikemodified-gravitytheo- r dominated phase followed by an accelerated expansion as a fi- aries.Thesimplestmodified-gravitycandidatesaretheso-called nal attractor. By using the parameters m¯ ≡ Rf /(1+ f ) and f(R)-theories,inwhichtheLagrangiandensityL = R+ f(R)is r ≡ −R(1+ f )/(R+ f), it can be shown that,RaRn earlyRmatter- anonlinearfunctionoftheRicciscalarR. R dominatedepochoftheUniversecanbeachievedifm¯(r≈−1)≈ Asiswellknown,metric f(R)-theoriescanbethoughtofasa 0+ and m¯/r(r ≈ −1) > −1. Furthermore, a necessary condition specialcaseofascalar-tensortheory;aBrans-Dickemodelwith foralate-timeacceleratedattractoris0<m¯(r≈−2)≤1. acouplingconstantω =0.Anacceleratedexpansionappears BD Thereareviable f(R)gravitytheoriesthatsatisfyallthecri- naturallyinthesetheories.Theveryfirstinflationarymodel,pro- teriamentionedabove(Hu&Sawicki2007;Starobinsky2007; posedbyStarobinskymorethanthreedecadesago(Starobinsky Appleby & Battye 2007; Cognola et al. 2008; Linder 2009; O’Dwyeretal.2013).However,thereisagenericdifficultyfrom (cid:63) [email protected] (cid:63)(cid:63) [email protected] which all these “viable” f(R) theories (Thongkool et al. 2009) (cid:63)(cid:63)(cid:63) [email protected] suffer: the curvature singularity in cosmic evolution at a finite (cid:63)(cid:63)(cid:63)(cid:63) [email protected] redshift (Frolov 2008). It can be shown that this type of sin- Articlenumber,page1of10 A&Aproofs:manuscriptno.main gularity problem can be cured, for instance, by adding a high- 2. γgravityreview curvature term proportional to R2 (Appleby et al. 2010) to the Weinvestigatespatiallyflatcosmologicalmodelsinthecontext density Lagrangian. Therefore, it is not possible to have cos- ofγgravity(O’Dwyeretal.2013),aviable f(R)theorydefined micaccelerationwithatotallyconsistent f(R)theorymodifying bythefollowingansatz: gravity only at low curvatures. We remark that we do not ad- dress this problem here and only consider modifications at low αR (cid:34)1 (cid:32) R (cid:33)n(cid:35) curvatures. f(R)=− ∗γ , , (1) n n R Weconsiderthespecificcaseofaviable f(R)theorycalled ∗ γabglerafv(iRty)t(hOe’oDriwesy,esrtreutcatul.re2f0o1r3m).atGioenneimricpaolsley,sisnucahlmstorostngalclovni-- where γ(n,x) = (cid:82) xtn−1e−tdt is the incomplete Γ-function and 0 straintsontheparametersofthemodelsthattheeffectiveequa- α, n and R∗ are free positive constants. In reality, γ gravity can tion of state parameter cannot be distinguished from that of a be thought of as a simple generalization of exponential gravity cosmologicalconstant.Inγgravitythesteepdependenceonthe (Linder2009) RicciscalarRfacilitatestheagreementwithstructureformation. O’Dwyer et al. (2013) showed that, in principle, the parameter f(R)=−αR∗(1−e−R/R∗), (2) thatcontrolsthesteepnessinγgravityallowsmeasurabledevia- tionsfromΛCDM(Cosmologicalconstant,ColdDarkMatter)at obtained by fixing n = 1 in Eq.(1). We emphasize that γ grav- ity can satisfy all the stability and viability conditions. As dis- bothlinearperturbationandbackgroundlevels,whilestillcom- cussedinO’Dwyeretal.(2013),forfixedn,thereisaminimum patiblewithbothcurrentobservations.Themaingoalofthispa- peristostudytheeffectsofγgravityonthestructureformation value(αmin)oftheparameterαsuchthatforvaluesα > αmin a late-timeacceleratedattractorisachieved.Weconsiderthiscase at nonlinear scales for choices of parameters where the model hasobservablesignatures1onthebackgroundexpansionhistory throughout.FromEq.(1)weobtainthefollowingderivatives: ofourUniverse. Wegoonestepfurtherandanalyzethenonlinearevolutionof fR =−αe−(cid:16)RR∗(cid:17)n, (3) sthtriuscptuapreesrcisombapsuetdedonfriosmansluigmhetrmicaoldsifiimcautliaotnioonfs.ITSIhSec(Lodlienathreast fRR = αRn(cid:32)RR (cid:33)ne−(cid:16)RR∗(cid:17)n. (4) ∗ et al. 2014), which in turn is a modification of the RAMSES hy- drodynamic N-body code (Teyssier 2002). See also Zhao et al. WenotefromEq.(3)thatwithincreasingn,thesteepnessofthe (2011);Lietal.(2012);Oyaizuetal.(2008);Motaetal.(2008); f(R) function increases. Higher n means smaller |f |, and the R0 Li et al. (2011); Boehmer et al. (2010); Zumalacarregui et al. departuresfromGRwillbesmalleraccordingly. (2013);Puchweinetal.(2013)forothercodesthathaveimple- Althoughthereisnocosmologicalconstant, f(0)=0,itfol- mentedandperformedsimulationsof f(R)gravityand(Winther lowsfromEq.(1)thatGRwithΛisrecoveredathighcurvatures. etal.2015)forarecentcodecomparisonbetweenthesecodes. Therefore,forR(cid:29)R themodelsbehavelikeΛCDM.Sincewe ∗ We only consider the modifications made to implement the aremainlyinterestedinphenomenathatoccurredafterthebegin- γ gravity field equations in the modified gravity part of the N- ningofthematter-dominatedera,weneglectradiationandwrite bodycode.Formoredetailsontheimplementationofthescalar the effective cosmological constant (the cosmological constant fieldsandothertechnicalitieswerefertoLlinaresetal.(2014). ofthereferenceΛCDMmodel)as For this purpose, we focus on simple observables such as the matterpowerspectrum,halomassfunction,densityprofiles,and Λ˜ = αR∗Γ(1/n)=3H˜2(1−Ω˜ ). (5) velocity profiles to investigate modified gravity signatures that 2n 0 m0 were previously studied in Li et al. (2012); Hammami et al. (2015); Schmidt et al. (2009); Gronke et al. (2015); Lombriser In the equation above, Ω˜m0 denotes the present value of the etal.(2012,2014);Wintheretal.(2012);Terukina&Yamamoto matter density parameter that a ΛCDM model would have if it (2012); Shi et al. (2015); Pujol et al. (2014); Li et al. (2013); hadthesamematterdensitytoday(ρ¯m0)asthemodifiedgravity Hellwingetal.(2013);Heetal.(2013,2015);Schmidt(2010); f(R)model. H˜0 representstheHubbleconstantinthereference Braxetal.(2013);Tessoreetal.(2015);Gronkeetal.(2015);Li ΛCDMmodel.Therefore,wehavem2 ≡8πGρ¯m0/3=Ω˜m0H˜02 = etal.(2013)thatcanalsobeobserved(Bertietal.2015;Zivick Ω H2,whereΩ andH arethepresentvalueofthematteren- m0 0 m0 0 et al. 2015; Mak et al. 2012; Cai et al. 2014, 2015; Song et al. ergydensityparameterandHubbleparameterinthe f(R)model, 2015; Wilcox et al. 2015; Jain & Khoury 2010; Schmidt 2010; respectively.ItisusefultorewriteR as ∗ Hellwingetal.2014). R 6nd Thispaperisorganizedasfollows:inSect.2werevisittheγ ∗ = , (6) gravity model and show the main properties of the background m2 αΓ(1/n) andlinearperturbationevolution.Section3detailsthedynamics equationsofscalaronfield(fR)andparticlemovementequations whered =(1−Ω˜m0)/Ω˜m0.Tocomputethebackgroundevolution, forγgravity,whichmustbesolvedbyourcodeduringthesim- westartfromthe f(R)fieldequationforaFLRWmetric ulations. The method for solving these equations is briefly ex- Rf − f plained in Sect. 4, and the code is tested in Sect. 5. Finally our H2(1+ f +R(cid:48)f )− R =m2e−3y, (7) R RR resultsareshowninSect.7,andweconcludeinSect.8. 6 where(cid:48) ≡ d/dy(y = lna), H ≡ a˙/aistheHubbleparameter(a 1 Withobservablesignatureswemeanthattheequationofstatediffers dotdenotesthederivativewithrespecttocosmictime),whichis enoughfromtheΛCDMvalueofw = −1atlowredshiftsthatsucha relatedtoRby deviation could be detected by near-future experiments like WFIRST Spergeletal.(2015),forexample. R=12H2+6HH(cid:48). (8) Articlenumber,page2of10 M.VargasdosSantosetal.:N-bodySimulationsofγGravity Ω˜ n α m0 0.267 2 1.05 0.267 2 1.18 0.267 2 1.5 Table1.Overviewofthemodelparametersforγgravity. Fig.2.Fractionaldifferenceinthematterpowerspectrumwithrespect toΛCDMfordifferentvaluesofnandα,asindicatedinTable1.These resultsareusedtocomparewiththenonlinearpowerspectrumfromour numericalsimulations. In GR, f = Q = 0, and there is no scale dependence for RR thedensitycontrastinthelinearregime.ForΛCDMthegrowing modecanbeexpressedintermsofhypergeometricfunction F Fig.1.Effectiveequation-of-stateparameterw asafunctionofred- 2 1 DE as(Silveira&Waga1994) shiftzfortheparametersgiveninTable1.Thestrongestdeviationfrom −1islowerthan4%. (cid:34) (cid:35) 1 11 δ+ ∝e−y2F1 ,1, ,−e3yd . (16) Tosolvetheseequations,weintroducethenewvariables 3 6 H2 We solved Eq. (14) numerically and obtained the growing x (y)= −e−3y−d, (9) 1 m2 mode for the γ gravity. By using (16), we then obtained the R fractional change in the matter power spectrum P(k) relative to x2(y)= m2 −3e−3y−12(d+x1). (10) ΛCDM. Figure 2 shows ∆Pk/PΛ at z = 0 for the three choices ofparametersshowninTable1. Withthesedefinitionsweobtain x x(cid:48)(y)= 2, (11) 3. N-bodyequations 1 3 R(cid:48) f(R)modelsareequivalenttoascalar-tensortheory(Braxetal. x2(cid:48)(y)= m2 +9e−3y−4x2, (12) 2008), where the first derivative of the f(R) function, fR. This fieldpropagatesaccordingtheequation whereR(cid:48) isgivenbyEq.(7).Itisstraightforwardtoverifythat, asdefined, x and x arealwayszeroduringtheΛCDMphase. Weherefocu1sonthe2threecasessummarizedinTable1. (cid:3)f = ∂Veff = (1− fR)R+2f +κ2T, (17) Furthermore, in terms of x1(y) and x2(y), the effective dark R ∂fR 3 energyequationofstate(w )isgivenby, DE where κ2 = 8πG/c4 and T is the trace of energy-momentum w =−1− 1 x2 . (13) tensor, T = gµνTµν. In the quasi-static limit (see, e.g., Noller DE 9x +d etal.(2014);Boseetal.(2015);Llinares&Mota(2014,2013)) 1 thisequationbecomes Fortheconsideredmodels,theevolutionofw asafunctionof DE theredshiftzisshowninFig.1. 1 R−R(a) ∇2f = −m2a−3δ , (18) For a general f(R) model the differential equation for the a2 R 3 m matterdensitycontrast(δ )inthelinearregimeforsubhorizon m scalesisgivenby(Pogosian&Silvestri2008;Zhang2006;dela where Cruz-Dombrizetal.2008) R(a) =3(a−3+4d)+∆ (a), (19) δ(cid:48)(cid:48)+(cid:32)2+ H(cid:48)(cid:33)δ(cid:48) − 1−2Q 3H02Ω˜m0 e−3yδ =0, (14) m2 R m H m 2−3QH2(1+ f ) m and ∆ (a) = x (a)+12x (a). The Ricciscalar R in functionof R R 2 1 f isgivenbyinvertingEq.(3) R where (cid:32) α (cid:33)1/n Q(k,y)=− 2fRRc2k2 . (15) R=R∗log |f | . (20) (1+ f )e2y R R Articlenumber,page3of10 A&Aproofs:manuscriptno.main The geodesic equation, needed to update the particle positions, 4. ImplementationintheISIScode reads Implementingscalar-tensortheoriesofgravityinN-bodycodeis (cid:32) (cid:33) ratherstraightforwardbecausethescalar-tensortheoriesallcon- 1 f x¨+2Hx˙=− ∇ Φ− R , (21) tributeasafifthforceandbecauseRAMSES,whichISISisbased a2 2 on, has been widely used, thoroughly tested, and optimized. In thissectionwedescribehowtheequationsweneedtosolveare whereΦisthenewtonianpotential,whichthedynamicsisgiven implementedinISIS.Formoredetailssee(Llinaresetal.2014). bythePoissonequation To solve for f directly is not numerically stable since the R solution can potentially vary over several orders of magnitude 3m2δ when going from deep voids to massive clusters in our simula- ∇2Φ= m. (22) tion. We therefore introduce a field redefinition |f | = A(u) → 2 a ∇|f | = b(u)∇u where b(u) = dA(u)/du. The Rgeneral field R equation for f discretized on a grid with the field-redefinition When implementing these equations in the N-body code, we f˜ ≡ −a2f =R A(u), where u is the field we solve for, can be R R needtorewritethemincode-unitsgivenby writtenas Φa2 L(ui,j,k)=∇code·[b(u)∇codeu]i,j,k+Ωm0(H0B0)2× (31) x˜= x/B0, Φ˜ = (H0B0)2, ×a4R∗ log(cid:34) αa2 (cid:35)1/n−a(δ ) −a−a4(cid:34)4d+ ∆R(a)(cid:35), dt˜= H0dt, ∇ = B .∇. (23) 3m2 A(ui,j,k) m i,j,k 3  a2 code 0 where HereB0isthesizeofthesimulationbox.Intermsof f˜R =−a2fR, ∇ ·[b(u)∇ u] = (32) code code i,j,k theevolutionequationsbecomes bi+1/2,j,k(ui+1,j,k−ui,j,k)−bi−1/2,j,k(ui,j,k−ui−1,j,k)+ h2 d2x˜ 1 dt˜2 =−∇codeΦ˜ − 2(B H )2∇codef˜R, (24) + bi,j+1/2,k(ui,j+1,k−ui,j,k)−bi,j−1/2,k(ui,j,k−ui,j−1,k)+ 0 0 h2 3 ∇2codeΦ˜ = 2Ωm0aδm, (25) + bi,j,k+1/2(ui,j,k+1−ui,j,k)−bi,j,k−1/2(ui,j,k−ui,j,k−1), h2 ∇2 f˜ =Ω (H B )2a4× (26) code R m0 0 0 ×−3Rm∗2 log(cid:32)αfa˜2(cid:33)1/n+(cid:34)a−3+4d+ ∆R3(a)(cid:35)+a−3δm. wanhderwehherisetwheeghraivdesdpeaficinnegdabn(du)b≡i±1/d2Ad,(uju,k).≡ 21(b(ui±1,j,k)+b(ui,j,k)) R TheequationsaresolvedusingNewton-Gauss-Seidelrelax- ation(withmultigridacceleration).Themethodconsistsofgoing Thesearetheonlyequationsweneedtoimplementandsolvein throughthegridandupdatingthesolutionusing theN-bodycode. L(u ) For comparison we also need the linearized field equation. unew =u − i,j,k . (33) Simulationswiththisequationcomparedtothefull fR equation i,j,k i,j,k ∂L(ui,j,k)/∂ui,j,k isagoodmeasureoftheamountofscreeningthattakesplacein Forthiswealsoneed∂L(u )/∂u ,whichisgivenby themodel.Thelinearized f equationissimply i,j,k i,j,k R ∂L(u ) ∂∇ ·[b(u)∇ u] i,j,k = code code i,j,k+ (34) 1 ∇2δf =m2(a)δf −m2a−3δ , (27) ∂ui,j,k ∂ui,j,k a2 R φ R m −Ω (H B )2a4R∗ b(ui,j,k) log(cid:34) α (cid:35)1/n−1, whereδf = f − f (a)andm2(a)= 1 .Incodeunits,taking m0 0 0 3m2 A(ui,j,k) A(ui,j,k)  R R R φ 3fRR(a) u=− δfRa2 ,weobtain where 2(H0B0)2 ∂∇ ·[b(u)∇ u] Ω a code code i,j,k = (35) ∇2 u=[m (a)aB ]2u+δ m0 , (28) ∂ui,j,k code φ 0 m 2 − bi+1/2,j,k+bi−1/2,j,k − bi,j+1/2,k−bi,j−1/2,k+ h2 h2 andthegeodesicequationbecomes − bi,j,k+1/2+bi,j,k−1/2 + 1c (cid:34)ui+1,j,k+ui−1,j,k−2ui,j,k+ h2 2 i,j,k h2 d2x˜ dt˜2 =−∇codeΦ˜ −∇codeu. (29) +ui,j+1,k+ui,j−1,k−2ui,j,k + ui,j,k+1+ui,j,k−1−2ui,j,k(cid:35), h2 h2 Wehave and c = c(u ), where c(u) = db(u). Some problems arise i,j,k i,j,k du a2(H B )2R(a)(cid:34) R (cid:35)n related to how boundary conditions are handled on the refined m2(a)a2B2 = 0 0 ∗ e[R(a)/R∗]n. (30) gridsinthecode.ThisisdiscussedinLlinaresetal.(2014);Li φ 0 3αn H2 R(a) 0 etal.(2012). Articlenumber,page4of10 M.VargasdosSantosetal.:N-bodySimulationsofγGravity Our solver needs a starting guess, and for this we use the whereδ= ρin −1characterizesthedensitycontrastbetweenthe cosmologicalbackgroundsolution,thatis,thesolutionwewould insideandoρouuttsideofthesphere,ρ¯ isthemeandensity,Risthe expectwhentherearenomattersources, radiusofthesphere,and B thesizeofthebox.Thevalueofδ 0 chosen for the test is 5000. For the f(R) test we used R = 25 |fR(a)|=αe−[R(a)/R∗]n →u= A−1(|fR(a)|). (36) Mpc/hand B0 = 250Mpc/h.Asphericalsymmetricconfigura- tion is effectively one-dimensional, therefore the field equation This is only needed when the simulation is stared as otherwise reducestoanODE,whichwesolvedbyusingMathematicaand wecanusetheoldsolutionasourguess.Forγgravityourchoice usedthistocomparewith. forAandtherelatedexpressionsforbare Forthesecondtestweanalyticallycomputedadensityfield ρ(x,y,z) ≡ ρ(x) (i.e., a 1D configuration), using the field equa- A(u)≡ f˜R(u)=αa2e−eu, (37) tion,sothatthesolutionisgivenbyasine:u∝2+sin(2πx). b(u)= dA(u) =−αa2e−eueu, (38) Figure 3 shows the result of both these tests and the differ- du entcolorsdepictthedifferentrefinementlevels.Weseethatthe db(u) curvesaresmoothwhengoingfromoneleveltoanother,which c(u)= =αa2e−eu(eu−1)eu, (39) du demonstratesthatboundaryconditionsarehandledproperly.The tests were performed using the serial version of the code, and andthebackgroundvalueofuis bothtestsgivetheexpectedresults,whichdemonstratesthatthe   codeworksproperly. u= A−1(f˜R(a))=nloga−3+R4d/(+3H∆2R)(a)/3, (40) ∗ 0 6. Simulations intermsofuwehave Theγgravitysimulationswererunusing5123 darkmatterpar- L(u)=∇[b(u)∇u]+Ωm0(H0B0)2× (41) ticleswithaboxsizeof B0 = 250Mpc/h,andwerefinedcells (cid:34)a4R (cid:32) ∆ (a)(cid:33)(cid:35) whenever it had more than eight particles in it. The highest re- × ∗eu/n−aδ −a−a3 4d+ R , finement level reached in our simulation was eight. The back- 3m2 m 3 groundcosmologyusedforthesimulationwascomputedusing ∂ L(u)= ∂ ∇[b(u)∇u]+Ω (H B )2(cid:32)a4R∗eu/n(cid:33). (42) Eq.(7)withh=0.71,ΩΛ =0.733andΩm0 =0.267. ∂u ∂u m0 0 0 3m2 n ThemodelparameterswerepresentedinSect.2,withsome plotsofvariousbackgroundquantities.Wealsocomparedthese When the fifth force is implemented, we simply replace it with simulations with simulations of the Hu-Sawicki f(R) model aneffectiveforceFeff thatincludestheeffectsofmodifiedgrav- fromLlinaresetal.(2014). ity wherever the code normally works with the gravitational force,F , N 7. Results Feff = FN +Fφ. (43) 7.1. Powerspectrum The expression for the fifth force in code-units is given in The nonlinear matter power spectrum is an important observ- Eq.(24),andweusethesamefive-pointstencilasRAMSESuses ableandcouldbeusedtodistinguishamongdifferentmodelsof tocomputethegravitationalforce∇Φ˜ tocomputethefifthforce structureformation.Asweshowedabove,γ gravitycanhavea ∇f˜R. strong effect on the growth rate of the linear perturbations. We expect these signatures to be detectable in the nonlinear matter powerspectrum. 5. TestsoftheN-bodysolver To compute the power spectrum we used a public code, To verify that the solver is implemented correctly, we tested it POWMES (Colombi & Novikov2011), which uses folding meth- severaltimes.Wepresentsomeofthesetestsinthissection.The odstocomputethepowerspectrum. density is provided to the code through a distribution of parti- Figure 4 displays the difference of the matter power cles. The density estimation (CIC) and refinement criteria are spectrum with respect to ΛCDM, defined as ∆P/PΛCDM ≡ thesameasthoseusedforthecosmologicalsimulations.Wealso P(k)/PΛCDM(k)−1forγgravitytogetherwiththecorresponding have the option to set the density field analytically in the code. predictionsfromlinearperturbationstheoryandtheHu-Sawicki To test that the treatment of the boundary of the refinement is f(R)modelforcomparison.Wefocusonthepresent-dayepoch, correct, we included two levels of refinement. The model pa- whichcorrespondstoz=0. rametersusedforthetestsarethesameasmentionedpreviously Figure 4 shows that the differences from ΛCDM are lower inTable1. than5−10%forallofourrunsintherange0.05hMpc−1 (cid:46)k (cid:46) ThefirsttestcasecorrespondstoasphereofradiusRofcon- 0.1hMpc−1)andinagreementwiththelinearperturbationresult stantdensitylocatedinthecenteroftheboxandembeddedina seeninFig.(2).ThedeviationfromΛCDMapproacheszerofor uniformbackground: larger scales. This is because the range of fifth force is smaller than the horizon, therefore the modifications of gravity are not ρ(r)= 1+(143π+δδ(cid:16))BρR0(cid:17)3, r<R , (44) feeevffleetrc,otn∆oPfth/tehPeΛlaCfirDgftMehstcfoosrnccateilneausc.etsOs,tnoansgdmrowawlelewsr,eientholankrlgiinenerγadregvsrciaaavtliieotsyn,,sti.hnHecofouwnl--l  1+ 4πρδ(cid:16)R(cid:17)3, r>R tmernaehsctahntaocnetimhsmeenHits.u-iTnShapiwslaiiycskoyannmthioneddseiecl,aswtciaohlneerset,hwtahtheitchcheharsmecdreuelceeeonsnintgshcermepeeoncwihneagr- 3 B0 Articlenumber,page5of10 A&Aproofs:manuscriptno.main Fig. 3. Scalar field from tests. Different colors depict the different refinement levels. The left panel shows the field from a spherical density distribution,therightsideshowsa1Dsinefieldobtainedinthesecondtest. Fig. 5. Fractional difference of the halo mass function for γ gravity Fig. 4. Fractional deviation of the dark matter power spectrum with (data points) with respect to the ΛCDM model. For comparison we respect to the ΛCDM model. The linear predictions of γ gravity are alsoshowtheresultsfromsimulationswiththeHu-Sawicky f(R)model represented by dashed lines. For comparison we also show the re- with|f |={10−4,10−5,10−6}(dashedlines)fromLlinaresetal.(2014). R0 sults from simulations with the Hu-Sawicky f(R) model with |f | = R0 {10−4,10−5,10−6}fromLlinaresetal.(2014). nismdoesnotworkverywellinoursimulation.Wediscussthis the signatures are very similar to the |f | (cid:38) 10−5 Hu-Sawicki inmoredetailinSect.7.4. model.Forthemostmassivehalos(M/MR0 (cid:38) 1015)theeffectof (cid:12) screeningismuchmorepronouncedinHu-Sawickimodels,and themass-functionapproachesΛCDM.Thisdoesnotoccurforγ 7.2. Halomassfunction gravityandisagainanindicationthatthechameleonscreening Thehalomassfunction,thenumberdensityofhaloswithagiven mechanismdoesnotworkveryefficientlyinoursimulations. mass,isausefultoolforinvestigatingtheefficiencyofamodel informinghalosofdifferentmasses.Tolocatehalosinthesimu- For γ gravity, the number of halos increases significantly, lationoutputs,weusedtheAHF(AmigaHaloFinder(Knollmann especiallyatthehighmassend,byupto40−100%forcluster- &Knebe2009)).Wedeterminedthehalomassfunctionbybin- sized halos. For Hu-Sawicki models, when the value of the f R ninghalosinlogarithmicmassintervals. field becomes comparable to the cosmological potential wells, In Fig. 5 we show the fractional difference with respect to the chameleon effect starts to operate. This can be seen in the ΛCDM of halo mass function computed from our simulations. massfunction,wheredeviationsfromΛCDMapproachzeroin Our measurement of the halo mass function itself is limited by thehigh-massendformodelswith|f |=10−5and|f |=10−6. R0 R0 statisticsandtoalesserextent,bytheresolutioninthehighand lowmassend.However,wecanreducetheimpactofthesetwo Thelargedeviationswefind∆n/nΛCDM ∼0.5−1inthehigh- effects by considering the relative difference between the halo massendareprobablyalreadyruledoutbypresentclustercounts massfunctionsmeasuredinmodifiedgravityandΛCDMsimu- (Cataneoetal.2015;Maketal.2012).Thisdoesnotruleoutthe lationswiththesameinitialcondition. model per se, but requires choices of parameters where |f | is R0 Figure5showsthatwefindanexcessofhalosintherange muchsmallertoday,leadingtonosignaturesinthebackground 1013 (cid:46) M/M (cid:46)1015probedbyoursimulationinγgravity,and evolution(i.e.,intheHubblefactor)oftheUniverse. (cid:12) Articlenumber,page6of10 M.VargasdosSantosetal.:N-bodySimulationsofγGravity Fig.6.FractionaldifferenceinthehalodensityprofileswithrespecttoΛCDMforfourdifferentmassbins.Forcomparisonwealsoshowthe resultsfromsimulationswiththeHu-Sawicky f(R)modelwith|f |={10−4,10−5,10−6}(dashedlines)fromLlinaresetal.(2014). R0 7.3. Haloprofiles and by 5% in the |f | = 10−6 case for all the three halo mass R0 ranges analyzed. Only for the |f | = 10−5 parameters a mixed R0 Wealsostudiedhowmodificationsofgravitychangethedensity behaviorisseen,thatis,forthelowertwohalomassrangesthe andvelocityprofilesofdarkmatterhalos.Wefocusedonhalos boostis∼ 20%,butfortheheaviesthalosthereisnodeviation in the mass ranges 0.5−1×1013Mpc/h, 1−5×1013Mpc/h, fromΛCDMintheinnerpartsand∼ 15%highervelocitiesare 0.5−1×1014Mpc/h,and1−5×1014Mpc/h.Thedensitypro- foundintheouterparts.Forγgravity,thedifferencebetweenthe fileswerecalculatedbybinningdarkmatterparticlesinannular models we simulated is more expressive for less massive halos bins for each halo. Our calculated density profiles are averages (5×1012 < M/M < 1013 and 1013 < M/M < 5×1013), for (cid:12) (cid:12) of all density profiles of the proper size, ranging from 10% of most massive halos (1014 < M/M < 5×1014) only the case thevirializationradius,r =0.1Rvir,totentimesthevirialization α=1.5isdistinguishablefromthe(cid:12)others. radius,r=10R .Thisrangewaschosentoproperlyincludeall vir behaviors of the fifth force on the dark matter halos while also avoiding the inner regions of the halos, where the resolution of 7.4. Fifthforceandscreening oursimulationsisnotsufficient. Figure 6 shows the fractional difference with respect to Generalrelativityisverywelltestedonverysmallscales,espe- ΛCDM in the density profiles. We first note that the inner re- cially inside the solar system. To ensure that Fφ is not relevant atthesescales,weneedascreeningmechanismtosuppressthe gions (R < R ) of halos for γ gravity are significantly denser vir than in ΛCDM. This difference is compensated for in outer re- fifthforceatsmall-scale,veryhighcurvatureregimes. gions (R > R ). Moreover, the profiles between the different WhenγgravitywasproposedinO’Dwyeretal.(2013),the vir modelparametersdonotdifferappreciablyfromeachother,and authors explored the compatibility between the model and the thispatternrepeatsforallranges.However,thedensityprofiles solar system experiments using the chameleon mechanism for for Hu-Sawicki models in general show stronger clustering in screening,asanyother f(R)(Braxetal.2008). thelow-densityregionsintheoutskirtsofhalosthanintheinner Asweshowedintheresultsabove,thescreeningmechanism regions. doesnotseemtobeworkingveryefficientlyforthemodelssim- Inthevelocitiesprofiles,showninFig.7,weexpectthefifth ulatedhere.Tocheckhowmuchthefifthforceisscreenedinour forcetoincreasethevelocitydispersion.Thiseffectisverysim- simulations, we compared the magnitude of the fifth and New- ilar for both models, except for most massive halos, for which tonian force on the particles in our simulation box, as in Davis theHu-Sawickimodelsaremorescreened,causingasubstantial etal.(2012). decreaseincomparisontoγgravity.FortheHu-Sawickimodel Figure8showsascatterplotforthiscomparisonatredshift the velocities are boosted by ∼ 20% in the |f | = 10−4 case z = 0.Thedispersionforsmall F < 0.01isexpected(numer- R0 N Articlenumber,page7of10 A&Aproofs:manuscriptno.main Fig.7.FractionaldifferenceinthevelocitydispersionprofileswithrespecttoΛCDMforfourdifferentmassbins.Forcomparisonwealsoshow theresultsfromsimulationswiththeHu-Sawicky f(R)modelwith|f |={10−4,10−5,10−6}(dashedlines)fromLlinaresetal.(2014). R0 ical scatter), here the forces are tiny, so the scatter here has a models considered here, Earth and Sun are almost completely very weak effect on the simulation. The important part in this screenedandpasslocalgravityconstraintsassumingthegalaxy figure is the behavior for large F . If we have screening, then isscreened. N weshouldseeF (cid:28) F /3inhigh-densityregions(whichcorre- However,fortheparametervaluesconsideredinthispaper, φ N spondstohighvaluesof F ).Theresultwefindshowsthatthe the screening condition gives that the galaxy is not screened N forceratioisroughlyonethird-whichisthelinearprediction- which again implies that the Earth and the sun is not screened everywhere,meaningthatthereisverylittlescreeningpresentin either3.Theonlycaveattothisisthatscreeningmighthavesur- oursimulations. vived from earlier times. At early redshift |f | (cid:28) 1 and almost R Tounderstandthisbetter,werevisitthescreeningconditions. allobjectsarescreened.Theninthecosmologicalevolution|f | Considering a spherical symmetric body with constant density R very quickly evolves to high values |f | = O(0.01−0.1). The ρ embeddedinthebackgroundofconstantdensityρ ,thesolu- R c b fieldinsideourgalaxymighthavebeentrapped,ensuringscreen- tionstothefieldequation(seee.g.Braxetal.(2012))meanthat ing.Thisisnotexpected,buttorigorouslyruleout(orconfirm) thefifthforceisgivenby this possibility, we would need simulations beyond the quasi- Fφ ∼ 13∆RRGrM2 e−mbr, sbteaytiocnadptphreoxscimopaetioofntthhiastpwaepears.sumedinoursimulations.Thisis where ∆R |f − f | = Rc Rb , R 2Φ 8. Summaryandconclusions N is the so-called screening factor (also called the thin-shell) and Wehaveinvestigatedthenonlinearevolutioninγgravity,a f(R) fRb (fRc)isthevalueof fR inthebackground(insidethebody), theoryofgravitythatisaviablealternativetoΛCDM.Themod- where ρ = ρb (ρc). If ∆RR (cid:28) 1, the fifth force is screened and els we investigated use a screening mechanism to suppress the we recover General Relativity. If ∆R (cid:38) 1, we instead find that deviations from General Relativity at small (solar system) and R F ∝ 1F andgravityissignificantlymodified. large cosmological scales. Specifically, this is the chameleon φ W3henNthefieldislocatedintheminimumofitseffectivepo- screeningmechanism.Asaresultofthisscreeningmechanism, tential in the background2 , we can calculate the screening fac- the strongest signatures in these models are expected to occur tor analytically. Assuming this, we find that for the γ gravity atthenonlinearregimeofstructureformation.Therefore,toun- 2 Backgroundherereferstothesurroundingsforthebodyinquestion, 3 Thisisonlytruefortheparametersconsideredinthispaper.If|f |(cid:46) R0 notnecessarilythecosmologicalbackground. 10−5,thenthegalaxy,Earth,andourSunarescreened. Articlenumber,page8of10 M.VargasdosSantosetal.:N-bodySimulationsofγGravity tweenthedifferentγgravitymodelsislargerforlow-masshalos, 5×1012 < M/M <1013,andwecandistinguishthemfromeach (cid:12) other;thedifferencesreach∼30%closetotheboundaryregion, R ∼ R . This difference decreases with the mass of the halos, vir andforthemostmassivehalos1014 < M/M <5×1014 .Only (cid:12) thecaseα=1.5islowerthan∼10%fromtheothers,whichare practicallyidentical. The chameleon mechanism - the screening mechanism that makes f(R)gravityviable-isnotveryeffectivefortheparame- terchoicesconsideredinthispaper.Thisexplainsthelargede- viations from ΛCDM in the observables we considered. Espe- ciallytheclustercountsignaturesof40−100%inthehigh-mass endM (cid:38)1014.5M disagreewithcurrentobservations.Thisdoes (cid:12) notruleoutthemodelperse,butrequirechoicesofparameters where|f |ismuchsmallertoday,whichimpliesthatthemodel R0 has no observable signatures in the background evolution (i.e., Fig.8.ComparisonbetweenthemagnitudesofNewtonianforce(F ) intheHubblefactor)oftheUniverse. N andfifthforce(F ),thedashedlineindicatesF /F =1/3.Theforces φ N φ areinunitsofH /c2. 0 Acknowledgements MVSthankstheBrazilianresearchagencyCAPESandtheUni- veiltheimprintsofsuchtheoriesatastrophysicalscales,weran versityofOsloforsupportandMaxB.GrönkeandAmirHam- severalcosmologicalN-bodysimulations.Wecomparedmodels with ΛCDM and the Hu-Sawicki model and showed that sev- mami for useful discussions. HAW is supported by BIPAC and the Oxford Martin School. DFM would like to thank the Re- eralastrophysicalobservables(halomassfunction,densitypro- search Council of Norway for funding. The simulations were files, and power spectra) show the signatures of the model sig- performed on the NOTUR cluster HEXAGON, which is the nificantlystrongly. computingfacilityattheUniversityofBergen. For the matter power spectrum we found a small deviation, lower than 10%, on large scales (k (cid:46) 10−1 h/Mpc), which is consistentwiththepredictionsoflinearperturbationtheory.For References small scales (k (cid:38) 10−1 h/Mpc), on the other hand, we found a strongincreaseinpower,thelargestdeviation(α=1.18)reaches Ade,P.A.R.etal.2014,Astron.Astrophys.,571,A22 Amendola, L., Gannouji, R., Polarski, D., & Tsujikawa, S. 2007, Phys. Rev., ∼ 60%,andintheothercases(α = 1.05andα = 1.5)itreaches D75,083504 ∼ 50%.ThisisquitedifferentfromtheresultsfoundintheHu- Appleby,S.A.&Battye,R.A.2007,Phys.Lett.,B654,7 Sawicky model, where the screening of the fifth force is much Appleby,S.A.,Battye,R.A.,&Starobinsky,A.A.2010,JCAP,1006,005 Berti,E.,Barausse,E.,Cardoso,V.,etal.2015,ArXive-prints more active in suppressing the enhancement of power on small Boehmer,C.G.,Burnett,J.,Mota,D.F.,&Shaw,D.J.2010,JHEP,07,053 scales, the largest deviation is ∼ 30%, much lower than for γ Bose,S.,Hellwing,W.A.,&Li,B.2015,J.CosmologyAstropart.Phys.,2,34 gravity. Brax,P.,Davis,A.-C.,Li,B.,&Winther,H.A.2012,Phys.Rev.D,86,044015 Brax,P.,Davis,A.-C.,Li,B.,Winther,H.A.,&Zhao,G.-B.2013,J.Cosmology Forthehalomassfunctionwefoundthatallofourrunsgave Astropart.Phys.,4,29 verysimilarresultinthemassrange11.5<log(M/M )<14.7, (cid:12) Brax,P.,vandeBruck,C.,Davis,A.-C.,&Shaw,D.J.2008,Phys.Rev.D,78, and these results are very similar to what is found in the Hu- 104021 Sawickymodelfor|f |≥10−5;a10−60%increaseinthehalo Cai,Y.-C.,Padilla,N.,&Li,B.2014,ArXive-prints R0 abundance. However, for γ gravity we found an enhancement Cai,Y.-C.,Padilla,N.,&Li,B.2015,MNRAS,451,1036 Capozziello,S.,Cardone,V.F.,Carloni,S.,&Troisi,A.2003,Int.J.Mod.Phys., of 40−100% in the abundance for log(M/M ) > 14.5, which (cid:12) D12,1969 canbecomparedtotheHu-Sawickimodel,wherethehalomass Carroll,S.M.,Duvvuri,V.,Trodden,M.,&Turner,M.S.2004,Phys.Rev.,D70, functionapproachesthatofΛCDM. 043528 For halo density profiles we found that the different runs Cataneo,M.,Rapetti,D.,Schmidt,F.,etal.2015,Phys.Rev.D,92,044009 Cognola,G.,Elizalde,E.,Nojiri,S.,etal.2008,Phys.Rev.,D77,046009 with γ gravity were not distinguished well in any mass range; Colombi,S.&Novikov,D.2011,POWMES:MeasuringthePowerSpectrumin inallsituationswesawthattheinnerregions(R<R )ofhalos anN-bodySimulation,AstrophysicsSourceCodeLibrary vir for γ gravity are significantly thicker than the halos of ΛCDM, Davis,A.-C.,Li,B.,Mota,D.F.,&Winther,H.A.2012,Astrophys.J.,748,61 around10%forlessmassivehalos,5×1012 < M/M <1013,in- delaCruz-Dombriz,A.,Dobado,A.,&Maroto,A.L.2008,Phys.Rev.,D77, (cid:12) 123515 creasing with the halo mass until they reached ∼ 30% for the Frolov,A.V.2008,Phys.Rev.Lett.,101,061103 most massive halos, 1014 < M/M < 5 × 1014. This differ- Gronke,M.,Llinares,C.,Mota,D.F.,&Winther,H.A.2015,Mon.Not.Roy. (cid:12) enceiscompensatedforintheouterregions(R > R ).Forthe Astron.Soc.,449,2837 vir Hu-Sawicki model the same effect acts for less massive halos Gronke,M.,Llinares,C.,Mota,D.F.,&Winther,H.A.2015,MNRAS,449, 2837 5×1012 < M/M(cid:12) <1013,butitchangesformoremassivehalos, Hammami,A.,Llinares,C.,Mota,D.F.,&Winther,H.A.2015,Mon.Not.Roy. 5×1013 < M/M(cid:12) <1014and1014 < M/M(cid:12) <5×1014whenthe Astron.Soc.,449,3635 screeningmechanismsuppressesthefifthforce.Thisthenleads He,J.-h.,Li,B.,&Hawken,A.J.2015,ArXive-prints He,J.-h.,Li,B.,&Jing,Y.P.2013,Phys.Rev.D,88,103507 to a suppression in the accretion of mass concentration in the Hellwing,W.A.,Barreira,A.,Frenk,C.S.,Li,B.,&Cole,S.2014,Physical innerregions(comparedtothenon-screenedcase). ReviewLetters,112,221102 Finally, we computed the velocity profiles, and as we ex- Hellwing,W.A.,Cautun,M.,Knebe,A.,Juszkiewicz,R.,&Knollmann,S.2013, pected, the velocities are significantly enhanced in γ gravity JCAP,1310,012 comparedtoΛCDM.Theothercasesarecomparableinallmass Hu,W.&Sawicki,I.2007,Phys.Rev.,D76,064004 Jain,B.&Khoury,J.2010,AnnalsPhys.,325,1479 ranges, but with some important signatures. The difference be- Knollmann,S.R.&Knebe,A.2009,Astrophys.J.Suppl.,182,608 Articlenumber,page9of10 A&Aproofs:manuscriptno.main Li,B.,Hellwing,W.A.,Koyama,K.,etal.2013,MNRAS,428,743 Li,B.,Hellwing,W.A.,Koyama,K.,etal.2013,Mon.Not.Roy.Astron.Soc., 428,743 Li,B.,Mota,D.F.,&Barrow,J.D.2011,Astrophys.J.,728,109 Li,B.,Zhao,G.-B.,Teyssier,R.,&Koyama,K.2012,J.CosmologyAstropart. Phys.,1,51 Li,B.,Zhao,G.-B.,Teyssier,R.,&Koyama,K.2012,JCAP,1201,051 Linder,E.V.2009,Phys.Rev.,D80,123528 Llinares,C.&Mota,D.2013,Phys.Rev.Lett.,110,161101 Llinares,C.&Mota,D.F.2014,Phys.Rev.,D89,084023 Llinares,C.,Mota,D.F.,&Winther,H.A.2014,Astron.Astrophys.,562,A78 Lombriser,L.,Koyama,K.,&Li,B.2014,JCAP,1403,021 Lombriser,L.,Schmidt,F.,Baldauf,T.,etal.2012,Phys.Rev.,D85,102001 Mak,D.S.Y.,Pierpaoli,E.,Schmidt,F.,&Macellari,N.2012,Phys.Rev.D, 85,123513 Mota,D.F.,Shaw,D.J.,&Silk,J.2008,Astrophys.J.,675,29 Noller,J.,vonBraun-Bates,F.,&Ferreira,P.G.2014,Phys.Rev.D,89,023521 O’Dwyer,M.,Joras,S.E.,&Waga,I.2013,Phys.Rev.,D88,063520 Oyaizu,H.,Lima,M.,&Hu,W.2008,Phys.Rev.,D78,123524 Perlmutter,S.etal.1999,Astrophys.J.,517,565 Pogosian, L. & Silvestri, A. 2008, Phys. Rev., D77, 023503, [Erratum: Phys. Rev.D81,049901(2010)] Puchwein,E.,Baldi,M.,&Springel,V.2013,MNRAS,436,348 Pujol,A.etal.2014,Mon.Not.Roy.Astron.Soc.,438,3205 Riess,A.G.etal.1998,Astron.J.,116,1009 Schmidt,F.2010,Phys.Rev.D,81,103002 Schmidt,F.2010,Phys.Rev.,D81,103002 Schmidt,F.,Lima,M.V.,Oyaizu,H.,&Hu,W.2009,Phys.Rev.,D79,083518 Shi,D.,Li,B.,Han,J.,Gao,L.,&Hellwing,W.A.2015,Mon.Not.Roy.Astron. Soc.,452,3179 Silveira,V.&Waga,I.1994,Phys.Rev.,D50,4890 Song,Y.-S.,Taruya,A.,Linder,E.,etal.2015,Phys.Rev.D,92,043522 Spergel,D.,Gehrels,N.,Baltay,C.,etal.2015,ArXive-prints Starobinsky,A.A.1980,Phys.Lett.,B91,99 Starobinsky,A.A.2007,JETPLett.,86,157 Terukina,A.&Yamamoto,K.2012,Phys.Rev.,D86,103503 Tessore,N.,Winther,H.A.,Metcalf,R.B.,Ferreira,P.G.,&Giocoli,C.2015, J.CosmologyAstropart.Phys.,10,36 Teyssier,R.2002,Astron.Astrophys.,385,337 Thongkool,I.,Sami,M.,Gannouji,R.,&Jhingan,S.2009,Phys.Rev.,D80, 043523 Wilcox,H.,Bacon,D.,Nichol,R.C.,etal.2015,MNRAS,452,1171 Winther,H.A.,Mota,D.F.,&Li,B.2012,Astrophys.J.,756,166 Winther,H.A.,Schmidt,F.,Barreira,A.,etal.2015,ArXive-print Zhang,P.2006,Phys.Rev.,D73,123504 Zhao,G.-B.,Li,B.,&Koyama,K.2011,Phys.Rev.D,83,044007 Zivick,P.,Sutter,P.M.,Wandelt,B.D.,Li,B.,&Lam,T.Y.2015,MNRAS, 451,4215 Zumalacarregui, M., Koivisto, T. S., & Mota, D. F. 2013, Phys. Rev., D87, 083010 Articlenumber,page10of10