楊 鑫， 王冠博， 李潤(rùn)東， 劉漢剛， 王 侃

（1.核物理與化學(xué)研究所，四川綿陽(yáng) 621900；2.清華大學(xué)工程物理系，北京 100084）

離子輸運(yùn)蒙特卡羅模擬的步長(zhǎng)選取

楊 鑫1，2， 王冠博1，2， 李潤(rùn)東1， 劉漢剛1， 王 侃2

（1.核物理與化學(xué)研究所，四川綿陽(yáng) 621900；2.清華大學(xué)工程物理系，北京 100084）

介紹帶電離子在物質(zhì)中輸運(yùn)模擬的詳細(xì)歷史法，開發(fā)模擬程序ITR（ions transport and reaction），同時(shí)實(shí)現(xiàn)抽取離子飛行步長(zhǎng)的基本方法和脈沖近似法，計(jì)算結(jié)果與CORTEO及TRIM程序的結(jié)果比較符合很好.討論兩種步長(zhǎng)選取方法，對(duì)脈沖近似法的最大步長(zhǎng)進(jìn)行限制，分析CORTEO程序與TRIM程序計(jì)算結(jié)果差別的主要原因. ITR同時(shí)具有TRIM和CORTEO的優(yōu)點(diǎn)，采用脈沖近似法減少需要計(jì)算的核碰撞次數(shù)，通過(guò)索引（Index）方法對(duì)散射角的計(jì)算進(jìn)行加速，計(jì)算效率大幅提高.

ITR；離子輸運(yùn)；詳細(xì)歷史法；蒙特卡羅

0 引言

在離子束分析技術(shù)、中子深度分析技術(shù)［1］、材料輻照損傷等研究中，帶電離子在物質(zhì)中輸運(yùn)過(guò)程的模擬是關(guān)鍵問(wèn)題.國(guó)際上，帶電離子“詳細(xì)歷史法”（detailed history）是研究熱門，不斷有新的程序出現(xiàn)，如SRIM［2－3］，SIMNRA［4］，CORTEO［5］等.其中，SRIM是由Ziegler等人開發(fā)的一套計(jì)算帶電離子在物質(zhì)中的阻止本領(lǐng)與輸運(yùn)過(guò)程的蒙卡模擬程序，采用ZBL（Ziegler-Biersack-Littmark）理論，TRIM是其中的離子輸運(yùn)計(jì)算模塊，是目前公認(rèn)的標(biāo)準(zhǔn)的詳細(xì)歷史計(jì)算程序.CORTEO程序是加拿大蒙特利爾大學(xué)Fran?ois Schiettekatte開發(fā)的一個(gè)離子輸運(yùn)快速模擬程序，理論基礎(chǔ)與TRIM程序相同，但其將散射角提前計(jì)算并以列表形式儲(chǔ)存，需要時(shí)通過(guò)索引（index）得到，而非每一次都直接計(jì)算，在對(duì)計(jì)算精度影響很小的情況下提高了計(jì)算效率.

本文根據(jù)詳細(xì)歷史法的原理，開發(fā)離子輸運(yùn)模擬程序ITR（ions transport and reaction），研究抽取離子飛行步長(zhǎng)的基本方法和脈沖近似法［2］，對(duì)兩者的聯(lián)系和區(qū)別進(jìn)行討論.同時(shí)，引入索引（Index）方法對(duì)核碰撞散射角的獲取進(jìn)行加速，計(jì)算效率大幅提高.

1 詳細(xì)歷史法

當(dāng)一個(gè)帶電離子在固體中穿行時(shí)，其能量損失可以分為兩部分，一部分是靶原子核做反沖運(yùn)動(dòng)的能量，用核阻止本領(lǐng)（－d E／d x）n表示；另一部分是激發(fā)或電離靶原子核外電子的能量，用電子阻止本領(lǐng)（－d E／d x）e表示.通?？梢越频卣J(rèn)為這兩種能量損失過(guò)程是獨(dú)立的.核阻止是造成離子運(yùn)動(dòng)軌跡偏轉(zhuǎn)的主要原因，電子阻止對(duì)離子能量的損失貢獻(xiàn)更大.

1.1 核阻止

1.1.1 勢(shì)函數(shù)

原子間受引力和斥力的耦合作用，斥力勢(shì)的描述很復(fù)雜，原子間斥力一部分是由于兩原子電子云重疊引起的，另一部分是帶正電的原子核之間的斥力.根據(jù)原子間距離r的不同，描述原子間作用力的勢(shì)函數(shù)有多種表達(dá)方式：當(dāng)r稍小于點(diǎn)陣內(nèi)原子平衡間距時(shí)，用Born-Mayer勢(shì)來(lái)表示.當(dāng)兩個(gè)核距離非常近，可以用庫(kù)侖勢(shì)表示.當(dāng)介于兩者之間時(shí)，描述起來(lái)要復(fù)雜得多.一般用屏蔽庫(kù)侖勢(shì)來(lái)描述.

其中，Z1，Z2為原子序數(shù)，e為單位電荷，Φ（r／a）稱為屏蔽函數(shù)，a是屏蔽半徑.

關(guān)于屏蔽函數(shù)與屏蔽半徑，有多種表示方式，其中以普適勢(shì)（universal）最為精確，它是1982年Ziegler和Biersack等人總結(jié)前人工作的基礎(chǔ)上提出理論模型，并擬合500多對(duì)離子和靶的組合進(jìn)行系數(shù)修正得來(lái)的，與實(shí)驗(yàn)偏差不超過(guò)5%，形式如下

其中x是約化半徑，x＝r／a，a0為Bohr半徑.

1.1.2 碰撞過(guò)程

離子與靶核發(fā)生碰撞的過(guò)程用經(jīng)典力學(xué)的二體碰撞理論描述.

在質(zhì)心系內(nèi)碰撞，取p為碰撞參數(shù).體系滿足能量守恒和角動(dòng)量守恒，可以得到在質(zhì)心系中，入射核的散射角為

在TRIM程序中，散射角的計(jì)算采用Magic Formula［3］，可以得到比較準(zhǔn)確的結(jié)果，但每一次碰撞都要對(duì)Magic Formula進(jìn)行計(jì)算，在很大程度上限制了程序的計(jì)算速度.而CORTEO程序則根據(jù)不同的能量和碰撞參數(shù)，先對(duì)散射角進(jìn)行計(jì)算存儲(chǔ)，需要時(shí)采用索引的方法得到，在對(duì)精度影響極低的情況下極大地提高了計(jì)算效率.

1.2 電子阻止

當(dāng)一個(gè)高能量的運(yùn)動(dòng)離子進(jìn)入固體內(nèi)后，絕大部分能量消耗在電子阻止上.電子阻止本領(lǐng)計(jì)算有眾多模型，如Bethe-Bloch模型、采用波恩近似與一級(jí)微擾的量子力學(xué)模型、Lindhard介電模型等.Ziegler以Linhard介電模型為基礎(chǔ)，推導(dǎo)了氫離子、He離子以及其他重離子在固體中的電子阻止本領(lǐng)計(jì)算方法［2］.沿離子飛行路程對(duì)阻止本領(lǐng)進(jìn)行積分就可以得到電子阻止造成的能量損失，通常，還需要考慮每一步的能量歧離.目前的電子阻止本領(lǐng)都是通過(guò)SRIM的SRModule模塊計(jì)算得到，而能量岐離計(jì)算有Bohr模型，Chu模型等［5］.

1.3 模擬步驟

詳細(xì)歷史法的計(jì)算流程如圖1所示.

2 確定飛行步長(zhǎng)

2.1 基本方法

基本的詳細(xì)歷史法中，認(rèn)為離子飛行的每個(gè)步長(zhǎng)僅與靶的原子間距有關(guān)，而與離子能量無(wú)關(guān)，CORTEO在默認(rèn)情況下采用該方法.

根據(jù)隨機(jī)態(tài)近似（RPA），在［L，L＋d L］的間隔里，找到一個(gè)碰撞核的概率P（L）d L服從Poisson分布

圖1 離子輸運(yùn)過(guò)程模擬流程圖Fig.1 Flowchart of ion transport simulation

2.2 脈沖近似法（impulse approximation，IA）

上面所述的詳細(xì)歷史法模擬中，每步的步長(zhǎng)僅與靶原子間距有關(guān)，平均為原子間距L0.因此，在離子運(yùn)動(dòng)的軌跡上，需要計(jì)算大量的核碰撞過(guò)程，尤其在離子能量相對(duì)較高時(shí)絕大多數(shù)的碰撞影響微乎其微，耗費(fèi)較多的計(jì)算時(shí)間.若能針對(duì)不同能量選擇不同步長(zhǎng)，使得每次碰撞的散射角分布基本相同，這樣高能時(shí)的自由飛行距離，即步長(zhǎng)L就會(huì)長(zhǎng)很多，從而大大減少模擬的步數(shù)，提高計(jì)算效率．對(duì)于同樣大小的散射角，當(dāng)離子入射能量較高時(shí)，碰撞參數(shù)p較小，而入射能量較低時(shí)，碰撞參數(shù)p較大.同時(shí)，觀察式（10），可以將碰撞參數(shù)寫為

其中b＝p／a，f（b）是一個(gè)與碰撞選用中心勢(shì)有關(guān)的函數(shù)，ε＝Eca／（Z1Z2e），為約化能量參數(shù).

為了使每次的碰撞模擬有研究?jī)r(jià)值，要求每次碰撞的能量傳遞不小于一個(gè)閾值Tmin，一般來(lái)說(shuō)Tmin取離位能或表面結(jié)合能.當(dāng)給定碰撞傳遞的能量，T＝Tmin時(shí)，由式（5）可得

3 程序開發(fā)與校驗(yàn)

按照上面的原理，開發(fā)了帶電離子詳細(xì)歷史法模擬程序ITR［6］，在散射角的計(jì)算中，包含了CORTEO的data數(shù)據(jù)包，采用索引技術(shù)以提高計(jì)算速度.對(duì)步長(zhǎng)的選取，提供了基本的選取法與脈沖近似法，可根據(jù)需求選擇.

3.1 ITR程序校驗(yàn)

分別對(duì)初始能量為500 keV和1.472 1 MeV的4He離子沿x坐標(biāo)方向垂直射入Si靶進(jìn)行模擬，與TRIM和CORTEO的計(jì)算結(jié)果進(jìn)行比較.

500 keV的4He離子垂直射入1μm厚的Si靶后的比較見圖2、圖3、表1.其中，l，m，n分別為飛行方向與x，y，z方向的夾角余弦值，由于對(duì)稱性，m，n按兩者的和作圖.

1.472 1 Mev的4He離子垂直射入5μm厚的Si靶后的比較見圖4、圖5、表2.

圖2 500 keV的4He離子垂直入射并穿透1μm厚的Si時(shí)的能譜Fig.2 Energy spectrum of 500 keV vertically implanting4He ions transmiting 1μm Si

圖3 500 keV的4He離子垂直入射并穿透1μm厚的Si時(shí)的飛行方向Fig.3 Flight direction of 500 keV vertically implanting4He ions transmiting 1μm Si

可以看到，以上兩種情況下，采用脈沖近似法選取步長(zhǎng)的ITR（IA）與基本的步長(zhǎng)選取法ITR（basic）的結(jié)果分別與TRIM和CORTEO符合得很好，驗(yàn)證了程序的正確性.脈沖近似法需要計(jì)算的碰撞次數(shù)遠(yuǎn)小于基本的步長(zhǎng)選取方法，可以提高計(jì)算效率.從耗時(shí)情況看，TRIM每次碰撞的散射角都直接采用Magic Formula計(jì)算，耗費(fèi)大量時(shí)間，計(jì)算效率欠佳，CORTEO采用索引（Index）方法對(duì)散射角的計(jì)算進(jìn)行加速后，計(jì)算效率高于TRIM，而ITR（IA）由于同時(shí)采用了脈沖近似和Index技術(shù)，綜合了兩者的優(yōu)點(diǎn)，計(jì)算速度更快.

3.2 步長(zhǎng)選取

上面的比較中，CORTEO與TRIM的計(jì)算結(jié)果符合不是很好，兩個(gè)程序最大的區(qū)別在于步長(zhǎng)選取方法的不同.從ITR中對(duì)兩種方法的選取結(jié)果看出，由于脈沖近似法得到的飛行步長(zhǎng)較大，計(jì)算的碰撞次數(shù)少了很多，因此其計(jì)算效率增加了.下面，對(duì)ITR采用脈沖近似法時(shí)的最大步長(zhǎng)進(jìn)行限制，討論其結(jié)果與基本選取方法的差異.步長(zhǎng)限制即限定最大的單次飛行步長(zhǎng)Lmax＝f×d，其中d為靶厚，f為限制系數(shù)．若脈沖近似法的步長(zhǎng)大于 Lmax，則令該步長(zhǎng)為L(zhǎng)max.

圖4 1.472 1 MeV的4He離子垂直入射并穿透5μm厚的Si時(shí)的能譜Fig.4 Energy spectrum of1.472 1 MeV vertically implanting4He ions transmiting 5μm Si

圖5 1.472 1 MeV的4He離子垂直入射并穿透5μm厚的Si時(shí)的飛行方向Fig.5 Flight direction of 1.472 1 MeV vertically implanting4He ions transmiting 5μm Si

表1 500 keV的4He離子垂直入射并穿透1μm厚的Si時(shí)程序的比較Table 1 Calculation of 500 keV vertically implanting4He ions transm iting 1μm Si

表2 1.472 1 MeV的4He離子垂直入射并穿透5μm厚的Si時(shí)程序的比較Table 2 Calculation of 1.472 1 MeV vertically im p lanting4He ions transm iting 5μm Si

分別對(duì)500 keV和1.472 1 MeV的4He離子垂直射入Si靶進(jìn)行研究.

500 keV時(shí)，f分別取0.01，0.001，0.000 4的結(jié)果如圖6和表3所示.

1.472 1 MeV時(shí)，f分別取0.001，0.000 4，0.000 1的結(jié)果如圖7和表4所示.

隨著f的減小，單個(gè)離子的碰撞次數(shù)越多，耗費(fèi)的計(jì)算時(shí)間越長(zhǎng)，能譜形狀由與TRIM接近逐漸轉(zhuǎn)向與CORTEO接近，特別是當(dāng)4He離子能量為1.472 1 MeV時(shí)，此效果更為明顯，說(shuō)明步長(zhǎng)選取方法的不同是CORTEO與TRIM間計(jì)算產(chǎn)生差別的主要原因.

圖6 500 keV的4He離子垂直入射并穿透1μm厚的Si時(shí)f取不同值的能譜Fig.6 Energy spectrum of500 keV vertically implanting4He ions transmiting 1μm Siwith different f

圖7 1.472 1 MeV的4He離子垂直入射并穿透5μm厚的Si時(shí)f取不同值的能譜Fig.7 Energy spectrum of 1.472 1 MeV vertically implanting4He ions transmiting 5μm Siwith different f

表3 500 keV的4He離子垂直入射并穿透1μm厚的Si時(shí)f取不同值程序的比較Table 3 Calculation of 500 keV vertically im p lanting4He ions transm iting 1μm Siwith different f

表4 1.472 1 MeV的4He離子垂直入射并穿透5μm厚的Si時(shí)f取不同值程序的比較Table 4 Calculation of 1.472 1 MeV vertically implanting4He ions transm iting 5μm Siwith different f

由前面的分析知道，步長(zhǎng)L和碰撞參數(shù)p之間是相互耦合的，而散射角反映了L路程中所有碰撞造成飛行方向偏移的積分效應(yīng)，L越大，p越小，散射角越大．因此，L的選取有一定的隨意性，可以根據(jù)不同的需求選取L的大小，最終會(huì)通過(guò)p的變化進(jìn)而影響對(duì)散射角的計(jì)算．步長(zhǎng)選取的基本方法和脈沖近似法即為兩種不同選取標(biāo)準(zhǔn)，基本方法以靶原子間距為平均值對(duì)步長(zhǎng)進(jìn)行抽樣，雖然看上去與真實(shí)情況更為接近，但當(dāng)離子能量較低時(shí)，散射角受核碰撞的影響是極其敏感的，若抽取到的L值較大，會(huì)使得一些重要的核碰撞被忽略.而脈沖近似法考慮離子能量的影響，既在離子能量較高時(shí)忽略影響不大的核碰撞，又在離子能量較低時(shí)充分考慮了核碰撞對(duì)散射角造成的影響，這種方式在減少核碰撞次數(shù)的同時(shí)又兼顧了重要性更高的核碰撞，因此，其得到的模擬結(jié)果更好.TRIM仍然是目前公認(rèn)的詳細(xì)歷史模擬法的標(biāo)準(zhǔn)程序.

4 結(jié)論與展望

介紹了帶電離子輸運(yùn)的詳細(xì)歷史法的物理原理及模擬步驟，編寫了模擬程序ITR并驗(yàn)證其正確性.對(duì)步長(zhǎng)選取方法的研究，找到了CORTEO與TRIM間計(jì)算產(chǎn)生差別的主要原因，分析認(rèn)為脈沖近似法仍然是目前最好的步長(zhǎng)選取方法.ITR程序結(jié)合TRIM和CORTEO的優(yōu)點(diǎn)，同時(shí)采用了選取步長(zhǎng)的脈沖近似法及獲取散射角的索引技術(shù)，提高了計(jì)算效率.這使得諸如在中子深度分析中，完全采用蒙特卡羅方法進(jìn)行模擬成為可能，以往都是在模擬多種情況后，通過(guò)擬合數(shù)據(jù)的方式得到.今后將結(jié)合中子深度分析技術(shù)的物理特性及探測(cè)器的優(yōu)化模擬，編寫專用的模擬程序MCNDP（Monte Carlo code for neutron depth profilling）.

［1］Ziegler JF，Cole GW，Baglin JE E.Technique for determining concentration profiles of boron impurities in substrates［J］.Journal of Applied Physics，1972，43：3809－3815.

［2］Ziegler JF，Biersack JP，Ziegler MD.SRIM—The stopping and range of ions inmatters［M］.New York：Pergamon Press，2008.

［3］Ziegler JF，Ziegler MD，Biersack J P.SRIM—The stopping and range of ions in matter［J］.Nuclear Instruments and Methods in Physics Research Section B，2010，268（11－12）：1818－1828.

［4］Maire M.Ion beam analysis of rough thin film［J］.Nuclear Instruments and Methods in Physics Research Section B，2002，194（2）：177－186.

［5］Fran?ois Schiettekatte.Fast Monte Carlo for ion beam analysis simulations［J］.Nuclear Instruments and Methods in Physics Research B，2008，266：1880－1885.

［6］Wang G，Yang X，Liu H，et al.“Detail history”Monte Carlo method for ion simulation［J］.Chinese Journal of Computational Physics，2013，30（2）：303－308.

Flight Path Selection in Monte Carlo Code ITR for Ion Transport

YANG Xin1，2，WANG Guanbo1，2，LIRundong1，LIU Hangang1，WANG Kan2
（1.Institution ofNuclear Physics and Chemistry，Mianyang，Sichuan 621900，China；2.Department ofEngineering Physics，Tsinghua University，Beijing 100084，China）

Strategy of detailed history method is summarized，depending on which code ITR（ions transport and reaction）was developed.Both basicmethod and impulse approximationmethod to select flight path of ions are complemented respectively.Numerical results agree with CORTEO and TRIMvery well.Relation between two methods is discussed.We found reason of CORTEO results differ from TRIM's by limitingmaximum of flight length derived by impluse approximation method.ITR has advantages of TRIMand CORTEO.Impluse approximation method reduces nuclear collisions simulated.Index technique is used to calculate scattering angle，which significantly improves efficiency of ITR.

ITR；ions transport；detailed historymethod；monte carlo

date： 2012－11－27；Revised date： 2013－05－17

TL99

A

1001-246X（2014）04-0417-07

2012－11－27；

2013－05－17

中國(guó)工程物理研究院科學(xué)技術(shù)發(fā)展基金（2010B0103009）及國(guó)家自然科學(xué)基金（11175162）資助項(xiàng)目

楊鑫（1988－），男，貴州，博士生，主要從事反應(yīng)堆物理及應(yīng)用研究，E-mail：yangx05＠126.com