PIC_MCC程序手册（精编版）

1、PIC-MCC程序手册1. PIC_MCC的模拟方法和数值计算1.1 PIC_MCC的模拟原理1.2 PIC 模型1.3 MCC模型1.4 模拟中所涉及的放电粒子1.5 模拟中的碰撞2. 数据分析程序正常运行所需文件及文件意义主要输出文件的意义1PIC_MCC的模拟方法和数值计算的模拟原理目前在计算机模拟中大量采用低压射频放电模型来模拟材料的加工及改性。但在低压情况下， 粒子和中性气体的碰撞不足，无法使其达到平衡，在这种情况下我们认为放电所产生 的负离子及电子的速度已经偏离了Maxwellian分布。因此，流体力学的模拟方法已经无法 准确的解决此类问题。我们选用一种新的方法，即用运动分析的方法

2、来解决低压等离子反应 器中的物理和化学过程，用包含大量粒子的模型来解决Boltzmann 方程。 粒子模型不仅可以解决刻蚀在基板上粒子的能量问题，而且还可以很好的解决粒子刻蚀率和各向异性的问题。Monte Carlo算法与 PIC 模拟方法的有机结合就形成了PIC_MCC模拟方法。在PIC_MCC中， 我们假设中性气体是在时间与空间位置上的一种特定分布。PIC_MCC模拟中采用的是一维空间（ Zn），速度方向为三维（Vx， Vy，Vz）。PIC模型在 PIC 模型中（如图a），粒子在电场力的作用下运动。粒子模拟只能解决少量粒子存在的模型， 这个模型中的粒子数量远远小于真实情况下等离子体中的粒子

3、数目，模拟中的每69个粒子即超粒子表示10 10 个粒子。 在模拟中我们必须有足够多的超粒子，以减少粒子离散及噪声扰动。超粒子与所划分的网格点数之比必须大于1。在模拟中，我们主要解决Maxwell 方程以及F=ma=q（E+v*B）。电场可以通过Maxwall 方程解出。 粒子在电场和磁场中所受的力可通过Newton_Lorentz方程解出。图 1.2.1模拟中所用的系统模型（ a），系统模型 z 轴方向上网格点的划分。0 点为接射频电极极板端， ZN 为接地电极极板端。1.2.1 PIC模型中一个时间步循环的数值计算图 1.2.2 PIC一个时间步的循环运算（ 1）电荷密度 j被分配到每个网

4、格点j 上，这个过程称之为电荷分配。因此， 先从连续的节点Zpi然后再到离散的节点Zj来计算电荷密度j。电荷分配函数可用S（ Zj-Zpi）来表示， 包括零节点， 第一节点和最后节点。图 1.2.3中所描述的是第一节点的电荷分配方式。这种分配方式把Zj节点上的j单胞和 Zj+1节点上的 （ j+1 ）单胞这部分电荷看做带电粒子或电荷云。这种带电粒子看做是一种有限度的刚性电荷云，他们可以在通过彼此时不受束缚而自由运动，这种模型我们称之为cloud_in_cell或者 CIC.图 1.2.3如果电荷粒子的密度是一定的，j 和 j+1 之间的距离为Z，那么电荷粒子qpi分配给节点j 的电量为：分配给

5、节点j+1 的电量为：因此，在Zpi 上的电荷粒子qpi分配给 j 节点的电荷密度j 为：（ 2）电荷密度可以用来计算网格点上的电场E。在静电场模拟中， *E= - ? B/ ? t 0, 所以 E=- ， 由一维条件下的Poisson方程可以得到：电场可以由以下公式计算得出：（ 3）E 又按照函数S（Zj -Z pi ）分配给网格点上的粒子。在一维的静电场模型中， 电场分配各网格点的电场为：带电粒子所受的电场力为F=qE，一维静电场模型中：（ 4）运动方程可以计算出带电粒子新的位置和速度。在一维静电场模型中，可用以下方程代替上述运动方程：因为带电粒子的速度v 和位置 x 是不能同时确定的，所

6、以leap_frog算法要采用不同模式原则。图 1.2.2 leap_frog算法的网格点划分示意图。应注意到初始条件下带电粒子在时间 t=0 时的速度是需要改变的，把 v（0）处的速度 V 退回到 v（ - t/2 ）处，然后通过带电粒子所受的电场力还计算 t=0 时的速度。（ 5）检查边界条件，检查粒子是否附着在极板上。初始时间t=0 时粒子的位置和初始- t/2时的速度已经给出， 带电粒子的密度也可通过计算得出。图1.2.2中所描述的（ 1）到（ 5）只是重复的循环， 直到等离子体达到收敛。1.2.2 边界处等离子体粒子的模拟在射频放电产生等离子体的模拟中，不仅要考虑中心等离子体处的粒子

7、行为，也要模拟边界处，即鞘层处的粒子行为。位势方程的边界条件可通过 Gauss 法则得出：S等离子体区域和上下两极板的总面积A0下极板（接射频电压的极板）的表面区域面积AN上极板（接地电极极板）的表面区域面积 0下极板（接射频电压的极板）的带电粒子密度 N上极板（接地电极极板）的带电粒子密度网格点的电势可通过以下方程计算得出:j=1， 2， N-1， N 为所划分的网点。一维系统的边界条件为：1.3 MCC model1.3.1 无碰撞的模拟方法PIC 模拟方法是一种碰撞模型。即使在低电压情况下带电粒子和中性气体的碰撞也对维持放电起着非常重要的作用。碰撞可以将 PIC 和 MC两种方法结合起来

8、进行模拟运算。 PIC 模拟的是所有粒子在同一时间步的运动。 而 MC方法模拟的是在碰撞中一些随机粒子的行为，只对每个时间步中的部分粒子的行为进行模拟，我们称之为 MCC模拟方法。1.41.5 电子和中性粒子的碰撞1.5.1 碰撞截面的数据我们认为中性气体（Ar ，CF4 和 N2）的分布是均匀的，它们的速度分配在室温情况下（ Tg= eV 或者 300K）是服从麦克斯韦分布的。因此，中性气体粒子和电子（平均 Tg>2 eV ）相比所具有的能量很小，我们认为它们是静止的。碰撞截面我们用（ g）来表示， g 是粒子在碰撞之前的相对速度。在模型中所有中性粒子的碰撞截面数据和相应的阈值在表中给

9、出。Type of collisionReactionThreshold （ eV）ArElastic scatteringe+Ar e+ArTotalelectronic excitatione+Ar e+Ar*+Ionizatione+Ar 2e+ArCF4Momentum transfere+CF4e+CF4 Vibrational excitatione+CF4e+CF4Vibrational excitatione+CF4e+CF4Vibrational excitatione+CF4e+CF4Electronic excitatione+CF4e+CF4*-Electron att

10、achmente+CF4F+CF35-Electron attachmente+CF4F+CF3 5-+Dissociatione+CF4e+F +CF3 12Dissociativeionization+e+CF42e+F+CF316Neutral dissociatione+CF4e+F+CF312Neutral dissociatione+CF4e+2F+CF217Neutral dissociatione+CF4e+3F+CF18N2Momentum transfere+ N 2e+ N 2*aexcitatione+ N 2 e+ N 2 (Y)+bIonizatione+ N 22

11、+ N 2 (Y)+Ionizatione+ N 2 2+ N2(Y)(B )18.21.5.2 粒子碰撞后的速度计算方法一般情况下，我们所考虑的是两个均匀球体粒子的之间的弹性碰撞。粒子在碰撞之前，我们设它们的质量和速度分别为，m和 M，v 和 V，它们之间的相对速度 为 g=v-V 。在不考虑一般情况下的损失，我们认为在碰撞之前的系统条件为：V=0， v=g。碰撞之后的速度为v ， V， g = v V。因为碰撞前后系统的总动量守恒，在这种条件下，可以计算出质心的速度VCM （图）。因为在质心坐标系中两个粒子的初始动量大小相等，方向相反，所以，在系统中，两个粒子所受的力大小相等，方向相反。所

12、以，碰撞后的动量也是大小相等，方向相反。碰撞之后粒子的速度方向仍然是平行的，但是却偏离了一定角度，如图图 粒子在质心坐标系下的运动轨迹。为散射角要求出粒子碰撞之后的速度，就要先求出散射角。如果可以求出散射角，那么可以通过以下公式来求出粒子碰撞之后的速度，在质心坐标系中粒子运动轨迹所在的平面叫做碰撞平面。碰撞平面和参考平面之间的夹角为，参考平面是任意选取的，因此夹角为：R 是在特定分布 0,1之间的一个随机数。两个粒子碰撞之后的速度有两个主要影响因素，第一个是角度，第二个是b。图 1.1.2碰撞参数和b系统要维持通量守恒，所以通过环面2 db 和立体角2 sin d的总粒子数相等。上式中负号表示

13、如果增加b，会相应的减小。db/d 在这里去绝对值，因为其为负值。散射角与两粒子之间的电势和速度有关系，不同的散射角，碰撞截面不同。（）（）不同的碰撞截面取决于粒子之间的相互作用力和电势。例如，在电子和Ar 的碰撞中，在计算电势时可以忽略掉库伦电势的影响（弹性碰撞，激发，电离）。不同情况下的碰撞截面由以下公式给出，（） =E/E0 是无量纲的能量，E 是电子的相对能量，E0 原子的单位能量（E0 = eV ） .散射角与随机数R 及之间的关系是：（）这样就可以求出散射角，而粒子弹性碰撞之后的速度也可通过公式（）和（） 求出。角度可通过公式（）求出。多数情况下的电子和中性粒子碰撞中，公式（）和（

14、） 中的 M+mM， g v。a) e_Ar excitatione+Are+Ar*根据动量和能量守恒公式，（）（）Eth 是非弹性碰撞的阈值。M+m M， g vv为v 的标量，（）2E=mv/2是电子碰撞之前的能量。激发过程如果看做是一个弹性碰撞过程，碰撞前的速度为v 何 V。碰撞之后的速度可以由公式（）和（）求出。公式中所有的v用 v 来代替。b) electron_AionizationA 是质量为M的中性粒子。（）+A 为碰撞后的离子，e1 是初始电子， e2 是发射出的电子。能量守恒公式，（）Eth 为电离反应的阈值。因为入射离子与电子的质量比很大，我们可以认为电子的动量远小于中性

15、粒子的动量，电子与中性粒子碰撞之后，电子偏离，中性粒子变成离子， 碰撞之后的轨迹仍然是未被干扰的。这种假设一方面使得在碰撞之前中性粒子的速度V =V；另一方面，能量守恒公式表示为，（）等式的左边我们看做是电离能量的改变量E，需要找到一种方法把散射电子和发射电子区分开。（）（）（）2Einc =mv/2为初始电子的能量。1 和 0 单位是电子伏特。当E 区分开后，我们可以从公式（）中来求得v。（）类似于激发情况，（）碰撞之后散射电子的速度v可有公式（）求出。式中v 用 v 来代替。碰撞之前，发射出的电子的能量实际是不存在的，我们假设为，（）碰撞之后的速度可由公式（）可算出，用v ej 来代替 v

16、， v ej 来代替 v。由于很多反应的碰撞截面数据是无法在文献中查询到的。如果没有碰撞截面的数据，我们就无法求出相应的散射角。总的碰撞截面可由公式（）得出 T =4（ g），这样可以得到sin d d/4 。意味着碰撞之后的速度方向随机，我们称这种性质的散射为各向同性。因为 electron_CF4 和 electron_N2 碰撞的碰撞截面数据是没有的，因此我们认为它们之间的碰撞是各向同性的。各向同性的散射角在区间0, ，可表示为，（）角度可从公式 （）求出， 而 electron_CF4 和 electron_N2 弹性碰撞之后的速度可由公式（）和（）求出。与electron_Ar的弹性

17、碰撞相似，M+m M ， g v。1.6 Ion_neutral collisions1.6.1 Cross_section data(a)Collisions of Ar+ with neutrals+因为所有Ar参与的反应的碰撞截面数据均已给出，我们用null_collision方法来处理此类问题。1) 主程序读取文件，并计算gas density gas density=pertorr*pressure/gtemp2) 主程序打开 data/ 文件，读取 e 和 Ar 的碰撞截面数据3) 根 据 语 句 “ 来 判 断 主 程 序 读 取 哪 些 碰 撞 截 面 数 据 文件。. arg

18、on/CF4)4) 主程序读取 CF4ion 与 CF4反应所需能量的文件（70，file=“ data/”）（70，file=“ data/”）（70，file=“ data/”）5) Spec(i)所代表的粒子分别为：spec（ 1）-espec（2）-Arspec（3）-CF3+spec（4）-F-spec（5）-CF3-a) 判断语句 if (nflag_col_spec(1).，由于此时只有e 一种离子存在，因此调用子程序prob_e_ngas(),来模拟计算电子和所有中性气体粒子的碰撞反应。if(nflag_col_spec(1). then call prob_e_ngas()e

19、nd ifb) 对模拟中气体种类进行判断（） ，然后对不同的species以及spec（ i ）选择调用不同的子程序来进行碰撞模拟。if.'argon'.'argon/CF4'.'argon/CF4/*N2')then if(nflag_col_spec(2).thencall prob_Ar_Ar() end ifend ifccprobability of various ion-molecule collisions;beta_inf=3. polar_cf4= polar_Ar =if.'argon/CF4'.'a

20、rgon/CF4/N2') then dois=3,nspecccprobabilityofCF4ions+ CF4reactiveand ellasticcollisionsif(nflag_col_spec(is). thenprob_col_spec(is)=e*sqrt(polar_cf4*PI*(mass(is)+* mass(0)/mass(is)/mass(0)/EPS0)*beta_inf*beta_inf*ratio_gas(2)*gden*dtis(is) cc*t henprob_col_ArCF=e*sqrt(polar_cf4*PI*(mass(2)+mass

21、(0)/mass(2)*/mass(0)/EPS0)*beta_inf*beta_inf*ratio_gas(2)*gden*dtis(2)end ifccprobabilityofCF4ions+ CF4reactiveand ellasticcollisions ccin case of pure CF4 dischargeelse if.'CF4') thendo is=2,nspecif(nflag_col_spec(is). thenprob_col_spec(is)=e*sqrt(polar_cf4*PI*(mass(is)+* mass(0)/mass(is)/m

22、ass(0)/EPS0)*beta_inf*beta_inf*ratio_gas(2)*gden*dtis(is)end if end doend if cc*a rgon/CF4/N2') thencall prob_Ar_N2 call prob_N2_N2end ifccwrite all probabilitieswrite(*,*)'prob',(prob_col_spec(i),i=1,nspec) write(*,*)if.'argon/CF4'.'argon/CF4/N2')thenwrite(*,*)'probC

23、F4ions- Ar',(prob_col_CFAr(i),i=3,nspec)write(*,*)'prob Ar+ - CF4',prob_col_ArCFwrite(*,*)'probAr+-N2charge exchange',prob_col_ArN2write(*,*)'probN2+-Archexandel col',prob_col_N2Arwrite(*,*) end ifcc*thenif.'argon')thencall mc_e_ngas(-1,2,-1,-1,-1,-1)else if.'

24、argon/CF4')thencall mc_e_ngas(0,2,3,4,5,-1) else if.'argon/CF4/N2')then call mc_e_ngas(0,2,3,4,5,6)else if.'CF4')thencall mc_e_ngas(0,-1,2,3,4,-1) end ifend ifcc*a rgon'.'argon/CF4'.or.*.'argon/CF4/N2')thencc*t hencall mc_Ar_Ar()end ifend ifcc*a rgon/CF4/N2

25、9;) thencall mc_N2_N2 ()end ifcc*a rgon/CF4'.'argon/CF4/N2')then cc*t hencallmc_ion_neut (3,0,dif_engy_cf3p,n_reac_cf3p,* num_col_ij_3,extra_col_ij_3,rate_col_ij_3)call mc_cf_Ar(3,2,num_col_ij_3,rate_col_ij_3) end ifcc*t hencall mc_ion_neut(4,0,dif_engy_fm,n_reac_fm,* num_col_ij_4,extra_

26、col_ij_4,rate_col_ij_4)call mc_cf_Ar(4,2,num_col_ij_4,rate_col_ij_4) end ifcc*thencallmc_ion_neut (5,0,dif_engy_cf3m,n_reac_cf3m,* num_col_ij_5,extra_col_ij_5,rate_col_ij_5)call mc_cf_Ar(5,2,num_col_ij_5,rate_col_ij_5) end ifcc*C F4')then cc* thencallmc_ion_neut(2,0,dif_engy_cf3p,n_reac_cf3p,* n

27、um_col_ij_3,extra_col_ij_3,rate_col_ij_3)end ifcc* thencall mc_ion_neut(3,0,dif_engy_fm,n_reac_fm,* num_col_ij_4,extra_col_ij_4,rate_col_ij_4)end ifcc*thencall mc_ion_neut(4,0,dif_engy_cf3m,n_reac_cf3m,* num_col_ij_5,extra_col_ij_5,rate_col_ij_5)end ifcc*mv'MV 'mvMV1(v '2V ')2Eth1(vV

28、 )22Takingintoconsiderationthat M+m M and g v we obtain(v 'v22Eth mEthvv1EWhere E=mv2/2 is the electron energy before collision. Theexcitationprocessistreatedas ifitwere an elasticcollisionwithpre_collisionvelocitiesvandV.Thepost_collisionvelocities are given byin which all v s are replaced byv

29、.a) Ionization(electron_A)Where A denotes the neutral particle with a mass , which canbe assumedequaltothatofthe ionM.The processisrepresentedase1 (v)A(V )e1 (v ')e2 ( v '')A(V ')Where A+ denotes an ion,e1 is the incident electron , e2 is the ejected electron , and the symbols in the

30、 parentheses denotethe velocities .The energy balance equation ismv2MV 2mv '2mv '2MV'2ejE22222thWhere Eth is the threshold energy of the ionization . Becauseofthelargeion_to_electronmass ratio,we can assumethatthe momentumoftheincidentelectronismuch lowerthanthe momentum oftheneutralpart

31、icle,.theincidentelectron removes an electron from the neutral ,and the neutral becomesanion,continuingonitstrajectronundisturbed .On one hand ,this assumption means that the created ion takes thevelocityanddirectionoftheneutralparticlebeforecollision ,.V=V .On theother hand,the energy balancecanbe

32、rewrittenmv2ejEthmv'2mv'2222The left_hand side is known as the excess energy E after ionization and we need to find an algorithm how to divide itinto the scattered and ejected electrons. In general, whenthereisno publishedwork on thedivisionoftheexcess energy , it is divided randomly into tw

33、omv '2EscatRE2ejmv'2Eej(1R)E211For the electron_Ar ionzation ,however,we use the expression for sampling the Eej by a random number R,Eej0tan R(tantan10 )tan10 2WhereEinc =mv/2istheenergy ofthe incidentelectron.The units of 0 and 1are electronvolts .Once the excess energy isdivided, v is cal

34、culated from Eq.v 'v22(Eth mEthvv1Eej )EejEincThe post_collision velocity of the scattered electron v isobtained byEq inwhich v isreplacedbyv .The pre_collisionvelocity for the ejected electron , which does not exist in reality, is assumedvejv2Eejvmand the post_collision velocity is obtained aga

35、in from Eq. byreplacing v withvejand v with v ej .Formanytypesofcollisionnodataondifferential cross_sectionisassumed nottodepend tothedeflectionangle x.Thetotalcross_sectioncalculatedfrom T=4 （ g）.Therefore , the probabilityDefinedbybecomessinxdxd /4 .Thismeansthatthe post_collisionvelocitiestakeran

36、domdirection.The scatteringwiththispropertyiscalledtherearenosufficientdataforthedifferentialcross_sectionsof electron_CF4 and electron_N2 collisions we assume that thescatteringisisotropic.Theapproximation oftheinteraction potentialwiththescreenedCoulombpotentialmadefor collisions with atomic Ar is

37、 not valid for molecularcase ofisotropic scattering the deflection angle , which is in the interval0,is randomly sampled bycos12RTheangleisrandomlysampledfromandthe post_collisionvelocitiesforelectron_CF4andelectron_N2elastic collisions can be found fromto electron_Ar elastic collisions,these equati

38、ons can be somewhat simplified by use of M+mM and g v.Electron_CF4andelectron_N2inelasticcollisionsaretreatedsimilartoelectron_Arinelasticcollisions,taking intoaccounttheisotropicscatteringaftercaseof ionization,the excess energy is divided by a random numberusingEqs .andbetweenthescatteredand eject

39、edelectron.Incase of attachment the incident electron is removed from thecalculationand thecreatednegativeiontakesthevelocityand direction of the neutral particle before the collision.a、Isotropic scatteringi. Cross_section data ii.Reactione+CF4e+CF4e+ N2e+ N2iii. Calculation of the post_collision ve

40、locities（）（）（）iv. subroutine newvel_e_isotrsubroutinenewvel_e_isotr(sengy,svel,svz,svy,svx,n_flag,is1,is2) include"param_p1d"integern_flag,is1,is2real*8sengy,svel,svx,svy,svz,rd,alphareal*8phi1,cosphi,sinphi,coschi,sinchi,up1,up2,up3,*smag,r11,r12,r13,*r21,r22,r23,r31,r32,r33coschi=（）sinch

41、i=dsqrt*coschi)phi1=*PI*rand_gen()（）cosphi=cos(phi1) sinphi=sin (phi1)cc*mass(is1)*/mass(is2)（）svel =svel*dsqrt（）sengy=sengy*（）end if r13=svz r23=svy r33=svxifthen up1= up2=up3= else up1= up2= up3= end ifr12=r23*up3-r33*up2 r22=r33*up1-r13*up3 r32=r13*up2-r23*up1smag=dsqrt(r12*r12+r22*r22+r32*r32)r1

42、2=r12/(smag+ r22=r22/(smag+ r32=r32/(smag+r11=r22*r33-r32*r23r21=r32*r13-r12*r33 r31=r12*r23-r22*r13svz=svel*(r11*sinchi*cosphi+r12*sinchi*sinphi+r13*coschi)svy=svel*(r21*sinchi*cosphi+r22*sinchi*sinphi+r23*coschi)svx=svel*(r31*sinchi*cosphi+r32*sinchi*sinphi+r33*coschi)cc*mass(is1)*/mass(is2)（）svel

43、 =svel*dsqrt（）sengy=sengy* end ifr13=svz r23=svy r33=svxifthen up1= up2= up3= else up1= up2= up3=end ifr12=r23*up3-r33*up2 r22=r33*up1-r13*up3 r32=r13*up2-r23*up1smag=dsqrt(r12*r12+r22*r22+r32*r32)r12=r12/(smag+ r22=r22/(smag+ r32=r32/(smag+r11=r22*r33-r32*r23 r21=r32*r13-r12*r33 r31=r12*r23-r22*r13

44、svz=svel*(r11*sinchi*cosphi+r12*sinchi*sinphi+r13*coschi)svy=svel*(r21*sinchi*cosphi+r22*sinchi*sinphi+r23*coschi) svx =svel*(r31*sinchi*cosphi+r32*sinchi*sinphi*+r33*coschi)cc*vx( is,i)cc*c csubroutine maxwvel(svx,svy,svz,ip,np,nflag) include"param_p1d"real*8svx,svy,svz,svr,stemp,theta,ra

45、ndom real*8tvx,tvy,tvz,phiintegerip,np,nflagreal*8sinphi,cosphi,sinthe,costhe,vel_trm cc*PI*rand_gen()random=rand_gen()do whilerandom=rand_gen() end dovel_trm=dsqrt(e*stemp/mass(ip) svr = vel_trm*dsqrt*log(random)svx = svr*dcos(theta) svy = svr*dsin(theta)random=rand_gen()do whilerandom=rand_gen() end dosvz = vel_trm*dsqrt*log(random)*sin*PI*rand_gen()cc*mv'MV 'mvMV1(v '2V ')2Eth1(vV )22Takingintoconsiderationthat M+m M and g v we obtainv 'v22Ethm(Ethvv1EWhere E=mv2/2 is the electron energy before collision. Theexcitationproce