银河系除了会发光的恒星之外，星际空间被自身不发光的星际介质占据。星际介质由气体和尘埃组成，尘埃吸收星光能量，并将其转化为红外波段的热辐射，同时通过光电作用加热气体。通过多波段的观测数据研究恒星等辐射天体和星际介质，需要了解辐射在通过介质时的转移情况，以获得辐射源和介质的几何、物理和化学性质。由于这个辐射转移过程相当复杂，一般通过构建数值模型来模拟和研究，即在一定的模型输入参数下，将数值模型的模拟结果与实际观测数据相比较，再通过拟合观测数据反演数值模型的输入参数，从而掌握与输入参数相关的天体物理化学性质。

尘埃的分布为三维不均匀分布，为了更真实地模拟尘埃的辐射转移过程，需采用三维辐射转移模型。然而，尘埃三维辐射转移问题有多方面的难点^[1]。首先它是一个非本地非线性的六维问题，包括了辐射场的坐标（三维）、传播方向（二维）和辐射波长这6个参数。其次，尘埃三维分布的复杂度越高，即模型所模拟的三维空间的几何分辨率越高，模型计算所需的存储空间越大，计算难度也越大。

三维辐射转移问题通常使用蒙特卡洛方法解决。虽然也有其它的技术和方法被开发来计算辐射转移，但基本是对辐射转移方程直接计算数值解。这种解方程的数值方法能有效地解决一维问题，但在二维和三维问题上变得十分复杂^[2]，更无法有效解决三维辐射转移问题，或者存在显著的缺陷^[1]（如忽略散射和尘埃热辐射等）。蒙特卡洛方法非常适用于任意的尘埃三维密度分布，因为它不是直接计算辐射转移方程，而是利用对概率密度函数(Probability Distribution Function,PDF)的随机抽样，达到在三维空间中传播大量光子从而获得统计结果的目的^[2]。

本文第1节介绍了基于蒙特卡洛方法的尘埃三维辐射转移的基本模拟过程。针对不同物理环境的模拟，在第2节分析述评了6个开源代码，并概括了它们的模拟性能。

1 用蒙特卡洛方法模拟尘埃的三维辐射转移

在蒙特卡洛算法中，辐射场不断地发出大量的光包，每一个光包由一定数量的光子组成。光包沿着一定的路径穿过尘埃介质，直到被吸收或离开介质。在此过程中的每一阶段，光包的下一段路径的特征参数都是基于适当的概率密度函数随机抽样确定的，具体的抽样方法可参见文^[3]。通过发射足够数量的光包，并对所有光包的路径进行统计分析而得到辐射场最终的分布情况。可见，蒙特卡洛方法是通过模拟辐射转移过程代替解辐射转移方程。

1.1 基本模拟过程

在最简化的情况下，辐射源发射某个频率的光子，光子在具有一定密度分布的三维网格（即离散化的三维空间^①）中传播，每一个光子在被散射或吸收（随即被再发射）之前都会经过一定的光深，这些光子将继续在网格中传播直到它们离开网格^[2]。离开网格的光子可用来生成光谱能量分布(Spectral Energy Distributions,SEDs)和图像。通过发射大量的光子重复以上过程，将提高模拟所得的光谱能量分布和图像的信噪比。

①基于一定的坐标系，将要模拟的三维空间离散化，形成三维网格。将模型参数（如尘埃密度和辐射强度等）随网格离散化，且参数值在同一个单元格内保持不变。

1.1.1 光包的发射和传播

为了提高模拟效率，辐射源实际发射的是光包而不是单个光子。在有多个源的情况下，光包来自哪个源随机抽样决定。每个源发射的光包数目可以和源的光度成正比，也可以假定每个源都发射同样数目的光包，但不同源发射的光包的权重不同^[2]。所发射光包的传播方向、频率随机抽样决定。

光包发射后开始传播过程。对光包将要经过的光深进行抽样，并和光包当前路径上的实际光深比较，当抽样光深大于实际光深时，光包将离开网格；反之，光包将和尘埃相互作用，具体是被散射或被吸收则由尘埃的发射率决定。被吸收后再发射的光包频率通过随机抽样获得。光包被散射或被吸收再发射后通过随机抽样得到新的传播方向和新的光深值，并重复上述传播过程直到光包离开网格。

蒙特卡洛方法跟踪大量光包的传播轨迹，光包在其存在的每个阶段都有一系列的量来标注其属性，比如坐标、传播方向和波长，根据研究目标的不同，还可以额外记录诸如多普勒位移和偏振等信息^[4]。当光包和三维网格中的介质发生相互作用，光包的上述属性将发生改变，并通过对适当的概率密度函数进行随机抽样更新光包的属性。当光包离开网格，它的最终属性被记录下来。通过这种方式记录大量的光包信息，将得到被模拟系统的总体统计特性。

1.2 模拟结果的处理

上述被模拟系统的总体统计特性即为模拟结果。在不同模型参数设置下进行多次模拟，并将模拟结果和观测数据相比较可获知观测对象的性质，一般步骤如下^[1]：

(1) 将尘埃密度分布、辐射源性质、尘埃光学性质等进行参数化；

(2) 确定模型网格的设置，将步骤(1)中的参数在网格中离散化，使用辐射转移代码计算得到光谱能量分布或者图像等模拟结果；

(3) 比较模拟结果和观测数据；

(4) 评价步骤(3)中两者的差距，改变模型的参数和（或）模型的假设，重复步骤(2)、(3)、(4)，直到模拟和观测的差距最小化；

(5) 找到最优的模型参数设置。

其中，合理评价观测和模拟误差对步骤(3)中的比较十分重要，因为它直接影响步骤(4)中对模型参数和假设的选择，最终影响最优模型设置的获得。

2 尘埃三维辐射转移模型代码及其特点分析

计算天体物理学家对一维辐射转移模型代码的开发使用已有近30年的经验，而三维辐射转移模型却是近些年随着计算机运算能力的提升才发展起来的。到20世纪90年代中期，还没有辐射转移模型代码可以实现一般三维几何体中星际尘埃对辐射的吸收、散射和尘埃热辐射的计算^[1]。可用的代码不是只适用于一维或二维几何体，就是不能计算完整的三维辐射转移问题，如忽略散射或尘埃热辐射。 在过去的20年间，计算机运算能力的显著增强和更高效算法的开发，促使了三维尘埃辐射转移模型代码的产生，这些代码，具备了计算三维几何体中尘埃对辐射的吸收和各向异性散射的能力。这些代码的应用范围十分广泛，可适用于原恒星、活动星系核(Active Galactic Nucleus,AGN)、分子云和星系等。

这些代码一般用Fortran、C或C++语言编写，同时为了提高运算效率，提供了并行化运算功能。代码的基本输入参数包括：尘埃三维分布情况、尘埃的物理化学性质、辐射源坐标和物理化学性质以及控制代码运行和输出的一些参数。

三维尘埃辐射转移模型中，尘埃颗粒的最简单假设是单一种类的处于局部热平衡的尘埃颗粒。这种情况下，尘埃辐射为灰体辐射，每个单元格的尘埃的平衡温度由吸收和辐射平衡决定。第2种更符合现实的做法是，考虑处于热平衡状态的不同种类和尺寸的尘埃，即将不同种类和尺寸的尘埃颗粒分开处理，放入不同的 “容器”(bin)内，每一个容器都处于局部热平衡，尘埃的总辐射是一系列灰体辐射的总和，其中每个容器的温度由吸收和辐射平衡决定。第3种是最接近实际情况的做法，即考虑来自小尺寸尘埃和多环芳香烃(Polycyclic Aromatic Hydrocarbons,PAH)的非局部热平衡辐射^[5]。非局部热平衡辐射和局部热平衡辐射的主要区别是尘埃温度不是用单个的平衡温度表征，而是用温度分布表征。在实践中，要在辐射转移代码中包含尘埃非局部热平衡辐射的计算是一项艰难的工程。下文介绍的代码都可以计算尘埃的热平衡辐射，其中个别代码还可同时计算非局部热平衡辐射。

到目前为止，已公开发表的尘埃三维辐射转移模型代码有近30个^[1]，其中90%为基于蒙特卡洛方法开发的代码。本节介绍其中6个开源的蒙特卡洛辐射转移模型代码（见表 1），这些代码基于不同的开发目的能针对性地模拟不同的科学问题。在模拟研究中，合理选择模型代码有利于得到更真实的天体信息。

2.1 HO-CHUNK的模拟性能和应用范围

原恒星一般由年轻恒星和包围它的盘、下落的包层以及双极外流组成。HO-CHUNK可模拟原恒星从Class 0到Class Ⅲ阶段的一系列演化过程中的光谱能量分布、颜色、偏振性质和图像^[6]。

高分辨率和多波段的观测揭示了原恒星盘存在复杂的三维结构，例如弯曲、裂缝、旋涡状的尘埃密度分布等。HO-CHUNK 可模拟盘的裂缝和旋涡结构、内盘的弯曲变形、中央恒星上的吸积热点、尘埃密度分布的不规则团块结构、包层的裂缝，以及双极外流的多重空腔壁^[7]等三维结构，使模型更符合实际情况。

HO-CHUNK 采用球极坐标系，网格结构的设置允许一个单元格内有多种尘埃成分。辐射源包括3类，分别为恒星辐射、盘吸积产生的辐射和星际辐射场。尘埃盘分为内盘和外盘，可分别具有不同的几何结构和尘埃性质，尘埃密度分布可为团块结构。模型计算来自大尺寸尘埃的热平衡辐射，以及来自小尺寸尘埃和多环芳香烃的非局部热平衡辐射，但不计算尘埃的散射。该代码预定义了多种尘埃模型，使用者也可以构建自己的尘埃模型。

该代码为信息传递接口(Message Passing Interface,MPI)并行运算代码，可在一台计算机的多个处理器或者多台联网的计算机上运行。针对模拟结果，HO-CHUNK提供IDL程序对结果进行后续处理与分析。

2.1.1 应用和成果

不同演化阶段的原恒星具有不同的尘埃辐射特征。通过使用HO-CHUNK模拟发现，Class 0阶段原恒星中波长大于100 μm的红外辐射来自其包层，而Class I阶段原恒星的这部分红外辐射则主要来自尘埃盘^[6]，这和前人的观测及模拟结果一致^[8]。 和低质量的原恒星相比较，高质量的原恒星的数量较少且演化时标较短，导致难以区分实际观测数据所对应的演化阶段，对进一步分析和挖掘观测数据造成障碍。使用HO-CHUNK模拟不同演化阶段的高质量原恒星的尘埃辐射特征和形态特征，并和观测数据相比较，从而确定实际观测数据所对应的演化阶段，同时也为鉴别观测对象是否为原恒星提供有力的线索^[7]。使用上述方法，文^[9]从多波段巡天数据中鉴定出了约1 000个原恒星。

2.2 Hyperion的模拟性能和应用范围

Hyperion可模拟任意三维空间的尘埃连续辐射^[2]。该代码基本是通用的，能适用于多种科学目标，可计算尘埃温度、光谱能量分布和多波段图像。

该代码由两部分组成：一部分进行辐射转移的计算，由Fortran执行，这部分代码和具体科学问题无关；第二部分包括了一些Python库，用于便捷地构建输入文件，并包含一些预定义的常用密度分布结构，同时提供了模拟结果的后续处理工具。

Hyperion采用三维直角坐标系、球极坐标系和圆柱坐标系。网格中可分布任意位置和数量的点源、一定半径的球状源、弥散源、星际辐射场和平面源，尘埃的分布也为任意三维分布。Hyperion计算来自大尺寸尘埃的热辐射，并同时计算尘埃散射和偏振辐射。

Hyperion是信息传递接口并行代码，并提供Python库便于使用者对模拟结果进行后续处理与分析。

2.2.1 应用和成果

Hyperion基于通用目的开发，但目前主要应用于原恒星的研究。和HO-CHUNK对原恒星复杂三维结构的预定义相比，Hyperion对原恒星的几何结构的设置更为灵活，被用于原恒星候选体的鉴定^[10]和原恒星盘几何结构的确定^[11]。

Class 0阶段的原恒星的光谱能量分布明显偏向红端且在亚毫米波段的流量较大。文^[10]将观测所得原恒星候选体光谱能量分布和不同模型参数设置下Hyperion所预言的光谱能量分布相比较，获知原恒星候选体的物理化学特征，从而确定候选体极可能为原恒星。文^[11]则通过比较模型模拟结果和实际观测数据来搜寻Class 0类型原恒星的盘结构特征，发现Class 0阶段的原恒星已可以具备盘结构。

2.3 3D Mocassin的模拟性能和应用范围

光致电离环境的研究是很多天体物理学问题的基础。3D Mocassin用于模拟包含尘埃的光致电离星云^[12]，能计算任意几何结构下被辐射场强电离的低密度气态环境，辐射场波长范围可从射电波段到伽玛射线波段^[13]。

3D Mocassin采用三维直角坐标系，能处理任意的三维几何结构和密度分布。辐射源可同时包含恒星辐射和弥散辐射场，可处理任意数量和位置的点源以及位于 X-Z 平面的平面源。3D Mocassin计算大尺寸尘埃的热辐射和散射。该代码在计算星云的发射线光谱、电子温度和密度结构的同时，还可以输出星云的光谱能量分布。

3D Mocassin是信息传递接口并行代码。相比其余5个代码，3D Mocassin的特色是能同时处理气体和尘埃的辐射转移，更符合星云中气体和尘埃相互混合的真实情况。尘埃和气体之间的相互作用包含了相互碰撞和光电效应，每一个单元格内的物理状况通过解热平衡和电离平衡方程获得。根据不同的模拟目的，也可对尘埃和气体分别计算，即构建纯尘埃或纯气体模型。

2.3.1 应用和成果

3D Mocassin主要用于富含尘埃和气体的行星状星云和超新星的研究中。通过拟合光谱信息可以获知星云的温度、尺度和光度等信息^[14]。

在超新星的研究中 ，3D Mocassin被用于计算不同尘埃质量和分布情况下超新星的辐射转移情况。通过拟合超新星随时间演化的尘埃辐射，可以获知超新星爆发后尘埃形成的时间和位置。文^[15]在分析拟合了SN 1987A从爆发后64天到24年间的光谱能量分布数据后，发现大部分的尘埃形成于超新星爆发后的晚期。

2.4 SKIRT的模拟性能和应用范围

SKIRT通过构建细致的辐射转移模型来研究盘星系中尘埃的吸收和散射对星系动力学参数的影响^[4]。模拟尘埃星系的动力学参数时，模型需要记录光包的多普勒位移。

SKIRT采用三维直角坐标系、一维球坐标系和二维圆柱坐标系，并预定义了多种几何结构用于设置辐射源和尘埃的空间密度分布，包括解析的球对称、轴对称分布，以及非解析的不对称分布。SKIRT可同时计算尘埃的热平衡辐射和非局部热平衡辐射。它结合了DustEM^[16]这个开源的尘埃辐射的数值计算工具，使其可以快速准确地计算被瞬时加热的极小尘埃颗粒和多环芳香烃的非局部热平衡辐射^[5]。模型输出光谱能量分布、星系动力学参数和图像。

SKIRT被设计成用户友好架构，即通过和用户交互的方式构建模型模拟所需的配置文件，然后通过命令行方式根据配置文件运行模拟程序^[17]。SKIRT是信息传递接口并行代码。

2.4.1 应用和成果

侧向旋涡星系的尘埃分布往往呈现特有的窄带形态，是研究尘埃在星系中的径向和垂直方向分布情况的理想对象。通过使用SKIRT模拟NGC 5166和IC 4225这两个侧向旋涡星系在远红外波段尘埃辐射的光谱能量分布，并和观测所得光谱能量分布相比较，文^[18]发现年轻星族对加热星系尘埃起了至关重要的作用，缺少年轻星族的星系模型所预言的尘埃辐射要低于实际观测值的3倍。

SKIRT同样适用于正向旋涡星系的模拟。文^[19]结合红外、紫外和光学的多波段数据，使用SKIRT模拟了M 51中尘埃的远红外和亚毫米波段的辐射。通过比较模型模拟结果和观测数据，分析了年轻星族和年老星族对加热尘埃所作的贡献，其中年老星族提供了加热尘埃所需总能量的37%。

2.5 STOKES的模拟性能和应用范围

由于光子被自由电子和尘埃颗粒散射，活动星系核的光在很大的波长范围内是偏振的，这对辐射和散射区域的几何结构形成了重要的约束。STOKES计算活动星系核统一模型所描述的辐射转移，主要用于模拟活动星系核中尘埃造成的偏振^[20]，从而研究不同几何形状和光深的尘埃环对辐射偏振的影响。

STOKES采用三维直角坐标系，模拟区域分为源区和散射区两部分。源区可包括圆柱形的连续辐射源、圆柱形的宽线源和双圆锥形的窄线源；散射区的尘埃分布结构包括圆柱形、圆环、双圆锥、喇叭形和球形。STOKES计算大尺寸尘埃的热辐射和散射，输出的主要模拟结果为光谱能量分布、辐射流量和偏振信息。模型还可记录每一个光子从发射到离开网格所经历的时间，从而模拟和分析时变的偏振。

经过多年开发，STOKES已被扩展为可模拟宽线区发射线的偏振对该区域速度的依赖情况^[21]。宽线区的运动主要是旋转和明显的湍流，很可能还有内流。STOKES不能进行并行运算。

2.5.1 应用和成果

使用STOKES进行模拟研究发现，尘埃环的形状对造成光线偏振的程度轻重有显著的影响^[20]。通常认为，活动星系核最靠近中心的部分是由多种产生辐射和散射的介质组成的，这些介质在辐射转移过程中耦合。光谱偏振观测为这些区域的几何结构提供了重要的线索。通过构建模型模拟各成分间的辐射转移耦合获得活动星系核的偏振特性，并将模型模拟结果和偏振观测结果相比较，从而约束活动星系核的几何结构。

文^[22]使用STOKES模拟了Ⅱ型塞弗特星系NGC 1068的X波段的光谱偏振，分析了偏振特性随辐射转移发生区域和观测角度的变化情况，预言了X波段的偏振观测将对活动星系核内部几何结构的确定起到重要作用。

2.6 RADMC-3D的模拟性能和应用范围

RADMC-3D定位为多用途辐射转移工具，是基于通用目的开发的代码。RADMC-3D可用于模拟任意三维空间尘埃的辐射转移，比如可用于模拟原行星盘、活动星系核、分子云和行星大气等的辐射转移^②。

②http://www.ita.uni-heidelberg.de/~dullemond/software/radmc-3d/

针对不同的模拟对象，RADMC-3D预定义了一系列的模型设置作为样本，方便用户在样本基础上改进模型参数。值得一提的是，RADMC-3D的使用文档是本文介绍的代码中最为完善的，并且提供专门的论坛便于使用者和开发者之间的交流^③。

③http://radmc3d.ita.uniheidelberg.de/phpbb/index.php

根据到目前为止的开发情况，RADMC-3D采用三维直角坐标系和球极坐标系，辐射源包括任意数量和位置的点源、一定半径的球状源和星际辐射场。该代码计算大尺寸尘埃热辐射和散射，输出模拟区域的光谱能量分布和图像，并提供IDL程序便于使用者对模拟结果进行后续处理与分析。RADMC-3D是OpenMP并行代码。

2.6.1 应用和成果

作为一个新近开发的代码，目前已发表的RADMC-3D的具体应用较少。具有代表性的应用是在赫比格Ae型星HD 139614的盘尘埃分布的研究，文^[23]使用RADMC-3D模拟该恒星尘埃盘的辐射图像和光谱能量分布。通过拟合观测数据，确定了尘埃盘中裂缝的存在，并得到了盘中尘埃的物理化学性质。

3 总 结

本文介绍了用蒙特卡洛方法模拟尘埃三维辐射转移的基本过程，并列举了6个基于蒙特卡洛方法的尘埃三维辐射转移模型的开源代码。HO-CHUNK、3D Mocassin、SKIRT和STOKES这4个代码可分别用于模拟原恒星、光致电离环境、盘星系和活动星系核，Hyperion和RADMC-3D则为通用型代码，可模拟任意三维空间尘埃的辐射转移。这些开源代码的基本模拟性能特点已整理成表 2，以便于相关研究人员根据研究目标选择使用代码。