Table Of ContentMon.Not.R.Astron.Soc.000,1–??(2014) Printed29August2016 (MNLATEXstylefilev2.2)
The fully non-linear post-Friedmann frame-dragging vector
potential:
Magnitude and time evolution from N-body simulations
6
1
0 Daniel B. Thomas1,2(cid:63), Marco Bruni1, David Wands1
2
g 1 - Institute of Cosmology and Gravitation, University of Portsmouth, Dennis Sciama Building, Burnaby Road,
u Portsmouth, PO1 3FX, UK
A 2 - Department of Physics, University of Cyprus, Aglantzia, Nicosia, 2109
6
2 Accepted;Received;inoriginalform
]
O
ABSTRACT
C Newtonian simulations are routinely used to examine the matter dynamics on non-
. linear scales. However, even on these scales, Newtonian gravity is not a complete
h
descriptionofgravitationaleffects.Apost-Friedmannapproachshowsthattheleading
p
order correction to Newtonian theory is a vector potential in the metric. This vector
-
o potentialcanbecalculatedfromN-bodysimulations,requiringamethodforextracting
r the velocity field. Here, we present the full details of our calculation of the post-
t
s Friedmann vector potential, using the Delauney Tesselation Field Estimator (DTFE)
a
code. We include a detailed examination of the robustness of our numerical result,
[
includingtheeffectsofboxsizeandmassresolutionontheextractedfields.Wepresent
2 the power spectrum of the vector potential and find that the power spectrum of the
v vector potential is ∼ 105 times smaller than the power spectrum of the fully non-
9 linear scalar gravitational potential at redshift zero. Comparing our numerical results
9 toperturbativeestimates,wefindthatthefullynon-linearresultcanbemorethanan
7
order of magnitude larger than the perturbative estimate on small scales. We extend
0
theanalysisofthevectorpotentialtomultipleredshifts,showingthatthisratiopersists
0
overarangeofscalesandredshifts.Wealsocommentontheimplicationsofourresults
.
1 for the validity and interpretation of Newtonian simulations.
0
5 Key words: gravitation – cosmology: theory – cosmology: large-scale structure of
1 the universe.
:
v
i
X
r On the largest scales in cosmology, theoretical calcula- 2015;Flender&Schwarz2012;Hauggetal.2012;Koppetal.
a
tionscanbecarriedoutusingstandardcosmologicalpertur- 2014),aswellasexaminingtherelativisticinterpretationof
bation theory. These calculations fully encompass General thesimulations(Chisari&Zaldarriaga2011;Green&Wald
Relativity (GR) but are limited to scales where the pertur- 2012; Bruni et al. 2014; Adamek et al. 2013). These stud-
bations,inparticularthedensityperturbation,aresmall.On ieshavepredominantlyfocussedonwhetherthedynamicsof
smallerscales,wherethefocusisonnon-linearstructurefor- densitycontrastandscalarpotentialaccuratelymatchthose
mation,NewtonianN-bodysimulationsareused.Thesesim- of GR.
ulationsdonotrequirethatthedensitycontrastissmall,but
In this paper, we are mostly interested in another im-
they suffer from the limitations of being Newtonian rather
portant limitation of Newtonian simulations: Even if the
than GR simulations. There is an entire field in cosmology
matter dynamics are being computed correctly, there are
dedicated to developing, running and analysing these New-
cosmological quantities of interest on non-linear scales that
tonianN-bodysimulations.Therehasbeensporadicinterest
havenocounterpartinNewtoniantheory.Examplesofthese
inunderstandingtheuseofNewtoniantheoryincosmology
quantities include the difference between the two scalar po-
(Tomita 1991; Shibata & Asada 1995; Takada & Futamase
tentials, gravitational waves and the vector potential in the
1997; Matarrese & Terranova 1996; Carbone & Matarrese
metric,allofwhichmustexistonnon-linearscalesinaGR
2005; Hwang et al. 2008; Hwang & Noh 2013; Milillo et al.
universe.
Theseextraquantitieswouldnaivelybeexpectedtobe
(cid:63) E-mail:[email protected] small if the Newtonian simulations are a good approxima-
(cid:13)c 2014RAS
2 D. B. Thomas et al. 2014
tiontoaGRuniverse.However,explicitlycalculatingthese the pertinent details of the post-Friedmann formalism and
quantitieshasseveraladvantages:Tostartwith,itwouldbe show the equation governing the vector potential. We will
good to have a quantitative check of whether these quanti- also present our definitions and notation regarding vector
tiesaresmall,andindeedhowsmalltheyare.Inparticular, power spectra. In section 2, we explain how the relevant
asweentertheeraofprecisioncosmology,weneedtocheck fieldswereextractedfromN-bodysimulationsandexamine
that these quantities will not affect the observables at the the robustness of this extraction. In section 3, we show the
percent-level. Furthermore, checking that these quantities power spectrum of the vector gravitional potential and its
are negligible provides a quantitative check on the Newto- time evolution, as well as comparing it to the closest an-
nian approximation in a ΛCDM cosmology. alytical results in the literature. We conclude in section 4.
Wewillbeworkingwiththepost-Friedmannformalism AppendixAcontainssomedetailsaboutvectorpowerspec-
(Mililloetal.2015;Milillo2010).Thisgeneralisestocosmol- tra and in Appendix B we show results from some of the
ogy the weak-field (post-Minkowski) approximation, with numerical tests that were carried out. Additional plots are
a post-Newtonian style expansion (Chandrasekhar 1965; availableasonlinesupplementarymaterial,dividedamongst
Weinberg 1972; Poisson & Will 2014) in inverse powers of three files: Resolution and BoxSize Dependence.pdf (here-
the speed of light c of the perturbative quantities. These afterRB),GridSize and Binning Dependence.pdf(hereafter
expansions need to be performed differently in cosmology GB) and ConsistencyChecks.pdf (hereafter CC).
compared to in the Solar System due to the different sit-
uations and aims in the two cases. For example, the time-
time and space-space components of the metric need to be
treated at the same order in cosmology in order for the re- 1 POST-FRIEDMANN FORMALISM
sultingequationstobeaconsistentsolutionoftheEinstein
The post-Friedmann approach is developed in Milillo et al.
equations.
(2015)andMilillo(2010),seethereforthefulldetails.This
The post-Friedmann formalism, when linearised, cor-
approach considers a dust (pressure-less matter) cosmology
rectly reproduces conventional linear perturbation theory
withacosmologicalconstant.TheperturbedFLRWmetric,
and can thus describe structure formation on the largest in Poisson gauge, is expanded up to order c−5, keeping the
scales.Moreimportantly,theleadingorderequationsinthe
g and g scalar potentials at the same order:
00 ij
1/cexpansioncanbeexaminedandareexpectedtoyieldthe
(cid:20) (cid:21)
non-linear Newtonian equations. Note that in this “Newto- g = − 1− 2UN + 1 (cid:0)2U2 −4U (cid:1)
nian”regime,thedensitycontrasthasnotbeenassumedto 00 c2 c4 N P
besmall.Theequationsinthisregimewillbeshowninsec- aBN aBP
tion 1, essentially comprising the Newtonian equations, as g0i = − c3i − c5i (1)
expected,plusanadditionalequation.Thisadditionalequa- (cid:20)(cid:18) (cid:19) (cid:21)
tionshowshowthevectorpotentialinthemetric,thelowest gij = a2 1+ 2cV2N + c14 (cid:0)2VN2 +4VP(cid:1) δij+ hci4j
orderbeyond-Newtonianquantity,isgeneratedbythemat-
terdynamics.Thisvectorpotentialisthebeyond-Newtonian The g and g scalar potentials have been split into the
00 ij
quantitythatwewillexamineinthispaper,itisthecosmo- Newtonian (U , V ) and post-Friedmann (U , V ) com-
N N P P
logical manifestation of the ubiquitous relativistic effect of ponents. Similarly, the vector potential has been split up
frame dragging. This effect has been measured in the Solar into BN and BP. Since this metric is in the Poisson gauge,
i i
System by Gravity Probe B (Everitt et al 2011). the three-vectors BN and BP are divergenceless, BN = 0
i i i,i
In this paper, we present a calculation of this vector and BP = 0. In addition, h is transverse and tracefree,
i,i ij
potentialbasedonextractingthedensityandvelocityfields hi = h,i = 0. Note that at this order, h is not dynam-
i ij ij
fromN-bodysimulations.WeexpandontheresultsofBruni ical, so it does not represent gravitational waves. From a
et al. (2014), which was the first calculation of an intrin- post-Friedmann viewpoint, there are two different levels of
sically relativistic quantity on fully non-linear scales from perturbationsinthetheory,correspondingtotermsoforder
large scale cosmological matter fields, rather than from in- c−2 and c−3, or of order c−4 and c−5 respectively. Defining
dividual astrophysical occurences. The main focus in this “resummed”variables,suchasΦ=2U +c−2(cid:0)2U2 −4U (cid:1),
N N P
paper is to present the method used to extract this vector thencalculatingtheEinsteinequationsandlinearisingthem,
potential from N-body simulations. In particular, we exam- reproduceslinearGRperturbationtheoryinPoissongauge.
ine the robustness of the numerical extraction of the vector Thus, this approach is capable of describing structure for-
potentialandpresentthetestswecarriedouttoexaminethe mation on the largest scales.
numericaleffectsofsimulationparametersontheextraction, Forsmallerscales,inadustcosmology,weareinterested
which were not presented in Bruni et al. (2014). intheweakfield,slowmotion,sub-horizon,quasi-staticand
The main physical results of this paper are figures 3 negligiblepressureregime.Thisissimplyderivedbyretain-
and 8, showing the power spectrum of the vector potential ing only the leading order terms in the c−1 expansion and
atredshiftzeroanditsevolutionwithtimerespectively.Ad- upondoingsowerecoverNewtoniancosmology,albeitwith
ditionally,wehavepresentedtheratioofthevectorpotential a couple of subtleties. The first is that the space-time met-
power spectrum to that of the scalar potential in figures 5 ric is a well-defined approximate solution of the Einstein
and 9. The results on the magnitude and evolution of the equations. The second is that we have an additional equa-
power spectrum of the vector potential in this paper were tion, which is a constraint equation for the vector gravita-
usedinThomasetal.(2014)toexaminethepossibleweak- tional potential BN. The full system of equations obtained
i
lensing consequences of the vector potential. fromtheEinsteinandhydrodynamicequations(Mililloetal.
Thispaperislaidoutasfollows.Insection1,wepresent 2015), given the evolution of the background a(t), is as fol-
(cid:13)c 2014RAS,MNRAS000,1–??
The non-linear post-Friedmann vector potential 3
lows.
Table 1.Parametersforthesimulations.
dδ vi
+ ,i (1+δ)=0 (2) Boxsize Particle Massresolution Numberof Softening
dt a
dv a˙ 1 (h−1Mpc) number (108 M(cid:12)) Realisations (h−1kpc)
i + v = U (3)
dt a i a N,i 80 5123 3.97 8 6.25
1 4πG 80 5123 3.97 1 4.0
c2a2∇2VN =− c2 ρbδ (4) 140 7683 6.31 8 6.25
2 140 5603 16.3 8 6.25
c2a2∇2(VN −UN)=0 (5) 160 10243 3.97 3 6.25
(cid:20) (cid:21) 160 8803 6.26 3 6.25
1 2a˙ 2 1 8πGρ
U + V˙ − ∇2BN = b (1+δ)v(6) 160 6403 16.3 8 6.25
c3 a2 N,i a N,i 2a2 i c3 i 160 6403 16.3 1 5.0
160 3203 130 8 15.0
As expected, we have the Newtonian continuity and Euler
200 10243 7.76 2 6.25
equationsfromthehydrodynamicequationsaswellasPois- 240 9603 16.3 3 6.25
sons equation from the Einstein equations. Note that the 240 4803 130 8 15.0
time derivative here is the convective derivative, dA/dt = 320 6403 130 8 15.0
∂A/∂t+viA /a,foranyquantityA.TheEinsteinequations
,i
yieldtwoadditionalequations:Thefirstisanequationforc-
ingthescalarpotentialsV andU tobeequal,consistent
N N Table 2.Redshiftsusedtoprobetimeevolutionofquantities.
withtherebeingonlyonescalarpotentialinNewtonianthe-
ory. Note that some approaches consider the potentials to
ScaleFactor Redshift Colourontimeevolutionplots
be a priori equal at leading order whereas here we assumed
the full generality of GR and the equality of the potentials 0.33 2.0 black
arosenaturallyontakingtheNewtonianregime.Thesecond 0.4 1.5 red
0.5 1.0 magenta
additionalequationrelatestheleadingordervectorgravitia-
tionalpotential,BN,tothemomentumofthematter.Thus, 0.6 0.67 yellow
i 0.7 0.43 green
evenintheregimewherethematterdynamicsarecorrectly
0.8 0.25 cyan
described by Newtonian theory, the frame-dragging poten-
0.9 0.11 blue
tial BiN should not be set to zero; this would correspond 1.0 0.0 brown
to putting an extra constraint on the Newtonian dynamics.
Wenotethatthereisasimilarequationinseveralotherfor-
malismsintheliterature(Takada&Futamase1997;Hwang 2 SIMULATIONS
&Noh2013;Green&Wald2012).Wecanseefromequation
Our simulations have all been run using the publicly avail-
(6) that the potential BN is sourced by the vector part of
i able N-body code Gadget2 (Springel 2005). Many simula-
theenergycurrentρ(cid:126)v.Thisismadeapparentbytakingthe
tions have been run in order to quantify the effects of box
curl of this equation, which gives
size and mass resolution on the quantities that we are ex-
∇×∇2B(cid:126)N =−(cid:0)16πGρ a2(cid:1)∇×[(1+δ)(cid:126)v], (7) tracting, see table 1 for a full list of the simulations. All of
b
thesimulationswererunwithdarkmatterparticlesonly,as
where the source term on the right hand side splits up into the equation for the vector potential is derived for a pres-
threeterms:thevorticity∇×(cid:126)v andthentwofurtherterms, surelessmatterandcosmologicalconstantcosmology.Toal-
low comparison to previous studies of vorticity (Pueblas &
∇×[(1+δ)(cid:126)v]=∇×(cid:126)v+δ∇×(cid:126)v+∇δ×(cid:126)v. (8)
Scoccimarro 2009), the simulations were run with the cos-
It is equation (7) that will be used for the rest of the mological parameters Ωm = 0.27, ΩΛ = 0.73, Ωb = 0.046,
paper. Since the matter dynamics are not affected at this h=0.72,τ =0.088,σ8 =0.9andns =1.Allofthesimula-
order, i.e. they are described by the standard Newtonian tions started at redshift 50 and had their initial conditions
equations (2-4), the density and velocity fields sourcing the created using 2LPTic (Crocce et al. 2006). Our final result
vector potential are Newtonian and can be extracted from for the vector potential is taken from the three 160h−1Mpc
N-body simulations. Using the definitions of vector power simulationswith10243 particles,thesewillbereferredtoas
spectra in Appendix A, the power spectrum of the vector the high-resolution (HR) simulations.
potential is given by
(cid:18)16πGρ a2(cid:19)2 1 2.1 Tesselation
P (k)= b P (k), (9)
B(cid:126)N k2 k2 δv Toextractthenecessaryfieldsfromthesimulations,theDe-
launey Tesselation Field Estimator (DTFE) code was used
with
(Cautun&vandeWeygaert2011).Standardmethodsofex-
P =P (k)+P (k)+P (k) (10) tracting fields from N-body simulations, such as Cloud-In-
δv ∇×(cid:126)v δ∇×(cid:126)v (∇δ)×(cid:126)v
Cells (CIC) (Hockney & Eastwood 1981) work well for the
+P (k)+P (k)+P (k).
(∇δ×(cid:126)v)(∇×(cid:126)v) (∇δ×(cid:126)v)(δ∇×(cid:126)v) (δ∇×(cid:126)v)(∇×(cid:126)v)
densityfield,astheparticles,bydefinition,sampletheden-
Unless stated otherwise, all plots of the gravitational po- sity field well. However, these methods have several short-
tentials show the dimensionless power spectrum ∆(k), see comings when applied to the extraction of velocity fields:
Appendix A for conventions. One is that the field is only sampled where there are par-
(cid:13)c 2014RAS,MNRAS000,1–??
4 D. B. Thomas et al. 2014
ticles, so in a low density region the velocity field is arti- sis here as smaller box sizes have systematically less power.
ficially set to zero. In addition, the extracted field will be SeeAppendixB6fortheresultsfromthesesimulationsand
amass-weighted,ratherthanvolume-weightedfield.Acon- howtheycomparetothelargerboxsizes.Forfurtherresults
sequence of these shortcomings is that, as the grid size is regardingtheeffectsofasmallsimulationboxoncosmologi-
increased,thevelocityfieldwillnotconverge.Infact,itwill calquantities,see(Bagla&Prasad2006;Bagla&Ray2005;
become zero in an increasing proportion of the grid cells as Gelb & Bertschinger 1994).
the grid size increases. Several authors have looked at us-
ing the Delauney tesselation (Schaap & van de Weygaert
2.2.1 A note on error bars
2000; van de Weygaert & Schaap 2009; Bernardeau & van
de Weygaert 1996) for astrophysical applications including SincewehaveonlythreerealisationsofourHRsimulations,
theexaminationofvelocityfields.SeealsoPueblas&Scoc- we cannot compute meaningful error bars. Thus, we have
cimarro (2009) for comparisons of extracting velocity fields not included any error bars on the majority of our plots.
with tesselations rather than more standard methods. The Instead,infigures1,2,5and6,wehaveplottedtheresults
DTFEcodeconstructstheDelauneytesselationofthesetof from the three individual realisations, in order to illustrate
particles, consisting of tetrahedra whose nodes are located by how much the results vary. Unless stated otherwise, the
at the particles positions. The tetrahedra are constructed results shown in the other plots show the average over the
such that the circumsphere of each tetrahedron does not realisations.Weexplicitlyexaminethevariationamongstre-
contain any of the particles except for the particles located alisationsinAppendixB4forseveralquantities,notablythe
atthenodesofthetetrahedroninquestion.Thismakesthe vorticity and vector potential. In particular, we note there
tesselationunique.Theparticles’velocitiesarethenlinearly thatwhenconsideringthevectorpotential,cosmicvariance
interpolatedacrosseachtetrahedron,yieldingavalueforthe on the largest scales affects smaller scales, as explained by
smoothed velocity field and its gradients at every point in a perturbative analysis (Lu et al. 2008; Hui-Ching Lu et al.
thesimulationvolume.AregularN3 gridislaiddownand 2009).SeeB4formorediscussionofthis.Wealsonotethere
grid
thecodesamplesNsamples pointsatrandomineachgridcell thatthevariationofthevorticityamongstrealisationsseems
andaveragesthefieldoverthesepoints,givingavalueforthe tobelargerthanthevariationofthedensity,althoughthere
smoothedfieldineachgridcell.Oncethefieldsareobtained seems to be no discussion of this in the literature.
on the regular grid, the power spectra are calculated using
the standard process of averaging the modulus-squared of
2.2.2 Mass Resolution
theFouriercoefficientsoveragivenrangeofk.Fortheanal-
yseshere,weusedNgrid/4bins,althoughvaryingthisvalue We have examined the dependence of the density, veloc-
does not affect the results (see Appendix B7). ity divergence, vorticity and vector potential on the mass
resolution of the simulations. For the density and velocity
divergence there is evidence for a mild dependence on mass
2.2 Convergence and Tests resolution for both of these fields on smaller scales. This
is likely to be due to the DTFE window function, which
It is important to ensure that our numerical result for the
cannot be compensated for, rather than a mass-resolution
vectorpotentialisrobustandindependentofthesimulation
dependence of the field itself. There is no evidence of any
parameters.Inthissubsectionwewillpresenttheresultsof
mass-resolution dependence of these fields on larger scales.
our examination into the effects of different simulation pa-
The variation of the density and velocity divergence with
rametersontheextractedvectorpowerspectrum.Sincethe
mass resolution can be seen in figures 1 and 2 of file RB.
velocityanddensityfieldsbothcontributetothesourcefor
The effect of the small-scale mass-resolution dependence is
the vector potential, we will examine the density, vorticity
negligible for our HR simulations, as seen when compar-
and velocity divergence spectra too: We will examine their
ing to alternative methods of calculating the density power
behaviour individually, compare them to other studies and
spectrum.
methods of extraction and also consider the consistency of
The dependence of the vorticity power spectrum with
the extracted fields through the relations
mass resolution is shown in figure B1. The power spectrum
k2P (k) = P (k) shows spurious additional power when the mass resolution
δ ∇δ
k2P (k) = P (k)+P (k). (11) is insufficient. However, once the resolution is sufficient, of
(cid:126)v ∇·(cid:126)v ∇×(cid:126)v order 109M , there is no evidence for any systematic de-
(cid:12)
The box size and mass resolution of the simulation are pendenceonmassresolution.Thisdependenceonmassres-
thetwomainparameterswhoseeffectontheextractedfields olution,followedbyconvergencearound∼109M ,matches
(cid:12)
needs to be examined. In addition, we have examined the previous findings, notably those of Pueblas & Scoccimarro
effect of varying the grid size and N , which are both (2009).
samples
internal DTFE parameters. The parameters of the differ- InfigureB2,weshowthedependenceofthevectorpo-
ent simulations used are in table 1. We chose the soften- tential on mass resolution. There is a clear dependence of
inglengthsoftheN-bodysimulationstobeconsistentwith thevectorpotentialonmassresolution,similartothatseen
Pueblas & Scoccimarro (2009) in order to recreate their for the vorticity. However, there are several differences. In
studyofthevelocitydivergenceandvorticity,howevervary- particular, the mass-resolution dependence seems to be less
ing the softening length did not influence the results, see importantforsmallerscales,wherethereisagreaterdepen-
Appendix B5. dence on box size (see later). In addition, the dependence
Although we did run some simulations with a box size on mass-resolution is still apparent around 109M . How-
(cid:12)
below140h−1Mpc,wehavenotincludedtheseintheanaly- ever, once there mass resolution has improved to around
(cid:13)c 2014RAS,MNRAS000,1–??
The non-linear post-Friedmann vector potential 5
6×108M ,thereisnoevidenceofamassresolutiondepen- thefactorofk inequation12takentobethevaluedefining
(cid:12)
dence of the vector potential. the centre of the bin. For the red curve, the exact k-value
To show this further, figure B8 shows the higher res- foreachmodeisusedwhencomputingthesumineachbin.
olution simulations in more detail, complete with the indi- For small bins, or fields where the values vary slowly as a
vidual realisations of the HR simulations. The y-axis here function of k, these two should agree and indeed they do
is k2P (k) in order to show the variance more clearly over for smaller scales where our (logarithmic) bins are smaller.
B(cid:126)
the range of scales being considered. The cyan line shows There is a difference between the methods for the largest
the simulation with the worst resolution (16.3×108M ) of scales in our simulations, this will be discussed below for
(cid:12)
thoseinthisplotandindeedthissimulationshowsasystem- each test.
atic deviation on the largest scales. The better resolution For the density field, the two methods for calculating
simulations show better convergence, with the 140h−1Mpc theratiodogivedifferentanswers.However,forbothmeth-
simulationwith7683 particlesbeingconsistentwiththeHR ods,thedeviationiswithin2%foreverybinexceptthefirst.
simulationsforessentiallytheentirerangeunderconsidera- Thus,thisconsistencycheckforthedensityfieldiswellsat-
tion.ThisconvergenceisexaminedfurtherinAppendixB4. isfied for all scales k(cid:62)0.2hMpc−1.
The consistency check for the velocity field is less well
satisfied: there is a sharp divergence in the power spectra
2.2.3 Box Size on the smallest scales, such that the check is not satisfied
within 10% at k ≈ 8h−1Mpc. This shows the effect of the
We have also considered the effect of varying the box size
DTFEwindowfunctionontheextractedfields.Wewillnot
on the extracted power spectra. As expected, there is no
consider the extracted vector potential for k larger than
evidence for any systematic dependence of the density, vor-
k≈8h−1Mpcwhenpresentingourresults.Furthermore,the
ticityandvelocitydivergencepowerspectraontheboxsize
twomethodsshowverydifferentbehaviour:themethodus-
of the simulations. This can be seen in figures 3, 4 and 5 of
ing the average k-value for each bin causes the consistency
fileRB.Notethat,forsufficientlysmallboxes,asystematic
test to fail on large scales. However, with the more exact
deviation can arise, see Appendix B6.
method,theconsistencycheckisverywellsatisfiedonallof
Figure B3 shows the box size dependence of the vector
the largest scales. This suggests that the dominant contri-
potential. As mentioned above, the vector potential does
butiontothebinsonthelargestscalescomesfromthelow-k
show some dependence on box size. The vector potential
end of each bin, hence the overestimation of k2P (k) when
showssignsofadependenceontheboxsizeonscalesbelow (cid:126)v
1h−1Mpc, however this is difficult to entangle from the ef- the average k-value for each bin is used. The strong effect
hereispartlycausedbytherelativelysteepslopeoftheve-
fects of mass resolution and the window function. For box
sizes below 200h−1Mpc, there is no systematic dependence locity power spectrum. We note that this effect would also
come into play when calculating the dimensionless veloc-
of the vector potential power spectrum with box size.
ity power spectrum for binned data. Nonetheless, the good
In Appendix B4, we examine the variation between realisa-
agreement of the consistency check when using the second
tionsforthevectorpotential,andrelateittothebehaviour
methodisstrongevidencethatthederivativesofthevelocty
that might be expected from perturbative arguments. In
field are being calculated correctly.
particular, figure B8 shows how the variation between re-
A further check that we can perform is to extract the
alisations is larger than the effects of box size and mass
resolution for simulations with box sizes below 200h−1Mpc complete momentum field, p(cid:126) = (1+δ)(cid:126)v, and decompose it
andmassresolutionofatleast6×108M .Thus,weexpect into its vector and scalar parts directly rather than dealing
(cid:12)
withderivatives.Thepowerspectrumofthevectorpotential
numerical effects from the simulation parameters to be a
can then be calculated from the vector part of the momen-
sub-dominant source of error as long as the parameters are
tum field, p(cid:126)v, using
within this range.
(cid:18)16πGρ a2(cid:19)2
PB(cid:126)N(k)= k2b Pp(cid:126)v(k). (14)
2.2.4 Consistency Checks
In figure 1 we show the ratio of the vector power spectrum
Thereareafewconsistencychecksthatcanbeperformedon
calculated using the two methods, with the different lines
the different fields that we are interested in. The quantities
correspondingtodifferentindividualrealisations.Thevector
that are used for the vector potential include the density
potentialscalculatedfromthetwomethodsarebroadlycon-
field and its gradients as well as the velocity field and its
sistent, within 20% for most of the range under considera-
gradients. There are two relations between these fields and
tion,andagreeingtowithinafactorof2fork(cid:62)0.2hMpc−1.
their derivatives,
Weareunsurewhatthecausesofthedifferencebetweenthe
k2P (k) = P (k) (12) twomethodsare.Inparticular,wecheckedforwhetherthere
δ ∇δ
k2P (k) = P (k)+P (k). (13) isaneffectcomingfromtheuseofkaveragedoverthebin,as
(cid:126)v ∇·(cid:126)v ∇×(cid:126)v
inthevelocityfieldconsistencycheck,howeverthiseffectis
Wehaveextractedthequantitiesontheleftandrightsides negligibleforthegravitomagneticpotential1.Thedifference
of these relations from our HR simulations and compared
them, see figure 1 in file CC for the ratio P (k)/k2P (k)
∇δ δ
and figure 2 in file CC for the ratio k2P(cid:126)v(k)/(P∇·(cid:126)v(k)+ 1 As an aside, we note that we also calculated the momentum
P∇×(cid:126)v(k)).Inbothcases,twocurvesareplotted,correspond- field by extracting the velocity field and density field separately
ing to two different methods of calculating the ratio. The ateachgridpoint,beforemultiplyingthemtogether.Thepower
blue line shows the ratio exactly as suggested above, with spectrum calculated from this field agrees well with that calcu-
(cid:13)c 2014RAS,MNRAS000,1–??
6 D. B. Thomas et al. 2014
between the methods is larger than the variation amongst
realisations for either method.
We can also extract the momentum field directly us-
ing a standard cloud-in-cells (CiC) approach (Hockney &
Eastwood 1981), and compare this to the momentum field
extracted using the DTFE code. The ratio between these
fieldsisshowninfigure2.Thereisgoodagreementbetween
the two methods of computing the momentum power spec-
trumonlargerscales,butwithadivergencebetweenthetwo
methodsonsmallerscales.Itisunclearwhichmethodwould
be expected to be more accurate on these smaller scales:
the DTFE method suffers from having a window function
that cannot be deconvolved, however the CiC method will
have cells with a zero momentum field, due to the lack of
nearbyparticles,forasufficientlylargegrid.Infact,theCiC
method does not converge as the grid size is increased. We
useda5123 gridfortheCiCcode,althoughwecheckedthat
changingthisto256or1024doesnotsignificantlyaffectthe
results.UnliketheDTFEmethod,derivativescannotbedi-
rectly extracted with the CiC method, so the consistency
Figure 1. The ratio of the vector potential power spectra com-
checks performed earlier for the DTFE method cannot be
puted using the vector part of the momentum field and the curl
applied to the CiC method. This also means that the first
of the momentum field. The blue, magenta and red curves show
method of extracting the vector potential, using the curl of
theratioforthethreerealisationsoftheHRsimulations,andthe
the momentum field, cannot be carried out with the CiC black (dashed) curve shows the average over these three. There
method. is reasonable agreement between the two power spectra for the
We present the vector power spectrum from both the smaller scales, however the two methods diverge for the largest
momentumfieldandthecurlmethodintheresultssection. scalesandthereisadifferenceofafactorof5atthelargestscales.
Wenotethatthelevelofagreementbetweenfigures1and2 Formostoftherangeofkunderconsideration(k(cid:62)0.2hMpc−1),
suggeststhatourvectorpotentialpowerspectrumisrobust thetwovectorpowerspectraagreetowithinafactorof2.
andcorrecttowithinafactorof2.Itisuncleartouswhich
methodshouldbetrustedmore;whilstthemomentumfield
method is simpler, the derivative method allows us to ex-
amine the different components, notably the vorticity, and
check that it behaves as expected. The differences between
thetwomethodsdonotaffecttheobservabilityofthevector
potential, see Bruni et al. (2014) and Thomas et al. (2014).
2.2.5 Comparison to previous findings
There are several works in the literature to which we can
compare our findings on the velocity field and its compo-
nents. As mentioned above, the vorticity and velocity di-
vergencepowerspectrawereextractedfromN-bodysimula-
tions in Pueblas & Scoccimarro (2009) using an alternative
implementation of the Delauney tesselation. They found a
strong dependence on resolution of the extracted vorticity
powerspectrumandanapproximatescalingofthevorticity
powerspectrumwiththeseventhpowerofthelineargrowth
factor.
The vorticity and velocity divergence power spectra in
Figure 2. The ratio of the vector potential power spectra com-
Pueblas&Scoccimarro(2009)areconsistentwiththespec-
putedusingthevectorpartofthemomentumfieldcalculatedus-
tra extracted for this paper and we found the same res- ingtheCloud-in-CellsmethodandtheDTFEmethod.Theblue,
olution dependence of the vorticity power spectrum (see magentaandredcurvesshowtheratioforthethreerealisations
above). However, as detailed in Appendix B2, we do not oftheHRsimulations,andtheblackcurveshowstheaverageover
findthesamescalingofthevorticityspectrumwiththesev- thesethree.Thetwomethodsagreeverywellonlargerscales,but
enth power of the linear growth factor (D ): Although this divergeforthesmallestscales.
+
lated by extracting the momentum field as a single field. The
sameagreementisnotobtainedwhenextractingthefieldδ2 and
comparing to squaring the density field, when using either the
DTFEcodeoraCiCmethod.
(cid:13)c 2014RAS,MNRAS000,1–??
The non-linear post-Friedmann vector potential 7
scaling seems to hold at low redshift, it no longer holds at multiple realisations with the same resolution, in order to
redshift one and beyond. At these earlier times, the power compare our findings. In addition, there is no examination
spectrum is smaller than expected from the growth factor of the time dependence and thus no confirmation or rejec-
to the seven scaling, so the vorticity power spectrum must tion of the D7 scaling of the vorticity spectrum at higher
+
have grown by less at redshift two than expected. redshifts.
Tworecentpublications(Zhengetal.2013;Kodaetal.
2013) have examined the velocity field from the point of
view of redshift space distortions. In these works, a differ-
3 RESULTS
ent method of extracting velocity fields is used, the nearest
particlemethod.Inthismethod,thevelocityatagridpoint In this section we present the power spectrum of the post-
is given by the velocity of the nearest particle to that grid FriedmannvectorpotentialascalculatedfromN-bodysim-
point. See those works for comments on the differences be- ulations. We show the power spectrum at z = 0 and the
tweenthenearestparticleandDelauneytesselationmethods different components of the source, as well as the evolution
ofextractingthevelocitypowerspectra.Here,wenotethat of the power spectrum between z = 2 and z = 0. In addi-
there appear to be pros and cons to both methods, with tion,weshowtheratiobetweenthevectorandscalarpower
no clear “better” method. It would be interesting to exam- spectra, and examine the time evolution of this quantity
ine how close the agreement between the vector potentials aswell.Thepowerspectraplottedforthescalarandvector
extracted by the DTFE and nearest particle methods is. gravitationalpotentialsarethedimensionlesspowerspectra.
Nonetheless, there are some general observations that The closest analytic result to our calculation is the second
canbecomparedbetweentheseworks.Notably,themagni- orderperturbativevectorpotentialcalculatedinHui-Ching
tudeofthevelocityandvorticityspectraisfoundtobesim- Lu et al. (2009). We will compare our results to theirs at
ilar, considering the differences in cosmological parameters. redshift z=0, as well as comparing the time evolution.
Also, the onset of non-linearity is found to occur at lower
k for the velocity divergence than for the density. In addi-
3.1 Results at redshift zero
tion, Zheng et al. (2013) finds a strong dependence of the
curl component of the velocity field on the resolution, simi- Infigures3and4,weshowthepowerspectrumofthepost-
larlytoboththispaperandPueblas&Scoccimarro(2009). Friedmann vector potential as well as the standard Newto-
Theyalsofindatimedependenceofthiscomponentthatis nian scalar potential, at z = 0, for the curl and momen-
approximately D7 up to z = 2, although this relationship tum field methods of extraction respectively. As expected,
+
breaks down by up to a factor of two for certain redshifts both methods show that the scalar potential is small over
andscales.Asmentionedabove,whilstoursimulationsalso all scales and the vector potential is subdominant. There is
findthistimedependenceofthevorticityatlowredshift,we a quantitative difference between the two methods on the
findthattherelationshipbreaksdownforz>1.Thereisno largestscales,butthisdifferenceisnotsufficienttoalterthe
examination of multiple realisations in Zheng et al. (2013) expected qualitative behaviour. Notably, the effect of the
and, similarly to the comments in Appendix B2 regarding vectorpotentialonweak-lensingpowerspectra,asexamined
Pueblas & Scoccimarro (2009), the difference between our inThomasetal.(2014),willremainnegligible,regardlessof
realisationsissufficienttoexplainthedifferencebetweenour which method is used to calculate the vector potential. We
results and those of Zheng et al. (2013). have been unable to determine the reason for this discrep-
The broad agreement between different methods, in- ancy and it is unclear to us which method should, a priori,
cludingagreementregardingresolutiondependenceandcon- be expected to be more accurate.
vergence, is promising. Details of the vorticity field and its Infigures5and6,weshowtheratiobetweenthepower
evolution require further study, but the vorticity is a sub- spectra of the vector and scalar gravitational potentials at
dominantcontributiontothevectorpotential.Asthesimu- redshift zero, for the two methods of extracting the vec-
lations and snapshots used in the papers mentioned in this tor potential. We plot the ratios for all three individual re-
section are different to ours, it is not possible to compare alisations of the HR simulations. For the curl method, as
the methods and extracted fields any more precisely. We shown in Bruni et al. (2014), this ratio is approximately
note that the three works mentioned here do not have mul- 2.5×10−5. This ratio does not vary significantly over the
tiple realisations of their high resolution simulations, so we rangeofscalesconsidered,althoughthereisaslightincrease
areunabletodetermineifthevariationinvorticitybetween towards smaller scales. However, for the momentum field
realisations found by us is reproduced (see Appendix B4). method, the ratio is not approximately constant due to the
As this manuscript was being prepared, Hahn et al. decreased power on large scales. We will compare this be-
(2014) appeared on the arxiv. This paper investigates the haviourtotheanalyticsecond-orderperturbativebehaviour
properties of velocity divergence and vorticity and confirms shortly, here we just note that the curl method produces
many of the findings of Pueblas & Scoccimarro (2009). In qualitative behaviour that is closer to the analytic predic-
particular,theyagreewithourresultsregardingtheconver- tion.
gence of the DTFE code for sufficient mass resolution and In figure 7, we show the power spectra of the three
our finding of a resolution dependence of the velocity di- sources of the vector potential using the curl method, see
vergence, which did not appear in Pueblas & Scoccimarro equation (8). The power spectra plotted here are given by
(2009). They use a different method to compute the vortic- P(k)/(cid:0)f2H2(2π)3(cid:1), where H is the conformal time Hubble
ityandvelocitydivergencepowerspectra,whichagreeswith constant and f =dlnD/dlna is the logarithmic derivative
the DTFE code for sufficient resolution. However, as with of the linear growth factor D. These units are chosen such
the previous papers, there seems to be no examination of that the power spectrum of the velocity divergence agrees
(cid:13)c 2014RAS,MNRAS000,1–??
8 D. B. Thomas et al. 2014
Figure 3. The scalar (dashed red line) and vector (solid blue Figure5.Theratioatredshiftzerobetweenthevectorpotential,
line) gravitational potential power spectra at redshift zero, with calculated using the curl method, and the scalar potential. The
the vector potential calculated using the curl method. The lin- three curves show the ratio for the three realisations of the HR
eartheoryscalarpotentialisshownforcomparison(dottedblack simulations.
line).
Figure6.Theratioatredshiftzerobetweenthevectorpotential,
Figure 4. The scalar (dashed red line) and vector (solid blue calculatedusingthemomentumfieldmethod,andthescalarpo-
line) gravitational potential power spectra at redshift zero, with tential.Thethreecurvesshowtheratioforthethreerealisations
thevectorpotentialcalculatedusingthemomentumfieldmethod. oftheHRsimulations.
Thelineartheoryscalarpotentialisshownforcomparison(dotted
blackline).
of the relationship between Newtonian simulations and GR
on fully non-linear scales. The small magnitude of the vec-
with the density power spectrum on linear scales and have tor potential suggests that running Newtonian simulations
the same units as the matter power spectrum, following issufficientlyaccurateforcosmologicalpurposes,whereasa
Pueblas&Scoccimarro(2009).Thevorticity,althoughoften larger calculated value for the vector potential would sug-
ignoredinperturbationtheory,istheonlyoneofthesethree gestthattheapproximationstakeninderivingthefullynon-
quantities that is linear in perturbations. This figure shows linearNewtonianequationsdonotholdsufficientlywell.As
thatitisnegligiblecomparedtotheothertwocomponents, far as relating Newtonian and relativistic cosmologies goes,
sothevectorpotentialisbeingpredominantlygeneratedby in the language of Green & Wald (2012), the smallness of
non-linear effects. this vector potential allows the use of the abridged dictio-
SincethisvectorpotentialisthefirstcorrectiontoNew- nary in Chisari & Zaldarriaga (2011), rather than the dic-
toniantheory,thiscalculationisthefirstquantitativecheck tionaryproposedinGreen&Wald(2012).Wenotethatthe
(cid:13)c 2014RAS,MNRAS000,1–??
The non-linear post-Friedmann vector potential 9
Figure 7. The power spectra of the three source terms for the Figure 8.Theevolutionofthevectorpotentialforsixdifferent
vectorpotentialinequation(8),thevorticity(dashedblackline), wavenumbers. From top to bottom, these are k = 0.23hMpc−1
∇δ×(cid:126)v (dot-dashed blue line) and δ∇×(cid:126)v (solid red line). The (brown), k = 0.55hMpc−1 (black), k = 0.79hMpc−1 (cyan),
powerspectraplottedherearegivenbyP(k)/(cid:0)f2H2(2π)3(cid:1),such k = 1.01hMpc−1 (blue), k = 2.51hMpc−1(magenta) and k =
that the power spectrum of the velocity divergence agrees with 5.03hMpc−1 (red).
thedensitypowerspectrumonlinearscalesandensuringthatall
of the power spectra have the same units, following Pueblas &
Scoccimarro(2009).Thelinearmatterpowerspectrumisshown
asadottedmagentalineforcomparison.
analysis here is for a ΛCDM cosmology, further work is re-
quired to determine the validity of Newtonian simulations
in general dark energy cosmologies.
3.2 Time evolution
In this section we will examine the time evolution of the
vectorpotential,anditsratiotothescalarpotential,forthe
redshiftslistedintable2.Thevectorpotentialisthissection
has been computed using the curl method. In figure 8, we
plot the ratio of the vector potential to the scalar potential
as a function of redshift. The different curves in this plot
show the evolution for different wavenumbers. We can see
that individual k-modes do not exhibit significant growth
over time, although the more non-linear scales do exhibit
slightlymorevariationintime.Similarlytothescalargrav-
itational potential, the vector potential at a fixed scale is Figure 9. The ratio of the vector potential to the scalar po-
tential plotted for six different wavenumbers. From bottom to
not monotonic over time on non-linear scales.
top (at redshift=1), these are k = 0.23hMpc−1 (brown), k =
In figure 9, we plot the ratio of the vector potential to
0.55hMpc−1 (black), k =0.79hMpc−1 (cyan), k =1.01hMpc−1
the scalar potential as a function of redshift. The different
(blue),k=2.51hMpc−1(magenta)andk=5.03hMpc−1 (red).
curves in this plot show the same wavenumbers as in figure
8. The ratio stays fairly constant over time, varying by less
thanafactoroftwoforagivenscale.Acrosstheentirerange
As a perturbative analysis, it is unclear how large a value
of times and scales under consideration, the ratio varies by
of k this calculation should be extended to. Here we will
less than a factor of 4. The ratio between the gravitational
assumeitisvalidonallofthescalesofoverlapbetweenthis
potentials is also not monotonic over the redshift range un-
method and ours.
der consideration for a given scale.
For the curl method of computing the vector power
spectrum,thereissimilarqualitativebehaviourbetweenthe
twomethods,withtheratioofthepowerspectraofthevec-
3.2.1 Comparison to perturbative calculation
tor and scalar potentials being fairly constant and of or-
InHui-ChingLuetal.(2009),ananalyticcalculationofthe der 10−5 in both methods. The difference between the two
vector potential was performed using perturbation theory. methods being that the ratio in Hui-Ching Lu et al. (2009)
(cid:13)c 2014RAS,MNRAS000,1–??
10 D. B. Thomas et al. 2014
isbetweenthevectorandthelineartheoryscalarpotential, evidencethatthevectorpotentialissensitivetothesoften-
whereas the our ratio is between the vector and the fully ing length, binning, number of samples (an internal DTFE
non-linearscalarpotential.Thismeansthatdespitethissim- parameter) or the grid size used in the analysis. There is a
ilarqualitativebehaviour,thepowerspectrumofthevector reasonable agreement between the different methods (curl
potential in Hui-Ching Lu et al. (2009) underestimates the and momentum field) of extracting the vector potential,
fully non-linear value on these scales by up to two orders although there is an unresolved discrepancy between the
of magnitude, the same factor by which the linear theory two methods on the largest scales. We do however note the
scalar potential power spectrum underestimates the power importance of the variation of the vector potential between
spectrum of the fully non-linear scalar potential. realisations, this issue is discussed more fully in Appendix
The momentum field method of calculating the vector B4.
powerspectrumresultsinless-similarqualitativebehaviour. Figures 3 and 8 comprise the main physical results of
It is unclear how well the gravitomagnetic potential would this paper, showing the magnitude of the vector potential
be expected to match the perturbative prediction on these powerspectrumatredshiftzeroanditsevolutionwithtime
scales as the velocity field differs from the linear theory at respectively. The magnitude of the vector potential power
larger scales than the density. spectrum can also be expressed in terms of its ratio to
The power spectrum of the perturbative vector poten- the power spectrum of the scalar potential, as shown in
tial is given in Hui-Ching Lu et al. (2009) as figures 5 and 9. We have shown that the power spectrum
of the vector potential is around 105 times smaller than
(cid:18)2∆ (cid:19)4(cid:18)3g[g(cid:48)+H}](cid:19)2
P (k)= R k2Π(u2) (15) the power spectrum of the scalar potential, over a range
s 5g∞ ΩmH2 of scales and redshifts. These values were used in Bruni
where P is the dimensionless power spectrum of the vec- et al. (2014) and Thomas et al. (2014) when examining
s
tor potential, ∆ is the primordial power of the curvature the observability of the vector potential, showing that it
R
perturbation,g isthegrowthfactorforthescalarpotential, is neglgible for currently planned weak-lensing surveys.
g is a normalisation parameter chosen so g(0)=1, Π is a The small magnitude of the vector potential found here is
∞
function of the transfer function, Ω is the time dependent the first quantitative check of the validity of Newtonian
m
matterdensityandHistheconformalHubbleconstant.The simulations compared to GR on fully non-linear scales and
second term in parentheses contains all of the time depen- supports the use of Newtonian simulations for computing
denceofthevectorpotentialpowerspectrumandessentially cosmological observables. In terms of interpreting the
acts as the growth factor for the vector potential. We have simulations, the small value of this vector potential seems
comparedthisperturbativepredictionforthegrowthfactor to justify the use of the abridged dictionary in Chisari &
ofthevectorpotentialtothegrowthmeasuredinthesimu- Zaldarriaga (2011), rather than the dictionary proposed
lations(seefigure5infileCC).Thisshowsthattheanalytic in Green & Wald (2012), for relating GR and Newtonain
predictionisnotthemainsourceofthetimeevolutionofthe cosmologies.
vector potential. The work carried out so far considers a ΛCDM cosmology,
sothisconclusionmaynolongerbetrueforadarkenergyor
modified gravity cosmology. The post-Friedmann approach
would need to be expanded to include modified Einstein
4 CONCLUSION AND DISCUSSION
equations and/or a fluid with pressure in order to examine
In this paper we have presented the post-Friedmann alternative cosmologies and determine whether the use
frame-dragging vector potential calculated on non-linear of Newtonian-type N-body simulations is still valid in
scales from N-body simulations. We have presented this thosecosmologies.Thepost-Friedmannexpansionhasbeen
vector potential at redshift zero, as well as examining its applied to f(R) gravity and the vector potential calculated
evolution with redshift. We have also presented the tests from f(R) simulations in Thomas et al. (2015). The vector
we have performed in order to establish the robustness of potential in f(R) was found to be larger than in General
our result, including tests of simulation parameters and Relativity. We hope that this, and further extensions to
different methods of extracting the source of the vector the work in this paper, will allow us to understand how
potential. generic the findings in this paper are, and thus justify
We have shown that our density, velocity divergence and one of the most widely used tools in cosmology, N-body
vorticityspectraareconsistentwiththeliteratureandshow simulations. Whilst this manuscript was being prepared for
similar behaviour regarding convergence tests, particularly submission, Adamek et al. (2014) appeared on the ArXiv.
mass resolution. We do not see the vorticity scaling with Their preliminary results seem to agree with the results of
the seventh power of the linear growth factor D (Pueblas this work. It will be interesting to perform a more in-depth
+
& Scoccimarro 2009) beyond z=1, however the differences comparison once the details of their work are available.
between our results and others’ are within the variance
between realisations. We have noted a larger variation of Acknowledgements We thank Marius Cautun for help
the vorticity than the density and velocity divergence fields withthepubliclyavailableDTFEcodeandHectorGilMarin
between different realisations, a result that does not seem for provision of, and help with, a Cloud-in-Cells code. We
to have been studied in the literature. also thank Marc Manera for useful discussions and tech-
We have shown that there is no evidence for a systematic nical assistance. Some of the numerical computations were
dependenceofthevectorpotentialspectrumonboxsizefor done on the Sciama High Performance Compute (HPC)
boxes smaller than 200h−1Mpc, or on mass resolution with cluster which is supported by the Institute of Cosmology
mass resolution better than 6×108M . There is also no and Gravitation (ICG) and the University of Portsmouth.
(cid:12)
(cid:13)c 2014RAS,MNRAS000,1–??
