Modeling and extraction of residual stresses in some MEMS multilayer systems

The intrinsic stresses introduced during the MEMS multilayer deposition process and the thermal stresses introduced during operation can significantly impact the performance, reliability, and lifetime of the product. The mechanical properties of the MEMS multilayer system are analyzed. An accurate closed-form solution is obtained for the multilayer system, including both intrinsic stresses and thermal stresses, and validated by numerical FEM simulation. Thin films are deposited on silicon substrates, and wafer bending measurements are carried out. Based on this formalism, the residual stresses in various multilayer systems are extracted and analyzed using parameter estimation theory.

1 Introduction

Residual stress is the stress that exists in an equilibrium solid without external loads. Physically, it is the sum of interatomic forces and the momentum flux intercepted by a (microscopically big) area element, divided by its area (Subramaniyan and Sun, 2008). There are intrinsic and extrinsic residual stresses. Intrinsic residual stresses include stresses induced by lattice mismatch, plastic deformation, grain growth, and other processes. Extrinsic residual stresses are induced by changes in external environmental conditions, such as temperature or humidity, and are reversed when these conditions return to their original values. The generation and evolution of residual stresses have very different origins at many different scales. They can have microscopic origins, such as point defects in solids, or mesoscopic origins, such as dislocations, grain boundaries, grains, and inclusions. They can also have macroscopic origins, such as inhomogeneity.

Residual stresses in MEMS structures can affect the reproducibility, product lifetime, and reliability of the devices because they can modify the structural, electrical, and optical properties of many MEMS materials. Residual stresses can also affect the yield of the MEMS production process because they can distort the wafers and chips. The measurement of residual stresses is therefore very important (Dutta and Pandey, 2021). Measurement methods include the X-ray method, Raman spectroscopy method, microscopic FIB method, and wafer bow method (Abadias et al., 2018; Huff, 2022). The recent cross-sectional nano-diffraction method uses an X-ray beam with a diameter of 50 nm, which can resolve the stress gradients and measure the stress in thin films (Zeilinger et al., 2016). The X-ray method measures residual strain. To calculate the residual stress, Young's modulus is also used, like in the wafer bow method (Freund and Suresh, 2010). One limitation of this method is that it is not suitable for amorphous films. Raman spectroscopy can extract the residual stress from the Raman shift of the specific lattice vibration spectral line. The deformation potential involved in the extraction must be calculated or calibrated (Ma et al., 2021). This method is not useful for amorphous films (Huff, 2022). MEMS technology involves many thin films. The most frequently used films, such as thermal oxide, CVD-deposited oxide, and nitride films, are amorphous. Therefore, the wafer bow measurement method is widely used for MEMS technology due to its non-destructive nature and ease of integration into the manufacturing process. Wafer bow measurements can be taken using the 3D laser line scanning method (Schlarp et al., 2018) or a simple two-laser beam set-up (Janssen, 2007).

The kinetic evolution of residual stresses in multicrystalline films has been systematically investigated (Engwall et al., 2016, and Chason et al., 2012). The results show a general tendency toward tensile stress due to the grain growth in the initial stage and compressive stress due to the net flux of atoms into the top of the grain boundary during film growth in the latter stage (Chason et al., 2018). The formation and evolution of residual stresses in amorphous films are more complex and dependent on the growth process and history (Murry, 1994 and Chen et al., 2003). Therefore, a convenient residual stress-monitoring and evaluation method would be useful.

While the mechanical properties of multilayer systems can be numerically simulated through finite-element analysis, the modeling of realistic structures, such as large wafers composed of multiple thin layers with high aspect ratios, requires substantial computational resources. Its application to process optimization is even more demanding. Hence, an analytical solution is desirable. Stoney's formula (Stoney, 1909) of the bending of a double metal layer due to residual stresses from manufacturing was the first development in this direction and has led to the application of this geometrical method in many fields. For extrinsic residual stresses, i.e., thermal stresses, this theory has been extended to multilayer structures (Obata, 2014; Shi et al., 2022). For intrinsic residual stresses, an equivalent reference temperature (ERT) technique has been developed and combined with laminate theory to model multilayer structures (Nejhad et al., 2003). In this paper we extend these formalisms to more general cases, including both intrinsic residual stresses and thermal stresses in multilayers simultaneously. The equilibrium conditions and continuity of strains in multilayer systems are used to derive the solution. Instead of the general anisotropic constitutive equation used in Nejhad et al. (2003), a linear isotropic elastic property is assumed. This leads to an accurate closed-form solution. This approach can be conveniently used to extract residual stresses in multilayer systems and optimize the MEMS production process flow.

Assuming that the residual stresses of underlying layers will not change during the deposition or etching of new layers, a methodology is proposed to extract residual stresses from multiple wafer bow measurements. Various multilayer structures are realized through successive deposition or successive etching. The wafer bows are measured using 3D laser line scanning, and the residual stresses in the different layers are extracted using the developed formalism. Unlike the X-ray and Raman spectroscopy methods, which can extract stresses without layer geometry information, this model requires layer thickness and wafer size to be known. Because possible errors exist in the involved layer thicknesses, wafer bow measurements, and layer properties, we apply a parameter estimation theory to our formalism to study error propagation and parameter precision. Finally, the formalism, the residual stress extraction methodology, its potential applications, and the residual stresses extracted for the various layers are discussed.

2 The mechanical properties of a multilayer system

2.1 General multilayer system

We consider a beam with a multilayer structure of n layers with Young's modulus E_i, Poisson's ratio ν_i, thermal expansion coefficient (TEC) α_i, intrinsic residual stress σ_i, and homogeneous thickness h_i for layer i, where $i=\mathrm{1},\mathrm{2},\mathrm{\dots },n$. Due to the residual stresses presented in the layers, the multilayer structure may experience a bending with radius R to balance the residual stresses. The origin of the y axis is at the bottom of the first layer (Fig. 1). The system is supposed to be at a temperature ΔT above the reference temperature T_r. Then, in the isotropic linear elastic approximation, the elastic stress σ_xi at layer i can be expressed as

$$\begin{array}{}\text{(1)}& {\mathit{\sigma }}_{xi}\left(y\right)=E_i({\mathit{ϵ}}_{\mathrm{0}}-{\mathit{\alpha }}_i\mathrm{\Delta }T)+E_i\mathit{\kappa }y+{\mathit{\sigma }}_i(i=\mathrm{1},\mathrm{2},\mathrm{\dots },n),\end{array}$$

where $\mathit{\kappa }=\frac{\mathrm{1}}{R}$ is the curvature of the bent structure, and ϵ₀ is the net strain. For a circular wafer with cylindrical symmetry, the biaxial modulus $E_i^{\prime }=\frac{E_i}{\mathrm{1}-{\mathit{\nu }}_i}$ is used for E_i. According to the virtual work principle, the system satisfies the equilibrium conditions of the total force and bending moment (Obata, 2014). Solving the equilibrium equations, we obtain the following two parameters:

$$\begin{array}{}\text{(2)}& {\displaystyle }\mathit{\kappa }={\displaystyle \frac{P_{\mathrm{2}}{I}_{\mathrm{1}}-P_{\mathrm{1}}{I}_{\mathrm{2}}}{I_{\mathrm{1}}{I}_{\mathrm{3}}-I_{\mathrm{2}}^{\mathrm{2}}}},\text{(3)}& {\displaystyle }{\mathit{ϵ}}_{\mathrm{0}}={\displaystyle \frac{P_{\mathrm{1}}{I}_{\mathrm{3}}-P_{\mathrm{2}}{I}_{\mathrm{2}}}{I_{\mathrm{1}}{I}_{\mathrm{3}}-I_{\mathrm{2}}^{\mathrm{2}}}},\end{array}$$

where

$$\begin{array}{}\text{(4)}& {\displaystyle }{I}_{\mathrm{1}}=\sum _{i=\mathrm{1}}^n{E}_i{h}_i,\text{(5)}& {\displaystyle }{I}_{\mathrm{2}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\sum _{i=\mathrm{1}}^n{E}_i(y_i^{\mathrm{2}}-y_{i-\mathrm{1}}^{\mathrm{2}}),\text{(6)}& {\displaystyle }{I}_{\mathrm{3}}={\displaystyle \frac{\mathrm{1}}{\mathrm{3}}}\sum _{i=\mathrm{1}}^n{E}_i(y_i^{\mathrm{3}}-y_{i-\mathrm{1}}^{\mathrm{3}}),\text{(7)}& {\displaystyle }{P}_{\mathrm{1}}=\sum _{i=\mathrm{1}}^n(E_i{\mathit{\alpha }}_i\mathrm{\Delta }T-{\mathit{\sigma }}_i)h_i,\text{(8)}& {\displaystyle }{P}_{\mathrm{2}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\sum _{i=\mathrm{1}}^n(E_i{\mathit{\alpha }}_i\mathrm{\Delta }T-{\mathit{\sigma }}_i)(y_i^{\mathrm{2}}-y_{i-\mathrm{1}}^{\mathrm{2}}).\end{array}$$

2.2 Simplification in special cases

For a two-layer system, n=2, one obtains the following at the reference temperature, ΔT=0:

$$\begin{array}{}\text{(9)}& \mathit{\kappa }={\displaystyle \frac{-\mathrm{6}(E_{\mathrm{1}}{\mathit{\sigma }}_{\mathrm{2}}-E_{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{1}})(h_{\mathrm{1}}{h}_{\mathrm{2}}^{\mathrm{2}}+h_{\mathrm{1}}^{\mathrm{2}}{h}_{\mathrm{2}})}{(E_{\mathrm{1}}{h}_{\mathrm{1}}^{\mathrm{2}}+E_{\mathrm{2}}{h}_{\mathrm{2}}^{\mathrm{2}}{)}^{\mathrm{2}}+\mathrm{4}{E}_{\mathrm{1}}{E}_{\mathrm{2}}{h}_{\mathrm{1}}{h}_{\mathrm{2}}(h_{\mathrm{1}}^{\mathrm{2}}+h_{\mathrm{1}}{h}_{\mathrm{2}}+h_{\mathrm{2}}^{\mathrm{2}})}}.\end{array}$$

Furthermore, for a single thin film deposited on a substrate at the reference temperature, if the top layer (i=2) is much thinner than the substrate (i=1), i.e., h₂≪h₁, and the substrate is free of residual stress, σ₁=0, the above formula is reduced to the Stoney formula.

$$\begin{array}{}\text{(10)}& \mathit{\kappa }={\displaystyle \frac{-\mathrm{6}{\mathit{\sigma }}_{\mathrm{2}}{h}_{\mathrm{2}}}{E_{\mathrm{1}}{h}_{\mathrm{1}}^{\mathrm{2}}}}\end{array}$$

A situation of interest is a thin multilayer on a thick substrate. If n thin layers are deposited on a residual-stress-free substrate with Young's modulus E_s and thickness h_s, such that $h_{\mathrm{s}}\gg {\sum }_{i=\mathrm{1}}^n{h}_i$, then

$$\begin{array}{}\text{(11)}& \mathit{\kappa }=-{\displaystyle \frac{\mathrm{12}{M}_{\mathrm{t}}+\mathrm{6}{F}_{\mathrm{t}}{h}_{\mathrm{s}}}{E_s{h}_{\mathrm{s}}^{\mathrm{3}}}}.\end{array}$$

The average elastic stress in layer i is approximately

$$\begin{array}{}\text{(12)}& \begin{array}{rl}{\mathit{\sigma }}_{xi}& ={\displaystyle \frac{\mathrm{6}{M}_{\mathrm{t}}-\mathrm{4}{F}_{\mathrm{t}}{h}_{\mathrm{s}}}{E_{\mathrm{s}}{h}_{\mathrm{s}}^{\mathrm{2}}}}{E}_i+{\mathit{\sigma }}_i-E_i({\mathit{\alpha }}_i-{\mathit{\alpha }}_{\mathrm{s}})\mathrm{\Delta }T\\ & (i=\mathrm{1},\mathrm{2},\mathrm{\dots },n),\end{array}\end{array}$$

where F_t and M_t are the total linear force density and the total linear moment density of the thin-layer system:

$$\begin{array}{}\text{(13)}& {\displaystyle }{F}_{\mathrm{t}}=\sum _{i=\mathrm{1}}^n({\mathit{\sigma }}_i-E_i({\mathit{\alpha }}_i-{\mathit{\alpha }}_{\mathrm{s}}\left)\mathrm{\Delta }T\right)h_i,\text{(14)}& {\displaystyle }{M}_{\mathrm{t}}=\sum _{i=\mathrm{1}}^n({\mathit{\sigma }}_i-E_i({\mathit{\alpha }}_i-{\mathit{\alpha }}_{\mathrm{s}}\left)\mathrm{\Delta }T\right)h_i{y}_{i\mathrm{c}},\end{array}$$

with y_ic being the center position of thin layer i. This can be used in the sensor design to optimize the chip deformation and elastic stress in the function layer. It can also be used to design optical thin-film interference filters. By adjusting the thickness and residual stress of the individual films, a flat filter can be ensured. This is because the alternative tensile and compressive residual stresses can be compensated for to achieve a small curvature.

Figure 1 The bent multilayer structure of n layers with corresponding elastic parameters E_i,ν_i, residual stresses σ_i, the thermal expansion coefficient α_i, and thickness h_i. The bending radius is R.

2.3 Validation of the theory

To verify the theory, we simulate a double-layer cantilever composed of a 10 µm silicon layer and a 1 µm nitride layer with COMSOL. The cantilever is 200 µm in length. Assuming a residual stress of 1×10⁸ Pa in the nitride layer, we obtain the elastic distribution at position x=90 µm, as shown in Fig. 2. According to the Saint-Venant principle, the effects of the boundary constraints can be neglected in a region sufficiently far from the boundary (von Mises, 1945). This elastic stress distribution can be compared with the results calculated using our theory (Fig. 2), demonstrating good agreement between the calculated and simulated results.

Figure 2 The elastic stress distribution in the double layer simulated with this theory (red line) and with COMSOL (dotted blue line).

3 Application to MEMS technology

The equibiaxial residual stress problem in a deposited circular wafer is axisymmetric. It can be shown that the 2D-axisymmetric problem is equivalent to the beam problem from the previous section, with the symmetric center at at the fixed end of the beam (Timoshenko and Goodier, 2010). The same formalism (Sect. 2.1) for the beam can be applied, with the axial residual stress interpreted as the value of the equibiaxial stress. In MEMS technology, many implantations, depositions, and annealing processes are involved. These processes induce residual stresses, resulting in wafer bending and chip deformation. This formalism can be used to ensure that the wafer bow remains within the process tolerance throughout the entire process flow by optimizing the film thickness and process parameters. To do so, the process-dependent residual stresses must be known. This can be achieved through the inverse problem of the formalism by starting with the measured wafer curvature and deducing the residual stresses.

3.1 Extraction of residual stresses of the multilayer system

If the residual stresses are kept constant after deposition, they can be extracted with wafer bow measurements. We consider a MEMS structure involving n layers with l unknown residual stresses that need to be extracted. Then, l≤n, because the residual stress in some layers may be known or assumed to be stress free. One designs an experiment with m wafers with different layer structures from the same process flow. These different film-stacking structures can be achieved by inserting or taking out wafers from the same process flow. This is the successive deposition approach, where the deflection of each deposited wafer is measured. The same can be achieved after the entire process by etching the deposited wafer layer by layer. This is the successive-etching approach, where the deflection of each etched wafer is measured. This yields a data set of m measurements. The flat wafer bends due to the residual stress of the deposited thin films. The bending profile can be approximated by a spherical cap with curvature κ. The wafer center shifts a deflection g (wafer bow) relative to the wafer rim. The deflections of the wafers are measured using 3D laser line scanning based on the principle of optical triangulation (Schlarp et al., 2018). For a wafer of diameter d, one can obtain the deflection from the bending curvature (κ) as

$$\begin{array}{}\text{(15)}& g=-{\displaystyle \frac{d^{\mathrm{2}}}{\mathrm{8}}}\mathit{\kappa }\end{array}$$

if the bending radius is much bigger than the wafer size. One can thus get m equations for the l unknown residual stresses. The residual stress of a given layer is assumed to be the same when it appears in different combinations of layer stacks because it is formed in the same process and under the same conditions. With an optimization method, one can extract the l residual stresses if m is bigger than or equal to l. A Python program was written for data analysis. The fsolve and minimize functions of the scipy module optimize are used for stress extraction in the cases of m=l and m>l, respectively.

This formalism can be generalized to analyze inhomogeneous residual stresses in a multilayer system. In this case, small cantilever test structures should be implemented on the wafer. In Huber et al. (2026), cantilevers were fabricated, and the local residual stress in a single overlay was analyzed using the Stoney formula. Our formalism can be applied to a multilayered cantilever by replacing the wafer radius with the cantilever length. Instead of using the biaxial modulus for the wafer bow, the original Young's modulus should be used for the cantilever deflection.

3.2 Parameter error estimation

The precision of the fitted residual stress depends not only on the fitting residuals themselves, but also on the measurement errors of deflection and thickness, as well as the parameter errors of Young's modulus and the thermal expansion coefficient of the films. Here, we neglect the effect of possible uncertainties in Poisson's ratios. These errors can propagate to the uncertainty in the extracted residual stresses. For measurement j, the theoretical deflections are

$$\begin{array}{}\text{(16)}& \begin{array}{rl}{g}_j& =f({\mathit{\sigma }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\sigma }}_l,h_{\mathrm{1}},\mathrm{..},h_n,E_{\mathrm{1}},\mathrm{\dots },E_n,{\mathit{\alpha }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\alpha }}_n)\\ & (j=\mathrm{1},\mathrm{2},\mathrm{\dots },m).\end{array}\end{array}$$

The total standard deviation of deflection δ_t(g_j) is therefore

$$\begin{array}{}\text{(17)}& \begin{array}{rl}{\mathit{\delta }}_{\mathrm{t}}^{\mathrm{2}}\left(g_j\right)& ={\mathrm{res}}^{\mathrm{2}}\left(g_j\right)+{\mathit{\delta }}^{\mathrm{2}}\left(g_j\right)+\sum _{i=\mathrm{1}}^n{\left({\displaystyle \frac{\partial g_j}{\partial h_i}}\mathit{\delta }\left(h_i\right)\right)}^{\mathrm{2}}\\ & +\sum _{i=\mathrm{1}}^n{\left({\displaystyle \frac{\partial g_j}{\partial E_i}}\mathit{\delta }\left(E_i\right)\right)}^{\mathrm{2}}+\sum _{i=\mathrm{1}}^n{\left({\displaystyle \frac{\partial g_j}{\partial {\mathit{\alpha }}_i}}\mathit{\delta }\left({\mathit{\alpha }}_i\right)\right)}^{\mathrm{2}}\\ & (j=\mathrm{1},\mathrm{2},\mathrm{\dots },m),\end{array}\end{array}$$

where res(g_j) is the fitting residual, which directly determines the fitting error; δ(g_j) is the measurement error of the deflection of the jth measurement; δ(h_i) is the thickness error of the ith layer; and δ(E_i) and δ(α_i) are the parameter errors of Young's modulus and the thermal expansion coefficient of the ith layer, respectively. They determine the propagated error of the fitted residual stresses. One can obtain the uncertainties in the l fitted unknown residual stresses from sensitivity analysis (Cowan, 1998). The m×l Jacobian matrix J, also known as the sensitivity matrix, can be calculated numerically with the following elements:

$$\begin{array}{}\text{(18)}& {\mathbf{J}}_{jk}={\displaystyle \frac{\partial g_j}{\partial {\mathit{\sigma }}_k}}(j=\mathrm{1},\mathrm{2},\mathrm{\dots },m)(k=\mathrm{1},\mathrm{2},\mathrm{\dots },l).\end{array}$$

The l×l Fisher information matrix of the fitted parameters is approximately (Ober et al., 2003)

$$\begin{array}{}\text{(19)}& \mathbf{F}={\mathbf{J}}^T\mathbf{GJ},\end{array}$$

where

$$\begin{array}{}\text{(20)}& \mathbf{G}=\mathrm{diag}\left({\displaystyle \frac{\mathrm{1}}{{\mathit{\delta }}_{\mathrm{t}}^{\mathrm{2}}\left(g_{\mathrm{1}}\right)}},{\displaystyle \frac{\mathrm{1}}{{\mathit{\delta }}_{\mathrm{t}}^{\mathrm{2}}\left(g_{\mathrm{2}}\right)}},\mathrm{\dots },{\displaystyle \frac{\mathrm{1}}{{\mathit{\delta }}_{\mathrm{t}}^{\mathrm{2}}\left(g_m\right)}}\right).\end{array}$$

The inversion of Fisher's information matrix is the covariance matrix C of fitting parameters $\mathbf{C}={\mathbf{F}}^{-\mathrm{1}}$. The standard deviations, the correlation matrix, and the confidence intervals of the fitted residual stresses can be calculated from the covariance matrix (Tarantola, 2005). The relative error of the kth fitted residual stress is

$$\begin{array}{}\text{(21)}& {\displaystyle \frac{\sqrt{{\mathbf{C}}_{kk}}}{{\mathit{\sigma }}_k}}(k=\mathrm{1},\mathrm{2},\mathrm{\dots },l).\end{array}$$

With nonzero errors, δ(g_j), δ(h_i), δ(E_i), and δ(α_i), one can see the error propagation.

3.3 Residual stress extraction with successive deposition method

We designed an experiment (experiment I) to study the residual stresses of four thin films – oxidation layer (FOX), MoSi, CVD oxide, and stress-reduced LPCVD nitride – to investigate the generation and evolution of residual stresses in different layers involved in a MEMS production process. Since the MoSi layer will be annealed, we designed a five-step experiment using five starting wafers of (100)-Si with a thickness of 390 µm. After the thermal oxidation process, wafer 1 was removed from the flow. A MoSi layer was then deposited on the remaining wafers, after which wafer 2 was removed from the flow. The remaining three wafers were annealed, and wafer 3 was removed from the flow. LPCVD oxide deposition was performed on the remaining two wafers, after which wafer 4 was removed from the flow. Finally, a stress-reduced LPCVD nitride deposition was performed on wafer 5. The process steps and the wafers used are shown in Table 1. The wafer bows of the five wafers with different layer stacking were measured using 3D laser line scanning. As an example, Fig. 3 shows the wafer scans for wafer 2 (with the MoSi layer as deposited) and wafer 3 (after high-temperature annealing of the deposited MoSi layer). The wafer bending profiles aligned well with the arc of a circle, except for the rim of the wafer. This agrees with the Saint-Venant principle. This demonstrates that the curvature can be easily calculated from wafer deflection using Eq. (15). The results indicate that the stress in the MoSi layer changed drastically after annealing. The mechanical properties of the layers involved in this paper are shown in Table 2. For stress extraction, the wafer bows shown in Table 1, which have been corrected with the initial wafer bow, are used. Using the error analysis theory (Sect. 3.2), we calculated the relative fitting errors with ideal experimental parameters and the total relative errors, which took into account the propagated error from the 1 % relative error in wafer deflection, and the same relative error in thickness, Young's modulus, and the thermal expansion coefficient of all involved layers. The layer thicknesses, the extracted residual stresses, the relative fitting errors, and the total relative errors of fitted residual stresses are shown in Table 3. The results show that the residual stress of oxide depends on the preparation method and preceding layers. The exact kinetics are still unknown and must be studied using more systematic technological parameters. Furthermore, the correlation matrix of the fitted parameters is shown in Table 4, which also includes the propagated errors. A large anticorrelation coefficient is observed between the residual stresses of annealed MoSi and oxide. Although the total relative error of the fitted parameters is small, this correlation distorts the confidence ellipse from an ideal circle and will increase the size of the profile confidence interval (Jeyaratnam, 1992).

Figure 3 The wafer profile measured using laser scanning (dotted red line) and the fitting to an arc (blue line) for wafer 2 (left) and wafer 3 (right). Wafer 2, with an as-deposited MoSi layer, exhibits tensile stress. Wafer 3, with an annealed MoSi layer at 1000 °C for 30 min, shows a compressive stress.

Table 1 The process flow and wafer bows in µm of experiment I.

Table 2 Mechanical properties.

Table 3 Layer thicknesses, extracted residual stresses, relative fitting errors, and total relative errors of the extracted residual stresses of experiment I.

Table 4 Symmetrical correlation matrix of the extracted residual stresses of experiment I.

3.4 Residual stress extraction with successive-etching method

For experiment II, we start from a 100 mm diameter n-type (100) silicon wafer, which was oxidized with 60 nm silicon dioxide and deposited with 1 µm silicon nitride with LPCVD. The silicon nitride layer was deposited with a high DCS(dichlorosilane):NH₃ flow ratio to effectively reduce the residual stress (Zheng et al., 2013). To extract the residual stress in the thin layers, the layer thickness and wafer bow are measured at different stages of multilayer etching. The first bow measurement, D1, was made for the initial processed wafer without any etching. The second bow measurement, D2, was made after the backside nitride layer was removed. The third bow measurement, D3, was made after the backside oxide layer was removed. The fourth bow measurement, D4, was made after the topside nitride layer was removed. Finally, the fifth bow measurement, D5, was made after the topside oxide layer was removed. The experimental results, the extracted residual stresses, the relative fitting errors, and the total relative errors in the fitted residual stresses (as in Sect. 3.3) are shown in Table 5. The results show that the relative fitting error is very small. The fitted residual stress of the nitride film, 1.3×10⁸ Pa, agrees well with the value of 1.35×10⁸ Pa presented in Zheng et al. (2013). This also confirms the effectiveness of the high DCS:NH₃ ratio in reducing the residual stress of LPCVD nitride films. The residual stresses of normal LPCVD nitride films is 1.3×10⁹ Pa (Sect. 3.5), while the residual stress of nitride prepared in experiment I (Sect. 3.3) is only reduced to 4.9×10⁸ Pa.

Table 5 Layer thicknesses, wafer bows, extracted residual stresses, relative fitting errors, and total relative errors of the extracted residual stresses of experiment II.

3.5  Extraction of thermal stresses and intrinsic residual stresses

For experiment III, two 100 mm diameter n-type (100) silicon wafers, both 498 µm thick, are deposited with 300 nm of silicon dioxide and 100.7 nm of silicon nitride, respectively. The deposited film on the backside is etched away. The Tencor FLX2900 system is used to perform wafer curvature measurements with temperatures scanned from 50 to 350 °C. For stress extraction, it is assumed that the intrinsic residual stresses and thermal expansion coefficient (TEC) of the films remain constant during wafer heating measurements. Figure 4 shows the measured and simulated wafer curvature as a function of temperature. The fitted residual stress of silicon dioxide film is $-\mathrm{1.3}\times {\mathrm{10}}^{\mathrm{8}}$ Pa, and that of silicon nitride film is 1.3×10⁹ Pa. The relative fitting errors are $\mathrm{8.4}\times {\mathrm{10}}^{-\mathrm{5}}$ and $\mathrm{4.3}\times {\mathrm{10}}^{-\mathrm{4}}$ for silicon dioxide and silicon nitride films, respectively. The corresponding total relative errors are $\mathrm{1.3}\times {\mathrm{10}}^{-\mathrm{3}}$ and $\mathrm{2.3}\times {\mathrm{10}}^{-\mathrm{3}}$, respectively, if the propagated errors are taken into account like in Sect. 3.3. The thermal stresses are induced by the TEC differences between the film and substrate. Because the experiments contain errors in both the curvature measurement and the temperature reading, the extracted residual stresses also contain propagated uncertainty. Using these parameters, one can calculate the average elastic stress in the thin layer, as illustrated in Fig. 5.

It has been demonstrated that this formalism can address both intrinsic and extrinsic residual stresses in a layered structure. This is a useful extension of established approaches. For example, the widely used Stoney formula (Chason et al., 2018) can be applied to a single layer deposited on substrate. Thermal stresses can be analyzed in a layered structure with any number of layers (Shi et al., 2022), but intrinsic stresses induced by manufacturing are not incorporated. The equivalent reference temperature method (Nejhad et al., 2003) for dealing with intrinsic residual stresses can also model multilayer systems; however, no closed-form solution can be obtained. Therefore, it is difficult for them to be used in the reverse problem of general residual stress extraction. Furthermore, a systematic methodology for parameter extraction is presented using either successive deposition or etching methods, and a parameter estimation method is proposed. Additionally, this closed-form formalism can be incorporated into established device design flows and online stress-monitoring software as a software gadget to improve device performance and process yield.

Figure 4 Experiment IV – the measured (dot) and simulated (line) curvatures at different temperatures for wafer S5 (SiO₂ film, blue) and wafer S6 (Si₃N₄ film, red).

Figure 5 Experiment IV – the average elastic stresses in thin film at different temperatures for wafer S5 (SiO₂ film, blue) and wafer S6 (Si₃N₄ film, red).

4 Discussion, summary, and outlooks

For properly designed experiments with either successive etching or successive deposition and wafer bow measurements, the residual stresses in different layers can be extracted. Except for a few special layers, the total relative errors (including all possible propagated errors) of the fitted residual stresses are reasonably small. It can be seen that, except for nitride and the as-deposited MoSi layer, all layers have compressive residual stress. The polycrystalline MoSi film exhibits compressive stress after high-temperature annealing, which agrees with the model (Chason et al., 2018). We investigated how different processes affect residual stress, especially for LPCVD nitride, and confirmed the effectiveness of the procedure suggested in Zheng et al. (2013) for reducing tensile stress.

Theoretical expressions have been derived for wafer bending of multilayer systems caused by isotropic, homogeneous residual stresses stemming from both process-induced and thermal stresses at arbitrary temperature. This is an extension of the thermal stress theory for multilayer systems (Shi et al., 2022) and the equivalent reference temperature technique for the intrinsic stresses in layered structures (Nejhad et al., 2003). A Python-based extraction procedure has been developed. Using this procedure, the residual stress in each layer can be determined from a set of wafer bow measurements. Error propagation and the total relative error are estimated for the extracted residual stresses by applying the parameter estimation theory to our cases. Although many error factors are involved in wafer bending, it has been shown that error propagation can be kept small. To improve parameter estimation, the layer structure should be chosen to minimize parameter correlation and error propagation for residual stress extraction. Once the residual stresses associated with the different process steps have been determined, the process flow and layer structure can be designed using the theory presented in Sect. 2.1. This allows the process to simultaneously achieve optimal chip performance and ensures a trouble-free flow. This increases the chance of success in a single run for the development of new devices.

In MEMS technology, layers are often patterned in periodic structures or a porous layer is applied. These form a type of metastructure. The equivalent elastic properties (Jhung et al., 2009) can be used in this formalism to carry out the wafer bow analysis. Some MEMS devices, such as microemitters, operate at high temperatures. The additional thermal stresses can accelerate electrochemical migration, influencing the stability and lifetime of the devices. Our measurements and calculations of thermal stress show that an embedded prestress can alleviate the problem. Thus, the design of microsensors and microactuators can be optimized to minimize chip bowing and maintain low elastic stresses in the functional layers under all operating conditions. This approach ensures a higher production yield and the long-term reliability of sensors and actuators.

Modern MEMS technology involves many layers in the process. To guarantee optical, electrical, and mechanical performance, process conditions are modified and optimized. This often results in the generation of residual stresses. These stresses can be modified or relaxed in subsequent processes. All of these factors make grasping residual stresses challenging. The induced wafer bow can affect process precision because lithographic steps must be performed on a deformed wafer. After processing, a deformed wafer can make the automated wafer-level test difficult. Therefore, managing residual stress generation and evolution throughout the entire process is important. On the other hand, wafer bending is a direct indicator of residual stress. The wafer bow test is a reliable, non-destructive method of measuring residual stress. The presented closed-form formalism can facilitate its application to multilayer systems. In the future, the kinetics of residual stress generation and evolution in amorphous films should be systematically investigated. In addition to the successive deposition and successive-etching procedures used for residual stress extraction, the production line should properly incorporate online measurements with quick and reliable residual stress evaluation.