PIC_MCC程序手册X.doc

PIC-MCC 程序手册1. PIC_MCC的模拟方法和数值计算1.1 PIC_MCC的模拟原理1.2 PIC模型 1.3 MCC模型1.4 模拟中所涉及的放电粒子1.5 模拟中的碰撞2. 数据分析2.1程序正常运行所需文件及文件意义2.2主要输出文件的意义1 PIC_MCC的模拟方法和数值计算1.1PIC_MCC的模拟原理 目前在计算机模拟中大量采用低压射频放电模型来模拟材料的加工及改性。但在低压情况下，粒子和中性气体的碰撞不足，无法使其达到平衡，在这种情况下我们认为放电所产生的负离子及电子的速度已经偏离了Maxwellian分布。因此，流体力学的模拟方法已经无法准确的解决此类问题。我们选用一种新的方法，即用运动分析的方法来解决低压等离子反应器中的物理和化学过程，用包含大量粒子的模型来解决Boltzmann方程。粒子模型不仅可以解决刻蚀在基板上粒子的能量问题，而且还可以很好的解决粒子刻蚀率和各向异性的问题。Monte Carlo算法与PIC模拟方法的有机结合就形成了PIC_MCC模拟方法。在PIC_MCC中，我们假设中性气体是在时间与空间位置上的一种特定分布。PIC_MCC模拟中采用的是一维空间（Zn），速度方向为三维（Vx，Vy，Vz）。1.2 PIC模型在PIC模型中（如图a），粒子在电场力的作用下运动。粒子模拟只能解决少量粒子存在的模型，这个模型中的粒子数量远远小于真实情况下等离子体中的粒子数目，模拟中的每个粒子即超粒子表示106109个粒子。在模拟中我们必须有足够多的超粒子，以减少粒子离散及噪声扰动。超粒子与所划分的网格点数之比必须大于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 被分配到每个网格点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的电量为： 分配给节点j+1的电量为： 因此，在Zpi 上的电荷粒子qpi 分配给j节点的电荷密度j 为： （2） 电荷密度可以用来计算网格点上的电场E。在静电场模拟中，*E=- B/t0,所以E=-，由一维条件下的Poisson方程可以得到： 电场可以由以下公式计算得出：（3） E又按照函数S（Zj-Zpi）分配给网格点上的粒子。在一维的静电场模型中，电场分配各网格点的电场为： 带电粒子所受的电场力为F=qE，一维静电场模型中： （4） 运动方程可以计算出带电粒子新的位置和速度。 在一维静电场模型中，可用以下方程代替上述运动方程： 因为带电粒子的速度v和位置x是不能同时确定的，所以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 边界处等离子体粒子的模拟 在射频放电产生等离子体的模拟中，不仅要考虑中心等离子体处的粒子行为，也要模拟边界处，即鞘层处的粒子行为。 位势方程的边界条件可通过Gauss法则得出： S等离子体区域和上下两极板的总面积 A0下极板（接射频电压的极板）的表面区域面积 AN上极板（接地电极极板）的表面区域面积 0下极板（接射频电压的极板）的带电粒子密度 N上极板（接地电极极板）的带电粒子密度 网格点的电势可通过以下方程计算得出: j=1，2，N-1，N为所划分的网点。 一维系统的边界条件为： 1.3 MCC model1.3.1 无碰撞的模拟方法 PIC模拟方法是一种碰撞模型。即使在低电压情况下带电粒子和中性气体的碰撞也对维持放电起着非常重要的作用。碰撞可以将PIC和MC两种方法结合起来进行模拟运算。PIC模拟的是所有粒子在同一时间步的运动。而MC方法模拟的是在碰撞中一些随机粒子的行为，只对每个时间步中的部分粒子的行为进行模拟，我们称之为MCC模拟方法。1.41.5 电子和中性粒子的碰撞1.5.1 碰撞截面的数据我们认为中性气体（Ar，CF4和N2）的分布是均匀的，它们的速度分配在室温情况下（Tg=0.026 eV或者300K）是服从麦克斯韦分布的。因此，中性气体粒子和电子（平均Tg2 eV）相比所具有的能量很小，我们认为它们是静止的。碰撞截面我们用（g）来表示，g是粒子在碰撞之前的相对速度。在模型中所有中性粒子的碰撞截面数据和相应的阈值在表1.1中给出。Type of collisionReactionThreshold（eV）ArElastic scatteringe+Are+ArTotalelectronic excitatione+Are+Ar*11.5Ionizatione+Ar2e+Ar+15.8CF4Momentum transfere+CF4e+CF4Vibrational excitatione+CF4e+CF40.108Vibrational excitatione+CF4e+CF40.168Vibrational excitatione+CF4e+CF40.077Electronic excitatione+CF4e+CF4*7.54Electron attachmente+CF4F- +CF35Electron attachmente+CF4F+CF3-5Dissociatione+CF4e+F- +CF3+12Dissociative ionizatione+CF42e+F+CF3+16Neutral dissociatione+CF4e+F+CF312Neutral dissociatione+CF4e+2F+CF217Neutral dissociatione+CF4e+3F+CF18N2Momentum transfere+ N2e+ N2excitatione+ N2e+ N2* (Y)aIonizatione+ N22+ N2+(Y)b15.6Ionizatione+ N22+ N2+(Y)+ (B2)18.1.5.2 粒子碰撞后的速度计算方法 一般情况下，我们所考虑的是两个均匀球体粒子的之间的弹性碰撞。粒子在碰撞之前，我们设它们的质量和速度分别为，m和M，v和V，它们之间的相对速度为g=v-V。在不考虑一般情况下的损失，我们认为在碰撞之前的系统条件为：V=0，v=g。碰撞之后的速度为v ，V，g= v V。因为碰撞前后系统的总动量守恒，在这种条件下，可以计算出质心的速度VCM （图1.11）。 (2.31)(2.32) (2.34)因为在质心坐标系中两个粒子的初始动量大小相等，方向相反，所以， (2.36)在系统中，两个粒子所受的力大小相等，方向相反。所以，碰撞后的动量也是大小相等，方向相反。碰撞之后粒子的速度方向仍然是平行的，但是却偏离了一定角度，如图1.11图1.11 粒子在质心坐标系下的运动轨迹。为散射角要求出粒子碰撞之后的速度，就要先求出散射角。如果可以求出散射角，那么可以通过以下公式来求出粒子碰撞之后的速度，(2.38) (2.39)， (1.17)在质心坐标系中粒子运动轨迹所在的平面叫做碰撞平面。碰撞平面和参考平面之间的夹角为，参考平面是任意选取的，因此夹角为： (2.40)R是在特定分布0,1之间的一个随机数。两个粒子碰撞之后的速度有两个主要影响因素，第一个是角度，第二个是b。图1.1.2 碰撞参数和b系统要维持通量守恒，所以通过环面2db和立体角2sind的总粒子数相等。(2.41)上式中负号表示如果增加b，会相应的减小。 (2.42)db/d在这里去绝对值，因为其为负值。(2.43)散射角与两粒子之间的电势和速度有关系，不同的散射角，碰撞截面不同。（2.44）（2.45）不同的碰撞截面取决于粒子之间的相互作用力和电势。例如，在电子和Ar的碰撞中，在计算电势时可以忽略掉库伦电势的影响（弹性碰撞，激发，电离）。不同情况下的碰撞截面由以下公式给出，（2.46）=E/E0 是无量纲的能量，E是电子的相对能量，E0 原子的单位能量（E0 =27.21 eV）.散射角与随机数R及之间的关系是： （2.47） 这样就可以求出散射角，而粒子弹性碰撞之后的速度也可通过公式（1.15）和（1.16）求出。角度可通过公式（1.18）求出。多数情况下的电子和中性粒子碰撞中，公式（1.15）和（1.16）中的M+mM，gv。a) e_Ar excitatione+Are+Ar*根据动量和能量守恒公式，（2.48）（2.49）Eth 是非弹性碰撞的阈值。M+mM，gv(2.50)v为的标量，（2.51）E=mv2/2是电子碰撞之前的能量。激发过程如果看做是一个弹性碰撞过程，碰撞前的速度为何V。碰撞之后的速度可以由公式（1.15）和（1.16）求出。公式中所有的v用来代替。b) electron_A ionizationA是质量为M的中性粒子。（2.52）A+ 为碰撞后的离子，e1 是初始电子，e2 是发射出的电子。能量守恒公式，（2.53）Eth 为电离反应的阈值。因为入射离子与电子的质量比很大，我们可以认为电子的动量远小于中性粒子的动量，电子与中性粒子碰撞之后，电子偏离，中性粒子变成离子，碰撞之后的轨迹仍然是未被干扰的。这种假设一方面使得在碰撞之前中性粒子的速度V=V；另一方面，能量守恒公式表示为，（2.54）等式的左边我们看做是电离能量的改变量E，需要找到一种方法把散射电子和发射电子区分开。（2.55）（2.56）（2.57）Einc =mv2/2为初始电子的能量。1和0 单位是电子伏特。当E区分开后，我们可以从公式（1.33）中来求得v。（2.58）类似于激发情况，（2.59）碰撞之后散射电子的速度v可有公式（1.15）求出。式中v用来代替。碰撞之前，发射出的电子的能量实际是不存在的，我们假设为，（2.60）碰撞之后的速度可由公式（1.15）可算出，用ej 来代替v，v ej 来代替v。由于很多反应的碰撞截面数据是无法在文献中查询到的。如果没有碰撞截面的数据，我们就无法求出相应的散射角。总的碰撞截面可由公式（1.21）得出T =4（g），这样可以得到sindd/4。意味着碰撞之后的速度方向随机，我们称这种性质的散射为各向同性。因为electron_CF4和electron_N2碰撞的碰撞截面数据是没有的，因此我们认为它们之间的碰撞是各向同性的。各向同性的散射角在区间0, ，可表示为，（2.61）角度可从公式（1.18）求出，而electron_CF4和electron_N2 弹性碰撞之后的速度可由公式（1.15）和（1.16）求出。与electron_Ar的弹性碰撞相似，M+mM ，gv。1.6 Ion_neutral collisions1.6.1 Cross_section data(a) Collisions of Ar+ with neutrals因为所有Ar+ 参与的反应的碰撞截面数据均已给出，我们用null_collision方法来处理此类问题。1) 主程序读取inp1d_mc.dat文件，并计算gas densitygas density=pertorr*pressure/gtemp2) 主程序打开data/e_Ar_eltot.dat文件，读取e和Ar的碰撞截面数据3) 根据语句“species.eq.XX”来判断主程序读取哪些碰撞截面数据文件。(species.eq.argon/CF4)4) 主程序读取CF4ion与CF4反应所需能量的文件（70，file=“data/dif_engy_cf3p.dat”）（70，file=“data/dif_engy_fm.dat”）（70，file=“data/dif_engy_cf3m.dat”）5) Spec(i)所代表的粒子分别为：spec（1）-e spec（2）-Ar spec（3）-CF3+ spec（4）-F- spec（5）-CF3-a) 判断语句if (nflag_col_spec(1).eq.1)，由于此时只有e一种离子存在，因此调用子程序prob_e_ngas(),来模拟计算电子和所有中性气体粒子的碰撞反应。if (nflag_col_spec(1).eq.1) thencall prob_e_ngas()end ifb) 对模拟中气体种类进行判断（(species.eq.XX），然后对不同的species以及spec（i）选择调用不同的子程序来进行碰撞模拟。if(species.eq.argon.or.species.eq.argon/CF4.or.species.eq.argon/CF4/* N2) thenif (nflag_col_spec(2).eq.1) then call prob_Ar_Ar() end if end ifcc*cc probability of various ion-molecule collisions; beta_inf=3. polar_cf4=2.9d-30 polar_Ar =1.642d-30 if (species.eq.argon/CF4.or.species.eq.argon/CF4/N2) then do is=3,nspeccc*cc probability of CF4 ions + CF4 reactive and ellastic collisions if (nflag_col_spec(is).eq.1) then prob_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*/cc probability of CF4 ions - Ar ellastic collisions prob_col_CFAr(is)=e*sqrt(polar_Ar*PI*(mass(is)+ * mass(2)/mass(is)/mass(2)/EPS0)*beta_inf*beta_inf* * ratio_gas(1)*gden*dtis(is) end if end docc*cc probability of Ar+ + CF4 ellastic collisions if (nflag_col_ij_2(3).eq.1) then prob_col_ArCF=e*sqrt(polar_cf4*PI*(mass(2)+mass(0)/mass(2) * /mass(0)/EPS0)*beta_inf*beta_inf*ratio_gas(2)*gden*dtis(2) end ifcc*cc probability of CF4 ions + CF4 reactive and ellastic collisions cc in case of pure CF4 discharge else if (species.eq.CF4) then do is=2,nspec if (nflag_col_spec(is).eq.1) then prob_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) end if end do end ifcc*/cc probability of (N2+) + N2 and (N2+) + Ar charge exchange if (species.eq.argon/CF4/N2) then call prob_Ar_N2 call prob_N2_N2 end ifcc*cc write all probabilities write(*,*)prob,(prob_col_spec(i),i=1,nspec) write(*,*) if (species.eq.argon/CF4.or.species.eq.argon/CF4/N2) then write(*,*)prob CF4 ions - Ar,(prob_col_CFAr(i),i=3,nspec) write(*,*)prob Ar+ - CF4,prob_col_ArCF write(*,*)prob Ar+ - N2 charge exchange,prob_col_ArN2 write(*,*)prob N2+ - Ar ch ex and el col,prob_col_N2Ar write(*,*) end ifcc*/cc end of initialization end ifcc*cc ELECTRON IMPACT COLLISIONS if (nflag_col_spec(1).eq.1) then if (species.eq.argon) then call mc_e_ngas(-1,2,-1,-1,-1,-1) else if (species.eq.argon/CF4) then call mc_e_ngas(0,2,3,4,5,-1) else if (species.eq.argon/CF4/N2) then call mc_e_ngas(0,2,3,4,5,6) else if (species.eq.CF4) then call mc_e_ngas(0,-1,2,3,4,-1) end if end ifcc*/cc ION IMPACT COLLISIONS if (nt_ion_i.eq.nt_ion0.or.nt_ion_i.eq.0) then if (species.eq.argon.or.species.eq.argon/CF4.or. * species.eq.argon/CF4/N2) thencc*/cc Argon ion impact collisions with Ar, CF4 or N2 if (nflag_col_spec(2).eq.1) then call mc_Ar_Ar() end if end ifcc*/cc N2+ + N2 and N2+ + Ar charge exchange if (species.eq.argon/CF4/N2) then call mc_N2_N2() end ifcc*/cc if Ar/CF4 discharge if (species.eq.argon/CF4.or.species.eq.argon/CF4/N2) thencc*/cc CF3+ ion impact collisions with CF4 and Ar if (nflag_col_spec(3).eq.1) then call mc_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*/cc F- ion impact collisions with CF4 and Ar if (nflag_col_spec(4).eq.1) then call mc_ion_neut(4,0,dif_engy_fm,n_reac_fm, * num_col_ij_4,extra_col_ij_4,rate_col_ij_4) call mc_cf_Ar(4,2,num_col_ij_4,rate_col_ij_4) end ifcc*/cc CF3- ion impact collisions with CF4 and Ar if (nflag_col_spec(5).eq.1) then call mc_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*/cc if species are only CF4 else if (species.eq.CF4) thencc*/cc CF3+ ion impact collisions if (nflag_col_spec(2).eq.1) then call mc_ion_neut(2,0,dif_engy_cf3p,n_reac_cf3p, * num_col_ij_3,extra_col_ij_3,rate_col_ij_3) end ifcc*/cc F- ion impact collisions if (nflag_col_spec(3).eq.1) then call 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*/cc CF3- ion impact collisions if (nflag_col_spec(4).eq.1) then call 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*/ end if end ifcc*/ return endcc*/cc end of MCMC1、Electron_neutral collisions1）Elastic collisione+Are+Are+CF4e+CF4e+ N2e+ N22)Inelastic collisiona) Excitationi. The mian reactionTotalelectronic excitation e+Are+Ar* 11.5Vibrational excitation e+CF4e+CF4 0.108Vibrational excitation e+CF4e+CF4 0.168Vibrational excitation e+CF4e+CF4 0.077Electronic excitation e+CF4e+CF4* 7.54Ratational excitation e+ N2e+N2 0.02Vibrational excitation e+ N2e+N2 0.291Vibrational excitation e+ N2e+N2 0.59Vibrational excitation e+ N2e+N2 0.88Vibrational excitation e+ N2e+N2 1.17Vibrational excitation e+ N2e+N2 1.47Vibrational excitation e+ N2e+N2 1.76Vibrational excitation e+ N2e+N2 2.06Vibrational excitation e+ N2e+N2 2.35Vibrational excitation e+ N2e+N2 0.290 e+ N2e+N2 * 6.17 6.785 7.35 7.36 7.80 8.16 8.40 8.55 8.89 11.03 11.88 12.35 13.00 ii. The main formula (2.48) (2.49)Where Eth is the threshold energy of the inelastic collision, i.e.the excitation in this case . Taking into consideration that M+mM and gv we obtain ( 2.50) (2.51)Where E=mv2/2 is the electron energy before collision. The excitation process is treated as if it were an elastic collision with pre_collision velocities and V. The post_collision velocities are given by Eqs.(2.37)-(2.39) in which all vs are replaced by . b) Ionization(electron_A)Where A denotes the neutral particle with a mass , which can be assumed equal to that of the ion M .The process is represented as (2.52)Where A+ denotes an ion,e1 is the incident electron , e2 is the ejected electron , and the symbols in the parentheses denote the velocities .The energy balance equation is (2.53)Where Eth is the threshold energy of the ionization . Because of the large ion_to_electron mass ratio ,we can assume that the momentum of the incident electron is much lower than the momentum of the neutral particle,i.e. the incident electron removes an electron from the neutral ,and the neutral becomes an ion,continuing on its trajectron undisturbed . On one hand ,this assumption means that the created ion takes the velocity and direction of the neutral particle before collision ,i.e. V=V .On the other hand , the energy balance can be rewritten (2.54)The left_hand side is known as the excess energyE after ionization and we need to find an algorithm how to divide it into the scattered and ejected electrons. In general, when there is no published work on the division of the excess energy , it is divided randomly into two (2.55) (2.56)For the electron_Ar ionzation ,however,we use the expression for sampling the Eej by a random number R,(2.57)Where Einc =mv2/2 is the energy of the incident electron. The units of 0 and1 are electronvolts .Once the excess energy is divided, v is calculated from Eq.(2.54) (2.58) (2.59)The post_collision velocity of the scattered electron v is obtained by Eq(2.37) in which v is replaced by .The pre_collision velocity for the ejected electron , which does not exist in reality, is assumed (2.60)and the post_collision velocity is obtained again from Eq.(2.37) by replacing v with and v with vej . For many types of collision no data on differential cross_section is assumed not to depend to the deflection angle x .The total cross_section calculated from Eq.(2.43)gives T=4（g）.Therefore , the probabilityDefined by (2.44)becomes sinxdxd/4.This means that the post_collision velocities take random direction. The scattering with this property is called isotropie.Since there are no sufficient data for the differential cross_sections of electron_CF4 and electron_N2 collisions we assume that the scattering is isotropic .The approximation