モーメント法によるニュートリノ輸送の数値解法
少ない計算量で影、散乱、伝播を扱える近似解法

花輪 知幸、菅野 裕次 (千葉大学)

モーメント法: 拡散近似の改良
Intensity
輻射強0th:1st:Iν (r, t, n)Iν (r, t, n)!

1Eν (x, t) ≡I (x, t, n) dnc!
F ν (x, t) ≡n I(x, t, n) dn!
←→
1
P ν (x, t) ≡nn I(x, t, n) dnc!

1Eν (x, t) ≡I (x, t. n dn) dn
2ndc!
拡散近似
F" ν (x, t)E #≡ M1 modeln I(x,(E,t, F)n) dn
1∂+ n · ∇ I (r, t,Pは推定n) = − σν Iν (r, t, n)
Fは推定!
c ∂t
←
→1!

Moment!2

"←→1P ν (x, t) !≡c
1nn I(x, t, n) dn

輻射輸送方程式
I (x, t. n dn) dnc#

Eν (x, t) ≡(2)

輸送吸収1"∂#F ν (x, t) ≡n I(x, t, n) dn
(3)
+n·
∇
I
(r,t,n)=−σI(r,t,n)νν1∂c ∂t+ n←·→∇ I (r, 1t,! n) = − σν Iν (r, t, n)
P ν (x, t) ≡nn I(x, t, n) dn !
(4)c ∂tc!
!
+ +σνσSSν (r,t,n)+σg(n,n)I(r,!
! t, !nν,sν
(r, t, n) + σ
g (n, n ) I (r, t, n ) dn!

ν νν,sν"#厳密
1∂+ n · ∇ I (r, t, n) = − σν Iν (r, t, n)c ∂t!
+ σν Sν (r, t, n) + σν,s g (n, n! ) Iν (r, t, n! ) dn! (5)

放射散乱

∂Eν
+ ∇ · F ν = σν (4πSν − cEν )∂t∂Eν+ ∇ · F ν =←σν (4πSν − cEν )→∂F∂t ν+ ←→∇ · P ν = − c (σν + σν,s ) F ν∂F ν+ ∇ · P ν = − c (σν + σν,s ) F ν∂t∂t

M1Closure relation (終結条件)←→Pν ="→1 − χ←3χ − 1I +nn
2
2#Fν
3 + 4f 2|F ν |$Eν , n =, χ =, f =Fν
Eν
5 + 2 4 − 3f 2

∆tnnn∆xF−FU n+1−Ui+1/2,j,ki−1/2,j,ki,j,kni,j,k+= S i,j,k ∆t∆xnn+1n− U i,j,k= USi,j,ki,j,k

風上差分 (特性曲線)
F−Fni+1/2,j,kni−1/2,j,k+∂Eν∆xcEν )+ ∇∆t· F ν = σν (4πSν −(simple) HLLn∂t=Si,j,k,j,k − λL F i+1,j,k + λR λL (U i+1,j,k − U i,j,k )λR F i,j,k −∂F
λL Fλ→U i,j,k )i+1,j,k + λR←
L (U i+1,j,k −(8)(HLL)ν
F i+1/2,j,k = λR − λL+λR∇− ·λLP ν = − c (σν + (8)σν,s ) F ν∂t
λR = c.
λL = −c(9)λR = c, λL = −c(9)安全だが拡散大
信号伝播の最大値,λR F i,j,k − λL F i+1,j,k + λR λL (U i+1,j,k − U i,(HLL)
F i+1/2,j,k =Godunov(特性曲線と伝播モード)"#λR − λL←→→1 − χ←3χ − 1Fν
3(定義も曖昧）P ν 拡散を抑えられるが、計算が面倒=I +nn E, χ =ν, n =
λR = c, λL = −c
22
Fν
5+

再構成 (本研究)Eν Fx,ν ,U =  Fy,ν Fz,ν



∂



2 3



3E
ν (1 − β )−4 ν −Fσ(4πSx,ννIν (n) = 2(1−β·n)2
 c P8π −c (σν + σν,3+βxx,ν
Fx = ,S=3f
 c2 Pxy,ν  −c (σνF + σν,!
β =, β = β|+ 4 − 3f 2c2 P2xz,ν
−c (σ|Fν + σν,Closure∂ relation∂ を満たす輻射
∂!4U +Fx +Fy +Fz = S

M1 モデルは影が扱える
入射吸収体

Diffusion
M1 (HLL)
M1(characteristics)

Gonzalez et al. ‘06

 , ∂F
+∇·F−cEν == σν (4πSν+ν)4 −22F= F=,Sν52xν←→2∂t Fy,ν cP−c(σ+σ)+ ∇· P ν =−c(σ+σ)
Fxy,ννν,sνν,sν←→∂F ν∂t2+ ∇ · P ν = − −cc (σν(σ+ σν,s)Fν )Fz,νcP+σxz,ννν,s∂tM1モデルの数値積分
"#EFσ(4πS−"#νx,ννν
2cEν )
保存形←→→ ←1 −∂χ←3χ ∂−11 − χ ←
F1ν ∂
3 +F
4fν
→→
3χ
−∂
2$σ , )F FP=I +nn EIν , +n =,nnχ−c=,= fν P=E,n=χcP(σ+ννx,νxx,ννν,sx,ν
2F=SF
 , 2F x U= + 22 F x2 + , FSy 2+z
ν
5+2
4−3f
F= ν
∂t∂y∂z Fy,ν  ∂x
c Pxy,ν−c (σν + σν,s )Fy,ν


2
 Fx,ν−c (σν + σν,s)FFz,ν cPExz,ν
σν z,ν(4πSνEF
σ (4πS − cE )x,ννν 2 ν Fx,νcP−c (σν +xx,ν 2 Fx,ν 
c Pxx,ν
U = ,  F x =   2−c (σν+, σν,sS)F=x,νU =  ∂  , F∂x =  Fy,νS =c Pxy,ν−c (σν +2 ∂ ,∂Fy,ν
c Pxy,ν
2 −c (σν + σν,s )Fy,ν
Fz,ν
cP−c (σν +U+F+F+F=Sxz,ν
2xyzFz,νc P∂yxz,ν∂t∂x∂z −c (σν + σν,s )Fz,ν
∂∂∂∂U +Fx +Fy +Fz = S
1∂t∂x
∂y∂z
ν

F x,i+1/2,j,k = F x,i+1/2,j,k (U x,i,j,k , U x,i+1,j,k )F x,i+1/2,j,k = F x,i+1/2,j,k (U x,i,j,k , U x,i+1,j,k )(i, j, k)(i+1, j, k)

Numerical Fluxnnn
F−FU n+1−U
i+1/2,j,ki−1/2,j,ki,j,ki,j,k+
∆t
∆x1n
=
Si,j,k1x

1!6

Pν,xz,i+1/2,j,k =「風上」性を考慮した再構成∗
Iν,i+1/2,j,k (n)=(+)
Pν,xz,i+1/2,j,k#+(−)
Pν,xz,i+1/2,j,kIν,i,j,k (n)(nx > 0)Iν,i+1,j,k (n) (nx < 0)nnnF−FU n+1−Ui+1/2,j,ki−1/2,j,ki,j,ki,j,k+∆t∆x= S ni,j,kλR F i,j,k − λL F i+1,j,k + λR λL (U i+1,j,k − U i,j,k )k =λR − λL
λR = c,λL = −c(8)(9)

23Eν (1 − β 2 )3−4Iν (n) =(1−β·n)8π 3 + β 23fF!
β =, β = β|F |2 + 4 − 3f 2

解析関数

(10)

M1 closure と整合する輻射場

セル内での吸収・放射・散乱境界のF!(+)
Fν,i+1/2,j,k減衰率
=中心のF

放射!

"−(∆τa,i +∆τs,i )/2 (+)
eFν,i+1/2,j,k + 1 − e−∆τa,i,j,k /2 πBν,i,j,k!(+)−(∆τa,i +∆τs,i )/2 !(+)!
Fν,i+1/2,j,k = e " cEν,i,j,k Fν,i+1/2,j,k +−∆τs,i,j,k /2−(∆τa,i,j,k )/2+e1−e!

1−"! 散乱4−(∆τa,i,j,k )/2! 1 − e−∆τs,i,j,k"/2!(+) !(+) !(+)!(−) !(−)+eFν,x,i+1/2,j,k+F−(∆τa,i +∆τs,i )/2 (+)= F=ν,x,i+1/2,j,k+F
k1/2,j,kFν,i+1/2,j,kν,x,i+1/2,j,k= eν,x,i+1/2,j,k Fν,i+1/2,j,k + 1 − e−∆τa,i,j,k /2 πBν

" cE! "!
"!
Sν,i,j,kS−∆τ/2ν−(∆τ
)/2 ,/2(+) 光学的に厚いν−∆τ(+)−∆τ/2−∆τi /2−∆τ(∆τ∆τ
"1)
ところで良い近似i /2i−s e+1e+ai 1=eF
1−e=eF+−e!(+)(+)
4
k1/2,j,kν,x,i+1/2,j,kPν,i+1/2,j,k=ν,x,i+1/2,j,ke−(∆τa,i +∆τs,i )/2 Pν,i+1/2,j,k 4
4!
" 2πB " S!
"!
Sνν
ν,i,j,k(+)
(+)
−∆τ/2
−∆τi+1/2−∆τ/2−∆τi+1 /2 + 1 − ea,i,j,ki+1 −∆τi+1 /2= e= Fe FF’− 1−− e1 − eFν,x,i+1/2,j,kk1/2,j,k
ν,x,i+1/2,j,k!(+)−(∆τ +∆τ 3c)/2 (+)
44
Pν,i+1/2,j,k = e!
Pν,i+1/2,j,k" E!(+) !(+)!(−) a,i,j,k−(∆τ)/2−∆τs,i,j,k /2 " ν,i,j,k!(−)!
1 − e −∆τ= P=ν,xx,i+1/2,j,k++Peν,xx,i+1/2,j,kPν,xx,i+1/2,j,k+ Pν,xx,i+1/2,j,k
ν,i,j,kk1/2,j,k/2 2πB
6+ 1−e3c!
"S!
"!
" Eν
(+)
S−∆τ
/2−∆τ/2ν(+)
ii−∆τ/2−∆τ/2
ν,i,j,k
光学的に厚い(∆τ,i∆τs " 1) ところで良い近似
−(∆τ
)/2 i
−∆τ
/2=e
P
+
1−ea
Pν,x,i+1/2,j,k + +1e − e 6 1 − ek1/2,j,k = e
ν,x,i+1/2,j,k66!
"2次精度は不安定
S!
ν " Sν
(−)−∆τ/2−∆τ/2i+1i(−)
光学的に厚い(∆τ,∆τ
"1)
ところで良い近似−∆τP−∆τi /2a+s1 − ei+1 /2=ePν,x,i+1/2,j,k + 1 − e 6
k1/2,j,k = e
ν,x,i+1/2,j,k6a,i,j,ks,i,j,ka,is,ia,i,j,ka,i,j,ks,i,j,k時間進化
1次精度
移流は単純な 陽解法
放射・吸収・散乱は「局所平衡」へ指数関数的に近づける (陽解法)
E (t+Δt) = η E(t) + (1 - η) Eeq
η= exp ( - Δt / τ)影
simple HLL

this work!10

ビーム
テスト
(1)
光の伝播
1次精度
Δx = 0.1!11∂H= J −S = 0J∂τ= 3H [τ + q (τ )]∂Kq閃光テスト" =2/3 H(1)∂τ
E =Econst.
flim
F
F → F=0|F |J = 3H [τ + q (τ )]> flim Eq "|F |2/3

(14)(15) (13)
(16)

(14)
(17) (15)ならflim E F
F →|F |

(16)

極端に細いビーム3は避ける!12閃光(2)
等方性とビームの鋭角性は二律背反!13

閃光 (3)E = const.
F = (0.9, 0, 0)
矢印: f = F/
HLLでは「拡散」「偏心」!14

t, n) dn(3)

SEDも計算できる(4)(x, t, n) dn

ex. 照射された原始惑星円盤からの赤外線Iν (r, t, n)
赤外線!
g (n, n! ) Iν (r, t, n! ) dn! (5)可視(紫外)4πSν − cEν )(σν + σν,s ) F ν

cf. Calvet et al. ’91
butFν
3 + 4f 2|F ν |$n =, χ =, f =FνEν
5 + 2 4 − 3f 2!15

高温層の場所は入射角度に依る!16

まとめ• M1 model は 影, 散乱, 伝播が扱える• 相対論的な天体では陽解法で計算可能• 再構成法は「良い」数値流束を与える• 安定: E > |F| > 0 が保たれる (Δt<Δx/c)• 光学的に厚いところで拡散近似に収束• 「影」がシャープ!17