Ar射頻放電特性隨時(shí)間演化的PIC/MCC模擬

張浚源 陳 峰 孫偉中 呂曉丹 賀平逆 茍富均,3

1（四川大學(xué)原子核科學(xué)技術(shù)研究所，輻射物理及技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室 成都 610064）

2（貴州大學(xué)PSI研究所MEMS課題組 貴陽(yáng) 550025）3（荷蘭皇家科學(xué)院等離子體所 荷蘭 2300）

射頻濺射制膜裝置中產(chǎn)生的非平衡態(tài)等離子體中的電子平均速度遠(yuǎn)大于離子，到達(dá)射頻極板表面的電子數(shù)遠(yuǎn)多于離子。這些電子一部分與表面發(fā)生復(fù)合，大部分則聚集在極板表面，導(dǎo)致極板附近相對(duì)于等離子體區(qū)呈負(fù)電勢(shì)，此負(fù)電勢(shì)阻擋電子向極板處運(yùn)動(dòng)，并吸引正離子，導(dǎo)致極板附近區(qū)域離子密度與電子密度分離而形成鞘層[1–4]。等離子體中的離子經(jīng)鞘層電場(chǎng)加速后轟擊靶材使其發(fā)生濺射，濺射出的中性靶原子或基團(tuán)大部分以擴(kuò)散的方式沉積在基片上，但有些在輸運(yùn)過(guò)程中與等離子體碰撞，形成離子而轟擊基片。這表明射頻濺射制膜工藝是個(gè)非常復(fù)雜的過(guò)程，它包括等離子體放電行為、濺射產(chǎn)物的輸運(yùn)以及等離子體與基片表面相互作用[5,6]。為理解薄膜制備過(guò)程中的各種微觀機(jī)理行為，有必要研究射頻濺射工藝中等離子體放電行為。然而，常用的等離子體診斷手段僅能探測(cè)等離子體區(qū)域和預(yù)鞘層區(qū)域[4]，用實(shí)驗(yàn)手段了解射頻放電行為有一定的局限性。

計(jì)算機(jī)模擬方法可有助于了解等離子體放電過(guò)程中的動(dòng)力學(xué)行為。邱清泉[7]等模擬了平面直流磁控濺射放電等離子體，發(fā)現(xiàn)隨著氣壓的升高，鞘層厚度基本不變，鞘層電勢(shì)降會(huì)增大。Lee J K[8,9]等用一維等離子體粒子模擬法(particle-in-cell/Monte Carlo collision, PIC/MCC), 研究了容性耦合等離子體(Capacitively Coupled Plasma, CCP)源中的氬氣放電行為，發(fā)現(xiàn)Ar+離子能量分布中的高能峰明顯隨著射頻頻率增大。Nanbu等[10]發(fā)現(xiàn)，相對(duì)于磁場(chǎng)而言，氣壓對(duì)電場(chǎng)的影響較小。

本文用PIC/MCC法利用一維模型模擬CCP源氬氣放電過(guò)程中等離子體動(dòng)力學(xué)行為。主要考察射頻極板處和腔體中心處粒子流量、電流以及鞘層厚度隨時(shí)間的演化關(guān)系。從理論上探討射頻等離子體放電過(guò)程中的物理和化學(xué)機(jī)制。

1 PIC/MCC模擬方法

PIC/MCC通過(guò)求解牛頓運(yùn)動(dòng)方程追蹤大量超粒子(106–109個(gè)電子或離子)在外加電場(chǎng)及自洽電場(chǎng)中的運(yùn)動(dòng)過(guò)程，同時(shí)應(yīng)用蒙特卡洛碰撞模型得到超粒子碰撞之后的運(yùn)動(dòng)狀態(tài)，得到模擬系統(tǒng)中每個(gè)超粒子在求解時(shí)間內(nèi)的運(yùn)動(dòng)軌跡[11-12]?？捎媒y(tǒng)計(jì)平均法求解等離子體密度、溫度、電勢(shì)、電場(chǎng)和粒子通量等放電特性。通過(guò)跟蹤不同時(shí)刻等離子體特性參量的變化，可以實(shí)現(xiàn)分析等離子體特性隨時(shí)間的演化過(guò)程。粒子模擬不需要電子和離子處于平衡態(tài)的假設(shè)，因此可以對(duì)鞘層區(qū)域進(jìn)行精確求解。

1.1 模擬模型

本文選用Georgieva和Bogaerts等[13]的CCP源模型。模擬過(guò)程中假設(shè)只有電場(chǎng)，并無(wú)磁場(chǎng)，位置方向?yàn)橐痪S，速度方向?yàn)槿S。模擬中腔體長(zhǎng)度zn=2.5 cm。zn=0為接射頻電壓極板，zn=2.5 cm為接地極板，腔體內(nèi)工作氣體為氬氣，氣體初始溫度為300 K，氣壓為26.66 Pa，初始電子和離子溫度分別為2 eV和0.043 eV。模擬過(guò)程中，射頻源電壓為f=V0sin(wt) (V0=100 V,w=2pf)，本次模擬中頻率f選為 13.56 MHz，并將腔體一維長(zhǎng)度方向劃分為100個(gè)格點(diǎn)，選用的時(shí)間步長(zhǎng)為3.7′10–11s(電子)和 9.25′10–10s(離子)。

1.2 模擬流程

圖1 PIC-MCC模擬單個(gè)循環(huán)流程Fig.1 PIC-MCC computational cycle, one time-step.

1.2.1 模擬初始化

初始值設(shè)置包括空間位置初始化和速度初始化。模擬開始前，電子和Ar+離子隨機(jī)均勻的分布在每個(gè)網(wǎng)格點(diǎn)上，帶電粒子速度分布服從麥克斯韋分布。

1.2.3 電場(chǎng)和電勢(shì)的計(jì)算

用粒子云分室法(CIC)[12]計(jì)算每個(gè)格點(diǎn)的電荷密度，其代入泊松方程，求得它們的電勢(shì)和電場(chǎng)。

1.2.4 帶電粒子受力

用CIC插值法求解每個(gè)帶電粒子所處位置的電場(chǎng)大小。設(shè)帶電粒子所處的位置為zpi，相鄰兩格點(diǎn)的位置和電場(chǎng)大小分別為zj,zj+1和Ej,Ej+1. 則zpi處的電場(chǎng)大小為：

由F=Eq得到此處帶電粒子所受的電場(chǎng)力為：

1.2.5 帶電粒子在電場(chǎng)中的運(yùn)動(dòng)

帶電粒子在電場(chǎng)力作用下運(yùn)動(dòng)，服從牛頓定律，因此，

用 leap_frog法[13]得到每個(gè)帶電粒子的新位置和速度。

1.2.6 邊界條件

模擬中所選模擬區(qū)域?yàn)樯漕l極板與接地極板之間，邊界為兩極板(zn=0,zn=2.5 cm)。射頻極板處(Zn=0)的電勢(shì)隨時(shí)間呈周期性變化f=Vsin(wt)，接地極板處(zn=2.5 cm)電勢(shì)f= 0。

1.2.7 MCC模型

用MCC模擬粒子在運(yùn)動(dòng)過(guò)程中的碰撞過(guò)程，判定是否發(fā)生碰撞，發(fā)生何種碰撞。在一個(gè)時(shí)間步長(zhǎng)Dt內(nèi)，第i個(gè)粒子發(fā)生碰撞的概率為：

式中vi為帶電粒子速度，sT為碰撞截面，Ei為帶電粒子能量，ni為 Ar原子密度。求解Pi值即可判斷粒子運(yùn)動(dòng)一段時(shí)間后是否會(huì)發(fā)生碰撞。取隨機(jī)數(shù)R，若R>Pi，則不發(fā)生碰撞。若RPi，則發(fā)生碰撞。

1.3 碰撞反應(yīng)

在Ar放電中，本文考慮以下5種反應(yīng)[13-14]：

表1 模擬中碰撞反應(yīng)的類型和閾值Table 1 All collisions considered in the model, with the corresponding thresholds for the cross-section data

2 結(jié)果與討論

圖2為射頻極板處鞘層平均厚度與時(shí)間關(guān)系。由模擬初始條件可知，未加射頻電壓前，極板間粒子均勻分布且作無(wú)規(guī)則熱運(yùn)動(dòng)。施加射頻電壓后，在極短時(shí)間內(nèi)射頻極板處形成厚0.28 cm的鞘層。初始階段鞘層厚度隨時(shí)間增加，經(jīng)過(guò)2.65′10–4s趨于穩(wěn)態(tài)，穩(wěn)態(tài)鞘層平均厚度約為0.41 cm。

圖2 射頻極板處鞘層厚度隨時(shí)間的變化關(guān)系Fig.2 Sheath thickness as a function of time at the powered electrode.

圖3 (a)為電子和Ar+離子在射頻極板處的平均粒子通量隨模擬時(shí)間的變化，Ar+離子轟擊射頻極板上的粒子通量隨時(shí)間增加而迅速減少。經(jīng)2.65′10–4s后，極板處的Ar+離子和電子平均粒子通量趨于平衡。平衡后 Ar+離子的通量(6.8′1017m–2·s–1)大于電子的通量(6.4′1017m–2·s–1)。圖 3(b)為電子和Ar+離子在腔體中心處的平均粒子通量隨模擬時(shí)間的變化，腔體中心處的粒子通量隨時(shí)間的增加迅速減少，1.03′10–4s后趨于平衡。4.12′10–4s后腔體中心處的電子平均通量與 Ar+離子平均通量幾乎相等。比較圖3(a)、(b)可知，極板處Ar+離子的平均粒子通量大于腔體中心處。

圖3 電子和Ar+離子在射頻極板上(a)和在腔體中心處(b)的粒子通量隨模擬時(shí)間的變化關(guān)系Fig.3 Flux of electron and Ar+ as a function of time, at the powered electrode (a) and middle of the reactor (b).

圖4 為射頻極板處和腔體中心處平均電流密度隨時(shí)間的變化。由圖 4(a)，射頻極板處的平均電流密度隨時(shí)間的增大而減小，2.65′10-4s后平均電流密度趨于平衡。當(dāng)射頻極板處平均電流密度趨于穩(wěn)定時(shí)，電子的平均電流密度幾乎為零，而Ar+離子的平均電流密度為 3.0′10–3A/m2。圖4(b)是腔體中心處的電子和 Ar+離子平均電流密度隨時(shí)間的變化，在腔體中心處Ar+離子平均電流密度基本保持不變(7.0′10–6A/m2)，而電子平均電流密度隨時(shí)間的增加逐漸減小，在1.03′10–4s之后趨于平衡，達(dá)到0.32 A/m2。從圖中還可以看出，腔體中心處電子的總電流密度大于Ar+離子。

圖4 射頻極板處(a)和腔體中心處(b)電流密度隨模擬時(shí)間的變化Fig.4 Current density as a function of time, at the powered electrode (a) and middle of the reactor (b).

圖5 所示為一個(gè)放電周期內(nèi)，電場(chǎng)強(qiáng)度在不同相位(i=0, p/2, p,3p/2)時(shí)的空間分布。在等離子體區(qū)域(0.75–1.75 cm)，電場(chǎng)強(qiáng)度很小且不隨時(shí)間發(fā)生變化。在鞘層區(qū)域和預(yù)鞘層區(qū)域(分別為射頻極板附近0.75–0 cm和接地板附近1.75–2.5 cm)，電場(chǎng)隨時(shí)間發(fā)生變化。在i=0時(shí)，z=0 和z=2.5 cm處，電場(chǎng)分別為–14387 V/m，15433 V/m 。在i=p/2時(shí),z=0和z=2.5 cm 處，電場(chǎng)分別為–26255 V/m，3512.5 V/m。在i=p時(shí)，z=0 和z=2.5 cm處，電場(chǎng)分別為–15190 V/m，14575 V/m。在i=3p/2時(shí)，z=0 和z=2.5 cm處，電場(chǎng)分別為–3296.2 V/m，26466 V/m。

圖5 單個(gè)射頻放電周期內(nèi)電場(chǎng)在不同相位的空間分布Fig.5 Electric field at different phases of an RF cycle.

射頻放電中極板間的電流由位移和傳導(dǎo)電流組成。位移電流是由隨時(shí)間變化的電場(chǎng)產(chǎn)生的，傳導(dǎo)電流的主要載流子為電子[15]。鞘層內(nèi)(0–0.5 cm 和2–2.5 cm)電子密度很低(圖7a)[1,3,7]，而電場(chǎng)強(qiáng)度隨時(shí)間變化(圖5)，則此區(qū)域內(nèi)電流密度主要由位移電流貢獻(xiàn)，故Ar+離子的平均電流密度大于電子的平均電流密度(圖 4a)。在等離子體區(qū)域(0.75–1.75 cm)，電場(chǎng)強(qiáng)度幾乎不隨時(shí)間發(fā)生變化(圖5)，這意味著位移電流為零，則該區(qū)域的電流主要由傳導(dǎo)電流貢獻(xiàn)，故此處的Ar+離子平均電流密度小于電子的平均電流密度(圖4b)。

在 e+Ar→2e+Ar+電離反應(yīng)中，有大量 Ar+離子產(chǎn)生，同時(shí)模擬中也有大量Ar+離子轟擊到極板上而損失掉。圖6(a)、(b)和(c)分別為腔體內(nèi)產(chǎn)生和損失的Ar+離子數(shù)密度以及腔體內(nèi)總的Ar+離子數(shù)密度隨時(shí)間的變化關(guān)系。由圖可知，隨模擬時(shí)間的增加，參加電離反應(yīng)產(chǎn)生的Ar+離子數(shù)密度增多，而損失的Ar+離子數(shù)密度和腔體內(nèi)的Ar+離子數(shù)密度減小。1.03′10–4s時(shí)，產(chǎn)生和損失的Ar+以及腔體內(nèi)的離子數(shù)密度趨于穩(wěn)定。

圖6 離化反應(yīng)產(chǎn)生的(a)、損失的(b)和腔體內(nèi)的(c)Ar+離子數(shù)密度Fig.6 Ar+ ion densities generated from ionization reaction(a), lost at the boundaries(b), and in the plasma reactor(c).

圖7給出了等離子體達(dá)到平衡后等離子體特性參量對(duì)時(shí)間平均后的空間分布。圖 7(a)是等離子體達(dá)到平衡時(shí)電子和Ar+離子數(shù)密度的空間分布，在中心等離子體區(qū)(0.75–1.75 cm)和預(yù)鞘層區(qū)(0.5–0.75和 1.75–2.0 cm)，電子和離子密度相等，而鞘層區(qū)域(0–0.5 cm和2–2.5 cm)電子密度明顯小于離子密度。圖7(b)為等離子體達(dá)到平衡時(shí)由空間電荷形成的電勢(shì)在空間的分布，鞘層區(qū)域相對(duì)中心等離子體區(qū)域呈負(fù)電勢(shì)，中心等離子體區(qū)電勢(shì)基本保持不變，而鞘層區(qū)域電勢(shì)迅速下降。由式(1)，由于中心等離子體區(qū)凈電荷密度r=0(圖7a)，因此電勢(shì)基本保持不變；而在鞘層區(qū)域r>0(圖7a)，極板間凈電荷為正電荷，因此電勢(shì)迅速下降。圖 7(c)為等離子體達(dá)到平衡時(shí)電場(chǎng)強(qiáng)度的空間分布，在腔體中心處電場(chǎng)強(qiáng)度幾乎為零，而在鞘層區(qū)域電場(chǎng)強(qiáng)度變化顯著。由式(2)，在中心等離子體區(qū)域電勢(shì)無(wú)變化，因此電場(chǎng)強(qiáng)度很小且基本不變；而在鞘層區(qū)域，電勢(shì)的迅速下降導(dǎo)致了電場(chǎng)強(qiáng)度的劇烈變化。圖7(d)為平衡時(shí)電子和 Ar+離子能量的空間分布，Ar+離子的能量在等離子體區(qū)域僅0.04 eV左右。因此，在等離子體區(qū)域，粒子的熱運(yùn)動(dòng)占優(yōu)勢(shì)。離子進(jìn)入鞘層區(qū)域后受到鞘層電場(chǎng)的加速作用(圖 7c)，其能量明顯升高，最高可達(dá)4.58 eV(在兩個(gè)極板處)。電子則在剛進(jìn)入鞘層階段能量有所增加，最高能量為3.41eV(0.25和2.25cm處)，當(dāng)電子進(jìn)入鞘層之后，由于鞘層的阻礙作用其能量迅速下降。

圖7 Ar+離子和電子的平均密度(a)、電勢(shì)(b)、電場(chǎng)強(qiáng)度(c)、平均能量(d)Fig.7 Density(a), potential(b), electric field(c) and averaged energy (d) of electron and Ar+ ion.

3 結(jié)論

本文通過(guò) PIC-MCC方法利用一維模型模擬了CCP源氬氣放電過(guò)程中等離子體的動(dòng)力學(xué)行為，通過(guò)對(duì)模擬結(jié)果分析，得到以下結(jié)論：

1) 射頻極板處鞘層在極短時(shí)間內(nèi)形成，其厚度隨著時(shí)間的增加而增厚，2.65′10–4s時(shí)鞘層厚度基本穩(wěn)定，厚度為0.41 cm。

2) 射頻極板處粒子通量隨時(shí)間的增加而逐漸減小，之后射頻極板處粒子通量達(dá)平衡，整個(gè)模擬中極板處離子通量略大于電子通量。腔體中心處的粒子通量也隨時(shí)間的增加逐漸減少最后逐漸趨于平衡且電子平均通量幾乎與Ar+離子平均通量相等。

3) 等離子體區(qū)域電流密度主要受傳導(dǎo)電流影響，其電子平均電流密度大于Ar+離子平均電流密度；鞘層區(qū)域主要受位移電流影響，其電子平均電流密度小于Ar+離子平均電流密度。

1 趙化僑. 等離子體化學(xué)與工藝[M].合肥: 中國(guó)科學(xué)技術(shù)大學(xué)出版社. 1993.ZHAO huaqiao. Plasm chemistry and technology[M].Hefei: University of Science and Technology of China Press. 1993

2 王德真, 張建紅, 宮野. 氣體放電陰極鞘層氬離子的自洽蒙特卡洛模擬[J]. 計(jì)算物理. 1995, 12(4): 483–488 WANG Dezhen, ZHANG Jianhong, Miyano. A self-consistent monte carlo simulation of ions in cathode sheath of argon gas discharges[J]. Computatinal physics,1995, 12(4):483–488

3 WANG Shuai, XU Xiang, WANG Younian. Frequency matching effects on characteristics of bulk plasmas and sheaths for dual-frequency capacitively coupled argon discharges: one-dimensional fluid simulation[J]. Plasma Sci.Technol, 2008, 10: 57–60

4 WANG Shuai, XU Xiang, WANG Younian. Numerical investigation of ion energy distribution and ion angle distribution in a dual-frequency capacitively coupled plasma with a hybrid model[J]. Physics of Plasmas, 2007,14: 113501

5 嚴(yán)心一. 薄膜技術(shù)[M] .北京: 兵器工業(yè)出版社, 1994 YAN Xinyi. Thin film technology[M]. Beijing: Weapons industry press, 1994

6 楊邦昌, 薄膜物理與技術(shù)[M]. 成都: 電子科技大學(xué)出版社. 1994 YANG Bangchang. Thin film physics and technology[M].Chengdu: University of Electronic Science and Technology of China Press. 1994

7 邱清泉, 勵(lì)慶孚. 工作參數(shù)對(duì)平面直流磁控濺射放電特性的影響[J]. 核聚變與等離子體物理. 2009, 29(2):182–187 QIU qingquan, LI qingfu. Influence of operating parameters on discharge characteristics of planar DC magnetron[J]. Unclear fusion and plasm physics, 2009.29(2): 182–187

8 Lee J K, Manuilenko O V, Babaeva N Y,et al. Ion energy distribution control in single and dual frequency capacitive plasma sources [J]. Plasma Sources Science &Technology, 2005. 14(1): 89–97

9 Lee J K, Babaeva N Y, Kim H C,et al. Simulation of capacitively coupled single- and dual-frequency RF discharges[J]. Ieee Transactions on Plasma Science, 2004,32(1): 47–53

10 Nanbu K, Kondo S. Jpn J. Analysis of three-dimensional dc magnetron discharge by the partial-in-cell/Monte Carlo method[J].App1.Phys.Part l, 1997, 36(7B): 4808–4813

11 Tao W F, Cai D S, Lembege B. et.al. Computer Physics Communications[J],2008. 179(12): 855-864.

12 Birdsall C K, Langdon A B. Plasma Physics via Computer Simulation[M].McGraw- Hill, New York, 1985

13 Georgieva V, Bogaerts A, Gijbels R. Numerical study of Ar/CF/N discharges in single-and dual-frequency capacitively coupled plasma reactors[J]. Journal of Applied Physics, 2003, 93(5): 2369–2379

14 楊津基. 氣體放電[M]. 北京:科學(xué)出版社, 1983 YANG jinji. Gas discharge[M]. Beijing: Science Press,1983

15 Liebermanallan M A, Lichtenberg J. Principles of plasma discharges and materials processing[M]. New Jersey,2005