Influence of displacement constraints to the surface reconstruction of stressed bicrystal thin films

In this work, surface morphology evolution of bicrystal thin films under the combined action of grain boundary and surface diffusion is investigated by considering different mechanical constraints. 2D surface topographies of thin films, that are (a) freestanding, (b) strongly bonded to its substrate and (c) strongly bonded to its substrate and one of sidewalls, are simulated using a numerical implementation of an irreversible thermo-kinetics model. Relationships which give the groove depth as a function of time are obtained. Results show that mechanical loading conditions are effective in determining the morphology and kinetics of grooving. For the three scenarios that had been investigated, it was found that the groove depth evolves linearly with different tip velocities under the same level of uniaxial tension. In freestanding films groove tip evolves faster; i.e. as the film gets constrained from its substrate and/or one of its sidewalls, the tip velocity slows down. It was also observed that high triple junction mobilities at low levels of applied stress hinder the effects of displacement constraints to groove shape, even in the case of asymmetric stress distributions inside the film. On the other hand, low triple junction mobilities at moderate applied stresses allow formation of asymmetric grain boundary grooves due to the induced asymmetry in the driving force for surface diffusion with respect to the grain boundary.


Introduction
The reliable operation of the building blocks of modern technology, e.g. microelectronic and optoelectronic circuits, data storage devices and micro electromechanical systems depend on the morphology, property and operational behavior of polycrystalline thin metal films. These films are generally in bamboo polycrystalline structure (thickness of thin film less then grain size) and exposed to external or internal stress fields during production (deposition on a substrate, lithography, and encapsulation) and operation. Deformation processes through diffusional creep may operate in order to relax these stresses [1]. With time, flat thin film surfaces are self-reconstructed and reorganized with the action of stress gradient induced atomic migration (together with other driving forces) to include voids and hillocks of several shapes. In bamboo films, the most prominent evolution occurs at the regions around grain boundaries that intersect free film surfaces. The triple junction region, where the two grains and void phase meet, starts to form a groove/slit/wedge/pit/crack or even a hillock depending on the driving forces and active mass transport mechanisms. Undoubtedly, the state of stresses throughout the films is important although they are not structural components [2][3][4][5][6][7].
For a non-coherent film/substrate interface, film surfaces other than the sidewalls are traction free and the film may move freely on the substrate. Previous works of Thouless [8] and Genin et al [9] on bicrystal and polycrystalline films assumed this condition and couple capillary driven surface diffusion and stress driven grain boundary diffusion to model the effect of stress to grain boundary grooving. Diffusional flux was taken to be constant over the grain boundary and the thermal grooving problem with a source term was studied. Sekerka et al [10] generalized this model by removing the assumption regarding how the grain boundary flux is split onto the two adjacent sides of the grain boundary groove. In contrast to the pinning effect of thermal grooves, they Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
showed that the presence of a compressive stress can promote grain boundary migration. In a coherent film/ substrate interface, the interface imposes a mechanical constraint and diffusion trough the interface may not be possible due to strong binding. Gao et al [11] considered this type of mechanical boundary condition at the film/ substrate interface and developed a model for grain-boundary diffusion induced evolution of triple junction region. Assuming that the surface diffusion is sufficiently fast, they observed formation of crack-like grain boundary wedges. The model predicts an exponential decay of stress in the thin film. Later, Zhang and Gao [12] extended the model further by coupling surface and grain boundary diffusion in a way that one acting to other as matter source or sink. The other limiting case, where grain boundary diffusion is fast, as well as fully coupled case produced same results in terms of relaxation of grain boundary traction. In the aforementioned studies, stress distribution along the grain boundary is either approximated with functions that sustain a constant flux or assumed to be generated by a network of edge dislocations. Liu and Yu [13] on the other hand, studied the constrained film problem by coupling the surface and grain boundary diffusion by solving the elastic boundary value problem using boundary element method (BEM) to obtain the stress field, whose gradient acts as a driver for grain boundary diffusion. Since the stress field was calculated for evolving surface profile they observed significant changes in stress relaxation of thin films with every parameter change effecting the profile evolution such as the ratios of surface to grain boundary mobilities and energies. Liu and Yu [14] studied another set of mechanical boundary condition in which both substrate/film interface and one of the film sidewall have displacement constraints. This is rather an interesting case, because it leads to an asymmetric stress distribution around the grain boundary.
The crucial role of elastic dipole tensor interactions (EDTI) on driving force for surface diffusion was considered previously, under tensile [15] and compressive [16] uniaxial stresses applied to freestanding thin films. As the dimensions of ultrathin films approached, small thicknesses can lead to great stress gradients, which act as a driving force for grain boundary diffusion [17]. In this paper, grain boundary diffusion is first coupled to surface diffusion by considering the EDTI interactions. Then, the effects of applied stress system as a driving force for surface diffusion through EDTI on morphology evolution of bicrystal thin films, by the coupled mass transport mechanism, is investigated. The computational complexity of the calculation of stress distribution at the interior collocation points (other than the boundaries), to compute grain boundary fluxes, is overcome by extending the indirect boundary element method (IBEM) for the elastostatic problem. Finally, the influence of displacement constraints to the process is studied by considering films that are strongly bonded to their substrates and/or one of their sidewalls, besides freestanding ones. The results showed that the existence of displacement constraints might change the groove shape and kinetics.

Methodology
In a two dimensional model of a thin film with two grains of finite lengths are separated by a single vertical grain boundary, discrete set of calculation points or nodes can be grouped into two: ordinary points on the surface other than the triple junction (TJ) singularity and the TJ points. The velocity of ordinary points, V , ord the longitudinal TJ velocity, V , g lonḡ and the boundary conditions (BC) at the TJ in terms of right and left side surface fluxes, J , ō  may be described in terms of dimensionless, scaled and normalized quantities as follows [18,19]: In the above expressions, the bar sign over the symbols is used to indicate scaled and normalized quantities. l is the arc length in 2D space, k is the local curvature (positive for a concave surface), and s h is the hoop stress which is normalized with respect to the nominal stress applied at the sidewalls of the specimen, s . o The dimensionless parameter S corresponds to the relative intensity (referred to the capillarity) of the elastic strain energy density (ESED) contribution on the stress-driven surface drift diffusion. In the case of plain strain condition, ESED takes the following form: s e n s  - / Here s̲ ̲ and e̲ ̲ are the stress and strain tensors, where tensorial quantities are indicated by double bars under the symbols, n is the Poisson's ratio and E is the Young's modulus. Similarly, X designates the elastic dipole tensor interaction (EDTI) parameter, which measures the relative intensity of the elastic interaction between the elastic dipole tensor, ̲ ̲  and the surface stress system generated by the applied tractions and/or body forces acting on the bulk medium. For isotropic mobile defects, i.e. defects that have elastic dipole tensor with vanishing deviatoric parts  º Tr I 1 3 , where W a is the mean atomic volume of adatoms, and colon (:) designates the double inner product, and Tr corresponds to the trace of a tensor. In equations (2) and (3) d , ā and h ḡ are the interatomic distance and the thickness of the GB layer, respectively. W ḡ and W s are mean atomic volumes of chemical specie in the GB and the surface layer s, respectively. q + and qare right and left dihedral angles, M lonḡ and M trans are the longitudinal and transverse generalized mobilities of TJ.
is the wetting parameter which is assumed to be isotropic in the present study, where g g and g s are the specific surface Gibbs free energy densities associated with the GB and the surface layers, respectively. The normalized grain boundary diffusion flux (J ḡ ) due to stress field gradients elaborated in equation (3) is given by equation (4), where M ḡ is the normalized grain boundary mobility. The dependence of diffusion coefficients on stress is neglected.
Normalization and scaling are performed according to the following relations: The mobility terms described in above equations are all normalized with respect to surface mobility associated with the mass flow at the surface layer: Here h s is the thickness of the surface layer and D s is the surface diffusion coefficient.
The stress field is calculated by assuming plane strain condition and solving the elastostatic problem at each time step using the IBEM that employs straight constant line elements [20]. Neumann (i.e., traction boundary condition) BCs are utilized along the free (top) surface of the test modulus, and at the sidewalls at which constant uniaxial tensile stresses are applied. Traction free BC is imposed along the lower film surface (film/substrate interface) in the case of a freestanding film. For strongly bonded (displacement free) films to their substrates and/or sidewalls Drichlet BCs are imposed.
Interior collocation points are generated dynamically within each grain, as it evolves with time, using the grain center of mass as a reference. A similar approach was used to compute the temperature field inside a single crystal thin film with an internal void by Ogurtani and Akyildiz [21]. In order to illustrate the meshing strategy and verify the IBEM elasticity computations, the calculated displacements for a copper film (Young's modulus: E=129.8 GPa; Poisson's ratio: ν=0.343 [22]) subjected to biaxial tension are presented in comparison with that of commercial finite element method (FEM) code ABAQUS (figure 1). The imaginary film has a hole in its center, where the distribution of the nodal points and elements used in IBEM and FEM calculations are shown in figures 1(a), (b). These calculations show that the IBEM code developed approximates well to the displacements calculated by the commercial FEM code.
After calculating the stress field, the local curvatures required in equation (1) at each surface node is calculated discrete geometric relations and time integration of the evolution equations has been carried out using explicit Adams-Bashford formula of second order that employs an adaptive time stepping [18,23]. A shape preserving parametric cubic spline interpolation which keeps concavity of the data is employed to reset the node positions after each calculation step, so that a constant segment size is forced.  [22]. Assuming that the normalized hoop stress is equal to 1 (s h =1) at a fictive point on the surface of the thin film, this plot presents the relative contributions to equation (1). When the applied stress is around σ o =180 MPa (s = 1.0 ō ), one can see that EDTI contribution is about four orders of magnitude bigger that the ESED contribution. Keeping this in mind, the ESED contribution was omitted for the rest of the study. Figure 3(a) shows the initial configuration of a film exposed to uniaxial tension s = 0.745 ō applied from its sidewalls, which corresponds to an EDTI value of X = 0.1. This figure is particularly useful in order to visualize the applied mechanical boundary conditions. The components of displacement vector at each node are given in figure 3b for a freestanding film on its substrate. Since the film is centered at the origin and there is no constraint, upon pulling from both sidewalls, the upper film surface (node index i = [0, 60] in figures 3(b), (c) is displaced along the -y direction while the lower surface is displaced along the +y direction. Similarly, the upper film surface is displaced linearly along the -x direction from the left edge to the triple junction, where the displacement is zero, and then displaced along the +x direction from the triple junction to the right edge, and vice versa for the lower surface. The left and right sidewalls are displaced with the same amounts along the -x and long trans¯) as well as the grain boundary mobility (M ḡ ) are set to unity in these simulations to couple surface and grain boundary diffusion mass transport in these simulations. The groove takes the form of a narrow and deep grain boundary pit with fast tip evolution along the vertical grain boundary. Hoop stress distribution during this process around the TJ is presented in figure 4(b).

Results and discussion
The initial configuration of a strongly bonded film to its substrate is given in figure 5(a). The zero displacement constraint put on the film/substrate interface (node index i= [80, 140]) can be seen clearly in figures 5(a)-(c). The upper film surface is displaced in a similar way along the x-axis as in the previous case, but the displacement profile along the y-axis is now concave upward in the negative region rather than being constant. This in turn results in a concave downward tensile hoop stress distribution along the free surface. Corresponding displacement fields, normal and shear components of stress tensor and trace of the stress tensor inside the film are presented in figures 5(d)-(i). In the presence of a substrate constrained, the grain boundary region on the upper surface evolves in time as shown in figure 6. Groove shape does not seem to be affected from the substrate constrained. Hoop stress distribution around the TJ is still symmetric with respect to the GB as presented in figure 6(b), yet the magnitude of tip stress decreases as a result of substrate constraint. Since the magnitude of stress is lower, it takes longer time to reach the same depth. Figure 7(a) shows the initial configuration of the thin film with substrate and sidewall constraints exposed to uniaxial tension σ 0 applied from its left sidewall. As shown in figure 7(b) the upper film surface is displaced along the -y and the -x directions from the right to the left edge, where there is no displacement at the right edge and at the substrate as dictated by the boundary conditions. Similarly, the left sidewall is displaced along -y and the -x directions from the substrate to the left edge corner. The corresponding hoop stresses that only exist the free  surfaces are given in figure 7(c). The maps presented in figures 7(d)-(i) clearly demonstrate that the displacements and corresponding stress components are asymmetric in this case. Figure 8 shows the evolution of the grain boundary region under double constrained condition. Although the stress distribution with respect to the GB is asymmetric as shown in figure 8(b), its effect on the groove profile is negligible. The decrease in the magnitude of tip stress is now more prominent and termination of the process takes longer. The similarity in surface morphologies is also illustrated in figure 9(a), which shows the terminal profiles for each condition. Based on these data, the motion of the triple junction throughout the process is also tracked and the change in groove depth as a function of time for all three mechanical boundary conditions is plotted in figure 9(b).  Normalized groove depth versus time data is found to be a linear function of normalized time as shown in figure 9(b) with a decreasing slope from freestanding film to substrate constrained and to double constrained film. In other words, the existence of displacement constraints slow down the tip velocity.
In order to see the effect of lowering the TJ mobility, a simulation with M=0.1 was performed on the double constrained thin film at the same level of stress, X = 0.1. Since the TJ motion is slower, the width of the groove is wider, i.e. more mass is depleted from the groove shoulders by surface diffusion as shown in figure 10(a). The groove profile is still symmetric. Increasing the level of stress, i.e. the source of asymmetry, is first investigated using M=1.0 as shown in figure 10(b). Increasing the EDTI (X = 0.4) leads to shorter times for attaining the same depth, yet the groove profile is still symmetric. On the other hand lowering the TJ mobility  Previously, asymmetric formations at the grain boundaries due to anisotropic diffusivity and/or surface energy in the form of facets or slits were reported [24][25][26][27][28]. The simulation results presented here showed that the asymmetry in driving force for surface diffusion with respect to a grain boundary due to external loadings might also cause an asymmetric groove profile evolution. In this sense, intrinsic stresses, coming from grain growth or existing defects in the film, which concentrate on grain boundary regions to form high-defect areas or hot spots [29,30], may also induce asymmetric grooves. Since the material system investigated here comprised of a single grain boundary, it is hard to compare the simulation results with an experimental work. However, with the development of novel procedures to produce thin film bicrystals in micro and nanoscale dimensions [31] and advanced materialographic sample preparation methods [32] it will be possible to validate the phenomenon presented herewith.

Conclusion
A coupled diffusion/elasticity problem was solved, and the time-dependent stress field and surface morphology in bicrystaline thin films subjected to different displacement constraints were determined. Since the value of the applied stress is the same in both cases (not originated from a mismatch strain), this study provides means to reveal the effect of the displacement constraints directly. The following conclusions can be drawn from the study: • The existence of displacement constraints change the groove shape and grooving kinetics.
• In freestanding films, groove evolution is faster than substrate and/or sidewall constrained films.
• High TJ mobilities and low levels of applied stress makes it harder to observe the effects of displacement constraints to groove shape, even in the case of asymmetric stress distributions inside the film.  • Low TJ mobilities at moderate applied stresses allow formation of asymmetric groove profiles by superimposing with the asymmetric elastic stress field with respect to the GB.
These results may give clues for the strain engineering of thin films.