Linear perturbative theory of the discrete cosmological N-body problem

5 MB·English

by  B. Marcos

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].