2017/10/21-23
第6回観測的宇宙論ワークショップ
＠弘前大学

Approaching and going beyond theshell crossing of structure formation
with perturbation theory
樽家 篤史（京都大学 基礎物理学研究所）
Stéphane Colombi (IAP)

嵯峨承平 (京大基研)

内容
大規模構造の摂動論とその問題
１次元宇宙での考察
３次元への拡張に向けた課題
まとめCosmologic

宇宙大規模構造と宇宙論
銀河サーベイの重要ターゲット：

観測スケール

50

100

40
80

0

150

1

ss? [h[h 1Mpc]
Mpc]

FIG. 4: Measured power spectra for the full LRG and main galaxy samples. Errors are uncorrelated and full window functions are shownin Figure 5. The solid curves correspond to the linear theory ΛCDM fits to WMAP3 alone from Table 5 of [7], normalized to galaxy bias
b = 1.9 (top) and b = 1.1 (bottom) relative to the z = 0 matter power. The dashed curves include the nonlinear correction of [29] for
A = 1.4, with Qnl = 30 for the LRGs and Qnl = 4.6 for the main galaxies; see equation (4). The onset of nonlinear corrections is clearly> 0.09h/Mpc (vertical line).
visible for k ∼

120
80

50
50

50

150
150

0

100

100

50

150
150

100

600,000 gals

50

100

50

150

k [h Mpc-1]

0

0

100

Alam et al. (’16)

100

BOSS DR12 - 0.5 < z < 0.75

BOSS DR9(SDSS-III)

1 1
s s[h
Mpc]
[h
Mpc]
k

BOSS DR1250s [h 1 Mpc]

P(k) [(h-1Mpc)^3]

100

k<0.2–0.3 h/Mpc @ z~0–2 （重力進化の線形〜弱非線形領域）
s (s , s ) [h Mpc ]Our Fourier convention is such that the dimensionless
power ∆2 of [77] is given by ∆2 (k) = 4π(k/2π)3 P (k).

B.

Clustering evolution

80

40

0

2

40

80

2

2

120

s2 (s , s ) [h 2 Mpc2]

L. Anderson et al.

BOSS DR12 - 0.5 < z < 0.75

150

40

7

150

14

S. Alam et al.

s [h 1 Mpc]

10

加速膨張の診断
重力のテストニュートリノ質量和0• バリオン音響振動• 赤方偏移ゆがみ• フリーストリーム減衰

摂動論による大規模構造の記述
精密観測データと比較する理論テンプレートとして、摂動論による解析的計算手法が威力を発揮
大規模構造の非線形重力進化＋ 赤方偏移空間ゆがみ・銀河バイアスの非線形性
最近• くりこみや高次摂動計算などのテクニックの発展• 高速計算可能なパブリックコード(Crocce & Scoccimarro '08, AT & Hiramtsu'08, Matsubara ’08, AT et al. '12, ...)

宇宙論データ解析に応用(WiggleZ, SDSS BOSS, …)

M (+baryon) particle.
ation
ation
&
'
ation
3
f (x, p) = n a {1 + δm (x)} δD &p − m a v(x)'.
(4.25)
3ニュートン重力 & ゆらぎの波長 << ハッブル半径
弱重力→
f (x, p) = n a {1 + δm (x)} δD &p − m a v(x)'.
(4.25)
e zeroth and 3first velocity moments of Eq. (4.23) yields
f (x, p) = n a {1 + δm (x)} δD p − m a v(x) .
(4.25)
e zerothandfirst velocity momentsof Eq.圧力ゼロの渦なし流体(4.23) yields
冷たい暗黒物質(CDM)
+
バリオン
∂δ
1
m
+
∇
[(1
+
δ
)v]
=
0,
(4.26)
e zerothandfirstvelocity
momentsofEq.
(4.23)
yieldsm∂δ1
質量密度場∂tm + a
速度場
∇
[(1
+
δ
)v]
=
0,
(4.26)m∂δ
1Juszkiewicz (’81),Vishniac (’83), Goroff et∂tm + a∇ [(1 + δm )v] = 0,(4.26)al. (’86), Suto & Sasaki (’91),
Makino, Sasaki
∂t 1a
∂v
1 ∂Ψ
& Suto (’92), Jain & Bertschinger (’94), ...
+ 1 (v · ∇) v = − 1 ∂Ψ ,(4.27)∂v∂t + a (v · ∇) v = − a ∂x ,(4.27)∂v
1
1 ∂Ψ
∂t + aa∂x(v · ∇) v = −,(4.27)
∂t
a ∂x
1 2a
∇
Ψ
=
4π
G
ρ
δ
.
(4.28)
mm
a12ニュートンポテンシャル
2
∇
Ψ
=
4π
G
ρ
δ
.
(4.28)mm
2
a1 ∇2 Ψ = 4π G ρ δ .
(4.28)
m m
標準摂動論
2
a
| |
1
パワースペクトル

大規模構造の摂動論=(1)+ (2) + (3) + · · ·(k; t) (k ; t) = (2 )3 D (k + k ) P (|k|; t)perturbative
for longitudinalandtransverse modes(i.e.,Thisψk,k isfinal step issolutions
to explicitly
construct thedisplacement
field itself.
final
explicitly construct
theFourier
displacement
field [46]).
itself. This is no
butstep
canisbetosystematicallydone inspace (e.g.,
Collisionless(Vlasov-Poisson
system)but can be Boltzmannsystematicallyequation
done in Fourier
space (e.g., [46]).

無衝突自己重力多体系
4.3 (Eulerian) Perturbation theory
#

$

∂
p ∂
∂Ψ ∂
+Perturbation−mf (x, p) = 0,
4.3 (Eulerian)theory2∂x ∂pより基本的な記述： ∂t ma ∂xCollisionlessBoltzmannequation (Vlasov-Poisson system)
supplementedwith
thePoissonequation:
CDM+ バリオン＝無衝突粒子からなる自己重力多体系
CollisionlessBoltzmann equation #(Vlasov-Poisson
system)
$
%
#
$
2
2 m
3
∇ Ψ(x)
=
4π
G
a
d
p
f
(x,
p)
−
ρ
.
m
∂
p
∂
∂Ψ
∂
分布関数
3
#
$ f (x, p) = 0,a+−m
無衝突ボルツマン方程式∂ ∂t p ma∂2 ∂x ∂Ψ∂x
∂ ∂p
+
−m
f (x, p) = 0,
2
Here,
m
is
the
mass
of
CDM
(+baryon)particle.
∂t ma ∂x∂x ∂p
（ヴラソフ方程式）
supplemented with the Poisson equation:$
supplemented with the Poisson equation: # %
m
#
$ρ .
2
2 %
3
∇
Ψ(x)
=
4π
G
a
d
p
f
(x,
p)
−
ポアソン方程式
m2
2 m a3 3Single-stream approximation∇ Ψ(x) = 4π G a
d p f (x, p) − ρm .
3aHere, m is the mass of CDM (+baryon) particle.
デルタ関数
&
'
初期設定
Here,
m is the mass of CDM (+baryon)
particle.
3
Ansatz f (x, p) = n a {1 + δm (x)} δD p − m a v(x) .
（シングルストリーム）

速度場
質量密度場
With thisansatz,takingthezerothandfirstvelocity
moments
of
Eq.
(4.23
Single-stream approximation
Single-stream approximation∂δm 1
シングルストリームの仮定では圧力ゼロの流体系に帰着+ ∇ [(1 + δm )v] = 0,&∂ta3
Ansatz f (x, p) = na {1 + δm (x)} &δD p − m a v(x)'3
ただし、この仮定はいずれ破れる（→
Ansatz f (x, p) = n a {1 +摂動計算の破綻）
δ (x)} δ p − m a v(x) .

一次元重力系での例

位置 x

位置 x

速度 v

速度 v

速度 v

速度 v

１次元重力系の時間発展 （宇宙膨張なし）

位置 x

初期条件 シェルクロッシング マルチストリーム
領域の発達
の形成

シングルストリームとして扱えるのはここまで

位置 x

ハローの形成宇宙論的状況だとハロー外部にシングルストリーム領域が広がるスプラッシュバック半径

Diemer et al.
＝シングルストリームとマルチストリームの境界（白いマル）

15 Mpc/h

The Astrophysical Journal, 843:140 (16pp), 2017 July 10

Figure 1. Comparison of conventional “virial” and splashback radii (Rvir and Rsp, shown as orange and white circles). The image shows the projected density througha slice in the L0125 simulation that is 30 h-1 Mpc wide and deep and 15 h-1 Mpc tall. The density ﬁeld is visualized using the GOTETRA code (P. Mansﬁeld et al.
2017, in preparation). Radii are shown for all halos with N200m . 1000 (equivalent to a mass of 1.4 ´ 1011 h-1 M: ), and the mass of the central halo is
M200m = 1.2 ´ 1014 h-1 M: (corresponding to almost a million particles). The splashback radii shown are deﬁned as Rsp87% , which corresponds most closely to thedensity drop measured by the SHELLFISH code (see Section 4.2). For a small fraction of halos, SPARTA could not determine a splashback radius because they had

30 Mpc/hDiemer et al. (’17)

摂動論の UV 問題
ハロー形成が起こる小スケールに目をつぶればシングルストリーム近似にもとづく摂動計算は問題ない？
z ! 0.375

P$k%"Pnowiggle$k%

1.4

COSMIC PROPAGATORS AT TWO-LOOP ORDER

1ループ補正1.3
Blas et al. (’14)

1.1
!!!
!
1.0

!
!

N体計算

k=0.1 [h/Mpc]!
!

２ループ補正

!!!

! !
!
!
!
!!
!
!
!!
!
!
! ! !
!
!
0.05
0.10
0.15

線形理論

小スケールからの寄与大

３ループ補正
0.20

0.25

0.30

k !h" Mpc#

z ! 1.75

JCAP0

0.9
0.00

D 89,
(2014)
Pn-loop (k)PHYSICAL
d lnREVIEWq Kn-loop(k,023502
q) P0 (q)

!
!

!

!

!

!
!!!

z=0.375

1.2

!

!!

!

FIG. 11 (color online). The shapeof the kernelAT
functions
Bernardeau,&
高次計算すると大スケールのゆらぎが小ス
P ðqÞKðk; qÞ (blue solid line), P ðqÞKðk; qÞ (green
FIG. 10 (color online). Regular parts of the density propagator
1.4

Reg Gp−loop ðkÞ at one-, two-, and three-loop order with, respec1þ

1-looptively, solid, dashed, and dotted lines. The calculations
are done
ケールと強くカップリング→
摂動論が破綻
!

$k%

1.3

2-loopNishimichi(’14)
dashed line) for k ¼ 0.1 h=Mpc andP0 ðqÞK 3-loop ðk; qÞ(red dotted line) as a function of q for z ¼ 0.5.
00

FIG. 8. Diagrammatic representation of the leading diagrams in
Definition and methodology.— Whatis the when
response
one can -q
encounterq ⇠ k. In the Eikonal approximation thek' -q -q
k'-q k to -qof the nonlinear power spectrum at wavenumber
at wavenumber q? At linear+=
qthe linear power spectrum-q-q k'level, it is simply a Dirac-delta function since each-qFourier -q=+ & AT+ (’16,’17)
mode evolves+ independently in+ standard cosmological
Nishimichi, Bernardea
N体シミュレーションでは、
scenarios.
Here we wish to introduce a well-definedkernelrepresentation of the two-loopk'co
FIG. 9. Diagrammaticfunctionand investigate k'
it at fully nonlinear level. We
大スケールに対する小スケールの影響は抑制されている
FIG.a
9. Diagrammaticconsider the nonlinear power spectrum asfunctionalrepresentation of the two-loop co
nl
nl
lin
of the出力
linearpower
spectrum,
i.e.,
P
=
P
[P
and 入力 (線形)(非線形) Pnl (k) =d ln q K(k, q) P],0 (q)define the kernel function as its functional derivative:
ntation of the two-loop contribution
to the resummed
multipoint-1
nl
k=0.1025hMacδP (k; z).
(1)応答関数 K(k, q; z) = q lin
δP (q; z)

大規模構造の応答関数

2-loopWe omit the explicit dependence
on z from the arguments
2-loopin what follows. 減衰項入りモデル
The normalization for K is chosen suchthat a small variation in P nl is related to that of P lin as!
k=0.1025 h Mac-1
δP nl (k) =d ln q K(k, q)δP lin (q).
(2)

This relation provides us a simple way to measure the kernel function from simulations. In order to do so, we pre-

摂動論のトラブル：小まとめ
従来の摂動論にひそむ問題点
シングルストリームの破れ → 高次の摂動計算が破綻（大スケールでも破綻）対処方法✓ 高次の計算をしない
適用範囲が限定される✓ 有効場理論のアプローチ

Baumann et al. (’12)
Carrasco, Herzberg & Senatore (’12), …

圧力ゼロ流体からのずれを表す有効ストレステンソルを導入
N対シミュレーションで較正理論の予言性が失われる動機
シングルストリームを超える取り扱いは摂動計算で可能か？
計算精度や適用範囲は従来の方法に比べて改善するか？

１次元宇宙での考察

AT & Colombi (’17)

ポストコラプス摂動論

ラグランジェ描像にもとづく新
しい摂動論

適応フィルタリング

マルチストリームの影響を低減

両者の組み合わせで小スケールまでシミュレーションを再現

Zel’dovtion givesInanone-dimensional
exact solutioncase,for thetheso-called
dynamics
of mtionthegives
an exact
solution forthe dynamic
one-dimensional
so-called
approximadimensionalIncase,the so-called
Zel’dovich
approximabefore case,shell-crossing.
TheZel’dovich
Zel’dovich
solution
also
3.1 Post-co
before shell-crossing.
The ofZel’dovich
solutiongives forana exact
solutionforthedynamics
masssheet
es an exacttion
solutionthedynamics
of
mass
sheet
natural basis for the cold initial condition. Thea naturalbasis
forsolution
the cold also
initialcondition.
Thebasic form
beforeThe
Zel’dovich
provideshell-crossing.
Theshell-crossing.
Zel’dovich
solution
also
providesgiven byby
(Zel’dovich’70)
naturalbasis forthegivencold The
initialcondition.
solutionStartinisal basis for athe
cold initialcondition.
solutionis Thefollows.
dD
+
we
first
followgivenbyyx(q; τ ) = x(q;q + τψ(q)D++ψ(q)(τ ), D+ (τv(q;
τ
)
=
ψ(q)
)
=
q
),
v(q;
τ ) = ψ(q): 線形成長因子
シングルストリームdτZel’dovichsoldD+ (τ )dD+ (τ )の厳密解x(q;τ )),= q +Here,ψ(q)), thev(q;
)=
ψ(q)
(20)
the
function
D+τ.isthegrowthfactorfa= q + ψ(q)
D
v(q;
τ )D
=+ (τ
ψ(q)
(20)
shell-crossing,
Here,functionD
is thelinear. growth
+ (τ
: 変移場
+lineardτdτ
the following
equation:
the following
equation:compute the bD+» is 2thelinearsatisfying
he function Here,D+ isthethefunction
linearfactorsatisfying» 2 growth–growth– factor
anexplicit
fun
d
3
位相空間
d
3
2
2 Ω
the following equation:owing equation:
−
H
a(τ
)) D
(τthe
) =shell-crossi
0.
m,0
+0.
0
−
Ω
H
a(τ
)
D
(τ
=
m,0
+
2
0
2
» 2–
dτ 2
2 –dτ
シェルクロッシングd3exerted at eac
3v
2
2
that
in
terms
of
the
cosmic
time
t,
Eq.
Ωm,0
a(τ )Note
D
(τ
)
=
0.
(21)
+
Ωm,0 H0 a(τ )2 −
D+ (τ
)Note
=H0.
(21)
0 thatin terms of the cosmic timet, Eq.
(21)
from
Eq.
(19)
dτ
2
2
to the standard form of the linear evolution
to the standard
form of the linear evolution
equa
tain
the
correc»–2 cosmic
2 is reduced
Note
that
in»terms
oft,thetime
t,3 Eq.
(21)
at in termsofthecosmic
time
Eq.
(21)
is
reduced–ddΩ
Hm,0 0 from Eq. (16).
xx2x2+linear
2H(t)−H0 3 equation:
D+ (t) = 0.
d
d
3
Ω
m,0
to
the
standard
form
of
the
evolution
2
tandardform of the linear evolution
dt equation:dt
2 aD(t)+ (t) = 0.
+
2H(t)
−
密度プロファイル
2
3 (t)
» 2
–
dt
dt
2
a
–
2
2
d
d
3
Ω
H
m,0 0
d
3 Ωm,0H
1 The shell-cros
0
+
2H(t)
−
D
(t)
=
0.
(22)
シェルクロッシング
+
2H(t)
−
D
(t)
=
0.
(22)
+
2
3
MNRAS
000,
1–21
(2015)
δdt
2 a (t)dt
2dt a3 (t)from Lagrangia

ゼルドビッチ解

0.04

0.04

0.04

0.02

0.02

0.02

0.00

0.00

0.00

!0.02

!0.02

!0.02

!0.04

!0.04

!0.04

0.0

0.2

0.4

0.6

0.8

1.0

0.0

0.2

0.4

0.6

0.8

1.0

0.0

0.2

0.4

0.6

0.8

1.0

10

10

10

8

8

8

6

6

6

4

4

4

2

2

MNRAS 000, 1–21 (2015)

0

0

0

2

(2015)
000, 1–21(2015) 000, 1–21
x MNRAS
x
0.0

0.2

0.4

0.6

0.8

1.0

0.0

0.2

0.4

0.6

0.8

1.0

0.0

0.2

0.4

x

0.6

0.8

1.0

まで厳密解

Let us now consider a more rwith random initial conditions.
setup in 1D, a relevant initial co
3D case may be given by the G
with the initial power spectrum
初期条件 2 Planck ΛCDMk
P1D (k) =P3D (k)2π

シミュレーション vs ゼルドビッチ

with P3D being the matter po
無次元パワースペクトル
we computed with the transfer

は３次元と同じ

MNRAS 000, 1–21 (2015)

シングルストリーム近似の限界がクリアにL = 1, 000 Mpc# of particles (sheets) : 200, 000# of runs：50

パブリックコード： Vlafroid
http://www.vlasix.org/uploads/Main/froid1D.1.5.tar.gz

>:

rized in Eqs. (53)F (q0 , τ ) = −ψ(q0 ) D+ (τ ), (65) for inner partandandC2.are defined b
where the quantities A, B, andF (q0 , τ ) = −ψ(q0 ) D+ (τ ),(30). Note that in (41)deriving Eq. (39), we havesystemfollows
Zel’dovich solution at |Q|
where the quantities A, B, and C are thedefinedby Eqs.
(28)–
AT & Colombi (’17)
resultant
expressions
are written in term(30). Note that in deriving Eq. (39), thewe haveassumedthatquantities
densityVelocity
peak at anp
ゼルドビッチ解をもとに、シェルクロッシング後のb c . Since the3.3.1
thesystem follows Zel’dovich solution
at |Q| >characterizing
Qshell-crossing
time τ0 , Eq. (39) is still app
the resultant expressions are writtenthe
in termsof the local
Let us first conside
マルチストリーム領域を扱う新しい摂動計算手法
casesatinposition
which there
appear other shell-crossi
quantities characterizing the density peakqand0(ii). In this case, t
bc .
|Q|
>
Q
the shell-crossing time τ0 , Eq. (39) is still applicable to the
計算の概要
cases in which there appear other shell-crossing regions
at変移場∆vout (Q; τ, τbc ) =b|Q| > Qc .
3.3 Corrections to the Zel’dovich flow
; |Q| < Qc ,ポストコラプス摂動論

ゼルドビッチ解1.
2.

xZel (q; ⌧ ) = q + (q; ⌧ ) q：ラグランジュ座標

Provided the explicit expression for the force in
3.3 Corrections to the Zel’dovich
flowwe now compute
シェルクロッシング
region,the corrections to t3(q; ⌧ ) = A(q0 ; t) + B(q0 ; ⌧ ) (q q0 ) + C(q
;
⌧
)
(q
q
)
+ · ·formal·flow
solution
Eqs.
0 based on
0 theRecallinginthe
fact(1
Provided the explicit expression for the force in multi-stream 周りでテイラー展開
which
give
the
approximate
expression
relevant
above
integralsarregion, we now compute the corrections to the Zel’dovich
D [seeF (x(q;⌧ )) = Appendix
rx (x(q;⌧ ))E
マルチストリーム領域の「力」を計算：
flow based on the formal solution invalued
Eqs. region:(16)
andZ (17),τ
the ′first term, we
which give the approximate expression relevant at the multi-′∆v(Q; τ, τq ) =dτ F (x(Q,τ the)), short perio
forvalued region:「力」を積分してゼルドビッチ解に対する反作用を求める：
τq
Z τ
Z τ
Z τ
∆v(Q; τ, τq ) =dτ ′ F (x(Q, τ ′ )), ∆x(Q; τ, τq ) = (42)dτ ′ ∆v(Q; τ ′ , τq ).
dτ ′ a(τ ′ ) Jτqτqτbc (Q)Z τ"
on the position in Lag
ラグランジュ座標
Q=q-q(43)0 の7次の多項式∆x(Q; τ, τq ) =dτ ′ ∆v(Q;τ ′ , τq ). Notice that depending

Taruya and S. Colombi

10

A. Taruya and S. Colombi

ポストコラプス摂動：孤立ハロー

Taruya and S. Colombi

AT & Colombi (’17)
シミュレーション
ombi
Taruya and10S. Colombi
A. Taruya最初のシェル
and10S. Colombi
A. Taruya and10S. Colombi
A. Taruya and
10
S.
ColombiA.
Taruyaand
S. Colombi
２度目のク
ゼルドビッチ
ombi
mbi
位相空間
ポストコラプス摂動
クロッシングロッシング

密度プロファイル

napshots of phase-space structure (upper inset) and density profile (lower
Figureinset)
2. forSnapshots
the single-clusterof phase-spaceformation
structurein Einstein(upper inset) and density profile (lower inset) for theverse. For the initial density contrast given in Eq. (68), results of N -bodyde Sitter
simulations
universe.are
Fordepicted
the initial
as density
red lines,
contrast
while thegiven in Eq. (68), results of N -body simulation
lts with Zel’dovich solution are shown in green dotted lines. The blueanalytic
solid lines
results
are thewithprediction
Zel’dovich
withsolutionbasic arepost-collapse
shown in green dotted lines. The blue solid lines are th
Snapshots of phase-space structure (upper inset) and density profile (lower
Figureinset)
2. Snapshots
for the single-clusterof phase-spaceformation
structurein Einstein(upper inset) and density profile (lower inset) for thnt.
PTtreatment.
niverse. For the initial density contrast given in Eq. (68), results of Nde-body
Sittersimulations
universe. For
are thedepicted
initialasdensity
red lines,
contrast
while given
the in Eq. (68), results of N -body simulation
ults with Zel’dovich solution are shown in green dotted lines. The blue
analytic
solid lines
results
are withthe prediction
Zel’dovichwith
solutionbasicarepost-collapse
shown in green dotted lines. The blue solid lines are t
ent.
PT treatment.

Snapshots of phase-space structure (upper inset) and density profile (lower inset) for the single-cluster formation in Einsteiniverse. For the initial density contrast given in Eq. (68), results of N -body simulations are depicted as red lines, while theultswith Zel’dovich
solution
are shownin inset)
green anddottedlines.profile
The blue
solid
lines
arethethesingle-cluster
prediction
with
basic
post-collapse
napshotsof phase-spacestructure(upperdensity(lowerFigureinset)
2. forSnapshotsof phase-spaceformation
structurein Einstein(upper inset) and density profile (lower inset) for thnt.

２度目のクロッシング後の近似は悪くなるものの、密度プ
ロファイルの形状はシミュレーションをそれなりに再現

iverse. For the initial density contrast given in Eq. (68), results of N-bodyde Sitter
simulations
universe.are
Fordepicted
the initial
as density
red lines,contrast
while thegiven in Eq. (68), results of N-body simulation
ults with Zel’dovich solution are shown in green dotted lines. The blueanalytic
solid lines
results
are thewithprediction
Zel’dovichwith
solutionbasic arepost-collapse
shown in green dotted lines. The blue solid lines are tnt.
PT treatment.

ucture (upper inset) and density profile (lower inset) for the single-cluster formation in Einstein-Let us now consider a more rwith random initial conditions.
setup in 1D, a relevant initial co
3D case may be given by the G
AT & Colombi (’17)
with the initial power spectrum
Planck ΛCDM
2k
P1D (k) =P3D (k)2π位相空間with P3D being the matter po
we computed with the transfer

ポストコラプス摂動：CDM的初期条件（無次元）パワースペクトルz=0

LinearMNRAS 000, 1–21 (2015)

Zel’dovich(analytic)

シミュレーション
ゼルドビッチポストコラプス摂動

ハロー同士の相互作用の影
響大→摂動的記述が悪化

適応ポストコラプス摂動
ハロー個々の構造に興味はない → スムージング
ただし、ハローの成長・合体過程は環境に依存
初期密度ピークに対して「適応フィルタリング」
フィルターされた初期密度場を用いてポストコラプス摂動

v

peak

２度目のクロッシング時刻に対応する臨界密度値

Rs
大フィルタースケール

小x

適応ポストコラプス摂動
ハロー個々の構造に興味はない → スムージング
ただし、ハローの成長・合体過程は環境に依存
初期密度ピークに対して「適応フィルタリング」
フィルターされた初期密度場を用いてポストコラプス摂動

v

peak

２度目のクロッシング時刻に対応する臨界密度値

Rs
大フィルタースケール

小
x

適応

ポストコラプス摂動：CDM的初期条件
適応フィルタリング（無次元）パワースペクトル
Linearz=0

初期密度ピークに適用(with filter scales determinedby first-barrier crossing)
位相空間

シミュレーション
ゼルドビッチポストコラプス摂動

AT & Colombi (‘17)

考察
シングルストリームを超える取り扱いは摂動計算で可能か？
フリーパラメーターなし
ポストコラプス摂動論
マルチストリームのラグランジュ的取り扱い

適応フィルタリング
初期密度ピークに応じたハローの粗視化

３次元にも有効（なはず）
ただし、３次元ではゼルドビッチ解はあくまで「近似」
シェルクロッシング前の記述も正確にできるか非自明宇宙論的ヴラソフコードの発展

ophysical Journal, 762:116 (18pp), 2013 January 10

c)doi:10.1088/0004-637X/762/2/116d)American Astronomical Society. All rights reserved. Printed in the U.S.A.

b. 323 + two level dynamic adaptive refinement

DIRECT INTEGRATION OF THE COLLISIONLESS BOLTZMANN EQUATION
IN SIX-DIMENSIONAL PHASE SPACE: SELF-GRAVITATING SYSTEMS

2013

1

Kohji Yoshikawa , Naoki Yoshida

2,3, and Masayuki Umemura

1646

Figure 14. The initial conditions for the “ripple-wave” test problem
1 Center for Computational Sciences, University of Tsukuba, 1-1-1 Tennodai,
305–8577,
kohji@ccs.tsukuba.ac.jp(cf. Sec. 4.2).Tsukuba,
Shown areIbaraki
the particle
locations Japan;(panel a),the density
2 Department of Physics, The University
field
using
the
tetrahedral
phase
space
elements(panelb), usingofTokyo,
Tokyo
113-0033,
Japan
A
Mon. Not. R. Astron. Soc. 000, 000–000 (0000)
Printed 12 January 2015
(MN L TEX style file v2.2)
3 Kavli Institute for the Physics and Mathematics of the Universe,
tri-linear
elements(panelc)
andusing
tri-quadraticelements(panel JapanThe University of Tokyo, Kashiwa, Chiba 277-8583,
d). The linear elements are discontinuous at element boundaries,
Received 2012 June 18; accepted 2012 November 23; published 2012 December 20
while the quadratic is continuous.
c. 5123 N-body

An adaptively refined phase-space
ABSTRACTelement method forcosmologicalsimulations
andcollisionlessdynamics
We present
a scheme for numericalsimulations ofcollisionlessself-gravitating
systems which directly integrates thetri-quadratic reconstructed from N-body 323

Vlasov–Poisson equations in six-dimensional phase space. Using the results from a suite of large-scale numericalsimulations, we demonstratethat the present scheme
can simulate collisionless self-gravitating systems properly.
?12
Oliver
Hahn
and
Raul
E.
Angulo†The integration
schemeisETH
based
on the positive flux conservation method recently developed in plasma physics.
1 Department
case:
sine
waves (phasespaceevolution)of Physics,
Zurich, CH-8093 Zürich, Switzerland
We test2theaccuracy
ofdeourcode
by performing
severaltest1,calculations,the stability of King spheres, theCentrode Estudios
Fı́sicadel Cosmos
de Aragón, PlazaSan Juan
Planta-2, 44001,including
Teruel, Spain.
gravitational instability, and the Landau damping. We show that the mass and the energy are accurately conserved for
all the test cases we study. The results are in good agreement with linear theory predictions and/or analytic solutions.
submitted to MNRAS Jan. 8, 2015
The distribution function keeps the property of positivity and remains non-oscillatory. The largest simulations are
tri-quadratic 32 self-consistent
Figure 17. The ripple wave collapse test with dynamic adap
ColDICE:parallelsolverusing
moving adaptive simplicial
run on 646 grids. The computation speed scales
well withathe
number Vlasov-Poissonof processors, and thus
our code
performs
refinement. The 323 runs use the same initial conditions as in Fi
ABSTRACTefficiently on massively parallel supercomputers.
tri-quadratic elements and one (top, panel a), and two (middle, pFigure 16. Comparison between a reconstruction of the tritessellation
Dark matter numerical simulations
anddensity
the N field-bodymethodessential
for understanding
2 standardb) of dynamic adaptive refinement. The bottom panel shows
quadraticfromthe 32are
N -bodyrun (top
structureforms–andevolvesinthe and
Universe.
However, evolution
the discrete
natureof N -body solution of a high-resolution N -body run using 5123 particles a
Key words: galaxies: kinematicshow
anddynamics
methods:
numerical
half-panel)
the self-consistentof thetri-quadratic

Cold initial condition
2016

stro-ph.CO] 8 Jan 20153

Cold initial condition

2016

ep 2015

simulations can a↵ect its accuracy
modelling
collisionlesssystems.
same 2563 PM force resolution. On clearly sees how adding melementswhen(bottom
half-panel).
One clearly seesthat N -body particle
a,b,c,⇤
Online-only material: color figures
supportinga points approaches the high-resolution N -body solu
significantly
perturbs
thesolution,in particular,caustics
areWe introduce a new approachnoise
to simulatethegravitational
evolutionof cold
collisionlessThierry
Sousbie, Stéphane
Colombi
not persistent.
fluids by solving the Vlasov-Poisson
equations in terms of adaptively refineable “Lagrangian Still, the top two panels have significantly fewer degrees of free
than the N -body run.
phase space elements”. Thesea Institut
geometricalelementsdearepiecewisesmoothbetweend’Astrophysique
Paris,CNRS UMR
7095mapsand UPMC,
98bis, bd Arago, F-75014 Paris, France
three-dimensional Lagrangian space and six-dimensional
Eulerian phase space and apb Department of Physics,
The
Universityof Tokyo, Tokyo 113-0033, Japanusing
refinement
in
Figure17,
comparingonce
more
proximate the continuum structureofthedistributionfunction.
allow
foragainst
dynamicalc Research3 Center for the Early Universe,They
Schoolof Science,
The University
Tokyobut
113-0033,
Japan speak for themse
these oftoyTokyo,problems
let the images
thefollow
512
particle
high-resN -bodysolution
the
same
force
adaptive splitting to accurately
theonevolutioneven
ingravitational
regions
ofatverystrongmixing.
movethesmoothpotential
calculatedfromthe
1. INTRODUCTION
resolution. We only
considerdescription
the tri-quadraticelements in this cold and perform a quantitative convergence study of refinem
The elements thus permit a deterministic
non-linearof self-gravitating
overall
distributionoftheparticles.
Despiteofin these
attractive
the next
section, where we apply the Lagrangian elem
althoughlimit.
the linear elements also perform reasonably
well.
and collisionless fluids in the case,continuous
3
methodto cosmological
We started
withthree-dimensional
thesamemajor
32 initialconditions
as inwhichthefeatures,thedisadvantage
of fixed
theSCF
methodis its structure formation.
tational interaction is one of the mostimportantphysicalWe discussin detailvarious one-,
twoandtest
problems
demonresolution
testofshown
Figurebutnow
employed
theforce hasso that the lowest
inflexibility
thatthe15,
basis
setthatmust
chosenstrate the ofcorrectness
and performance
our inmethod.
We
show
ourbe
methodes in the dynamics and the formation
astrophysrefinement criterion with a threshold of 0.1 to dynamicallyseveral advantages
to standard
N -bodyalgorithmstheby i)global
explicitly
trackingof
thethe systems underorder terms
reproduce
ects, such as star clusters, galaxies,
and the comparedlarge- split
elements
if required(the results using
velocitystructurefine-grainedAbstract
distribution function,ii) naturally
representing
caustics,iii)refinement
providing an 5 A FIRST APPLICATION: COSMOLOGICAL
are however
not significantlydi↵erent).
The solution
allowing
investigation(Weinberg
1999).
otherwords,SCF methodOF A WARM DM UNIVERS
ucture of the universe. Stars andarbitrarily
dark matterthese field
regularindensity
thatis defined everywherein space,iv)Ingivingdirectlya theSIMULATION
for one additional
level of refinement
is shown
inthe top panel,
Resolving
numerically
Vlasov-Poisson
equationsfor
initiallycoldsystemscanbe reduced
can be applied
to thethesymmetric
gravitational
collapseor to following the evolution
vitating systems are essentially collisionless,exceptforsmooth and regulargravitational
potential
field, thusonly
eliminatingneed for any
type of

解析計算とシミュレーションの詳細な比較ができるようになってきた

the one for two levels in the middle panel, and the reference

We now apply our Lagrangian phase space element method

a=0.023

a=0.033

a=0.025

a=0.048

2D collapse with sinewave initial condition
Initial displacement @ a=0.01:✓
2⇡L
0.4 sin L qx
(q) =
2⇡
0.3
sin
2⇡
L qy

◆

projected
densitya=0.088a=0.255vx
yx

http://www.vlasix.org/index.php?n=Main.ColDICE

a=0.023

a=0.033

a=0.025

2D collapse with sinewave initial condition
Initial displacement @ a=0.01:✓
2⇡L
0.4 sin L qx
(q) =
2⇡
0.3
sin
2⇡
L qy

a=0.048◆

projected
densitya=0.088a=0.255
y

vx

vx

y
x

x

vy is color coded
http://www.vlasix.org/index.php?n=Main.ColDICE

３次元シェルクロッシングの記述
嵯峨くんの講演
ラグランジュ的摂動論にもとづく記述 （LPT・Q1D）(e.g., Matsubara ’15, Rampf & Frisch ’17)

ヴラソフシミュレーションとの比較0.8（準１次元的コラプス）

0.4

0.79

0.2

0.77

1
✓ ✏x sin qx2⇡
L ⌘ @ 0.4 sin L Aqx(1)(q)⌘)==aeinit ✏y sin qy2⇡
LPT (q,
0.3 sin L qy
2⇡✏z sin qz

0.6vx

変移ベクトル0

a =a =0.035
0.035

パラメーター：ainit = 0.0005(✏x , ✏y , ✏z ) = ( 24, 4, 3)

◆0
-0.2
-0.4
-0.6
-0.8

0.82

y=z=0

0.81
0.8
0.78
0.76
-0.12

0.8
0.6
0.4
0.2
0
-0.2
-0.4
-0.6
-0.8
-0.04-0.03-0.02-0.01 0 0.010.020.030.04

-0.4

-0.2

-0.1

-0.08

LPT 5thQ1D 1stQ1D 2nd
LPT growing 10th

0
x

0.2

0.4

-0.06

３次元シェルクロッシングの記述
嵯峨くんの講演
ラグランジュ的摂動論にもとづく記述 （LPT・Q1D）(e.g., Matsubara ’15, Rampf & Frisch ’17)

ヴラソフシミュレーションとの比較a =a =0.040
0.035
0.035
0.040
0.8
0.8
0.6
0.6

（準１次元的コラプス）
1✓ ✏x sin qx2⇡
L ⌘ @ 0.4 sin L Aqx(1)(q)⌘)==aeinit ✏y sin qy2⇡
LPT (q,
0.3 sin L qy
2⇡✏z sin qz

y=z=0

0.4
0.4
0.2

vx

変移ベクトル0

0.82
0.91
0.9
0.81
0.89
0.8
0.88
0.79
0.87
0.78
0.86
0.77
0.85
0.76
0.84
-0.12
-0.1

パラメーター：ainit = 0.0005(✏x , ✏y , ✏z ) = ( 24, 4, 3)

◆0
-0.2
-0.4
-0.6
-0.6
-0.8
-0.8

0.8
0.6
0.6
0.4
0.4
0.2
0.2
0
-0.2
-0.2
-0.4
-0.4
-0.6
-0.8
-0.6
-0.04-0.03-0.02-0.01 0 0.010.001
0.020.030.002
0.04-0.002-0.001

-0.4

-0.2

-0.1
-0.08

-0.08
-0.06

LPT 5th
5th
LPT
Q1D 1st
1st
Q1D
Q1D 2nd
2nd
Q1D
LPT growing
growing 10th
10th
LPT

0
x

0.2

0.4
0.4

-0.06
-0.04

３次元シェルクロッシングの記述
嵯峨くんの講演
ラグランジュ的摂動論にもとづく記述 （LPT・Q1D）(e.g.,proﬁleMatsubara& Frisch’17)
Density(εx, ’15,εy, εRampf
)
=
(-24,
-4,
-3)
z

ヴラソフシミュレーションとの比較

0.8
��
0.8
0.6
0.6

（準１次元的コラプス）
1✓ ✏x sin qx2⇡
L ⌘ @ 0.4 sin L Aqx(1)(q)⌘)==aeinit ✏y sin qy2⇡
LPT (q,
0.3 sin L qy
2⇡✏z sin qz

0.89
0.8
0.88
0.79
0.87
0.78
0.86
0.77
0.85
0.76
0.84
-0.12
-0.1

��
0.4
0.4

0.2
��

vx

変移ベクトル0

a =a =0.040
0.035
0.035
シェルクロッシングが起こる
0.040
0.82
0.91
a = 0.03814
直前の密度プロファイル
0.9
0.81
y=z=0

パラメーター：ainit = 0.0005(✏x , ✏y , ✏z ) = ( 24, 4, 3)

◆0
�

-0.2
-0.4
�
-0.6
-0.6
�

-0.8
-0.8

0.8
0.6
0.6
0.4
0.4
0.2
0.2
0
-0.2
-0.2
-0.4
-0.4
-0.6
-0.8
-0.6
-0.04-0.03-0.02-0.01 0 0.010.001
0.020.030.002
0.04-0.002-0.001

LPT full 4th
LPT full 5th
Q1D 1stQ1D 2nd
LPT growing 10th�� × ��-�

-0.4

-0.2

�� × ��-��� × ��-�

-0.1
-0.08

-0.08
-0.06

-0.06
-0.04

LPT 5th
5th
LPT
Q1D 1st
1st
Q1D
Q1D 2nd
2nd
Q1D
LPT growing
growing 10th
10th
LPT

0
x

����

0.2����

0.4
0.4

����

３次元シェルクロッシング後
シェルクロッシング後の１次元ポストコラプス摂動論と比較（2次のQ1Dに対する力の反作用を求めた）
ヴラソフシミュレーション
ポストコラプス摂動
a=0.045

パラメーター：ainit = 0.0005(✏x , ✏y , ✏z ) = ( 24, 4, 3)

a=0.060
a=0.080

準一次元的なコラプスでも、２度目以降のシェルクロッシングが早まり、１次元の摂動計算とずれる（多次元の効果）

まとめ
シングルストリーム近似にもとづく従来の計算手法をこえる大規模構造の摂動論の開発・発展
シングルストリーム近似にもとづく摂動計算の問題点
高次摂動で破綻

シェルクロッシングを超える（１次元）：• ポストコラプス摂動論
マルチストリームのラグランジュ的取り扱い• 適応フィルタリング
初期密度ピークに応じたハローの粗視化
３次元への拡張に向けた研究も進展中 → 嵯峨くんの講演