Table Of Content

Cumulant Correlators in 2D and 3D Scale Free Simulations Dipak Munshi1 and Adrian L. Melott2 1 Astronomy Unit, Queen Mary and Westfield College, London E1 4NS, United Kingdom. 2 Department of Physics and Astronomy, University of Kansas, Lawrence, Kansas 66045,U.S.A. Email: [email protected], [email protected] Abstract Shape dependence of higher order correlations introduces complication in direct determination of these 8 quantities. For this reason theoretical and observational progress has been restricted in calculating one 9 9 point distribution functions and related moments. Methods based on factorial moments of two point count 1 probabilitydistribution(cumulantcorrelators)wererecentlyshowntobeefficientinsubtractingdiscreteness effects and extracting useful information from galaxy catalogs. We use these cumulant correlators Q to n NM a study clustering in scale free simulations both in two and three dimensions. Using this method in the J highly nonlinear regime we were able to separate hierarchal amplitudes Q ,R and R associated with 3 a b 3 different tree graphs contributing to third and fourth order correlation functions. They were found to increase with power on large scales. Results based on factorial moments of one point cumulants show very 1 good agreement with that on cumulant correlators. Comparison of simulation results with perturbation v theory and extended perturbationtheory were found to be in reasonableagreement. Comparisonswere also 1 1 made against predictions from hierarchal models for higher order correlations. We argue that finite volume 0 corrections are very important for computation of cumulant correlators. 1 0 SubjectHeadings:largescalestructureoftheuniverse-galaxies: statistics-methods: dataanalysis -methods: 8 9 simulation / h 1. Introduction p - o There is growing evidence which suggest that the large scale structure in the universe was formed by the r t gravitational amplification of small inhomogeneities. Statistical characterization of clustering is important s a for understandingthe dynamics ofgravitationalclustering. Correlationfunctions andcountin cellstatistics : v are among the oldest and most widely used statistics to quantify clustering. Evaluation of higher order i correlationfunctions with full shape dependence is difficult given limited size of presentday galaxycatalogs X and numerical simulations (Fry & Peebles 1978, Peebles 1980). Much progress has therefore been made r in estimation of volume averages of these quantities (Gaztanaga 1992, Bouchet et al. 1993, Gaztanaga & a Frieman 1994, Colombi et al. 1996a, Szapudi et al. 1996), which can be directly related to moments of one pointcountprobabilitydistributionfunction. Analyticalmethodsbasedontreelevelperturbationtheoryand loopcorrectionswere used to evaluate these quantities in the quasi-linearregime (Peebles 1980,Juszkiewicz et al. 1993, Bernardeau 1992, Bernardeau 1994, Munshi et al. 1994, Colombi et al. 1992, Scoccimarro & Frieman 1996a-b,). In the highly nonlinear regime, although no full theory exists, scaling models are often used to predict their behaviour (Peebles 1980, Balian & Schaeffer 1989). Computational methods were also developed based on factorial moments to subtract Poisson noise from discrete data. Errors associated with these one point cumulants and other related quantities such as void probability function have also been estimated and correction procedure have been developed (Colombi et al. 1995, Szapudi & Colombi 1996, Munshi et al. 1997). Although one point quantities carry useful information about dynamics and background geometry av- eraging associated with such quantities causes a significant loss of information. Alternative methods use correlationof pairs of cells (Szapudi et al. 1992, Meiksin et al. 1992,Szapudi et al. 1995). Cumulant corre- latorsprovideanaturalgeneralizationofonepointcumulants. Thehierarchalansatzinthestrongclustering regime and perturbation theory in the weak clustering regime have definite predictions for these quantities. Using cumulant correlators it is possible to separate amplitudes associated with different tree topologies up 1 tofifthorderinthecorrelationhierarchy. Szapudi&Szalay(1997)haveanalysedAPMdatausingcumulant correlators. Weusethesestatisticstostudyscale–freesimulationsintwoandthreedimensions(2Dand3D). In the next section we outline the theoretical framework necessary for using cumulant correlators, predictions from perturbation theory, hierarchal ansatz and extended perturbation theory. In section §3 we describe our numerical simulations and techniques to evaluate cumulant correlators from simulation data. Wediscussimplicationsofourresultsinsection§4andcompareitwithexistingresultsfromgalaxycatalogs. 2. Theoretical Predictions We define factorial moment correlatorsof a pair of cells of volume l3 separated by a distance r as h(N ) (N )i−h(N ) ih(N )i 1 k 2 l 1 k 2 l w (l,r)= ; k 6=0,l6=0, (1) kl hNik+l and the normalized one point factorial moment h(N ) i 1 k w (l)= , (2) k0 hN ik 1 where we have used the notation (N) = N(N −1)...(N −k +1), and the angular braces h...i denote k average over all possible positions of the cells. Itispossibletorelatew andw tonormalizedonepointcumulantsS =hδNi/hδ2i(N−1) andcumulant k0 kl N correlators C = hδN(x )δM(x )i/hδ(x )δ(x )ihδ2i(N+M−2). However we find it to easier to work with NM 1 2 1 2 Q =S /N(N−2)andQ =C /N(N−1)M(M−1). Whilecentralmomentscanalsobeusedtocomputed N N NM NM these quantities, factorial moments are known to be more suitable for subtracting Poissonshot noise. One can introduce the generating function for the factorial moments in terms of the cumulants Q . N ∞ W(x)=exp Γ xNQ (3) X N N N=1 where we have defined Γ =NN−2ξ N−1/N! and ξ is the volume average of ξ over volume of the cell l3. N s s 2 The generating function can be linked with the one point void probability distribution function. In generalQ parametersshowscaledependence,increasingfromquasi-linearphasetohighlynonlinearphase. N Hierarchal form of correlation functions demand Q to be independent of scale in the nonlinear regime. N We can also construct generating function of factorial moments W(x,y) = ∞W xMyN/m!n! which P0 mn can be related to generating function of cumulant correlatorsQ(x,y) ∞ Q(x,y)=ξ xMyNQ Γ Γ NM, (4) l X NM M M M=1,N=1 by the following equation W(x,y)=W(x)W(y)(expQ(x,y)−1). (5) Expanding equation (5) one can compute cumulant correlatorsQ from factorial moments W . NM NM Q 2Γ Γ ξ =w /2−ξ 12 1 2 r 12 r Q 3Γ Γ ξ =w /6−w /2−w /2+ξ 13 1 3 r 13 12 20 r Q Γ24ξ =w /4−w +ξ −ξ2/2 (6) 22 2 r 22 12 r r 2 As it is clear from equations (6), cumulant correlators are determined by two different length scales, length scale correspondingto ξ and scale correspondingto ξ . We use the terms quasi-linear and nonlinear s r depending upon values of ξ . s 2.1 Nonlinear Regime The extreme nonlinear stage of gravitational clustering is believed to be well described by a power law solution for correlationfunction ξ =r−γ, where γ can be expressedas a function of initial power spectrum 2 (Davis & Peebles 1977, Balian & Schaeffer 1989. All higher order correlations exhibit a self similar scaling in strong clustering regime ξ (λr ,...λr )=λ−γ(N−1)ξ (r ,...r ). (7) N 1 N N 1 N MoreexplicitfunctionalformsforN-pointcorrelationsareconstructedbysummingoverproductsofN−1 two-pointcorrelationfunctionscorrespondingtodifferenttopologies,eachofwhichrepresentingatreegraph spanning N-points with amplitudes T n,α (N−1) ξ (r ,...r )= T ξ . (8) N 1 N X N,α X Y (ri,rj) α,N−trees labelings edges Tree graphsspanning three points can haveonly one topology. The amplitude associatedwith this graph Q contributes through three different configurations of these points, so that cumulant correlators to have 3 following form: hδ(x )2δ(x )i=2Q ξ ξ =Q (2ξ ξ +ξ2). (9) 1 2 12 s r 3 s r r Higher order cumulants get contributions from different types of tree diagrams. At fourth order trees connectingfourpointshavetwodifferenttopologies,knownasthe snake(withitsamplitudedenotedbyR ) a and the star topology (corresponding amplitude denoted conventionally by R ). Using this notation we can b write hδ (x )3δ (x )i=9Q ξ ξ2 =6ξ ξ2R +3ξ ξ2R +6ξ2ξ R +ξ3R (10) 1 1 2 2 13 r s r s a r s b r s a r b hδ(x )2δ(x )2i=4Q ξ ξ2 =4ξ ξ2R +4ξ2ξ R +4ξ2ξ R +4ξ3R . (11) 1 2 22 r s r s a r s b r s a r b ThesesimultaneousequationsarelinearinR andR andcanbesolvedonceQ andQ aredetermined. a b 22 13 Whenξ <<1 whichwillbe the casewhentwocellsarefaraway,onecandefine linearcumulantcorrelators r (linear in ξ ) by considering only terms linear in ξ . The linear solution to equation (11) are R =Q and r r a 22 R =3Q −2Q . In such a linear accuracy one can show that b 13 22 Q ≈Q (12) NM N+M Althoughamplitudesoftreetermswithdifferenttopologiescanbeestimatedbycumulantcorrelators,the situation becomes more complicated at higher order as number of degenerate correlators become less than number of topologies, making the system indeterministic. Other interesting questions regarding scaling of generating functionss Q(x,y) and related question about the nature of Q in the highly nonlinear regime NM can be addressed when bigger N-body simulations become available. 2.2 Quasi-linear Regime Inquasi-linearregimewhenξ issmallerthanunity itis possibletoexpandδ =δ(1)+δ(2)+δ(3)...where s perturbation expansion is valid as long as the series is convergent. Using perturbative calculation of two point quantities it was possible for Bernardeau (1996) to express C at lowest order. Perturbative terms pq contributing to the expansion of hδp(x )δq(x )i can be written as 1 2 3 p q hδp(x )δq(x )i≈ h δ(pi)(x ) δ(qi)(x )i (13) 1 2 X Y 1 Y 2 decompositions i=1 j=1 where sum is taken over all possible decompositions and p and q satisfies i i p q p + q =2(p+q)−2. (14) X i X j i=1 j=1 In lowest order in perturbation theory it is possible to write (Bernardeau 1996), hδp(x )δq(x )i =C hδ2(x)ip+q−2hδ(x )δ(x )i=C ξp+q−2ξ . (15) 1 2 c pq 1 2 pq s r Using method of generating function it was shown that C can be decomposed into following relation pq (Bernardeau, 1996), C =C C . (16) pq p1 q1 In 3D the lowest order C can be expressed in terms of spectral index n, pq C =68/21−1/3(n+3) (17) 21 C =11710/441−61/7(n+3)+2/3(n+3)2. (18) 31 It should be noted that such factorization is possible only in case of tree-level diagrams also calculations were done using large separation limit. Given the whole hierarchy of C it is possible to compute bias NM associated with gravitationalclustering in the quasi-linear regime (Bernardeau 1996). 2.3 Connecting different regimes: Extended Perturbation Theory TreelevelperturbationtheorycanpredictS parametersfortop-hatsmoothing,itwasused(asoutlined N above) to compute cumulant correlatorsQ to arbitrary order in large separation limit. The expressions N,M derivedforcumulantsdependsonindexofinitialpowerspectrum. Howeveritwaslaterrealizedthatthesame quasi-linear expressions can describe evolution of S parameters from quasi-linear regime to the nonlinear N regime by allowing spectral index n to vary with scale (Colombi et al. 1996b). It is interesting that whole hierarchyofS canbedescribedbysuchavariationofeffectivespectralindexn (unfortunatelynostraight N eff forwardrelationexistsbetweentruenonlinearpowerspectraandn ). Anextensionofperturbationtheory eff to the non–perturbative regime was considered by Szapudi & Szalay (1997) for cumulant correlatorsQ . N,M Tree level perturbation theory as mentioned before predicts C = C C . Using same logic as extended pq p1 q1 perturbationtheoryonecanexpectsucharelationtoholdevenforsmallcellsizesatlargedistancesalthough separately each of these term may vary significantly. 3. Simulations and Data Analysis Thesimulationsusedherearenumericalmodelsforthegravitationalclusteringofcollisionlessparticlesinan expanding background. We study evolution of initial Gaussian perturbations in Ω=1 universe. All the 2D simulationsaredonewithaparticle-mesh(PM)codewith5122particleswithanequalnumberofgridpoints and in 3D 1283 particles with 1283 grid points The code has at least twice the dynamical resolution of any other PM code with which it has been compared. The 2D simulations are described in detail in Beacom et. al (1991), with a video of the evolution in Kauffmann and Melott (1992). The 3D simulations are described in Melott and Shandarin (1993). Both sets, which consist of multiple realizations (different random seeds) for a variety of power spectra, have been widely used for comparative studies of various statistical methods and dynamical approximations. 4 Fig. 1.— Lowerpaneldisplays factorialmomentcorrelatorsw as a function of separationr in 2D in units kl of grid scale. Degenerate parallel lines correspond to increasing values of N +M from 2 to 6. Left panels display results for scale free spectra with power law index n = −2, middle panel for n = 0 and right panel n = 2 in two dimension. Middle and upper panels plot cumulant correlators of different orders. Dashed curves in middle panels represent Q Q for different values of N and M, closest solid curves represent N1 M1 Q forsamevaluesofN andM. Accordingtopredictionsfromextendedperturbationtheorytheyshould NM overlap Q = Q Q . Top most panels compare predictions of hierarchal ansatz Q = Q . NM N1 M1 NM N+M Adjacent curves correspondto different N and M with same N +M. Error bars were calculated by finding scatter in results with different realization of same initial power spectra. 5 Fig. 2.— Same as Figure - 1 but results from only one realization are plotted to show level of agreement with theoretical predictions in individual realizations. 6 In this paper we analyze a subset of the simulations with featureless power-law initial spectra of the general form, P(k) ∝ kn for k ≤k , (19) c = 0 for k >k . (20) c We haveanalyzedpower-lawmodelswith n=2,0,−2in2Dandn=1,0,−1,−2in3Dwithacutoffineach case at the Nyquist wave number: k = 256k for 2D and k = 64k for 3D where k = 2π/L is the c f c f f box fundamental mode associated with the box size. We choose σ(k ), the epoch when the scale 2π/k is going nonlinear as a measure of time. NL NL 1 σ(k )=RkkfNyP(k)kk.2 (21) NL RkkfNLP(k)kk. The first scale to go nonlinear is the one corresponding to the Nyquist wave number. This happens, by definition,whenthevarianceσisunity. Ofcourseasσincreasessuccessivelargerscalesenterinthenonlinear regime. Thesimulationswerestoppedatλ =2l , 4l , 8l ,...., L /2. Inourstudy in2Dwetook NL grid grid grid box the epoch when L /16 for analysis of n = −2 spectrum and for n = 0 and2 we studied the epoch when box L /8 is going nonlinear. In 3D we had less dynamic range so choose the epoch when L /4 for analysis box box of n = −2 spectrum and for n = −1,0, and1 we studied the epoch when L /2 is going nonlinear. The box separation r that we study is much less compared to the length scale going nonlinear so it is unlikely that our results will be affected much by boundary conditions. The growth rate of various modes in the linear regime were studied by Melott et al. (1988) for this PM code. The results at λ= 3l are equivalent to the ones obtained by a typical PM code at λ= 8l , due grid grid to the staggered mesh scheme. So we expect that our code performs well at the wavelength associatedwith four cells and since the collapse of 4l -size perturbations will give rise to condensations of diameter 2l grid grid or less, the smallest cell size that can be safely resolved is 2l . grid Computations of count probability distribution function (CPDF) P (n) and two-point count probability l distributionP (n,m)weredonebylayingdownagridofmeshspacinglandcountingoccupationnumberin l,r eachcellto computethe probabilityoffinding nobjectsincellsize l andalsothe jointprobabilityoffinding n and m objects in two cells separated by a distance r. Statistics were improved by perturbing the grid in each orthogonal direction and keeping the mesh undistorted while repeating the counting process. Both in twoandthreedimensionsweconsideredonlycellsofsizel . Whilethisparticularchoiceofcellsizeisopen grid to question, we will show that our results match remarkably well with our earlier studies done using bigger cell sizes. With this cell size we could reachprobabilities as few times 10−6 in 2D and 3D. Computations of factorialmomentsandfactorialmomentcorrelatorsweredoneusingcomputedvaluesofP (n)andP (n,m). l l,r Computedvaluesofξ¯ forcellsusedwerefoundtobe 8.98,27.87,60.39in2Dforspectralindex−2,0, and2 2 respectively. In 3D these values were found to be 25.46,88,125.5, and 171 respectively. Different spurious effects such as shot noise and finite volume correctionsare important while computing countprobability distribution functions. For smallcell sizes Poissonnoise starts dominating,whereas larger cellsizesaredominatedby finite volumecorrections. We find w parametersto be dominatedby shotnoise kl with increasing separation of cells, this effects starts to be severe with higher order cumulants. This effect starts dominating even for small separation for spectra with less power on larger scales. Using large cells produces effect similar to smoothing the density field and hence reducing the effective correlation length scale between them. These two dominant effects reduce the range of separation which can be probed for studying cumulant correlators. Finite volume corrections are more difficult to quantify. Methods based on scaling of count in cell statistics were shown to be effective in correcting finite volume effects (Munshi et al. 1997),similar arguments canbe usedto developcorrectionsfor cumulantcorrelators. The validity of sucha method depends on correctnessof the scaling model, which we test in this paper. Developments of methods to correct Q for finite volume effect are left for future work. NM Computedw for2DsimulationsareplottedinFigure-1asafunctionofcellseparationr. Corresponding kl results for 3D are presented in Figure - 5. Parallel degenerate lines correspond to different values of N + 7 Fig. 3.— Lower order hierarchal amplitudes Q , R , R calculated from the fully nonlinear cumulative 3 a b correlators in 2D are plotted against separation r measured in grid units. Solid curve represent estimates of Q lower and upper dashed lines represent amplitudes of fourth order snake and star graphs R and R 3 a b respectively. Open triangles and square correspond measurement of Q and Q respectively from factorial 3 4 moments of counts in cell statistics after finite volume corrections were taken into account (from Munshi et al. 1997). Filled triangles and squares are measurements of same quantities using cumulant correlators without finite volume correction. Fig. 4.— Each panel displays factorial moment correlatorsw as a function of separation r in unit of grid kl scale for different initial scale–free power law spectra in 3D. Scatter in each plot is computed from scatter in four different realizations of same power spectra. 8 Fig. 5.— Nonlinear cumulant correlatorsis plotted against separationr in 3D as measured in units of grid scale. UppermostcurvecorrespondtoQ ,lowestcurvetoQ andthedashedonetoQ2 . Solidanddashed 31 22 21 straight lines correspond to predictions for Q , and Q2 from perturbations theory at large separation. 31 21 Fig. 6.— Hierarchal amplitudes Q , R , R calculated from the nonlinear cumulative correlators in 3D 3 a b are plotted against separation r measured in grid units. Solid curve represent estimates of Q lower and 3 upperdashedlinesrepresentamplitudesoffourthordersnakeandstargraphsR andR respectively. Open a b trianglesandsquarecorrespondtomeasurementofQ andQ respectivelyfromfactorialmomentsofcounts 3 4 in cell statistics after finite volume corrections were taken into account (from Munshi et al. 1997). Filled trianglesandsquaresare measurementsofsame quantities using cumulantcorrelatorswithout finite volume correction. 9 M increasing from bottom to top in each panel. Middle and top panels show variation of Q with r. NM Middle panels test the validity of extended perturbation theory. Dashed lines correspondto Q Q while N1 M1 neighboringsolid lines representQ for same values of N andM. Errorbars are computed by estimating NM the scatterin resultsfromfourdifferentrealizationofeachspectra. It isinteresting tonote that whileQ NM shows an increasing trend with increasing N +M for spectra with large scale power n = −2 in 2D, the trend is reversed for n = 0 and 2. The results presented in 2D are for N +M = 4,5 and 6 respectively. In 3D meaningful computation of Q were possible only for N +M < 4. We find that computed values of N,M Q are more stable towards fluctuation at large separationif one of the indices of cumulant correlatorsis NM larger than the other (e.g. w is more stable at larger separation as compared to w ). Close associations 31 22 of dashed curves with solid curves proves the validity of extended perturbation theory. It is also easy to notice that agreement is better in case of 2D compared to 3D. Larger simulation size might be a possible explanation for such an effect. Given that we restrict our result to highly nonlinear regime, it is intersting to note that computed values of C are not completely different from predictions of perturbation theory NM at large separation. In Figure - 5, dashed and solid straight lines are predicted values of Q2 and Q from 21 31 equation(17)andequation(18). ToppanelsinFigure-1andFigure-5teststhevalidityofequation(12)in 2Dand3Drespectively. Nearestneighborcurvescorrespondto differentvalues ofN andM whichproduces same N +M. Agreement in this case is similar to that of extended perturbation theory. Using cumulant correlators it is possible to compute amplitudes of different tree topologies contributing to four and five point correlationfunctions. We have computed amplitudes Q, R and R . Figure - 3 shows a b thattheseamplitudesarealmostconstantin2Dinthehighlynonlinearregime. In3DFigure-6showslarge fluctuationsincomputedvaluesofR althoughingeneraltheyarefairlyconstantwithinthelimitedrangeof b nonlinearitystudiedby us. Opentrianglesandsquaresdenote measuredvalues ofQ andQ usingfactorial 3 4 moments of CPDF. On the other hand solid triangles and squares denotes measured values of Q and Q 3 4 from cumulant correlators. Q was calculated by using the relation 16Q = 12R +4R , from computed 4 4 a b averagevalues of R andR . The remarkableagreementofboth the methods increasesconfidence in results a b basedoncumulantcorrelators. Itis alsoto be notedthatwhereascomputationofonepointcumulants were done by taking volume correctionsinto account, sucheffects were neglectedin the computationof cumulant correlators. OurresultsshowthatQ ismuchclosertoR whichmightindicatethatstartopologiesingeneral 4 b dominateoversnaketopologies. However,moresystematicstudiesarenecessarywhichwillincorporatefinite volume effects and have more dynamic range than this study. 6. Discussion Measurementsofcumulantcorrelatorsandlowerorderone pointcumulants havealreadybeendone for Lick catalog, SDES and the APM survey. Analyzing cumulant correlatorsin the APM survey, Szapudi & Szalay (1997)demonstratedthevalidityofthehierarchalansatztounprecedentedaccuracy. TheyfoundQ =1.15, 3 R = 5.3 and R = 1.15. These results were also shown to be in good agreement with computation of a b one point cumulants based on factorial moments Q = 1.7 and Q = 4.17. Measurements form SDES give 3 4 Q = 1.16 and Q = 2.96. It was suggested that slightly low value of Q form SDES was lack of nonlinear 3 4 4 correctionsin their computations. Other measurements fromAPMproduces Q =1.6Szapudi etal. (1996) 3 and Q = 1.7 by Gaztanaga (1994). The values of fourth order cumulants found by these authors were 3 Q =3.2 and Q =3.7 respectively. 4 4 Use of methods based on cumulant correlators by Szapudi and Szalay has produced lower values for R b comparedtoR . Inoursimulationsbothin2Dand3Dwefindthatrelativepositionofthesetwoamplitudes a dependoninitialpowerspectra. InallcasesR wasfoundtobe eitherequaltoorgreaterthanR ,although b a the separation was found to decrease with increasing n. For power spectra such as n = 1 in 3D and n = 2 in 2D which have more power at smaller scales these two quantities are found to almost coincide with each other. Although this is in disagreement with findings of Szapudi and Szalay (1997), our results seems to be closer to that of by Fry & Peebles (1978) who found R = 2.5±0.6 and R = 4.3±1.2 from their a b analysis of Lick catalog. At any rate, they are analyzing data in projection, and we a simulation with full information. We do not know to what extent projection and the possible inclusion of non–gravitational physics may influence the results. 10