Modeling of stress and strain fields induced during the Smart-Cut process on silicon. Influence of different couplings for diffusion of hydrogen at a microscopic scale

The Smart-Cut technology consists in the increasing of external temperature and pressure imposed by the diffusion of hydrogen ions within silicon substrates leading to a wafer splitting. In the present work, we studied the evolution of the stress field in the crystalline lattice of silicon, the diffusion of hydrogen ions that can lead to the growth and coalescence of cavities. Meanwhile, we test several models and simulate these phenomena by a numerical approach, in order to compare the results to experimental observations.


Introduction.
The Smart-Cut technology is an industrial process widely used in the manufacturing field of microelectronic components with the semiconductor materials, to produce high quality silicon-onisolator (SOI) in high volume [1]. This technology takes place in the background of improving and optimizing the production methods, allowing for a more efficient and more economic process, and providing for a more controllable finished product quality. Because of the great interests of Smart-Cut technology, it is necessary to investigate the basic mechanisms underlying the thin wafer separation phenomena during the annealing treatment after hydrogen implantation [1,2]. Furthermore, there is insufficient knowledge about the formation of penny-shape microcracks and about the mechanisms of substrate fracture, as well as its spreading occurring in the high hydrogen ion concentration layer. Consequently, a solution can be found through the use of multiscales modeling and corresponding numerical simulations of such processes [3,4]. Therefore, we are presently interested in understanding theoretically the mechanisms of those phenomena, especially the behaviors in the substrate during the first stage of the splitting process [5,6]. In this paper, we studied particularly different geometrical parameters and different initial conditions imposed on hydrogen penny-shape microcracks. We implemented this model in finite element analysis software (Abaqus) to calculate the concentration fields of hydrogen. Finally, a weakly coupled model (taking into account simultaneously diffusion and mechanical effects and its coupling) is assembled to simulate the phenomena during the first stage of the splitting process. Total duration is between 1 and 3 hours. Thus, we manage to calculate with numerical methods the mechanical stresses developing around hydrogen penny-shape microcracks/defects and the geometric evolution of those defects in the wafer.

Elements of the modeling.
Physical aspects. The modeling is composed of two parts: a diffusional aspect and a mechanical aspect. The system is composed of two phases: the substrate of silicon Si and the gas-filled H 2 penny-shape microcracks. It is assumed that no silane ( y x H S ) can form within the system. Only the mechanical behavior of the silicon is really investigated. It means that the gas inside the pennyshape microcracks is described simply using a pressure and a temperature. Temperature is supposed homogeneous in all the system, whereas pressure is not necessarily homogeneous especially along the gas/silicon interface. Moreover, the latter is considered with a very slow speed (quasi-static domain). Mechanical behaviors. Description of the mechanical phenomena takes place at two separate scales. At microscopic scale, diffusion of hydrogen interacts with mechanical behavior of Si leading to an evolution of penny-shape microcracks. This effect can include both growth and coalescence of the defects, which should be treated through the continuous damage mechanics. At macroscopic scale, the results of the coalescence of penny-shape microcracks along a preferential direction lead to the propagation of cracks. In the present article, we focus only on the microscopic scale. The silicon has a ductile-brittle transition with temperature [7]. For the considered temperature (700K) during the annealing treatment, the silicon is ductile. Ductile damage meaning coupling with elastoplasticity should thus be taken into account [8]. However, for the first stage of the splitting process, it is supposed that hardening and damage are not too significant. This strong assumption is used to simplify the modeling, because we focus on the study of the coupling between mechanics and diffusion at the early stage of the splitting process. In other words, the silicon is elastic and perfectly plastic. For elasticity, an anisotropic behavior is considered since the substrate is a monocristal of orientation [100]. Moreover, according to Varshni et al. [9], the influence of temperature on the different components of the silicon stiffness tensor can be calculated through the following equation: (1) where T is the temperature, 2 1 , , a a C ij ref are constants given in the crystallographic frame, whose numerical values are given in [9]. Secondly, for the considered orientation [100], the elastic modulus has to be calculated in the laboratory frame [10], which can be done with: (2) The same kind of relations can be given for the Poisson ratio [11]. Eventually, when the calculated stress will reach a threshold, then plasticity is supposed to take place and the present simulations (without hardening and damage) will be no longer valid. At 773K, Nakao et al. [7] evaluate the silicon yield stress of 1,68 GPa, whereas the fracture stress is only of 1.78 GPa. Because of the gas within the penny-shape microcracks, a pressure occurs along the gas/silicon interface that acts on the silicon as an internal boundary condition. As a first approximation, this pressure could be calculated using the perfect gas state law. For the volume of the penny-shape microcracks (initially around 8 nm of diameter), for a pressure at the yield stress and for a temperature of 700K, we obtain a number of hydrogen molecules H 2 of 4 10 9 . 8 × . It leads to a hydrogen concentration of . This result is absolutely unrealistic. Strictly, because of the substrate, Feng et al. [12] propose a method to calculate the critical concentration to obtain the fracture of the penny-shape microcracks based on the perfect gas law and taking into account elasticity of the substrate. The calculation leads to an hydrogen concentration of 3

/ 109
nm H and the corresponding density is 3 / 368 m kg . Once again, at this density, the gas is no more "perfect" because of the frequent interactions between hydrogen molecules. Thus, we propose to consider a Van der Waals gas given by equation:

708
Residual Stresses IX where A N is the Avogadro constant, R is the gas constant. V is the volume of the penny-shape microcracks. P is the internal pressure of the penny-shape microcracks. . By extending the model of Feng et al. [12], and using now Eq. 3, a more complex relation has been determined: where γ is the critical surface energy needed to fracture the silicon (based on Griffith criteria); h is the height of the penny-shape microcracks considering it as a cylinder. It is interesting to note that Eq. 4 depends on the radius of this penny-shape microcrack only through its volume. It leads to a critical value for the concentration of hydrogen molecules concentration of , which satisfy the chosen assumptions. Diffusional behaviors. Description of the diffusional phenomena is based on Fick law. Especially, it is necessary to take into account the Ostwald effect as observed by different authors [13,14]. Indeed, the growth of cavities corresponds to the transport of hydrogen from penny-shape microcracks to others. Big microcracks grow using indirectly small ones. As the number of hydrogen molecules increases, the pressure increases and opens the penny-shape microcracks (regulated by the volume increase). If two penny-shape microcracks meet, their volume merges, corresponding to a coalescence phenomenon. This is the mechanism proposed by Ostwald. Different approaches have suggested in order modeling this effect (phase fields method [15]…). In the present paper, we directly solve the Fick law: where H C is the concentration of hydrogen atoms in the silicon; ) (T D is the diffusion coefficient of hydrogen in the substrate, depending on the temperature. The Fick law is supposed to be applied on a 2D plane where the failure takes place. Based on the random walk model, the diffusion coefficient can be calculated according to: where nm a Si 543 . 0 = is the lattice parameter of the silicon (diamond-like structure); 1 13 0 10 − ≈ s ν is the jump frequency at 500K. The energy in Eq. 6 depends on the kind of considered defects. For interstitial migration, it is 0.45 eV. For lacunar migration, it is 0.55 eV. However, due to the lack of accuracy on the jump frequency, a mean value for the energy has been considered that leads to a value for the diffusion coefficient at 700K of 1  ), by using the Sievert relation [16], given by: where ) (T K is the reaction constant, which depends on temperature. Such an evolution is proposed by different authors based on the collisions of hydrogen [17]. The parameter H v is the dissociation . The reaction energy K E is evaluated to 0.29 eV. The parameter 0 K is a constant whose value is a priori 20.13 USI. Geometries. For the present article, two geometries are considered. It corresponds to 2D surfaces defined as 100nm*100nm squares. Such a choice is made because of the splitting process, which takes place only in a plane perpendicular to the free surface. On this surface, the cylinder pennyshape microcracks are circles. The number of circles on the surface, as well as their diameter (average of 8 nm), has been chosen according to the experimental observations of Personnic et al. [18]. Fig. 1 and Fig. 2 show the geometries. The first one has been chosen for the purely diffusional simulation (without taking into account of any mechanical coupling). When coupling with mechanics, calculation time increases drastically, so that we reduce the number of penny-shape microcracks, as shown in Fig. 2. We consider then only four penny-shape microcracks with respectively 4nm, 10nm, 20nm and 40nm diameters. It is supposed that each cylindrical penny-shape microcrack is h = 1nm high. The effect of the initial geometry distribution is not presented here.  Fig. 2 : Geometry for coupled simulation (diffusion with mechanics) Boundary and initial conditions. External and internal conditions have been defined for both the coupled (diffusion with mechanics) and non-coupled model (diffusion). On the external surface, the boundary condition is homogeneous, whatever the quantity (concentration or displacement). On the internal surface (corresponding to the gas/silicon interface), because of the diffusion phenomena, the boundary conditions are not a priori homogeneous. By using Eq. 7, after a first step, the pressure P due to hydrogen within the penny-shape microcracks is calculated. All the relations (3) and (7) have been implemented into Abaqus® FEM software using the subroutine VUINTER. The lateral surface S(t) and volume V(t) of this penny-shape microcracks are evaluated using the geometry evolution at each step. Then a calculation of the average pressure is done using: where "i" means for the ith node around the penny-shape microcracks. The volume and the average pressure are used to calculate the number of hydrogen molecules 2 H N by using Eq. 3 for the Van der Waals gas. According to the previous assumptions, the reaction of H is limited to a transformation in H 2 , and vice-versa. Moreover, the interface is considered in a quasi-static domain. Then the hydrogen flux around the penny-shape microcracks can be roughly calculated with: This average hydrogen flux is used as boundary condition for the next step of calculation. Initial conditions have been used to test either the influence of geometry or the influence of hydrogen 710 Residual Stresses IX concentration in the penny-shape microcracks (chosen around an average value of 10 atom/nm 3 ), whereas the hydrogen concentration in silicon is initially chosen at 2.45 atom/nm 3 . Numerical simulations. The simulations have been performed with Abaqus FEM software. Because it does not exist any predefined unit system in Abaqus, we have customized an adapted unit system, before launching the calculations. By considering the different natures of simulations, we defined the type of mesh element as linear standard (CPE3RT triangular strain plane element of Abaqus) for the simulation without coupling and as linear explicit coupled concentrationdisplacement for the simulation of coupled phenomena. Due to the irregular shape of the geometries, we choose free structured triangle elements for both simulations and set 10 seconds for both simulation step times. Results without coupling. By considering only the diffusional models defined by Eq. 5 and Eq. 6 with initial and boundary conditions, the results of the corresponding simulations are presented in Fig. 3 and   Fig. 3 presents the results of hydrogen concentration. We can observe that hydrogen concentration between two penny-shape microcracks with high initial concentration is higher than in other parts of the system, even if the flux is not significant. One can suggest that growth then coalescence of penny-shape microcracks may occur following such a preferential way (high concentration area), because of an increase of pressure in these areas. The penny-shape microcracks will then propagate to these areas so that the two penny-shape microcracks ends of these two areas merge together. Consequently, the coalescence effect would occur, what we try to check taking into account a mechanical coupling with diffusion. Fig. 4 presents the results of hydrogen flux fields. We can observe that this is more important between penny-shape microcracks of different concentration. Moreover, this flux is more important between penny-shape microcracks and the surface of the square domain. This is because of the difference of concentration. The higher is the difference of hydrogen concentration, the higher is the corresponding flux.
Results with coupling. By considering now the diffusional models coupled to the mechanical part, defined by Eqs. 1 to 9, the results of the corresponding simulations are presented in Fig. 5 and Fig 6. Fig. 5 presents the hydrogen concentration field for the coupled model. We can observe that concentration of area embedded by penny-shape microcracks is smaller than the initial concentration of other areas in the silicon. It means that hydrogen seems to be absorbed by the penny-shape microcracks all around. The bigger penny-shape microcrack absorbs more hydrogen; therefore the hydrogen concentration is smaller at the interface than for other gas/silicon interface. Fig. 6 presents the Von Misès equivalent stress in the substrate at the beginning of the splitting process. In red area, the stress is around 39 GPa, and is around 36 GPa in the orange area. In these two areas, this stress is higher than the yield stress or the fracture stress. We can thus imagine that the failure will grow in the red area, and then it will propagate along those areas. The failure would grow until the penny-shape microcracks merge all together.

Conclusions.
In the present work, investigations have been performed on the modeling and simulation of the Ostwald mechanism. It has been shown that correct concentrations in the penny-shape microcracks need the use of Van der Waals gas. Moreover, the diffusion is necessary to model this mechanism. The boundary condition at gas/silicon interface is of primary importance and leads to a complex modeling. The coupling with mechanics has been performed in the beginning of the splitting process, to see the corresponding stress (pressure, deviatoric stress…). Consequently, an evolution is seen with high residual stresses even at the beginning of the splitting process and indicating the need to take into account other relaxation mechanisms (plasticity…) in the early stage.