Analysis of Anisotropic Etching Characteristics and Morphology Simulation of Crystal Silicon

doi:10.16553/j.cnki.issn1000-985x.2025.0211

Abstract

Abstract: The anisotropic wet etching characteristics of crystal silicon are complex and easily affected by etching conditions and mask shapes，which makes it difficult to accurately predict and control the evolution process and morphology of etching structural. This study was based on experimental data on the all crystal plane etching rate and mask etching structure of crystal silicon，and analyzed in detail the anisotropic etching mechanism and mask etching forming process. A simulation model of crystal silicon etching morphology （Si-LST） was constructed using the Level-Set interpolation method，which achieved accurate interpolation of the all crystal plane etching rate using a small amount of crystal plane etching rate and accurate simulation of etching morphology under any mask. The results show that the Si-LST has high simulation accuracy for concave，convex and composite masks，which can provide efficient process design assistance for the processing and surface quality control of crystal silicon micro-structures.

Key words: crystal silicon; wet etching; anisotropic; etching mechanism; morphology simulation

CLC Number:

O793

ZHANG Hui, QIAN Jun. Analysis of Anisotropic Etching Characteristics and Morphology Simulation of Crystal Silicon[J]. Journal of Synthetic Crystals, 2026, 55(2): 241-252.

Figures/Tables 16

Fig.1 Experiment results of wet etching crystal silicon hemisphere. （a） Etching rate top view distribution of crystal planes on crystal silicon hemisphere （red is high etching rate crystal plane，blue is low etching rate crystal plane）；（b） position distribution of extreme rate crystal planes on crystal silicon hemisphere；（c） etching rate curve of all crystal silicon crystal planes and etching surface morphology of main crystal planes

Fig.2 Evolution process of forming sidewall etching at mask boundary

Fig.3 Evolution process of crystal planes at mask boundary intersection

Fig.4 Forming process of etching structural surface on sidewall of Si（100） concave window mask. （a） Shape and etching structure of concave window mask；（b） etching progress of sidewall structure plane at different stages

Fig.5 Forming process of convex corner structure of Si（100） convex mask. （a） Shape and etching structure of convex mask；（b） convex corner etching structure at initial etching stage；（c） dynamic adjustment process of convex mask boundary

Fig.6 Crystal silicon unit cell （a），symmetry characteristics of different crystal orientations （b）~（d），and three typical silicon crystal plane structures （e）~（g）

Fig.7 Surface atomic structure types on Si （110） crystal plane

Table 1 Types and proportions of surface atomic structures on crystal planes with different etching rates

Etching rate	Crystal plane	Si原子配位类型：P%
Low	$(111)$	（0，3，9）：100.0%
Middle	$(100)$	（0，2，6）：100.0%
Middle	$(221)$	（0，3，9）：50.0% （2，1，7）：50.0%
High	（311）	（1，1，5）：30.0% （1，2，7）：30.0% （0，3，9）：30.0%
High	$(110)$	（2，1，7）：100.0%

Table 1 Types and proportions of surface atomic structures on crystal planes with different etching rates

Etching rate	Crystal plane	Si原子配位类型：P%
Low	$(111)$	（0，3，9）：100.0%
Middle	$(100)$	（0，2，6）：100.0%
Middle	$(221)$	（0，3，9）：50.0% （2，1，7）：50.0%
High	（311）	（1，1，5）：30.0% （1，2，7）：30.0% （0，3，9）：30.0%
High	$(110)$	（2，1，7）：100.0%

Fig.8 Step etching flow model of crystal silicon and schematic diagram of etching Si（111） crystal plane

Fig.9 Non-orthogonal three-dimensional coordinate system （a，b，c） constructed based on four crystal plane groups ｛111｝，｛110｝，｛100｝，｛311｝

Table 2 Transformation relationships of coordinate systems in six non-orthogonal three-dimensional coordinate regions

Coordinate matrix	Scaling matrix	Coordinate transformation	Normalization factor
$100110311$	$100010003$	a=x-y-2z b=y-z c=3z	1/x
$100101311$	$100010003$	a=x-2y-z b=z-y c=3y	1/x
$111110311$	$100010003$	a=（y-x）/2+z b=y-z c=3/2·（x-y）	1/x
$111101311$	$100010003$	a=（2y-x+z）/2 b=z-y c=3/2·（x-z）	1/x
$111101113$	$100010003$	a=（x+2y-z）/2 b=x-y c=3/2·（z-x）	1/z
$001101113$	$100010003$	a=-x-2y+z b=x-y c=3y	1/z

Table 2 Transformation relationships of coordinate systems in six non-orthogonal three-dimensional coordinate regions

Coordinate matrix	Scaling matrix	Coordinate transformation	Normalization factor
$100110311$	$100010003$	a=x-y-2z b=y-z c=3z	1/x
$100101311$	$100010003$	a=x-2y-z b=z-y c=3y	1/x
$111110311$	$100010003$	a=（y-x）/2+z b=y-z c=3/2·（x-y）	1/x
$111101311$	$100010003$	a=（2y-x+z）/2 b=z-y c=3/2·（x-z）	1/x
$111101113$	$100010003$	a=（x+2y-z）/2 b=x-y c=3/2·（z-x）	1/z
$001101113$	$100010003$	a=-x-2y+z b=x-y c=3y	1/z

Table 3 Interpolation calculation formula of etching rate in six non-orthogonal three-dimensional coordinate regions

Unit vector

Sphere scope

Interpolation rate calculation formula

［100］［110］［311］

$0 ≤ θ ≤ a r c t a n (1 / 3)$

$0 ≤ w ≤ a r c t a n (s i n θ)$

$R (x, y, z) = (1 / x) (x - y - 2 z) R 100 + (y - z) R 110 + 3 z R 311$

$R θ, w = 1 - t a n θ - 2 t a n w c o s θ R 100 + t a n θ - t a n w c o s θ R 110 + 3 t a n w c o s θ R 311$

$a r c t a n (1 / 3) ≤ θ ≤ (π / 4)$

$0 ≤ w ≤ a r c t a n (c o s θ - s i n θ) / 2$

［100］［101］［311］

$0 ≤ θ ≤ a r c t a n (1 / 3)$

$w ≥ a r c t a n (s i n θ)$

$w ≤ a r c t a n (c o s θ - 2 s i n θ)$

$R (x, y, z) = (1 / x) (x - 2 y - z) R 100 + (z - y) R 110 + 3 y R 311$

$R θ, w = 1 - 2 t a n θ - t a n w c o s θ R 100 + t a n w c o s θ - t a n θ R 110 + 3 t a n w c o s θ R 311$

［111］［110］［311］

$a r c t a n (1 / 3) ≤ θ ≤ (π / 4)$

$w ≥ a r c t a n (c o s θ - s i n θ) / 2$

$w ≤ a r c t a n (s i n θ)$

$R (x, y, z) = 1 x {[12 (y - x) + z] R 111 + (y - z) R 110 + 32 (x - y) R 311}$

$R θ, w = 12 t a n θ + t a n w c o s θ - 12 R 111 + t a n θ - t a n w c o s θ R 110 + 32 1 - t a n θ R 311$

［111］［101］［311］

$0 ≤ θ ≤ a r c t a n (1 / 3)$

$w ≥ a r c t a n (c o s θ - 2 s i n θ)$

$w ≤ a r c t a n (c o s θ)$

$R (x, y, z) = 1 x [12 (2 y - x + z) R 111 + (z - y) R 110 + 32 (x - z) R 311]$

$R θ, w = 12 (2 t a n θ + t a n w c o s θ - 1) R 111 + t a n w c o s θ - t a n θ R 110 + 32 1 - t a n w c o s θ R 311$

$a r c t a n (1 / 3) ≤ θ ≤ (π / 4)$

$w ≥ a r c t a n (s i n θ)$

$w ≤ a r c t a n (c o s θ)$

［111］［101］［113］

$0 ≤ θ ≤ (π / 4)$

$w ≥ a r c t a n (c o s θ)$

$w ≤ a r c t a n (c o s θ + 2 s i n θ)$

$R (x, y, z) = 1 z [12 (x + 2 y - z) R 111 + (x - y) R 110 + 32 (z - x) R 311]$

$R θ, w = 12 (c o s θ t a n w + 2 s i n θ t a n w - 1) R 111 + c o s θ t a n w - s i n θ t a n w R 110 + 32 1 - c o s θ t a n w R 311$

［001］［101］［113］

$0 ≤ θ ≤ (π / 4)$

$w ≥ a r c t a n (c o s θ + 2 s i n θ)$

$w ≤ π / 2$

$R (x, y, z) = 1 z [(- x - 2 y + z) R 100 + (x - y) R 110 + 3 y R 311]$

$R θ, w = - c o s θ t a n w + 2 s i n θ t a n w + 1 R 100 + c o s θ t a n w - s i n θ t a n w R 110 + 3 s i n θ t a n w R 311$

Table 3 Interpolation calculation formula of etching rate in six non-orthogonal three-dimensional coordinate regions

Unit vector

Sphere scope

Interpolation rate calculation formula

［100］［110］［311］

$0 ≤ θ ≤ a r c t a n (1 / 3)$

$0 ≤ w ≤ a r c t a n (s i n θ)$

$R (x, y, z) = (1 / x) (x - y - 2 z) R 100 + (y - z) R 110 + 3 z R 311$

$R θ, w = 1 - t a n θ - 2 t a n w c o s θ R 100 + t a n θ - t a n w c o s θ R 110 + 3 t a n w c o s θ R 311$

$a r c t a n (1 / 3) ≤ θ ≤ (π / 4)$

$0 ≤ w ≤ a r c t a n (c o s θ - s i n θ) / 2$

［100］［101］［311］

$0 ≤ θ ≤ a r c t a n (1 / 3)$

$w ≥ a r c t a n (s i n θ)$

$w ≤ a r c t a n (c o s θ - 2 s i n θ)$

$R (x, y, z) = (1 / x) (x - 2 y - z) R 100 + (z - y) R 110 + 3 y R 311$

$R θ, w = 1 - 2 t a n θ - t a n w c o s θ R 100 + t a n w c o s θ - t a n θ R 110 + 3 t a n w c o s θ R 311$

［111］［110］［311］

$a r c t a n (1 / 3) ≤ θ ≤ (π / 4)$

$w ≥ a r c t a n (c o s θ - s i n θ) / 2$

$w ≤ a r c t a n (s i n θ)$

$R (x, y, z) = 1 x {[12 (y - x) + z] R 111 + (y - z) R 110 + 32 (x - y) R 311}$

$R θ, w = 12 t a n θ + t a n w c o s θ - 12 R 111 + t a n θ - t a n w c o s θ R 110 + 32 1 - t a n θ R 311$

［111］［101］［311］

$0 ≤ θ ≤ a r c t a n (1 / 3)$

$w ≥ a r c t a n (c o s θ - 2 s i n θ)$

$w ≤ a r c t a n (c o s θ)$

$R (x, y, z) = 1 x [12 (2 y - x + z) R 111 + (z - y) R 110 + 32 (x - z) R 311]$

$R θ, w = 12 (2 t a n θ + t a n w c o s θ - 1) R 111 + t a n w c o s θ - t a n θ R 110 + 32 1 - t a n w c o s θ R 311$

$a r c t a n (1 / 3) ≤ θ ≤ (π / 4)$

$w ≥ a r c t a n (s i n θ)$

$w ≤ a r c t a n (c o s θ)$

［111］［101］［113］

$0 ≤ θ ≤ (π / 4)$

$w ≥ a r c t a n (c o s θ)$

$w ≤ a r c t a n (c o s θ + 2 s i n θ)$

$R (x, y, z) = 1 z [12 (x + 2 y - z) R 111 + (x - y) R 110 + 32 (z - x) R 311]$

$R θ, w = 12 (c o s θ t a n w + 2 s i n θ t a n w - 1) R 111 + c o s θ t a n w - s i n θ t a n w R 110 + 32 1 - c o s θ t a n w R 311$

［001］［101］［113］

$0 ≤ θ ≤ (π / 4)$

$w ≥ a r c t a n (c o s θ + 2 s i n θ)$

$w ≤ π / 2$

$R (x, y, z) = 1 z [(- x - 2 y + z) R 100 + (x - y) R 110 + 3 y R 311]$

$R θ, w = - c o s θ t a n w + 2 s i n θ t a n w + 1 R 100 + c o s θ t a n w - s i n θ t a n w R 110 + 3 s i n θ t a n w R 311$

Fig.10 Interpolation results of etching rate of all crystal silicon crystal planes under etching condition of 71 ℃，40%KOH. （a） Top view etching rate in <111> crystal orientation family；（b） top view etching rate in <100> crystal orientation family；（c） comparison of interpolation results with experimental data

Table 4 Etching rates of thirteen typical crystal planes

Crystal plane	（100）	（110）	（210）	（211）	（221）	（310）	（311）	（320）	（331）	（530）	（540）	（111）	（411）
Etching rate at 40%KOH/（μm·min^-1）	0.60	1.25	1.26	0.50	0.54	0.95	0.90	1.10	0.80	1.16	1.18	0.01	0.76

Fig.11 Etching rates of all crystal plane obtained by interpolation of etching rates of thirteen crystal planes. （a） Top view etching rate in <111> crystal orientation family；（b） top view etching rate in <100> crystal orientation family；（c） comparison of interpolation results with experimental data

Fig.12 Comparison of experimental results （a）~（c） and simulation results （d）~（f） of mask etching

References 17

[1]	张　辉，钱　珺，洪莉莉. 蓝宝石晶体湿法刻蚀各向异性研究与机理分析［J］. 人工晶体学报，2023，52（6）：1128-1135.
	ZHANG H，QIAN J，HONG L L. Anisotropy and mechanism of wet etching of sapphire crystal［J］. Journal of Synthetic Crystals，2023，52（6）：1128-1135 （in Chinese）.
[2]	张　辉. 石英各向异性湿法刻蚀机理及工艺模型研究［D］. 南京：东南大学，2018.
	ZHANG H. The etching mechanism and process model research for quartz anisotropic wet etching［D］. Nanjing：Southeast University，2018 （in Chinese）.
[3]	蔡鹏鹏. 石英各向异性湿法刻蚀特性及模拟研究［D］. 南京：东南大学，2016.
	CAI P P. The character research and simulation of anisotropic wet etching quartz［D］. Nanjing：Southeast University，2016 （in Chinese）.
[4]	GROSSE A，GREWE M，FOUCKHARDT H. Deep wet etching of fused silica glass for hollow capillary optical leaky waveguides in microfluidic devices［J］. Journal of Micromechanics and Microengineering，2001，11（3）：257-262.
[5]	CHEN Y C，LIN B W，HSU W C，et al. The formation of smooth facets on wet-etched patterned sapphire substrate［J］. ECS Journal of Solid State Science and Technology，2014，3（2）：R5-R8.
[6]	GAO Z D，WANG Q J，ZHANG Y，et al. Etching study of poled lithium tantalate crystal using wet etching technique with ultrasonic assistance［J］. Optical Materials，2008，30（6）：847-850.
[7]	GOSÁLVEZ M A，CHENG D，NIEMINEN R M，et al. Apparent activation energy during surface evolution by step formation and flow［J］. New Journal of Physics，2006，8（11）：269.
[8]	WANG J，GUO L W，JIA H Q，et al. Fabrication of patterned sapphire substrate by wet chemical etching for maskless lateral overgrowth of GaN［J］. Journal of the Electrochemical Society，2006，153（3）：C182.
[9]	GOSÁLVEZ M A，FERRANDO N，XING Y，et al. Simulating anisotropic etching of silicon in any etchant：evolutionary algorithm for the calibration of the continuous cellular automaton［J］. Journal of Micromechanics and Microengineering，2011，21（6）：065017.
[10]	FERRANDO N，GOSÁLVEZ M A，COLÓM R J. Evolutionary continuous cellular automaton for the simulation of wet etching of quartz［J］. Journal of Micromechanics and Microengineering，2012，22（2）：025021.
[11]	ZHANG H，XING Y，GOSÁLVEZ M A，et al. Removal probability function for kinetic Monte Carlo simulations of anisotropic etching of silicon in alkaline etchants containing additives［J］. Sensors and Actuators A：Physical，2015，233：451-459.
[12]	谢海波，傅　新，刘　玲，等. 基于玻璃湿法刻蚀的微流控器件加工工艺研究［J］. 机械工程学报，2003，39（11）：123-129.
	XIE H B，FU X，LIU L，et al. Research on glass wet etching for micro fluidic devices［J］. Chinese Journal of Mechanical Engineering，2003，39（11）：123-129 （in Chinese）.
[13]	AGOS LVEZ M，NIEMINEN R M. Surface morphology during anisotropic wet chemical etching of crystalline silicon［J］. New Journal of Physics，2003，5：100.
[14]	GOSÁLVEZ M A. Atomistic modeling of anisotropic etching of crystalline silicon［D］. Finland：Laboratory of Physics Helsinki University of Technology，2003.
[15]	XING Y，GOSÁLVEZ M A，SATO K，et al. Evolutionary determination of kinetic Monte Carlo rates for the simulation of evolving surfaces in anisotropic etching of silicon［J］. Journal of Micromechanics and Microengineering，2012，22（8）：085020.
[16]	RADJENOVIĆ B，RADMILOVIĆ-RADJENOVIĆ M. 3D simulations of the profile evolution during anisotropic wet etching of silicon［J］. Thin Solid Films，2009，517（14）：4233-4237.
[17]	王振海. Level Set算法研究及其在拓扑优化设计中的应用［D］. 西安：西北工业大学，2006.
	WANG Z H. Study on Level Set methods and its applications in topology optimization［D］. Xi’an：Northwestern Polytechnical University，2006 （in Chinese）.