Dark matter statistics for large galaxy catalogs: power spectra and covariance matrices

1.2 MB·

by  Anatoly Klypin

Mon.Not.R.Astron.Soc.000,000–000(0000) Printed23January2017 (MNLATEXstylefilev2.2)
Dark matter statistics for large galaxy catalogs: power
spectra and covariance matrices
Anatoly Klypin1(cid:63) and Francisco Prada2
1 Astronomy Department, New Mexico State University, Las Cruces, NM, USA
7
2 Instituto de Astrof´ısica de Andaluc´ıa (CSIC), Glorieta de la Astronom´ıa, E-18080 Granada, Spain
1
0
2

n
a 23January2017
J

0
2
ABSTRACT

Upcomingandexistinglarge-scalesurveysofgalaxiesrequireaccuratetheoreticalpre-
]
O dictions of the dark matter clustering statistics for thousands of mock galaxy cata-
logs. We demonstrate that this goal can be achieve with our new Parallel Particle-
C
Mesh (PM) N-body code (PPM-GLAM) at a very low computational cost. We run
h. about 15,000 simulations with ∼2billion particles that provide ∼1% accuracy of the
p dark matter power spectra P(k) for wave-numbers up to k ∼ 1hMpc−1. Using this
- large data-set we study the power spectrum covariance matrix, the stepping stone
o
for producing mock catalogs. In contrast to many previous analytical and numer-
r
t ical results, we find that the covariance matrix normalised to the power spectrum
s C(k,k(cid:48))/P(k)P(k(cid:48)) has a complex structure of non-diagonal components. It has an
a
[ upturn at small k, followed by a minimum at k ≈ 0.1 − 0.2hMpc−1. It also has
a maximum at k ≈ 0.5 − 0.6hMpc−1. The normalised covariance matrix strongly
1
evolves with redshift: C(k,k(cid:48)) ∝ δα(t)P(k)P(k(cid:48)), where δ is the linear growth factor
v
and α ≈ 1 − 1.25, which indicates that the covariance matrix depends on cosmo-
0
logical parameters. We also show that waves longer than 1h−1Gpc have very little
9
6 impact on the power spectrum and covariance matrix. This significantly reduces the
5 computationalcostsandcomplexityoftheoreticalpredictions:relativelysmallvolume
0 ∼(1h−1Gpc)3 simulations capture the necessary properties of dark matter clustering
. statistics.Allthepowerspectraobtainedfrommanythousandsofoursimulationsare
1
publicly available.
0
7
Key words: cosmology:Largescalestructure-darkmatter-galaxies:halos-meth-
1
ods: numerical
:
v
i
X
r
a 1 INTRODUCTION non-linear matter power spectrum P(k) (e.g., Smith et al.
2003; Heitmann et al. 2009; Schneider et al. 2016), the co-
Accuratetheoreticalpredictionsfortheclusteringproperties
variance matrix is a much less studied quantity that is re-
ofdifferentgalaxypopulationsarecrucialforthesuccessof
quiredfortheanalysisoftheobservationaldata.Mostofthe
massiveobservationalsurveyssuchasSDSS-III/BOSS(e.g.,
timethemainattentionshiftstothelaststep,i.e.,theplace-
Alam et al. 2016), SDSS-IV/eBOSS (Dawson et al. 2016),
ment of “galaxies” in the density field with all the defects
DESI(DESICollaborationetal.2016)andEuclid(Laureijs
and uncertainties of estimates of the dark matter distribu-
et al. 2011) that have or will be able to measure the posi-
tionandvelocitiesbeingoftenincorporatedintothebiasing
tionsofmillionsofgalaxies.Thepredictionsofthestatistical
scheme (e.g., Chuang et al. 2015). However, one needs to
propertiesofthedistributionofgalaxiesexpectedinmodern
accurately reproduce and to demonstrate that whatever al-
cosmologicalmodelsinvolveanumberofsteps:itstartswith
gorithm or prescription is adopted to mock galaxies, also is
the generation of the dark matter clustering and proceeds
able to match the clustering and covariance matrix of the
withplacing“galaxies”accordingtosomebiasprescriptions
dark matter field.
(see Chuang et al. 2015; Monaco 2016, for a review).
Whileasubstantialefforthasbeenmadetoestimatethe The covariance matrix C(k,k(cid:48)) is the second-order
statistics of the power spectrum: C(k,k(cid:48)) = (cid:104)P(k)P(k(cid:48))(cid:105)−
(cid:104)P(k)(cid:105)(cid:104)P(k(cid:48))(cid:105),where(cid:104)...(cid:105)impliesaveragingoveranensam-
(cid:63) E-mail:[email protected] ble of realizations. The power spectrum covariance and its
(cid:13)c 0000RAS
2 Klypin & Prada
cousinthecovarianceofthecorrelationfunctionplayanim- thecomputationalvolumesshouldbelargeenoughtocover
portantroleinestimatesoftheaccuracyofmeasuredpower theentireobservationalsample.Indeed,thiswillbetheideal
spectrum,andtheinversecovariancematricesareusedines- case,ifweweretomeasuresomestatisticsthatinvolvewaves
timatesofcosmologicalparametersdeducedfromthesemea- as long as the whole sample. The problem with producing
surements (e.g., Anderson et al. 2012; Sa´nchez et al. 2012; extra large simulations is their computational cost. For ac-
Dodelson&Schneider2013;Percivaletal.2014).Thepower curate predictions one needs to maintain the resolution on
spectrum covariance matrix measures the degree of non- small scales. With a fixed resolution the cost of computa-
linearity and mode coupling of waves with different wave- tions scales with volume. So, it increases very quickly and
numbers. As such it is an interesting entity on its own. becomes prohibitively expensive.
Because it is expensive to produce thousands of sim- However, most of the observable clustering statistics
ulations (e.g., Taylor et al. 2013; Percival et al. 2014) us- rely on relatively small scales, and thus would not require
ingstandardhigh-resolutionN-bodycodes,therewererela- extra large simulation boxes. For example, the Baryonic
tivelyfewpublicationsonthestructureandevolutionofthe AcousticOscillations(BAO)areonscalesof∼100h−1Mpc
covariance matrix. Most of the results are based on simu- and the power spectrum of galaxy clustering is typically
lations with relatively small number of particles (16 million measured on wave-numbers k>0.05hMpc−1. In order to
∼
ascomparedwithour∼2billion)andsmallcomputational make accurate theoretical predictions for these scales, we
volumesof500−600h−1Mpc(Takahashietal.2009;Lietal. need relatively small ∼1h−1Gpc simulation volumes.
2014;Blotetal.2015).Therewasnosystematicanalysisof Whatwillbetheconsequencesofusingsmallcomputa-
the effects of mass, force resolution and volume on the co- tional volumes? One may think about few. For example,
variance matrix. long-waves missed in small simulation boxes may couple
So far, there are some uncertainties and disagreements with small-scale waves and produce larger power spectrum
evenontheshapeofthecovariancematrix.Neyrinck(2011); (and potentially covariance matrix). This is the so called
Mohammed&Seljak(2014);Carronetal.(2015)arguethat SuperSampleCovariance(SSC,Gnedinetal.2011;Lietal.
once the Gaussian diagonal terms are removed, the covari- 2014; Wagner et al. 2015; Baldauf et al. 2016). Another ef-
ance function should be a constant. This disagrees with fect is the statistics of waves. By replicating and stacking
numerical (Li et al. 2014; Blot et al. 2015) and analytical small boxes (to cover large observational volumes) we do
(Bertolini et al. 2016; Mohammed et al. 2016) results that not add the statistics of small waves, which we are inter-
indicatethatthenon-diagonalcomponentsofthecovariance ested in. For example, the number of pairs with separation
matrix have quite complex structure. of say 100h−1Mpc will be defined only by how many inde-
Therearealsoissuesandquestionsrelatedwithnumer- pendentpairsareinasmallbox,andnotbythemuchlarger
ical simulations. How many realizations are needed for ac- numberofpairsfoundinlargeobservationalsamples.These
curatemeasurementsofthecovariancematrix(Tayloretal. concernsarevalid,butcanberesolvedinanumberofways.
2013;Percivaletal.2014)andhowcanitbereduced(Pope SSCeffectsdependonthevolumeandbecomeverysmallas
&Szapudi2008;Pearson&Samushia2016;Joachimi2017)? thecomputationalvolumeincreases.Statisticsofwavescan
Whatresolutionandhowmanytime-stepsareneededforac- bere-scaledproportionallytothevolume(e.g.,Mohammed
curate estimates? How large the simulation volume should & Seljak 2014; Bertolini et al. 2016). We investigate these
be to adequately probe relevant scales and to avoid defects issuesindetailinourpaper,togetherwiththeperformance
relatedwiththeboxsize(Gnedinetal.2011;Lietal.2014)? ofthepowerspectrumdependingonthenumericalparame-
The goal of our paper is to systematically study the ters of our PPM-GLAM simulations. All the power spectra
structure of the matter covariance matrix using a large set obtainedfrommanythousandsofoursimulationsaremade
of N-body simulations. The development of N-body cos- publicly available.
mological algorithms has been an active field of research In Section 2 we discuss the main features of our PPM-
for decades. However, the requirements for the generation GLAMsimulationcode.Moredetaileddescriptionandtests
of many thousandsof high-quality simulations areextreme. are presented in the Appendixes. The suite of simulations
Existing codes such as GADGET, RAMSES, or ART are usedinthispaperispresentedinSection3.Convergenceand
powerful for large high-resolution simulations, but they are accuracy of the power spectrum are discussed in Section 4.
notfastenoughformediumqualitylarge-numberofrealisa- Theresultsonthecovariancematrixofthepowerspectrum
tionsrequiredforanalysisandinterpretationoflargegalaxy are given in Section 5. We campare our results with other
surveys. New types of codes (e.g., White et al. 2014; Tas- works in Section 6. Summary and discussion of our results
sev et al. 2013; Feng et al. 2016) are being developed for are presented in Section 7.
this purpose. Here, we present the performance results of
ournewN-bodyParallelParticle-MeshGLAMcode(PPM-
GLAM), which is the core of the GLAM (GaLAxy Mocks)
2 PARALLEL PARTICLE-MESH GLAM CODE
pipelineforthemassiveproductionoflargegalaxycatalogs.
PPM-GLAM generates the density field, including peculiar There are a number of advantages of cosmological Particle-
velocities, for a particular cosmological model and initial Mesh (PM) codes (Klypin & Shandarin 1983; Hockney &
conditions. Eastwood 1988; Klypin & Holtzman 1997) that make them
The optimal box size of the simulations is of particu- useful on their own to generating a large number of galaxy
larinterestforproducingthousandsofmockgalaxycatalogs mocks(e.g.,QPM,COLA,FastPM;Whiteetal.2014;Tas-
andforthestudiesoflarge-scalegalaxyclustering.Withthe sevetal.2013;Fengetal.2016),orasacomponentofmore
upcomingobservationalsamplestendingtocoverlargerob- complexhybridTREE-PM(e.g.,Gadget2,HACC;Springel
servational volumes (∼ 50Gpc3), one naively expects that 2005; Habib et al. 2014) and Adaptive-Mesh-Refinement
(cid:13)c 0000RAS,MNRAS000,000–000
Dark matter statistics for large galaxy catalogs 3
Table 1. Numerical and cosmological parameters of different simulations. The columns give the simulation identifier, the
size of the simulated box in h−1Mpc, the number of particles, the mass per simulation particle mp in units h−1M(cid:12), the
mesh size Ng3, the gravitational softening length (cid:15) in units of h−1Mpc, the number of time-steps Ns, the amplitude of
perturbationsσ8,thematterdensityΩm,thenumberofrealisationsNr andthetotalvolumein[h−1Gpc]3
Simulation Box particles mp Ng3 (cid:15) Ns σ8 Ωm Nr TotalVolume
PPM-GLAMA0.5 5003 12003 6.16×109 24003 0.208 181 0.822 0.307 680 85
PPM-GLAMA0.9 9603 12003 4.46×1010 24003 0.400 136 0.822 0.307 2532 2240
PPM-GLAMA1.5 15003 12003 1.66×1011 24003 0.625 136 0.822 0.307 4513 15230
PPM-GLAMA2.5 25003 10003 1.33×1012 20003 1.250 136 0.822 0.307 1960 30620
PPM-GLAMA4.0 40003 10003 5.45×1012 20003 1.250 136 0.822 0.307 4575 292800
PPM-GLAMB1.0a 10003 16003 2.08×1010 32003 0.312 147 0.828 0.307 10 10
PPM-GLAMB1.0b 10003 13003 1.78×1010 26003 0.385 131 0.828 0.307 10 10
PPM-GLAMB1.5 15003 13003 3.88×1010 26003 0.577 131 0.828 0.307 10 33
PPM-GLAMC1a 10003 10003 8.71×1010 30003 0.333 302 0.828 0.307 1 1
PPM-GLAMC1b 10003 10003 8.71×1010 40003 0.250 136 0.828 0.307 1 1
BigMDPL1 25003 38403 2.4×1010 – 0.010 – 0.828 0.307 1 15.6
HMDPL1 40003 40963 7.9×1010 – 0.025 – 0.828 0.307 1 64
Takahashietal.2 10003 2563 1.8×1012 – – – 0.760 0.238 5000 5000
Lietal.3 5003 2563 2.7×1011 – – – 0.820 0.286 3584 448
Blotetal.4 6563 2563 1.2×1012 – – – 0.801 0.257 12288 3469
BOSSQPM5 25603 12803 3.0×1011 12803 2.00 7 0.800 0.29 1000 16700
WiggleZCOLA6 6003 12963 7.5×109 38883 – 10 0.812 0.273 3600 778
References:1Klypinetal.(2016),2Takahashietal.(2009),3Lietal.(2014),4Blotetal.(2015),
5Whiteetal.(2014),6Kodaetal.(2016)
codes (e.g. ART, RAMSES, ENZO; Kravtsov et al. 1997; sity field used in the Poisson equation is obtained with the
Teyssier 2002; Bryan et al. 2014). Cosmological PM codes Cloud-In-Cell(CIC)schemeusingthepositionsofdarkmat-
are the fastest codes available and they are simple. terparticles.Oncethegravitationalpotentialisobtained,it
We have developed and thoroughly tested a new par- isnumericallydifferentiatedandinterpolatedtotheposition
allel version of the Particle-Mesh cosmological code, that of each particle. Then, particle positions and velocities are
provides us with a tool to quickly generate a large num- advanced in time using the second order leap-frog scheme.
ber of N-body cosmological simulations with a reasonable The time-step is increased periodically as discussed in Ap-
speed and acceptable resolution. We call our code Parallel pendix A. Thus, a standard PM code has three steps that
Particle-Mesh GLAM (PPM-GLAM), which is the core of are repeated many times until the system reached its final
the GLAM (GaLAxy Mocks) pipeline for massive produc- moment of evolution: (1) Obtain the density field on a 3D-
tion of galaxy catalogs. PPM-GLAM generates the density mesh that covers the computational volume, (2) Solve the
field, including peculiar velocities, for a particular cosmo- Poissonequation,and(3)Advanceparticlestothenextmo-
logical model and initial conditions. Appendix A gives the ment of time.
detailsofthecodeandprovidestestsfortheeffectsofmass ThecomputationalcostofasinglePPM-GLAMsimula-
andforceresolutions,andtheeffectsoftime-stepping.Here, tiondependsonthenumberoftime-stepsN ,thesizeofthe
s
we discuss the main features of the PPM-GLAM code and 3D-meshN3,andtheadoptednumberofparticlesN3.The
g p
provide the motivation for the selection of appropriate nu- CPUrequiredtosolvethePoissonequationismostlydeter-
merical parameters. mined by the cost of performing a single 1D-FFT. We in-
The code uses a regularly spaced three-dimensional corporateallnumericalfactorsintoonecoefficientandwrite
mesh of size N3 that covers the cubic domain L3 of a sim- the CPU for the Poisson solver as AN3. The costs of den-
g g
ulation box. The size of a cell ∆x = L/N and the mass sityassignmentandparticledisplacement(includingpoten-
g
of each particle m define the force and mass resolution re- tialdifferentiation)scaleproportionallytoN3.Intotal,the
p p
spectively: CPU time T required for a single PPM-GLAM run is:
tot
m = Ω ρ (cid:20) L (cid:21)3 = (1) Ttot =Ns(cid:2)ANg3+(B+C)Np3(cid:3), (4)
p m cr,0 N
p where B and C are the coefficients for scaling the CPU es-
= 8.517×1010(cid:20)Ωm (cid:21)(cid:20)L/h−1Gpc(cid:21)3h−1M , (2) timate for particle displacements and density assignment.
0.30 N /1000 (cid:12) Thesenumericalfactorswereestimatedfordifferentproces-
p
(cid:20)L/h−1Gpc(cid:21) sors currently used for N-body simulations and are given
∆x = h−1Mpc, (3) in Table A1. For a typical simulation analysed in this pa-
N /1000
g
per (N = 2400, N = N /2) the CPU per time-step is
g p g
where N3 is the number of particles and ρ is the critical ∼0.5hoursandwall-clocktimeperstep∼1−3minutes.The
p cr,0
density of the universe at present. total cost of 1000 PPM-GLAM realizations with N = 150
s
PPM-GLAM solves the Poisson equation for the gravi- is 75K CPU hours, which is a modest allocation even for a
tational potential in a periodical cube using a Fast Fourier small computational cluster or a supercomputer center.
Transformation (FFT) algorithm. The dark matter den- Memory is another critical factor that should be con-
(cid:13)c 0000RAS,MNRAS000,000–000
4 Klypin & Prada
sidered when selecting the parameters of our simulations. this is the largest set of simulations available today. Power
PPM-GLAM uses only one 3D-mesh for storing both den- spectraandcovariancematricesarepubliclyavailableinour
sityandgravitationalpotential,andonlyonesetofparticle Skies and Universes site1.
coordinates and velocities. Thus, for single precision vari- The PPM-GLAM simulations labeled with letter A
ables the total required memory M is: are the main simulations used for estimates of the power
tot
spectrum and the covariance matrix. Series B and C are
designed to study different numerical effects. In particu-
M = 4N3+24N3 Bytes, (5)
tot g p lar, C1a are actually four simulations run with the differ-
(cid:18) N (cid:19)3 (cid:18) N (cid:19)3 ent number of steps: N = 34,68,147,302. There are also
= 29.8 g +22.3 p GB, (6) s
2000 1000 four C1b simulations that differ by the force resolution:
(cid:18) N (cid:19)3 Ng =1000,2000,3000,4000.
= 52 p GB, for N =2N . (7) WecomparetheperformanceofourPPM-GLAMsimu-
1000 g p
lationswiththeresultsobtainedfromsomeoftheMultiDark
Thenumberoftime-stepsN isproportionaltothecom- simulations2runwithL-Gadget2:BigMDPLandHMDMPL
s
putational costs of the simulations. This is why reducing (seefordetailsKlypinetal.2016).Theparametersofthese
the number of steps is important for producing a large set large-boxandhigh-resolutionsimulationsarealsopresented
of realisations. White et al. (2014) and Koda et al. (2016) in Table 1. For comparison we also list very large number
use just ∼ 10 time-steps for their QPM and COLA simu- of low resolution simulations performed by Takahashi et al.
lations. Feng et al. (2016) and Izard et al. (2015) advocate (2009),Lietal.(2014)andBlotetal.(2015)withtheGad-
usingN ≈40stepsforFast-PMandICE-COLA.Theques- get2,L-Gadget2andAMRcodes,respectively,tostudythe
s
tionstillremains:whatoptimalnumberoftime-stepsshould power spectrum covariances. Details of the QPM (White
be adopted? However, there is no answer to this question et al. 2014) and COLA (Koda et al. 2016) simulations that
without specifying the required force resolution, and with- were used to generate a large number of galaxy mocks for
out specifying how the simulations will be used to generate theBOSSandWiggleZgalaxyredshiftsurveysarealsogiven
mock galaxies. in Table 1. Note that the QPM simulations have very low
In Appendix B we provide a detailed discussion on the forceresolution,whichrequiressubstantialmodelingonhow
effectsoftime-stepping.Wearguethatforthestabilityand darkmattershouldbeclusteredandmovingonthescaleof
accuracy of the integration of the dark matter particle tra- galaxies.
jectoriesinsidedense(quasi-)virialisedobjects,suchasclus- WeestimatethepowerspectrumP(k)ofthedarkmat-
tersofgalaxies,thetime-step∆tmustbesmallerenoughto ter density field in all our 10 simulation snapshots for each
satisfy the constraints given by eqs. (A20) and (A22). For realisation, but in this paper we mostly focus on the z =0
example, FastPM simulations with 40 time-steps and force results. For each simulation we estimate the density on a
resolutionof∆x=0.2h−1Mpc(seeFengetal.2016)donot 3D-mesh of the size indicated in Table 1. We then ap-
satisfy these conditions and would require 2-2.5 times more ply FFT to generate the amplitudes of the Fourier har-
time-steps. However, a small number of time-steps mani- monics δ in phase-space. The spacing of the Fourier
i,j,k
fests itself not in the power spectrum (though, some de- harmonics is equal to the length of the fundamental har-
clineinP(k)happensatk∼1h−1Mpc).Itseffectismostly monic∆κ≡2π/L.Thus,thewave-vectork correspond-
i,j,k
observed in a significantly reduced fraction of volume with ing to each triplet (i,j,k), where i,j,k = 0,...N −1, is
g
largeoverdensitiesandrandomvelocities,whichpotentially k = (i∆κ,j∆κ,k∆κ). Just as the spacing ∆x = L/N
i,j,k g
introduces undesirable scale-dependent bias. in real-space represents the minimum resolved scale (see
Because our main goal is to produce simulations with Sec.2),thefundamentalharmonic∆κistheminimumspac-
the minimum corrections to the local density and peculiar ing in Fourier-space, i.e. one cannot probe the power spec-
velocities, we use N ≈ 100−200 time-steps in our PPM- trum below that scale. To estimate the power spectrum we
s
GLAM simulations. This number of steps also removes the use a constant bin size equal to ∆κ. This binning results
need to split particle displacements into quasi-linear ones in very fine binning at high frequencies, but preserves the
and the deviations from quasi-linear predictions. Thus, in phase resolution at very small frequencies (long waves).
this way we greatly reduce the complexity of the code and A correction is applied to the power spectrum to com-
increaseitsspeed,whilealsosubstantiallyreducethemem- pensate the effects of the CIC density assignment: P(k) =
ory requirements. P (k)/(cid:2)1−(2/3)sin2(πk/2k )(cid:3), where the Nyquist fre-
raw Ny
quency of the grid is k = (N /2)∆κ = π/∆x. The same
Ny g
number of grid points is used for estimates of the power
spectrum as for running the simulations. We typically use
3 SIMULATIONS
resultsonlyfork<(0.3−0.5)k .Nocorrectionsareapplied
Ny
We made a large number of PPM-GLAM simulations – forthefinitenumberofparticlesbecausethesearesmallfor
about 15,000 – to study different aspects of the clustering the scales and particle number-densities considered in this
statisticsofthedarkmatterdensityfieldintheflatΛCDM paper.
Planck cosmology. The numerical parameters of our simu- Similar to CIC in real space, we apply CIC filtering
lations are presented in Table 1. All the simulations were in Fourier space. For each Fourier harmonic ki,j,k the code
started at initial redshift z = 100 using the Zeldovich
init
approximation.Thesimulationsspanthreeordersofmagni-
tudeinmassresolution,afactorofsixinforceresolutionand 1 http://projects.ift.uam-csic.es/skies-universes/
differbyafactorof500ineffectivevolume.Toourknowledge 2 http://www.multidark.org
(cid:13)c 0000RAS,MNRAS000,000–000
Dark matter statistics for large galaxy catalogs 5
Figure 1.Powerspectraofdarkmatteratredshiftzero.Thelinearpowerspectrumisshownasadottedline.Left:ThePPM-GLAM
simulations used for the plot are A0.5 (red full curve), A1.5 (black dot-dashed) and A2.5 (blue dashed). They closely reproduce the
clusteringofthehigh-resolutionMultiDarksimulations(blackfullcurve)uptok≈1hMpc−1withexactdeviationsonlargerkdepending
on the force resolution. Right: Zoom-in on the region of the BAO peaks. The power spectra were multiplied by k5/4 and arbitrarily
normalised to reveal more clearly the differences between the simulations. Because there are only two realisations of the MultiDark
simulations(BigMDPL&HMDPL),statisticaldeviationsduetocosmicvarianceareseenatdifferentk (e.g.,k≈0.045,0.07hMpc−1).
Cosmic variance of the simulations (errors of the mean) are nearly negligible because there are thousands of PPM-GLAM realisations.
ThedeviationsseenfortheA0.5points(redcircles)atsmallk areduetothelargebinwidth,ink-space,forthissimulation.
finds left and right bins by dividing the length of the wave- PPM-GLAM simulations. As a reference we also show the
vectorbythefundamentalharmonic∆κandthenbytaking linear power spectrum (dotted line).
the integer part INT(|k |/∆κ). The contributions to the
i,j,k
leftandrightbinsarefoundproportionallytothedifference Oursimulationscloselyreproducewelltheclusteringof
between harmonic and bin wave-numbers. This procedure thehigh-resolutionMultiDarksimulations(blackfullcurve)
reduces the noise in the power spectrum by ∼ 30% at the both for long-waves (as may have been expected) and even
cost of introducing dependencies in power spectrum values for larger wave-numbers up to k ≈ 1hMpc−1 with exact
inadjacentbinsinFourierspace.Effectsofthisfilteringare deviations on larger k depending on force resolution. The
included in estimates of the covariance matrix: they mostly lack of force resolution results in the decline of the power
change (reduce) diagonal components. spectrumatlargewave-numbers:asresolutionincreasesthe
power spectrum becomes closer to the MultiDark results.
This is clearly seen in the left panel of Figure 2 where we
show the ratios of the PPM power spectra P(k) of A0.5
(reddot-dashedcurve),A1.5(bluefullcurve)andA2.5(red
4 POWER SPECTRUM
dashedcurve)tothepowerspectrumobtainedfromtheMul-
Westudytheconvergenceperformanceofthepowerspectra tiDarksimulations.Wealsolabelintheplottheforceresolu-
obtained with PPM-GLAM by comparing our results with tionforeachofthePPM-GLAMsimulations.Theresultsfor
thosedrawnfromthehigh-resolutionMultiDarksimulations theA0.5simulationsarepresentedonlyfork>0.15hMpc−1
listed in Table 1. Specifically, we average power spectra of because the bin smearing becomes visible (∼ 2%) at lower
BigMDPLandHMDPLsimulationsweightedproportianally frequencies due to the small volume of each individual re-
to the volume of each simulation. Left panel in Figure 1 alisation. We also plot the average of the B1.0a and B1.0b
showsthepowerspectrainlog−log scalefortheA0.5(red simulations(blacklong-dashedcurve).Again,thedeviations
fullcurve),A1.5(blackdot-dashed)andA2.5(bluedashed) are less than 1% on large scales and they start to increase
(cid:13)c 0000RAS,MNRAS000,000–000
6 Klypin & Prada
Figure 2. Convergence of power spectra in real-space at redshift zero. Left: Ratios of the real-space power spectra in PPM-GLAM
simulationsto the power spectrum of the MultiDarksimulations. Simulations used for the plot areindicated in the plot.Forthe black
long dashed curve we use the average of B1.0a and B1.0b simulations. The lack of force resolution results in the decline of the power
spectrum at large wave-numbers. As the resolution increases, the power spectrum becomes closer to the MultiDark results. At low
valuesofk<0.3hMpc−1 deviationsat∼0.5%levelarecausedbycosmicvarianceintheMultiDarksimulations.Otherwisethereareno
∼
systematicsrelatedwiththefiniteboxsize.Right:ConvergenceofpowerspectrainPPM-GLAMsimulationswithdifferentboxsizesand
resolutions.WeplottheratiosofpowerspectraofA0.9,A1.5,A2.5andA4.0simulationstothecombinedpowerspectrumPPM thatis
foundbyaveragingbestsimulationsfordifferentrangesofwave-numbers(seetext).Heretheeffectsofcosmicvariance(seenontheleft
panel)arenegligiblebecauseoftheaveragingoverthousandsofrealisations.Missingwaveslongerthanthesimulationboxeshavelittle
effectforL>1h−1Mpcandk>0.05hMpc−1.
∼
as we go to larger k with the magnitude of the error de- This is clearly seen in the left panel of Figure 2 were
pending on the force resolution. Note that the ratios of the we show the ratios of GLAM power spectra to P(k) in the
PPM-GLAM results to those in the MultiDark simulations MultiDarksimulations.ResultsfortheA0.5simulationsare
are the same at long-waves with k < 0.1hMpc−1. This is presentedonlyfork>0.15hMpc−1 becausethebinsmear-
related with the cosmic variance present in the MultiDark ingbecomesvisible(∼2%)atlowerfrequencies.Again,the
P(k)sincethereareonlytworealisations,i.e.BigMDPLand deviations are less than 1% on large scales and start to in-
HMDPL. crease as we go to large k with the magnitude of the error
The right panel in Figure 1 zooms-in on the relevant depending on the force resolution. Note that the ratios of
domaink≈0.07−0.2hMpc−1oftheBAOpeaks.Inthisplot the PM results to those in the MultiDark simulations are
thepowerspectrumP(k)ismultipliedbythefactork5/4 to the same at k < 0.1hMpc−1. This is related with the cos-
reducethedynamicalscaleallowingustoseethedifferences mic variance in the MultiDark P(k) – there are only two
assmallasafractionofpercent.Thecosmicvarianceofthe realizations of MultiDark.
PPM-GLAM simulations (errors of the mean) are nearly In order to test the effects of force resolution (cid:15) and fi-
negligible because there are thousands of realisations. The nite box size L we construct the combined power spectrum
observed deviations of the A0.5 points at small k (e.g. k ≈ by taking the average of the best PPM-GLAM simulations
0.08hMpc−1) are due to the large size of the binning in k in separate ranges of frequency: (1) for k < 0.1hMpc−1
space defined by the width of the fundamental harmonic we average P(k) of all A1.5, A2.5, and A4.0 realisations,
∆k = 2π/L = 0.0125hMpc−1. If we consider simulations (2) for the range 0.1hMpc−1 < k < 0.2hMpc−1 we take
with small binning, i.e. simulation box L>1h−1Gpc, then the average of the A2.5, A1.5, A0.9 simulations, (3) for
∼
thedeviationsfromtheMultiDarksimulationsarelessthan 0.2hMpc−1 < k < 0.4hMpc−1, A0.9 and A0.5 simulations
1 per cent on the large scales. are used, and (4) for larger wave-numbers we consider only
(cid:13)c 0000RAS,MNRAS000,000–000
Dark matter statistics for large galaxy catalogs 7
Figure3.Comparisonofpowerspectrainredshift-spaceatredshiftzero.Monopole(leftpanel)andquadrupole(rightpanel)werescaled
with different power of k to reduce the dynamical range of the plots. The solid black curves correspond to the BigMDPL simulation
estimatedusing5%ofallparticles.ThedashedcurvesanderrorbarsarefortheB1.5simulations.
the A0.5 realisations. We show in Figure 2 (right panel) withhighresolution,missingwaveslongerthanthesescales
theratiosofP(k)ofA0.9,A1.5,A2.5andA4.0simulations are deeply in the linear regime, and thus the SSC effects
to the combined power spectrum. The deviations of each are expected to be small, which the right panel in Figure 2
simulation from this combined P (k) spectrum is a mea- clearly demonstrates.
PM
sureoftheerrorsineachsimulation.Notethattheseerrors
are broadly consistent with those of MultiDark except for The SSC effects should manifest themselves as an in-
the very long waves where we now do not have artificial creaseinP(k)atsmallkinsimulationswithlargeLascom-
deviations due to the cosmic variance. This plot gives us pared with simulations with smaller L. Indeed, we see this
means to estimate what resolution is needed if we are re- effect at very long-waves. For example, at k=0.03hMpc−1
quired to provide some specified accuracy at a given scale. the power spectrum in A0.5 simulations was below that of
For example, if the errors in P(k) should be less than 1% A4.0 by 4%. However, the effect becomes much smaller as
at k < 0.5hMpc−1, then the resolution of the simulation the box size increases. The error becomes less than 0.2%
shouldbe∆x=0.62h−1Mpc.For,1%atk<1hMpc−1 the for the A1.5 simulations. It is also much smaller for shorter
resolution should be ∆x=0.2h−1Mpc. waves. For example, the error is less than 0.5% for A0.9
simulations at k >0.05hMpc−1. This can be understood if
TherightpanelinFigure2alsogivesveryusefulinfor- oneestimatestheamplitudeofdensityfluctuationsinwaves
mation on the effects due to the finite box size. The size of longerthanthesimulationbox.ForL=500h−1Mpcandfor
the computational box is an important factor that affects thePlanckcosmologythermsdensityfluctuationδ(>L)is
the total CPU time, the statistics of the large-scale fluc- relatively large: δ(> 500h−1Mpc) = 0.027. It is nearly ten
tuations, and possibly the non-linear coupling of long- and times smaller for A1.5 simulations: δ(> 1500h−1Mpc) =
short waves. Non-linear coupling effects are of some con- 0.0036.
cern (e.g., Gnedin et al. 2011; Li et al. 2014; Wagner et al.
2015; Baldauf et al. 2016) because the long-waves missed While the main interest of this paper is in the cluster-
in small-box simulations (called Super Sample Covariance ing in real-space, we also tested dark matter clustering in
or SSC) can affect both the power spectrum and the co- redshift-space and compared that with the MultiDark sim-
variance matrix. The magnitude of the SSC effect critically ulations. We plot in Figure 3 both monopole (left panel)
depends on the size of the computational box. Because our andquadrupole(rightpanel)forPPM-GLAMB1.5andBig-
main target is relatively large boxes L ≈ (1−1.5)h−1Mpc MDPL, which shows a remarkable agreement.
(cid:13)c 0000RAS,MNRAS000,000–000
8 Klypin & Prada
4.0 1.0
1.0 3.6 1.0 0.9
3.2 0.8
0.8 0.8
1k (hMpc)′−00..46 11222.....26048log(Cov[k,k])′ 1k (hMpc)′−00..46 00000.....345671/2Cov[k,k]/(Cov[k,k]Cov[k,k])′′′
0.2 0.2
0.8 0.2
0.4 0.1
0.2 0.4 0.6 0.8 1.0 0.2 0.4 0.6 0.8 1.0
k (hMpc−1) k (hMpc−1)
0.0 0.0
(cid:112)
Figure4.CovariancematrixC(k,k(cid:48))onlogarithmicscale(leftpanel)andcovariancecoefficientC(k,k(cid:48))/ C(k,k)C(k(cid:48),k(cid:48))(rightpanel)
ofthedarkmatterpowerspectrumatredshiftzeroforthePPM-GLAMA1.5simulations.Withtheexceptionofthenarrowspikeatthe
diagonalofthematrix,thecovariancematrixisasmoothfunction.Horizontalandverticalstripesseeninthecorrelationcoefficientat
smallk areduetocosmicvariance.
5 COVARIANCE MATRIX. where the numerical factor α is equal to unity for the
Nearest-Grid-Point(NGP)assignmentinFourier-spaceand
ThecovariancematrixC(k,k(cid:48))isdefinedasareducedcross
α = 2/3 for the CIC assignment used in this paper. Note
product of the power spectra at different wave-numbers for
that for a fixed bin width ∆k the number of harmon-
thesamerealisationaveragedoverdifferentrealisations,i.e.
ics, and thus, the amplitude of the Gaussian term scales
C(k,k(cid:48))≡(cid:104)P(k)P(k(cid:48))(cid:105)−(cid:104)P(k)(cid:105)(cid:104)P(k(cid:48))(cid:105). (8) proportional to the computational volume Nh ∝ L3 with
C (k)∝1/L3.
Gauss
The diagonal and non-diagonal components of the covari- There are two ways of displaying the covariance ma-
ance matrix have typically very different magnitudes and trix. One can normalise C(k,k(cid:48)) by its diagonal compo-
evolve differently with redshift. Their diagonal terms are nent:r(k,k(cid:48))≡C(k,k(cid:48))/(cid:112)C(k,k)C(k(cid:48),k(cid:48)).Thisquantityis
larger than the off-diagonal ones, but there are many more calledthecorrelationcoefficient,andbydefinition,r(k,k)≡
off-diagonal terms making them cumulatively important 1.Thecovariancematrixcanalsobenormalisedbythe”sig-
(Taylor et al. 2013; Percival et al. 2014; O’Connell et al. nal”, i.e. the product of power spectra at the two involved
(cid:112)
2016). Off-diagonal terms are solely due to non-linear clus- wave-numbers: C(k,k(cid:48))/P(k)P(k(cid:48)).
teringeffects:inastatisticalsensetheoff-diagonaltermsare Figure 4 shows the covariance matrix C(k,k(cid:48)) and the
equal to zero in the linear regime. The diagonal component correlationcoefficientforthePPM-GLAMA1.5simulations
C(k,k)canbewrittenasasumofthegaussianfluctuations at z=0. With the exception of a narrow spike at the diag-
due to the finite number of harmonics in a bin and terms onal of the matrix, the covariance matrix is a smooth func-
that are due to non-linear growth of fluctuations: tion. Horizontal and vertical stripes seen in the correlation
coefficientatsmallkareduetocosmicvariance.Theygrad-
C(k,k)≡C (k)+C (k,k), (9)
Gauss non uallybecomeweakerasthenumberofrealisationsincreases
(e.g., Blot et al. 2015).
where the Gaussian term depend on the amplitude of the
power spectrum P(k) and on the number of harmonics N : The diagonal terms of the covariance matrix are pre-
h
sented in Figure 5. In the left panel we compare the results
2 4πk2∆k of various simulations at z = 0 with different box sizes. In
C (k)=α P2(k), N = , (10)
Gauss Nh h (2π/L)3 ordertodothat,werescaletheindividualC(k,k)tothatof
(cid:13)c 0000RAS,MNRAS000,000–000
Dark matter statistics for large galaxy catalogs 9
Figure 5.DiagonalcomponentsofthecovariancematrixforPPM-GLAMsimulationswithdifferentboxsizesatz=0(leftpanel),and
fordifferentredshiftsfortheA0.9simulations(rightpanel).Allresultswererescaledtothesimulationvolumeof1.5h−1Gpc.Thedotted
linesshowtheGaussiancontributionCGauss givenbyeq.(10)withα=2/3,whichgivesagoodapproximationuptok∼<0.2hMpc−1.
At larger wave-numbers the covariance matrix is substantially larger than the Gaussian term due to the non-linear coupling of waves.
Thediagonaltermsevolvewithredshiftandsensitivelydependontheforceresolution.
the volume of the A1.5 simulation with (1.5h−1Gpc)3. Up reasonswhyoneneedsthousandsofrealisationstoestimate
to k<0.2hMpc−1 the covariance matrix is well described the covariance matrix reliably.
∼
by the Gaussian term, but at larger wave-numbers it sub- However, with the exception of the diagonal compo-
stantiallyexceedstheGaussiancontributionduetothenon- nentes,thecovariancematrixisasmoothfunctionthatcan
linear coupling of waves. The force resolution plays an im- be approximated analytically using a relatively small num-
portant role here. ber of parameters. The following approximation gives 3%
accuracyforthez=0covariancematrixatk>0.1hMpc−1
∼
TherearenoindicationsthatSSCwaves(modeslonger and ∼10% at smaller wave-numbers:
than the simulation box) affect the diagonal components.
C(k,k(cid:48)) = 8.64×10−6L−3P(k)P(k(cid:48))
If present, SSC effects should results in enhanced C(k,k)
in simulations with very large simulation boxes (Li et al. × (cid:2)1+g(k)+g(k(cid:48))(cid:3)2 (11)
2014;Baldaufetal.2016).Forexample,A4.0resultsshould (cid:34) (cid:35)2
have a larger covariance matrix as compared with that of × 1−αe−(k−2σk2(cid:48))2 +(1+α)e−2(k(0−.1kσ(cid:48)))22 ,
the A0.9 simulations. However, the right panel of Figure 5
clearlyshowstheoppositeeffect:A0.9resultsareaboveA4.0
presumablyduetothebetterforceresolutionthatproduces where
larger non-linear effects.
(cid:34) (cid:18) k (cid:19)2(cid:35) 1.1k0.60
g(k)=1.4exp − + , (12)
Figures 6 and 7 demonstrate that the non-diagonal 0.07 1+1.2k3
terms of the covariance matrix have a complex structure.
Theydependonbothkandk(cid:48) inanon-trivialwayandthey and
evolve with redshift in a complicated fashion. The bottom α=0.03/k(cid:48), σ=4κ=4(2×10−3π/L). (13)
leftpanelinFigure6highlightsoneofthemainissueswith
the non-diagonal terms: each term is small ∼ 3×10−3 as HeretheboxsizeLisinunitsofh−1Gpcandwave-numbers
comparedwiththesignal,buttherearemanyofthem.The are in units of hMpc−1. In spite of the fact that the fit has
fact that individual components are so small is one of the 10 free-parameters, it is still a useful approximation.
(cid:13)c 0000RAS,MNRAS000,000–000
10 Klypin & Prada
Figure 6. Slices through the z = 0 covariance matrix C(k,k(cid:48)) at different values of k(cid:48) for ∼ 4500 realisations of PPM-GLAM
A1.5. Left: Covariance coefficient C(k,k(cid:48))/[C(k,k)C(k(cid:48),k(cid:48))]1/2 (top panel) and covariance matrix normalised to the power spectrum
[C(k,k(cid:48))/P(k)P(k(cid:48))]1/2(bottompanel)fork(cid:48)=0.12,0.2,0.3,0.5hMpc−1.Notethechangeinscaleofthey-axes.Thecovariancematrix
normalised by the signal is very small ∼ 3×10−3 and relatively flat as compared with the large differences seen in the covariance
coefficientsonthetoppanel.Right:Detailedviewsofthecovariancematrix.SolidcurveswiththeerrorbarsarefortheA1.5simulations.
Results for the A0.9 simulations scaled to the 1.5h−1Gpc volumes are shown with dashed curves. Dotted curves show the analytical
approximationgiveninEqs.(11-13).Askincreasesthecovariancematrixfirstdecreases,reachestheminimumatk≈(0.1−0.2)hMpc−1
andthenhasamaximumatk≈(0.5−0.6)hMpc−1.Inaddition,ithasadiponbothsidesofthediagonalcomponents.
The approximation for C(k,k(cid:48)) is so complicated be- The evolution of the covariance matrix with redshift
cause the covariance matrix has a complex dependence on is equally complex as that illustrated in Figure 7 that
k and k(cid:48). For a fixed k(cid:48) the covariance matrix declines showsC(k,k(cid:48))fortheA1.5simulationsatdifferentredshifts.
with increasing k. For example, the covariance matrix de- Curves on the right panels were scaled up with the linear
clines by a factor of ∼ 2.5 from k = 0.03hMpc−1 as com- growthfactorδ(t).Resultsindicatethatonverylongwaves
pared to k = 0.1hMpc−1. It reaches a flat minimum at k<0.07hMpc−1 the covariance matrix grows very fast as
∼
k ≈ (0.1−0.2)hMpc−1 and then increases by factor 1.5 C ∝ δ(t)P(k,t)P(k(cid:48),t) ∝ δ5(t). At the intermediate scales
andreachesamaximumatk≈(0.5−0.6)hMpc−1.Inaddi- 0.1hMpc−1 < k < 0.5hMpc−1 the growth is even faster:
tion, it has a dip on both sides of the diagonal components C ∝δ5.25(t)
(see Figure 6), which is approximated by the last terms in We also conclude that the covariance matrix decreases
eq. (12). with the increasing computational volume. We already saw
Theapproximationgivenbyeqs.(12-13)areonlyforthe itforthediagonalcomponents.Thesameistrueforthenon-
non-diagonalcomponentsofC(k,k(cid:48)).Forthediagonalcom- diagonal terms as illustrated by the left panels in Figure 8.
ponentsthefallowingapproximationgivesa∼2%-accurate The scaling of the whole covariance matrix with volume is
fit to the A0.9 simulations: trivial: C ∝L−3. It is not an approximation or a fit. It is a
(cid:20) (cid:21)
2 scaling law. Right panels in the figure show the results for
C(k,k)=P2(k) α +A2 , (14)
N C rescaled to the (1.5h−1Gpc)3 volume of the A1.5 simula-
h
tions.Thisisanimportantscaling,whichisoftenforgotten.
where parameter A is a slowly increasing function of wave-
Errorsintheestimatesofthepowerspectrumoffluctuations
number:
A= 5.5×10−3 (cid:34)1+(cid:18) k (cid:19)1/2− k(cid:35), k<1hMpc−1. (baenedntthouoslathrgeeeirfroonrseiwnecroesmtoosltoagcickatlopgaertahmeretmerasn)ywsomualdllhsiamve-
L3/2 1.5 4 ulations to mimic large observational volume and forget to
(15) re-scale the covariance matrix.However,whentherescaling
(cid:13)c 0000RAS,MNRAS000,000–000