Second-order Lagrangian perturbation theory initial conditions for resimulations

0.48 MB·English

by  Adrian Jenkins

Mon.Not.R.Astron.Soc.000,000–000 (0000) Printed11January2010 (MNLATEXstylefilev2.2)
Second-order Lagrangian perturbation theory initial
conditions for resimulations.
Adrian Jenkins⋆
Institute forComputational Cosmology, Department of Physics, Universityof Durham, South Road, Durham, DH1 3LE, UK
0 11January2010
1
0
2 ABSTRACT
We describe and test a new method for creating initial conditions for cosmologicalN-
n
body dark matter simulations based on second-orderLagrangianperturbation theory
a
J (2lpt). The method can be applied to multi-mass particle distributions making it
suitable for creating resimulation, or ‘zoom’ initial conditions. By testing against an
1
1 analytic solution we show that the method works well for a spherically symmetric
perturbationwithradialfeaturesrangingovermorethanthreeordersofmagnitudein

] linearscaleandelevenordersofmagnitudeinparticlemass.Weapplythemethodand
O repeat resimulations of the rapid formation of a high mass halo at redshift ∼ 50 and
C the formation of a Milky-Way mass dark matter halo at redshift zero. In both cases
. we find that many properties of the final halo show a much smaller sensitivity to the
h
startredshiftwiththe2lptinitialconditions,thansimulationsstartedfromZel’dovich
p
initial conditions. For spherical overdense regions structure formation is erroneously
-
o delayedinsimulationsstartingfromZel’dovichinitialconditions,andwedemonstrate
r for the case of the formation the high redshift halo that this delay can be accounted
t
s for using the sphericalcollapse model.In addition to being moreaccurate,2lptinitial
a conditions allow simulations to start later, saving computer time.
[

Key words: cosmology:theory – methods: N-body simulations
2
v
8
5
2 1 INTRODUCTION tially with time, it became feasible in the 1990s to sim-
0
ulate the formation of single virialised halos in the CDM
. Computer simulations of cosmological structure formation
0 modelwithenoughparticlestobeabletoprobetheirinter-
havebeencrucialtounderstandingstructureformationpar-
1 nal structure (e.g. Navarro et al. 1996, 1997; Ghigna et al.
ticularly in the non-linear regime. Early simulations e.g.
9 1998; Moore et al. 1999). These more focused simulations
0 Aarseth et al.(1979)usedonlyaroundathousandparticles required new methods for generating the initial conditions.
: tomodelalargeregionoftheUniverse,whilerecentsimula-
v The algorithm of Hoffman & Ribak (1991) for setting up
tionsofstructureformationhavemodelledoverabillionpar-
i constrained Gaussian random fields was one method which
X ticles in just a single virialised object (Springel et al. 2008;
could be applied to setting up initial conditions for a rare
Stadel et al. 2009).SincetheadventoftheCold DarkMat-
r peak where a halo would be expected to form. The essence
a ter (CDM) model (Peebles 1982; Daviset al. 1985), most
ofthistechniqueisthatitallows selection ofaregion based
computational effort has been expended modelling struc-
onthepropertiesofthelineardensityfield.However,theob-
ture formation in CDM universes. Early work in the 1980s
jects we actually observe in the Universe are the end prod-
focused on ‘cosmological’ simulations where a representa-
ucts of non-linear structure formation and it is desirable, if
tiveregionoftheuniverseismodelled,suitableforstudying
wewanttounderstandhowthestructureweseeformed, to
large-scale structure (Davis et al. 1985). The starting point
beabletoselectobjectsonthebasisoftheirfinalproperties.
fortheseCDMsimulations,theinitialconditions,areGaus-
sianrandomfields.Thenumericaltechniquesneededtocre- This requirement led to an alternative method for
ate the initial conditions for cosmological simulations were producing initial conditions based first on selecting ob-
developedinthe1980sandaredescribedinEfstathiou et al. jects from a completed simulation e.g. Katzet al. (1994),
(1985). Navarro& White(1994).Theinitial conditionsforthefirst
As the algorithms for N-body simulations have im- or parent simulation, were created using the methods out-
proved, and the speed of computers increased exponen- linedbyEfstathiou et al.(1985).Thedensityfieldiscreated
out of a superposition of plane waves with random phases.
Having selected an object at the final (or any intermedi-
⋆ [email protected] ate) time from the parent simulation a fresh set of initial
(cid:13)c 0000RAS
2 A. Jenkins
conditions with highernumericalresolution in theregion of In this Section we give a brief summary of Lagrangian
interest, which we will call ‘resimulation’ initial conditions perturbation theory (Subsection 2.1), concentrating on the
(alsocalled‘zoom’initialconditions)canbemade.Particles second-order approximation, and illustrate the difference
of different masses are laid down to approximate a uniform theymakeforthesphericaltop-hatmodelinSubsection2.2.
mass distribution, with the smallest mass particles being The method we describe could be developed further to cre-
concentratedintheregionfromwhichtheobjectforms.The ate third order Lagrangian perturbation initial conditions.
new initial conditions are made by recreating the the same However,Scoccimarro(1998)concludedthatthegainsfrom
plane waves as were present in parent simulation together going from 2nd to 3rd order were relatively modest while
with theaddition of newshorter wavelength power. the complexity of computing the third order terms is con-
An alternative technique for creating resimulation ini- siderable.
tialconditionswasdevisedbyBertschinger(2001)basedon
earlierworkbySalmon(1996)andPen(1997)whereaGaus-
2.1 Summary of second-order Lagrangian
sian random field with a particular power spectrum is cre-
perturbation theory
atedstartingfromawhitenoisefield.Recentlyparallelcode
versions using this method has been developed to generate We give only a bare-bones account of second-order pertur-
verylargeinitialconditions(Prunet et al.2008;Stadel et al. bation theory here. Detailed discussions of Lagrangian per-
2009).Afeaturecommontoboththesetechniques,todate, turbation theory in general can be found in the literature
is that the displacements and velocities are set using the (e.gBuchert1994;Buchert et al.1994;Bouchet et al.1995;
Zel’dovich approximation (Zel’dovich 1970), where the dis- Scoccimarro1998).WefollowthenotationusedinAppendix
placements scale linearly with the linear growth factor. D1of Scoccimarro (1998).
It has long been known that simulations starting from In Lagrangian perturbation theory the Eulerian final
Zel’dovichinitialconditionsexhibittransientsandthatcare comoving positions x are related to the initial positions q
must be taken choosing a sufficiently high start redshift so via a displacement field,Ψ(q):
thatthesetransientscandecaytoanegligibleamplitude(e.g
x=q+Ψ. (1)
Efstathiou et al. 1985). A study of the behaviour of tran-
sientsusingLagrangianperturbationtheorybyScoccimarro It can beshown forsecond-orderLagrangian perturba-
(1998), showed that for initial conditions based on second- tion theory (2lpt) that the displacement field is given, in
order Lagrangian perturbation theory, the transients are terms of two potentials φ(1) and φ(2), by:
both smaller and decay more rapidly than first-order, or
x=q−D ∇ φ(1)+D ∇ φ(2), (2)
Zel’dovich, initial conditions. Scoccimarro (1998) gave a 1 q 2 q
practicalmethodforimplementingsecond-orderLagrangian where D is the linear growth factor, and D the second
1 2
initial conditions. The method has been implemented in order growth factor, which is given by D ≈ −3D2/7. The
2 1
codes for creating cosmological initial conditions by Sirko subscripts q refer to partial derivatives with respect to the
(2005) and Crocce et al. (2006). These codes are suitable Lagrangian coordinates q. Similarly, the particle comoving
for making initial conditions for cosmological simulations, velocities v are given to second order by:
targeted at large-scale structure, but they do not allow the
v=−D f H∇ φ(1)+D f H∇ φ(2), (3)
creationofresimulationinitialconditions,thefocusofmuch 1 1 q 2 2 q
current work on structureformation. where H is the Hubble constant, the fi = dlnDi/dlna
In this paper we describe a new method for creat- (a is the expansion factor). For flat models with a non-
ingsecond-orderLagrangian initialconditionswhichcanbe zero cosmological constant, the following relations apply
used to make resimulation initial conditions. The paper is f1 ≈Ω5/9,f2 ≈2Ω6/11 (Bouchet et al.1995),whereΩisthe
organised as follows: in Section 2 we introduce 2lpt theory matter density. In practice, these approximations for f1,f2
initial conditions and motivate their advantages for study- andD2areextremelyaccuratewhenappliedtomakingmost
ing structure in the non-linear regime by applying them to actualΛCDMinitialconditions,asthestartredshiftishigh
the spherical top-hat model; in Section 3 we describe how enough that Ω is very close tounity.
Zel’dovichresimulation initialconditionsaremade,andthe The potentials φ(1) and φ(2) are obtained by solving a
newmethodforcreating2lptinitialconditions,andtestthe pair of Poisson equations:
method against an analytic solution for a spherically sym- ∇2φ(1)(q)=δ(1)(q), (4)
metric perturbation; in Section 4 we apply themethod and q
analysetheformationofadarkmatterhaloathighandlow where δ1)(q) is the linear overdensity,and
redshift for varying starting redshifts for both Zel’dovich
∇2φ(2)(q)=δ(2)(q). (5)
and 2lpt initial conditions; in Section 5 we summarise and q
discuss the main results; in the Appendix we evaluate two We will refer to the term δ(2)(q), as the ‘second-order
alternateinterpolationschemeswhichhavebeenusedinthe overdensity’. In fact, this quantity is related to the linear
process of making resimulation initial conditions. overdensity field by thefollowing quadratic expression:
δ(2)(q)= φ(1)(q)φ(1)(q)−[φ(1)(q)]2 , (6)
,ii ,jj ,ij
Xi>j(cid:16) (cid:17)
2 THE MOTIVATION FOR USING SECOND
where we use φ ≡∂2φ/∂q ∂q for short.
ORDER LAGRANGIAN PERTURBATION ,ij i j
The most common technique used in making cosmo-
THEORY INITIAL CONDITIONS
logical initial conditions, to date, is the Zel’dovich approx-
. imation (Zel’dovich 1970). Zel’dovich initial conditions can
(cid:13)c 0000RAS,MNRAS000,000–000
2lpt initial conditions 3
equivalently be called first order Lagrangian initial condi- equations of motion of for all the particles forward in time.
tions,andaregivenbyignoringthetermswithφ(2) inequa- For a real simulation the accuracy of the subsequent solu-
tions (2) and (3). A variant of Zel’dovich initial conditions tion would depend not only on the accuracy of the initial
was advocated by Efstathiou et al. (1985). In this case the conditions, but also on numerical aspects such as the time-
velocitiesarenotsetproportionaltothedisplacements,but integrationschemeandtheforceaccuracy.Forourpurposes,
proportionaltothegravitationalforceevaluatedforthedis- wewill assume that thesenumericalsources of errorcan be
placedparticledistribution.Thismethodforsettingtheve- overcome, and will only consider the effects of using either
locities givesweakertransients,asdiscussed in Scoccimarro Zel’dovich or 2lpt initial conditions on the subsequent evo-
(1998),althoughthescalingofthetransientswithredshiftis lution of thetop-hat.
unaffected. In this paper we use the term Zel’dovich initial The Zel’dovich and 2lpt initial conditions needed to
conditionstomeanthevelocitiesareproportionaltothedis- start a cosmological simulation of a top-hat are most easily
placements.Wenote,andwewillexploitthefactlater,that obtainedusingtheratherlesselegantpower-seriesrepresen-
the individual terms on the right-hand side of equation (6) tationof thesphericaltop-hatsolution. Thepower-seriesin
arethederivativesofthefirstorderZel’dovichdisplacements t corresponding toequation (7) is:
which are linear in thepotential φ(1).
R 12πt 23 1 12πt 34 3 12πt 2
r= 0 − − ... ,
4 "(cid:18) tc (cid:19) 20(cid:18) tc (cid:19) 2800(cid:18) tc (cid:19) #
2.2 The Spherical top hat model.
(9)
Thespherical top-hatmodel(Peebles 1980) providesacon- and thecorresponding velocity, r˙ is:
venient way to illustrate the relative merits of Zel’dovich 2πR 12πt −31 1 12πt 13 9 12πt
and 2lpt initial conditions. Apart from the significant ad- r˙ = 0 − − ... .
vantage of mathematical convenience, there are other rea- tc "(cid:18) tc (cid:19) 10(cid:18) tc (cid:19) 2800(cid:18) tc (cid:19) #
sons for considering this special case: (10)
The leading order term in equation (9) represents the
(i) The first collapsed structures in a cosmological sim- expansionofaperturbationwithmeancosmicdensitywhich
ulation will form from the regions with the highest linear grows in time in the same way as the scale factor for an
overdensity. In the limit of high peak height, the regions
Einstein-de Sitter universe. The second to leading term
around the maxima in a smoothed linear overdensity field gives the linear growth of an overdense spherical grow-
becomespherical (Bardeen et al.1986).Aswewill see, it is ing mode perturbation and it is easy to verify from equa-
the first structures which are most sensitive to whether the tion (9) that the linear overdensity at t = t , when the
c
initial conditions are Zel’dovich or 2lpt. non-linear solution collapses to a point, is given by δ =
c
(ii) It is easy to show that the second-order overdensity, 3(12π)2/3/20) ≃ 1.686, familiar from Press-Schechter the-
δ(2), is related to the first order overdensity, δ(1), through ory (Press & Schechter1974).
the inequality: δ(2) ≤ 13(δ(1))2, with equality occurring WeobtaintheZel’dovichinitialconditionforthespher-
for isotropic compression or expansion. For a fixed linear ical top-hat by taking just the first two leading terms in r
overdensity, the second order overdensity is maximised for and r˙ in equations (9) and (10). The equivalent 2lpt initial
isotropic compression, so the spherical top-hat is the limit- conditionscorrespondtotakingthefirstthreetermsinstead.
ing case where the affects of second-order initial conditions We can follow the top-hat solution starting from these ini-
are maximised. tialconditions,upuntilthepointwhenthespherecollapses
to zero radius, by integrating the equation of motion equa-
For convenience we will consider, for now, the top-hat
tion (8) forward in time from time, t, using the values of r
in an Einstein-de Sitter universe. In this case the spherical
and r˙ at time t as theboundary conditions.
top-hatsolutionismostelegantlyexpressedinaparametric
In Figure 1 the black solid line shows the parametric
form as follows. For an overdense spherical top-hat which
top-hat solution given by equation (7). The red solid curve
beginsexpandingatt=0andreachesamaximumphysical
shows a spherical top-hat model integrated forward from
radius of expansion, R , and collapses back to zero radius
0
the arbitrarily chosen time, t = t /27, starting from the
at time, t , theparametric solution is: c
c
Zel’dovich initial condition (marked with a red dot). The
R t
r= 0 (1−cosη) ; t= c (η−sinη), (7) cyan dashed curve shows the corresponding result starting
2 2π with the 2lpt initial condition at t = t /27. The solution
c
where0≤η≤2π.Theseequationscanbederivedbysolving using the 2lpt initial condition is very much closer to the
theequation of motion for thetop-hat: analyticsolutionthanthecurvestartedfromtheZel’dovich
initial condition.
π2R3 1
r¨=− 0 , (8) Thecollapsetimeforthetop-hatforboththeZel’dovich
2t2 r2
c and2lptinitialconditionsoccurslaterthantheanalyticso-
and applying theappropriate boundary conditions. lution. This delay, or timing error, ∆t, can be quantified as
If we were to solve the spherical top-hat problem in- follows. Multiplyingequation (8) byr˙ andintegrating gives
steadbythetechniquesusedincosmological numericalsim- an integral of the motion: E = r˙2/2−π2R3/2t2r. For the
0 c
ulations, we would begin by creating initial conditions con- analytic solution this can easily be evaluated at the radius
sisting ofaset ofparticles with givenmasses, positions and of maximum expansion which gives E = −π2R2/2t2. The
0 c
velocities,atatime,t,soonafterthestartoftheexpansion. corresponding values of E for the Zel’dovich and 2lpt solu-
The solution would be obtained by evaluating the gravita- tions can be obtained by substituting the first two or three
tionalinteractionsbetweentheparticlesandintegratingthe termsrespectivelyfrom equation(9)andequation(10)into
(cid:13)c 0000RAS,MNRAS000,000–000
4 A. Jenkins
Similarly,the2lptinitialconditionfractionaltimingerroris
given by:
∆t 23δ2 t 4/3
= c s Ω−0.23. (12)
t 210 t c
c (cid:18) c(cid:19)
Alternatively, these results can be expressed in terms
of the start redshift of a numerical simulation, z . Relative
s
to an analytic spherical top-hat with collapse occurring at
redshift, z , where z ≪ z , the fractional timing error is
c c s
given for Zel’dovich initial conditions by:
∆t 0.51 1+z
≈ c , (13)
t Ω0.36 1+z
c c (cid:18) s(cid:19)
and for 2lpt initial conditions it is:
∆t 0.31 1+z 2
≈ c , (14)
t Ω0.49 1+z
c c (cid:18) s(cid:19)
where we have made use another approximate relation:
t2/3(1+z) ∝ Ω0.13 relating redshift, z, to the age of the
Figure 1. Evolution of a top-hat spherical perturbation. The universe, t and the value of Ω all at time, t, which is accu-
blackcurveshowstheradiusofasphericaltop-hatperturbationas rate to better than one percent in the range 0.25 < Ω ≤ 1
afunctionoftimeforanEinstein-deSitteruniverse.Theanalytic for a flat universewith a cosmological constant.
solutionstartsexpanding at t=0andcollapsesatatimet=tc Ignoring the weak Ω dependence, the timing error for
and has a maximum physical radius of expansion of R0. The the2lptinitialconditionshowsaquadraticscalingwiththe
red and dashed cyan curves show the spherical top-hat solution ratiooftheexpansionfactoratthestarttothatatthetime
obtained by integrating the equation of motion for the top-hat
of collapse, while for Zel’dovich initial conditions, the scal-
starting from t = tc/27 and respectively initialising the radius ing is just linear. These results for spherical collapse are
and its velocity with the Zel’dovich and 2lpt initial condition.
consistent with the more general behaviour of the decay
The result for the 2lpt initial condition is more accurate than
of transients discussed in Scoccimarro (1998). Theabsolute
for the Zel’dovich initial condition where the time of collapse is
significantlydelayed. value of the timing error, ∆t, for a fixed start redshift, ac-
tuallydecreaseswithincreasingcollapsetimefor2lptinitial
conditions. For Zel’dovich initial conditions, the timing er-
rorgrows slowly with increasing collapse time.However,by
theexpressionforE.Theexpansiontime,t ,ofthetop-hat startingatsufficientlyhighredshift,thetimingerrorforthe
c
over the complete cycle, by analogy with Keplerian orbits, Zel’dovichinitialconditionscanbemadesmall.Itiscommon
isrelatedtoE by:dlnt /dln|E|=−2 forboundmotionin practice to start simulations from Zel’dovich initial condi-
c 3
anEinstein-deSitteruniverse.Wecangeneralise thisresult tionsatahighredshiftforthisreasonsothatthetransients
to flat universes with a cosmological constant by numeri- havesufficient time to decay toa negligible amplitude.
cally integrating equation (8) with an added term for the Equating the fractional timing errors for equations
cosmological constant.Inthiscasewefind,byfittingtothe equation (13) and equation (14) gives an objective way of
results of numerical integration, an equivalent approximate defininganequivalentstartingredshiftforZel’dovichor2lpt
relation: dlnt /dln|E| = −2Ω0.23, valid to an fractional initialconditionsforagivenchoiceofcollapseredshift.Tobe
c 3 c
accuracy of less than 1.2 percent for 0.25 < Ω ≤ 1, where conservative,itwouldbenaturaltochoosethecollapsered-
c
Ω is the matter density in units of the critical density at shifttobethefirstoutputredshiftofthesimulationwhichis
c
the time of collapse, provided t ≪ t so the cosmological ofscientificinterest.Togiveaconcreteexample,ignoringthe
s c
constant is negligible at t . Ωdependence,ifthe2lptinitialconditionsarecreatedafac-
s
Knowing dlnt /dln|E| and assuming again the condi- torof10inexpansionbeforethefirstoutput,theequivalent
c
tionthatthestartingtimeoftheinitialcondition,t ,obeys, Zel’dovich initial conditions would need to start at a factor
s
t ≪t , allows the fractional timing error, ∆t/t , to be ob- ofabout160inexpansionbeforethefirstoutput.Using2lpt
s c c
tained, to leading order in t /t , for both types of initial initialconditionsinsteadofZel’dovichinitialconditionscan
s c
condition,from theratioofE fortheanalyticsolution with therefore lead to a significant reduction in computer time
respectivelywiththevalueofE forZel’dovichor2lptinitial because the simulations can be started later. A later start
conditions. hasfurtheradvantagesasforceerrors,errorsduetothedis-
Applying these results we find after some algebra that creteness in the mass distribution and time integration er-
the timing error for Zel’dovich initial conditions to lowest rorsallaccumulateastheratioofthefinalexpansionfactor
order is: totheinitialexpansionfactorincreases.Allofthesereasons
make 2lpt initial conditions attractive for the simulator. In
∆t 3δ t 2/3 the next section we describe the new method to make 2lpt
= c s Ω−0.23. (11)
t 10 t c initial conditions.
c (cid:18) c(cid:19)
(cid:13)c 0000RAS,MNRAS000,000–000
2lpt initial conditions 5
3 MAKING SECOND ORDER LAGRANGIAN In this paper we will not be concerned with steps (i)
RESIMULATION INITIAL CONDITIONS and(ii),whicharecommon toresimulations ingeneral,but
willconcentrateonadaptingsteps(iii)and(iv)toallowthe
In this Section we explain what resimulation initial condi-
creation of 2lpt initial conditions.
tions are and describe how they can be made. We will de-
Thedisplacement andvelocity fieldsneededtoperturb
scribethemethodcurrentlyinusetomakeinitialconditions
theparticle load can in practice be efficiently created using
for the Virgo Consortium for Supercomputer Simulations1.
Fourier methods. Adopting periodic boundary conditions
The Virgo resimulation code to date generates Zel’dovich
when setting up the perturbations forces the Fourier rep-
initial conditions. InSubsection3.1 we describetheprocess
resentation of the density, displacement and velocity fields
for making these Zel’dovich resimulation initial conditions.
tobecomposedofdiscretemodesinFourierspace.Theam-
In Subsection 3.2 we describe the new method to generate
plitudes and phases of these modes are set as described in
2lptinitialconditions.InSubsection3.3wedemonstrateby
Efstathiou et al. (1985) to create a Gaussian random field
comparingwithananalyticsolution thatthemethodworks
withtheappropriatepowerspectruminFourierspace.Start-
in practice.
ing from theFourier representation of the density field it is
straightforward to use discrete Fourier transforms to gen-
3.1 Setting up Zel’dovich resimulation initial erate each of the components of the displacement field in
conditions realspace.ThediscreteFouriertransformgivesthedisplace-
mentson a cubiclattice.
ThecurrentmethodusedbytheVirgocodeisbasedonthe Inpractice,becausecomputermemoryislimited,more
ideasfirstdescribedinNavarro & White(1994).Theactual than one Fourier grid is required to be able to complete
implementationhaschangedconsiderablyovertimewithin- step (iv).These grids havedifferent physicalsizes, with the
creasingly up-to-date accounts being given in Power et al. largest grid covering the entire periodic cubic volume, as is
(2003), Navarro et al. (2004) and Springel et al. (2008). the case when making cosmological initial conditions. The
Theneedforresimulation initial conditionsstemsfrom additional grid(s) are given a smaller physical size, which
thedesiretorepeat aparticularsimulation ofstructurefor- means even with limited computer memory, that they can
mationwith(usually)highermassresolutioninoneormore generateadisplacementfieldwithpowerathigherwavenum-
regions than was present in the original, or parent, simu- bersthanthelargest grid,butonlyoverarestricted partof
lation. The parent simulation may have been run starting the simulation volume. By nesting sufficient Fourier grids
from cosmological initial conditions or it may be a resim- about the high resolution region it is possible to generate
ulation itself. The following list gives the steps needed to power all they way down to the particle Nyquist frequency
create a set of resimulation initial conditions from a parent inthehighresolutionregion.Adescriptionofthecreationof
simulation. theinitialconditionsfortheAquariushalos,whichrequired
two Fourier grids, is given in Springel et al. (2008).
(i) Identify theregion(s) of interest in theparent simula-
For Zel’dovich initial conditions, thedisplacement field
tion where higher mass resolution is required.
atanylocationisformedbycoaddingthedisplacementcon-
(ii) Create a force-free multi-mass particle distribution
tributions from the one or more of these overlapping cubic
which adequately samples the high resolution region(s) of
Fouriergrids, andthecorresponding velocity field issimply
interest identified in (i) with low mass particles. This par-
proportionaltothedisplacementfield.ForeachFouriergrid
ticle distribution needs to be a good approximation to an
only a subset of the possible modes are assigned non-zero
isotropic and homogeneous mass distribution. The initial
values.Ink-spacethese‘populated’modesaresplitbetween
densityfluctuationsarecreated byperturbingthepositions
thedifferentgridsinanestedfashionaroundtheorigin.Each
and assigning velocities to each of theparticles. The region
gridpopulatesadisjointregionofk-spaceconcentriconthe
exterior to the high resolution region can often be sampled
origin,andtakentogetherthemodesfromallthegridsfilla
more crudely than theparent simulation, using larger mass
spherical region in k-space around the origin. The smallest
particles, as its only purpose is to ensure that the gravita-
physical grid, which typically fits around just the high res-
tional tidal field in the high resolution region is accurately
olution region,isusedtorepresent thehighest wavenumber
reproduced.Typicallythemassesintheexteriorregionscale
modes in k-space. The highest wavenumber possible is lim-
roughlyasthecubeofthedistancefromthehighresolution
ited by the particle Nyquist frequency which is determined
region.Wewillrefertothisspecialparticledistributionfrom
bythe interparticle separation.
now on as the‘particle load’.
A displacement field is generated on each Fourier grid.
(iii) Recreate the displacement and velocity fields for the
Thecompletedisplacementfieldiscalculatedbyaddingthe
parent simulation and apply the displacements and assign
thevelocities to theparticle load created in (ii). thecontributionsfromthevariousgridstogether.Todothis
it isnecessary todefineaset of common positions at which
(iv) Add extra small-scale fluctuations to the region of
toco-addthedisplacements.Themostnaturalcommonpo-
interestonly,toaccountforthefact thatthemassdistribu-
sitions to use are those of the particle load itself, as the
tionintheresimulationisnowmorefinelysampledthanthe
displacements need to be known at these locations to gen-
parentsimulationandthereforehasahigherspatialNyquist
eratetheZel’dovichinitialconditions.Thedisplacementsat
frequency.ForZel’dovich initial conditions thishigh spatial
theparticlepositionsarecomputedbyinterpolatingtheval-
frequencypowercanbeaddedlinearly tothelowfrequency
uesfromthenearestFouriergridpoints.Adiscussionofthe
power copied from theparent simulation.
interpolation methods currently in use and their associated
errorsisgivenforcompletenessintheAppendix.Theinter-
1 http://www.virgo.dur.ac.uk/ polationerrorsonlybecomesignificantforfluctuationswith
(cid:13)c 0000RAS,MNRAS000,000–000
6 A. Jenkins
a wavelength comparable to the Nyquist wavelength of the sitions. Computationally this entails storing six additional
Fourier grid. quantitiesper particle in computer memory.Weuse theQI
For the high resolution region a second scale, the par- interpolation scheme, described in the Appendix, to inter-
ticle Nyquist wavelength also becomes important. Even in polate the derivatives to the positions given by theparticle
the absence of interpolation errors, the discreteness of the load, as this scheme helps to minimise theinterpolation er-
mass distribution affects the growth of linear fluctuations rors,givenlimitedcomputermemory.Wehaveimplemented
with wavelengths close to the particle Nyquist frequency the algorithm to compute δ(2) for each particle in the par-
(Joyce et al. 2005). We will not be concerned about either ticle load in a Fortran/C code called IC 2LPT GEN. This
of these sources of error in themain part of the paper. codealso computestheZel’dovichdisplacementsandveloc-
ities for each particle at thesame time.
Thefinalstagetomakethe2lptresimulationinitialcon-
3.2 Generating 2lpt resimulation initial conditions
ditions is to solve the Poisson equation (5). At thispoint it
In Subsection 3.1, we outlined a method for making isnolongertrivialtosolvetheequationbyFouriermethods
Zel’dovich initial conditions for resimulations. In this Sub- because thesecond-order overdensity field is known only at
section we will take this method for granted, and concen- the positions of the particle load, which do not in general
trateonhowtoaugmenttheZel’dovichinitialconditionsto form a regular grid. However it is relatively simple to solve
produce 2lpt initial conditions. This is achieved, as can be the Poisson equation using an N-body code Poisson solver.
seen from equations (2) and (3), by computing the com- It is natural, given that the purpose of creating the initial
ponentsofthefield∇ φ(2).Thesecond-ordercorrectionsto conditions in the first place is to integrate them forward in
q
the Zel’dovich initial conditions for both the positions and time with an N-body code, to use the same N-body code
thevelocities are proportional to this quantity. tosolveboth equation(5) andtointegrate theequationsof
Thesecond-orderpotential,∇ φ(2),isobtainedbysolv- motion.
q
ing the Poisson equation (5), where the source term is the As mentioned in Subsection 3.1, in point (ii), the La-
second-order density field, δ(2), given by equation (6). The grangian positions of the particle making up the particle
right-handsideof thisequation isa non-linearcombination load are special. The particles are distributed so that the
ofsixdistinctterms,eachofwhicharesecondderivativesof massdensityfieldisanexcellentapproximationtoauniform
thefirstorderpotential,φ(1).Thissamefirstorderpotential mass density,on scales larger than theinterparticle separa-
isusedtogeneratetheZel’dovichdisplacementsandveloci- tion.Furthermorethegravitationalforceonanyparticledue
ties. For this reason, as pointed out by Scoccimarro (1998), totheotherscancelsouttoaverygood approximation.We
the calculation of the second-order overdensity field, fits in can exploit the properties of the particle load in a simple
verywellwiththegenerationofZel’dovichresimulation ini- waytoreconstructthesecond-orderdensityfield.Todothis
tialconditions,fromacomputationalpointofview.Thesix wecreatea‘2lpt particle load’ bytakingthepositions from
second derivatives can be calculated using Fourier methods the particle load, but replacing the particle masses with a
in the same way as is used to compute the first derivatives ‘2lpt’ mass, given by:
of φ(1) needed for theZel’dovich approximation. m =m(1+λδ(2)), (15)
2LPT
For initial conditions that can be made using a single
Fourier grid, the six derivatives can be combined to give whereλisaconstantnumericalfactor,whosevaluefornow
the value of δ(2) at each Fourier grid point. The Poisson we will taketo be unity.
equation (5), can then be solved immediately using stan- With minimal modification, an N-body code can then
dard Fourier methods (Hockney & Eastwood 1988) to give be used to compute the gravitational force due to the 2lpt
the components of ∇ φ(2) at the grid points. Choosing a particle load acting on itself. By construction this amounts
q
suitable cubic grid arrangement for the particle load en- to solving equation (5), with the resulting gravitational
suresthattheLagrangian positionsoftheparticlescoincide acceleration giving, ∇ φ(2), at the positions given by the
q
exactly with Fourier grid points. This means the particle particle load. Provided the N-body code also reads in the
2lpt displacements and velocities can be evaluated without Zel’dovichvelocitiesfortheparticleload,the2lptinitialcon-
theneed for any spatial interpolation. This is theapproach dition positions andvelocities can beevaluated using equa-
adopted by the 2lpt codes made available Sirko (2005) and tions (2) and (3) respectively. Discarding the 2lpt masses,
Crocce et al. (2006). andreplacingthemwith thetrueparticle massescompletes
The approach just described does not work when more the creation of the 2lpt resimulation initial conditions. The
than one Fourier grid is used to generate the Zel’dovich N-bodycodecanthenbeusedinthenormalwaytointegrate
initial conditions, because δ(2) is a non-linear function of theequations of motion for all of the particles.
thefirst-orderpotential,φ(1),unliketheZel’dovichdisplace- The amplitude of the second-order overdensity field,
mentsandvelocitieswhicharelinearfunctions.Thesolution and similarly ∇ φ(2) depends on the redshift of the initial
q
to this problem is simply to compute all six of the deriva- conditions.Byintroducingthescalingfactorλintothe‘2lpt’
tives,whicharelinearfunctionsofφ(1),foreachFouriergrid, masses, and rescaling the resulting ∇ φ(2) by 1/λ to com-
q
andco-addthem.Thenoncethesummationiscompletefor pensate,itispossibletomakethecalculation performedby
all six derivatives for a given particle, compute δ(2) using theN-bodycodeindependentoftherequestedredshiftofthe
equation (6). initialconditions.Itisdesirabletoimplementthisscalingif
Aswith theZel’dovich displacements,thevaluesof the one wants to calculate ∇ φ(2) accurately for high redshift
q
six derivatives need to be interpolated to a common set of initial conditions, as otherwise the inevitable force errors
positions before theycan be co-added. It provesconvenient from the N-body force calculation will begin to dominate
again to use the particle load to define these common po- over the small contribution from the second-order overden-
(cid:13)c 0000RAS,MNRAS000,000–000
2lpt initial conditions 7
sityfield.Theprecisevalueofλisunimportant,providedit malisation and thesame cosmological parameters) thenthe
islargeenough.Asuitablechoiceisgivenbysettingλtobe outputof thetwo codes agrees extremely well.
as large as possible while avoiding any of the ‘2lpt’ masses Testing the 2lpt resimulation initial conditions gener-
becoming non-positive. ated usingthemethod described in theprevioussection for
We have implemented the final step to make the 2lpt ageneral case is problematic asit is non-trivialtocompute
initialconditionsinpracticewiththeP-Gadget3code,writ- thesecond ordertermsbyalternative means.However,it is
ten by Volker Springel. P-Gadget3 is an N-body/SPH code relativelysimpletotestthecodeforasphericallysymmetric
written in C. This code is similar to the public Gadget-2 perturbation for which one can calculate the properties of
code (Springel 2005), and has been used, for example, for the initial conditions analytically. The IC 2LPT GEN code
the Aquarius project (Springel et al. 2008). Operationally needs only minimal alteration to be able to create a spher-
theFortran/Ccode,IC 2LPT GEN,oncompletionoutputs ical perturbation. The only change needed is to replace the
a special format initial condition file(s) for the P-Gadget3 partofthecodewhichgeneratestheamplitudesandrandom
code to read in. The file is special in that the particle po- phasesforeachFouriermode,asrequiredforgeneratingran-
sitions are the Lagrangian positions, i.e. the particle load dom Gaussian fields, with a module that instead computes
positions,thevelocitiesaretheZel’dovichvelocities,andan theappropriateamplitudeandphaseforaspecifiedspherical
additional block of data is provided which gives every par- perturbationwithacentreofsymmetryatagivenlocation.
ticle a second mass - the ‘2lpt-mass’ described above. This In other respects the code is unaltered, which means the
required a modification to P-Gadget3 so that it first can sphericaltest remains apowerfultest of theoverallmethod
recognise the special file format required by the 2lpt initial and theparticular numerical implementation.
conditions through a new flag stored the file header. Hav- The test initial condition most closely resembles what
ing read in the initial conditions, P-Gadget3 first does a would be expected around a high overdensity peak and is
force calculation using the 2lpt masses for each particle as similar to the initial conditions in Section 4 for the highest
described above, and once the 2lpt initial conditions have resolution initial conditions used in the papers Gao et al.
been created, integrates theequationsof motion forward in (2005)andReed et al.(2005).Thesetwopaperswerebased
time as usual. on a simulation of a relatively massive and very rare high
The linear density field created by the resimulation redshift halo, a suitable location for finding one of the first
methodoutlinedaboveisbandwidthlimited(seePress et al. metal free stars.
(1992) for a discussion). The second-order overdensity field For our test we add a spherically symmetric pertur-
is similarly bandwidth limited but because equation (6) is bation, to a homogeneous and isotropic Einstein de-Sitter
quadratic,thesecond-orderoverdensityhastwiceaslargea universe. We adopt a simple analytic form for which the
bandwidth in each physical dimension. If the second-order second-orderoverdensityandthe2lptdisplacementsandve-
overdensity field is known at the Lagrangian particle posi- locitiescanbecomputedanalytically.Theperturbationwas
tions then it is only fully sampled if the linear density is formed by summing n component terms to give the linear
limited to a wavelength which is twice that of the particle overdensity field as follows:
Nyquist wavelength (this approach is implemented in the
Crocce et al. (2006) code). It is common practice in mak- n
ing initial conditions to use waves for the linear density δ(1)(r)= δi(1)(3−r2/σi2)exp(−r2/2σi2), (16)
fielddown totheNyquistwavelength itself. Inthiscase the i=1
X
second-order overdensity field at the positions of the parti-
cles undersamples the second-order density field. However, where r is the radius, δ(1) is the amplitude of the linear
i
asmentionedearlier,thegrowthrateoffluctuationsofvery overdensity,andσ thelinearsize, fortheithcomponentof
i
shortwavelengthsisalsocompromisedbythediscretenessof theperturbation.Theprefactortotheexponentialischosen
the mass distribution in any case whether using Zel’dovich so that: ∞r2δ(1)dr=0, meaning that the net mass of the
0
or 2lpt initial conditions so we will ignore this issue, but it perturbation is zero.
R
remainsthecasethatitisalwaysadvisabletodoresolution The choice of a radial Gaussian cut-off ensures that
tests to establish the limitations of thesimulations. the spherical perturbation is well localised which is a nec-
essaryrequirement:theuseofperiodicboundaryconditions
inmakingtheresimulation initialconditionsviolatesspher-
ical symmetry,and meansthechoice of thescale lengthsσ
3.3 Testing the 2lpt resimulation initial i
has to be restricted or the resulting perturbation will not
conditions
beclosetobeingsphericallysymmetric.Wehavechosenthe
We can test limited aspects of the IC 2LPT GEN code scale lengths, σi, so that each component of the perturba-
againstthecodemadepublicbyCrocce et al.(2006)forthe tionisessentiallygeneratedbyacorrespondingFouriergrid,
casesofbothZel’dovichand2lpt cosmological initial condi- withthevalueofσi beingnogreaterthanabout10percent
tions, where the second-order overdensity field can becom- of the side length of the grid and each grid is centred on
puted using a single Fourier grid. The two codes were writ- the centre of spherical symmetry. These requirements en-
tenindependentlyalthoughtheydosharesomesubroutines sure that each component of the perturbation is close to,
fromPress et al.(1992).WefindthatwhentheCrocce et al. although not exactly, spherical.
(2006) code is modified so as to set up precisely the same Allthatisneededtobuildtheperturbationwiththeini-
density field (i.e. exactly the same modes, with the same tialconditions codeis theequivalentFourier representation
phases and amplitudes for each mode, and the same shape of the real space perturbation equation (16). The Fourier
linearinputpowerspectrum,thesamepowerspectrumnor- representation,δˆ(k)foraperturbationwithoverdensityδ(x)
(cid:13)c 0000RAS,MNRAS000,000–000
8 A. Jenkins
is given by:
δˆ(k)= δ(x)exp(ik·x)d3x, (17)
Z
where k is the wavevector. Equation (16) has the following
Fourierrepresentationasafunctionofthemagnitudeofthe
wavevector k:
n
δˆ(k)=(2π)3/2 δ(1)k2σ5exp(−k2σ2/2). (18)
i i i
i=1
X
The solution of equation (4) is:
n
φ(1)(r)= δ(1)σ2exp(−r2/2σ2). (19)
i i i
i=1
X
Forageneralsphericalperturbation,thesourcetermforthe
Poissonequation(5)simplifiesinsphericalpolarcoordinates
to:
1 d dφ(1) 2
∇2φ(2) = r (20)
r2dr (cid:18) dr (cid:19) !
which can be integrated to give the second order displace-
ment:
Figure2. Comparisonbetweenquantitiesgeneratedbytheres-
dφ(2) 1 dφ(1) 2 imulation codes for a spherical perturbation and the exact ana-
= , (21) lyticalresult.ThetoppanelshowstheZel’dovichvelocities(v(1))
dr r dr
(cid:18) (cid:19) areaccuratelyreproducedbytheIC 2LPT GENcode.Themid-
where the constant of integration has been set to zero to dle panel shows the second-order overdensity (δ(2)), output by
avoid a source term at the coordinate origin. the IC 2LPT GEN code. This quantity is needed as an input to
For the spherical perturbation given by equation (16), theP-Gadget3codetogeneratethesecondorderLagrangiandis-
thepotential is given by: placements and velocities. The bottom panel shows the second
orderLagrangianvelocitiescomputedbytheP-Gadget3codeus-
n n
1 ingthevaluesofthesecond-orderoverdensitycomputedforeach
φ(2)(r)= δ(1)δ(1)σ2 exp(−r2/σ2), (22)
2 i j ij ij individual particle. The arrows markthe sidelengths of the five
Xi=1Xj=1 periodicFouriergridsthatareusedtogeneratetheperturbation
where 2/σ2 =1/σ2+1/σ2. intheresimulationcode.Thedottedverticallinesmarktheshort
ij i j wavelengthcut-offforeachoftheFouriergrids.
The second-order overdensity,given by equation (6) is:
n n 2r2
δ(2)(r)= δ(1)δ(1) 3− exp(−r2/σ2). (23)
i j 3σ2 ij is replaced by the condition kx2+ky2+kz2 < Kmax to en-
Xi=1Xj=1 (cid:18) ij(cid:19) sure a spherical cut-off in k-space to minimise small-scale
p
For the test, a set of initial conditions was made sum- anisotropy.
ming five terms in Eq. (16). The simulation volume was Following stage (ii) of the process described in Subsec-
chosen to be periodic within a cubic volume of side length tion3.1,aspecialparticleloadwascreatedforthetestcon-
100 units and the centre of symmetry was placed in the sisting of concentric cubic shells of particles centred on the
middle of the simulation volume. Five Fourier grids were centre of the volume. In detail, the cubic volume of side
used also centred on the centre of the volume. The param- length 100 was divided into n3 cubic cells. A particle was
eters chosen are given in Table 1. The Table gives the val- placed at thecentreof just those cells on thesurface of the
ues of δ(1) and σ for the five components making up the cube(atotalofn3−(n−2)3 cells).Themassoftheparticle
i i
spherical perturbation. For the purposes of the test the ac- was given bythecell volumetimes thecritical density.The
tual value of the overdensity is unimportant. We arbitrar- remaining(n−2)3 cellsdefineasmaller nestedcubeofside
ily chose the central overdensity to be unity. This is an or- length100×(n−2)/n.Thissmallercubicvolumewasagain
der of magnitude higher than what would typically occur dividedinton3 cells and theprocess of particlecreation re-
in a realistic simulation. The quantityL gives the phys- peated.Theprocesswasrepeated360times,takingn=51.
grid
ical dimension of each Fourier grid. The first grid, which Onthefinaliteration alln3 cells werepopulated,filling the
covers the entire cubic periodic domain, has a dimension hole in the middle. Provided the spherical distribution is
of 100 units on a side. The parameter N is the size centredonthecentreofthecubethisarrangement ofparti-
grid
of the Fourier transform used in the computations. The clesallowsasphericalperturbationcoveringawiderangeof
parameters K and K define which Fourier modes scalestowellsampledeverywhere.Theparticleloadhascu-
min max
are used to make the spherical perturbation. All modes bicratherthansphericalsymmetry,sothisarrangementwill
are set to zero except for the modes in each grid with generate some asphericity, but for sufficiently large n, this
wavenumber components, 2π/L (k ,k ,k ), which obey becomessmall. Similarresultswereobtainedwithasmaller
grid x y z
K ≤ max(|k |,|k |,|k |) < K . These limits apply to value,n=23.
min x y z max
grids1,2,3,4andgrid5,exceptforthehigh-klimit,which Figure 2 shows the results of the test. The top panel
(cid:13)c 0000RAS,MNRAS000,000–000
2lpt initial conditions 9
i δi(1) σi Lgrid Ngrid Kmin Kmax
1 0.1 15 100 512 1 80
2 0.15 2 15 512 12 80
3 0.2 0.2 1.5 512 8 80
4 0.25 0.024 0.18 512 10 80
5 0.3 0.003 0.0225 512 10 80
Table 1. Parameters for a spherical perturbation used to test
the method for generating 2lpt initial conditions. The index i
andquantitiesδi(1) andσi characterisetheperturbationusedfor
the test as given by equation (16). The remaining columns give
parameters of the Fourier grids used to realise the perturbation
andareexplainedinthemaintext.
shows the Zel’dovich velocities (v(1)) as a function of ra-
dius, scaled by the Hubble velocity relative to the centre.
The black line shows the analytic solution, while the red
dots show binned values measured from the output of the
IC 2LPT GEN code. The vertical arrows mark the edge
length of the five cubic Fourier grids, while the five ver-
tical dotted lines show the high spatial frequency cut-off Figure 3. The physical radius of a sphere containing a fixed
amountofmassinasimulationoftheformationofamassivehalo
(given by K ) of the imposed perturbations for the cor-
max
at high redshift. The sphere is centred on a density maximum
responding cubic grid. The first panel establishes that the
located by applying the shrinking sphere algorithm detailed in
IC 2LPT GEN code is able to generate the desired linear
the text. The four curves show simulations starting from either
perturbation.
Zel’dovich or 2lpt initial conditions, at a start redshiftof 399 or
ThesecondpanelofFigure2showsinblackthesecond- 599.Theagreementbetweenthetwo2lptinitialconditionsisvery
orderoverdensityoftheanalyticsolution,whilethereddots good. As expected from the spherical collapse model, discussed
showtheresultsfromtheoutputoftheIC 2LPT GEN.This inSubsection2.2,thecollapseoftheZel’dovichinitialconditions
shows that the second-order overdensity is successfully re- isdelayedwiththelargestdelayoccurringforzs=399.
produced by thecode.
The thirdpanelof Figure 2shows thesecond-orderve-
locity(v(2)).Theblackcurveistheanalyticsolutionandthe 4.1 Structure formation at high redshift
red dots show the binned second-order contribution gener-
The initial conditions for the spherically symmetric test in
atedbytheP-Gadget3code,usingthevaluesofthesecond-
Subsection3.3,wasintendedtoapproximatetheinitialcon-
orderoverdensitycomputedforeachparticlewhicharegen-
ditions in a series of resimulations presented in Gao et al.
eratedbytheIC 2LPT GENcode.Themethodworkswell,
(2005) to investigate the formation of the first structure in
the fractional error in the value of v(2) is below 1 percent
the Universe. In that paper, a cluster was selected from a
everywhere except close to the scale of the first grid where
simulation of a cosmological volume. A resimulation of this
spherical symmetry is necessarily broken and theperturba-
cluster was followed to z = 5, and the largest halo in the
tion has an extremely small amplitudeanyway.
high resolution region containing the proto-cluster was it-
Inconclusionthemethodisabletoworkwellforaper-
self resimulated and followed to z =12. The largest halo in
turbationthatspansaverylargerangeofspatialscalesand
the high resolution region at this redshift was then resimu-
particle masses. The 5th component has a length-scale of
lated and followed to z = 29. This was repeated one more
just 1/5000th ofthe1st component andtheparticlemasses
time ending with a final resimuation (called ‘R5’) yielding
associated with representing these components vary by the
cubeofthisnumber-morethanelevenordersofmagnitude. a halo of mass M200 = 1.2×105h−1M⊙ (that is the mass
within a spherical region with a mean density of 200 times
Most practical resimulations are less demanding, although
the critical density) at redshift 49, simulated with particles
theresimulationofGao et al.(2005),whichwereproducein
thenext Section, is similar. with a mass of just 0.545h−1M⊙. The original motivation
forstudyingsuchahaloisthathalos likethisoneareapo-
tentialsitefortheformation ofthefirstgenerationofmetal
free stars (Reed et al. 2005).
We have recreated the R5 initial conditions with
both Zel’dovich and 2lpt resimulation initial conditions. In
4 REAL APPLICATIONS
Gao et al. (2005), the starting redshift was 599, and we
In this Section we remake some resimulation initial condi- have repeated a simulation starting from this redshift and
tions with Zel’dovich and 2lpt initial conditions and run also z = 399, for both Zel’dovich and 2lpt initial condi-
s
themwiththeP-Gadget3code.Weinvestigatehowtheend tions. We ran each simulation with P-Gadget3 to z = 431
3
point of the simulation depends on the start redshift and and produced outputs at expansion factors of 0.01, 0.015,
thetypeof initial condition. 0.017, 0.019, 0.021, 0.0225. These initial conditions remain
(cid:13)c 0000RAS,MNRAS000,000–000
10 A. Jenkins
Figure 3 shows these five radii as a function of the ex-
pansionforthefourresimulations.Theradiiforthetwores-
imulations starting from 2lpt initial conditions at redshifts
399 and 599 are extremely close. The simulations starting
from Zel’dovich initial conditions show a significantly de-
layed collapse, with the lower start redshift showing the
largest delay.This is in accordance with what would beex-
pected from the predictions of the spherical collapse model
as shown in Subsection 2.2.
Thecollapse occurring intheR5simulation isonly ap-
proximately spherical as evidenced by Figures 3 and 4 in
Gao et al.(2005).Itisnonethelessinterestingtoseewhether
theresultsfrom thesphericaltop-hatmodeldoagreequan-
titatively as well as qualitatively. To do this we first esti-
mated an expansion factor of collapse for each of the five
2lpt z = 599 curves, by fitting a spherical top-hat to each
s
curve,andtakingavaluefortheexpansionfactorofcollapse
from thecorresponding top-hatmodel. Wechose thesepar-
ticular curves because the expected timing error for these
2lpt initial conditions should be very small, compared to
the simulations started from Zel’dovich initial conditions,
andsoshouldnotbiastheestimateofthecollapseexpansion
factor significantly. In Figure 4 the same radii are plotted
Figure 4. Same as Figure 3 except that the expansion factor
as in Figure 3. However the values of the expansion factors
plottedforeachofthefourcurveshasbeencorrectedtoaccount
for the time delay expected in the spherical collapse model and alongthex-axisare‘corrected’usingtheequations(13)and
givenbyequations (13)and(14).Thedensityinsideeachradius (14) to cancel out the expected timing errors for both the
exceeds ten times the critical density to the right of the black 2lpt and Zel’dovich initial conditions. The timing error is
dots.Thecorrectionisexpected toworkwellforasphericalcol- assumed constant along each curve. Strictly speaking this
lapse above an overdensity of around ten (see main text). Even correction applies at the moment the spherical top-hat so-
through the collapse is not perfectly spherical applying the cor- lution collapses to a point. However inspection of Figure 1,
rection does greatly improve the match between the four curves
togetherwith thefact thattop-hatsolution issymmetricin
particularly to the right of the black dots. The spacing of the
time,foranyinitialcondition,aboutthepointofmaximum
outputsisinsufficienttoresolvefeaturesarounda=0.018where
expansion (for an Einstein-de Sitter universe, an excellent
shellcrossingbecomesimportantinthelowertwocurves.
approximation at redshift ∼50.), shows that thetiming er-
ror mostly arises around the time of maximum expansion
the most complex that have been created by the Virgo and therefore the correction at the time of collapse, should
resimulation code, requiring seven Fourier grids. The very be applicable once the top-hat has collapsed by more than
large dynamic range also required that the simulation of ∼0.85 of theradius of maximum expansion. For thespher-
R5 itself was carried out with isolated boundary conditions ical top-hat this corresponds to the density exceeding the
following a spherical region of radius 1.25h−1Mpc which critical density by a factor of about ten, and this overden-
is much smaller than the parent simulation of side length sity is marked on each curve in Figure 4 by a black dot.
479h−1Mpc. We adopted the same isolated region for the To the left of the dots the density is less than ten so the
new simulations. correction is expected to be too strong, which it is. To the
The use of isolated boundary conditions gives the iso- right of the dots the shift works rather well with the lines
lated region as a whole a velocity which depends both on showingmuchbetteragreementthaninFigure3.Thespac-
the starting redshift and initial condition type. In order to ing of the outputs is too crude to follow the collapse of the
compare the different versions of R5 at particular redshifts innermost two radii very accurately, so it is not surprising
it is first necessary to establish a reference point or centre theagreementarounda=0.018islessgoodthanelsewhere.
from which to measure. Using a position determined from In conclusion, the results using 2lpt initial conditions
the main halo is one potential approach, but it was found are much less sensitive to the start redshift and it is possi-
in Gao et al. (2005) that using a group finder like friends- ble to start the simulation later. The formation of the halo
of-friends on this rather extreme region results in groups is delayed using Zel’dovich initial conditions and the mag-
thatare highlyirregular in shape,and arethereforeunsuit- nitude the delay agrees with what one would expect for a
able. Instead we have found a centre for each simulation spherical top-hat collapse.
output using the shrinking sphere algorithm described in
Power et al. (2003) starting with a sphere that encloses the
4.2 Aquarius halos from 2lpt initial conditions
entireisolatedregion.Havingfoundacentre,wehavegrown
spheresfrom thispoint tofindtheradiienclosing masses of Asafinalexample,welookattheeffectofusing2lptresim-
(104,4×104,16×104,64×104,256×104)h−1M⊙.Wechecked, ulationinitialconditionsontheredshiftzerointernalstruc-
foroneoutputredshift,bymakingdotplotsthatthiscentre tureofoneoftheAquariushalos(Springelet al.2008).The
can be identified visually as being in the same recognisable start redshift of the Aquarius halos was 127, so even for
lump. aspherical top-hat collapse thetiming error given byequa-
(cid:13)c 0000RAS,MNRAS000,000–000