Linear perturbative theory of the discrete cosmological N-body problem B. Marcos, T. Baertschiger Dipartimento di Fisica, Universit`a “La Sapienza”, P.le A. Moro 2, I-00185 Rome, Italy, & ISC-CNR, Via dei Taurini 19, I-00185 Rome, Italy. M. Joyce Laboratoire de Physique Nucl´eaire et de Hautes Energies, Universit´e Pierre et Marie Curie-Paris 6, UMR 7585, Paris, F-75005 France. A. Gabrielli ISC-CNR, Via dei Taurini 19, I-00185 Rome, Italy, & SMC-INFM, Dipartimento di Fisica, Universit`a “La Sapienza”, P.le A. Moro 2, I-00185 Rome, Italy. 6 F. Sylos Labini 0 “E. Fermi” Center, Via Panisperna 89 A, Compendio del Viminale, I-00184 Rome, Italy, 0 & ISC-CNR, Via dei Taurini 19, I-00185 Rome, Italy. 2 n Abstract a J We present a perturbative treatment of the evolution under their mutual self-gravity of particles 0 displaced off an infinite perfect lattice, both for a static space and for a homogeneously expanding 2 spaceasincosmologicalN-bodysimulations. Thetreatment,analogoustothatofperturbationstoa crystalinsolidstatephysics,canbeseenasadiscrete(i.e. particle)generalizationoftheperturbative 1 solution in the Lagrangian formalism of a self-gravitating fluid. Working to linear order, we show v explicitlythatthisfluidevolutionisrecoveredinthelimitthattheinitialperturbationsarerestricted 9 to modes of wavelength much larger than the lattice spacing. The full spectrum of eigenvalues of 7 the simple cubic lattice contains both oscillatory modes and unstable modes which grow slightly 4 faster than in thefluid limit. A detailed comparison of ourperturbativetreatment, at linear order, 1 with full numerical simulations is presented, for two very different classes of initial perturbation 0 spectra. Wefindthattherangeofvalidityissimilartothatoftheperturbativefluidapproximation 6 (i.e. uptocloseto“shell-crossing”), butthattheaccuracyintracingtheevolutionissuperior. The 0 formalism provides a powerful tool to systematically calculate discreteness effects at early times in / h cosmological N-bodysimulations. p - PACSnumbers: 98.80.-k,05.70.-a,02.50.-r,05.40.-a o r t s I. INTRODUCTION which solve numerically for the evolution of a system of a N particles interacting purely through gravity, with a : v softening at very small scales. The number of particles Xi The standard paradigm for formation of large scale N in the very largest current simulations [8] is 1010, ∼ structure in the universe is based on the growth of many more than two decades ago, but still many orders r a small initial density fluctuations in a homogeneous and of magnitude fewer than the number of realdark matter isotropic medium (see e.g. [1]). In the currently most particles ( 1080 in a comparable volume for a typical ∼ popular cosmologicalmodels, a dominant fraction (more candidate). The question inevitably arises of the accu- than80%) ofthe clusteringmatter inthe universeisas- racywithwhichthese“macro-particles”tracethedesired sumedtobeintheformofmicroscopicparticleswhichin- correlation properties of the theoretical models. This is teractessentiallyonlybytheirself-gravity. Atthemacro- theproblemofdiscretenessincosmologicalN-bodysimu- scopicscalesofinterestincosmologytheevolutionofthe lations. Itisanissuewhichisofconsiderableimportance distribution ofthis matter is then very welldescribed by as cosmology requires ever more precise predictions for the Vlasov or “collisionless Boltzmann” equations cou- its models for comparison with observations. pled with the Poisson equation (see e.g. [2]). A full Up to now the primary approach to the study of dis- solution, either analytical or numerical, of these equa- creteness in N-body simulations has been through nu- tions starting from appropriate initial conditions (IC) is merical studies of convergence (see e.g. [9, 10]), i.e., not feasible. There are,on the one hand, various pertur- one changes the number of particles in a simulation and bative approaches to their solution (for reviews see e.g. studies the stability of the measured quantities. Where [3, 4]), which allow one to understand the evolution in results seem fairly stable, they are assumed to have con- somelimited range(essentially ofsmallto moderateam- verged to the continuum limit. While this is a coherent plitude fluctuations). On the other hand, there are cos- approach, it is far from conclusive as, beyond the range mological N-body simulations (for reviews see [5, 6, 7]), of perturbation theory, we have no theoretical “bench- 2 marks” to compare with. Nor is there any systematic on a lattice interacting by a pure unscreened Coulomb theory of discreteness effects, e.g., we have no theoreti- force. At linear order (harmonic analysis) one simply cal knowledge of the N dependence of the convergence. solves a 3N 3N eigenvalue problem to determine the × GiventhattypicallyN isvariedoveraverymodestrange eigenmodes and eigenvalues of the displacements off the (typically one or two orders of magnitude) compared to crystal. This can be done at very low cost in computa- that separatingthe simulationfromthe model (typically tionalresourcesbecause ofthe symmetries ofthe lattice. 70 orders of magnitude) there is much room for error. Stable (i.e. dynamically oscillating) modes in one prob- Different mechanisms by which discreteness effects lem become unstable (growing and decaying) modes in may make the evolution of N-body simulations different the other problem, and vice versa. One consequence of to that of the fluid limit have been discussed in the lit- this which we will discuss briefly here, and more exten- erature. A very basic consideration is that of the dis- sively in [23], is that, for what concerns discreteness in cretenesseffectsintroducedalreadyintheIC,beforeany N-body simulations, there are qualitatively different fea- dynamical evolution. Indeed there are necessarily dis- tures on the simple cubic (sc) lattice compared to the crepancies between the correlation properties of the dis- body centered cubic (bcc) or face centered cubic (fcc) cretised IC and those of the input theoretical model, as lattice. there are intrinsic fluctuations associated with the par- Acrucialstepinouranalysisisevidentlythe compari- ticles themselves. For analysis and discussion of these sonofoursolutionsfortheevolutionwiththoseobtained effectssee,e.g.,[11,12,13,14,15,16]. Whatisprobably fromthetreatmentofthecontinuousself-gravitatingsys- the mostobvious effect of discreteness,andcertainly the tem. This latter problem can also be solved in a limited one most emphasised in the literature, is two body colli- range with a perturbative treatment of the equations of sionality: pairs of particles can have strong interactions a self-gravitating fluid, which are obtained by trunca- with one another, which is an effect absent in the colli- tion of the collisionless Boltzmann equation. Given that sionlesslimit. Foranalysisanddiscussionoftheseeffects ourperturbationschemeworkswiththedisplacementsof see, e.g., [17, 18, 19, 20, 21]. the particles, one might anticipate that the appropriate In this paper we present an approach which allows in perturbationschemetocomparewithisthatgiveninthe principle a systematic understanding of discreteness ef- Lagrangianformulationofthesefluidequations,inwhich fectsbetweenthesetworegimes,intheevolutionfromthe the evolution is described in terms of the trajectories of IC up to the time when two body collisions start to oc- fluid elements [29]. We show explicitly, at linear order, cur. We do so by developing a perturbative solution to that this is the case: taking the limit in which the per- the fully discrete cosmological N-body problem, which turbations to the lattice are of wavelengths much larger is valid in this regime. This essentially analytic solu- than the lattice spacing ℓ the evolution described by our tion can be compared to the analogous fluid (N ) schememapspreciselyontothatatthesameorderinthe → ∞ solution,andonecanunderstandexhaustivelythemodi- Lagrangian description of the fluid. The Zeldovich ap- ficationsintroduced,atagiventime andlengthscale,by proximation,whichissimplytheasymptoticformofthis thefinitenessofN. Whiletheusefulnessoftheapproach solution, can then be understood in very simple analogy is restrictedto the regimeof validity ofthis perturbative with the long-wavelengthcoherent “plasma oscillations” approach, we can gain considerable insights into the ef- in a unscreened charged plasma. fects of discretenessandhowthey introduce error. Some The paper is organized as follows. In the next section oftheessentialresultshavealreadybeenbrieflyreported wedescribetheperturbativetreatmentforperturbations in [22]. In this paper we describe in much greater detail off a perfect lattice, for the specific case of gravity. We the perturbative method used to describe the evolution, workfirstly,forsimplicity,inastaticEuclideanuniverse, and evaluate its regime of validity by extensive compari- givingtheexplicitexpressionsfortheevolutionfromgen- son with numerical simulations. In a forthcoming paper eral IC (i.e. any perturbation from the lattice). In the [23] we will discuss the application of this method to followingsection we explicitly solvethe evolutionfor the the study of discreteness effects in N-body simulations, case of a simple cubic lattice, and discuss in detail the providing precise quantifications of these effects in the structure of the spectrum of eigenvalues. In Sect. IV we regime in which our treatment is valid. generalize our results to the case of an expanding uni- The perturbative scheme we employ is one which is verse, and then show explicitly the recovery of the fluid well known and standard in solid state physics, as it is limit given by the solution at the same order of the La- that used in the treatment of perturbations of a crystal grangianformulationoftheequationsofaself-gravitating aboutthe localorglobalminimum ofits internalenergy. fluid. In the next sectionwe presenta comparisonof the Indeed the class of cosmological N-body simulations we evolution described by our approximation with that of considerarethosewhichstartbymakingverysmallper- fullnumericalN-bodysimulations. Weconsiderbothun- turbationstoparticlesinitiallyplacedonaperfectsimple correlatedinitial perturbations(a “shuffledlattice”)and lattice [6, 9, 10, 24, 25, 26]. Up to an overall change in a set of highly correlated perturbations (with a power sign,ourperturbativeschemeispreciselythatonewould spectrumofdensityfluctuations k−2),likethatincur- ∼ use for the analysis of the extensively studied Coulomb rent cosmological models. We verify that the agreement latticeorWignercrystal(seee.g. [27,28]),ofN particles isverygood,andthattheevolutionistracedwithconsid- 3 erablybetteraccuracythanbythefluidlimitatthesame Taking the infinite volume limit V , neither the → ∞ order. Both approximations break down when particles gravitational potential (1), nor the gravitational force start to approach one another (i.e. “shell-crossing” in (3), arewelldefined. In the firstcasethe resultdiverges, fluidlanguage),whichweparametrisethroughanappro- while in the second it may be finite or infinite, but its priate statistical measure. In the final section we sum- value depends on how the limit is taken1. marise our results and discuss various further develop- In Euclidean spacetime this behaviour in the infinite ments of the method presentedhere which could be pur- volume limit may be regulated by the introduction of a sued, notably the extension of the perturbation scheme negative background — the so-called Jeans swindle (see tohigherthanlinearorder. Thismaythrowfurtherlight e.g. [2, 30]) — so that the potential is defined as on how insights about discreteness gainedusing this for- 1 malism, which will be discussed at length in [23], extend φ(r) = G lim m (V,r′) into the regime of highly non-linear evolution. We also − V→∞ r r′ V (cid:20) r′6=r| − | discuss briefly the possible interest of solving the cos- X 1 mological N-body problem on the bcc lattice, as well as ρ d3r′ (V,r′) . (4) the possible extension of our method to IC generated on − 0 ZR3 |r−r′|V (cid:21) “glassy” configurations. This modifies the usual Poisson equation to 2φ(r)=4πG(ρ(r) ρ ). (5) 0 II. LINEARIZATION OF GRAVITY ON A ∇ − PERTURBED LATTICE The expression (4) is well defined2, provided (i) that the limit V is taken in a physically reasonable → ∞ In this section we start by defining and studying some way3, and (ii) that the fluctuations in the system are generalpropertiesofthegravitationalpotentialandforce sufficiently rapidly decaying at large scales4. In the cos- ofaninfinite systemofpointparticles. We thenconsider mologicalcontextthis negativebackgroundappears nat- the particular case of a perturbed infinite lattice in a urally as a consequence of the expansion of the universe static Euclidean space, the generalization to an expand- (see Sect. IV). ing universe being given in Sect. IV. Since the force is The simulations of self-gravitating systems we are in- zero in the unperturbed lattice, the dominant contribu- terested in are performed using a finite cubic simulation tion to the force in the perturbed case is linear in the box of side L and volume VB = L3, subject to peri- relative displacements of the particles. In the last sub- odicboundaryconditions. Theforceonaparticleisthus section, we consider the equations of motion resulting computednotonlyfromallthe otherparticlesinside the from this linearized force. simulation box, but also from all the copies of the par- ticles contained in the “replicas”. The reason for using these boundary conditions is that they introduce the in- A. Definition of the force and the potential evitable finite size effects without breaking translational invariance: every particle can be considered to be at the centre of the finite box and therefore sees the boundary Letusconsidercarefullyfirstthedefinitionofthegrav- in the same way. The infinite system we consider is thus itational force in an infinite system of point particles of an infinite number of replicas of a finite cubic box, with equal mass m. We will assume that this system (either a negative background as described above to make the stochasticordeterministic)ischaracterizedbyawellde- force well defined5. In this case the gravitationalpoten- fined mean number density n > 0, and mass density 0 tial may be written as ρ = mn . The gravitational potential of a particle, per 0 0 unitmass,atr,due totheparticlesinafinite volumeV, φ(r)= lim [φ (r)+φ (r)], (6) is: V→∞ b p 1 φ(r)= Gm (V,r′), (1) − r r′ V r′6=r| − | X 1 F(r)isaconditionallyconvergent series. where the sum is over all the particles contained in the 2 Foramoredetaileddiscussionofthegravitationalforceininfinite system,and (V,r)isthewindowfunctionforthevolume systemsseealso[31]. V, i.e., V 3 E.g.,takingtheinfinitevolumelimitincompact sets. 4 If P(k) is the power spectrum of density fluctuations, itis sim- 1 if r V, ple to show, using the modified Poisson equation Eq. (5), that (V,r)= ∈ (2) convergence ofthefluctuations inthegravitationalpotential re- V 0, otherwise. (cid:26) quires limk→0knP(k) = 0 for n > 1. For finite fluctuations in theforceonerequiresn> 1. The force per unit of mass (i.e. the acceleration), due 5 Note also that, because th−e system is just a lattice when con- to these same particles, is given by the gradient of the sidered at scales larger than the box size, the fluctuations are potential: alwayssufficientlysuppressedatlargescalessothatthegravita- tional forceiswell defined. Thus any possibledivergence inthe F(r)= φ(r). (3) fluctuations offorcewillberegulatedbytheboxsizeL. −∇ 4 where Note that the contribution from the background (10) is identically zero if one takes a window function with in- 1 φ (r)=Gρ d3r′ (V,r′) (7) version symmetry in r (e.g. a sphere or cube centred on b 0 R3 r r′ V r). Z | − | is the contribution from the background,and ∗ (V,r′+nL) φ (r)= Gm V (8) p − r r′ nL B. Linearization of the gravitational force n,r′ | − − | X thecontributionfromtheparticles. Herethesumoverr′ We consider the infinite lattice generatedby the repli- is restricted to the particles in the box, while the other cation of a sc lattice of volume VB of side L with N ele- sum, over the three integers n (i.e. over the images of ments, i.e., whose lattice vectors are R = (m1,m2,m3)ℓ r′),hasa“*”toindicatethatthetermr′ =risexcluded withmi [0,N1/3 1] N andℓ=L/N1/3 is the lattice when n=0. spacing6∈. This lat−tice∩(with a particle at each site) is The gravitationalforce is: now perturbed by applying displacements u(R) to each particle R, so thatthe new positions of the particles can F(r)= lim [Fb(r)+Fp(r)], (9) be written as V→∞ where r=R+u(R). (12) r r′ F (r)=Gρ d3r′ − (V,r′) (10) b 0 R3 r r′ 3V The “particle” term in the gravitational force [i.e. Z | − | Eq. (11)] can then be expanded order by order in Taylor and series about its value in the unperturbed lattice. At lin- F (r)= Gm ∗ r−r′−nL (V,r′+nL). (11) oeabrtaoirnder in the relative displacements u(R)−u(R′) we 6 Thepgenera−lizationno,rf′a|lrl t−her′ca−lcunlaLt|i3onVs presented here to any BravaislatticeissXtraightforward(seee.g. [32]). ∗ R R′+nL u(R) u(R′) [u(R) u(R′)] [R R′+nL] F (r) = Gm − + − 3 − · − (R R′+nL) p − R R′+nL3 R R′+nL3 − R R′+nL5 − n,R′(cid:26)| − | | − | | − | (cid:27) X (V,R′+nL). (13) ×V The first term in this sum The convergence criterion for each term of (13) is ∗ R R′+nL R R′ > u(R) u(R′). (15) Gm − (V,R′+nL) (14) | − | | − | − R R′+nL3V n,R′ | − | Note that the validity of the power expansion does not X depend on the displacement of the particle R on which has the poor infinite volume behaviour which is regu- we compute the force, but on relative displacements of lated, as discussed above, by the contribution coming the particles at the position R and R′. Under the ac- fromthe backgroundEq.(10). The totallinearizedforce tion of the gravitational interaction, the displacements is then also well defined, and given by the infinite vol- u(R) will typically grow so that the condition Eq. (15) ume limit of Eq. (13) summed with Eq. (10). In the is violated after some time. However when some pairs case that we choose to calculate using the infinite vol- of particles no longer satisfy condition (15), it may nev- ume limit of a volume V with inversion symmetry in r ertheless continue to apply for the rest of the particles (i.e. the displaced position of the particle), the full lin- and (13) may remain a sufficiently good approximation earized force is thus given by Eq. (13). If, however, we to the force. In order to have a precise characterization choose to sum in a volume with inversion symmetry in of the regime of validity of the approximationapplied to the lattice site R, it is simple to show that Eq. (14) is follow the dynamical evolution of a perturbed lattice, it identically zero. The backgroundterm then contributes, is necessary to compare the results with those obtained with the sum [(10)+(14)] remaining invariant. from evolution under full gravity. We will perform such 5 a comparison in Sect. V using N-body simulations. It is convenient to write the linearized force just dis- cussed in terms of the so-called dynamical matrix (R) D (see e.g. [32, 33]): F(r)= (R R′)u(R′). (16) D − R′ X This matrix has the following properties: it is a com- plete symmetric operator,i.e., (R)= ( R) with µν νµ D D − inversion symmetry, i.e., (R) = ( R). Further, µν µν D D − since the same displacement applied to all the particles produces no net force, we have (R) = 0. For R µν FIG.1: Computation ofthediagonal termsof thedynamical D anypairinteractionpotentialv(r) itis straighforwardto matrix at R=0. P show that it can be written as [32, 33] (R=0)=∂ ∂ w(R) (17a) C. Equations of motion in a static Euclidean µν µ ν D 6 universe (R=0)= ∂ ∂ w(R′) (17b) µν µ ν D − RX′6=0 In this section we derive the equations of motion of the particles in the linear approximation,and then solve where them. We treat first a static Euclidean space, giving the ∂2w(r) generalization to a cosmological expanding universe in ∂ ∂ w(r )= (18) µ ν 0 ∂r ∂r Sect. IV. (cid:20) µ ν(cid:21)r=r0 Using Newton’s second law and Eqs. (12) and (16) we and w(r) is the periodic function defined as can write the equation of motion of the particles as: u¨(R,t)= (R R′)u(R′,t), (22) w(r)= v(r+nL), (19) D − R′ n X X where the double dot denotes a double derivative with i.e.,thepotentialduetoasingleparticleandallitscopies. respect to time. The expression (22) is a system of vec- For gravity we have v(r) = Gm/r and Eq. (19) is im- torial coupled second order differential equations which − plicitlyunderstoodtoberegulatedasdiscussedatlength canbe reducedtoaneigenvalueproblem,usingstandard abovebythe additionofauniformnegativebackground. techniques. From Bloch’s theorem [32] it follows that We will describe below, and in App. A, how we use the Eq. (22) can be diagonalized by the following combina- well-knownEwaldsummationtechniquetoexplicitlyper- tion of plane waves: form this sum. Equation (17b) gives the force on a particle, at first u(R,t)= 1 u˜(k,t)eik·R, (23) order in the displacements, when it is displaced and all N k the others remain unperturbed (see Fig. 1). For gravity X it is straightforward[31] to show that where the sum over k is restricted to the first Brillouin zone, i.e., for a sc lattice to 4π (0)= Gρ δ , (20) 2π Dµν 3 0 µν k= n, (24) L i.e.,thelinearizedforcefs(r)onaparticledueonlytoits withn=(n ,n ,n )suchthatn [ N/2,N/2[ Z. We 1 2 3 i owndisplacementuwithrespecttotherestofthelattice denote by u˜(k,t) the Fourier tran∈sfo−rm of u(R,t∩): is u˜(k,t)= u(R,t)e−ik·R, (25) 4π fs(r)= Gρ0u(R). (21) R 3 X where the sum is restricted to the simulation box (i.e. The simplestwaytoderivethis resultis bysumming the without considering the replicas). Inserting Eq. (23) in force in spheres centred on the unperturbed position of Eq. (22), we obtain for each k: the displaced particle. In this case it is straighforward to show,by symmetry, that the linearizeddirect particle u¨˜(k,t)= ˜(k)u(k,t), (26) contributionEq.(13)iszeroandthefullforceisgivenby D the background term Eq. (10). The result follows then where ˜(k) is the FT of (R), defined analogously to D D simply fromGauss’law whichgivesthat the forcecomes (25). From the properties of (R) given above, it fol- D only from the region inside the sphere shown in Fig. 1. lows that ˜(k) is a real and symmetric operator which D 6 satisfies7 where the matrix elements of the “evolution operators” and are 4π P Q lim ˜ (k)= Gρ δ . (27) k→0Dµν 3 0 µν 3 (k,t)= U (k,t)(ˆe (k)) (ˆe (k)) , (35a) µν n n µ n ν P We can now solve Eq. (26) by diagonalizing the 3 3 n=1 matrix ˜(k). Foreachk,thisdeterminesthreeorthon×or- X3 D maleigenvectorsˆe (k) withthree associatedeigenvalues (k,t)= V (k,t)(ˆe (k)) (ˆe (k)) . (35b) n µν n n µ n ν ωn2(k) (n=1,2,3) satisfying the eigenvalue equation: Q nX=1 ˜(k)ˆe (k)=ω2(k)ˆe (k). (28) Theoperator thusevolvestheinitialdisplacementfield D n n n and the iniPtial velocity field. Q Wecandecomposeeachmodeu˜(k,t)inthebasis ˆe (k) n { } as III. DETERMINATION AND ANALYSIS OF 3 THE SPECTRUM OF EIGENVALUES OF ˜(k) D u˜(k,t)= ˆe (k)f (k,t). (29) n n nX=1 In this section we describe the determination of the eigenvectors and spectrum of eigenvalues of the dynam- UsingEqs.(26),(28)and(29)wegetthefollowingequa- ical matrix for gravity. We then discuss the physical tion for the coefficients f (k,t): n meaning of the results, notably identifying how the fluid limitisobtainedandhowcorrectionstothislimitmaybe f¨ (k,t)=ω2(k)f (k,t). (30) n n n calculated. In this discussion we will use extensively the strict analogy between the case we are treating and the Depending onthe signofω2(k), weobtaintwoclassesof n Coulomb lattice, or Wigner cystal, studied in condensed solutionsU (k,t)andV (k,t). Wechoosethem,without n n matter physics (see e.g. [28]). This is a system of posi- any loss of generality, satisfying tively charged particles embedded in a negative neutral- izingbackground. Theparticlesinteractwitharepulsive U (k,t )=1, U˙ (k,t )=0, (31a) n 0 n 0 1/r potential instead of the attractive 1/r potential of V (k,t )=0, V˙ (k,t )=1. (31b) Newtoniangravity. Thusallourresults−aremappedonto n 0 n 0 those for the corresponding Coulomb lattice by making The function U (k,t) is associated with initial displace- the formalsubstitution Gm2 e2, where e is the elec- n ments and V (k,t) with initial velocities. If ω2(k) 0 tronic charge8. →− n n ≥ then Un(k,t)=cosh(ωn(k)(t t0)), (32a) A. Numerical computation of the spectrum of ˜(k) − D V (k,t)=sinh(ω (k)(t t ))/ω (k). (32b) n n 0 n − The spectrum of the matrix ˜(k) must be computed If ω2(k)<0 numerically. The matrix (R)Dis constructed using the n D Ewald sum method [32, 33, 34, 35] to speed up the con- U (k,t)=cos( ω2(k)(t t )), (33a) vergenceofthesum. Wecontinuetoworkhereexplicitly, n | n | − 0 as above, with a sc lattice of side L, with lattice spacing Vn(k,t)=sin(p|ωn2(k)|(t−t0))/ |ωn2(k)|. (33b) ℓ and N elements9. To determine the dynamical matrix p p we use the Ewald method to evaluate w(r) as given in Whereas the modes (32) with positive eigenvalues cause Eq.(19),splittingitintotwopiecesusinganappropriate anexponentialgrowthofperturbationinthesystem,the damping function : modes (33) with negative eigenvalues leads to oscilla- C tions. The evolution of the displacement field from any w(r)= v(r+nL) (r+nL,α) initial state u(R,t0) is then given by the transformation C | | n X (36) + v(r+nL)[1 (r+nL,α)], 1 −C | | u(R,t)= (k,t)u˜(k,t0)+ (k,t)u˜˙(k,t0) eik·R Xn N P Q k Xh i (34) 8 Thepotentialwehaveusedhereforgravityhasbeendefinedper unit mass, i.e., in our notation v(r) = e2/mr for the Coulomb lattice. 7 But note that ˜µν(k = 0) = R µν(R) = 0, i.e., ˜(k) is 9 Thegeneralizationtoaparallelepipedbox,andtootherBravais discontinuousatDk=0. D D lattices,isstraightforward(seee.g. [32]). P 7 whereαisaarbitrary“dampingparameter”ofwhichthe 1.2 resultis independant. The function (r,α) is chosento [1,0,0] be equal to unity at r = 0 and rapidCly| d|ecaying to zero 1 [1,1,0] [1,1,1] as r goes to infinity. The first sum is then evaluated in | | 0.8 real space and the second one in Fourier space, making use of the Parseval theorem [36], being chosen so that the secondterminEq.(36)isanaClyticatr=0andthus ε(k)n0.6 rapidly convergent in Fourier space. A common choice 0.4 for a 1/r pair potential is 0.2 (r,α)=erfc(αr). (37) C | | | | 0 The expression for the function w is then: 0 0.5 1 1.5 k/k w(r)=w(r)(r)+w(k)(r). (38) N FIG.2: Spectrumof eigenvaluesforsimplecubiclattice with In the gravitationalcase 163 particles. The lines correspond to chosen directions in k space. 1 w(r)(r)= Gm erfc(αr+nL), (39a) − r+nL | | n | | X 4π 1 k2 B. Analysis of the spectrum of eigenvalues in a w(k)(r)= Gm exp | | cos[k r], simple cubic lattice − V k2 −4α2 · B k=6 0| | (cid:18) (cid:19) X (39b) We now describe the spectrum of eigenvalues of the dynamical matrix (R) for a sc lattice. As we have dis- where V is the volume of the box and the wavevectors D B cussedintheintroduction,thisisthelatticewhichisused k are as in Eq. (24), but with n ranging over all triple verywidelyinN-bodysimulationsofstructureformation integers (i.e. not restricted to the first Brillouin zone). in cosmology. There is no k = 0 term in the sum (39) because of the In Fig. 2 we plot the spectrum of a sc lattice, for N = presenceofthenegativebackground: whensummedover 163, obtained with the method outlined in the previous all the particles, this term is equal to subsection. We show the normalized eigenvalues kli→m0φ˜0(k)=−kli→m04πkG2ρ0, (40) εn(k)= 4ωπn2G(kρ) (43) 0 i.e.,thek=0modeofthepotential(calculatedfromthe asafunctionofthemodulusofthekvectors,normalized Poisson equation in Fourier space) which is cancelled by to the Nyquist frequency k = π/ℓ. With this normal- N the contribution coming from the negative background. isation the spectrum remains substantially the same as The Ewald sum for the dynamical matrix can then be we increase the number of particles: the only change is calculated directly using Eq. (17) and (39). The result, thattheeigenvaluesbecomedenserintheplot,fillingout as in Eq. (38), is divided in two parts: theapproximatefunctionalbehaviourswithmorepoints. Forourdiscussionherethereisnointerestinconsidering (R)= (r)(R)+ (k)(R), (41) a greater number of points than that we have chosen. D D D For each vector k there are three eigenvalues ω2(k), n for which the explicit expressions are given in App. A. n=1,2,3. Eachfamily of eigenvalues(i.e. with samen) For the results quoted here we have taken α = 2/L defines a surface, corresponding to the three branches of [37]. Using this numerical value of α, it is sufficient to thefrequency-wavevectordispersionrelation. Sectionsof sum for these surfaces are plotted for some chosen directions of the vector k in Fig. 2. 6π n 3 k . (42) | |≤ | |≤ L to obtain a well converged determination of the dynam- 1. An expression for ˜(k) and the Kohn sum rule D ical matrix. The diagonalization calculation involves es- sentially N operations (where N is the number of parti- Before proceeding further it is useful to derive some cles). It is perfectly feasible, with modest computer re- important results we will employ much in what follows. sources,toperformthisdiagonalisationforparticlenum- These are well known in the context of the application bers as large as those used in the largestcurrentN-body of this formalism in condensed matter physics (see e.g. simulations. [28]). First of all, we derive an analytical expression for 8 thedynamicalmatrixinFourierspace. Letusdecompose InthecontextoftheCoulomblatticethisisawell-known in Fourier modes the function w(r) defined in Eq. (19) result, the so-called Kohn sum rule. In this case the quantity which appears on the r.h.s. of the sum, instead w(r)= V1B k w˜(k)eik·r, (44) ofrfeq4uπeGnρc0y,. iWs −eωwp2ill=d−isc4uπses2nfu0r/tmherwbheelroewωtphiesstihgenipfilcaasnmcae X of this analogy. where the sum over k is performed over all k space, i.e., We can use these results and the above sum rule to not restricted to the first Brillouin zone and compute — in a different way than in Eqs. (20)–(21) — the R= 0 term of the dynamical matrix (R) (i.e. the w˜(k)= d3rw(r)e−ik·r. (45) term giving the force on a particle, at lineaDr orderin the ZVB relativedisplacements,whenitaloneisperturbedoffthe lattice). Using the Kohn sum rule (53), the trace of the The derivatives of the periodic potential are dynamical matrix is: 1 w (r)= k k w˜(k)eik·r. (46) µν µ ν tr[ (R)]=4πGρ . (54) −V 0 B k D X If the crystal has three equivalent orthogonal directions Using the definition of the dynamical matrix then the diagonal terms of the dynamical matrix will be ˜ (k)= (R)e−ik·R (47) equal. In the case of lattices with special symmetries µν µν D D (like the sc, bcc and fcc) it is simple to show that when R X a single particle is displaced along the direction of an and Eqs. (17) and (46) we obtain: axis, the force acting on it is parallel to the direction of displacement10. This implies that the non-diagonal ˜ (k)= 1 k′k′w˜(k′) eiR·(k′−k) eik′·R termsofthedynamicalmatrixarezero. Wecantherefore Dµν −VB kX′,R µ ν (cid:16) − (cid:17) conclude that (48) 4 wherewecanincludethetermR=0inthesumbecause µν(0)= πGρ0δµν. (55) D 3 it vanishes. Using the orthogonality relation, we have ei(k−k′)·R =N δk′,k+K, (49) 2. The branches of the dispersion relation and the fluid R K limit X X where the k are restricted to the first Brillouin zone and We have noted that the spectrum of eigenvalues has a K are the reciprocal vectors of R satisfying clearbranchstructure. Toidentifythedifferentbranches K=2k m, (50) it is useful to consider the k 0 limit keeping the in- N → terparticle distance ℓ constant. We expect this to corre- with m Z3. Substituting Eq. (49) in (48) we obtain spond to the fluid limit: a plane wave fluctuation eik·r ∈ finally the expression [28]: with k 1/ℓ has a variation scale much larger than the ≪ interparticle distance, and therefore does not “see” the ˜µν(k)= n0kµkνw˜(k) (51) particles. D − From Eq. (51) the limit for k 0 is straightforward n0 [(kµ+Kµ)(kν +Kν)w˜(k+K) KµKνw˜(K)], → − − as the contribution of the sum on the r.h.s. vanishes in KX6=0 this limit11 where n is the number density ofparticles. In the grav- itational0case, the integral (45) cannot be evaluated an- lim ˜µν(k)= n0kˆµkˆνw˜(k). (56) k→0D − alytically. However, neglecting finite size effects, this in- tegral can be computed over the whole space and the Using the eigenvalue equation (28) with Eqs. (51) and periodic potential w(r) is approximated by the interac- (52), it follows that the solutions in the fluid limit are tion pair potential v(r)= Gm/r, so that − 1. one longitudinal eigenvectorpolarized parallelto k w˜(k) v˜(k)= d3rv(r)e−ik·r = 4πGm. (52) with normalized eigenvalue ε1(k→0)=1 and ≃ R3 − k2 Z Using this it is straightforwardto show (see App. B) the following simple result: 10 Thiscanbeexplicitlyshowne.g. usingEq.(A2)(takingthelimit α 0andassumingthatthesumoverthereplicasconverges). 3 11 W→e have assumed that the sum in Eq. (51) is well defined — ωi2(k)=−n0k2w˜(k)=4πGρ0. (53) whichisthecaseforthegravitationalinteraction—sothatitis possibletotakethelimitbeforeperformingthesum. i=1 X 9 2. two transverse eigenvectors polarized in the plane 1 transverse to k with normalized eigenvalues [1,0,0] [1,1,0] ε (k 0)=0. 2,3 [1,1,1] → Asthespectrumofeigenvaluesε (k)isexactlythesame, n upto anoverallnegativemultiplicative constant,tothat 1| of the Coulomb lattice, we adapt the same terminology k)- 0.01 asinthiscontext. Thebranchofeigenvalueswhoseasso- ε(1 | ciated eigenvectors converges to the longitudinal eigen- vector as k 0 is called the optical or longitudinal → branch. The twoother brancheswhoseeigenvectorscon- vergeto the transverseeigenvectorsarecalledthe acous- tic branches. For finite k, the eigenvectors are not ex- 0.0001 actly parallel or perpendicular to kˆ for all k but belong 0.03 0.1 k/k 1 2 N nevertheless to one of the three branches, which define three-dimensionalhyper-surfacesinthefour-dimensional FIG.3: Opticalbranchfordifferentdirectionsofk. Thethick space (ω,k) space. line is proportional to k2. The appearance ofan optical branchin a monoatomic crystal is a characteristic feature of the 1/r interaction potential (at large r). In the case of a more rapidly de- Wecanestimateanalyticallythecorrectionstothislimit caying potential at large scales, i.e., 1/r1+α with α > 0, for small k (i.e. for large wavelengths) by expanding the it becomes a third acoustic branch. In the case of a po- eigenvaluesandeigenvectorsofthefulldynamicalmatrix tential that decays slower at large r, i.e., α < 0, the about k = 0. We note that this corresponds to calcu- optical branch diverges as k 0. The physical inter- lating the difference, at large wavelengths, between the → pretation of the optical branch is that it represents the evolutionoftheperturbedlatticewithafinitenumberof coherent excitation of the whole lattice with respect to particlesandthatofthe fluidlimit. These arethuswhat the background [27]. In a Coulomb crystal, the optical are,inthecontextofcosmologicalsimulations,“discrete- modeisproducedbythelatticemovingagainstthisback- ness effects” introduced by the modelling of the fluid by ground producing a “plasma oscillation”, at the plasma sucha system. We willdiscuss at lengththis application frequency ω defined above. This mode is, as we have of this formalism in [23]. p just seen, purely longitudinal, i.e., the perturbations are WhenexpandingthedynamicalmatrixinTaylorseries parallelto k, while the tranverse modes, i.e., the pertur- about the fluid limit k 0, it is simple to show that for → bations orthogonalto k have zerofrequency. The reason 1/rinteractionsthisseriesisinevenpowersofk,because forthisbehaviouroflongwavelengthdensityfluctuations (R) is real and ˜(k) analytic for k 0 (see [27, 38]). D D → can be easily understood. The density fluctuations are It is therefore possible to write the corrections to the related, in this fluid limit, to the displacements through eigenvalues of the optical mode as: the continuity equation: ω2(k) 4πGρ (1 b (kˆ)k2), (60) δρ u, (57) 1 ≃ 0 − 1 ∼∇· which implies in k space that where the expression for b1(kˆ) can be computed by di- agonalizing ˜(k) expanded up to (k2). The leading δρ˜ k u˜. (58) D O correction to the two acoustic modes may be written ∼ · Thustranversemodesdonotsourcedensityfluctuations, ω2(k) 2πGρ b (kˆ)k2, (61a) andtherefore(by the Poissonequation)they donotpro- 2 ≃ 0 2 duce a force. In the case of gravity,instead of oscillating ω2(k) 2πGρ b (kˆ)k2. (61b) as in a plasma, the longitudinal mode may be amplified 3 ≃ 0 3 orattenuated(depending onthe initialperturbation),in The Kohn sum rule implies that b (kˆ) = (b (kˆ) + 1 2 awaywhichisindependentofk. Aswewilldiscussinde- b (kˆ))/2. InFig.3weshowtheopticalbranch,invarious 3 tailbelow,thisisjustthewellknownlinearamplification different chosen directions. The approximation with the of density fluctuations in a self-gravitating fluid. leading term in the Taylor expansion is very good up to the Nyquist frequency. InFig.4weshowhowtheanisotropyoftheeigenvalues 3. Corrections to the fluid limit increasesasthemodulusofthewavevectorincreases(i.e. whenwelookatsmallerspatialscales). Weplot,forthree We have just seen that the fluid limit is obtained by rangesofvaluesofthemodulusofk,thevalueofthenor- taking the dynamical matrix as malized eigenvalues as a function of the angle θ between 4πGρ k and the axis that forms a minimal angle with it. As θ D˜(k)= k2 0kµkν. (59) increases (i.e. as cosθ decreaseswith 0<θ <π/2) there 10 1.02 0<k/k <0.1 N 0.1<k/k <0.2 N 0.2<k/k <0.3 N 1 k) 2ε( 0.98 (i) (ii) 0.96 FIG. 5: Schematic representation of (i) a mode collapsing faster than fluid limit and (ii) an oscillating mode. 0.6 0.7 0.8 0.9 cos(θ) FIG. 4: Variation of the value of the eigenvalues for various behaviour qualitatively different to that generically ex- rangesasafunctionof thecosineoftheanglebetween kand pectedbasedontheanalysisofthefluidlimit. Translated the axes of the lattice which forms a minimal angle with it. We see that the effects of anisotropy are strongest for the to the analagousCoulombsystem, the resultmeans that short-wavelength modes, and decrease as we go towards the a sc lattice is, in this case, unstable (as there are grow- fluid limit. ing modes). While this result has not apparently been shown in the literature, it is not an unexpected result in this context. It has been established [39, 40] that for is a clear trend of decrease in the eigenvalue, in each of the (classical) Coulomb lattice that the ground state is the three cases. The difference as a function of orienta- the bcc lattice. It has a lower binding energy than the tion of the vector k is, however, much more marked for fcc lattice, which in turn is a lower energy configuration larger k, i.e., at scales closer to the Nyquist frequency. than the sc lattice. Our result implies that the latter is Thisisnotunexpected: the effectsofanistropy(whichis not only a higher energy state, but that it is strictly un- completely absent in the fluid limit, in which the eigen- stable. Indeed we note that the specific modes we have valuesareindependentoftheorientationk)arenaturally considered above describe a “sliding” of adjacent places strongest for the short wavelength modes. inansc lattice whichdeformittowardsthe lowerenergy configuration represented by the fcc lattice. 4. Oscillatory modes The spectrum of the sc lattice Fig. 2 includes some modes [e.g. for k = (k ,0,0)] with eigenvalues on the IV. GENERALIZATION TO AN EXPANDING x optical branch larger than the fluid limit. For exam- UNIVERSE ple, this is the case for modes with initial displacement u(r,0) xˆexp(ik x), shown in Fig. 5-(i). Adjacent ∝ x In the previous section, we have described the gravi- planes collapse towards one another, faster than in the tational evolution of a perturbed lattice in a static Eu- fluid limit. The Kohn sum rule Eq. (53) states that the clidean universe. In the cosmological context, density sum of the three eigenvalues ω2(k) is equal to 4πGρ . n 0 fluctuations are a perturbation around an homogeneous Therefore, the existence of modes collapsing faster than andisotropicFriedmann-Robertson-Walker(FRW) solu- the fluid limit implies that there are other modes with tion of Einstein’s field equations of general relativity. In negative eigenvalues ω2(k), i.e., which oscillate. This n cosmologicalN-body simulations,since the regionsstud- is the case, e.g., of the mode with initial displacement iedaresmallerthanthe Hubble radiusandthe velocities u(r,0) yˆexp(ik x), shown in the Fig. 5-(ii). In this ∼ x are non-relativistic, one considers the limit in which the case, contiguous planes oscillate as indicated in the fig- equations of motion of the particles are strictly Newto- ure. We will study these modes in greater detail using nian in physicalcoordinates r [1]. These coordinates are numerical simulation in Sect. VC. related to the comoving coordinates x of the FRW solu- The existence of oscillating modes in a perturbed and tion by cold purely self-gravitating system (i.e. without any additional interaction or velocity dispersion giving rise to a restoring pressure12) is an unexpected curiosity, a r(t)=a(t)x(t), (62) 12 If there is a non negligible velocity dispersion, it known that fluctuations atscalessmallerthantheJeanslengthoscillate[2].