## Abstract

Direct wafer bonding is a manufacturing process that is used in the fabrication of electronic, optical and mechanical microsystems. The initial step in the process requires that the wafers are sufficiently smooth, flat and compliant such that short-range surface forces can elastically deform the wafers and bring the surfaces into complete contact. Analytical and computational mechanics models of this adhesion process as well as experiments that validate these models are presented in this work. An energy-based analysis is used to develop the models that allow acceptable limits of wafer-scale flatness variations to be predicted. The analytical models provide basic insight into the process while the finite-element-based model reported provides a method to analyse a broad range of realistic cases, including the bonding of wafers with etch patterns and arbitrary geometries. Experiments, in which patterned silicon wafers with different magnitudes of wafer-scale shape variations were bonded, were performed and the results demonstrate that the proposed models can accurately predict the size and shape of the bonded area based on the wafer geometry, elastic properties and work of adhesion.

## 1. Introduction

Direct wafer bonding is the joining of thin (0.4–1 mm), large diameter (50–300 mm) semiconductor substrates without the use of an intermediate bonding medium. A two-step procedure is used to join the wafers, in which the wafers are initially bonded at room temperature via macroscopically short-range surface forces and then heated to form strong primary bonds at the interface (Tong & Gosele 1999). The process is attractive, as it allows high-strength bonds that are stable at high temperatures to be formed without the introduction of a bonding interlayer or the application of clamping loads during thermal treatment. The numerous attributes of this bonding approach have resulted in it being employed in a range of applications, including the manufacture of semiconductor substrates, advanced micro- and optoelectronic devices and microelectromechanical systems (MEMS).

While direct bonding continues to be adopted in an even wider range of applications, a lack of quantitative models has limited the ability to design and control the process. This has made adoption of direct bonding for new applications a time-consuming and expensive endeavour and has restricted the use of pre-bond metrology. A large fraction of the manufacturing problems that arise result from incomplete bonding in the room-temperature contacting step. If the wafers are sufficiently smooth and flat and there is an adequate density of appropriate bonding species on the surfaces, a bond front will propagate outward from the point of initial contact and pull the wafers together. However, if the wafers have high surface roughness or different shapes due to large flatness variations, bonding may fail. Establishing the acceptable limits of surface roughness and flatness variations is essential to develop robust bonding processes.

Several authors have reported models that examine the role of surface waviness (Tong & Gösele 1995; Yu & Suo 1998) and surface roughness (Gui *et al*. 1999; Miki & Spearing 2003) in the room-temperature bonding step. The models that describe the effect of surface waviness use an energy-balance approach that compare the elastic strain energy required for the two wafers to deform to a common geometry to the reduction in surface energy that occurs when the bond is formed. The analyses given by Tong & Gösele (1995) and Yu & Suo (1998) provide insight into the importance of wafer geometry and give a first-order criterion for bonding. However, since they compare the reduction in total surface energy to the total gain in strain energy, the models cannot be used to predict the degree to which wafers will bond or the effect of local flatness irregularities or etch patterns. In recent work, Turner & Spearing (2002) proposed a bonding criterion that allows the extent of bond advance to be predicted by considering the energy accumulated in the wafers as the bond propagates. This analysis, which is similar to the familiar Griffith energy balance in fracture mechanics and the approach employed by Johnson *et al*. (1971) in their work quantifying the role of adhesion in the contact of compliant spheres, allows the size of the bonded area to be predicted and permits bonding difficulty to be evaluated locally across a wafer. Turner & Spearing (2002) presented the basic framework, analytical expressions and two-dimensional finite-element simulations to illustrate the use of this approach for assessing the bonding of bare and patterned wafers in which the out-of-plane shape of the wafers are defined by a simple curvature. The results presented illustrate the essentials of the modelling approach, but did not provide a route for applying the model to realistic wafers with arbitrary geometries or verify the models experimentally.

The current work presents computational mechanics models that extend the models developed in Turner & Spearing (2002) to a wider range of cases and reports experimental results that validate the models. The paper begins with a review of the general energy-based bonding criterion and the analytical models that were developed using this criterion, as they provide background for the development of the computational model and provide insight into the design of the experiments. Next, a finite-element-based model that uses techniques adapted from fracture mechanics to evaluate the bonding criterion and calculate the shape of the bond front in wafer pairs with arbitrary geometries is presented. The model is outlined and its capabilities demonstrated through a set of test cases. Then, experiments, which were conducted to test the validity of the model, are described. Finally, the model and experiments are compared and the results are discussed.

## 2. Mechanics framework

Wafers typically have a range of flatness variations, from surface roughness at the smallest scale, to surface waviness with millimetre-scale spatial wavelengths, to shape variations that extend across the entire wafer. All of these flatness variations can play a role in bonding and must be controlled in order to realize robust bonding processes. The focus in this work is the effect of wafer-scale shape variations, which are deviations of the midsurface of the wafer from a flat reference plane, as shown in figure 1*a*. These shape deviations can result from manufacturing variations in wafer production as well as the deposition of residually stressed thin films on the surface. Shape maps of two bare 100 mm diameter silicon wafers that are from the same manufacturing lot are shown in figure 1*b* and demonstrate the type of shape variations that are commonly observed on prime grade semiconductor wafers. While the wafers shown in figure 1*b* have moderate flatness variations with a peak-to-valley height of 10 μm, larger variations, up to 100 μm across a 100 mm wafer, are not uncommon on poorly manufactured wafers or wafers with a highly stressed film on one surface.

In order for bonding to be achieved, wafer-scale shape variations, such as those shown in figure 1, must be accommodated during the room-temperature bonding step through elastic deformation. This combination of adhesion and elastic deformation that results in bonding is not relegated to just direct bonding, but is also found in a range of natural and engineered systems. Adhesion that is mediated by short-range surface forces plays a role in the contact between soft solids (Johnson *et al*. 1971), the ability of certain animals to adhere to walls (Spolenak *et al*. 2005) and stiction problems in MEMS (de Boer *et al*. 2000). In all of these examples, an energy-based approach, in which the change in strain energy in the body is compared to the change in surface energy, has been used to model the adhesion processes. This approach, which was first proposed by Johnson *et al*. (1971), and later demonstrated to be equivalent to the Griffith energy balance in linear elastic fracture mechanics (Maugis & Barquins 1978), is well suited to examine the direct bonding process.

The application of this energy-based approach to direct bonding can be understood by considering the model bonding system illustrated in figure 1*a*. Two wafers that have different initial shapes must elastically deform to a common shape for bonding to occur. The bonding process is driven by short-range surface forces. The overall driving force may be quantified in terms of the surface energies, *γ*_{1} and *γ*_{2}, and the interface energy, *γ*_{12}, as the work of adhesion, . As the bonded area increases, the strain energy, *U*, and total interface energy increases, while the total surface energy decreases. Calculating the total system energy and differentiating to find the equilibrium bond area of the system yields the criterion for bond front advance (Turner & Spearing 2002)(2.1)The term d*U*/d*A*, which is a function of the geometry and elastic properties of the wafers, represents the strain energy accumulated in the wafers as the bond advances. While it may be helpful to think about d*U*/d*A* as a strain energy accumulation rate, the quantity is equivalent to the strain energy release rate, *G*, that is commonly used in fracture mechanics. The fact that d*U*/d*A*=*G* is important, as it allows numerous analytical and computational methods that have been developed for fracture mechanics to be employed for adhesion problems.

To employ the criterion given in equation (2.1) to predict success or failure in bonding processes, *G* must be calculated based on the wafer geometry and elastic properties and then compared to the work of adhesion. Routes for calculating *G* analytically and numerically to examine the effect of wafer-scale shape variations are demonstrated in §§3 and 4.

## 3. Analytical model

If the shapes of the initial wafers can be described by simple functions, analytical models can be developed to predict the final shape of the bonded pair and the extent of bond advance. A relevant simple shape to consider, as it is the shape that results when there is a uniform residually stressed film on one surface of a wafer, is a wafer which is parabolic with a shape described by a curvature, *κ*. The position of the midsurfaces of the two wafers, *w*, as a function of the radial coordinate, *r*, is given as(3.1)

The application of the bonding criterion given in §2 to this case was examined in Turner & Spearing (2002). The analysis used classical plate theory and gives the final curvature of the bonded pair,(3.2)and an expression for the strain energy release rate when bonding unpatterned wafers,(3.3)where *η* is the thickness ratio, *η*=*h*_{1}/*h*_{2}, *Σ* is the modulus ratio, *Σ*=*E*_{1}/*E*_{2} and *R* is the non-dimensional bond front radius, *R*=*c*/*b* (*c*, bond radius; *b*, radius of wafer). It is important to note that when the wafers are unpatterned and have a parabolic shape, d*G*/d*R* is always negative, suggesting that bond propagation will be unstable if the work of adhesion is constant across the wafer. Thus, wafers will either bond completely or not bond at all. As such, the criterion for bonding may be written by taking *R*=0(3.4)

In many practical applications, features are etched on the surfaces of the wafers prior to bonding and these features can affect the bond propagation. In Turner & Spearing (2002), the effects of shallow and deep features on bond propagation were considered. Here, the results for bonding with shallow features are reviewed. The simple pattern, shown in figure 2, was considered and is used in the experiments here. The strain energy release rate for a wafer pair with this spoke pattern at the interface is(3.5)

Depending on the pattern, bond propagation in patterned wafers may be stable or unstable. As an example, for the spoke pattern used in the experiments in this work, where *R*_{0}=0.14 and *C*_{t}=0.30, d*G*/d*R*>0 for *R* from *R*_{0} to 1, the bond front advance will be stable. Given that the propagation is stable, the bond front will arrest when *G*=*W*, which allows the radius of the bonded area to be calculated. For the pattern used in this work and assuming a Poisson ratio of 0.18, which is an in-plane average for a (100) silicon wafer, the bonded radius is(3.6)

Note that the size of the bonded area is a function of the work of adhesion, elastic modulus, wafer thickness and curvature difference.

## 4. Numerical model

The analytical model in §3 provides insight into the essential scaling in the bonding process and can be used to predict the final shape and bonded areas in pairs in which the wafers are elastically isotropic and have simple geometries and etch patterns. For more complicated cases, where the wafer shapes are not described by simple functions, the wafers are single crystal and have anisotropic elastic properties, or the etch patterns are complex, a more flexible model is required. The challenge in predicting the bonded area in these more complicated cases is that the shape of the bond front is not known *a priori* and must be determined by finding the bond area where the bonding criterion, *G*=*W*, is satisfied along the entire front. This is significantly more challenging than the axisymmetric case in which the bonded area was assumed to be circular and defined by a single bond radius. In order to examine this more general class of problems, a finite-element-based numerical model that implements the bonding criterion outlined in §2 was developed. Finite-element methods were used as they allow complex geometries and material properties to be defined and *G* to be calculated locally along the bond front. The finite-element model was coupled with a procedure to iterate the bond front position that allowed the size and shape of the bond area to be predicted.

### (a) Finite-element model

The wafers used in the experimental work, and hence used to define the geometry and elastic properties in the model, were 100 mm (100) silicon wafers. The mesh employed is shown in figure 3*a*. The wafer has a 100 mm diameter and a flat (18 mm long) that is oriented parallel to the [110] direction. The spoke pattern used in the experiments (*R*_{0}=0.14 and *C*_{t}=0.30) was incorporated into the mesh. The wafers were modelled as a linear elastic material with anisotropic properties (silicon: C_{11}=165.7 GPa, C_{12}=63.9 GPa, C_{44}=79.6 GPa). The out-of-plane geometry was defined by the wafer thickness and wafer shape (position of the wafer midsurface). For analysis of the experiments, the thickness and shape were imported directly from measurements of the wafers. The model includes both wafers in the bonding pair and the shape and the thickness of each are specified independently.

Figure 3*a* shows a mesh with a circular bond front, as was often used as the initial condition, while figure 3*b* shows a mesh after iteration in which the bonding criterion is satisfied along the entire bond front. The basic arrangement of the elements in plane shown in figure 3 was used throughout this study. There are four elements across the width of each spoke and three across the regions between the spokes. The mesh is refined around the bond front in the radial direction. The refined region is 5 mm wide and has four elements perpendicular to the front. The model was meshed with 20-node continuum solid elements with two elements through the thickness of each wafer. The mesh geometry shown was chosen through a series of convergence studies that ensured that the *G* values calculated were not a function of the mesh density.

The bonding process was modelled by calculating the distance between the two wafer surfaces from the shape and thickness data and then applying an appropriate displacement to close the gap at the interface. The displacements were applied between corresponding nodes on the opposite wafer surfaces and were specified in a way such that the gap was closed while maintaining a constant reaction force on corresponding nodes on each wafer. This is required as only the total displacement that is needed to close the gap is known and the relative deflections of the two wafers is governed by equilibrium considerations. The mesh was generated using a custom Matlab script and was exported to and solved with the commercial finite-element package Abaqus.

### (b) Local calculation of G

The quantity d*U*/d*A* that is compared to the work of adhesion to assess bonding is equivalent to the strain energy release rate. A host of techniques, including the elemental crack advance method, the *J*-integral and the virtual crack closure technique (VCCT), have been developed for calculating the strain energy release rate using finite-element methods. While any of these techniques can be adapted for adhesion analyses, the VCCT is used here as the implementation is straightforward and it allows *G* to be calculated locally along the entire bond front from one finite-element simulation. The basic principle of the VCCT is that the strain energy required to advance a crack may be estimated by calculating the work required to close the crack over a small area immediately ahead of the crack front. The work is calculated from the reaction forces at the crack tip and the displacements immediately behind the crack tip. *G* was calculated locally along the front using the standard equations for 20-node solid elements (Krueger 2002). The VCCT provides the ability to calculate the mode I and II components of the strain energy release rate. In cases where dissimilar materials are bonded, or the wafers have significantly different thicknesses, the bond will be subject to a mixed mode loading that may affect the behaviour and both components must be considered. However, in the present work, in which both wafers are silicon and have nearly the same thickness, the loading at the interface is nearly pure mode I and the mode II component does not need to be considered.

### (c) Iterative scheme

The equilibrium bond position is the position at which *G*=*W* is satisfied along the entire bond front. As the shape of the bond front is not known *a priori*, a bond front position is initially assumed. The finite-element model is solved and it is checked to see if the front is in equilibrium. If the bond shape does not satisfy equilibrium, the position of the bond front is adjusted and equilibrium is checked again. The adjustment to the bond front position is done in a systematic fashion using a quasi-Newton method with a rank-one Broyden update scheme (Hitchings *et al*. 1996).

The strain energy release rate at any node on the bond front can be considered to be a nonlinear function of the position of all the nodes along the bond front. The mesh is constructed in a manner such that the position of the nodes on the front is defined solely as a function of the radial coordinates of the bond front. This allows the strain energy release rate at any node to be defined as(4.1)where *i* indicates the node and *n* represents the total number of nodes along the front. Equilibrium is satisfied along the entire front when the strain energy release rate at every node along the front is equal to the work of adhesion. This condition may be written as a nonlinear set of equations(4.2)

The vector, * R*, (

*=[*

**R***r*

_{1}…

*r*

_{n}]

^{T}), fully specifies the bond front position. The strain energy release rates at each of the nodes may be calculated for a given

*from one finite-element simulation.*

**R**The * R* that satisfies equation (4.2) specifies the bond front position at equilibrium and may be solved iteratively. The approach used here to solve equation (4.2) was adapted from Hitchings

*et al*. (1996) and uses a quasi-Newton method to update the bond shape after each iteration. After each iteration, a new bond front position is calculated:(4.3)where

*k*indicates the iteration number. An approximate Jacobian,

*, is updated after each iteration using the rank-one Broyden update formula(4.4)where(4.5)*

**B**The iteration scheme calculates the new bond front position based on the previous position of the front, the difference between *G* and the specified work of adhesion and the approximate Jacobian. Given that *G* may vary significantly along the bond front and the exact Jacobian is not being used, equation (4.3) may predict a new bond front that is not physically possible (i.e. outside the boundary of the wafer) or would cause significant mesh distortion. To prevent this, the radial shifts, (Δ* R*=

**R**^{k+1}−

**R**^{k}), that are predicted in each iteration, are checked to ensure that they do not exceed a specified value, Δ

*R*

_{lim}. Values of 3–5 mm for Δ

*R*

_{lim}were used. If any of the calculated shifts along the front did exceed the specified Δ

*R*

_{lim}, the shifts along the entire front were scaled by a factor, , such that the maximum shift at any node along the front was Δ

*R*

_{lim}. The bond front position was iterated until

*G*was within ±1% of the specified work of adhesion value along the entire front. Convergence was typically achieved in 20–50 iterations. This iterative scheme and its integration with the other parts of the model is summarized in figure 4.

### (d) Results

The capabilities of the model were explored by examining the bonding of different combinations of wafers with the shapes shown in figure 5. The bonding of wafers B and E shown in figure 5 illustrates the basic operation of the model well. The shapes shown in the figure were imported into the model and an initial circular front with a radius of 27 mm was assumed (figure 3*a*). Figure 6 shows the distribution of *G* along the circular front. The variation in *G* that is observed results from the asymmetry in the shape of the wafers and the anisotropic elastic properties. The model required 32 iterations to determine the bond front in which *G* was within ±1% of the target of 25 mJ m^{−2}. The distribution of *G* along the converged bond front is shown in figure 6 and the shape of the equilibrium bond front is shown in figure 3*b*.

The bonded area shapes for six different wafer pairs at various values of work of adhesion are shown in figure 7. Each pair shown in figure 7 is labelled by a two-letter combination which indicates the shapes of the two wafers shown in figure 5 that were considered. The predicted fronts for the three pairs shown in the top row of figure 7 illustrate the effect of varying degrees of asymmetry on the bond front shape. Wafer A has the shape with strongest asymmetry and thus the bond front predicted for pair A–E is the most asymmetric in the set. The asymmetry of the underlying wafers decreases from left to right in the top row of figure 7 and a corresponding change in the shape of the bond front is seen. Comparison of the size of the bonded areas at a constant work of adhesion (20 or 25 mJ m^{−2}) across pairs A–E, B–E and C–E shows that the bonded area is smallest in pair C–E. As all the wafers have the same maximum curvature and pair C–E has the least asymmetry, this result suggests that a conservative estimate of the bond area could be made by using the axisymmetric analytical model with the maximum curvature difference.

Comparing pairs B–E and C–E to pairs B–D and C–D illustrates that the shapes of both wafers in the pair influence the shape of the final bonded area. It is seen that in pairs B–E and C–E, where the bottom wafer has an axisymmetric curvature that is comparable to the maximum curvature of the asymmetric wafer, the bonded areas are less asymmetric than when the bottom wafer is flat (pairs B–D and C–D). Pair B–F shows the bonding of two asymmetric wafers with their principal curvatures orthogonal to one another. It is seen that the bonded area is not significantly asymmetric, which can be understood by realizing that the gap at the interface (*w*_{B}−*w*_{F}) is axisymmetric. The final point that is illustrated in the predictions shown in figure 7 is the effect of elastic anisotropy. In all the predictions, the bond areas have sides that are slightly flattened as a result of the cubic symmetry of the silicon.

## 5. Experimental

### (a) Approach

The aim of the experiments in this work was to validate the mechanics-based models for wafer bonding presented. Figure 8 shows the basic arrangement used in the experiments. An unpatterned wafer with a tensile residually stressed film on the back surface was bonded to a wafer with a shallow etched spoke pattern on the bonding surface. A spoke pattern with *R*_{0}=0.14 and *C*_{t}=0.30 was used to achieve stable propagation such that the bond front would arrest when *G*=*W* and result in a measurable bond area.

The basic premise of the experiments is to bond wafers with different shapes, monitor the propagation and assess the ability of the model to predict the final shape and size of the bond front. While the elastic properties of the wafers are known and the geometry and thickness was measured, the work of adhesion was not known and could not readily be measured. The work of adhesion in this context is a lumped parameter that is a function of the type and density of the bonds formed at the interface as well as smaller-scale flatness variations (waviness and roughness). Given that the work of adhesion could not be measured independently, the approach used in these experiments was to bond multiple pairs with different shapes, calculate the work of adhesion for each pair from the bond position using the model and then compare the work of adhesion across different pairs. Given that the wafers were from the same manufacturing lot and were cleaned in the same manner prior to bonding, the expectation is that the extracted work of adhesion values should be constant across all the pairs if the model properly accounts for the different geometries of the wafers.

All wafers used in the bond propagation study were 100 mm (100) silicon wafers that were nominally 525 μm thick. The wafers were single side polished and were taken from the same manufacturing lot. A tensile residually stressed nitride film was deposited on the back surface of one wafer in the pair to add controlled wafer-scale shape variations. The thickness of the film was varied from 0.2 to 0.6 μm to achieve a range of wafer shapes. These wafers were fabricated by initially growing a 500 nm thermal oxide, then deposition of high stress (*σ*≈0.7–1.0 GPa) stoichiometric silicon nitride. The silicon nitride was removed from the bonding surface by masking the back with a deposited low-temperature oxide, and etching the nitride from the front in a standard hot-phosphoric etch solution (H_{3}PO_{4} (85%) at 165 °C). The mask oxide on the back surface and the 500 nm thermal oxide on the bonding surface were removed using a buffered oxide etch (7 : 1 NH_{4}F(40%) : HF(49%)). The spoke patterned wafers were fabricated by first growing a 1000 nm thermal oxide, then patterning the oxide using standard photolithography and a wet oxide etch. The pattern was transferred to the silicon using deep reactive ion etching. The etch depth on all wafers was between 2 and 3 μm, which is deep enough to prevent bonding, but sufficiently shallow not to alter the stiffness of the wafer. Following etching, the oxide mask was removed using a buffered oxide etch. A standard RCA clean was used immediately prior to bonding to ensure a hydrophilic surface and consistent surface chemistry.

During bonding, the wafers were supported at the centre on a 3 mm diameter pin and pressed into contact with another pin to initiate the bond. All experiments were performed at room temperature and pressure in a class 100 clean room at a relative humidity of approximately 37%. The bond front propagation was monitored using infrared (IR) transmission imaging and the bonded areas were measured directly from digitized images. After bonding, the pairs were separated and the geometry of each wafer was measured using an ADE 9900 capacitance gauge (ADE Corporation 2000) that has a 2 mm^{2} capacitance probe that is scanned across the wafer to obtain full wafer maps of the thickness and shape.

### (b) Analysis

The work of adhesion for each bonded pair was calculated based on the wafer geometry, elastic properties, etch pattern and measured bond area using the axisymmetric model described in §3 and the numerical model presented in §4. Employing the axisymmetric model (equation (3.6)) was not straightforward as the shapes of the wafers and the bonded areas for the most part were not axisymmetric. This made it difficult to define a single bond radius and wafer curvature difference for each pair difficult. As a result, for each pair, a minimum, average, and maximum bond radius and curvature difference were determined from the experiments and a minimum, average, and maximum *W* for each wafer pair were calculated using equation (3.6). The finite-element model presented a more direct route to calculate work of adhesion values for the pairs as their asymmetric geometries could be accounted for properly in the model. For each pair, the measured shapes and thicknesses of the individual wafers were imported into the model and the bonding behaviour was mapped out as a function of work of adhesion.

Figure 9 shows an example of this process for one of the pairs in the experiment. The measured shapes of the individual wafers that are imported into the model are shown in figure 9*a,b* and the bond fronts predicted using the finite-element model for a range of work of adhesion values are shown in figure 9*c*. The work of adhesion for each pair may be determined by selecting the predicted bond front that fits the experimentally observed front the best.

### (c) Results and discussion

Six wafer pairs were bonded in the study. The shapes of the individual wafers in each pair and the resulting bond radii are summarized in table 1. The wafer shape has been expressed as the peak-to-valley shape difference (P−V) across the centre 90 mm of the wafer. The sign next to the P−V values indicates if the wafer shape was concave up (+) or down (−) during bonding. The signed P−V value clearly does not fully describe the wafer shape, but is given as it provides an understanding of the relative shapes of the wafers in the experiments. Measurements of the wafer thickness showed that the mean thicknesses of the wafers in the set ranged from 523 to 528 μm and that the total thickness variation of each wafer was less than 4 μm.

The bond front position (and hence the work of adhesion) was observed in the experiments to be time dependent. The bond front was observed to advance at a rate of several millimetres per second in the first few seconds after contact and then continue to grow at a much slower rate for hours after contact. This time dependence is presumably associated with the dynamics of bond formation at the interface as well as the adsorption of molecules (i.e. water) on the surface and is likely to be a function of various factors such as temperature and humidity. No attempt was made to characterize this dependence fully or to further understand its origin as it is outside the focus of validating the mechanics models in this work. To minimize the effect of time in the extracted work of adhesion values, all images and measurements of the bond front were recorded 10 min after initial contact of the wafers. A period of 10 min was selected as the rate of bond advance was relatively slow at this point and the time period was short enough such that multiple pairs could be bonded sequentially without the ambient conditions (temperature and humidity) changing significantly.

The IR images of four of the bonded pairs listed in table 1 are shown in figure 10. Also shown in the images are the bond fronts predicted using the finite-element model. In each image there are two predictions, the smaller of the predicted bond areas corresponds to *W*=35 mJ m^{−2} and the larger to *W*=40 mJ m^{−2}. Pair P5 is not shown as the bond front propagated to the edge of the wafer. The sixth pair in the experimental set is not shown either as irregular bond front propagation, which is consistent with surface contamination, was observed.

The work of adhesion values extracted from the pairs in the bond propagation study is summarized in figure 11. Work of adhesion values calculated via the axisymmetric model and the finite-element model are both shown. For the axisymmetric analytical results, the markers in figure 11 correspond to the *W* calculated from the average bond radii and curvature difference of the pair and the limits of the error bars represent the *W* calculated from the minimum and maximum bond radii and curvature. Across pairs P1–P4, the average work of adhesion values extracted using the axisymmetric model range from 29 to 43 mJ m^{−2}, while those extracted with the numerical model range from 36 to 40 mJ m^{−2}. The variability in the results extracted with the axisymmetric model is not surprising given that clear asymmetries are present in the wafer shapes and observed bond fronts. Pair P5 is the pair that bonded completely and as such only a lower bound on the work of adhesion may be calculated from the test. The lower bound calculated from P5 is consistent with the work of adhesion values of pairs P1–P4. The irregular propagation observed in pair P6 prevented the model from being employed to compute a work of adhesion.

The results in figures 10 and 11 show that the finite-element model accurately predicts the shape of the bond area and that the work of adhesion values extracted from pairs with different geometries are the same. The magnitude of work of adhesion values extracted is reasonable. The work of adhesion values reported are ‘apparent’ work of adhesion values in that they include effects of smaller-scale flatness variations as well factors such as humidity. Previous measurements on polycrystalline and single crystal silicon surfaces with a thin native oxide have shown that work of adhesion values can range from 10 to 140 mJ m^{−2} (e.g. Miki & Spearing 2003; Maboudian & Carraro 2004), depending on the roughness and the relative humidity. The work of adhesion values of 36–40 mJ m^{−2} measured in this work fall in this range.

## 6. Conclusions

A comprehensive set of models for predicting bonded area based on the wafer geometry, elastic properties and work of adhesion have been presented and experimentally validated. The analytical model is useful for analysing wafers where the wafer-scale shape variations are parabolic and to make first-order estimates of bonding in asymmetric wafer pairs. The numerical model presented is significantly more flexible and allows the bonding of wafers with complicated shapes and etch patterns to be analysed. The excellent agreement between the predicted and observed bond areas shown in figure 10 and the consistency in extracted work of adhesion values (numerical results) shown in figure 11 indicate that the proposed model accurately captures the mechanics of the bonding process and has the ability to model the effect of wafer-scale shape variations in direct bonding.

The set of models presented here have significant potential to aid in the design and control of direct bonding processes. The analytical model provides a route for quickly estimating the ability to bond two wafers and clearly shows the effect of changing the geometry or material of the wafers in the process. The numerical model allows tolerances to be set on wafer geometry and enables the effect of specific etch patterns on bonding to be studied. Furthermore, the models presented can be used to aid in the design of novel wafer bonding processes, such as the bonding of wafers in a stressed configuration to produce uniaxially strained silicon.

Beyond wafer bonding, this work also demonstrates a unique approach of using computational fracture mechanics techniques to analyse surface-force-driven adhesion problems that arise when soft materials or thin compliant structures make contact. As investigations into contact in biological and micromechanical systems are pushed further, there will undoubtedly be a need to employ computational approaches to adequately understand the mechanics of the systems. The computational techniques demonstrated here are general and provide a set of tools that allow such problems to be addressed.

## Acknowledgments

The authors would like to acknowledge the MIT Microsystems Technology Laboratory (MTL), where all of the microfabrication and experimental work was performed, and the assistance of the ADE Corporation (P. Hester) in obtaining wafer shape measurements. This work was supported by the Cambridge MIT Institute (CMI-059/P-IR(FT) MEMS).

## Footnotes

↵* Author and address for correspondence: Department of Mechanical Engineering, University of Wisconsin–Madison, 1513 University Ave, Madison, WI 53719, USA (kturner@engr.wisc.edu).

↵† Present address: School of Engineering Sciences, University of Southampton, Southampton SO17 1BJ, UK (spearing@soton.ac.uk).

- Received April 15, 2005.
- Accepted August 18, 2005.

- © 2005 The Royal Society