朱佳麗

(1. 中國科學(xué)院空間天文與技術(shù)重點(diǎn)實(shí)驗(yàn)室，北京　100012；2. 中國科學(xué)院國家天文臺(tái)，北京　100012)

?

塵埃三維蒙特卡洛輻射轉(zhuǎn)移模型述評(píng)*

朱佳麗1,2

(1. 中國科學(xué)院空間天文與技術(shù)重點(diǎn)實(shí)驗(yàn)室，北京100012；2. 中國科學(xué)院國家天文臺(tái)，北京100012)

通過觀測(cè)和模擬星際空間的輻射轉(zhuǎn)移過程，探索天體的物理化學(xué)性質(zhì)和分布演化過程，是天體物理學(xué)研究的一種重要手段。星際塵埃在星際輻射的產(chǎn)生和加工方面都扮演了重要角色。星際塵埃的分布為三維不均勻分布，需采用塵埃三維輻射轉(zhuǎn)移模型以更真實(shí)地模擬塵埃的輻射轉(zhuǎn)移過程。在研究分析塵埃三維輻射轉(zhuǎn)移模型有關(guān)文獻(xiàn)的基礎(chǔ)上，重點(diǎn)介紹了用蒙特卡洛方法模擬塵埃三維輻射轉(zhuǎn)移的基本圖景和對(duì)模擬結(jié)果的一般處理步驟。詳細(xì)梳理了適用于不同物理環(huán)境的6 種模擬塵埃三維輻射轉(zhuǎn)移的開源代碼，比較分析了它們的模擬對(duì)象、輻射源類型、塵埃分布和組成等要素，總結(jié)評(píng)析了各個(gè)代碼的特點(diǎn)和優(yōu)勢(shì)，便于需要利用數(shù)值模擬來進(jìn)行科學(xué)研究的相關(guān)人員根據(jù)科學(xué)目標(biāo)進(jìn)行鑒別和采用。

塵埃；蒙特卡洛方法；數(shù)值模擬；三維輻射轉(zhuǎn)移

銀河系除了會(huì)發(fā)光的恒星之外，星際空間被自身不發(fā)光的星際介質(zhì)占據(jù)。星際介質(zhì)由氣體和塵埃組成，塵埃吸收星光能量，并將其轉(zhuǎn)化為紅外波段的熱輻射，同時(shí)通過光電作用加熱氣體。通過多波段的觀測(cè)數(shù)據(jù)研究恒星等輻射天體和星際介質(zhì)，需要了解輻射在通過介質(zhì)時(shí)的轉(zhuǎn)移情況，以獲得輻射源和介質(zhì)的幾何、物理和化學(xué)性質(zhì)。由于這個(gè)輻射轉(zhuǎn)移過程相當(dāng)復(fù)雜，一般通過構(gòu)建數(shù)值模型來模擬和研究，即在一定的模型輸入?yún)?shù)下，將數(shù)值模型的模擬結(jié)果與實(shí)際觀測(cè)數(shù)據(jù)相比較，再通過擬合觀測(cè)數(shù)據(jù)反演數(shù)值模型的輸入?yún)?shù)，從而掌握與輸入?yún)?shù)相關(guān)的天體物理化學(xué)性質(zhì)。

塵埃的分布為三維不均勻分布，為了更真實(shí)地模擬塵埃的輻射轉(zhuǎn)移過程，需采用三維輻射轉(zhuǎn)移模型。然而，塵埃三維輻射轉(zhuǎn)移問題有多方面的難點(diǎn)[1]。首先它是一個(gè)非本地非線性的六維問題，包括了輻射場(chǎng)的坐標(biāo)(三維)、傳播方向(二維)和輻射波長這6個(gè)參數(shù)。其次，塵埃三維分布的復(fù)雜度越高，即模型所模擬的三維空間的幾何分辨率越高，模型計(jì)算所需的存儲(chǔ)空間越大，計(jì)算難度也越大。

三維輻射轉(zhuǎn)移問題通常使用蒙特卡洛方法解決。雖然也有其它的技術(shù)和方法被開發(fā)來計(jì)算輻射轉(zhuǎn)移，但基本是對(duì)輻射轉(zhuǎn)移方程直接計(jì)算數(shù)值解。這種解方程的數(shù)值方法能有效地解決一維問題，但在二維和三維問題上變得十分復(fù)雜[2]，更無法有效解決三維輻射轉(zhuǎn)移問題，或者存在顯著的缺陷[1](如忽略散射和塵埃熱輻射等)。蒙特卡洛方法非常適用于任意的塵埃三維密度分布，因?yàn)樗皇侵苯佑?jì)算輻射轉(zhuǎn)移方程，而是利用對(duì)概率密度函數(shù)(Probability Distribution Function, PDF)的隨機(jī)抽樣，達(dá)到在三維空間中傳播大量光子從而獲得統(tǒng)計(jì)結(jié)果的目的[2]。

本文第1節(jié)介紹了基于蒙特卡洛方法的塵埃三維輻射轉(zhuǎn)移的基本模擬過程。針對(duì)不同物理環(huán)境的模擬，在第2節(jié)分析述評(píng)了6個(gè)開源代碼，并概括了它們的模擬性能。

1　用蒙特卡洛方法模擬塵埃的三維輻射轉(zhuǎn)移

在蒙特卡洛算法中，輻射場(chǎng)不斷地發(fā)出大量的光包，每一個(gè)光包由一定數(shù)量的光子組成。光包沿著一定的路徑穿過塵埃介質(zhì)，直到被吸收或離開介質(zhì)。在此過程中的每一階段，光包的下一段路徑的特征參數(shù)都是基于適當(dāng)?shù)母怕拭芏群瘮?shù)隨機(jī)抽樣確定的，具體的抽樣方法可參見文[3]。通過發(fā)射足夠數(shù)量的光包，并對(duì)所有光包的路徑進(jìn)行統(tǒng)計(jì)分析而得到輻射場(chǎng)最終的分布情況?？梢?，蒙特卡洛方法是通過模擬輻射轉(zhuǎn)移過程代替解輻射轉(zhuǎn)移方程。

1.1基本模擬過程

在最簡化的情況下，輻射源發(fā)射某個(gè)頻率的光子，光子在具有一定密度分布的三維網(wǎng)格(即離散化的三維空間*基于一定的坐標(biāo)系，將要模擬的三維空間離散化，形成三維網(wǎng)格。將模型參數(shù)(如塵埃密度和輻射強(qiáng)度等)隨網(wǎng)格離散化，且參數(shù)值在同一個(gè)單元格內(nèi)保持不變。)中傳播，每一個(gè)光子在被散射或吸收(隨即被再發(fā)射)之前都會(huì)經(jīng)過一定的光深，這些光子將繼續(xù)在網(wǎng)格中傳播直到它們離開網(wǎng)格[2]。離開網(wǎng)格的光子可用來生成光譜能量分布(Spectral Energy Distributions, SEDs)和圖像。通過發(fā)射大量的光子重復(fù)以上過程，將提高模擬所得的光譜能量分布和圖像的信噪比。

1.1.1光包的發(fā)射和傳播

為了提高模擬效率，輻射源實(shí)際發(fā)射的是光包而不是單個(gè)光子。在有多個(gè)源的情況下，光包來自哪個(gè)源隨機(jī)抽樣決定。每個(gè)源發(fā)射的光包數(shù)目可以和源的光度成正比，也可以假定每個(gè)源都發(fā)射同樣數(shù)目的光包，但不同源發(fā)射的光包的權(quán)重不同[2]。所發(fā)射光包的傳播方向、頻率隨機(jī)抽樣決定。

光包發(fā)射后開始傳播過程。對(duì)光包將要經(jīng)過的光深進(jìn)行抽樣，并和光包當(dāng)前路徑上的實(shí)際光深比較，當(dāng)抽樣光深大于實(shí)際光深時(shí)，光包將離開網(wǎng)格；反之，光包將和塵埃相互作用，具體是被散射或被吸收則由塵埃的發(fā)射率決定。被吸收后再發(fā)射的光包頻率通過隨機(jī)抽樣獲得。光包被散射或被吸收再發(fā)射后通過隨機(jī)抽樣得到新的傳播方向和新的光深值，并重復(fù)上述傳播過程直到光包離開網(wǎng)格。

蒙特卡洛方法跟蹤大量光包的傳播軌跡，光包在其存在的每個(gè)階段都有一系列的量來標(biāo)注其屬性，比如坐標(biāo)、傳播方向和波長，根據(jù)研究目標(biāo)的不同，還可以額外記錄諸如多普勒位移和偏振等信息[4]。當(dāng)光包和三維網(wǎng)格中的介質(zhì)發(fā)生相互作用，光包的上述屬性將發(fā)生改變，并通過對(duì)適當(dāng)?shù)母怕拭芏群瘮?shù)進(jìn)行隨機(jī)抽樣更新光包的屬性。當(dāng)光包離開網(wǎng)格，它的最終屬性被記錄下來。通過這種方式記錄大量的光包信息，將得到被模擬系統(tǒng)的總體統(tǒng)計(jì)特性。

1.2模擬結(jié)果的處理

上述被模擬系統(tǒng)的總體統(tǒng)計(jì)特性即為模擬結(jié)果。在不同模型參數(shù)設(shè)置下進(jìn)行多次模擬，并將模擬結(jié)果和觀測(cè)數(shù)據(jù)相比較可獲知觀測(cè)對(duì)象的性質(zhì)，一般步驟如下[1]：

(1)將塵埃密度分布、輻射源性質(zhì)、塵埃光學(xué)性質(zhì)等進(jìn)行參數(shù)化；

(2)確定模型網(wǎng)格的設(shè)置，將步驟(1)中的參數(shù)在網(wǎng)格中離散化，使用輻射轉(zhuǎn)移代碼計(jì)算得到光譜能量分布或者圖像等模擬結(jié)果；

(3)比較模擬結(jié)果和觀測(cè)數(shù)據(jù)；

(4)評(píng)價(jià)步驟(3)中兩者的差距，改變模型的參數(shù)和(或)模型的假設(shè)，重復(fù)步驟(2)、(3)、(4)，直到模擬和觀測(cè)的差距最小化；

(5)找到最優(yōu)的模型參數(shù)設(shè)置。

其中，合理評(píng)價(jià)觀測(cè)和模擬誤差對(duì)步驟(3)中的比較十分重要，因?yàn)樗苯佑绊懖襟E(4)中對(duì)模型參數(shù)和假設(shè)的選擇，最終影響最優(yōu)模型設(shè)置的獲得。

2　塵埃三維輻射轉(zhuǎn)移模型代碼及其特點(diǎn)分析

計(jì)算天體物理學(xué)家對(duì)一維輻射轉(zhuǎn)移模型代碼的開發(fā)使用已有近30年的經(jīng)驗(yàn)，而三維輻射轉(zhuǎn)移模型卻是近些年隨著計(jì)算機(jī)運(yùn)算能力的提升才發(fā)展起來的。到20世紀(jì)90年代中期，還沒有輻射轉(zhuǎn)移模型代碼可以實(shí)現(xiàn)一般三維幾何體中星際塵埃對(duì)輻射的吸收、散射和塵埃熱輻射的計(jì)算[1]。可用的代碼不是只適用于一維或二維幾何體，就是不能計(jì)算完整的三維輻射轉(zhuǎn)移問題，如忽略散射或塵埃熱輻射。

在過去的20年間，計(jì)算機(jī)運(yùn)算能力的顯著增強(qiáng)和更高效算法的開發(fā)，促使了三維塵埃輻射轉(zhuǎn)移模型代碼的產(chǎn)生，這些代碼，具備了計(jì)算三維幾何體中塵埃對(duì)輻射的吸收和各向異性散射的能力。這些代碼的應(yīng)用范圍十分廣泛，可適用于原恒星、活動(dòng)星系核(Active Galactic Nucleus, AGN)、分子云和星系等。

這些代碼一般用Fortran、C或C++語言編寫，同時(shí)為了提高運(yùn)算效率，提供了并行化運(yùn)算功能。代碼的基本輸入?yún)?shù)包括：塵埃三維分布情況、塵埃的物理化學(xué)性質(zhì)、輻射源坐標(biāo)和物理化學(xué)性質(zhì)以及控制代碼運(yùn)行和輸出的一些參數(shù)。

三維塵埃輻射轉(zhuǎn)移模型中，塵埃顆粒的最簡單假設(shè)是單一種類的處于局部熱平衡的塵埃顆粒。這種情況下，塵埃輻射為灰體輻射，每個(gè)單元格的塵埃的平衡溫度由吸收和輻射平衡決定。第2種更符合現(xiàn)實(shí)的做法是，考慮處于熱平衡狀態(tài)的不同種類和尺寸的塵埃，即將不同種類和尺寸的塵埃顆粒分開處理，放入不同的 “容器”(bin)內(nèi)，每一個(gè)容器都處于局部熱平衡，塵埃的總輻射是一系列灰體輻射的總和，其中每個(gè)容器的溫度由吸收和輻射平衡決定。第3種是最接近實(shí)際情況的做法，即考慮來自小尺寸塵埃和多環(huán)芳香烴(Polycyclic Aromatic Hydrocarbons, PAH)的非局部熱平衡輻射[5]。非局部熱平衡輻射和局部熱平衡輻射的主要區(qū)別是塵埃溫度不是用單個(gè)的平衡溫度表征，而是用溫度分布表征。在實(shí)踐中，要在輻射轉(zhuǎn)移代碼中包含塵埃非局部熱平衡輻射的計(jì)算是一項(xiàng)艱難的工程。下文介紹的代碼都可以計(jì)算塵埃的熱平衡輻射，其中個(gè)別代碼還可同時(shí)計(jì)算非局部熱平衡輻射。

到目前為止，已公開發(fā)表的塵埃三維輻射轉(zhuǎn)移模型代碼有近30個(gè)[1]，其中90%為基于蒙特卡洛方法開發(fā)的代碼。本節(jié)介紹其中6個(gè)開源的蒙特卡洛輻射轉(zhuǎn)移模型代碼(見表1)，這些代碼基于不同的開發(fā)目的能針對(duì)性地模擬不同的科學(xué)問題。在模擬研究中，合理選擇模型代碼有利于得到更真實(shí)的天體信息。

表1　開源塵埃三維輻射轉(zhuǎn)移模型代碼

2.1HO-CHUNK的模擬性能和應(yīng)用范圍

原恒星一般由年輕恒星和包圍它的盤、下落的包層以及雙極外流組成。HO-CHUNK可模擬原恒星從Class 0到Class III階段的一系列演化過程中的光譜能量分布、顏色、偏振性質(zhì)和圖像[6]。

高分辨率和多波段的觀測(cè)揭示了原恒星盤存在復(fù)雜的三維結(jié)構(gòu)，例如彎曲、裂縫、旋渦狀的塵埃密度分布等。HO-CHUNK 可模擬盤的裂縫和旋渦結(jié)構(gòu)、內(nèi)盤的彎曲變形、中央恒星上的吸積熱點(diǎn)、塵埃密度分布的不規(guī)則團(tuán)塊結(jié)構(gòu)、包層的裂縫，以及雙極外流的多重空腔壁[7]等三維結(jié)構(gòu)，使模型更符合實(shí)際情況。

HO-CHUNK 采用球極坐標(biāo)系，網(wǎng)格結(jié)構(gòu)的設(shè)置允許一個(gè)單元格內(nèi)有多種塵埃成分。輻射源包括3類，分別為恒星輻射、盤吸積產(chǎn)生的輻射和星際輻射場(chǎng)。塵埃盤分為內(nèi)盤和外盤，可分別具有不同的幾何結(jié)構(gòu)和塵埃性質(zhì)，塵埃密度分布可為團(tuán)塊結(jié)構(gòu)。模型計(jì)算來自大尺寸塵埃的熱平衡輻射，以及來自小尺寸塵埃和多環(huán)芳香烴的非局部熱平衡輻射，但不計(jì)算塵埃的散射。該代碼預(yù)定義了多種塵埃模型，使用者也可以構(gòu)建自己的塵埃模型。

該代碼為信息傳遞接口(Message Passing Interface, MPI)并行運(yùn)算代碼，可在一臺(tái)計(jì)算機(jī)的多個(gè)處理器或者多臺(tái)聯(lián)網(wǎng)的計(jì)算機(jī)上運(yùn)行。針對(duì)模擬結(jié)果，HO-CHUNK提供IDL程序?qū)Y(jié)果進(jìn)行后續(xù)處理與分析。

2.1.1應(yīng)用和成果

不同演化階段的原恒星具有不同的塵埃輻射特征。通過使用HO-CHUNK模擬發(fā)現(xiàn)，Class 0階段原恒星中波長大于100 μm的紅外輻射來自其包層， 而Class I階段原恒星的這部分紅外輻射則主要來自塵埃盤[6]，這和前人的觀測(cè)及模擬結(jié)果一致[8]。

和低質(zhì)量的原恒星相比較，高質(zhì)量的原恒星的數(shù)量較少且演化時(shí)標(biāo)較短，導(dǎo)致難以區(qū)分實(shí)際觀測(cè)數(shù)據(jù)所對(duì)應(yīng)的演化階段，對(duì)進(jìn)一步分析和挖掘觀測(cè)數(shù)據(jù)造成障礙。使用HO-CHUNK模擬不同演化階段的高質(zhì)量原恒星的塵埃輻射特征和形態(tài)特征，并和觀測(cè)數(shù)據(jù)相比較，從而確定實(shí)際觀測(cè)數(shù)據(jù)所對(duì)應(yīng)的演化階段，同時(shí)也為鑒別觀測(cè)對(duì)象是否為原恒星提供有力的線索[7]。使用上述方法，文[9]從多波段巡天數(shù)據(jù)中鑒定出了約1 000個(gè)原恒星。

2.2Hyperion的模擬性能和應(yīng)用范圍

Hyperion可模擬任意三維空間的塵埃連續(xù)輻射[2]。該代碼基本是通用的，能適用于多種科學(xué)目標(biāo)，可計(jì)算塵埃溫度、光譜能量分布和多波段圖像。

該代碼由兩部分組成：一部分進(jìn)行輻射轉(zhuǎn)移的計(jì)算，由Fortran執(zhí)行，這部分代碼和具體科學(xué)問題無關(guān)；第二部分包括了一些Python庫，用于便捷地構(gòu)建輸入文件，并包含一些預(yù)定義的常用密度分布結(jié)構(gòu)，同時(shí)提供了模擬結(jié)果的后續(xù)處理工具。

Hyperion采用三維直角坐標(biāo)系、球極坐標(biāo)系和圓柱坐標(biāo)系。網(wǎng)格中可分布任意位置和數(shù)量的點(diǎn)源、一定半徑的球狀源、彌散源、星際輻射場(chǎng)和平面源，塵埃的分布也為任意三維分布。Hyperion計(jì)算來自大尺寸塵埃的熱輻射，并同時(shí)計(jì)算塵埃散射和偏振輻射。

Hyperion是信息傳遞接口并行代碼，并提供Python庫便于使用者對(duì)模擬結(jié)果進(jìn)行后續(xù)處理與分析。

2.2.1應(yīng)用和成果

Hyperion基于通用目的開發(fā)，但目前主要應(yīng)用于原恒星的研究。和HO-CHUNK對(duì)原恒星復(fù)雜三維結(jié)構(gòu)的預(yù)定義相比，Hyperion對(duì)原恒星的幾何結(jié)構(gòu)的設(shè)置更為靈活，被用于原恒星候選體的鑒定[10]和原恒星盤幾何結(jié)構(gòu)的確定[11]。

Class 0階段的原恒星的光譜能量分布明顯偏向紅端且在亞毫米波段的流量較大。文[10]將觀測(cè)所得原恒星候選體光譜能量分布和不同模型參數(shù)設(shè)置下Hyperion所預(yù)言的光譜能量分布相比較，獲知原恒星候選體的物理化學(xué)特征，從而確定候選體極可能為原恒星。文[11]則通過比較模型模擬結(jié)果和實(shí)際觀測(cè)數(shù)據(jù)來搜尋Class 0類型原恒星的盤結(jié)構(gòu)特征，發(fā)現(xiàn)Class 0階段的原恒星已可以具備盤結(jié)構(gòu)。2.33D Mocassin的模擬性能和應(yīng)用范圍

光致電離環(huán)境的研究是很多天體物理學(xué)問題的基礎(chǔ)。3D Mocassin用于模擬包含塵埃的光致電離星云[12]，能計(jì)算任意幾何結(jié)構(gòu)下被輻射場(chǎng)強(qiáng)電離的低密度氣態(tài)環(huán)境，輻射場(chǎng)波長范圍可從射電波段到伽瑪射線波段[13]。

3D Mocassin采用三維直角坐標(biāo)系，能處理任意的三維幾何結(jié)構(gòu)和密度分布。輻射源可同時(shí)包含恒星輻射和彌散輻射場(chǎng)，可處理任意數(shù)量和位置的點(diǎn)源以及位于X-Z平面的平面源。3D Mocassin計(jì)算大尺寸塵埃的熱輻射和散射。該代碼在計(jì)算星云的發(fā)射線光譜、電子溫度和密度結(jié)構(gòu)的同時(shí)，還可以輸出星云的光譜能量分布。

3D Mocassin是信息傳遞接口并行代碼。相比其余5個(gè)代碼，3D Mocassin的特色是能同時(shí)處理氣體和塵埃的輻射轉(zhuǎn)移，更符合星云中氣體和塵埃相互混合的真實(shí)情況。塵埃和氣體之間的相互作用包含了相互碰撞和光電效應(yīng)，每一個(gè)單元格內(nèi)的物理狀況通過解熱平衡和電離平衡方程獲得。根據(jù)不同的模擬目的，也可對(duì)塵埃和氣體分別計(jì)算，即構(gòu)建純塵?；蚣儦怏w模型。

2.3.1應(yīng)用和成果

3D Mocassin主要用于富含塵埃和氣體的行星狀星云和超新星的研究中。通過擬合光譜信息可以獲知星云的溫度、尺度和光度等信息[14]。

在超新星的研究中，3D Mocassin被用于計(jì)算不同塵埃質(zhì)量和分布情況下超新星的輻射轉(zhuǎn)移情況。通過擬合超新星隨時(shí)間演化的塵埃輻射，可以獲知超新星爆發(fā)后塵埃形成的時(shí)間和位置。文[15]在分析擬合了SN 1987A從爆發(fā)后64天到24年間的光譜能量分布數(shù)據(jù)后，發(fā)現(xiàn)大部分的塵埃形成于超新星爆發(fā)后的晚期。

2.4SKIRT的模擬性能和應(yīng)用范圍

SKIRT通過構(gòu)建細(xì)致的輻射轉(zhuǎn)移模型來研究盤星系中塵埃的吸收和散射對(duì)星系動(dòng)力學(xué)參數(shù)的影響[4]。模擬塵埃星系的動(dòng)力學(xué)參數(shù)時(shí)，模型需要記錄光包的多普勒位移。

SKIRT采用三維直角坐標(biāo)系、一維球坐標(biāo)系和二維圓柱坐標(biāo)系，并預(yù)定義了多種幾何結(jié)構(gòu)用于設(shè)置輻射源和塵埃的空間密度分布，包括解析的球?qū)ΨQ、軸對(duì)稱分布，以及非解析的不對(duì)稱分布。SKIRT可同時(shí)計(jì)算塵埃的熱平衡輻射和非局部熱平衡輻射。它結(jié)合了DustEM[16]這個(gè)開源的塵埃輻射的數(shù)值計(jì)算工具，使其可以快速準(zhǔn)確地計(jì)算被瞬時(shí)加熱的極小塵埃顆粒和多環(huán)芳香烴的非局部熱平衡輻射[5]。模型輸出光譜能量分布、星系動(dòng)力學(xué)參數(shù)和圖像。

SKIRT被設(shè)計(jì)成用戶友好架構(gòu)，即通過和用戶交互的方式構(gòu)建模型模擬所需的配置文件，然后通過命令行方式根據(jù)配置文件運(yùn)行模擬程序[17]。SKIRT是信息傳遞接口并行代碼。

2.4.1應(yīng)用和成果

側(cè)向旋渦星系的塵埃分布往往呈現(xiàn)特有的窄帶形態(tài)，是研究塵埃在星系中的徑向和垂直方向分布情況的理想對(duì)象。通過使用SKIRT模擬NGC 5166和IC 4225這兩個(gè)側(cè)向旋渦星系在遠(yuǎn)紅外波段塵埃輻射的光譜能量分布，并和觀測(cè)所得光譜能量分布相比較，文[18]發(fā)現(xiàn)年輕星族對(duì)加熱星系塵埃起了至關(guān)重要的作用，缺少年輕星族的星系模型所預(yù)言的塵埃輻射要低于實(shí)際觀測(cè)值的3倍。

SKIRT同樣適用于正向旋渦星系的模擬。文[19]結(jié)合紅外、紫外和光學(xué)的多波段數(shù)據(jù)，使用SKIRT模擬了M 51中塵埃的遠(yuǎn)紅外和亞毫米波段的輻射。通過比較模型模擬結(jié)果和觀測(cè)數(shù)據(jù)，分析了年輕星族和年老星族對(duì)加熱塵埃所作的貢獻(xiàn)，其中年老星族提供了加熱塵埃所需總能量的37%。

2.5STOKES的模擬性能和應(yīng)用范圍

由于光子被自由電子和塵埃顆粒散射，活動(dòng)星系核的光在很大的波長范圍內(nèi)是偏振的，這對(duì)輻射和散射區(qū)域的幾何結(jié)構(gòu)形成了重要的約束。STOKES計(jì)算活動(dòng)星系核統(tǒng)一模型所描述的輻射轉(zhuǎn)移，主要用于模擬活動(dòng)星系核中塵埃造成的偏振[20]，從而研究不同幾何形狀和光深的塵埃環(huán)對(duì)輻射偏振的影響。

STOKES采用三維直角坐標(biāo)系，模擬區(qū)域分為源區(qū)和散射區(qū)兩部分。源區(qū)可包括圓柱形的連續(xù)輻射源、圓柱形的寬線源和雙圓錐形的窄線源；散射區(qū)的塵埃分布結(jié)構(gòu)包括圓柱形、圓環(huán)、雙圓錐、喇叭形和球形。STOKES計(jì)算大尺寸塵埃的熱輻射和散射，輸出的主要模擬結(jié)果為光譜能量分布、輻射流量和偏振信息。模型還可記錄每一個(gè)光子從發(fā)射到離開網(wǎng)格所經(jīng)歷的時(shí)間，從而模擬和分析時(shí)變的偏振。

經(jīng)過多年開發(fā)，STOKES已被擴(kuò)展為可模擬寬線區(qū)發(fā)射線的偏振對(duì)該區(qū)域速度的依賴情況[21]。寬線區(qū)的運(yùn)動(dòng)主要是旋轉(zhuǎn)和明顯的湍流，很可能還有內(nèi)流。STOKES不能進(jìn)行并行運(yùn)算。

2.5.1應(yīng)用和成果

使用STOKES進(jìn)行模擬研究發(fā)現(xiàn)，塵埃環(huán)的形狀對(duì)造成光線偏振的程度輕重有顯著的影響[20]。通常認(rèn)為，活動(dòng)星系核最靠近中心的部分是由多種產(chǎn)生輻射和散射的介質(zhì)組成的，這些介質(zhì)在輻射轉(zhuǎn)移過程中耦合。光譜偏振觀測(cè)為這些區(qū)域的幾何結(jié)構(gòu)提供了重要的線索。通過構(gòu)建模型模擬各成分間的輻射轉(zhuǎn)移耦合獲得活動(dòng)星系核的偏振特性，并將模型模擬結(jié)果和偏振觀測(cè)結(jié)果相比較，從而約束活動(dòng)星系核的幾何結(jié)構(gòu)。

文[22]使用STOKES模擬了II型塞弗特星系NGC 1068的X波段的光譜偏振，分析了偏振特性隨輻射轉(zhuǎn)移發(fā)生區(qū)域和觀測(cè)角度的變化情況，預(yù)言了X波段的偏振觀測(cè)將對(duì)活動(dòng)星系核內(nèi)部幾何結(jié)構(gòu)的確定起到重要作用。

2.6RADMC-3D的模擬性能和應(yīng)用范圍

RADMC-3D定位為多用途輻射轉(zhuǎn)移工具，是基于通用目的開發(fā)的代碼。RADMC-3D可用于模擬任意三維空間塵埃的輻射轉(zhuǎn)移，比如可用于模擬原行星盤、活動(dòng)星系核、分子云和行星大氣等的輻射轉(zhuǎn)移*http://www.ita.uni-heidelberg.de/～dullemond/software/radmc-3d/。

針對(duì)不同的模擬對(duì)象，RADMC-3D預(yù)定義了一系列的模型設(shè)置作為樣本，方便用戶在樣本基礎(chǔ)上改進(jìn)模型參數(shù)。值得一提的是，RADMC-3D的使用文檔是本文介紹的代碼中最為完善的，并且提供專門的論壇便于使用者和開發(fā)者之間的交流*http://radmc3d.ita.uniheidelberg.de/phpbb/index.php。

根據(jù)到目前為止的開發(fā)情況，RADMC-3D采用三維直角坐標(biāo)系和球極坐標(biāo)系，輻射源包括任意數(shù)量和位置的點(diǎn)源、一定半徑的球狀源和星際輻射場(chǎng)。該代碼計(jì)算大尺寸塵埃熱輻射和散射，輸出模擬區(qū)域的光譜能量分布和圖像，并提供IDL程序便于使用者對(duì)模擬結(jié)果進(jìn)行后續(xù)處理與分析。RADMC-3D是OpenMP并行代碼。

2.6.1應(yīng)用和成果

作為一個(gè)新近開發(fā)的代碼，目前已發(fā)表的RADMC-3D的具體應(yīng)用較少。具有代表性的應(yīng)用是在赫比格Ae型星HD 139614的盤塵埃分布的研究，文[23]使用RADMC-3D模擬該恒星塵埃盤的輻射圖像和光譜能量分布。通過擬合觀測(cè)數(shù)據(jù)，確定了塵埃盤中裂縫的存在，并得到了盤中塵埃的物理化學(xué)性質(zhì)。

3　總　結(jié)

本文介紹了用蒙特卡洛方法模擬塵埃三維輻射轉(zhuǎn)移的基本過程，并列舉了6個(gè)基于蒙特卡洛方法的塵埃三維輻射轉(zhuǎn)移模型的開源代碼。HO-CHUNK、3D Mocassin、SKIRT和STOKES這4個(gè)代碼可分別用于模擬原恒星、光致電離環(huán)境、盤星系和活動(dòng)星系核，Hyperion和RADMC-3D則為通用型代碼，可模擬任意三維空間塵埃的輻射轉(zhuǎn)移。這些開源代碼的基本模擬性能特點(diǎn)已整理成表2，以便于相關(guān)研究人員根據(jù)研究目標(biāo)選擇使用代碼。

表2　開源代碼基本模擬性能概覽

[1]Steinacker J, Baes M, Gordon K D. Three-dimensional dust radiative transfer[J]. Annual Review of Astronomy and Astrophysics, 2013, 51: 63-104.

[2]Robitaille T P. HYPERION: an open-source parallelized three-dimensional dust continuum radiative transfer code[J]. Astronomy & Astrophysics, 2011, 536: A79-A97.

[3]Whitney B A. Monte Carlo radiative transfer[J]. Bulletin of the Astronomical Society of India, 2011, 39: 101-127.

[4]Baes M, Davies J I, Dejonghe H, et al. Radiative transfer in disc galaxies-III. the observed kinematics of dusty disc galaxies[J]. Monthly Notices of the Royal Astronomical Society, 2003, 343(4): 1081-1094.

[5]Baes M, Verstappen J, De Looze I, et al. Efficient three-dimensional NLTE dust radiative transfer with SKIRT[J]. The Astrophysical Journal Supplement Series, 2011, 196(2): 22-36.

[6]Whitney B A, Wood K, Bjorkman J E, et al. Two-dimensional radiative transfer in protostellar envelopes. II. an evolutionary sequence[J]. The Astrophysical Journal, 2003, 598: 1079-1099.

[7]Whitney B A, Robitaille T P, Bjorkman J E, et al. Three-dimensional radiation transfer in young stellar objects[J]. The Astrophysical Journal Supplement Series, 2013, 207: 30-49.

[8]Wolf S, Padgett D R, Stapelfeldt K R. The circumstellar disk of the butterfly star in Taurus[J]. The Astrophysical Journal, 2003, 588: 373-386.

[10]Stutz A M, Tobin J J, Stanke T, et al. A herschel and apex census of the reddest sources in orion: searching for the youngest protostars[J]. The Astrophysical Journal, 2013, 767(1): 1-36.

[11]Tobin J J, Looney L W, Wilner D J, et al. A sub-arcsecond survey toward class 0 protostars in perseus: searching for signatures of protostellar disks[J]. The Astrophysical Journal, 2015, 805(2): 125-182.

[12]Ercolano B, Barlow M J, Storey P J. The dusty MOCASSIN: fully self-consistent 3D photoionization and dust radiative transfer models[J]. Monthly Notices of the Royal Astronomical Society, 2005, 362(3): 1038-1046.

[13]Ercolano B, Young P R, Drake J J, et al. X-ray enabled MOCASSIN: a three-dimensional code for photoionized media[J]. The Astrophysical Journal Supplement Series, 2008, 175(2): 534-542.

[14]Wesson R, Barlow M J, Corradi R L M, et al. A planetary nebula around nova v458 vulpeculae undergoing flash ionization[J]. The Astrophysical Journal, 2008, 688: L21-L24.

[15]Wesson R, Barlow M J, Matsuura M, et al. The timing and location of dust formation in the remnant of SN 1987A[J]. Monthly Notices of the Royal Astronomical Society, 2015, 446(2): 2089-2101.

[16]Compiègne M, Verstraete L, Jones A, et al. The global dust SED: tracing the nature and evolution of dust with DustEM[J]. Astronomy & Astrophysics, 2011, 525: A103-A117.

[17]Camps P, Baes M. SKIRT: an advanced dust radiative transfer code with a user-friendly architecture[J]. Astronomy and Computing, 2015, 9: 20-33.

[18]De Geyter G, Baes M, De Looze I, et al. Dust energy balance study of two edge-on spiral galaxies in the Herschel-ATLAS survey[J]. Monthly Notices of the Royal Astronomical Society, 2015, 451: 1728-1739.

[19]De Looze I, Fritz J, Baes M, et al. High-resolution, 3D radiative transfer modeling I. the grand-design spiral galaxy M 51[J]. Astronomy & Astrophysics, 2014, 571: A69-A91.

[20]Goosmann R W, Gaskell C M. Modeling optical and UV polarization of AGNs I. imprints of individual scattering regions[J]. Astronomy & Astrophysics, 2007, 465: 129-145.

[21]Marin F, Goosmann R W, Gaskell C M, et al. Modeling optical and UV polarization of AGNs II. polarization imaging and complex reprocessing[J]. Astronomy & Astrophysics, 2012, 548: A121-A146.

[22]Goosmann R W, Matt G. Spotting the misaligned outflows in NGC 1068 using X-ray polarimetry[J]. Monthly Notices of the Royal Astronomical Society, 2011, 415: 3119-3128.

[23]Matter A, Labadie L, Augereau J C, et al. Inner disk clearing around the Herbig Ae star HD 139614: Evidence for a planet-induced gap[J]. Astronomy & Astrophysic, 2016, 586: A11-A27.

Commentary on Three-dimensional Dust Monte-Carlo Radiative Transfer Models

Zhu Jiali1,2

(1. Key Laboratory of Space Astronomy and Technology, Beijing 100012, China, Email: jlzhu@nao.cas.cn;2. National Astronomical Observatories, Chinese Academy of Sciences, Beijing 100012, China)

Cosmic dust plays an important role in producing and processing interstellar radiation. Dust grains absorb starlight and convert it into infrared thermal radiation. This radiative transfer process is very complicated. To study properties of the stellar and interstellar medium by using multi-band observation data, numerical models need to be constructed to simulate the complicated radiative transfer process. Comparing the simulation results with the observation data, we can invert the input parameters of the numerical models, so as to recover the astrophysical and chemical properties associated with the input parameters. The three-dimensional (3D) distribution of dust requires the use of 3D radiative transfer models so as to simulate radiative transfer more realistically. Monte-Carlo method is applicable to any 3D dust density distribution. Instead of solving radiative transfer equations directly, it uses probability distribution function for random samplings. Monte-Carlo method propagates a large number of photons in a 3D grid to obtain the statistical results. To simulate different physical environments, we introduce six 3D dust Monte-Carlo radiative transfer models, including HO-CHUNK, Hyperion, 3D Mocassin, SKIRT, STOKES and RADMC-3D. These codes are open source codes based on the Monte-Carlo method. HO-CHUNK, 3D Mocassin, SKIRT and STOKES are used to simulate the protostar, photo-ionization environment, disk galaxies and active galactic nucleus respectively. Hyperion and RADMC-3D are general-purpose codes that can be used to simulate any 3D geometry. All of them calculate the dust local thermal equilibrium (LTE) emission, while HO-CHUNK and SKIRT can simultaneously calculate the non-LTE emission from very small grains and polycyclic aromatic hydrocarbon. All of the codes calculate the scattering of dust grains except HO-CHUNK. To improve the computing performance, all of them have been parallelized except STOKES. Basic input parameters of all the codes are: 3D distribution of the dust, physical and chemical properties of the dust, 3D distribution of the radiation sources, physical and chemical properties of the sources, and parameters control the simulation process and outputs. General simulation performance of these codes has been summarized into a table so as to help researchers find appropriate code which is more suitable for their research objectives.

Dust; Monte-Carlo technique; Numerical simulation; 3D radiative transfer

國家自然科學(xué)基金 (11503052)；中國科學(xué)院國家天文臺(tái)青年人才基金項(xiàng)目資助.

2016-03-02；

2016-04-06

朱佳麗，女，助理研究員. 研究方向：星際介質(zhì)，數(shù)值模擬. Email: jlzhu@nao.cas.cn

N32; P115.2+9

A

1672-7673(2016)04-0392-08

CN 53-1189/PISSN 1672-7673