Radiative Transfer Calculation bySolving Moment Equations

Tomoyuki Hanawa, Yuji Kanno (Chiba U.)

Radiative Transfer is a Key Issue inMany Astrophysical Situations• Neutrino Emission in Core Collapse SNe• Neutrino Emission from Merging Binary NSs• Super Eddington Luminous BHs• Irradiated Protoplanetary Disks– and more

• M1 scheme achieves an intermediate angularresolution at low cost.
– Diffusion approximation is coarse while fullradiative transfer costs much.
!2Contents– M1 Moment Equations• Comparison with diffusion
approximation and full radiativetransfer• Reconstruction Method–with & without absorption– 1D Plane Parallel Equilibrium– 2D and 3D Examples• Shadow, Directional Characteristics, Propagation– Irradiated P.-P. Disk (２color）!3

Description of Radiation Field
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
2nd:c!
Diffusionmodelt, (E,F" ν (x,Approx.
t) #≡ E M1n I(x,n)F)dn1∂
estimate
Festimate
P+ n · ∇ I (r, t, !
n) = − σν Iν (r, t, n)
c ∂t4
←
→1!

Moment!

M1 Model can handle “shadow”.
Radiation

Absorber

Diffusion
M1 (HLL)
M1(characteristics)

Gonzalez et al. ‘06

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

EquationTransfer
E (x, t) ≡ of IRadiative(x, t. n dn) dn(2)c#νTransfer
Absorption
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 ∂t
c!
!
+ +σνσ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∂
Exact+ 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)"Emission∂Eν+ ∇ · F ν = σν (4πSν − cEν )∂t∂Eν+ ∇ · F ν =←σν (4πSν − cEν )→∂F∂t ν+ ←→∇ · P ν = − c (σν + σν,s ) F ν∂F ν+ ∇ · P ν = − c (σν + σν,s ) F ν∂t∂t

M1

closure relation←→
Pν =Scattering"→
1 − χ←
3χ − 1I +nn
22#Fν
3 + 4f 2|F ν |$Eν , n =, χ =, f =FνEν5 + 2 4 − 3f 2η = tan←→∂F ν+ ∇ · P ν = − c (σν + σν,s ) F ν∂tFx

Reconstructed Radiation Field←→Pν ="

2 3#3(1−β)
E2→
1 − χ←
3χ − 1I (θ, ϕ) =Fν3+4f|Fθνsin| ψ[1−βcosθ
cos
ψ
−βsin
2$I +nn Eν , n8π=(3 +, βχ ) =, f =222Fν
Eν
5 + 2 4 − 3f

3 (1 − β 2 )3 E−4I (θ) =(1−βcosθ)
8π (3 + β 2 )closure relationThis intensity defined by Eq. (14) reproduces the closure relatio
F.

Fz(+) ≡( π/2 ( 2πI (θ, ϕ) cos θ sin θ dθ dϕ( π/2 ( 2πI (θ, ϕ) cos2 θ sin θ dθ dϕ( π/2 ( 2πI (θ, ϕ) cos θ sin!72 θ cos ϕ dθ dϕ00Fz(+)
3 q + 6 β 2 cos2 ψ q −1 − β 4 cos2 ψq −3 + 8=E
4 (β 2 + 3)
1+Pzz
≡Pzz
E0303β cos ψq −1 + 3β cos ψq + 4β 2 c2 + 1 −=2 (β 2 + 3)+Pzx≡00+ ν∇−·FFcE− cE+ ∇ · F ν = σν (4πSν ν )= σν (4πSν+ν)4 −2ν52∂t∂t  c2 P Fx,ν −c(σ+σ)Fxx,ννν,sx,ν , F x ∂F

←
→
∂F
←→ν
ν
=
S=
2· P ν =−, c (σ+∇·σPc (σν++ σσν,s ) )FFν+ ∇+
 Fy,ν νν,s )νF=
ν−c−(σcPxy,ννν,sy,ν
∂t∂t
IntegrationM1 Model
Numericalof2FEz,ννc F"Px,ν−c(σν ν+ −σν,sxz,νcE)F#σν (4πSν )z,νConservationForm←→ 2 1 − #←→χ
3χ − 1 Fν" F←2, )FPI+nnE,n=χx,ν=cP−c(σ+σν3χ=
νx,ν
xx,ν
νν,s→
←→
1−
χ−
1F
3+4fν Pν 2 ∂ Eν,, nS=2∂Fν
,F==$=I+nn,χ=, fx2∂∂
 Fy,ν 2P(σ+σ)F
2 +F νF −cxy,ννν,sy,ν
5+24−3fUFc2+F+=Sxyz2∂t∂x c PExz,ν
∂y∂z Fx,ν−c (σν + σν,s )FFz,ν
σν z,ν(4πSν2 c2 Pxx,ν −c (σν + Fx,ν  
 , F x =  2 σν (4πS ,ν −ScE=ν ) U = EνFx,ν
 Fy,ν
cP−c (σν +xy,ν
2 ∂  F∂∂∂c Pxx,ν x,ν 
2−c (σν + σν,s )Fx,ν Fz,ν
cP−c (σν +U = , Fx F= x+, +S = 
U+FF=Sxz,ν
2yz
Fy,νc P∂y−c (σν + σν,s )Fy,ν
xy,ν
∂t
∂x
∂z
∂
∂ z,νFz,ν
c2 Pxz,ν ∂−c (σ∂ν + σν,s )F
U +Fx +Fy +Fz = S
∂t∂x
∂y∂z





F x,i+1/2,j,k = F x,i+1/2,j,k
1 (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 FluxnnnF−FU n+1−
Ui+1/2,j,ki−1/2,j,ki,j,ki,j,k+∆t∆x1= S ni,j,k1x1!8∆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

Upwind (Characteristics)

F ni+1/2,j,k − F ni−1/2,j,k+∂Eν
∆t· F ν = σν (4πSν −∆xcEν )+∇(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 )(HLL)
F i+1/2,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)ν= λR − λL+λR∇− ·λLP ν = − c (σν + (8)σν,s ) F ν

λR = c,∂tλR = c.
λL Signal= −c(9) Diffusive
λL = −c(9) Safe but
Max.
Speed,λ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 =(characteristicsand eigen modes)
λR − λLGodunov"#←→ diffusive→1 − butχ←3χ − 1(mean characteristicsF ν are
3Less
costs more.
Pν =I +nn Eν , n =, χ =λR = c, λL = −cnot well defined.）
22Fν
5+Reconstruction (this work)

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 2
c2 P2xz,ν
−c (σ|Fν + σν,

consistentwith
theclosurerelation∂∂∂∂!9
U +Fx +Fy +Fz = S

Pν,xz,i+1/2,j,k =UpwindReconstruction∗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)nnn
F−FU n+1−U
i+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

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

2Consistentwith M1(10)closure3f!
β =,2 + 4 − 3f 2Fβ = β|F |(+)
Fν,x,i+1/2,j,k(−)
Fν,x,i+1/2,j,kReconstructed Numerical Flux
Fν,x,i+1/2,j,k =(+)

Fν,x,i+1/2j,k =(−)
Fν,x,i+1/2j,k=Pν,xx,i+1/2,j,k =Pν,xy,i+1/2,j,k =Pν,xz,i+1/2,j,k =∗Iν,i+1/2,j,k(n)=""n x > 0∗nx Ii+1/2,j,k(n) dn∗
nx Ii+1/2,j,k(n) dn

nx < 0(+)
Pν,xx,i+1/2,j,k(+)
Pν,xy,i+1/2,j,k(+)
Pν,xz,i+1/2,j,k#+(−)+ Pν,xx,i+1/2,j,k(−)+ Pν,xy,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)without
absorption(11)

Numerical fluxes are explicit functionsof E and F.
11!

Numerical Fluxes Given by Reconstruction+Fz=!

3q + 6β c q − c β + 8βcq
42 2 2%4 4

3"#8 (3 + β ) q
2&

3$−1

3 c3
1β+2 22Pzz =+ 3βcq + 4β c + 1 − β
22 (3 + β )q2
3 + 4f!
χ=(12)
fx+
224−3fPxz= β 5F+zffy+= βFz Consistent with the M1 clopecificPyzIntensity
furefzc = cos ψ ≡f''2 )33E
(1−β222 + β 2−42q=
1−βs=
1−βc
I (θ, ϕ) =(1 − β cos χ)(13)
28π 3 + β3f!
β =(14)2 + 4 − 3f 2!12
f cos χ = fz cos θ + sin θ (fx cos φ + fy sin φ)(15)(26(27(28(29(30(31+),j,k=(+)
Fz,i,j,kπ/2=2πI (θ, ϕ) cos θ sin θ dθ dϕFluxes
EachCubeFace
I (θ, ϕ)fromcos θ sinθdθdϕ(26)I (θ, ϕ) cos θ sin θ dθ dϕ
F=! π/2 ! 2π00! 2π0 ! π/2 0(+)
z,i,j,k00(27)! π/2 ! 2π(+)(+)(−)+−Fz,i,j,k(+) = (−)(+)I (θ,(−)ϕ) cos θ sin θ dθ dϕ(+)
xz
FFFFPPPP
Ix(θ,ϕ)cossindθxxdϕ(26)
Fz,i,j,k−(+)(−)(+)z (−)(+)
(−)
+ θxx
− θ (+)
0 (−) zz
0 (+) P (−)
z =Pzzxzz
Pxx
Pxx
Pxz
FFFFPPPzz0 zzxxxxxxzzzz
0! π/2 ! 2π! π/2 ! 2π(27)(+)
I (θ, ϕ) cos θ sin θ dθ dϕFz,i,j,k =I (θ,ϕ) (−)
cos θ (+)
sin 0θ dθdϕ(26) (−)=0 (+)
(+)
(−)
(−)
(+)
),k (+)
(−)
(+)
(−)
(+)
(−)
0 z
Fz PzzFzPzz Fx Fx Pxx Pxx Pzz Pzz! π/2 !
Fx Fx 0 Pxxx
Pxx!
!
π/2! π/2 ! 2π2π

(26)
2π(+)
(27)

F)θ sin θ dθ(+)dϕI (θ, (E,ϕ) cos(26) Fz,i,j,k =I (θ, ϕ) cos θ sin θ dθ 0dϕ=00(+)(−)(+)(−) Fz,i,j,k(+)
(−)
(+)
(−)
z PF(−)F(+)
Fx Pxx Pxx P0zz P0 zz(+)
zz PF(−)xPx
P(27)
xxxxzzzz! π/2 ! 2π0Fz

I (θ, ϕ)
(26)
(27)

+)
(−)
(+)
(−)
PPPxxxzz x zzz F (+) F (−) F + F − P (+) P (−) P (+) P (−)
zzxxxx(+)xx (−)
zz
zz+x z Fz

(27)

(+)
(−)
(+)
(−)
Fx Fx− PxxPxxPzzPzz

Numerical Flux Evaluated by Reconstruction!13

Numerical Flux Modified by Absorption
and Emission
F= F+ Femission!(+)
ν,x,i+1/2,j,k

ν,x,i+1/2,j,k!(−)ν,x,i+1/2,j,k−∆τi /2(+)
Fν,x,i+1/2,j,kPν,xy,i+1/2,j,k =" Sν
e+ 1−e4!
"−∆τi+1 /2 (+)−∆τi+1 /2 Sν
eFν,x,i+1/2,j,k − 1 − e
4!(+)
!(−)
Pν,xx,i+1/2,j,k + Pν,xx,i+1/2,j,k!
"−∆τi /2 (+)−∆τi /2 Sν
ePν,x,i+1/2,j,k + 1 − e6!
"−∆τi+1 /2 (−)−∆τi /2 SνePν,x,i+1/2,j,k + 1 − e6!(+)!(−)
Pν,xy,i+1/2,j,k + Pν,xy,i+1/2,j,k!(+)
Pν,x,i+1/2,j,k= e−∆τi /2(+)
Pν,x,i+1/2,j,k= e−∆τi+1 /2!(+)
Fν,x,i+1/2,j,k=!(+)

Fν,x,i+1/2,j,k =Pν,xx,i+1/2,j,k =!(+)
Pν,x,i+1/2,j,k=!(−)

Pν,x,i+1/2,j,k =!(−)
Pν,x,i+1/2,j,k(−)
Pν,x,i+1/2,j,kabsorption!

−∆τi /2

Modification due toEmission and(12)Absorption!14

!(+) of Absorption −∆τ(+)by theEffects
Evaluated
i /2!(+)=!(−)FFe+
ν,x,i+1/2,j,k
ν,x,i+1/2,j,k
Fν,x,i+1/2,j,k= FormalFν,x,i+1/2,j,k+F
ν,x,i+1/2,j,kSolution!(+)
Fν,x,i+1/2,j,k!(+)=!
"(+)−∆τi /2−∆τi /2 Sν
eFν,x,i+1/2,j,k−e−∆τi+1+/2 1 (+)

Fν,x,i+1/2,j,k = e

Fν,x,i+1/2,j,k4 −!
"S!(+)(+)−∆τ/2−∆τ/2FluxatBoundary
Absorption
FluxatCenterEmissivity
i+1i+1!(+) !(+)!(−) F !(−)
F=e−1−e
ν,x,i+1/2,j,k
ν,x,i+1/2,j,k= F= Fν,x,i+1/2,j,k+ F + Fν,x,i+1/2,j,k,x,i+1/2,j,k!(+)!(−)
1/2,j,kP=P+P
ν,xx,i+1/2,j,k!(+)
!(−)
!
"!(+)
!(−)
! ν,xx,i+1/2,j,k"−∆τ
ν,xx,iSν/2 Sν(+)−∆τ(+)/2 (+)
ν,x,i+1/2,j,k

ν,x,i+1/2,j,k−∆τi /2 =−∆τi /2i FiFν,x,i+1/2,j,k=+FPν,xx,i+1/2,j,kP+P
ν,x,i+1/2,j,k
ν,x,i+1/2,j,k=eF+
1−eν,xx,i+1/2,j,k
ν,xx,i+1/2,j,k=eF+
1−e,x,i+1/2,j,k
ν,x,i+1/2,j,k
1/2,j,k
ν,x,i+1/2,j,kGood Approximation
at a"44!
S!
ν " Sν
!(+)
(+)
−∆τ/2
−∆τ
/2!(+)(+)
ii!(+)(+)−∆τ/2large!
"!
"−∆τ(+)−∆τiS/2ν
=eF+1−ei
i /2Sν
(+)(+) P FP−∆τi+1/2−∆τ/2−∆τ/2−∆τ/2
ν,x,i+1/2,j,kν,x,i+1/2,j,k=Fi+1e F =Pν,x,i+1/2,j,k+e 4
i+11 −i+1eP+=e−1−eν,x,i+1/2,j,k=e−1−eν,x,i+1/2,j,k1/2,j,k
ν,x,i+1/2,j,k,x,i+1/2,j,kν,x,i+1/2,j,k
ν,x,i+1/2,j,k6
4
4!
"Sν
!(+)
−∆τi+1 /2 (+)
−∆τ/2
i+1!
"!(+) !(+)!(−)
SFν,x,i+1/2,j,k=eF−
1−e!(−)
ν!(−)(−)−∆τ/2−∆τ/2
ν,x,i+1/2,j,ki+1i=
P+P=ν,xx,i+1/2,j,kPν,xx,i+1/2,j,k+ PPν,x,i+1/2,j,k1/2,j,k4
Pν,x,i+1/2,j,k=e+1−exx,i+1/2,j,k
ν,xx,i+1/2,j,k
ν,xx,i+1/2,j,k!(−)(−)−∆τ/2!(+)
!(−)
i+16!
"
P=P+PP=eP+S!
"
ν,xx,i+1/2,j,kν
(+)
S−∆τ
/2
−∆τ
/2
ν,xx,i+1/2,j,k
ν,xx,i+1/2,j,kν
(+)(+)
ii!(−) −∆τν,x,i+1/2,j,k−∆τν,x,i+1/2,j,ki /2 !(+)
i /2=eP+1−e= e =ν,x,i+1/2,j,kPν,xy,i+1/2,j,k++ P
1 ν,xy,i+1/2,j,k−e 6
1/2,j,kPν,xy,i+1/2,j,kP,x,i+1/2,j,k
ν,x,i+1/2,j,k!
" Sν!(+)
(+)

6P
=eP+
1−e!
"
ν,x,i+1/2,j,k
ν,x,i+1/2,j,k!(+)
!(−)
S!(+)
(+)
!
"ν(−)
−∆τ/2
−∆τ/2−∆τ/2Si6
i+1i(−) P = e−∆τ/2 (−)−∆τi /2+ νP=ePi+1P+1−eP=P=eP+1−eν,x,i+1/2,j,kν,x,i+1/2,j,k1/2,j,k
ν,x,i+1/2,j,k
ν,xy,i+1/2,j,k!15!
"
ν,xy,i+1/2,j,k,x,i+1/2,j,k
ν,x,i+1/2,j,k
6!(−)
(−)
−∆τi+1 /2−∆τi /2 6 Sν ν,xy,i+−∆τi /2−∆τi /2

∂K
∂τ

= H

∂H
Aboid Extremely Sharp
Beam= J −S = 0∂τ∂KJ = 3H [τ + q (τ )]= H∂τ q " 2/3
flim E FF →|F q| (τ )]3H[τ+If

J =q " 2/3|F | > flim E
flim E F
F →|F |!16

Beam Test
Reflection Boundary
Incident
Radation

S=0Δτ = 0

Initial(HLL) Free
StreamingΔτ = 5Reflection

Shadow
simple HLL

this work!18

Shadow
Test (2)!19

Beam Test
(1)
Propagation

First OrderAccurate
Δx = 0.1!20

Beam Test(2)f = 0.999

f = 0.700!21

BeamTest (3)!22Beam Test (5)Burst Test (1)E = const.
F=0!24

Burst(2)Sphericity and Sharpness
Conflict!25

Burst(3)

E = const.
F = (0.9, 0, 0) E

arrows: f = F/E
HLL: diffusive, bias

!26

(+) (−)−∆τ−∆τ
i /2i+1 /2eePν,x,i+1/2,j,kPν,x,i+1/2,j,ke−∆τi+1 /2(12)1D Plane Parallel
(12)

(−)
Pν,x,i+1/2,j,kClosure relation∂H= J −S = 0∂τ∂H= J −S = 0
∂K∂τ
= H
∂τ
∂K= H
∂τ

J = 3H [τ + q (τ )]q " 2/3

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

Boundary:
No incident from outside!27

t, n) dn(3)

IR from Irradiated(4)
Protoplanetary Disk(x, t, n) dn

Iν (r, t, n)IR!
g (n, n! ) Iν (r, t, n! ) dn! (5)Optical (UV)

4πSν − cEν )(σν + σν,s ) F ν

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

Temperature Distribution Depends on the Incident Angle!29

Summary• M1 model can handle shadow andscattering, both of which have important
effects in celestial bodies.
• Reconstruction gives us simple but good
numerical fluxes.
– If Δt<Δx/cm, |F| < E– Burst Test
– Light Propagation
Dullemond & Monnier ’10!30

Numerical Instability in2nd Order Accuarate FluxCharacteristics

ES a new RHD
code
Gonzalez+’07431

Characteristics: Averagevelocity!32