Logarithmic corrections to Newtonian gravity and Large Scale Structure
Salvatore Capozziello, Mir Faizal, Mir Hameeda, Behnam Pourhassan, Vincenzo Salzano
aa r X i v : . [ g r- q c ] F e b Eur. Phys. J. C manuscript No. (will be inserted by the editor)
Gravitational Nonlocality from Large Scale Structure
Salvatore Capozziello a,1,2,3,4 , Mir Faizal b,5,6,7 , Mir Hameeda c,8,9 , BehnamPourhassan d,10 , Vincenzo Salzano e,11 Department of Physics “E. Pancini”, University of Naples “Federico II”, Naples, Italy INFN Sez. di Napoli, Compl. Univ. di Monte S. Angelo, Edificio G, Via Cinthia, I-80126, Naples, Italy Scuola Superiore Meridionale, Largo San Marcellino 10, I-80138, Naples, Italy Laboratory for Theoretical Cosmology, Tomsk State University of Control Systems and Radioelectronics (TUSUR), 634050Tomsk, Russia Irving K. Barber School of Arts and Sciences, University of British Columbia - Okanagan, 3333 University Way, Kelowna,British Columbia V1V 1V7, Canada Department of Physics and Astronomy, University of Lethbridge, Lethbridge, Alberta, T1K 3M4, Canada Canadian Quantum Research Center 204-3002 32 Ave Vernon, BC V1T 2L7 Canada Department of Physics, Government Degree College Tangmarg, Kashmir 193402, India Visiting Associate, IUCCA, Pune, 41100, India School of Physics, Damghan University, Damghan, 3671641167, Iran Institute of Physics, University of Szczecin, Wielkopolska 15, 70-451 Szczecin, PolandReceived: date / Accepted: date
Abstract
Gravitational nonlocality effects at large scalecan be investigated using the cosmological structureformation. In particular, it is possible to show if andhow gravitational nonlocality modifies the clusteringproperties of galaxies and of clusters of galaxies. Thethermodynamics of such systems can be used to ob-tain important information about the effects of suchlager scale nonlocality on clustering. We will comparethe effects of such larger scale nonlocality with obser-vational data. It will be demonstrated that the obser-vations seem to point to a characteristic scale of nonlo-cality of the gravitational interactions at galactic scales.However, at larger scales such statistical inferences aremuch weaker.
The clustering of galaxies is the main mechanism toaddress the large scale structure of the Universe. Suchclustering mechanism can be studied confronting nu-merical simulations [1, 2] with observations [3, 4], usingthe local matter distributions of galaxies organized ingroups, filaments and clusters. However, to analyze nu-merically the large scale structure formation, it is possi-ble to approximate galaxies as point-like particles, and a e-mail: [email protected] b e-mail: [email protected] c e-mail: [email protected] d e-mail: [email protected] e e-mail: [email protected] then study the clustering adopting the standard formal-ism of statistical mechanics. This approximation is validbecause the distance between galaxies is much largerthan their proper size. In fact, it has been demonstratedthat the clustering mechanism can be studied consider-ing a quasi-equilibrium description because the macro-scopic quantities change very slowly as compared withthe local relaxation time scales [5–8]. Pressure, aver-age density and average temperature of clusters are themacroscopic quantities to be taken into account. Herethe temperature is the effective one obtained by the ki-netic theory of gases, for a gas of galaxies (each galaxybeing approximated as a point particle). So, by thisquasi-equilibrium description, the clustering of galax-ies can be dealt as a thermodynamic system [9, 10].Specifically, one has to adopt a gravitational partitionfunction [11, 10] where point-like galaxies gravitation-ally interact. Such a partition function can be also usedto study phase transitions for systems of interactinggalaxies [12–14].It has to be pointed out that such a gravitationalpartition function diverges because the extended struc-ture of galaxies has been neglected in this approxima-tion. However, these divergences can be removed byadding a softening parameter, which accounts for theextended structure of galaxies as reported in [15]. Ina recent study it has also been proved that the non-local gravity is divergence free [16]. Thermodynamicsfor systems of galaxies is obtained from the partitionfunction (regularized with a softening parameter), andthis can in turn be used to study the clustering proper- ties [17]. It has been observed that any modification ofgravitational partition function by a softening parame-ter, in turn, modifies the thermodynamic quantities ofa given system. Vice versa, because these thermody-namic quantities are related to the clustering parame-ter, it also changes the clustering parameter itself. Infact, it can be demonstrated that the modification tothe skewness and kurtosis of clustering systems occurdue to such a softening parameter [18]. This formalismhas also been used to demonstrated that galaxy clus-ters are surrounded by halos, and this has been doneadopting a wide range of samples [19, 20]. Furthermorethe isothermal compressibility gives also information onthe clustering of galaxies [21].Specifically, the clustering of galaxies occurs due tothe gravitational force between them. However, it isnow well-known by a wide range of cosmological ob-servations (e.g. Type Ia Supernovae, Baryon AcousticOscillations and Cosmic Microwave Background) thatthe Universe is undergoing an accelerating expansion[22–24]. This cosmic expansion is pulling the galaxiesaway from each other. So, the cosmic expansion seemsto oppose the gravitational force, and thus it is impor-tant to include its effect in the gravitational partitionfunction. It has been argued that the cosmological con-stant related to the accelerated expansion can produceimportant local effects on the scale of galaxy clusters[25–29]. These effects can be investigated using modelof interacting dark matter and dark energy [25, 26]. Itis also possible to constrain the dark energy from theclustering of galaxies [30]. The cosmic expansion canalso change the dynamics at the scale of galaxy clus-ters [31]. Such a modification of cluster dynamics hasalso been investigated using different groups of galax-ies [32]. Finally, one can state that the cosmologicalconstant has consequences for the formation of galaxyclusters [33]. In summary, it is important to study theeffects of cosmic expansion at the scale of clusters. Froma more technical point of view, it has to be noted thatthe Fisher matrix formalism has been used to analyzethe effects of cosmic expansion on clusters [34].Thus, it is possible to obtain information about thecosmological constant, and, in general, about dark en-ergy dynamics, from the local clustering of galaxies. So,it is important to incorporate the effects of cosmic ex-pansion in the statistical mechanical description of theclustering of galaxies. In fact, an accurate measurementof the statistics of galaxies can be used to constraint thevalue of the cosmological constant [35], and the distri-bution function of galaxies can be used to constraintthe amount of the dark energy [36].As such a statistical distribution function can be ob-tained from the gravitational partition function [9, 10], it is important to incorporate the effects of cosmic ex-pansion in the partition function. So, the modificationof the gravitational partition function by the cosmolog-ical constant term has also to be considered, and thismodified gravitational partition function can be usedto study the clustering process in an expanding Uni-verse [37]. This can be technically done by first analyz-ing the Helmholtz free energy and entropy for a systemof galaxies. Then the clustering parameter for this sys-tem can be obtained, and it has been observed that thisclustering parameter depends on the value of the cosmo-logical constant. The modification of the gravitationalpartition function with a time dependent cosmologicalconstant can be also used to study the effects of dy-namical dark energy on clustering of galaxies [38]. Ithas been observed that the correlation function for thissystem are consistent with observations.However, it is possible to obtain cosmic expansionwithout considering dark energy related to some newfundamental particle. In fact, large scale nonlocality re-lated to some effective gravitational potential seems arealistic effect of gravitational field at infrared scales asrecently reported in several papers, see for example [39–43]. It has been argued that such nonlocalities can occurdue to some quantum gravitational effects [39, 40, 44].They can be relevant also in gravitational waves modesand polarizations [45, 46]. It has been demonstratedthat string field theory produces nonlocal actions forthe component fields, and these nonlocal actions canalso be used to construct nonlocal cosmological models[47, 48]. It is also possible to construct models of non-local cosmology from loop quantum gravity [49, 50].This can be done adopting a sort of condensation ofstates in loop quantum gravity. In summary, it is possi-ble to obtain nonlocal cosmologies using quantum grav-itational effects at ultraviolet scales that propagate upto infrared scales.In other words, as nonlocal cosmology can explainthe cosmic expansion without dark energy [39–43], suchnonlocal corrections can also be used to incorporate theeffects of cosmic expansion in the gravitational parti-tion function. This one can then be used to analyze theeffects of nonlocality on clustering of galaxies. It maybe noted that even though general relativity is stronglyconstrained at the Solar System scales by several exper-iments and observations [51–53], nonlocal modificationto general relativity can occur at large scales startingfrom galactic one [54]. It has also been proposed thatsuch nonlocality could be constrained from the cluster-ing of galaxies [55]. Thus, it is important to analyze theeffect of nonlocal modification of gravity on the gravi-tational partition function and point out how it affectsthe clustering process. Such a modification can be an- alyzed using nonlocal corrections to the gravitationalpotential [54, 56–59]. In general, the nonlocal gravita-tional partition function can be used to analyze theeffects of nonlocality on clustering of galaxies.This paper is devoted to this issue.In Sec. 2, we discuss the clustering mechanism. Af-ter defining the nonlocal gravitational partition func-tion, we consider the related thermodynamics and spa-tial correlation function. Observational data analysisis described in detail in Sec. 3. Specifically, we takeinto account data from galaxies and clusters and studythe effects of nonlocality on the clustering mechanism.Results are discussed in Sec. 4, while conclusions aredrawn in Sec. 5. In general, a nonlocal modification of the Hilbert-Einsteinaction can have the following form [39, 40] S = 12 κ Z d x √− g (cid:2) R (cid:0) f ( (cid:3) − R ) (cid:1)(cid:3) , (1)where R is the Ricci scalar and f ( (cid:3) − R ) is an arbitraryfunction, called distortion function , of the nonlocal term (cid:3) − R , which is explicitly given by the retarder Green’sfunction G [ f ]( x ) = ( (cid:3) − f )( x ) = Z d x ′ p − g ( x ′ ) f ( x ′ ) G ( x, x ′ ) . (2)Setting f ( (cid:3) − R ) = 0 , the above action is equivalent tothe Einstein-Hilbert one. The nonlocality is introducedby the inverse of the d’Alembert operator. In order to"localize" the action, an auxiliary scalar field can beintroduced so that (cid:3) − R = φ and then formally (cid:3) φ = R . As a consequence, characteristic lengths related tothe nonlocal terms naturally come out [54].The weak-field limit of such a dynamics gives riseto corrections to the Newtonian potential which are in-teresting at large scales. These corrections can be poly-nomial or logarithmic and, in general, introduce char-acteristic scales. For a detailed discussion on this pointsee [54].In the framework of this theory, let us now model theclustering of galaxies using the gravitational partitionfunction which we are going to define. We will explicitlyadopt the nonlocal gravity, as discussed in [57–59], toobtain the corrections to such a gravitational partitionfunction. After analyzing the effects of nonlocality, thenonlocal gravitational partition function will be usedto analyze the clustering of galaxies. This can be doneby considering the thermodynamics of this system, and then relating it to the clustering parameter. Then wewill calculate the spatial correlation function.2.1 Non-local Gravitational Partition FunctionLet us take into account a system with a large numberof galaxies. This system can be analyzed using a quasi-equilibrium description, where the change in macro-scopic quantities is slower than the local relaxation timescales [5–8]. It is possible analyze such a system adopt-ing an ensemble of cells, with same volume V , or radius R and average density ρ . As the number of galaxiesand their total energy can vary between these cells, thesystem can be analyzed in the framework of the grandcanonical ensemble. Thus, it is possible to define a grav-itational partition function. In this picture, a system of N galaxies of mass m interacting with a nonlocal grav-itational potential can be written as [11, 10] Z ( T, V ) = 1 λ N N ! Z d N pd N r × exp (cid:18) − (cid:20) N X i =1 p i M + Φ nl ( r ) (cid:21) T − (cid:19) , (3)where p i are the momenta of different galaxies and Φ nl is the nonlocal gravitational potential energy. Here T isthe average temperature, which is obtained from the ki-netic theory of gases where each particle is representedby a galaxy. Now integrating on the momentum space,we can write this gravitational partition function as Z N ( T, V ) = 1 N ! (cid:18) πmTΛ (cid:19) N/ Q N ( T, V ) , (4)where Q N ( T, V ) , the configuration integral, can be ex-pressed as Q N ( T, V ) = Z .... Z Y ≤ i 2; 5 / − R /ǫ ) is the hypergeometric function. Now using R ∼ ρ − / ,we observe that Gm / R T = 3 Gm / ρ − / T = 3 G m ρ / T − / . So, using the scale invariance, ρ → λ − ρ , T → λ − T and r → λr , we can define x as x = 32 ( Gm ) ρT − = βρT − , (14)where β = 3( Gm ) / . From this expression for x ,we can obtain any general configuration integral. Us-ing the so called dilute approximation , we can assume ( Φ i,j ) nl /T very small and take only the first term of theexponential expansion, f ij = Gm T q r ij + ǫ − Gm T λ ln q r ij + ǫ + Gm T λ ln λ (15)Thus, calculating the general configuration integral, weobtain Q N ( T, V ) = V N (1 + αx ) N − (16)where α = 2 R λ ln λR + 3 λ R r ǫ R + 1 + ǫ R tan − Rǫ (17) − λǫ R ln q ǫ R + 1 + 1 ǫR − 12 ln (cid:18) ǫ R + 1 (cid:19) − ǫ R + 13 Now the gravitational partition function for this dilutegravitating system can be written as Z N ( T, V ) = 1 N ! (cid:0) πmTΛ (cid:1) N/ V N (cid:0) αx (cid:1) N − (18)This expression for the gravitational partition functioncan be used to analyze the effects of nonlocality on theclustering of galaxies.2.2 ThermodynamicsIt has to be noted that the gravitational partition func-tion has already been used to study the thermodynam-ics of systems of galaxies, see [9] and [10]. We can thusadopt the nonlocal gravitational partition function toanalyze the thermodynamics of this nonlocal system.The Helmholtz free energy F = − T ln Z N ( T, V ) is F = − T ln (cid:18) N ! (cid:0) πmTΛ (cid:1) N/ V N (cid:0) αx (cid:1) N − (cid:19) . (19)In the plots of Fig. 1, we analyze the dependence ofthe Helmholtz free energy on various parameters. Thedependence of Helmholtz free energy on x is representedin Fig. 1 (a). We can observe that the Helmholtz freeenergy is a decreasing function of x for large value ofthe nonlocality parameter λ and an increasing functionof x for smaller value of the nonlocality parameter λ .So, the value of nonlocality parameter changes the de-pendence of the Helmholtz free energy on x . In Fig. 1 ( a ) - xF ( b ) λ F ( c ) TF ( d ) NF Fig. 1 Helmholtz free energy dependence on: (a) x with: T = 1 , N = 50 , λ = 2 / (solid green), λ = 1 (dashed red), λ = ∞ (dotted blue); (b) λ with: x = 1 , N = 50 , T = 0 . (solid green), T = 0 . (dashed red), T = 0 . (dotted blue); (c) T with: x = 1 , N = 60 (solid green), N = 40 (dashed red), N = 20 (dotted blue); (d) N with: T = 1 , λ = 0 . ; x = 1 (solid green), x = 0 . (dashed red), x = 0 (dotted blue). We have assumed unit values for all the other parameters. (b), the Helmholtz free energy is plotted in terms of λ .We can see how there is a minimum for the Helmholtzfree energy, and this minimum seems to occur at a spe-cific value of λ . This minimum does not seem to changewith the temperatures of the system. This observationcan explain the change in the behavior of the Helmholtzfree energy with the change in the value of λ , as wasobserved in Fig. 1 (a). In Fig. 1 (c), we plot the depen-dence of the Helmholtz free energy on the temperatureof the system. Here the temperature corresponds to thetemperature of a gas of galaxies, with each galaxy act-ing as a particle analog. We observe that generally theHelmholtz free energy increases with the temperature,but the rate of this increase depends on the number ofgalaxies. As we have plotted the dependence on tem-perature by considering the system as a gas of galaxies,we also analyze the dependence of the Helmholtz freeenergy on the number of galaxies in Fig. 1 (d). We ob-serve that the Helmholtz free energy first decreases for arelatively small number of galaxies, becoming negative,and then increases as the number of galaxies becomesmore than a certain value. The entropy S can now be calculated from this Helmholtzfree energy, as S = − (cid:18) ∂F∂T (cid:19) N,V = N ln( ρ − T / ) + ( N − 1) ln (cid:0) αx (cid:1) − N αx αx + 52 N + 32 N ln (cid:0) πmΛ (cid:1) . (20)Now, for large N , using N − ≈ N , we obtain SN = ln( ρ − T / ) + ln (cid:0) αx (cid:1) − B l + S N , (21)where S = N + N ln (cid:0) πmΛ (cid:1) . where B l = αx αx . (22)This is the general clustering parameter for a system ofgalaxies interacting in presence of nonlocal gravity. Theinternal energy U = F + T S of a system of galaxies cannow be written as U = 32 N T (cid:0) − B l (cid:1) . (23)We have plotted the internal energy of a system ofgalaxies in Fig. 2. Its behavior is analyzed by investi-gating its dependence on various parameters. In Fig. 2 ( a ) xU ( b ) λ U ( c ) - - TU ( d ) NU Fig. 2 Internal energy dependence on: (a) x with: T = 1 , N = 50 , λ = 2 / (solid green), λ = 1 (dashed red), λ = 0 (dottedblue); (b) λ with: x = 1 , N = 50 , T = 0 . (solid green), T = 0 . (dashed red), T = 0 . (dotted blue); (c) T with: x = 1 , N = 60 (solid green), N = 40 (dashed red), N = 20 (dotted blue); (d) N with: T = 1 , λ = 0 . ; x = 1 (solid green), x = 0 . (dashed red), x = 0 (dotted blue). We have assumed unit values for all the other parameters. (a), we investigate how the internal energy of the sys-tem varies with x . We observe that it is an increasingfunction of x for smaller values of λ , and a decreas-ing function of x for larger values of λ . In Fig. 2 (b),we analyze the dependence of the internal energy onthe nonlocality parameter λ . We observe that there isa minimum for internal energy, and again this mini-mum does not seem to depend on the temperature ofthe system. We also investigate the dependence of theinternal energy on the temperate of the system in Fig.2 (c). It is observed that there is a discontinuity in thedependence of internal energy on temperature, whichoccurs independently of the number of galaxies. In Fig.2 (d), the dependence of internal energy on the numberof galaxies is plotted. It is observed that internal energyincreases with the increase in the number of galaxies.However, the rate of increase is larger for larger valuesof x .Similarly, we can write the pressure P and chemicalpotential µ for galaxies interacting through a nonlocal gravitational potential as P = − (cid:18) ∂F∂V (cid:19) N,T = N TV (cid:0) − B l (cid:1) ,µ = (cid:18) ∂F∂N (cid:19) V,T = T ln( ρT − / ) (24) − T ln (cid:0) αx (cid:1) − T 32 ln (cid:0) πmΛ (cid:1) − T B l . The probability of finding N galaxies can be written as F ( N ) = P i e NµT e − UnT Z G ( T, V, z ) = e NµT Z N ( V, T ) Z G ( T, V, z ) , (25)where Z G is the grand-partition function defined by Z G ( T, V, z ) = ∞ X N =0 z N Z N ( V, T ) . (26)and z is the activity. Thus, for a system of gravitation-ally interacting galaxies in nonlocal gravity , we canwrite exp (cid:20) N µT (cid:21) = (cid:18) ¯ NV T − / (cid:19) N (cid:18) B l (1 − B l ) (cid:19) − N · exp [ − N B l ] (cid:18) πmΛ (cid:19) − N/ . (27) Now the grand-partition function can be written as ln Z G = P VT = ¯ N (1 − B l ) , (28)and the distribution function can be expressed as F ( N ) = ¯ N (1 − B l ) N ! (cid:0) N B l + ¯ N (1 − B l ) (cid:1) N − · exp (cid:2) − N B l − ¯ N (1 − B l ) (cid:3) . (29)In Fig. 3, we analyze the behavior of the distributionfunction for galaxies, and its dependence on differentparameters used here. In general, the distribution func-tion has the universal feature of having a maximum,while its functional shape changes, and depends on theconsidered parameters.2.3 Spatial correlation functionIt is well know that the interaction between differentgalaxies will cause correlation in their positions. The in-tegral of the correlation function over a certain volumecan be expressed in terms of the mean square fluctua-tion of the total number of galaxies in a given volume.Thus we can write Z ξdV = < ( ∆N ) > ¯ N − (30)This can be further represented in terms of thermody-namic quantities as Z ξdV = − N TV (cid:18) ∂V∂P (cid:19) T − (31)Taking the volume derivative and using the equation ofstate for pressure P , we get ξ = − N TV ∂∂V (cid:18) ∂V∂P (cid:19) T + 2 N TV (cid:18) ∂V∂P (cid:19) T (32)The pressure can be written in terms of the clusteringparameter as P = N TV (cid:0) − B l (cid:1) , (33)We can write the change of volume with pressure atconstant temperature (cid:18) ∂V∂P (cid:19) T = − V N T (1 − B l ) (34)The expression for ξ is ξ = − B l V (1 − B l ) (35) where ∂B l /∂V = − ( x/V )( dB l /dx ) = B l (1 − B l ) /V and x = βρT − = βN T − /V . Now we can write Z ξdV = B l (2 − B l )(1 − B l ) (36)The volume integral of correlation function is the quan-tity which can be compared with observational data inorder to detect possible effects of nonlocal gravity. We will test the nonlocal gravitational clustering withthe cluster catalog in [60], containing , clustersof galaxies in the redshift range . ≤ z < . fromthe Sloan Digital Sky Survey III (SDSS-III). We willtake advantage of this catalog because it has all theneeded ingredients to analyze clustering properties ofboth galaxies and clusters (although an alternative morerecent version is in [61]), using two different approaches.3.1 GalaxiesThe analysis of the clustering of galaxies will be con-ducted with a sort of “smart” version of the standardcounts-in-cell procedure. In general, one would first needto define cells/volumes with a given size, and then tocount the gravitational structures of interest within them.Here, our “smart” volumes will be the clusters identifiedin the SDSS-III survey, and we will count the galaxieswhich are within each of them. From this procedure, wewill be able to calculate the observed distribution func-tion F ( N ) and to compare it with the theoretical ex-pectation in nonlocal gravity as given by Eq. (29). The“smartness” in this procedure is that the cells/volumeswe are going to analyze are physically realized systems,structures which exist and have been identified, con-trarily to the standard way of proceeding, where theexistence of clusters and interacting systems of galaxiesis not assured neither verified. On the other hand, per-forming the count in this way, we are going to inevitablymiss voids, which instead indirectly retain some infor-mation about the clustering properties.All the needed data are provided by the cataloguein [60]: the volume of the cells or, equivalently, the ra-dius of the clusters, r in Mpc, is defined as the radiuswithin which the mean density of a cluster is timesthe critical density of the Universe at the same redshiftand N is the number of member galaxy candidateswithin r of each cluster. The only caveat to take inmind is that the radii of the cells we are considering,i.e. the r , are much smaller than the scales where ( a ) x F ( N ) ( b ) xF ( N ) ( c ) NF ( N ) ( d ) NF ( N ) Fig. 3 The distribution function dependence on: (a) x with: N = 50 , ¯ N =10, λ = 10 (solid green), λ = 1 (dashed red), λ = 0 . (dotted blue); (b) x with: ¯ N =10, λ = 1 , N = 20 (solid green), N = 10 (dashed red), N = 5 (dotted blue); (c) N with: x = 1 , ¯ N = 10 , λ = 10 (solid green), λ = 1 (dashed red), λ = 0 . (dotted blue); (d) N with: ¯ N = 10 , λ = 1 ; x = 10 (solid green), x = 5 (dashed red), x = 1 (dotted blue). We have assumed unit values for all the other parameters. the quasi-equilibrium clustering should be more effec-tive [62, 63]; we will discuss this point when presentingour results.We have divided the full catalogue in groups by red-shift, with bin widths of ∆z = 0 . , and by radius, withbin widths of ∆r = 0 . Mpc, and we have eventuallyselected only the groups with a sufficiently large num-ber of clusters to enable a strong and reliable statisti-cal analysis. The final groups with which we will workwill have radii in the three ranges, . < r < . , . < r < . and . < r < . Mpc, and red-shifts in the interval . ≤ z ≤ . . We need to re-mind here that the catalogue is complete only up to aredshift z ∼ . , for clusters with an estimated mass M > · M ⊙ ; at higher redshift, a bias towardsmaller clusters (lower masses) and possibly to lowercounts-in-cell is possible [60].In Table 1, we show the results derived from usingEq. (29) with the chosen data. The parameters involvedin the analysis of Eq. (29) have been rewritten as:- R = r , the radius of the cell (cluster), which en-ters Eq. (29) through Eq. (17). We have assumed itto be constant within the three bins we have iden-tified; in the following, it will act as a scaling factorfor some parameters; - N = N , the number of galaxies within each cell(cluster);- the clustering parameter b , defined as b = x x , (37)and which enters Eq. (29) through: B l = b α b ( α − , (38)with α given by Eq. (17);- the dimensionless softening parameter, ˜ ǫ = ǫ/R ;- the dimensionless nonlocal characteristic length, ˜ λ = λ/R .The only parameter which fully characterizes the non-local gravity model we are focusing on is thus ˜ λ . Allother parameters are standard in the sense that theywould appear also whether standard GR would be con-sidered. In Table 1 we also report the number of clustersin each redshift bin, N cl , although this is not a fittingparameter.When it comes to fit the data, we have defined andcompared results from two different statistical tools. Wehave performed a least-square minimization using the χ defined as χ ls = X i (cid:2) F theo ( N i ) − F obs ( N i ) (cid:3) , (39)where F obs ( N i ) is the real counts-in-cell distributionextracted from the data, and F theo ( N i ) is the theoret-ical counts-in-cell calculated from the nonlocal gravitymodel by using Eq. (29). The index i derives from thefact that in each radius bin, the clusters have a vari-able amount of galaxies within them, ranging from someminimal number N min to a maximal one N max . Thus,the index i selects the finite natural values of this range,where we evaluate the distribution function. We havealso defined a second χ as χ jk = X i (cid:2) F theo ( N i ) − F obs ( N i ) (cid:3) σ i , (40)where the σ i are the errors on F obs ( N i ) which we havederived from a jackknife-like procedure as exposed be-low.For each redshift and radius bin:1. we cut a variable fraction F (randomly selected inthe range [10 , ) from the total bin population;2. for each cut sample we derive the counts-in-cell dis-tribution function ∼ times;3. we thus obtain a distribution of F obs ( N i ) from the ∼ sets of N i ;4. from each distribution, which looks very close to astandard Gaussian, we derive the standard devia-tion and assume such value as the error σ i .The data points and the errors on the distribution func-tion F obs ( N i ) which we finally obtain are shown respec-tively as black dots and bars in Figs. (5) - (6) - (7). Thebest fitting F ( N ) distributions obtained from the min-imization of χ ls are shown in red; those ones derivedfrom χ jk are in green.The minimization of the defined χ is performed byusing a Monte Carlo Markov Chain (MCMC) approach,running chains with points and using the uninfor-mative and very general priors: ¯ N ≥ ; ≤ b ≤ ; ≤ ˜ ǫ ≤ ; ˜ λ > .3.2 ClustersConcerning the count-in-cells where the counted objectswithin volumes are the clusters of galaxies, we follow aprocedure which is very similar to what is described in[62]. In order to be as much as clear as possible, we aregoing to enumerate all the steps in the following list.The initial set up consists in the coming steps: (1a) we choose the size of the spherical cells/volumes weare going to analyze. We have selected four differentphysical lengths, R = 10 , , , and Mpc;(2a) we divide the clusters in three groups by redshift, . ≤ z < . , . ≤ z < . and . ≤ z < . . These values have been chosen in order tospan approximately the same comoving volume ineach group: assuming a baseline Planck cosmology[24], with H = 67 . km s − Mpc − and Ω m =0 . , the volume is ≈ . · Mpc . The numberof clusters from the SDSS catalog falling into eachgroup are, respectively, , and (thetotal sums up to , because we are taking onlyclusters with z ≤ . , for which the catalogue iscomplete);(3a) all clusters in the catalog are provided with red-shift z and celestial coordinates, the J . rightascension β and declination δ . These quantities areconverted into physical Cartesian coordinates usingthe standard formulae, x = D C ( z ) cos δ cos β ,y = D C ( z ) cos δ sin β , (41) z = D C ( z ) sin δ cos β , (42)where D C is the comoving distance calculated as-suming the same fiducial cosmology as in the previ-ous step. In Fig. 4 we show both the angular pro-jection of the SDSS survey area (top panel), andthe “physical length” distribution of the catalogue,with the three bins shown in different colors (bot-tom panel);(4a) we define a Cartesian grid covering the SDSS sur-vey area. Each point in this grid will be the centerof a cell which will be used for the count-in-cell pro-cedure. In order to guarantee a full coverage of thearea, considering that the volumes are spheres, thedistance between each point (i.e. each center of eachcell) is √ R ;(4a) the SDSS survey does not cover the sky in a uni-form way; thus, we are not going to include in ouranalysis all the points from the previously definedCartesian grid. Actually, converting the Cartesiancoordinates of each point back to celestial coordi-nates, we only select those ones which fall in therange defined by the initial cluster catalog. This stepwill acquire more significance as further steps in ourprocedure will be highlighted in the next points.Having defined this initial set up, we have performedthe following operations on/with it:(1b) we randomize the grid defined in the previous point (3 a ) , by shifting its origin by a random number in - α (°) δ ( ° ) Fig. 4 SDSS survey area. (Top panel.) Full clusters catalog in Celestial coordinates, right ascension α and declination δ . (Bottom panel.) Full clusters catalog in physical distances from three different viewpoints. The different colors represent thethree redshift bins we have described in the text, with: . ≤ z < . as dark grey; . ≤ z < . as middle grey; and . ≤ z < . as light grey. the interval [0 , √ R ] . Actually, we apply this shiftindependently to each Cartesian axis of the grid;(2b) after each shift, we only retain cells which satisfythe condition described in the previous point (5 a ) .This assures us that any cell which should fall outof the SDSS survey will not be considered in thefollowing analysis;(3b) additionally, we apply a jackknife cut of of cellsselected randomly;(4b) the previous steps, (1 b ) − (2 b ) − (3 b ) , are realized times, each time having a different shift in the originof the grid system (point (1 b ) ), a possibly differentnumber of cells falling out of the SDSS survey are(point (2 b ) ), and a different set of cells cut frompoint (3 b ) .After all these steps are performed, we eventually haveat our disposal the grid of centers of cells/volumes, foreach cell size R provided at point (1 a ) , and divided inthe three subgroups following the redshift criterion de-scribed at point (2 a ) , which will serve to perform thecount-in-cell with the cluster catalog data. Basically, foreach point in the grid, we will retain only the clusters whose distance from a given grid point (i.e. cell cen-ter) will be lower than the size R taken into account.Because of point (4 b ) above, we basically have dif-ferent grid sets to compare with the original catalogue,namely, we will have F ( N ) spanning the ranges [ N min , N max ] where, as before, N min and N max are theminimum and maximum numbers of objects (clusters)within each cell/volume.The median of such ensemble of values for F ( N ) will serve as observational estimation for this quantityin the χ ls analysis, as defined in Eq. (39). Instead, themedian and the standard deviation, assumed as error σ i , will be used when performing the χ jk analysis, asdefined in Eq. (40).3.3 ResultsThe complete outcomes of our statistical analysis arepresented in Table 1 for galaxies and Table 2 for clustersof galaxies. As a first comment, valid for both cases, letus notice that the softening parameter ˜ ǫ is not presentin the tables because it is basically unconstrained, ex- hibiting a uniform distribution all over the given priorrange [0 , . Thus, we can infer that its role in the fit ismarginal. For what concerns galaxies, we can see how the param-eter ¯ N is very well constrained, and there is no differ-ence between using χ ls or χ jk . This statement is validwhatever radius R and redshift bin is considered. Ac-tually, we must remember that the catalog is completeonly up to z = 0 . , after which we do have a bias to-ward smaller clusters, namely, clusters with a lower ¯ N .And this is effectively seen in the tables, where ¯ N risesslightly with redshift, becomes practically constant, andthen starts to decrease for z > . .The same consideration is basically true also for theclustering parameter b , even if in this case we can seehow there is a trend with varying R , while no peculiarchange with redshift z can be detected. In general, thisparameter lies in the range [0 . , . for all cases; how-ever, while for the smaller R we can see how it is verywell constrained, for growing R , we are able to put onlyupper limits on it, which means that also small valuesof b are compatible with data. As described in [63], weknow that smaller values of b correspond to clusterswhich have not yet fully virialized. It is worth noticingthat such trend is on the opposite side with respect to[19, 20, 63] where smaller cells volumes correspond tosmaller values of b ; although we must also stress that in[19, 20, 63], the considered physical volumes are muchlarger than those we have defined here as well as theredshift bins, and the number of objects in each bin ismuch lower. In our case the bins are much finer, bothin radius and in redshift.For what concerns the only fully characterizing pa-rameter of the nonlocal gravitational potential, i.e. ˜ λ ,we can see how it is remarkably very well constrainedin a specific and well defined range. Statistically speak-ing, we can see how it approximately lies in the range [0 . , . , which, in physical units, corresponds roughlyto [0 . , . Mpc. We can detect a slight decrease forvarying redshift, with smaller values at larger z ; whileit seems to be more pronounced a trend with varying R ,with larger ˜ λ for larger clusters. It is quite striking tohave such values for this parameter, which are actuallytotally consistent with the size of gravitational struc-tures/volumes under scrutiny. At the same time, theyare also very different from a stable GR limit, whichwould be achieved when λ → ∞ . Actually, from Eq. 7,one can see how it would be possible to recover an un-stable GR-like limit even for λ → r ij . Finally, from Figs. 5 - 6 - 7, we can also visuallycheck the quality of the fits: the larger is the cell volume,the better is the fit. In the case of clusters, results are a bit more fuzzy.First of all, we focus on ¯ N . We can see here a trendboth in redshift and in volume size R : ¯ N decreases withredshift and increases with R . The latter correspon-dence is of course expected, as large cells can containmore clusters; the former one is also somehow expected,because we might predict an intrinsic lack of clusters athigher redshifts, as they are caught in an epoch of for-mation, while a larger number of them is observable atsmaller redshifts. Notice that our bins all have the samecomoving volumes, so that any geometry-related issueshould not be effective in this case.The behavior of the clustering parameter b is nowfully consistent with literature [19, 20, 63], with smallervalues for smaller volumes, as well as for higher red-shifts.The nonlocal characteristic length ˜ λ , instead, in thiscase is quite unconstrained. We are unable to put strongboundaries on it, and we can also see how there is atendency to prefer large values, thus approaching theGR limit. Nonlocality effects can occur at large distances as theconsequence of modified gravitational potential. Specif-ically, nonlocal gravity, besides the effects at cosmolog-ical scale as a possible engine for accelerated expansion[39, 42], is relevant also for large scale formation. Inparticular, it is capable of triggering the clustering ofgalaxies and clusters of galaxies.In view of this statement, we have analyzed the ef-fects of nonlocal gravitational potential on the struc-ture formation comparing the model with observationaldata. This has been done by treating both galaxiesand clusters of galaxies (in two distinct analysis) aspoints in a statistical mechanical system, and then an-alyzing their clustering using a gravitational partitionfunction. This function has been modified accordingto the nonlocal gravitational potential and its influ-ence on the properties of the gravitational clusteringhas been scrutinized. Finally, the model has also beencompared with observational data, and we have demon-strated that nonlocality is consistent with the observa-tions, at least with those concerning the galactic scale,i.e. within volumes corresponding to clusters of galaxies. Table 1 Galaxies: results from fitting Eq. (29) with data from the cluster catalog from [60]. Data and fits are divided inredshift bins (first column) and cell radius bin (as indicated by R intervals). N cl is not a fitting parameter, but the number ofclusters in each redshift and radius bin. All the fitting parameters are described in the text. . < R < . Mpc . < R < . Mpc . < R < . Mpc z N cl ¯ N b ˜ λ N cl ¯ N b ˜ λ N cl ¯ N b ˜ λχ ls [0 . , . 1] 291 9 . +0 . − . . +0 . − . . +0 . − . 390 11 . +0 . − . . +0 . − . . +0 . − . 298 14 . +0 . − . < . 30 0 . +0 . − . [0 . , . 15] 1052 10 . +0 . − . . +0 . − . . +0 . − . . +0 . − . < . 36 0 . +0 . − . 949 15 . +0 . − . < . 26 0 . +0 . − . [0 . , . 2] 2343 10 . +0 . − . . +0 . − . . +0 . − . . +0 . − . < . 37 0 . +0 . − . . +0 . − . < . 23 0 . +0 . − . [0 . , . 25] 3069 10 . +0 . − . . +0 . − . . +0 . − . . +0 . − . < . 35 0 . +0 . − . . +0 . − . < . 30 0 . +0 . − . [0 . , . 3] 3318 10 . +0 . − . . +0 . − . . +0 . − . . +0 . − . < . 39 0 . +0 . − . . +0 . − . < . 33 0 . +0 . − . [0 . , . 35] 4024 10 . +0 . − . . +0 . − . . +0 . − . . +0 . − . < . 39 0 . +0 . − . . +0 . − . < . 34 0 . +0 . − . [0 . , . 4] 4379 10 . +0 . − . . +0 . − . . +0 . − . . +0 . − . < . 39 0 . +0 . − . . +0 . − . < . 35 0 . +0 . − . [0 . , . 45] 4925 10 . +0 . − . . +0 . − . . +0 . − . . +0 . − . < . 41 0 . +0 . − . . +0 . − . < . 37 0 . +0 . − . [0 . , . 5] 3982 9 . +0 . − . . +0 . − . . +0 . − . . +0 . − . < . 44 0 . +0 . − . . +0 . − . < . 37 0 . +0 . − . [0 . , . 55] 2838 8 . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . < . 38 0 . +0 . − . [0 . , . 6] 1971 8 . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . < . 40 0 . +0 . − . [0 . , . 65] 1046 8 . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . 977 11 . +0 . − . < . 43 0 . +0 . − . χ jk [0 . , . 1] 291 9 . +0 . − . . +0 . − . . +0 . − . 390 11 . +0 . − . < . 38 0 . +0 . − . ∗ . +0 . − . < . 02 47666(0 . . , . 15] 1052 10 . +0 . − . . +0 . − . . +0 . − . . +0 . − . < . 34 0 . +0 . − . 949 15 . +0 . − . < . 23 0 . +0 . − . [0 . , . 2] 2343 10 . +0 . − . . +0 . − . . +0 . − . . +0 . − . < . 40 0 . +0 . − . . +0 . − . < . 15 0 . +0 . − . [0 . , . 25] 3069 10 . +0 . − . < . 37 0 . +0 . − . . +0 . − . < . 35 0 . +0 . − . . +0 . − . < . 26 0 . +0 . − . [0 . , . 3] 3318 10 . +0 . − . . +0 . − . . +0 . − . . +0 . − . < . 37 0 . +0 . − . . +0 . − . < . 33 0 . +0 . − . [0 . , . 35] 4024 10 . +0 . − . . +0 . − . . +0 . − . . +0 . − . < . 37 0 . +0 . − . . +0 . − . < . 32 0 . +0 . − . [0 . , . 4] 4379 10 . +0 . − . . +0 . − . . +0 . − . . +0 . − . < . 35 0 . +0 . − . . +0 . − . < . 34 0 . +0 . − . [0 . , . 45] 4925 10 . +0 . − . . +0 . − . . +0 . − . . +0 . − . < . 42 0 . +0 . − . . +0 . − . < . 37 0 . +0 . − . [0 . , . 5] 3982 9 . +0 . − . . +0 . − . . +0 . − . . +0 . − . < . 43 0 . +0 . − . . +0 . − . < . 35 0 . +0 . − . [0 . , . 55] 2838 9 . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . < . 39 0 . +0 . − . [0 . , . 6] 1971 9 . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . < . 35 0 . +0 . − . [0 . , . 65] 1046 8 . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . 977 11 . +0 . − . < . 38 0 . +0 . − . Thus, it is possible to conclude that the nonlocal grav-itational effects can be revealed at such scales. On theother hand, the clustering of larger volumes of space,in which the points are the clusters of galaxies, do notshow a clear evidence in favor of a nonlocal modificationto gravity. For a more general discussion, it has be noted thatother modifications of gravitational potential, explain-ing dark matter or dark energy effects, have been stud-ied considering modifications of the gravitational parti-tion function. For example, it is well-known that f ( R ) gravity modifies the large distance behavior of gravita-tional potential [64–67]. Thus, the clustering of galax- Table 2 Clusters: results from fitting Eq. (29) with data from the cluster catalog from [60]. Data and fits are divided inredshift bins (first column) and cell radius bin (as indicated by R intervals). N cl is not a fitting parameter, but the number ofclusters in each redshift and radius bin. All the fitting parameters are described in the text. R = 10 Mpc R = 20 Mpc R = 30 Mpc R = 40 Mpc z N cl ¯ N b ˜ λ ¯ N b ˜ λ ¯ N b ˜ λ ¯ N b ˜ λχ ls [0 . , . ∗ . +0 . − . < . 20 482 (370) 0 . +0 . − . . +0 . − . < 544 1 . +0 . − . . +0 . − . . +0 . − . . +0 . − . +684 − [0 . , . . +0 . − . < . 29 312 (14) 0 . +0 . − . < . 16 161 (9) 1 . +0 . − . . +0 . − . 28 (17) 2 . +0 . − . . +0 . − . . , . 42] 20634 0 . +0 . − . < . 27 1430 (1022) 0 . +0 . − . < . 10 0 . 97 (1 . 04) 0 . +0 . − . . +0 . − . 301 (0 . 71) 2 . +0 . − . . +0 . − . 617 (506) χ jk [0 . , . ∗ . +0 . − . . +0 . − . 37 (30) 0 . +0 . − . . +0 . − . . +0 . − . < . 31 2 . 43 (0 . 66) 3 . +0 . − . . +0 . − . . , . . +0 . − . . +0 . − . 58 (5 . 6) 0 . +0 . − . . +0 . − . 178 (905) 1 . +0 . − . . +0 . − . 527 (201) 2 . +0 . − . . +0 . − . . . , . 42] 20634 0 . +0 . − . < . < . 65 0 . +0 . − . < . < . 33 0 . +0 . − . . +0 . − . . +0 . − . . +0 . − . 81 (5 . < R < < z < N p ( N ) < R < < z < N p ( N ) < R < < z < 10 15 200.000.050.100.15 N p ( N ) < R < < z < N p ( N ) < R < < z < N p ( N ) < R < < z < N p ( N ) < R < < z < N p ( N ) < R < < z < N p ( N ) < R < < z < N p ( N ) < R < < z < N p ( N ) < R < < z < N p ( N ) < R < < z < N p ( N ) Fig. 5 Comparison between data and theoretical expectation Eq. (29) for cell size bin . < R < . Mpc . Black points: data;black bars: jackknife-like observational errors. Solid red line: best fit from minimization of χ ls ; solid green line: best fit fromminimization of χ jk .4 < R < < z < N p ( N ) < R < < z < 10 15 20 250.000.020.040.060.080.100.120.14 N p ( N ) < R < < z < 10 15 20 250.000.020.040.060.080.100.120.14 N p ( N ) < R < < z < 10 15 20 250.000.020.040.060.080.100.12 N p ( N ) < R < < z < N p ( N ) < R < < z < 10 15 200.000.050.100.15 N p ( N ) < R < < z < 10 15 20 250.000.020.040.060.080.100.120.14 N p ( N ) < R < < z < N p ( N ) < R < < z < N p ( N ) < R < < z < N p ( N ) < R < < z < N p ( N ) < R < < z < N p ( N ) Fig. 6 Comparison between data and theoretical expectation Eq. (29) for cell size bin . < R < . Mpc . Black points: data;black bars: jackknife-like observational errors. Solid red line: best fit from minimization of χ ls ; solid green line: best fit fromminimization of χ jk . ies, interacting through effective potentials derived from f ( R ) gravity, has also been studied using a modifiedgravitational partition function, and it has been ob-served that also these models are consistent with ob-servations for large samples of galaxies and clusters[68, 69], although the corresponding statistical infer-enec is a bit weaker that the nonlocal scenario we haveconsidered here. It is worth noticing that the clusteringof galaxies in f ( R ) gravity is also consistent with theconstraints coming from the Planck data [70].The MOND gravitational partition function has alsobeen used to study the thermodynamics and the clus-tering properties of systems of galaxies [13, 71]. Furthermore, the gravitational phase transition canalso be analyzed using the gravitational partition func-tion for systems of galaxies [14]. In fact, a first orderphase transition occurs due to the clustering of galax-ies from a homogeneous phase. In this framework, it ispossible to take into account the Yang-Lee theory forsystems of galaxies, and use the complex fugacity to an-alyze the phase transition [72]. A forthcoming step willbe to extend the Yang-Lee theory for nonlocal gravi-tational potential, and analyze the phase transition forthe related gravitational partition function.In particular, it has to be noted that the clusteringof galaxies has also been investigated using the cos-mic energy equation [73]. As the cosmic energy equa- < R < < z < 10 15 20 250.000.020.040.060.080.100.12 N p ( N ) < R < < z < 10 15 20 250.000.020.040.060.080.10 N p ( N ) < R < < z < 10 15 20 25 300.000.020.040.060.080.10 N p ( N ) < R < < z < 10 15 20 25 300.000.020.040.060.080.10 N p ( N ) < R < < z < 10 15 20 25 300.000.020.040.060.080.10 N p ( N ) < R < < z < 10 15 20 25 300.000.020.040.060.080.10 N p ( N ) < R < < z < 10 15 20 25 300.000.020.040.060.080.10 N p ( N ) < R < < z < 10 15 20 25 300.000.020.040.060.080.100.12 N p ( N ) < R < < z < 10 15 20 250.000.020.040.060.080.100.12 N p ( N ) < R < < z < 10 15 20 250.000.020.040.060.080.100.120.14 N p ( N ) < R < < z < 10 15 20 250.000.050.100.15 N p ( N ) < R < < z < N p ( N ) Fig. 7 Comparison between data and theoretical expectation Eq. (29) for cell size bin . < R < . Mpc . Black points: data;black bars: jackknife-like observational errors. Solid red line: best fit from minimization of χ ls ; solid green line: best fit fromminimization of χ jk . tion is derived by approximating galaxies as points ina statistical mechanical system, a softening parame-ter can modify the cosmic energy equation [74]. It hasbeen demonstrated that a large scale modification ofthe gravitational potential also modifies the cosmic en-ergy equation [75].As final remark, we can say that, according to theresults reported here, nonlocal gravity affects the clus-tering process. In other words, the size and the distri-bution of large structures as galaxies and clusters ofgalaxies, may contain a sort of intrinsic “signature” ofnonlocality at infrared scales. In some sense, the non-locality invoked at ultraviolet scales to match gravityand quantum mechanics at fundamental level [44] may have a cosmological counterpart capable of triggeringthe structure formation. Acknowledgements SC acknowledgs the support of the Istituto Nazionaledi Fisica Nucleare (INFN), sezione di Napoli, iniziativespecifiche QGSKY and MOONLIGHT-2. References 1. C. S. Frenk et al. , Astrophys. J. , 554 (1999),arXiv:astro-ph/9906160 . R = 10 Mpc0.05 < z < N p ( N ) R = 10 Mpc0.281 < z < N p ( N ) R = 10 Mpc0.361 < z < N p ( N ) R = 20 Mpc0.05 < z < N p ( N ) R = 20 Mpc0.281 < z < N p ( N ) R = 20 Mpc0.361 < z < N p ( N ) R = 30 Mpc0.05 < z < N p ( N ) R = 30 Mpc0.281 < z < N p ( N ) R = 30 Mpc0.361 < z < N p ( N ) R = 40 Mpc0.05 < z < N p ( N ) R = 40 Mpc0.281 < z < N p ( N ) R = 40 Mpc0.361 < z < N p ( N ) Fig. 8 Comparison between data and theoretical expectation Eq. (29) for cell size bin . < R < . Mpc . Black points:data; black bars: jackknife-like observational errors. Solid red line: best fit from minimization of χ ls ; solid green line: best fitfrom minimization of χ jk . 2. D. J. Barnes, S. T. Kay, Y. M. Bahé, C. DallaVecchia, I. G. McCarthy, J. Schaye, R. G.Bower, A. Jenkins, P. A. Thomas, M. Schaller,R. A. Crain, T. Theuns, and S. D. M. White,Mon. Not. Roy. Astron. Soc. , 1088 (2017),arXiv:1703.10907 [astro-ph.GA] .3. T. J. Ponman, A. J. R. Sander-son, and A. Finoguenov,Mon. Not. Roy. Astron. Soc. , 331 (2003),arXiv:astro-ph/0304048 .4. G. M. Voit, Rev. Mod. Phys. , 207 (2005),arXiv:astro-ph/0410173 .5. W. C. Saslaw and A. Yang, arXive-prints , arXiv:0902.0747 (2009), arXiv:0902.0747 [astro-ph.CO] .6. W. C. Saslaw and F. Fang,Astrophys. J. , 16 (1996).7. W. C. Saslaw and R. K. Sheth,Astrophys. J. , 504 (1993).8. H. Rahmani, W. C. Saslaw, andS. Tavasoli, Astrophys. J. , 1121 (2009),arXiv:0902.1103 [astro-ph.CO] .9. M. Itoh, S. Inagaki, and W. C. Saslaw,Astrophys. J. , 476 (1993).10. W. C. Saslaw, S. M. Chitre, M. Itoh, and S. Ina-gaki, Astrophys. J. , 419 (1990).11. W. C. Saslaw, Astrophys. J. , 11 (1986). 12. W. C. Saslaw and F. Ahmad,Astrophys. J. , 1246 (2010).13. S. Upadhyay, B. Pourhassan, and S. Capozziello,Int. J. Mod. Phys. D , 1950027 (2018),arXiv:1809.03579 [gr-qc] .14. M. S. Khan and M. A. Malik,Mon. Not. Roy. Astron. Soc. , 2629 (2012).15. F. Ahmad and Hameeda,Astrophysics and Space Science , 227 (2010).16. M. Hameeda, B. Pourhassan, M. C. Rocca, andA. B. Brzo, Eur. Phys. J. C , 146 (2021).17. W. C. Saslaw and A. J. S. Hamilton,Astrophys. J. , 13 (1984).18. F. Ahmad, W. C. Saslaw, and N. I. Bhat,Astrophys. J. , 576 (2002).19. G. R. Sivakoff and W. C.Saslaw, Astrophys. J. , 795 (2005),arXiv:astro-ph/0503509 [astro-ph] .20. H. Rahmani, W. C. Saslaw, andS. Tavasoli, Astrophys. J. , 1121 (2009),arXiv:0902.1103 [astro-ph.CO] .21. F. Ahmad, W. C. Saslaw, and M. A. Malik,Astrophys. J. , 940 (2006).22. A. G. Riess et al. (Supernova SearchTeam), Astron. J. , 1009 (1998),arXiv:astro-ph/9805201 .23. S. Perlmutter et al. (Supernova Cosmol-ogy Project), Astrophys. J. , 565 (1999),arXiv:astro-ph/9812133 .24. N. Aghanim et al. (Planck),Astron. Astrophys. , A6 (2020),arXiv:1807.06209 [astro-ph.CO] .25. V. G. Gurzadyan, A. A. Kocharyan, andA. Stepanian, Eur. Phys. J. C , 24 (2020),arXiv:2001.02634 [astro-ph.CO] .26. V. G. Gurzadyan and A. Stepa-nian, Eur. Phys. J. C , 869 (2018),arXiv:1810.09846 [physics.gen-ph] .27. V. G. Gurzadyan and A. Stepa-nian, Eur. Phys. J. C , 169 (2019),arXiv:1902.07171 [physics.gen-ph] .28. V. G. Gurzadyan and R. Pen-rose, Eur. Phys. J. Plus , 22 (2013),arXiv:1302.5162 [astro-ph.CO] .29. V. G. Gurzadyan and A. Stepa-nian, Eur. Phys. J. C , 568 (2019),arXiv:1905.03442 [astro-ph.CO] .30. M. Le Delliou, R. J. F. Mar-condes, and G. a. B. L. Neto,Mon. Not. Roy. Astron. Soc. , 1944 (2019),arXiv:1811.10712 [astro-ph.CO] .31. P. Boonserm, T. Ngampitipan, A. Simpson,and M. Visser, Phys. Rev. D , 024050 (2020), arXiv:1909.06755 [gr-qc] .32. S. Peirani and J. A. D. F. Pacheco,Astron. Astrophys. , 845 (2008),arXiv:0806.4245 [astro-ph] .33. I. T. Iliev and P. R. Shapiro,Mon. Not. Roy. Astron. Soc. , 468 (2001),arXiv:astro-ph/0101067 .34. A. Stark, C. J. Miller, andD. Huterer, Phys. Rev. D , 023543 (2017),arXiv:1611.06886 [astro-ph.CO] .35. M. Aguena and M. Lima,Phys. Rev. D , 123529 (2018),arXiv:1611.05468 [astro-ph.CO] .36. D. Wen, A. J. Kemball, and W. C.Saslaw, Astrophys. J. , 160 (2020),arXiv:2001.06119 [astro-ph.CO] .37. M. Hameeda, S. Upadhyay,M. Faizal, and A. F. Ali,Mon. Not. Roy. Astron. Soc. , 3699 (2016),arXiv:1610.07990 [gr-qc] .38. B. Pourhassan, S. Upadhyay,M. Hameeda, and M. Faizal,Mon. Not. Roy. Astron. Soc. , 3166 (2017),arXiv:1704.06085 [gr-qc] .39. S. Deser and R. P. Woodard,Phys. Rev. Lett. , 111301 (2007),arXiv:0706.2151 [astro-ph] .40. J.-C. Ding and J.-B.Deng, JCAP , 054 (2019),arXiv:1908.11223 [astro-ph.CO] .41. S. Bahamonde, S. Capozziello, and K. F. Di-alektopoulos, Eur. Phys. J. C , 722 (2017),arXiv:1708.06310 [gr-qc] .42. S. Bahamonde, S. Capozziello, M. Faizal, andR. C. Nunes, Eur. Phys. J. C , 628 (2017),arXiv:1709.02692 [gr-qc] .43. F. Bajardi, S. Capozziello, andD. Vernieri, Eur. Phys. J. Plus , 942 (2020),arXiv:2011.01317 [gr-qc] .44. L. Modesto and L. Rachwał,Int. J. Mod. Phys. D , 1730020 (2017).45. S. Capozziello, M. Capriolo, and S. No-jiri, Phys. Lett. B , 135821 (2020),arXiv:2009.12777 [gr-qc] .46. E. Belgacem, Y. Dirian, A. Finke, S. Foffa,and M. Maggiore, JCAP , 022 (2019),arXiv:1907.02047 [astro-ph.CO] .47. A. S. Koshelev, JHEP , 029 (2007),arXiv:hep-th/0701103 .48. A. S. Koshelev and S. Y. Ver-nov, Eur. Phys. J. C , 2198 (2012),arXiv:0903.5176 [hep-th] . 49. S. Gielen, D. Oriti, and L. Sin-doni, Phys. Rev. Lett. , 031301 (2013),arXiv:1303.3576 [gr-qc] .50. S. Gielen, JCAP , 013 (2019),arXiv:1811.10639 [gr-qc] .51. N. Yunes and X. Siemens,Living Rev. Rel. , 9 (2013),arXiv:1304.3473 [gr-qc] .52. C. Chicone and B. Mashhoon,Class. Quant. Grav. , 075005 (2016),arXiv:1508.01508 [gr-qc] .53. S. X. Tian and Z.-H. Zhu,Phys. Rev. D , 124059 (2019).54. K. F. Dialektopoulos, D. Borka, S. Capozziello,V. Borka Jovanović, and P. Jo-vanović, Phys. Rev. D , 044053 (2019),arXiv:1812.09289 [astro-ph.GA] .55. L. Amendola, Y. Dirian, H. Nersisyan,and S. Park, JCAP , 045 (2019),arXiv:1901.07832 [astro-ph.CO] .56. F. W. Hehl and B. Mashhoon,Phys. Rev. D , 064028 (2009),arXiv:0902.0560 [gr-qc] .57. F. W. Hehl and B. Mash-hoon, Phys. Lett. B , 279 (2009),arXiv:0812.1059 [gr-qc] .58. H.-J. Blome, C. Chicone, F. W. Hehl, andB. Mashhoon, Phys. Rev. D , 065020 (2010),arXiv:1002.1425 [gr-qc] .59. C. Chicone and B. Mashhoon,J. Math. Phys. , 072501 (2016),arXiv:1510.07316 [gr-qc] .60. Z. L. Wen, J. L. Han, and F. S. Liu,Astrophysical J. Supplement , 34 (2012),arXiv:1202.6424 [astro-ph.CO] .61. Z. L. Wen and J. L. Han,Astrophys. J. , 178 (2015),arXiv:1506.04503 [astro-ph.GA] .62. A. Yang and W. C. Saslaw,Astrophys. J. , 123 (2011),arXiv:1009.0013 [astro-ph.CO] .63. A. Yang and W. C. Saslaw,Astrophys. J. , 113 (2012),arXiv:1203.3823 [astro-ph.CO] .64. S. Capozziello, Int. J. Mod. Phys. D , 483 (2002),arXiv:gr-qc/0201033 .65. S. Capozziello and M. De Lau-rentis, Phys. Rept. , 167 (2011),arXiv:1108.6266 [gr-qc] .66. T. P. Sotiriou and V. Faraoni,Rev. Mod. Phys. , 451 (2010),arXiv:0805.1726 [gr-qc] . 67. S. Nojiri, S. D. Odintsov, and V. K.Oikonomou, Phys. Rept. , 1 (2017),arXiv:1705.11098 [gr-qc] .68. S. Capozziello, E. De Filippis, and V. Salzano,Mon. Not. Roy. Astron. Soc. , 947 (2009),arXiv:0809.1882 [astro-ph] .69. S. Capozziello, M. Faizal, M. Hameeda,B. Pourhassan, V. Salzano, and S. Upadhyay,Mon. Not. Roy. Astron. Soc. , 2430 (2018),arXiv:1711.06630 [astro-ph.CO] .70. I. De Martino, M. De Laurentis,F. Atrio-Barandela, and S. Capozziello,Mon. Not. Roy. Astron. Soc. , 921 (2014),arXiv:1310.0693 [astro-ph.CO] .71. M. Hameeda, B. Pourhassan, M. Faizal,C. P. Masroor, R. U. H. Ansari, andP. K. Suresh, Eur. Phys. J. C , 769 (2019),arXiv:1911.01739 [gr-qc] .72. M. S. Khan and M. A. Malik,Astrophysics and Space Science , 211 (2013).73. A. Wahid, F. Ahmad, and A. Nazir,Astrophysics and Space Science , 241 (2011).74. F. Ahmad, A. Wahid, M. A. Malik, and S. Masood,International Journal of Modern Physics D , 119 (2009).75. M. Hameeda, S. Upadhyay, M. Faizal,A. F. Ali, and B. Pourhas-san, Phys. Dark Univ.19