Optimization of the Monte Carlo Simulation for Sapphire in Wet Etching

Abstract

In this paper, the Monte Carlo simulation for sapphire in wet etching is optimized, which improves the accuracy and efficiency of simulated results. Firstly, an eight-index classification method is proposed to classify the kinds of surface atoms, which can make assigned removal probabilities more accurately for surface atoms. Secondly, based on the proposed classification method of surface atoms, an extended removal probability equation (E-RPE) is proposed, which makes the errors between simulated and experimental rates smaller and greatly improves the accuracy of the simulated result of the etch rate distribution under the experimental condition (H₂SO₄:H₃PO₄ = 3:1, 236 °C). Thirdly, a modified removal probability equation (ME-RPE) considering the temperature dependence is proposed based on the error analysis between the simulated and experimental rates under different temperature conditions, which can simulate etch rates under the different temperature conditions through a group of optimized energy parameters and improve the simulation efficiency. Finally, small errors between the simulated and experimental rates under the different temperature conditions (H₂SO₄:H₃PO₄ = 3:1, 202 °C and 223 °C) verify the validity of the ME-RPE for temperature change. The optimization methods for the Monte Carlo simulation of sapphire in wet etching proposed in this paper will provide a reference for the simulation of other crystal materials.

1. Introduction

Sapphire is widely used in microelectromechanical systems (MEMS) because of its good light transmittance, insulation and chemical stability [1,2,3,4,5,6,7,8,9]. Wet etching is the preferred method for preparing patterned sapphire substrates. Etchants are commonly mixed solutions of H₂SO₄ and H₃PO₄ with various volume ratios, and the etching temperature is commonly above 200 °C [10,11,12,13,14,15,16,17]. The higher the temperature is, the more intense the etching reaction is. Considering the stability of the etching reaction and the requirement for etching devices, the etching temperature is usually in the range of 200–300 °C [11,13]. The simulation for sapphire in wet etching can improve wet etching technology and reduce the development cost of sapphire microstructures [18,19].

The Monte Carlo method is an atomic simulation method, which can effectively simulate etch rates and etch structures [18,20,21,22,23,24]. According to the simplified atomic structure, the literature [18] proposed a six indices classification method ($n_{1s}^G,n_{1b}^G,n_{1s}^R,n_{1b}^R,n_{2s},n_{2b}$) to classify the kinds of surface atoms. The simplification of the atomic structure of sapphire is shown in Figure 1; the co-edge connection and co-plane connection among octahedrons are simplified as the grey-bond connection and red-bond connection among Al atoms. The six indices in the six indices classification method ($n_{1s}^G,n_{1b}^G,n_{1s}^R,n_{1b}^R,n_{2s},n_{2b}$) represent the amounts of first-grey-bond surface atoms, first-grey-bond bulk atoms, first-red-bond surface atoms, first-red-bond bulk atoms, second surface atoms and second bulk atoms, respectively. The proposed removal probability equation (RPE) of surface atoms is as follows:

$$P\left(n_{1s}^G,n_{1b}^G,n_{1s}^R,n_{1b}^R,n_{2s},n_{2b}\right)=P_0\times P_1^G\left(n_{1s}^G,n_{1b}^G\right)\times P_1^R\left(n_{1s}^R,n_{1b}^R\right)\times P_2\left(n_{2s},n_{2b}\right)$$

(1)

Among them,

$$P_0=\left(1+e^{-\beta E_1}\right)\times \left(1+e^{-\beta E_2}\right)\times \left(1+e^{-\beta E_3}\right)$$

(2)

$$P_1^G\left(n_{1s}^G,n_{1b}^G\right)=1/(1+e^{\beta \times F_1^G\left(n_{1s}^G,n_{1b}^G\right)}) \mathrm{and} F_1^G\left(n_{1s}^G,n_{1b}^G\right)={\epsilon }_1\times n_{1s}^G+{\epsilon }_2\times n_{1b}^G-E_1$$

(3)

$$P_1^R\left(n_{1s}^R,n_{1b}^R\right)=1/(1+e^{\beta \times F_1^R\left(n_{1s}^R,n_{1b}^R\right)}) \mathrm{and} F_1^R\left(n_{1s}^R,n_{1b}^R\right)={\epsilon }_3\times n_{1s}^R+{\epsilon }_4\times n_{1b}^R-E_2$$

(4)

$$P_2\left(n_{2s},n_{2b}\right)=1/(1+e^{\beta \times F_2\left(n_{2s},n_{2b}\right)}) \mathrm{and} F_2\left(n_{2s},n_{2b}\right)={\epsilon }_5\times n_{2s}+{\epsilon }_6\times n_{2b}-E_3$$

(5)

Here, $\beta =K_{B}T$, parameter $K_B$ is the Boltzmann constant, and parameter $T$ is the absolute temperature (K). Energy parameters ${\epsilon }_1$ (${\epsilon }_3$) and ${\epsilon }_2$ (${\epsilon }_4$) are average energies between target atoms and their first-grey (red)-bond surface atoms, as well as first-grey (red)-bond bulk atoms, respectively; energy parameters ${\epsilon }_5$ and ${\epsilon }_6$ are average energies between target atoms and their second surface atoms and second bulk atoms, respectively; energy parameters $E_1$, $E_2$ and $E_3$ are critical energies between target atoms and their first-grey-bond atoms, first-red-bond atoms and second atoms, respectively.

Based on the proposed six indices classification method and the corresponding RPE, the literature [18] performed the Monte Carlo simulation for sapphire in wet etching under the experimental condition (236 °C, H₂SO₄:H₃PO₄ = 3:1). The simulated results show that there are still errors in the simulation for the etch rate distribution, mainly in that the maximum rate zone fails to split into two zones. On this basis, the Monte Carlo simulation of sapphire will be optimized to improve the accuracy of the simulation. Furthermore, based on the role of the temperature variable in the removal probability equation (the higher the temperature is, the higher the removal probabilities of surface atoms are, and the higher the etch rates are), a modified removal probability equation (ME-RPE) is proposed by analyzing the errors between simulated and experimental rates under different temperature conditions (H₂SO₄:H₃PO₄ = 3:1, 202–304 °C). Thus, etch rates under different temperature conditions can be accurately simulated according to a group of energy parameters, and the simulation efficiency can be improved.

Therefore, the rest of the paper can be divided into the following parts. In the second part, etch rates at different temperatures are obtained through etching experiments of sapphire hemispheres. In the third part, an eight indices classification method is proposed, which can classify kinds of surface atoms in more detail, so that the removal probabilities of surface atoms can be assigned more accurately. By comparing the simulated and experimental results of the etch rate distribution, the effectiveness of the eight indices classification method of surface atoms is verified. In the fourth part, the errors between simulated rates and experimental rates under different temperature conditions are compared, and a ME-RPE is proposed, so that the etch rates under different temperature conditions can be accurately simulated according to a set of energy parameters, and the simulation efficiency can be improved. The validity of the ME-RPE is verified by comparing simulated and experimental rates at different temperatures. In the fifth part, this article is summarized.

2. Experimental

Etch rate distributions under different temperature conditions can be obtained through etching experiments of sapphire hemispheres [18,25,26,27]. The radius of the sapphire hemisphere selected in this paper is 21.50 mm, and its upper plane toward the <0 0 0 1> direction. The mixture of H₂SO₄ (98 wt%) and H₃PO₄ (87 wt%) in a volume ratio of 3:1 is selected as the etching solution for sapphire. On this basis, three groups of etching temperatures are selected, namely 202 °C, 223 °C and 236 °C. The etching times of sapphire hemispheres at three temperatures are 32 h, 28 h and 24 h, respectively. Three groups of experimental parameters of sapphire hemispheres are shown in

Table 1. Three groups of experimental parameters of sapphire hemispheres.

Group	Volume Ratio	Etching Temperature (°C)	Etching Time (h)
1	H2SO4:H3PO4 = 3:1	202	32
2	H2SO4:H3PO4 = 3:1	223	28
3	H2SO4:H3PO4 = 3:1	236	24

.

The procedures of obtaining etch rates through etching experiments of sapphire hemispheres can be summarized as follows [28,29,30]:

(1)

Coordinate points on the three sapphire hemispherical surfaces are measured every two degrees. The measurement range is from −10° to 90° (latitude) and from 0° to 358° (longitude).

(2)

According to the experimental parameters shown in Table 1, three groups of etching experiments of sapphire hemispheres are carried out without stirring.

(3)

The coordinate points on the three sapphire hemispherical surfaces are measured again after etching. The step is the same as (1).

(4)

According to the ratio of the etching distances along the radial direction to the etching times of crystallographic planes, etch rates of sapphire under different temperature conditions (H₂SO₄:H₃PO₄ = 3:1, 202 °C, 223 °C and 236 °C) are obtained.

Etch rate distributions of sapphire at three different temperatures (H₂SO₄:H₃PO₄ = 3:1, 202 °C, 223 °C and 236 °C) are shown in Figure 2a–c. Etch rate distributions at three different temperatures all show the three-fold rotational symmetry. It can be seen from rate curves in the <−1 1 0 0> crystallographic zone (Figure 2d) and <1 1 −2 0> crystallographic zone (Figure 2e) that etch rates increase with the increase in the temperature, but rate inflection points and rate extreme points on the etch rate curves do not change. Etch rates of main crystallographic planes under three different temperature conditions are shown in

Table 2. Etch rates of main crystallographic planes under three different temperature conditions.

Temperature (°C)	(0 0 0 1) (μm/h)	(1 1 −2 18) (μm/h)	(1 −1 0 11) (μm/h)	(1 −1 0 4) (μm/h)
202	0.67	1.01	1.34	1.28
223	1.38	2.01	2.45	2.35
236	1.77	2.42	3.11	2.97

.

3. Results and Discussion

3.1. Accuracy Optimization of the Monte Carlo Simulation of Sapphire

The literature [18] simulated etch rates at the experimental condition (H₂SO₄:H₃PO₄ = 3:1, 236 °C) by selecting eleven key planes as the basis for checking the energy parameters in the removal probability equation (RPE) with the evolutionary algorithm. The simulation of the etch rate distribution obtained based on the optimized energy parameters is shown in Figure 3a. The etch rate distribution simulated in the literature [18] has three maximum rate zones, but the maximum rate zone cannot be split into two zones (zone A and zone B). It shows that the removal probabilities assigned to surface atoms on planes in the maximum rate zone is not accurate enough. Therefore, it is necessary to further classify the kinds of surface atoms in more detail. As shown in Figure 3c,d, based on the six indices classification method of surface atoms, the eight indices classification method of surface atoms $\left(n_{1s}^G,n_{1b}^G,n_{1s}^R,n_{1b}^R,n_{2s}^G,n_{2b}^G,n_{2s}^R,n_{2b}^R\right)$ is proposed. Based on the six indices classification method, the eight indices classification method further classifies the second atoms $2s$ and $2b$ into second-grey-bond surface atoms $2Gs$, second-grey-bond bulk atoms $2Gb$, second-red-bond surface atoms $2Rs$, and second-red-bond bulk atoms $2Rb$. Based on this, the original removal probability equation (RPE) is optimized to the expanded removal probability equation (E-RPE):

$$P\left(n_{1s}^G,n_{1b}^G,n_{1s}^R,n_{1b}^R,n_{2s}^G,n_{2b}^G,n_{2s}^R,n_{2b}^R\right)=P_0\times P_1^G\left(n_{1s}^G,n_{1b}^G\right)\times P_1^R\left(n_{1s}^R,n_{1b}^R\right)\times P_2^G\left(n_{2s}^G,n_{2b}^G\right)\times P_2^R\left(n_{2s}^R,n_{2b}^R\right)$$

(6)

Among them,

$$P_0=\left(1+e^{-\beta E_1}\right)\times \left(1+e^{-\beta E_2}\right)\times \left(1+e^{-\beta E_3}\right)\times \left(1+e^{-\beta E_4}\right)$$

(7)

$$P_1^G\left(n_{1s}^G,n_{1b}^G\right)=1/(1+e^{\beta \times F_1^G\left(n_{1s}^G,n_{1b}^G\right)}) \mathrm{and} F_1^G\left(n_{1s}^G,n_{1b}^G\right)={\epsilon }_1\times n_{1s}^G+{\epsilon }_2\times n_{1b}^G-E_1$$

(8)

$$P_1^R\left(n_{1s}^R,n_{1b}^R\right)=1/(1+e^{\beta \times F_1^R\left(n_{1s}^R,n_{1b}^R\right)}) \mathrm{and} F_1^R\left(n_{1s}^R,n_{1b}^R\right)={\epsilon }_3\times n_{1s}^R+{\epsilon }_4\times n_{1b}^R-E_2$$

(9)

$$P_2^G\left(n_{2s}^G,n_{2b}^G\right)=1/(1+e^{\beta \times F_2^G\left(n_{2s}^G,n_{2b}^G\right)}) \mathrm{and} F_2^G\left(n_{2s}^G,n_{2b}^G\right)={\epsilon }_5\times n_{2s}^G+{\epsilon }_6\times n_{2b}^G-E_3$$

(10)

$$P_2^R\left(n_{2s}^R,n_{2b}^R\right)=1/(1+e^{\beta \times F_2^R\left(n_{2s}^R,n_{2b}^R\right)}) \mathrm{and} F_2^R\left(n_{2s}^R,n_{2b}^R\right)={\epsilon }_7\times n_{2s}^R+{\epsilon }_8\times n_{2b}^R-E_4$$

(11)

Similar to the RPE, eight energy parameters ${\epsilon }_1$,${\epsilon }_2$, ${\epsilon }_3$, ${\epsilon }_4$, ${\epsilon }_5$, ${\epsilon }_6$, ${\epsilon }_7$ and ${\epsilon }_8$ are average energies; in addition, four energy parameters $E_1$, $E_2$, $E_3$ and $E_4$ are threshold energies.

Based on the proposed eight indices classification method and the extended removal probability equation (E-RPE), the removal probabilities assigned to surface atoms will be more accurate. Here, eleven key planes located in the <−1 1 0 0> and <1 1 −2 0> crystallographic zones are also selected as the basis for checking the energy parameters in the E-RPE. The software used in the paper is C++, and the computer employed consists of an Intel(R) Core(TM) i5-10400F CPU @ 2.90GHz processor, a RAM (DDR4) of 16 GB and an AMD Radeon RX 580 2048SP graphics card (the graphics card is not related to the Monte Carlo simulation). Therefore, a group of energy parameters at the experimental condition (H₂SO₄:H₃PO₄ = 3:1, 236 °C) are obtained: ${\epsilon }_1$ = 0.14, ${\epsilon }_2$ = 0.09, ${\epsilon }_3$= 0.57, ${\epsilon }_4$= 0.31, ${\epsilon }_5$ = 0.18, ${\epsilon }_6$ = 0.03, ${\epsilon }_7$ = 0.09, ${\epsilon }_8$ = 0.11, $E_1$ = 0.28, $E_2$ = 0.15, $E_3$ = 0.54 and $E_4$ = 66 (eV).

Smaller errors between simulated rates and experimental rates as shown in Figure 4a,b show that the energy parameters have been well optimized. Figure 4c shows the different periods of simulations of the etch rate distribution at the experimental condition (H₂SO₄:H₃PO₄ = 3:1, 236 °C). From the early period, the middle period to the later period of the simulation, the characteristics of the simulated etch rate distribution are more and more obvious, and the simulation in the later period is the best. However, in the last period, the simulation of the etch rate distribution begins to deteriorate. It indicates that the simulation time of wet etching of the hemisphere is too long, and the spherical surface has been over etched. Figure 4d shows the comparison between the etch rate distribution under the experimental condition, the previous simulation based on the six indices classification method in the literature [18] and the current simulation based on the eight indices classification method. The small errors between the experimental result and current simulation show that the proposed eight indices classification method makes the removal probabilities assigned to surface atoms more accurately, and can simulate the obvious characteristics of the etch rate distribution (the maximum rate zone is split into two zones).

3.2. Efficiency Optimization of the Monte Carlo Simulation of Sapphire

According to the extended removal probability equation (E-RPE), the higher the temperature variable $T$ is, the higher the removal probabilities of surface atoms are, and the higher the etch rates are. Theoretically, based on a group of energy parameters obtained under one temperature condition, etch rates under different temperature conditions can be simulated by changing only the temperature variable $T$. Therefore, the simulated rates in the temperature range of 202–304 °C will be obtained according to a set of energy parameters (${\epsilon }_1$ = 0.14, ${\epsilon }_2$ = 0.09, ${\epsilon }_3$ = 0.57, ${\epsilon }_4$ = 0.31, ${\epsilon }_5$ = 0.18, ${\epsilon }_6$ = 0.03, ${\epsilon }_7$ = 0.09, ${\epsilon }_8$ = 0.11, $E_1$ = 0.28, $E_2$ = 0.15, $E_3$ = 0.54 and $E_4$ = 1.66 (eV)) based on the eight-index classification method and compared with the corresponding experimental rates.

According to the experimental rates under the existing three groups of temperature conditions (H₂SO₄:H₃PO₄ = 3:1, 202 °C, 223 °C and 236 °C), the experimental rates within the range of 202–304 °C can be obtained through the Arrhenius formula [31,32]:

$$v=A{e}^{-E_a/K_{B}T}$$

(12)

Here, $v$ is the etch rate (μm/h), $A$ is the pre-exponential factor, $E_a$ is the activation energy of the crystallographic plane (eV), $K_B$ is the Boltzmann constant and $T$ is the absolute temperature (K).

Equation (12) can be transformed into Equation (13). According to Equation (13), the logarithm of $v$ is linear with the reciprocal of $T$.

$$\ln v=-\frac{E_a}{K_B}\times \frac{1}{T}+\ln A$$

(13)

For example, the linear equation of (0 0 0 1) plane can be fitted according to the experimental rates (0.67, 1.38 and 1.77 μm/h) of (0 0 0 1) plane at 202 °C, 223 °C and 236 °C, as shown in Figure 5a. Thus, the experimental rates of (0 0 0 1) plane in the range of 202–304 °C are obtained, as shown in Figure 5b. The experimental rates of other crystallographic planes in the range of 202–304 °C can be obtained in the same way.

In this paper, low rate plane (−1 1 0 4), medium rate plane (0 0 0 1) and high rate planes (1 −1 0 4) and (1 −1 0 11) located in the <1 1 −2 0> crystallographic zone are selected to analyze the errors between their simulated and experimental rates. The comparisons between simulated and experimental rates of these four planes are shown in Figure 6.

It can be seen from Figure 6 that the simulated rates increase with the increase in the temperature. It indicates that the effect of the temperature variable $T$ in the extended removal probability equation (E-RPE) is sensitive and positive. However, there are larger errors between simulated and experimental rates in the range of 202–304 °C. When 202 °C ≤ t < 236 °C, simulated rates are higher than experimental rates; when 236 °C < t ≤ 304 °C, simulated rates are lower than experimental rates. In addition, the error degree of etch rates within 202 °C ≤ t < 236 °C is less than that within 236 °C < t ≤ 304 °C.

Therefore, the extended removal probability equation (E-RPE) needs to be modified to reduce the errors between simulated and experimental rates. The modified removal probability equation (ME-RPE) must also satisfy that when the temperature is 236 °C, the ME-RPE is equal to the E-RPE.

According to the analysis for errors between simulated and experiment rates in the temperature range of 202–304 °C, the decreasing and increasing properties of the exponential function can be used to modify the extended removal probability equation (E-RPE). It can be seen from the expression of the exponential function $y=a^x$ that the value range of $a$ is $a$ > 0 and $a$ ≠ 1. When 0 < $a$ < 1, $y=a^x$ is a decreasing function; when $a$ > 1, $y=a^x$ is an increasing function. Therefore, the errors that simulated rates are higher than experimental rates in the temperature range of 202–236 °C can be modified by the decreasing property of $y=a^x$ (0 < $a$ < 1), and the errors that simulated rates are lower than experimental rates in the temperature range of 236–304 °C can be modified by the increasing property of $y=a^x$ ($a$ > 1). Figure 7 is the schematic diagram of the extended removal probability equation (E-RPE) modified by the properties of the exponential function, so that the absolute temperature $T_0=\left(236+273.15\right)K$ corresponds to the Celsius temperature $t_0$ = 236 °C, the ratio of $T$ ($t$ < 236 °C) to $T_0 (t_0$ = 236 °C) corresponds to $0<a<1$, and the ratio of $T$ ($t$ > 236 °C) to $T_0 (t_0$ = 236 °C) corresponds to $a>1$.

Therefore, temperature-modified factor $B_1={\left(\frac{T}{T_0}\right)}^m$ is proposed in this paper, and the modified temperature can be expressed as $T=T\times B_1$. Among them, $T_0=\left(236+273.15\right)K$, $T=\left(t+273.15\right)K$ and $t$ is the Celsius temperature (202–304 °C). Parameter $m$ is the modified index ($m$ > 0); the larger m is, the greater the modified range of rates are. The effect of the temperature-modified factor $B_1$ is shown in Figure 8a. When 202 °C ≤ t < 236 °C, $B_1$ < 1, the value of the modified $T$ decreases, resulting in a decrease in the etch rates $v$. When 236 °C < t ≤ 304 °C, $B_1$ > 1, the value of the modified $T$ increases, resulting in an increase in the etch rates $v$. It can be seen from Figure 8b that with the increasement of the modified index $m$, etch rates $v$ decrease continuously in the range of 202 °C ≤ t < 236 °C, while etch rates $v$ increase continuously in the range of 236 °C < t ≤ 304 °C. Since the error degrees between simulated and experimental rates are different in the two ranges of 202 °C ≤ t < 236 °C and 236 °C < t ≤ 304 °C, when $m$ reaches a certain value, the rate errors in the range of 236 °C < t ≤ 304 °C are small, and the rate errors in the range of 202 °C ≤ t < 236 °C are still large, as shown in the orange curve and black curve in Figure 8b.

Therefore, in order to reduce the errors within the range of 202 °C ≤ t < 236 °C, temperature-modified factor $B_1$ is further optimized to $B_2=B_1{}^{{\left(\frac{T}{T_0}\right)}^n}$, that is, $B_2={\left(\frac{T}{T_0}\right)}^{m\times {\left(\frac{T}{T_0}\right)}^n}$. Among them, $T_0=\left(236+273.15\right)K$, $T=\left(t+273.15\right)K$, m and n are modified indices (m > 0 and n > 0). The modified temperature can be expressed as $T=T\times B_2$. The effect of the temperature-modified factor $B_2$ is shown in Figure 8c. The modified result is shown in Figure 8d. When m and n reach certain values, the errors between simulated and experimental rates in the two ranges of 202 °C ≤ t < 236 °C and 236 °C < t ≤ 304 °C are both small. In addition, when t = 236 °C, the ME-RPF is equal to the E-RPF.

Therefore, low rate plane (−1 1 0 4), medium rate plane (0 0 0 1) and high rate plane (1 −1 0 4) and (1 −1 0 11) are selected as the basis for checking the modified indices m and n. After several tests on the values of m and n, when m = 2.8 and n = 4.2, the errors between simulated and experimental rates in the temperature range of 202 °C ≤ t ≤ 304 °C are small, as shown in Figure 9. Therefore, the modified removal probability equation (ME-RPE) is obtained. Among them, $B={\left(\frac{T}{T_0}\right)}^{2.8\times {\left(\frac{T}{T_0}\right)}^{4.2}}$, $T_0=\left(236+273.15\right)K$, $T=\left(t+273.15\right)K$ and 202 °C ≤ t ≤ 304 °C.

In order to verify the effectiveness of the modified removal probability equation (ME-RPE) for temperature change, etch rates under the other two temperature conditions (H₂SO₄:H₃PO₄ = 3:1, 202 °C and 223 °C) are simulated by changing the temperature variable $T$. Small errors between simulated rates and experimental rates under the two temperature conditions (H₂SO₄:H₃PO₄ = 3:1, 202 °C and 223 °C) shown in Figure 10 and Figure 11 verify the effectiveness of the modified removal probability equation (ME-RPE) on temperature change.

4. Conclusions, Limitations, and Future Research

In this paper, the Monte Carlo simulation for sapphire in wet etching is effectively optimized, which improves the accuracy and efficiency of the simulation. On the basis of the six indices classification method, this paper classifies the kinds of surface atoms in more detail and proposes an eight indices classification method $\left(n_{1s}^G,n_{1b}^G,n_{1s}^R,n_{1b}^R,n_{2s}^G,n_{2b}^G,n_{2s}^R,n_{2b}^R\right)$. Accordingly, the E-RPE is established. Based on this, the accuracy of the simulated etch rate distribution has been greatly improved (the maximum rate zone can be split into two zones). Furthermore, based on the role of the temperature variable in the E-RPE, the errors between simulated and experimental rates under different temperature conditions are analyzed. Based on the decreasing and increasing properties of the exponential function, the E-RPE is modified, and the ME-RPE is obtained. Based on this, the small errors between simulated rates and experimental rates under different temperature conditions verify the effectiveness of the ME-RPE for temperature change. The ME-RPE significantly improves the efficiency of the Monte Carlo simulation. The optimization of the Monte Carlo simulation for sapphire in wet etching in this paper can also provide a reference for the simulation of wet etching of other crystal materials, such as GaN, SiC and LiNbO₃. However, the optimization for the Monte Carlo simulation in this paper is still based on the simplified atomic structure of sapphire. In the future, the Monte Carlo simulation will be carried out based on the original atomic structure of sapphire to further improve the accuracy of the simulation.