離子輸運蒙特卡羅模擬的步長選取

楊 鑫， 王冠博， 李潤東， 劉漢剛， 王 侃

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

離子輸運蒙特卡羅模擬的步長選取

楊 鑫1，2， 王冠博1，2， 李潤東1， 劉漢剛1， 王 侃2

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

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

ITR；離子輸運；詳細歷史法；蒙特卡羅

0 引言

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

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

1 詳細歷史法

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

1.1 核阻止

1.1.1 勢函數(shù)

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

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

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

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

1.1.2 碰撞過程

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

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

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

1.2 電子阻止

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

1.3 模擬步驟

詳細歷史法的計算流程如圖1所示.

2 確定飛行步長

2.1 基本方法

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

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

圖1 離子輸運過程模擬流程圖Fig.1 Flowchart of ion transport simulation

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

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

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

為了使每次的碰撞模擬有研究價值，要求每次碰撞的能量傳遞不小于一個閾值Tmin，一般來說Tmin取離位能或表面結(jié)合能.當給定碰撞傳遞的能量，T＝Tmin時，由式（5）可得

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

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

3.1 ITR程序校驗

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

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

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

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

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

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

3.2 步長選取

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

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

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

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

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

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

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

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

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

圖6 500 keV的4He離子垂直入射并穿透1μm厚的Si時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時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時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時f取不同值程序的比較Table 4 Calculation of 1.472 1 MeV vertically implanting4He ions transm iting 5μm Siwith different f

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

4 結(jié)論與展望

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

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

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