x - CfCA

プラズマ過程から迫る
高エネルギー
天体現象
東京大学RESCEU
茂山研D2 鈴木昭宏
自己紹介
鈴木昭宏
東京大学理学系研究科天文学専攻D2
ビッグバンセンター茂山研
専門は、プラズマ天体物理
特に非熱的放射の起源に興味がある
対象天体はいろいろ
→超新星とその残骸,GRBs,AGN,microQSO
Outline
PART I : プラズマ天体物理の基礎
1.Introduction
2.Vlasov-Maxwell system
3.Plasma Instability
4.Wakefield Acceleration
PART II : 最近の興味
1.Bulk comptonization
2.Supernova shock breakout
3.XRF080109/ SN 2008D
Outline
PART I : プラズマ天体物理の基礎
1.Introduction
2.Vlasov-Maxwell system
3.Plasma Instability
4.Wakefield Acceleration
PART II : 最近の興味
1.Bulk comptonization
2.Supernova shock breakout
3.XRF080109/ SN 2008D
Outline
PART I : プラズマ天体物理の基礎
1.Introduction
2.Vlasov-Maxwell system
3.Plasma Instability
4.Wakefield Acceleration
PART II : 最近の興味
1.Bulk comptonization
2.Supernova shock breakout
3.XRF080109/ SN 2008D
plasma......
plasma......
Collisionless plasma
構成粒子同士の衝突が無視できるプラズマ
Maxwell-Boltzmann分布からずれた速度分布(非熱的成分など)が実現さ
れる
宇宙の至る所に存在する(地球磁気圏,SNR,GRB,AGN jet,ICMなど)
3CALE ,ENGTH -EAN &REE 0ATH
%LECTRON
Collisionless shocks
無衝突プラズマ中で起こる衝撃波
shock frontでの散逸過程は粒子と電磁場との相互作用が担う
(cf. 流体衝撃波での散逸過程は粒子衝突による粘性)
粒子加速のサイトと考えられている
Ex.) 地球磁気圏でのBow shock, SNR, GRB, PWN, AGN Collisionless shocks
Bow shock, SNR : non-relativistic
GRB shock : Γ∼100
AGN jet : Γ∼10
Crab Nebula : Γ∼10
6
Plasma parameters
Plasma frequency
Inertial length
Gyrofrequency
Outline
PART I : プラズマ天体物理の基礎
1.Introduction
2.Vlasov-Maxwell system
3.Plasma Instability
4.Wakefield Acceleration
PART II : 最近の興味
1.Bulk comptonization
2.Supernova shock breakout
3.XRF080109/ SN 2008D
Plasma kinematics
プラズマ= 荷電粒子の集まり + 電磁場
特に、荷電粒子間の相互作用はdebye shieldingによって無視できる
collisionも無視することで、粒子は電磁場中の運動方程式に従う
電磁場はMaxwell方程式で決まる
dp
p
=q E+
B
dt
mγ
· E = 4πρ
·B =0
1 B
c t
1 E
B=
+ 4πj
c t
E=
: particle
Debye Shielding
v (θ) =
2Eexp
2 ∆r ρ
4πrin
in in
in
プラズマ= 荷電粒子の集まり + 電磁場
点電荷のポテンシャルが逆の電荷を持つ粒子の遮蔽によって、変更され
Eexp
る
q
φ=
r
デバイ長λDまでしか、粒子間相互作用は効かない
但し、場からは外力を受ける
vin (θ) =
2Eexp
2 ∆r ρ
4πrin
in in
Eexp = 1051 erg
dp
p
=q E+
B
dt
mγ
· E = 4πρ
·B =0
1 B
c t
1 E
B=
+ 4πj
c t
E=
: particle
q
φ=
r
51
= 10 erg
q
φ = ex
r
1 + α cos(2θ)
1 − 2α/3 + 7α2 /1
ρin
q
r
φ = exp −
r
λD
Fluid vs plasma
流体 'HQVLW\
6KRFNIURQW
Q
7
(OHFWURQ
,RQ
X
YHORFLW\
無衝突プラズマ 'HQVLW\
6KRFNIURQW
(OHFWURQ
,RQ
無衝突プラズマの扱いでは速度(位相)空間が必要
Vlasov-Maxwell system
me (vy2 + vz2 )
me vx2
f (vx , vy , vz ) ∝ exp − 2
exp −
2 T2
T
2k
m
(v
m2k
v
B
⊥
B vz )
e
e x
y +
f (vx , vy , vz ) ∝ exp −
exp −
plasmaの構成粒子の位相空間
2kB T⊥
2kB2T
2
me (vx + V )
me (vx − V )
での分布関数の時間発展を記
f (vx , vy , vz ) ∝ exp −
+ exp −
exp −
2
2
2k
T)
2k
T)
B
BV
m
(v
+
V
m
(v
−
e
x
e
x
述する方程式
eq.
f (vx , vy , vz ) ∝ expVlasov
−
+ exp −
exp −
2kB T m∂f
2
2
∂fj 2kB T∂fj
Zj e
me v1x2
e (vjy + vz )
f(x,p,t): 位相空間中の微小要素
+f (v
vx·, vy , vz+
=0
E+ v×
) ∝ exp −
expB − ·
2k
T
2k
T
∂t
∂x Z emj
∂vB
B c⊥
∂f
∂f
∂f
1
j
j
j
j
[x,x+Δx] [p,p+Δp]に何個粒
+v·
+
=0
E+ v×B ·
2
∂t
∂x me (vm
c
∂vV )2
j+ V )
m
(v
−
x
e
x
子がいるか
eqs.+ exp −
f (vx , vy , vz ) ∝ exp Maxwell
−
exp −
2kB T
2kB T
∇ · E =4πδ
Source termsを介して粒子と
∂fj
∂fj
Zj e
∂fj
1
+v·
+1 ∂E E + v × B ·
=0
∂t
∂x
m
c
∂v
電磁場が相互作用する
j + 4πj
∇×B=
c ∂t
δ=
dp
p
=q E+
B j
dt
mγ
·B =0
1 B
E=
c t
1 E
B=
+ 4πj
c t
: particle
j
E =4πδ
∞ ∇ ·∞
∞
∞ ∞ ∞∞
∂E v)dvx dvy dvz
Zj e
fj1(x,
∇
×
B
=
+ 4πj
Source
−∞ terms
−∞ −∞ c ∂t
Zj δe =
· E = 4πρ
j=
∞
j
−∞
j=
Zj e
−∞
−∞
Zj e
j
∞
fj (x, v)dvx dvy dvz
vf
(x,
v)dvx dvy dvz
j
−∞ −∞
−∞
−∞
∞
∞
∞
−∞
−∞
vfj (x, v)dvx dvy dvz
Outline
PART I : プラズマ天体物理の基礎
1.Introduction
2.Vlasov-Maxwell system
3.Plasma Instability
4.Wakefield Acceleration
PART II : 最近の興味
1.Bulk comptonization
2.Supernova shock breakout
3.XRF080109/ SN 2008D
Plasma oscillation
e-p plasmaのシート
electronを動かす
静電場が発生
復元力が働き、振動が励起される
(OHFWULFILHOG
(OHFWURQ
(OHFWURQ
,RQ
,RQ
(OHFWURQ
(OHFWURQ
,RQ
,RQ
Counter-stream
As the initial condition, we consider a counter-streaming plasma that is homogeneous
in space;
fj0 =
nj
2π 3/2 vj3
vx2 + vy2 + (vz + Vj )2
vx2 + vy2 + (vz − Vj )2
exp −
+ exp −
2
vj
vj2
,
(14)
where nj , vj , and Vj are the number density, the thermal velocity, and the bulk velocity of
IY
species j. We assume that they are constants and ni = ne = n0 for the charge neutrality.
In our simulation, the values of these parameters are as follows; n0 = 1, ve = 0.05, vi =
0.05 me /mi , and Ve = Vi = 0.2. The mass ratio is assumed to be mi /me = 1 and 16. The
initial configuration of the electromagnetic field is
9]
9 x sin y,
φ = Ax = Ay = 0,9 Az = − sin
(15)
which is equivalent to
Ex = Ey = Ez = 0,
Bx = sin x cos y,
By = − cos x sin y,
Bz = 0,
(16)
where we have introduced a small parameter (= 10−5 ). In other words, we treat the above
magnetic field as a perturbation to the initial distribution of particles with no electromagnetic
field. Next, we address the boundary conditions. The simulation domain consists of the
spatial intervals given by x, y ∈ [−π, π] and the velocity intervals given by vx , vy ∈ [−0.4, 0.4]
and vz ∈ [−0.6, 0.6]. The periodic boundary condition is imposed in the x and y direction,
Two-stream Instability
counter-streaming plasma
%X
%LECTRON %LECTRON longitudinal electric ﬁeld のperturbation
XAXIS
me (vy2 + vz2 )
me vx2
%LECTRON , vz ) ∝ exp
−
exp
−
(1)
ある方向へ進むelectronが 2kB T⊥
2kB T
「渋滞」する
me (vy2 + vz2 )
me (vx + V )2
me (vx − V )2
(2)
+ exp −
exp −
2kB T
2kB T
2kB T
∂fj
Zj e
∂fj
1
+v·
+
=0
E+ v×B ·
∂x
mj
c
∂v
Maxwell eqs.
%X
∇ · E =4πδ
Zj e
j
Zj e
∞
∞
−∞
−∞
−∞
∞
−∞
∞
−∞
∞
−∞
fj (x, v)dvx dvy dvz
%LECTRON %LECTRON (4)
XAXIS
%LECTRON %LECTRON (5)
vfj (x, v)dvx dvy dvz
%LECTRIC FIELD
(3)
1 ∂E
∇×B=
+ 4πj
c ∂t
∞
%LECTRON Two-stream Instability
counter-streaming plasma
%X
%LECTRON %LECTRON longitudinal electric ﬁeld のperturbation
XAXIS
me (vy2 + vz2 )
me vx2
%LECTRON , vz ) ∝ exp
−
exp
−
(1)
ある方向へ進むelectronが 2kB T⊥
2kB T
「渋滞」する
me (vy2 + vz2 )
me (vx + V )2
me (vx − V )2
(2)
+ exp −
exp −
2kB T
2kB T
2kB T
∂fj
Zj e
∂fj
1
+v·
+
=0
E+ v×B ·
∂x
mj
c
∂v
Maxwell eqs.
%X
∇ · E =4πδ
1 ∂E
∇×B=
+ 4πj
c ∂t
Zj e
j
Zj e
∞
∞
∞
−∞
−∞
−∞
∞
−∞
∞
−∞
∞
−∞
fj (x, v)dvx dvy dvz
(3)
(4)
ELECTRIC CURRENT
(5)
vfj (x, v)dvx dvy dvz
%LECTRON %LECTRIC FIELD
Weibel Instability
z-axis
Electron 2
conter-streaming plasma
y-axis
Electron 1
磁場のperturbationが存在
me (vy2 + vz2 )
me vx2
−
exp −
z ) ∝ expﬁlamentの形成
2kB T⊥
2kB T
(vx + V )2
me (vx − V )2
+ exp −
2kB T
2kB T
exp −
∂fj
Zj e
∂fj
1
v·
+
=0
E+ v×B ·
∂x
mj
c
∂v ZAXIS
Maxwell eqs.
∇ · E =4πδ
1 ∂E
∇×B=
+ 4πj
c ∂t
Zj e
Zj e
∞
∞
∞
−∞
−∞
−∞
∞
−∞
∞
−∞
∞
−∞
(1)
Electron 3
me (vy2
+
2kB T
vz2 )
: Electron
(2): Magnetic field
: Orbit
(3)
%LECTRON (4)YAXIS
%LECTRON fj (x, v)dvx dvy dvz
(5)
%LECTRON vfj (x, v)dvx dvy dvz
%LECTRON XAXIS
Electron 4
x-axis
Weibel Instability
z-axis
Electron 2
conter-streaming plasma
y-axis
Electron 1
磁場のperturbationが存在
me (vy2 + vz2 )
me vx2
−
exp −
z ) ∝ expﬁlamentの形成
2kB T⊥
2kB T
(vx + V )2
me (vx − V )2
+ exp −
2kB T
2kB T
exp −
(1)
Electron 3
me (vy2
+
2kB T
vz2 )
∂fj
Zj e
∂fj
1
v·
+
=0
E+ v×B ·
∂x
mj
c
∂v ZAXIS
Maxwell eqs.
∇ · E =4πδ
Zj e
Zj e
∞
∞
−∞
−∞
−∞
∞
−∞
∞
−∞
∞
−∞
fj (x, v)dvx dvy dvz
(2): Magnetic field
: Orbit
(3)
(4)YAXIS
1 ∂E
∇×B=
+ 4πj
c ∂t
∞
: Electron
%LECTRIC CURRENT
(5)
vfj (x, v)dvx dvy dvz
XAXIS
Electron 4
x-axis
me (vy + vz )
me vx
f (vx , vy , vz ) ∝ exp −
exp −
2kB T⊥
2kB T
Vlasov Equation (2D3V)
me (vx + V )2
me (vx − V )2
f (vx , vy , vz ) ∝ exp −
+ exp −
exp −
2kB T
2kB T
Vlasov eq.
splittin schemeを用いると、
∂fj
∂fj
Zj e
∂fj
1
∂f
∂f
+v·
+
v×B ·
=0
E + (16)
∂f
+ vx ∂f
=0
vlasov方程式は5本の1次元移
∂t
∂x
mj
c (16) ∂v
+ ∂x
vx
=0
∂t
∂t
∂x
流方程式に分解される
∂f
∂f
∂f
+ vy ∂f
=0
+ ∂y
vy
=0
∂t
1次元移流方程式なら色んな解
∂t
∂y
∂f
∂f
き方がある
∂f= 0
∂f
+ (Ex + vy Bz )
+ (Ex + vy B∂v
=0
∂t
z) x
→僕が使っているのは
∂t
∂vx
∂f
2nd-order
van Leer ∂f
scheme
∂f
∂f= 0
+ (Ey − vx Bz )
+ (Ey − vx B∂v
=0
∂t
z) y
∂t
∂vy
∂f
∂f
∂f
∂f= 0
+ (vx By − vy Bx )
=0
+ (vx By − vy B∂v
∂t
x) z
∂t
∂vz
∂g
∂g
∂g
+ aζ ∂g
=0
+ ∂ζ
aζ
=0
∂t
∂t
∂ζ
g(t + ∆t, ζ − aζ δt) = T [ζ, ∆t]g(t, ζ)
g(t + ∆t, ζ − aζ δt) = T [ζ, ∆t]g(t, ζ)
(17)
(17)
∂f
∂f
(18)
+ vx (18)
=0
∂t
∂x
∂f
∂f
+ vy (19)= 0
∂t
∂y(19)
∂f
∂f
+ (Ex + vy(20)
Bz )
=0
∂t
∂vx
(20)
∂f
∂f
+ (Ey − vx(21)
Bz )
=0
∂t
∂vy
(21)
∂f
∂f
(22)
=0
+ (vx By − vy B(22)
x)
∂t
∂vz
∆t
∆t
∆t
∆t
∆t
∆t
]T
[y,
]T
[x,
]T
[v
,
]T
[v
,
]T
[v
,
]T
∆t
∆t
∆t
∆t
∆t
∆t [vz , ∆t]
x
y
x
4
2
4
4
2
4
∂g
∂g
[x,
]T
[y,
]T
[x,
]T
[v
,
]T
[v
,
]T
[v
,
x 4
y 2
x 4 ]T [vz , ∆t]
(23)
4
2
4
+ aζ (23)
=0
∆t
∆t
∆t
∆t
∆t
∆t
∂ζ
× T×[vTx[v
, 4, ]T
v)
∆t[vy , 2 ]T
∆t[vx , 4 ]T
∆t[x, 4 ]T
∆t [y, 2 ]T
∆t [x, 4 ]f
∆t(t, x,∂t
]T
[v
,
]T
[v
,
]T
[x,
]T
[y,
]T
[x,
]f
(t,
x,
v)
x 4
y 2
x 4
4
2
4
g(t + ∆t, ζ − a δt) = T [ζ, ∆t]g(t, ζ)
∂f
∂f= 0
+ (Ey − vx Bz )
(19)
+
(E
−
v
B
)
=
0
(19)
∂t
y
x ∂v
z y
∂f
∂f
∂t
∂vy
+ vx
=0
∂t
∂x
∂f
∂f
∂f
∂f= 0
(20)
+ (vx By − vy Bx )
∂f
∂f
=0
(20)
+ (vx By − vy B∂v
∂t
x) z
+
v
=
0
y
∂t
∂vz
∂t
∂y
∂g
∂g
∂g
+ aζ ∂g
=0
(21)
∂f
∂f
+
a
=
0
(21)
∂t
∂ζζ
+ (Ex + vy Bz )
=0
∂t
∂ζ
∂t
∂vx
g(t + ∆t, ζ − aζ δt) = T [ζ, ∆t]g(t,
(22)
∂f ζ) ζ)
∂f
g(t + ∆t, ζ − aζ δt) = T [ζ, ∆t]g(t,
(22)
+ (Ey − vx Bz )
=0
∂vy
t
∆t
∆t
∆t
∆t ∂t
∆t
]T
[y,
]T
[x,
]T
[v
,
]T
[v
,
]T
[v
,
]T
[v
,
∆t]
2 ∆t ]T [x,
4 ∆t ]Tx[vx4, ∆t ]Ty[vy2, ∆t ]Tx[vx4, ∆t ]Tz[vz , ∆t]
, ∆t
]T
[y,
(23)
4
2
4
4
2
4
∂f
∂f
(23)
Vz依存性をΔt解く
∆t
∆t
∆t
∆t = 0
+]T
(vx[y,
By∆t
− ]T
vy B
x)
× T×[vTx[v
, ∆t
]T
[v
,
]T
[v
,
]T
[x,
[x,
]f
(t,
x,
v)
∆t
∆t
∆t
∆t
∆t
∆t
y
x
4
2
4
∂v
x4, 4 ]T [vy2, 2 ]T [vx4, 4 ]T∂t[x,
4 ]T [y, 2 ]T [x, z4 ]f (t, x, v)
Vlasov Equation (2D3V)
∂g
∂g
+ aζ
=0
∂t
∂ζ
x,y依存性をΔt/2だけ解く Vx,Vy依存性をΔt/2だけ解く
×T
[vx , ∆t
4 ]T
Vx,Vy依存性をΔt/2だけ解く
[x,
∆t
4 ]T
(18)
(19)
(20)
(22)
∆t
∆t
∆t
∆t
∆t
f (t + ∆t, x, v) =T [x, ∆t
]T
[y,
]T
[x,
]T
[v
,
]T
[v
,
]T
[v
,
x
y
x
4
2
4
4
2
4 ]T [vz , ∆t]
[vy , ∆t
2 ]T
(17)
(21)
g(t + ∆t, ζ − aζ δt) = T [ζ, ∆t]g(t, ζ)
[vx , ∆t
4 ]T
(16)
[y,
∆t
2 ]T
[x,
∆t
4 ]f (t, x, v)
x,y依存性をΔt/2だけ解く
ここで、Maxwell eqs.を解いている
(23)
S[x, y, ∆t] ≡ T [x, ∆t/2]T
[x,
∂fj [y, ∆t]T
∂fj
Zj e ∆t/2].
1
∂t
+v·
∂x
+
mj
∂fj
=0
E+ v×B ·
c
∂v
(10)
Maxwell Equations
· E =4πδterms in Equations (3) are
If the distribution function at t + ∆t is given, the ∇source
1 ∂Ethe time evolution of the
calculated by Equation (6). The wave equations
(3)
governing
∇
×
B
=
+ 4πj
Source terms
c ∂t
electromagnetic
are solved
by using 2the 2Fourier transformation,
which imposes the
∞
∞
∞
2
Sourcefields
termsは右の式から計
m
(v
+
v
)
me vx
e y δ =z
Zj e
fj (x, v)dv
x dvy dvz
f (vx , vboundary
−
exp
−
(1)
periodic
conditions
in
the
x
and
y
directions
implicitly;
y , vz ) ∝ exp
−∞ −∞ −∞
算する
∝
2kB T⊥
2kB T
Lx /2
exp
j
∞
Ly /2
∞
∞
j=
Zj e
2
2 vfj (x, v)dvx dvy dvz
2
2
m−∞
+ v−∞
ˆ Veqs.はFFTでフーリ
me (vx +
)kx , ky ) =
me (vx − V )exp[−i(k
e (v
−∞
Maxwell
yy)]h(t,
z ) x, y)dxdy,
h(t,
k
j xx +
y
−
(2)
+ exp −
exp −
−L
/2
−L
/2
x 2kB Ty
2kB T
2kB T
(11)
エ変換して、4th-order
Runge
∂fjKutta
∂fj
where∂fLjx (L
the Z
length
of 1the computational
domain in the x(y) direction and h is a
je
y ) is
+v·
+
=0
(3)
E+ v×B ·
∂t of (t, x,
∂xy). The
mj inversec Fourier transformation
∂v
function
is
defined
by
Maxwell eqs.
∇ · E =4πδ
Ly /2
Lx /2
1
ˆ kx ,(4)
h(t, x, y) =
exp[i(kx x + ky y)]h(t,
ky )dxdy.
(12)
12∂E
∇ × B(2π)
= Lx L+y 4πj
−Lx /2 −Ly /2
c ∂t
∞
∞
∞
These transformations
are computed
by the standard Fast-Fourier-Transformation scheme
δ=
Zj e
fj (x, v)dvx dvy dvz
(see, e.g. Press
et−∞
al 2007)
numerically. Substitution of these relations for φ, A, ρ, j into
−∞ −∞
j
(5) equations;
∞ following
∞
Equation (3) leads ∞to the
four second-order ordinary differential
j=
vfj (x, v)dvx dvy dvz
Zj e
j
−∞
2 ˆ −∞
dφ
=
2
dt
−∞
−(kx2 +
ky2 )φˆ − 4π ρˆ,
ˆ
d2 A
2
2 ˆ
ˆj,
=
−(k
+
k
)
A
−
4π
x
y
dt2
(13)
ˆ
which are solved by the Runge-Kutta method of order 4. To ensure that the obtained φ(t)
e
vy , vz−
)+∝4πj
expy − z
exp −
f (vx , vy , vz ) ∝ exp ∇
− × Bfx(v
=x ,exp
2k
T
B
⊥
2kB T⊥ c ∂t
2kB T
Procedure
f (vx , vy , vz ) ∝
∞
∞
2kB T
∞
me (vx + V )22
me (v
Vv)22)
x2 −
m
(v
+
m
(v
+
V
)
m
(v
−
V
)
δ =e f (v
, vjye, vz ) ∝ exp −2e fxj (x, v)dv+
exp
−z 2 e y
x xZ
z
2y dv
x dv
m
(v
+
v
)
m
v
exp −
+
exp
−
exp
−
2k
T
2k
T
e
e
y
z
B
B
x
−∞ −∞ −∞
2
f (vx , vy2k
,jvBz T
) ∝ exp −
2kB T
me (vy2 + vz2 )
exp −
(2) 2kB T
∞
∂fj
Zj e 2kB T 1
∂fj
j
+ vvf
· (x,∂f
+v)dv
v
×
B
·
=0
E
+
∂fjj = ∂fZjj e Zj e
1
jm x dvy dv
j ·
∂t v × B ∂x
cz
∂v 2
j 0
+v·
+
=
(3)
E
+
2
2 −∞ 2 −∞
22
2
m
(v
+
v
−∞
∂t
∂x
m
c
∂v
mj e (vx + V j) me vx
me (v
−yV+) vz )
e y
z)
mxe (v
f (vx , vy , vz ) ∝ fexp
− − ∇ · E =4πδ exp −
(vx , v−y , vz ) ∝ exp − + exp exp
(1)
2k
T
2k
T
2k
T
B
B
2k
T
2k
T
1. x積分(1/2ステップ)B∇ · E∂f=4πδ
B
j B ⊥ ∂fj
+v·
= 0∇ × B = 1 ∂E + 4πj
(4)
∂E ∂x 1
2
2
c
∂t
2∂tZ1j e
2 ∂f
∂fj me (vx ∇
∂f
j
j
m
(v
+
v
)
+×
V )B =
me (vx −BV ) ·
e y
z
++
4πj
E
∞
∞ = 0−
f (vx , vy , vz ) ∝ exp ∂t
− + v · ∂x + m
(2)
+c exp
− c v ×∞
exp
∂t
∂v
j δ=
2kB T
2k
T
2k
T
B y dvz
Zj e B
fj (x, v)dvx dv
∞
∞
2. 電荷・電流密度を計算
3.
4.
2k∂fB T∞
⊥
∞
2kexp
BT −
∞
−∞
−∞
+v
·
−∞
−∞
∂t
∂x
(1)
(5)
=0
feed back
電磁場 E,B
分布関数 f
feed back
(2)
(3)
(2)
(6)
(4)
(3)
−∞
j v)dv dv dv
δ=
Zj e
f
(x,
j
x
y
z
∇
·
E
=4πδ
∂fj
∂f
Z
e
1
∞ ∂f
∞
∞
j
j
j
j + v · −∞ +−∞ −∞
= 0 vf (x, v)dv dv dv
(3)
Ej =
+ vZ×eB ·
j
j
x
y
z (5)
∂t
∂x
m
c
∂v
∞
∞ j ∞ 1 ∂E
−∞ −∞ −∞
j + 4πj
∇ × B =vfj (x,
j=
Zj e
v)dv
x dvy dvz
c
∂t
Maxwell eqs.
· E =4πδ
−∞ ∇−∞
−∞
j
∞
∞
∞
(4)
1
∂E
δ=
Zj e∂f
f4πj
j (x, v)dvx dvy dvz
∇j ×+B
=∂fj = +
v
·
0
(6)
分布関数 f feed back
c ∂t
−∞
j
∂t−∞ −∞
∂x
∞∞ ∞ ∞ ∞ ∞
v積分
∂f
Z
e
∂fv)dv
1
j
δj== j +ZZ
j eje E + v × Bfj (x,
x0dvy dvdv
z dv feed back 電磁場(7)
·
= v)dv
E,B
vf
j −∞ −∞ −∞
j (x,
x
y
z
∂tj
mj −∞ c−∞ −∞ ∂v
j
(5)
∞
∞
∞
j=
Zj e
vfjj (x, v)dvx dvy dvz分布関数 f feed back
∂fj
∂f
5. x積分 (1/2ステップ)
−∞
j
(1)
(1)
(5)
(4)
(5)
(6)
me (vy2 + vz2 )
me vx2
f (vx , vy , vz ) ∝ exp −
exp 2−
2
2
2kBm
T⊥
vz ) 2kApJ
BT
mSuzuki&Shigeyama
e (vy +(2009)
e vx
f (vx , vy , vz ) ∝ exp −
exp −
2kB T⊥ 2
2kB T
me (vx + V )
me (vx − V )2
f (vx , vy , vz ) ∝ exp −
+ exp −
exp −
2kB T m (v − V )2
2kB T me (v 2 + v 2
me (vxVlasov
+ V )2 eq.
e x
y
z
物理空間2次元,速度空間3次元
f (vx , vy , vz ) ∝ exp −
+ exp −
exp
−
2
2
2
me (vy + v2k
m
v
e
z ) BT
x
2kB T ∂fj f (vx , vy ∂f
2k
T
B
Zj e−
1 exp − ∂fj
, vzj) ∝ exp
のVlasov-Maxwell system
2k
T
+v·
+
EB+⊥ v × B · 2kB T = 0
∂v
∂fj
∂fj ∂t Zj e ∂x 1 mj2
∂fc j
2
+
v
·
+
v
×
B
·
=
0
E
+
m
m
(v
+
V
)
m
(v
−
V
)
e
x
e
x
Vlasov eq.とMaxwell eqs.を
+ exp∂v
−
exp −
∂tf (vx , vy , v∂x
mj −
c
z ) ∝ exp
2kB T
Maxwell eqs.2kB T
交互に時間発展させていく
∇ · E∂f=4πδ
∂fj
Zj e
∂fj
1
j
Simulation
∂t
Source termsを介して粒子と
電磁場が相互作用する
∇×B=
δ=
Zj e
−∞
j
j=
Zj e
j
∞
∞
−∞
+v·
1 ∂E∂x
c ∂t
∞
∞
−∞
−∞
Source
δ∞=
−∞
j=
j
mj
+ 4πj
=0
E+ v×B ·
c
∂v
∇ · E =4πδ
1 ∂E
∇
×
B
=
+ 4πj
fj (x, v)dvxc dv
terms
∂ty dvz
Zj e
∞
∞
∞
∞
fj (x, v)dvx dvy dvz
vfj−∞
(x, −∞
v)dv−∞
x dvy dvz
−∞
Zj e
j
+
∞
∞
∞
−∞
−∞
−∞
vfj (x, v)dvx dvy dvz
∂y
∂f
∂f ∂t
=0
+– (v
x–By − vy Bx )
7
∂f
∂f
∂t 2
∂v
2
2
2
2 z+ (E + v 2 B )
=
vx + vy + (vz + Vj )
vx + vy + (vz −
x Vj )y z
nj
∂t
∂vx (14
+
exp
−
,
fj0 = 3/2 3 exp −
∂g
∂g
Suzuki&Shigeyama
2π vj
vj2
+ aζ
= 0 vj2 (2009) ApJ
∂f
2.4. Simulation
setups
∂t
∂ζ ∂f
+ (Ey − vx Bz )
=
∂t∆t]g(t,
∂vy o
where nj , vj , and Vj are the number density,
theζ thermal
and the
g(t + ∆t,
− aζ δt) velocity,
= T [ζ,
ζ) bulk velocity
As
counter-streaming
plasma
that isneutrality
homoge
Plasmaの初期条件：
species
j. the
We initial
assumecondition,
that they we
are consider
constantsa and
ni = ne =∂f
n0 for
the charge
∂f
∆t
∆t
∆t
∆t
∆t
∆t
f (t + ∆t,the
x, v)
=T
[x,
]T [vas
[v
,ev4=
]T0.05,
=
+,n(v
B1,
−
) z , ∆t]
space;
x , follows;
x]T
yx
yB
x[v
In in
our
simulation,
values
of these
parameters
v
v
4 ]T [y,
2 ]T [x, 4 are
4 ]T [vy
一様なconter-stream
plasma
02 =
i =
∂t
∂vz
∆t
∆t
∆t
∆t = 1 and
∆t 16. The
∆t
0.05 me /mi , and Ve = Vi = 0.2.2 The
ratio
is
assumed
to
be]Tm
i /me ]T
×
[v
,
]T
[v
,
]T
[v
,
[x,
[y,
]T
[x,
2 Tmass
2
2
2
2
x
y
x
4 Vj )
2
4
vx + vy + (vz +
v4x + vy∂g
+4(vz −∂g
Vj2)
nj
initial configuration
field is + exp −5−
expelectromagnetic
−
fj0 = 3/2 3of the
+
=0 ,
2
2 aζ
2π vj
vj
∂tvj
∂ζ
= 10
φ = Ax = Ay = 0, Az = − sin g(t
x sin+y,∆t, ζ − aζ δt) = T [ζ, ∆t]g
(15
ne =number
ni = 1 density,
ve = 0.05
me /mi and
Ve =the
Vi bulk
= 0.2veloc
where nj , vj , and Vj are the
thevthermal
i = 0.05 velocity,
∆t ni =∆t
∆t
∆t
species
j. We assume
that
arex,constants
ne]T=[x,n0∆tfor
the
charge
neut
電磁場の初期条件
which
is equivalent
to
f (tthey
+ ∆t,
v) =T [x,and
]T
[y,
]T
[v
,
]T
[v
,
x
y
4
2
4
4
2
In our simulation, the values of these parameters are as follows;
n0 ∆t
= 1, ve =
∆t
∆t 0.05,
×
T
[vcos
]T y,
[vy , B2 ]T
[v
[x
x , x4 sin
x , 4 ]T(16
E
=
E
=
E
=
0,
B
=
sin
x
cos
y,
B
=
−
=
0,
x i , and
y
x
z e = 1 and 16
0.05 me /m
Vez = Vi = 0.2.
The mass ratio yis assumed to be mi /m
initial configuration of the electromagnetic field −5
is
= 10−5
where we have introduced a small parameter (= 10 ). In other words, we treat the above
magnetic field as a perturbation
toA
the
initial distribution
of particles
with no electromagnetic
φ=
Az = −
sin x sin y,
x = Ay = 0,
field. Next, we address the boundary conditions. The simulation domain consists of the
spatial
given by
whichintervals
is equivalent
to x, y ∈ [−π, π] and the velocity intervals given by vx , vy ∈ [−0.4, 0.4
and vz ∈ [−0.6, 0.6]. The periodic boundary condition is imposed in the x and y direction
= Ez =
Bx = sinfunction
x cos y, assumed
By = −tocos
x sin for
y, |vBx |,z |v
=y0,
x = Ey space,
while in the E
velocity
the0,distribution
vanish
| > 0.4
and |vz | > 0.6. The number of zones in the physical space
is 32 × 32 and that in the velocity
−5
where
we×have
a small parameter (= 10 ). In other words, we treat the
space
is 40
40 ×introduced
60.
in space;
Conditions
Result
初期条件は右図のような磁場
の擾乱を与える
x,y方向には周期境界条件
Nx=Ny=32
Nu=Nv=40
Nw=60
Suzuki&Shigeyama (2009) ApJ
Result
初期条件は右図のような磁場
の擾乱を与える
x,y方向には周期境界条件
Nx=Ny=32
Nu=Nv=40
Nw=60
Suzuki&Shigeyama (2009) ApJ
–8–
–8–
Result
-6*157
ies j in each direction is defined by
Suzuki&Shigeyama (2009) ApJ
kinetic
energy
of
species
j
in
each
direction
is
defined
by
∞
∞
∞
π
π
mj
vx2 fj dxdydvx dvy dvz ,π π ∞ ∞ (17)
=
∞
mj
2 −π −π −∞ −∞ −∞
2
v
(17)
K
=
jx
x fj dxdydvx dvy dvz ,
π
π
∞
∞
∞
2
mj
−π −π −∞ −∞ −∞
2
エネルギー変化
=
v fj dxdydvx dvmy dvz ,π π ∞ ∞ (18)
∞
j
2 −π −π −∞ −∞ −∞ y K
2
=
v
fj dxdydvx dvy dvz ,
(18)
jy
y
π
π
∞
∞
∞
2
mj
−π −π −∞ −∞ −∞
線形成長を再現 vz2 fj dxdydvx dvmy dvz ,π π ∞ ∞ (19)
=
∞
j
2 −π −π −∞ −∞ −∞
– 18 – x dvy dvz ,
Kjz =
vz2 fj dxdydv
(19)
2 −π −π −∞ −∞ −∞
the magnetic energies are defined by
while
the electric and the magnetic energies are defined by
π
π
1
π! π
(Ex2 + Ey2 + Ez2 )dxdy,
(20)
Ee =
1
2 −π −π
(Ex2 + Ey2 + Ez2 )dxdy,
(20)
Ee =
π
π
2
−π −π
1
2
2
2
π
π
Em =
(Bx + By + Bz )dxdy.
(21)
1
"&"!
2
2 −π −π
Em =
(Bx + By2 + Bz2 )dxdy.
(21)
2 −π −π
d Kjy evolve in the same way, because of the symmetry in the vx -vy
"&"""!
Thetwo
values
of Kthe
Kjyphase
evolveand
in the
the saturation
same way, because of the symmetry in the vx -vy
jx and
ution is divided into
phases,
linear
The time
evolution
is divided
into two phases,
phase, the electric plane.
and magnetic
energies
increase
exponentially.
The the linear phase and the saturation
phase.
In the(Califano
linear phase,
electric
and magnetic energies increase exponentially. The
#%
ent with a linearized
analysis
et al. the
1998).
After
!" the linear
growth
is consistent
with
a linearized
analysisto(Califano et al. 1998). After the linear
cts become significant
and rate
the magnetic
energy
becomes
comparable
)*+
phase,
effectstakes
become
significant
and
the magnetic energy becomes comparable
to
gy of electrons, while
thenon-linear
electric energy
relatively
small
values.
#$
)*,
!"
the total kinetic energy of electrons, while
the electric energy takes relatively small
values.
ween the two phases, Kex increases and Kez decreases. This is one of
-.*/012/
In the transition between the two phases, Kex increases and Kez decreases. This
is one of
3456*02/
es of the Weibel instability, the anisotropy of the motions is mediated
#!"
the remarkable features of the Weibel!"instability,
the anisotropy of the motions is mediated
en particles and electromagnetic fields.
"
'"
!""
!'"
(""
by interactions between particles and electromagnetic
fields.
028*
olution of the ion kinetic energy in each direction and the electric and
From the time evolution of the ion kinetic energy in each direction and the electric and
PIC simulation
Particle-In-Cell : 粒子法の一種
→plasma業界ではスタンダードな手法
最近は2D,3Dの相対論的無衝突衝撃波の
シミュレーションが流行っている
▷Silva et al (2003)
▷NIshikawa et al (2003)
▷Frederiksen et al (2004)
▷Kato (2007)
▷Spitkovsky (2007)
itkovsky
密度
磁場
空間
PIC simulation
Particle-In-Cell : 粒子法の一種
→plasma業界ではスタンダードな手法
最近は2D,3Dの相対論的無衝突衝撃波の
シミュレーションが流行っている
▷Silva et al (2003)
▷NIshikawa et al (2003)
▷Frederiksen et al (2004)
▷Kato (2007)
▷Spitkovsky (2007)
density through a transverse cut at z = 100, with a small inset showing the ion current
00 = 30δi , with the small inset now instead showing the electron current. The arrows
0
h
e
e
3
e
l
2
J. Trier Frederiksen, C. B. Hededal, T. Haugbølle, Å. Nordlund
Density distribution
x方向
F IG . 2.— Electron (top) and ionz方向
(bottom) currents, averaged over the xdirection, at time t = 1200.
Outline
PART I : プラズマ天体物理の基礎
1.Introduction
2.Vlasov-Maxwell system
3.Plasma Instability
4.Wakefield Acceleration
PART II : 最近の興味
1.Bulk comptonization
2.Supernova shock breakout
3.XRF080109/ SN 2008D
Particle Acceleration
統計加速(Stochastic Acceleration) →衝撃波加速(Diﬀusive Shock Acceleration): 粒子がshock frontを往
復してエネルギーを得る -2
→energy spectrum N(ε) ε を再現可能 →但し、効率よくDSAを起こすには始めから高エネルギー粒子が必要
直接加速：電磁場のエネルギーを直接粒子に与える → surﬁng加速
→ リコネクション加速 → wakeﬁeld加速 →ドリフト加速
▷ DSAの seed particle を作れる
Particle Acceleration
統計加速(Stochastic Acceleration) →衝撃波加速(Diﬀusive Shock Acceleration): 粒子がshock frontを往
復してエネルギーを得る -2
→energy spectrum N(ε) ε を再現可能 →但し、効率よくDSAを起こすには始めから高エネルギー粒子が必要
Berezhko&Voelk(2006)
RXJ1713.7-3946
Particle Acceleration
統計加速(Stochastic Acceleration) →衝撃波加速(Diﬀusive Shock Acceleration): 粒子がshock frontを往
復してエネルギーを得る -2
→energy spectrum N(ε) ε を再現可能 →但し、効率よくDSAを起こすには始めから高エネルギー粒子が必要
Berezhko&Voelk(2006)
RXJ1713.7-3946
DENSITY
DENSITY
3HOCK FRONT
3HOCK FRONT
XAXIS
η.
∂ρ
+ ∇ · (ρu) = 0 etc
∂t
XAXIS
η. ^ .
∂f CR
∂f
1
∂fCR
+ u·
= ∇[κ∇] fCR + (∇ · u)p ·
∂t
∂x
3
∂p
$IFFUSION COEFFICIENT
Particle Acceleration
統計加速(Stochastic Acceleration) →衝撃波加速(Diﬀusive Shock Acceleration): 粒子がshock frontを往
復してエネルギーを得る -2
→energy spectrum N(ε) ε を再現可能 →但し、効率よくDSAを起こすには始めから高エネルギー粒子が必要
直接加速：電磁場のエネルギーを直接粒子に与える → surﬁng加速
→ リコネクション加速 → wakeﬁeld加速 →ドリフト加速
▷ DSAの seed particle を作れる
Wakefield Acceleration
元々、レーザー物理業界で提案された加速機構
→Tajima & Dawson (1979)
ultra-intense laserをplasmaへ照射することで電子が加速される →レーザーによる粒子加速器への応用 →Esarey et al. (1996), Mourou et al. (2006)
最近、天体物理への応用が考えられはじめた →Chen et al. (2002) : 宇宙線加速への応用を提案
→Lyubarsky (2006) : relativistic shockで実現されることを示す →Hoshino (2008) : 粒子法によるシミュレーション実験
Wakefield Acceleration
,RQ
(\%] /LJKWSXOVH
(OHFWURQ
light pulseが通過 [D[LV
radiation pressure
(ponderomotive force)
によって電子の空間分布が
変わる
(\%]
/LJKWSXOVH
[D[LV
(\%]
/LJKWSXOVH
([
x方向に静電場が発生
[D[LV
PIC simulation
Particle-In-Cell : 粒子法の一種
→plasma業界ではスタンダードな手法
Hoshino (2008) : 相対論的無衝突衝撃波
(Γ=10)でwakeﬁeld加速が起こることを
シミュレーションで示した
dp
p
=q E+
B
dt
mγ
· E = 4πρ
·B =0
1 B
c t
1 E
B=
+ 4πj
c t
E=
: particle
Simulation
Vlasov eq.
Maxwell eqs.
物理空間1次元,運動量空間2次
元(x,p,q)のVlasov-Maxwell
system
Vlasov eq.とMaxwell eqs.を
交互に時間発展させていく
Source termsを介して粒子と
電磁場が相互作用する
Source terms
First, I divide the
with
of [xmin , xxmax ] ×
minphase
max space
min
maxthe range
min max
p [pmin
q , pmax ] × [qmin , qma
has the volume
∆x∆p∆q.
Thus, ∆x,
∆p,
pmax
− pand
qmax − qmin
xmax
− xmin 5.2.1
discretization
min ∆q are
Thus,
∆x, ∆p, and ∆q
are
, which
∆p = has the volume
, ∆q ∆x∆p∆q.
=
.
(5.14)
∆x = cells each of
N−
N−
pmaxthe
pmin space
qmax
qminof [x p, x ]−×p[p , p ] × [q q , q −] into
xmaxN−x xmin First, I divide
p phase
qrange
x
−
x
with
the
min max
max
min (5.14)
max
minmax
max qmin
min
max
min
, ∆p =
,
∆q
=
.
∆x =
=∆x, ∆p, and ∆q, are∆q =
Np∆x
Nq, ∆p
cells each of which
has =
the volume
∆x∆p∆q.
Thus,
at (x, p, q) = (xi ,N
pjx, qk ), where
Nx
Np
Nq
Semi-Lagrange Scheme
RATEGY
INTEGRATION
xmax − xmin
at (x, p, q)FOR
= (xxiNUMERICAL
,i p=j ,xqmin
where
k ), +
∆x(i − 1/2) for 1 ≤
i≤
,
∆x
= Nx ,
位相空間をmeshに区切る
Its center is located at (x, p, q) = (xi , pj , qNkx),
∆p =
where
33
pmax − pmin
qmax − qmin
,(5.15a)
∆q =
.
Np
Nq
pxji =
(5.15b)
= pxmin
+∆p(j
∆x(i−
−1/2)
1/2) for
for 11≤≤ji≤
≤N
Npx, ,
(5.15a)
min +
there exists a reliable scheme
for the
integration
equations,
the Buneman-Boris
method
[10].
This
Its center
is located
at (x,ofp,these
q)x=
(x
, pmin
), where
ix
j , qk+
=
∆x(i
−
1/2)
for
1
≤
i
≤
N
,
i
x
=
for
11≤≤
kj ≤≤NNqp. ,of motion of a relativistic particle
(5.15c)
= qpmin
+∆q(k
∆p(j −
−
1/2)
for the
(5.15b)
min +
is widelyあるcell(i,j,k)の中にある粒子
used inqpkjPIC
simulations
to1/2)
integrate
equation
equivalent
to
x
=
x
+
∆x(i
−
1/2)
for
1
≤
i
≤
N
,
i
min
x
p≤
1 ≤ j ≤ Np ,
j =
qk =
+ ∆q(k
1/2)
1 ≤t kas,
Nqp. min + ∆p(j − 1/2) for (5.15c)
number
(5.22).of数,エネルギーを定義
particles
of qamin
species
s in−the
cell for
at time
pj = pmin + ∆p(j − 1/2) for 1 ≤ j ≤ Np ,
rst,
using
the
Buneman-Boris
method,
one
obtains
the
orbit
of a+particle
located
the center
q
=
qmin
∆q(k −
1/2) atfor
1 ≤ kof≤each
Nq . cell
xi +∆x/2
+∆q/2
k
j +∆p/2
umber of sparticles of
a species s pin
the cell atqktime
t as,
= qmin +of∆q(k
− 1/2) for 1 ≤
≤ Ncalculated
s qklocation
q.
Iijk
define
the
∆tk are
(t) =the coordinates
dx (Xi , Pj , Q
dpk ) to specify
dqf
(x, p, q,
t), the particle at t +(5.16)
k ) at t.N
pnumber
i +∆x/2
j +∆p/2
k +∆q/2
Next,xixI−∆x/2
define
the
ofqkq−∆q/2
particles
ofs aofspecies
cellatattime
time
−∆p/2
Next,
Ipjdefine
the number
of particles
a speciess sin
in the
the cell
t as,t as,
s
Nijk
(t) =
dx
dp
dqf (x, p, q, t),
(5.16)
xi +∆x/2
+∆p/2
qk +∆q/2
xi +∆x/2
pjpj+∆p/2
qk +∆q/2
articles
contained in
the cell at t;
xi −∆x/2
pj −∆p/2
qk −∆q/2
t+∆t
t+∆t
t+∆t
s
s
s q N (t) =
s q, t),
p
p
dx
dp
dqf
(x, p,
ijk
s
⊥
s
⊥
⊥
N
(t)
=
dx
dp
dqf
(x, p, q, t)
q
+∆q/2
p
+∆p/2
x
+∆x/2
j
i the
kijk B
dt,
P
=
p
+
R
E
+
dt,
Q
=
q
+
R
E
−
B
dt.
(5.23)
articles
contained
in
cell
at
t;
i+
j
j
k
k
x
−∆x/2
p
−∆p/2
q
−∆q/2
q
q
i
j
k
s Γs
qk −∆q/2
Γs dqΓsxfis−∆x/2
Γs
(x, p, q, t). ptj −∆p/2
(5.17)
dp
dx t
=
t Eijk (t)
+∆p/2
jenergy
i +∆x/2and the
pjp−∆p/2
xix−∆x/2
of particles
contained in the cell at t;
k +∆q/2
qkq−∆q/2
s s
s
and
the
energy
of
particles
contained
in
cell
at t;pj +∆p/2
dqΓ
f the
(x,
p,
q, t).
(5.17)
dp
dx
(t)
=
ssume Ethat
the
other
particles
in
qk +∆q/2
xi +∆x/2
ijk
I discretize the electromagnetic
them
only at the positions xi ,
s
qk −∆q/2
pjfields
−∆p/2by defining
xi −∆x/2
dqΓs f s (x, p, q, t).
Eijk
(t) x
=i +∆x/2 dx pj +∆p/2 dp
qk +∆q/2
e cell move with the
same velocp,q
p,q
qk −∆q/2
pj −∆p/2
x⊥i −∆x/2
s s
s
⊥ fields ⊥
⊥
Ie particle
discretize
the
electromagnetic
by
defining
them
only
at
the
positions
x
,
Ei having
(t) = Ebeen
(xi , located
t), Ei at
(t) = E (xE
,
t),
B
(t)
=
B
(x
,
t).
(5.18)
dqΓ
f (x, p, q, t
dp
dx
(t)
=
i
i ijk
i
i
On the other hand, I discretize the
electromagnetic
fields by defining
them 2only at the positio
qk1−∆q/2
pj −∆p/2
xi −∆x/2
er, which is a good approxima⊥
⊥
⊥
⊥
present aEimethod
defined by equations (5.16)(t) = Eto(xcalculate
Eithe
(t)time
= E evolution
(xi , t), of
Bithe
(t)quantities
= B (xi , t).
(5.18)
i , t),
⊥
⊥
sufficiently small
cell.
The
intuEi 1(t)
= 2E (x3i , t), Ei (t)
= Eby
(xidefining
, t), Bi⊥ (t)
= B ⊥only
(x3i , t).
On the other hand, I discretize the
electromagnetic
fields
them
at the
lanation
the scheme
is shown
present
a for
method
to calculate
the time evolution of the quantities defined by equations (5.16)Suzuki&Shigeyama
4
5
In the following, I present a method
to
calculate
the
time
evolution
of
the
quantities
⊥
⊥
⊥
⊥defined
6
5
4
e 5.1.
In each
panel, the hori(2010a)
J.Comp.Phys.
Ei (t) = E (xi , t), Ei (t) = E (xi , t), Bi (t) = B (xi , t)
(5.18).
7
xis represents the x-axis and the
6
7 to8calculate
9
In
the
following,
I
present
a
method
the
time
evolution
of
the
quantities d
axis represents the p- and q-axes.
9
8
h I draw the (5.18).
phase space as two
onal, actual calculations are prein the three dimensional phase
, p, q). The procedure to calcux
x
number of particles in the cell
5.2. STRATEGY FOR NUMERICAL INTEGRATION
there exists a reliable scheme for the integration of these equations
Advection because
part
method is widely used in PIC simulations to integrate the equation of motio
Semi-Lagrange Scheme
5.2.3
FOR NUMERICAL INTEGRATION
33the Buneman-Bori
because
there
exists
a
reliable
scheme
for
the
integration
of
these
equations,
equation
(5.22).
For the integration of the advection part, I make use of the characteristics of the Vlasov equa
first, using to
theintegrate
Buneman-Boris
method,ofone
obtains
orbit of a pp
method is widely used in PICAt
simulations
the equation
motion
of the
a relativistic
RATEGY FOR NUMERICAL INTEGRATION (xdx
q method
pthe33
at t.dp
I defines the coordinates
Pj , QkThis
)sto specify
location
⊥ (Xi ,dq
⊥
⊥
i , pj , qkp) the
equation
(5.22).
a reliable
scheme
for
the
integration
of
these
equations,
Buneman-Boris
[10].
=
,
=
R
E
+
B
,
=
R
E
−
B
,
q
q
cellの中心の粒子の軌跡を計算
s
s
s
as,
dt
Γ
dt
Γ
dt
Γ
At first, the
using
the Buneman-Boris
one obtains
the equivalent
orbit of a particle
located at th
ed in PIC simulations to integrate
equation
of motion ofmethod,
a relativistic
particle
to
(xi , pjfor
, qk )the
at t.
I define theofcoordinates
(Xi , Pj the
, Qk )Buneman-Boris
to specifyt+∆t
the location
of [10].
the particle
する(t
→ t+Δt)
there exists
a reliable
scheme
integration
these equations,
method
This at t
t+∆t
p
q
s
⊥
X
dt,
P
=
p
+
R
E
+
B
dt, toQk = q
as,
i = xi + of motion
j
j
is
widely
used
in
PIC
simulations
to
integrate
the
equation
of
a
relativistic
particle
equivalent
q
s
s
he Buneman-Boris
method,
one
obtains
the
orbit
of
a
particle
located
at
the
center
of
each
cell
Γ
Γ
←Buneman-Boris法
t
t
nefine
(5.22).
the coordinates (Xi , Pj , Qk ) to specify
of the
particle at t +
are calculated
t+∆t
t+∆t
t+∆t
p the location
q ∆t
I
then
assume
that
the
other
particles
in
s
⊥
s
⊥
rst, using the Buneman-Boris
one obtains
the
orbit
of
a
particle
located
at
the
center
of
each
cell
Xi =method,
xi +
dt,
P
=
p
+
R
E
+
B
dt,
Q
=
q
+
R
E
j
j
k
k
q
q
s
s velocthe same cell move
with the same
Γ
Γ
cellの中にN
ijk
(t)個の粒子が一
t
t
t
k ) at t. I define the coordinates (Xi , Pj , Qk ) to specify the location of the particle at t + ∆t are calculated
ity as the particle
having been located at
t+∆t
t+∆t
p
q that
p ⊥
様に分布し、中心の粒子と同
I then
assume
in is a good
the particles
center, which
approximas
⊥ the other
s
⊥
dt, Pj = pj +
Rq E + s B dt, Qk = qk +
Rq E − s B dt. (5.23)
tion
for
sufficiently
small
cell. TheΓintuthe
same
cell
move
with
the
same
velocΓs
Γ
t
t
p,q
p,q
t+∆t
t+∆t
t+∆t
じ軌跡で運動すると近似
p
q
p
itive⊥explanation
the scheme is shown
ity+
as the particle
located
at for
s having been
s
⊥
⊥
dt,
P
=
p
R
E
+
B
dt,
Q
=
q
+
R
E
−
B
dt. (5.23)
i+
j
j
k
k
q
q
s
sin Figure
s
5.1. In each panel, the hori1
the
center,
which
is
a
good
approximathe other
particles
in
Γ
Γ
Γ
t
t
t
zontalThe
axisinturepresents the x-axis and the 2 ijk(t+Δt)とする
tion for sufficiently small cell.
with thecell中の粒子全体が移動した結果、cellに入ってきた粒子数を数えてN
same veloc- p,q
1
3
p,q
vertical
axis
represents
the
pand
q-axes.
ssume
that located
the other
in
itive explanation
for the scheme is shown
aving been
at particles
Although I draw the phase space as two
6
4 2 5
es cell
move
with
the
same
velocin
Figure
5.1.
In
each
panel,
the
hori1
a good approximap,q
p,q are predimensional, actual calculations
E
ijk
(t)も同様に時間
zontal
axis
represents
the
x-axis
and
the
esmall
particle
been located at
cell.having
The intuin the three dimensional phase
3
2 formed
1
3
7
81
9
2
vertical
axis
represents
the
pand
q-axes.
is is
a shown
good approximaspace (x, p, q). The procedure to calcu発展させる
rer,
thewhich
scheme
4
5
Although
I draw the phase
space
as
two
sufficiently
small
cell.
The
intulate
the
number
of
particles
in
the
cell
3
2
3
61
5
each panel, the hori- dimensional, actual4 calculations
pre- at the center of the sur(cell 5) are
located
lanation
for
the
scheme
is
shown
Suzuki&Shigeyama
4
5
7
ts the x-axis and the formed in the three dimensional
6
phase
rounding
cells
at
t
+
∆t
is
as
follows;
(1)
6
5
e 5.1.
In and
each
panel, the hori7
8
94
(2010a)
J.Comp.Phys.
nts
the pq-axes.
Figure 95.1: Schematic views of
calculate
the orbit of a particle located
space (x, p, q). The procedure
to calcu7
8
xis
represents
the
x-axis
and
the
x6
e phase space as two late the number of particles
at the
center of each cell (the left panel)
in
7 the8 cell 9
axis
representsare
theprep- and q-axes.
usingofthe
method. (2)
9
calculations
(cell 5) located at the center
theBuneman-Boris
sur8
he Idimensional
draw the phase
spacerounding
as two cells at t + ∆t is count
the number
phase
as follows;
(1) of particles which enter the original position of cell 5 und
bution of particles inside
each
In other words,
thethe
number
of particles
onal,
actual calculations
are pre- the orbit of a particle
Figure
5.1:cell.
Schematic
views of
integration
of t
located
procedure
to calcu- calculate
particles located
in the gray zones in the right panel of Figure
5.1. Therefor
x
x
inparticles
the three
dimensional
phase
at
the
center
of
each
cell
(the
left
panel)
in the cell
Ap becomes
,he
p, center
q). The
procedure
to
calcuof the sur- using the Buneman-Boris method. (2)
x
x
j+1
number
of
particles
in
the
cell
i+1
k+1the assumption
count
the
number
of
particles
which
enter
the
original
position
of
cell
5
under
+ ∆t is as follows; (1)
|Xi −
s
s
The of
L.H.S
this equation
thethe
advection
the rest andas
kinetic
energy ofbetween
particlesthe
alo
the of
Vlasov
equationrepresents
(5.6), while
R.H.S isofinterpreted
the exchange
of the
Vlasov
equation (5.6), energy.
while the R.H.S is interpreted as the exchange between the kine
and
the where
electromagnetic
and theOn
electromagnetic
energy.
s )2 + p2 + q 2 .
the other hand,
introducing the following variables,
Γs = (Rm
On the other hand, introducing the following variables,
⊥
⊥
⊥
⊥
E
(x)
−
B
(x)
E
(x)
+
B
(x)
s
s
⊥
⊥
⊥
⊥
Here Rm and Rq are dimensionless
constants
defined
by
,
H(x)
=
,
G(x)E=
E (x) − B (x)
(x) + B (x)
2 , H(x) =
2 ,
G(x) =
Transverse成分に関しては、
2
2
⊥ qs
⊥
m
s
one
can
rewrite
the
Maxwell
equations
for
the
perpendicular
components
E
and
B
as
s
s
G(x,t),H(x,t)という関数を定義
⊥
⊥
,
R
=
,
R
=
q
m
one can rewrite the Maxwell equations for the perpendicular components
E and B
as
m
⊥
⊥ e
e
∂G ∂G ⊥ J
∂H
∂H ⊥ J
+
=
−
,
−
.
∂G ∂G
J
∂H
∂H
J= −
∂t =
∂x− Maxwell
+
, 2 equations
−∂t
=∂x−
.to2the followin
すると、source termがある
respectively. On the other
hand,
the
lead
∂t
∂x
2
∂t
∂x
2
In the following,
I integrate theINTEGRATION
above equations instead of the equations for the compon
5.2.InSTRATEGY
FOR
NUMERICAL
線形移流方程式になるので、
the following, I integrate the above equations
instead of∂E
the
⊥ equations
∂E
∂B ⊥ for the components
∂B ⊥
Maxwell eqs.
+
⊥
s
∂t
∂t
∂x
∂t
where the
dimensionless
electric for
current
densities J andintegration
J are expressed in terms of f (x, p, q, t) a
5.2
Strategy
numerical
これを解けばよい
5.2
= −J ,
+
Strategy for numerical integration
∞
∞
p
= −J ⊥ ,
s
=
Rqs numerical integration
f
(x, p, q, t)dpdq,
In this section, I describeJa method
for
of the dimensionless Vla
s
Γ
−∞
−∞
In this section, I describe a method fors numerical integration of the dimensionless Vlasov-M
(5.11) introduced in the previous section. ∞ ∞
(5.11) introduced in the previous
section. s
q s
⊥
J
=
Rq
f (x, p, q, t)dpdq.
s
−∞ −∞ Γ
s
5.2.1 discretization
5.2.1
discretization
These two relations close the system.
First, I divide the phase space with the range of [xmin , xmax ] × [pmin , pmax ] × [qmin , qmax
First, I divide the phase space with the range of [xmin , xmax ] × [pmin , pmax ] × [qmin , qmax ] into
cells each of which has the volume ∆x∆p∆q. Thus, ∆x, ∆p, and ∆q are
cells
of which has the of
volume
∆x∆p∆q. Thus, ∆x, ∆p, and ∆q are
5.1.3 each
Transformation
equations
Suzuki&Shigeyama
pmax − pmin
qmax − qmin
xmax − xmin
p
−
p
q
x
−
x
,
∆p
=
,
∆q
=
.
∆x
=
max
min
max − qmin
max I transform
min
(2010a) J.Comp.Phys.
For convenience of the following
the
equations
introduced
above.
,
∆p
=
,
∆q
=
.
∆xsections,
=
Nx
Np
Nq
s
N
N
N
x some algebraic manipulations
p
q
Multiplying equation (5.6) by Γ (p, q) and
lead to the following
eq
Its center is located at (x, p, q) = (xi , pj , qk ), where
s s
s s
s s
Its
center
), where
∂(Γ f ) isplocated
f )
p
∂(Γ
f )
∂(Γs f s )at (x,sp, q) = (xq i , p⊥j , qk∂(Γ
s
s pE + qE
⊥
⊥
+ s
+ Rq E + s B xi = xmin +
+ ∆x(i
Rq E− 1/2)
− s Bfor 1 ≤ i ≤ =
NR
,q
∂t
Γ
∂x
Γxi = xmin∂p
Γ
∂q
Γs
+ ∆x(i − 1/2) for 1 ≤ i ≤ Nx , x
pj = pof
+rest
∆p(j
−kinetic
1/2) energy
for 1 of
≤ particles
j ≤ Np , along the
min
The L.H.S of this equation represents thepadvection
the
and
=
p
+
∆p(j
−
1/2)
for
1
≤
j
≤
N
j
min
p,
q
+ ∆q(k
− 1/2)
for between
1 ≤ k ≤the
N kinetic
.
of the Vlasov equation (5.6), while the R.H.Sqkis =
interpreted
as the
exchange
ener
qk = qmin +min
∆q(k − 1/2) for 1 ≤ k ≤ Nq . q
and the electromagnetic energy.
Next, I define the number of particles of a species s in the cell at time t as,
Next,
definehand,
the number
of particles
of a species
s in the cell at time t as,
On theI other
introducing
the following
variables,
Simulation Setup
初期条件
▷plasmaは静止している ▷thermal dispersionも無い
▷電磁場も存在しない
Initial condition
Boundary condition
境界条件 ▷電磁パルスを入射させる
(\%]
/LJKWSXOVH
,RQ
(OHFWURQ
[D[LV
Energy spectrum
位相空間(x,p,q)の運動量の角
度依存性と空間依存性を積分
Lorentz factorに対する粒子
数
Power-law成分が現れる
(index -2)
前半のまとめ
無衝突プラズマの基礎
Vlasov-Maxwell方程式系: 位相空間中でのプラズマのダイナミクス
但し、定性的な振る舞いは荷電粒子の運動で理解できる
two-stream instability, Weibel instability, wakeﬁeld加速
非熱的放射の起源
粒子加速
衝撃波のbulk kinetic energy
→高エネルギー粒子 (Fermi加速など)
→非熱的放射 (synchrotron, 逆コンプトン)
ここまでの話
(無衝突プラズマ)
光子ープラズマ相互作用
衝撃波のbulk kinetic energy
→非熱的放射 (bulk comptonization)
ここからの話
(衝突プラズマ,流体)
Outline
PART I : プラズマ天体物理の基礎
1.Introduction
2.Vlasov-Maxwell system
3.Plasma Instability
4.Wakefield Acceleration
PART II : 最近の興味
1.XRF080109/ SN 2008D
2.Supernova shock breakout
3.Bulk comptonization
XRF 080109/SN 2008D
渦巻き銀河NGC 2770に出現
したIb/c型超新星
Swiftが超新星の誕生の瞬間に
付随していると思われるX線
放射を観測
X-ray spectrumは
single power-law/powerlaw+black bodyでﬁtできる
→photon index Γ -2.3
→Tph 0.1keV
Soderberg+(2008)
XRF 080109/SN 2008D
VLBIでの電波観測で、ejecta
を空間分解
→ ejecta velocity < 0.75c
250
4
200
MilliARC SEC
2
20
98bw
(broad Ic)
04gq (Ib)
15
02ap
(broad Ic)
10
08D (Ib)
5
100
05kl (Ic)
50
-4
0
5 10 15
Bietenholtz et al.(2008)
-6
6
4
2
0
MilliARC SEC
-2
-4
-6
-1
Doppler velocity (10 km s
0
60o
o
90o
3
-2
10
30o
5
03jd
(broad Ic)
0
-15 -10 -5
15
70
150
0
[O I] 6300
Model BP8
0o
Normalized flux + const
6
20
Normalized flux + const
nebular phaseの可視光観測
でdouble-peaked [OI] を発
見 → side-viewed bipolar
ﬂow
[O I] 6300
0
-15 -10 -5
0
5 10 15
Doppler velocity (103 km s-1)
Tanaka et al.(2009)
Outline
PART I : プラズマ天体物理の基礎
1.Introduction
2.Vlasov-Maxwell system
3.Plasma Instability
4.Wakefield Acceleration
PART II : 最近の興味
1.XRF080109/ SN 2008D
2.Supernova shock breakout
3.Bulk comptonization
重力崩壊型超新星
*VYLJVSSHWZL
Fe
massive stars(M>8-10M◉ )の
進化の最終段階
Fe coreの光分解で、星の内部
に衝撃波が発生
*VYLIV\UJL
⁵⁶Fe+γ→13⁴He+4n
⁴He+γ→2p+2n
p+e-→ν+n
e
:OVJRWYVWHNH[PVU
:OVJRIYLHRV\[
衝撃波が星の外層を突き抜ける
ことで、ejectaを放出
但し、重力崩壊から衝撃波の発
生と伝播まではよく分かってい
ない
ρc =3 10¹⁴[g/c.c.]
:5L_WSVZPVU
UV/X-ray ﬂash
post shock ∼0.1keV
Optical
Shock Breakout
超新星爆発において、星のcore
付近で発生した衝撃波がその星
の外層をつき抜ける現象
重力崩壊型超新星では、星の外
層を伝搬する衝撃波は
radiation-mediated shock →衝撃波の前方のoptical
depthがτ c/Vsとなるときに
shockからphotonがもれだす
つまり、shock breakout以降
に超新星は観測可能になる
GHQVLW\
9V
τ
[
HOHFWURQ
SKRWRQ
photon diﬀusion velocity Vdiﬀ=c/τ
radiation temperature Tr∼106K
SNLS-04D2DC
Supernova Legacy Survey
GALEX衛星によってUV
ﬂashが検出され、その後超新
星を観測
Schawinski+(2008)
Outline
PART I : プラズマ天体物理の基礎
1.Introduction
2.Vlasov-Maxwell system
3.Plasma Instability
4.Wakefield Acceleration
PART II : 最近の興味
1.XRF080109/ SN 2008D
2.Supernova shock breakout
3.Bulk comptonization
Comptonization
photonが電子に散乱されるこ
とでspectrumが歪む
thermal compton=SZ効果
Kν
bulk compton
(Blandford&Payne 1981)
→電子の、圧縮のある流れがあ
るときに起こる
comptonization(e.g. 衝撃波,
降着流)
Kν
HOHFWURQ
SKRWRQ
速度大
速度小
Kν
Kν
HOHFWURQ
SKRWRQ
Comptonization
GHQVLW\
Fermi加速とのアナロジー
shock front付近での散乱:
→ エネルギーA倍 (A>1)
shock frontからの脱出確率:
→ 1-P
Pはoptical depthから決まる
(OHFWURQ
/RJ1(
[
1
31
31
31
31
31
($($($($($(
(QHUJ\(
Comptonization
GHQVLW\
Fermi加速とのアナロジー
shock front付近での散乱:
→ エネルギーA倍 (A>1)
shock frontからの脱出確率:
→ 1-P
(OHFWURQ
/RJ1(
Pはoptical depthから決まる
[
1
31
31
31
31
31
本研究の狙い
Shock breakoutでのbulk comptonization
でSN 2008Dのnon-thermal emissionを説
明するにはどんなパラメータが必要か？
($($($($($(
(QHUJ\(
Shock proﬁle
GHQVLW\
ejecta中のphotonの伝播を計
算するには、scattering
source(電子)の速度,密度が必要
→shock proﬁle
8 N[
L
ρ=k x A
plane-parallel
Sakurai (1960) : 星の外層を
shockが通過するときの自己相
似解(α=3: radiative
envelope)
atmosphere
[
!&
)
!"
'%"&
&
)
!'%"&
012345-6%7489:;
/0123,4'567+7+8
,-.(+!
,-.%+&
,-.$+!
,-.#+$
!%
!"
!$
!"
!#
!"
!"
%'!"
!"
&'!"
!"
('!"
!"
)'!"
!"
*'!"
!!
!'!"
/32,9170':;<=',>0'2?;:970'57=8
!+!'!"
!!
)
!$%"&
)
!*%"&
)
!(%"&
"&
!"%"&
!"#'%"&
"&
!"#$%"&
"&
-.!,#"
-.!+#*
-.!$#"
-.!/#$
"&
+%"&
"&
*%"&
"&
,%"&
"&
(%"&
"&
)%"&
""
"%"&
<5:-=>41%?@38%-A1%:B@?=41%748;
"#"%"&
""
Photon propagation
monte-carlo法でphotonの伝
播を解く
各タイムステップごとにshock
front からseed photonを入射
させる
density proﬁleからoptical
depthを計算
←Thomson断面積を仮定
GHQVLW\
cΔt
[
scattering probability: 1-exp(-neσcΔt)
散乱される場合は、散乱する電
子の速度をvelocity proﬁleか
ら決めて、photonの散乱後の
エネルギーと散乱角を計算する
スペクトル
GHQVLW\
LQLWLDO
shock proﬁleのパラメータ:
Vi=0.3c,
11
Xi=10 [cm]
ILQDO
9L
9I
shock breakoutが起こる瞬
間から計算を開始する
→τ=c/Vi → 1
[L
!&
)
!"
'%"&
&
)
!'%"&
012345-6%7489:;
/0123,4'567+7+8
,-.(+!
,-.%+&
,-.$+!
,-.#+$
!%
!"
!$
!"
!#
!"
!"
%'!"
[
!"
&'!"
!"
('!"
!"
)'!"
!"
*'!"
!!
!'!"
/32,9170':;<=',>0'2?;:970'57=8
!+!'!"
!!
)
!$%"&
)
!*%"&
)
!(%"&
"&
!"%"&
!"#'%"&
"&
!"#$%"&
"&
-.!,#"
-.!+#*
-.!$#"
-.!/#$
"&
+%"&
"&
*%"&
"&
,%"&
"&
(%"&
"&
)%"&
""
"%"&
<5:-=>41%?@38%-A1%:B@?=41%748;
"#"%"&
""
スペクトル
GHQVLW\
seed photonは0.1keVの
black bodyで衝撃波面から入
射させる
LQLWLDO
ILQDO
9L
9I
高エネルギー側にpower-law
tail: Γ -2.4
7(8+&5/09/3:0*0;'/<,"("=
#!!!
[L
%&'()*
+),-./+012
304&56),4
#!!
#!
#
!"#
!"!#
#!
#!!
#!!!
3:0*0;/&;&5>2/<&?=
$
#!
[
衝撃波速度への依存性
衝撃波速度を変えて計算
赤: Vi=0.1c
青: Vi=0.2c
緑: Vi=0.3c
黒: Vi=0.4c
non-thermal tailが顕著になる
には、衝撃波速度 Vs>0.3cが必
要
但し、Vs>0.3cではindexはほ
とんど変わらずΓ -2
後半のまとめ
XRF 080109/SN 2008DのX線
観測で得られたphoton indexを
再現(Vs>0.3c)
7(8+&5/09/3:0*0;'/<,"("=
SN shock breakoutにbulk
comptonizationを考慮すること
でnon-thermal tailが作れる
#!!!
#!!
#!
#
!"#
!"!#
#!
13
R=10 [cm]を仮定することで、
energeticsとduration( 500s)を
説明できる
%&'()*
+),-./+012
304&56),4
#!!
#!!!
3:0*0;/&;&5>2/<&?=
$
#!
非熱的放射の起源
粒子加速
衝撃波のbulk kinetic energy
Suzuki&Shigeyama (2008) Phys.Plasmas
Suzuki (2008) Phys.Plasmas
Suzuki&Shigeyama (2009) ApJ
Suzuki&Shigeyama (2010a) J.Comp.Phys.
→高エネルギー粒子 (Fermi加速など)
→非熱的放射 (synchrotron, 逆コンプトン)
光子ープラズマ相互作用
前半の話
(無衝突プラズマ)
後半の話
衝撃波のbulk kinetic energy
(衝突プラズマ,流体)
→非熱的放射 (bulk comptonization)
Suzuki&Shigeyama (2010b) ApJ submitted

Download Report