硅各向异性腐蚀的仿真模型研究

严一峰，方玉明

(南京邮电大学 电子科学与工程学院，江苏 南京 210023)

硅各向异性腐蚀的仿真模型研究

严一峰，方玉明

(南京邮电大学 电子科学与工程学院，江苏 南京 210023)

为了在硅衬底上加工出多样的微结构，满足微机电系统的发展，硅各向异性腐蚀的仿真模型被大量研究。从衬底模型、腐蚀规则的制定、腐蚀速率、时间步长的角度出发，介绍了这些要素对几何模型和原子模型所包含多种仿真方法的模拟精度和模拟效率的影响，归纳了各种仿真方法的适用场合和特点，并为提出更优化的方案提供了方向。其中连续CA方法和MC方法揭示了腐蚀的本质，模拟精度高，三维处理能力强，被广泛应用。而原子模型算法简单，运算速度快，在对结构简单的MEMS仿真中也得到应用。所有仿真方法都可在权衡精度和效率的利弊后做出改进。

硅；各向异性腐蚀；仿真方法；模拟精度和模拟效率

因为简单的蚀刻设备和出色的能力来制造特殊的微观结构，硅的各向异性湿法化学刻蚀技术被广泛用于制造各种微电子机械系统(MEMS)。硅的各向异性湿法化学腐蚀的过程很复杂，从宏观上理解，不同晶向腐蚀速率的差异是由腐蚀液类型、晶向、温度及浓度等因素决定的，从微观上理解，考察的是不同的原子类型，差别甚至达到几百几千倍[1]。因而各种用于硅各向异性腐蚀模拟的模型被提出用于优化各向异性腐蚀过程和提高MEMS设计的效率。

近年来研究者们提出的腐蚀模型主要分为宏观上理解腐蚀过程的几何模型与微观上理解腐蚀过程的原子模型两大类。分析了两类模型中不同仿真方法模拟精度和模拟效率的影响要素，归纳了各种仿真方法的适用场合和特点。为以后研究者们在进行硅各向异性腐蚀模拟时根据需要选取合适的模型提供了参考，也为提出更优化的方案提供了方向。

1 硅各向异性腐蚀模型

几何模型和原子模型中每种仿真方法模拟效率和模拟精度的影响要素是不同的，相同要素影响的程度也不同，从这些角度去研究各向异性腐蚀模型，可以更深入的理解每种方法的特点，因而在选用与改进上都能取得更好的效果。

1.1 几何模型

几何模型不同于原子模型，衬底模型被看作一个连续的整体，不同晶面沿着其法线传播，不同的几何方法用不同的几何规则模拟这一过程。几何规则的选取决定了模拟的效果，主要的几何方法分为以下几种：

1.1.1 Wulff-Jaccodine方法

Wulff-Jaccodine方法中认为腐蚀前端面被看作一个集合的平面沿着其法线向内推进进行腐蚀，一定时间内推进的距离是晶面腐蚀速率与时间的乘积。但由于新出现平面难以准确确定，模拟结果不够精确。姜岩峰等人用腐蚀衍化的办法找到了新出现的平面，这种近似的方法不能满足对高模拟精度的要求[2]。

1.1.2 Slowness方法

Slowness方法采用速率的倒数计算腐蚀前端面点或线的轨迹来确定任意时刻的腐蚀形状。根据尖角的轨迹对尖角处插入的小的线段进行演化来确定这些线段被去除的时间，这是一个迭代的过程。由于三维空间中尖角处由多个面形成，在三维空间中应用Slowness方法的过程很复杂。Sequin用这种分析几何变化的方法进行了模拟，由于缺乏腐蚀速率数据，没能准确的模拟出腐蚀轮廓，也没能在腐蚀前端面穿透出晶片后继续进行计算[3]。

1.1.3 Eshape方法

Eshape方法中定义了一个E矢量，定义为从初始处到每个单位时间结束后两个相邻处正切交点处所构成的矢量，如图1所示。

图1 E形状中所定义的E矢量示意图

Eshape方法结合了Wulff-Jaccodine方法和Slowness方法的优点，不仅在模拟精度和适用结构上有所改进，还因为考虑的是整体效果，不是局部领域的点，可实现反向模拟。

但Eshape方法可以被许多其它模型共享，当多个模型交叉形成新图形时需要更多全局的细节计算。

1.1.4 Level Set方法

Level Set方法通过Level Set函数零水平集的计算追踪腐蚀前端曲面的运动。图2所示的是二元函数z=φ(x,y)所表示的三维曲面与xoy坐标面的交线。这个交线是随时间运动的曲线，用z=φ(x,y)的零水平集可以描述交线的运动轨迹[4]。

Level Set方法解决了上述几种几何方法与原子方法相比模拟精度不高和三维处理能力弱的问题，同时具有几何方法共有的算法简单、运算速度快的优点。然而，还有非常多的工作需要去完成来提升到目前商业元胞自动机仿真软件的性能，比如通过与运动波理论的结合也能像原子模型一样分析原子级的微观腐蚀原理。

图2 Level Set方法示例

计算结果的精度归根结底由Level Set方程的数值解法决定，即采用有限差分法时所构造差分格式的精度，分为离散化方案和自适应方案。申伟采用高阶精度的ENO和WENO格式用以离散空间导数来提高模拟精度，但这样也增加了求解过程的复杂性[5]，Radjenovi等人使用自适应方案中延伸的稀疏场方法，该方法用了一个近似的距离函数，使得它重新计算零水平集每个时间步长的邻居变得可行，不仅增加了精度还减少了计算工作量[6]。

几何模型中硅衬底的描述是宏观的，不同于原子模型中的元胞阵列，几何模型将硅衬底看作一个连续的整体。基于这种衬底制定的腐蚀规则，不能反应腐蚀发生的本质过程，模拟精度低，但算法简单，运算速度快。腐蚀规则中也只能采用宏观的面腐蚀速率，通常采用实验测定后进行插值计算的方法获取关键方向上的腐蚀速率，比如测定硅半球在给定腐蚀条件下前后的外形变化，一个关键晶向需要测定不同温度时的腐蚀速率来获取Arrhenius定律中所需参数，完整腐蚀速率库可在应用Cubic(立方插值法)、v4(双谐样条插值法)、Hubbard等不同的插值方法后获取[5,7]。

1.2 原子模型

在原子方法中，硅衬底是晶格上所有原子组成的阵列，通过在腐蚀时间内不断判断腐蚀前端面原子的去留完成仿真。硅各向异性腐蚀模拟的MC方法与CA方法是分别基于蒙特卡罗方法和元胞自动机思想的。

1.2.1 CA方法

元胞自动机思想中元胞的状态在时间上的进化是元胞周围环境作用的结果。Than和Buttgenbach等人提出的随机CA方法将同一晶面的原子看作一次随机过程，因此腐蚀表面模糊不清，而且阈值的确定方式对模拟结果的影响很大。Zhu和Liu等人提出的连续CA方法通过腐蚀速率的计算去除原子的“生命值”，获得了高模拟精度，但由于每个腐蚀步长计算的对象是三维阵列中的所有元胞，降低了模拟效率。Than和Buttgenbach等人提出的动态CA方法在计算机仿真时，只把腐蚀前端面的原子作为处理对象，,程序的简化实现了模拟效率的大幅提高，因此动态CA方法被广泛应用。

碰撞理论与CA方法结合计算腐蚀速率是最早从腐蚀原理上入手的计算方法，该理论认为腐蚀液的温度以及活化粒子的浓度决定了OH根离子与硅原子发生有效碰撞的几率，描述了腐蚀液浓度和温度等对腐蚀速率的影响，局限性在于只能定义出(100)，(110)，(111)3个主要晶面的腐蚀速率。跃迁态理论同样可以融合进连续CA方法解释硅各向异性腐蚀，不同于前者，计算的是微观的原子腐蚀速率，如图3所示。表面类型原子根据不同配置可分为14种关键类型，不同类型原子的微观激活能不同，即转变态能量减去初始态能量的不同，这是由OH根离子对不同配置原子总化学键的削弱不同造成的，腐蚀速率可用波尔茨曼方程求解[8-10]。

图3 一步反应的微观激活能

高米勒指数晶面上原子种类大多数情况下超过两个，如果用实验测定的宏观腐蚀速率，连续CA方法只能在有限范围内精确测定。跃迁态理论与连续CA方法的结合，成功把宏观的面腐蚀速率转化为微观的原子腐蚀速率，加速了高指数晶面上新颖MEMS器件的研究。台阶流动模型曾被引入连续CA方法解决腐蚀速率不同的解释问题，极大增强了连续CA方法的性能[11]。

动态CA方法相比较于连续CA方法只需考虑腐蚀表面原子，可采用更高效的衬底存储方式进一步提高模拟效率。如图4所示，JPD等人给出了(100)晶面<110>晶向中表面晶胞4个不等价层中晶体中原子和自动机元胞的对应关系[12]。动静态结合存储中，一个硅晶胞中的原子被抽象成4层原子，表面的4层原子可循环表示整个硅衬底，这种存储方式在原子的存储数量以及单个原子所需存储的数据量都要优于静态存储。以(100)晶面1 000×1 000晶胞的晶片为例，一个晶胞平均一层有两个硅原子，模拟腐蚀深度为500层原子，静态存储存储量为40×1 000×1 000×2×500 Byte，约为40 GB,而动静态结合存储所需的存储量为4×1 000×1 000×2×4 Byte，约为32 MB[13]。CA方法在大尺寸模拟时，一个元胞如果代表一个原子的话，过多的原子数目都会对计算机造成巨大压力，掩膜图形分辨率决定了元胞代表的原子数量，这样虽然提高了模拟效率，却降低了模拟精度[14]。

图4 (100)晶面<110>晶向上四个不等价层的单元晶胞

连续CA方法的模拟精度和模拟效率也受到时间步长的影响。Zhu Z等人根据连续CA方法的特点提出了时间补偿的概念,将剩余时间补偿到下一时间步长,模拟精度得到保证，但仿真效率会降低一个数量级[15]。陈劲源等人针对传统动态CA方法中串行步骤过多的问题，在腐蚀步骤直接并行的基础上提出了超级步长的概念，腐蚀算法中串行步骤的工作被有效减少，实现了十几倍的加速[16]。

1.2.2 MC方法

蒙特卡罗方法是解决随机性问题的，在硅各向异性腐蚀模拟中的MC方法，把原子分为固定和脱落两种状态，MC方法预先确定原子间键断裂的概率，而原子间键则决定原子是否被移除，表面的所有原子都经过测试，以检验它们是否能从晶体中除去。与CA方法相比，两者的三维处理能力都很强，但MC方法在小尺度仿真时能发挥更大的作用，模拟精度更高。

跃迁态理论融合进MC法中时，不同类型原子在不同概率的情况下发生能级跃迁。Gosalves等人曾在上述模型基础上通过计算测试得到与实验吻合较好的概率方程，函数本身取决于参数，如第一和第二相邻原子的数目，第一邻居和第二邻居之间的平均相互作用能。MC方法与跃迁态理论的结合计算出了主要类型原子的腐蚀概率，相对于由几个基本晶面速率确定的腐蚀规则，模拟精度从根本上得到提高[12,17-18]。台阶流动模型也能融入MC方法，如图5所示。基于台阶流动模型的刻蚀速度可以被看作由台地(terrace)的速度和与其垂直的台阶(step)的速度组成，与台地原子相比，台阶处的原子刻蚀速度更高，因此表面的刻蚀基本只移除台阶处的高腐蚀速度原子，与此同时台地的刻蚀在缓慢的进行，形成了不断产生的扭折，不同类型原子在台地和台阶上所占的比例不同，造成了腐蚀速率的差异。该模型的特征是通过原子的离散化系统提供了了解系统演化的可能性。(110)面锯齿状条纹表面产生的过程，(100)表面金字塔状的凸起结构，证明这种模型解释各向异性腐蚀是有效的[19]。与CA方法相同的是，从原子级的微观腐蚀原理得出的腐蚀概率能在更大的仿真范围取得精确的结果，但MC方法的计算效率受到整体上表面原子配置转移概率大小的影响，如(111)面中，计算会变的非常缓慢。

图5 台阶结构

MC方法在进行衬底存储时，是将单元内的原子

定义为实际原子的。在大尺寸模拟时，将计算得到的结果选择一个合适的比例映射到实际尺寸和时间规模上。这样做的优势是大尺寸模拟时可以达到与CA方法一样的存储效果，而小尺寸模拟时可以最大化MC方法小尺度仿真的能力。MC方法的腐蚀步长有确定的算式，与CA方法一样提升了精度，但不会牺牲模拟效率。

原子遍历机制也决定着MC方法拥有更高的模拟精度。MC方法可以研究表面细微结构和粗糙度等特征，尤其突出于研究腐蚀过程中腐蚀液不同局部环境以及金属离子掺杂等影响下的接触面的微观腐蚀形态，这是因为CA法并行删除原子造成的，MC法中同一晶面的原子会在不同时间被腐蚀掉[20]。Peng等人分析了硅在HF和双氧水中被腐蚀时金属粒子的运动行为，(100)面的金属粒子通过原位双极电化学反应产生的能量转化为推动其运动的机械能，产生了垂直于(100)面的纳米线和纳米孔，这种微观腐蚀形态的展现采用MC方法更好[21]；而在系统结构的初步设计 中，CA方法能发挥更大的作用。总的来说，原子模型揭示了各向异性腐蚀过程的本质，注重的是局部原子之间互相作用，相比较于几何模型，模拟精度高，三维处理能力强。

2 不同仿真方法的归纳

原子模型和几何模型所包含的仿真方法在硅各向异性腐蚀模拟中应用的归纳如表1所示。

表1 原子模型和几何模型所包含的仿真方法在硅各向异性腐蚀模拟中应用的归纳

3 结束语

原子模型由于三维处理能力强，模拟精度高等优点受到青睐，几何模型由于算法简单，模拟效率高的优点在简单的二维结构模拟中也得到应用。原子模型中连续CA方法适用于系统结构的初步设计，而MC方法在小尺度仿真时能发挥更大作用；几何模型中，Level Set方法和Eshape方法具有一定的三维处理能力，其中Level Set方法通过数值算法、腐蚀速率求解中插值方法的改进能达到原子模型的精度，同时具备几何模型高效的特点，有待进一步研究。模拟精度和模拟效率是选取硅各向异性腐蚀模型的两大因素，可以从衬底从衬底模型、腐蚀规则的制定、腐蚀速率、时间步长的角度权衡这两方面的利弊，对腐蚀模型进行改进。

[1] 余丹铭,梁利华,许杨剑.微电子机械技术的研究和发展趋势[J].电子机械工程,2005(1):5-9.

[2] Sequin C H.Computer simulation of anisotropic crystal etching[C].CA,USA:International Conference on Solid- State Sensors and Actuators,1991.

[3] 姜岩峰,黄庆安.硅各向异性腐蚀三维计算机模拟系统的建立[J].机械强度,2001,23(4):484- 487.

[4] 杨思燕.Level Set方法在图像处理中的应用[J]. 电子科技,2015,28(7):153-156.

[5] 申伟.基于Level Set的硅湿法刻蚀模拟研究[D]. 南京:东南大学,2012.

[6] Radjenovi B,Radmilovi-Radjenovi M,Mitri M.Level set approach to anisotropic wet etching of silicon[J].Sen- sors,2010,10(5):4950-4967.

[7] 张佩君,黄庆安.硅各向异性腐蚀速率图的模拟[J]. 固体电子学研究与进展,2002,22(1):89-92.

[8] Gosalvez M A,Nieminen R M,Kilpinen P,et al. Anisotropic wet chemical etching of crystalline silicon: atomistic Monte-Carlo simulations and experiments[J]. Applied Surface Science,2001,178(1):7-26.

[9] Gosálvez M A,Foster A S,Nieminen R M.Atomistic simulations of surface coverage effects in anisotropic wet chemical etching of crystalline silicon[J].Applied Surface Science,2002,202(3):160-182.

[10] Zhou Z F,Huang Q A,Li W H.An atomic level model for silicon anisotropic etching processes: Cellular automaton simulation and experimental verification[J]. Applied Physics Letters,2007,91(17):174101- 174103.

[11] Rezvankhah M A,Merati A R.Step flow model in continuous cellular automata method for simulation of anisotropic etching of silicon[J].Journal of Micro/ Nano- lithography Mems & Moems,2013,12(2):23004-23009.

[12] Oliveira J P D,Carreo M N.Simulation of aniso- tropic etching of silicon using a celluar automata model [J].ECS Transactions,2008,14 (1):89-93.

[13] 尹宝林,韩秋菊,康建初,等.单晶硅各向异性腐蚀模拟中的存储结构和算法[C].天津:全国微米/纳米技术学术会议,2000.

[14] 周丽玲.基于CA模型的单晶硅各向异性腐蚀微观模拟[D].大连:大连理工大学,2006.

[15] Zhu Z,Liu C.Simulation of anisotropic crystallineetching using a continuous cellular automata algorithm [J].Journal of Computer Modeling in Engineering and Science,2000,1(1):11-17.

[16] 陈劲源,李建华,郭卫斌.改进的硅各向异性腐蚀GPU并行模拟[J].计算机应用,2013,33(12):3317-3320.

[17] 朱鹏,幸研,易红.基于原子模型的各向异性腐蚀模拟方法对比分析[J].传感技术学报,2007,20(12): 2540-2544.

[18] Gosálvez M A,Foster A S,Nieminen R M. Multiscale modeling of anisotropic wet chemical etching of crystalline silicon[J].EPL,2007,60(3):467-473.

[19] 朱鹏,幸研,易红,等. Metropolis蒙特卡罗模拟湿法腐蚀算法研究及应用[J].半导体学报,2008,29(1):183-188.

[20] 朱鹏.基于蒙特卡罗的各向异性湿法腐蚀模拟研究[D].南京:东南大学,2008.

[21] Peng K,Lu A,Zhang R,et al.Motility of metal na- noparticles in silicon and induced anisotropic silicon etching[J].Advanced Functional Materials,2010,18(19):3026- 3035.

Study on Simulation Model of Anisotropic Etching of Silicon

YAN Yifeng,FANG Yuming

(School of Electronic Science and Engineering,Nanjing University of Posts and Telecommunications,Nanjing 210023,China)

In order to produce various microstructures on silicon substrate and satisfy the development of MEMS, the simulation model of silicon anisotropic etching has been studied extensively. From the perspective of substrate model, formulation of corrosion rules, corrosion rate, time step,introduces the influence to simulation accuracy and efficiency of simulation of various simulation methods contained in the geometric model and atomic model from these elements,summarizes the application and characteristics of various simulation methods, and provides the direction for putting forward a more optimal solution . Among them, the continuous CA method and the MC method reveal the essence of corrosion. It has high simulation accuracy and strong three-dimensional processing ability, and has been widely used. The atom model algorithm is simple and fast, and it is also applied to the simulation of MEMS with simple structure. All simulation methods can be improved after weighing the pros and cons of accuracy and efficiency.

silicon;anisotropic etching;simulation method;simulation accuracy and simulation efficiency

2017- 06- 25

江苏省自然科学基金(BK20131380);南京邮电大学大学生创新训练计划项目(XYB2017057)

严一峰(1992-)，男，硕士研究生。研究方向：微机电系统设计。

TN773

A

1007-7820(2018)02-056-05