SPH粒子データを直接用いた
Lya輻射輸送計算コードの開発

安部牧人（筑波大）!

共同研究者
長谷川賢二（名古屋大）、Benoit Semelin (パリ天文台)、矢島秀伸（東北大）、梅村雅之（筑波大）!

第4回 銀河進化研究会2017@大阪大学 2017/6/8

Lya輝線とLya輝線銀河•

Lyman alpha (Lya) 輝線••

星形成銀河の強力なプローブ(Partridge & Peebles 67ʼ、大質量星によるガスの電離 & 再結合で光る)

Lya輝線で非常に明るい銀河（Lyman alpha emitter, LAE）•

high-zで多く同定 、銀河進化の観点で重要な天体•

脱出光子スペクトルには銀河内部の情報（星形成率、ガスの速度場の情報）

asymmetric proﬁle (shifted to bluer)•

LAE: 理論的な理解が今後の課題•

Lya光子の脱出過程-•

Lya emissivity-•(Yajima+15ʼ)Lya輻射輸送計算（取り扱いが難しい）ISMの状態を正確に求める必要

詳細な銀河モデル + Lya輻射輸送計算（RT）

SPH計算による銀河モデル•

銀河形成シミュレーションにおける流体計算•⇢

SPH法（粒子法）、広いダイナミックhレンジを取り扱うことができる•

近傍粒子の重ね合わせ

SPH+UV輻射輸送計算（輻射流体力学）による銀河形成シミュレーション•

流体計算(SPH)+自己重力(DM)+非平衡化学進化+輻射輸送を同時に解き、ダイナミクスを追うこと（輻射流体計算）が可能になってきた•

SPH粒子を輻射輸送計算のグリッドに•

上流 j粒子とi粒子の情報からΔτijの計算
輻射流体計算ray-tracing in SPH (e.g., Susa 06ʼ)•

宇宙論的銀河形成シミュレーション(Hasegawa&Semelin 13ʼ)星からのUV feedback込みの計算

Lya輻射輸送計算を組み合わせた計算が課題

Lya輻射輸送（Radiative transfer, RT）
Voigt proﬁle✓ Lya輝線：2P→1Sの遷移

Doppler coreφH(x)•104-52ΔνD

A係数が非常に大きい（6.26×108 s-1）
➡ 共鳴散乱過程

Lorentz wing•

多重散乱して系外に脱出✓ ラインプロファイル（Voigt proﬁle）•
Doppler幅ΔνDで規格化した振動数 xv||Hd~0νout

•

~vv?

macroscopic⌫atom = ⌫ind~νin✓ 振動数シフト

microscopic

1⌫out = ⌫atom / 1~v · d~c~v · d~0c

coreで非常に光学的に厚くてもwingの光子は系から脱出できる振動数の取り扱い（振動数シフト）が重要~ · d~
Vc!

~ · d~0Vc!

•

原子静止系：coherent散乱•

実験室系：散乱前後で光子の方向が変わると、Doppler shiftにより振動数が変わる✓ ダストがある場合は、ダストによる吸収/散乱Monte Carlo法によるLya輻射輸送計算
Lya輻射輸送を輻射輸送方程式から解くのは大変（散乱+振動数シフト）••Monte Carlo法がよく採用される•

Lya光子がガスと相互作用せずに透過する確率は e-τ➡ 光子が散乱される場所（散乱するまで進めるτ）を確率的に与える

手続き•(1)target τ(光子が進む距離)をランダムに決定する
τtarget = -lnξ (ξ：一様乱数)(2)τ=τtargetになるまでtauを積分（ray-tracing）

τtargetτtarget(3)τ=τtargetで散乱（方向はランダムに決める）➡ 振動数シフトを評価(4)あらたなtarget τを決め、(1)に戻る•

光源から多数の光子パケット(~105)を飛ばし、Lya光子の散乱過程を追う

SPH based Lyα輻射輸送コードの開発
これまでのLya RT Monte Carloコードはメッシュベース•

mesh-assignによってSPH計算の情報を人工的に落とす

•

Lya光子

mesh-based : SPHの情報をmesh-assign
meshの物理量からLya emissivity決定、Lya RT計算
•SPH-based：高密度領域を高解像度で分解

SPH計算の情報をそのまま使ってLya輻射輸送を解きたい（mesh-free）
SEURAT: SPH scheme Extended with UV line RAdiative Transferを開発•

SPH粒子を輻射輸送計算の際のgridとして用いる➡

高密度領域を高解像度で分解したまま、Lya輻射輸送計算を行うことができる

SPH based ray-tracing in SEURAT✓

τの計算 in meshfree

振動数シフト•光線上に一番近い”下流”のSPH粒子を探してくる••

物理量をray上にassign

τ > τtarget になるまでrayに近い粒子の探索 & τ積分•

2粒子の中間点でtauを評価： ⌧i+ 21 = ⌧i 12 + xi ni (si+ 2112)

i 粒子の物理量でtau積分!
•si

散乱位置の評価：⌧target = ⌧i 12 +xi ni (stargetsi12) （図だとi = 3）

•starget (τ = τtarget)で散乱パケットの新たな進行方向を決めて再びray-tracing•ray-tracing中に下流粒子が見つからなければ光子は系から脱出したとみなす

Meshfree Lya radiative transfer•

SEURATで標準的なテスト問題を解く：
1)emerging spectrum from static sphere2)emerging spectrum from expanding sphere3)escape fraction from dusty slab

643 SPH粒子、光子パケット数=105

Emerging spectrum from static sphere
脱出光子スペクトルの解析解が知られている
98p=x2⇡ <ihpJ(x) = p:24a⌧0 1 + cosh2⇡ 3 /27(|x3 |/a⌧0 ) ;(Dijkstra+06ʼ)0.025

o0 = 104

0.02

J(x)

0.015
o0 = 1050.01τ0 = 104, 105, 106

o0 = 106

0.005

T= 10 K
0-60-40-20

x = (⌫0x⌫ )/0204060⌫D✓

多重散乱に伴う振動数シフトによってwing部分から系外に脱出（double peak）✓

光学的に厚い系ほど広がったスペクトル✓

解析解とよく一致

Emerging spectrum from expanding sphere
脱出光子スペクトルの解析解は無いが、よくテスト計算に用いられる✓bulk motionによってblue側(x > 0)の光子がガス静止系でline centerに➡

red側(x <0)から脱出してくるような非対称な形のemerging spectrum
SEURAT

LICORICE (mesh-based, Semelin+07ʼ)

0.2

2×1018cm20.15

J(x)

NH
0.12×1019cm22×1020cm20.05

V (r) =✓rrmax◆0-60vmax

NH = 2×1018-20 cm2-50-40x = (⌫-30⌫x0 )/-20⌫D-100mesh-based codeの結果とよく一致Dust scattering/absorption•dust : 脱出光子スペクトルの形に影響•

dust込みのoptical depth：

多重散乱するほどdust吸収の確率が高くなる

散乱

d⌧x =

吸収
22nds+⇡rQnds+⇡rx HId s,x dd Qa,x nd ds

d⌧x,HInd d ds = d⌧x,dQs,x , Qa,x : ダスト散乱/吸収のQ-value, nd = nH

✓mHmd◆✓

MdMH◆

光子パケットが中性水素に散乱される確率p=•⌫ nHI /( ⌫ nHI + nd d )コード内でのdustによる吸収の取り扱い
1. stochastic elimination of photon packets (光子パケットを消す)•

Dustによって散乱する確率： $ = Qs /(Qs + Qa )➡

1

$ で光子パケットを消滅させる

2. attenuation of ﬂux•

Lya光子がdustに吸収されない確率は e ⌧d,abs (τd,abs: dust吸収のτ)➡

光子がすすむ度に、packetの光子数を減衰させる(Laursen+09ʼ)Escape fraction from a dusty slab•

dusty slabからのLya光子脱出確率•

解析解が知られている（Neufeld 90ʼ）0τ0

T= 10 Klog fesc-1
elimination
attenuation
analytic solution-21ihfesc =pcosh ⇠ 0 (b⌧0 )1/3 ⌧d,abs•0⇠ =p3/(0.525⇡ 5/12 )-3-10
log (bo0)1/3od,abs

1“stochastic elimination”, “attenuation of ﬂux”ともに解析解をよく再現

Meshfree Monte Carlo based Lya RT•

ray-tracingで下流粒子が見つからない→光子パケットを脱出と判断➡

•

うまく粒子を見つけないと光子の脱出を誤って評価、正しくLya RT計算できない

meshfree Monte Carlo計算が困難になるケース：

1.

ローカルな光学的厚みが非常に大きい•Δτ>>1のとき、τ = τtargetは1粒子の中で満たされる パケットが1粒子の周りを多重散乱•

パケットの位置によっては、neighbor listから下流粒子を見つけられない➡

neighborのneighborまで粒子探索2.

ガス密度の非一様性が大きい•

高密度領域の粒子は低密度領域の粒子をneighbor listに含めない•

低密度領域の粒子は、高密度領域の粒子を近傍粒子リストに含めている場合がある

hiij➡ scatter list (kをneighbor listに含めている粒子)から粒子探索kn

Adaptivity of SEURATSEURATは銀河形成シミュレーションに適用可能か••一様密度といった単純な密度構造ではない中で、近傍粒子探索が正しく行われるか（meshfree Monte Carlo計算が可能か）検証

Ray-tracing in the clumpy cloud•

密度コントラスト; nclump /nf = 100, 1000, 10000( nclump /nf = 1000 の時、clumpのサイズとﬁeld粒子の平均間隔が同じくらい)1

100.5z/rcloud

01-0.50.1-1-1-0.5

0x/rcloud 0.5

0.5

1 -1

0-0.5 y/rcloud

1n[cm-3]••

cloud中心に光源を配置•

100,000の光子パケットを飛ばす•“脱出した”（ray-tracingをやめた）光
子パケットの位置を検証

Spatial distribution of “escaped” photon
光子パケットが“脱出した”と判断した位置の分布
5510510with scatterwith scatterw/o scatternclump/n = 102

00.2

0.4

0.6

count

count

count

310

104

310

nclump/n = 1031020.81101

00.2

r/rcloud•0.40.6

310

nclump/n = 104

102

0.8

1101

r/rcloud

00.2

0.4

0.6

0.8

r/rcloud“neighbor of neighbor” のみの場合 (w/o scatter particle list)•

高密度領域の粒子が低密度領域の粒子を見つけられない••

w/o scatter

104

102

with scatterw/o scatter

104

101

10

光子パケットの”脱出”を誤って評価

scatter listを考慮した場合-

密度コントラストによらず、全パケットが正しく系外に脱出-

scatter listを近傍粒子探索の際に考慮することで、密度構造が複雑な系に対してもmeshfreeでLya輻射輸送計算が可能

1

まとめ•

LAEの理論的理解には銀河形成シミュレーションにLya RTを組み合わせた計算が必要•

銀河形成シミュレーションにはSPH法が広く用いられるが、Lya RT計算コー
ドはメッシュベースで開発されてきた•SPH粒子を輻射輸送計算の際のグリッドとして用いる、meshfreeのMonte Carlo
Lya RTコードを開発した
SEURAT: SPH scheme Extended with Ultraviolet line RAdiativeTransfer•

SEURATでは、水素原子によるLya光子の散乱に加えダスト吸収/散乱過程を実装、テスト計算の結果、meshfreeでこれらの過程が正しく解けることを確認•“neighbor of neighbor”、“scatter particle list”を用いた近傍粒子探索により、光学的に非常に厚く、また密度の非一様性が大きい系でのray-tracingが可能