A sharpinterface model and its numerical approximation for solidstate dewetting with axisymmetric geometry
Abstract
Based on thermodynamic variation, we rigorously derive the sharpinterface model for solidstate dewetting on a flat substrate in three dimensions in the form of cylindrical symmetry. The governing equations for the model belong to fourthorder geometric partial differential equations, with proper boundary conditions such that the total volume of the system is conserved and the total energy is dissipative during the time evolution. We propose a weak formulation for the sharpinterface model and then apply parametric finite element method to solve it efficiently. Extensive numerical simulations results are presented lastly to demonstrate the morphological characteristics for solidstate dewetting.
keywords:
Solidstate dewetting, axisymmetry, thermodynamic variation, surface diffusion, contact line migration, parametric finite element method (PFEM).sort&compress
1 Introduction
A drop of water on a leaf or a soap bubble in the air tend to form nearspherical shape, as the sphere has the smallest surface area to volume ratio. This phenomenon, known as surface tension, results form the cohesive forces among the liquid molecules, and drives the system to move towards state with lower surface energy. In micro and nanoscale, surface tension effects can also be observed in solids. The solid thin films supported by rigid substrates are typically unstable. They could dewet or agglomerate to form small islands on the substrate, and exhibit complex morphological evolutions even below their melting temperature (1); (2); (3); (4); (5); (6); (7). During this process, the thin film remains in the solid state, thus it is called solidstate dewetting (3). Compared to liquid dewetting, solidstate dewetting can be strongly influenced by the anisotropy of surface energy (1); (8); (7); (9); (10), and the mode of mass transport for solidstate dewetting is dominated by surface diffusion (11); (12) rather than fluid dynamics.
In recent year, the solidstate dewetting problems have been widely studied and different mathematical models have been proposed to investigate its kinetics. From these models, surface diffusion flow and contact line migration have been recognized as the two main kinetics. The first sharp interface model for solidstate dewetting was proposed by D.J. Srolovitz and S.A. Safran (11) to investigate the hole growth under the assumption of isotropic surface energy and cylindrical symmetry. This model was then generalized to both the twodimensional case (13) and threedimensional case (14) in Lagrangian representation. Moreover, the the ”marker particle” method was simultaneously proposed for solving these models. Some other models have been proposed to include the anisotropy, such as the discrete model by Dornel (15), kinetic Monte Carlo models (16); (17) and model via crystalline method (8); (9). More recently, the twodimensional solidstate dewetting has been fully studied via sharp interface models (18); (19); (20); (21). In these models, the interface between the thin film and vapor is assumed to be an open curve with two end points(the contact points) attached on the axis(flat substrate). The open curve evolves via surface diffusion flow while the positions of the two contact points are determined by the relaxed contact angle conditions. Compared to previous models, these model are rigorously derived via the energy variational method, thus account for the full anisotropic free energy of the system and represent a completely mathematical description. The governing equations for the model belongs to fourthorder (for weak anisotropy) or sixthorder (for strong anisotropy) geometric partial differential equations (PDEs) with boundary conditions at the two contact points. These geometric PDEs were then numerically solved by the novel parametric finite element method (PFEM) (21); (22), which show great advantages over the previous ”marker particle” method. Precisely, ”the marker particle” method is an explicit finite difference scheme, thus posing a very severe restriction on time step for the numerical stability. Besides, at each time step, this algorithm requires remeshing to redistribute the mesh points for the numerical stability. On the contrary, the PFEM overcomes these drawbacks by a semiimplicit mixed form finite element method. By introducing a particular tangential motion for the mesh points, the PFEM has very good properties with respect to the distribution of the mesh points. Although extensive numerical examples have been presented in these papers (18); (19); (20); (21), and a generalized Winterbottom construction was proposed to study the equilibrium problem, yet these numerical simulations are restricted to two dimensions and can not fully demonstrate the complicated geometric evolution during the process of solidstate dewetting.
In this work, we study the threedimensional solidstate dewetting problem by imposing the axissymmetry to both geometry of the thin film and the surface energy density. The assumption of this symmetry could greatly reduce the surface evolution to curve evolution by merely considering the crosssection profile of the thin film. As illustrated in Fig. 1, a cylindrical thin film lies on a fat substrates. The interface that separates the vapor and thin film is represented by a open surface with boundaries given by two closed curved and .
Due to the cylindrical symmetry, the open surface can be parameterised as follows
(1.1) 
Here in these expressions, is the radial distance, is the azimuth angle, is the film height, and represents the arc length along the axial direction curve (the generatrix, see Fig. 1(b)). and represent the inner contact line denoted by and outer contact line denoted by , respectively.
The total free energy of the system for solidstate dewetting problems can be written as
(1.2) 
Here denotes the area enclosed by the two contact lines, and are constants and represent the surface energy densities of film/substrate and vapor/substrate respectively. The surface energy density of the film/vapor(surface ) is given by , which is dependent on the unit out normal vector of the surface. The cylindrical symmetry can reduce this dependence to the orientation of curve in the axial direction. Thus we use to denote the surface energy density of film/vapor, with
(1.3) 
If we let and represent the radius of outer contact line and inner contact lines respectively on the substrate, then the total energy can be rewritten as
(1.4) 
The main objectives of this paper are as follows: (1) to derive the sharp interfacemodel for solidstate dewetting under the assumption of cylindrical symmetry, (2) to develop a parametric finite element method (PFEM) for solving the derived model based on a good weak formulation, (3) to demonstrate some geometric features during the evolution for solidstate dewetting. So the rest of the paper is organized as follows. In section , we show rigorously the derivation of sharp interface model via thermodynamic variation and present the equilibrium equation. In section , we present parametric finite element method for solving the geometric PDEs. In section , extensive numerical results will be presented to show the morphological characteristics for solidstate dewetting.
2 The sharp interface model
2.1 Thermodynamic variation
Given the parametrization of the surface in Eq. (1.1), the two tangential vectors of the surface can be calculated directly as follows
(2.1) 
The unit outer normal vector of the surface is then given by
(2.2) 
We consider as the first and as the second parameter for the surface, then the first fundamental form is given by
(2.3) 
with Or it can be written in the metric tensor notation as
(2.4) 
Consider a small axissymmetric perturbation of the surface and denote the perturbed surface by :
(2.5) 
Due to the symmetry, this perturbation can also be viewed as a corresponding perturbation of the curve in radial direction. As illustrated in Fig. 2, it depicts the perturbation of an initial toroidal thin film (Fig. 2(a)) or island film (Fig. 2(b)). we can rewrite the perturbation in coordinates as follows
(2.6) 
where is a small perturbation parameter and . Note here if we still take and as the two parameters of the surface , the corresponding matrix tensor .
Based on Eq. (1.4), the total free energy of the new system with perturbed surface can be written
(2.7)  
Denote and as the unit tangential and outer normal vector of the curve respectively, we expand the following terms at and obtain
(2.8) 
Now in order to calculate the first variation of the total energy of the system, we can write the total energy as
(2.9) 
Here is the coefficient of the term, substituting Eq. (2.8) into Eq. (2.7), we get
(2.10)  
Integration by parts, we obtain
(2.11)  
Let be the curvature of curve . Moreover, for simplicity, we assume that the curve has two contact points (see Fig. 2(a)). Denote as the inner and outer contact angle of curve respectively, we have the following expressions:
In addition, we requires that the two contact lines are always on the substrate after the small cylindrical perturbation, that is to say
Then we have
(2.12)  
From the above equations, we immediately obtain the variation of the total energy with respect to the surface and the two contact lines as
(2.13)  
(2.14)  
(2.15) 
Note here the variation is calculated for the case when the thin film has two contact lines with toroidal geometry. For island film, as depicted in Fig. 2(b), there is no inner contact line, thus the after the perturbation, we require at the boundary the following equations should be satisfied
(2.16) 
Thus integration by parts in Eq. (2.11) will not produce boundary terms at .
2.2 The model and its properities
From the anisotropic GibbsThomson relation (23), the chemical potential can be defined as
(2.17) 
with representing the atomic volume. Subsequently the normal velocity of the surface is given by surface diffusion (12)
(2.18) 
In these expressions, is the surface diffusivity, is the thermal energy, is the number of diffusing atoms per unit area, is the surface gradient. The two contacts lines moves on the substrates, and the velocities are given by the timedependent GinzburgLandau kinetic equations (18); (19),
(2.19)  
(2.20) 
with denoting the contact line mobility.
Now take the characteristic length scale and characteristic surface energy scale as and respectively, by choosing the time scale as with , and the contact line mobility is scaled by . Note that , thus if we let represent the curve in the radial direction, the sharp interface model for solidstate dewetting with anisotropic surface energy in three dimensions with cylindrical symmetry can be described in the following dimensionless form
(2.21a)  
(2.21b)  
(2.21c) 
where represents the total arc length of the curve , all the variables are in dimensionless form and we still use the same notation for brevity. The initial condition is given as
(2.22) 
The above governing equations are together with the following boundary conditions

contact line condition
(2.23) 
relaxed (or dissipative) contact angle condition
(2.24) 
zeromass flux condition
(2.25)
where are the (dynamic) contact angles at the inner contact line and outer contact line, respectively, denotes the dimensionless contact line mobility, and is defined as
(2.26) 
In the case when the thin film does not have an inner contact line, is fixed at and the symmetry also require that the Neumann boundary condition for z at , see Eq. (2.23).
In the following, we will show that the total volume of the thin film (the volume enclosed by the surface and the substrate) is conserved and the total energy of the system is dissipative. We introduce a new parameter to parameterize the evolution curves such that is a family of open curves in the plane, where , and . It should be noted that the relationship between the parameter and the arc length can be given as , and then we can obtain that .
Proposition 2.1 (Mass conservation).
Let and satisfies the PDEs (2.21a)(2.21b) with boundary conditions (2.23)(2.25), let denote the volume of the thin film, we then have:
Proof.
Proposition 2.2 (Energy dissipation).
In the proof of the energy dissipation, we only show the case when there exists the inner contact line. Note that when there is no inner contact line, the energy dissipation is still satisfied, as the boundary term at will no longer exist.
2.3 Equilibrium shape
The equilibrium shape for solidstate dewetting problem is a minimization of the total surface energy under the constraint such that the total volume of the thin film is fixed. The problem can be mathematically expressed in the following
(2.28) 
Based on the first variation of the dimensionless total surface energy, we know that the equation for the equilibrium shape can be obtained by vanishing the first variation of the total surface energy
(2.29) 
Assume the curve of the equilibrium shape in the axis direction is given by , then in order to ensure the first variation of the surface energy to be zero, it is equivalent to the following two conditions
(2.30a)  
(2.30b) 
Following the Winterbottom construction (24); (20), we know the equation for the equilibrium shape for solidstate dewetting under the cylindrical symmetry can be given as follows:
(2.31) 
When , in view of Eq. (1.3), we obtain , this implies that no inner contact line exists in the equilibrium shape. is the static contact angle for the outer contact line, and it satisfies , which immediately gives .
3 Parametric finite element method
3.1 Variational formulation
To formulate the variational formulation for the geometric PDEs: Eqs. (2.21a),(2.21b),(2.21c) with boundary conditions Eqs. (2.23),(2.24),(2.25), we need to introduce some functional spaces with respect to the evolution curve . First, we can define the space on the domain as
equipped with the inner product for any scalar (or vector) functions . Moreover, we define the the following two special functional spaces:
(3.1)  
(3.2) 
with , are given as two parameters related to the radius of the two movin contact lines. It is obvious that from the above definition, we have .
Let , , then the variational formulation for solidstate dewetting in three dimensions with cylindrical symmetry can be stated as follows: given , for any time , we want to find with coordinates positions of two contact lines and , such that:
(3.3a)  
(3.3b)  
(3.3c) 
Eq. (3.3b) is obtained directly by multiplying for Eq. (2.21b) and then integrating over . For Eq. (2.21a), firstly we multiply to cancel the fraction term, and then multiply the test function . By noting the boundary condition Eq. (2.25), integrating by parts over immediately gives Eq. (3.3a). Moreover, by reformulating Eq. (2.21c) as , and then multiplying the vector valued test function , integrating by parts, Eq. (3.3c) can be derived in view of the boundary condition Eq. (2.23). For the cases or , the variational formulation could be different as the functional space for the test function depends on .
3.2 Fully discrete scheme
To present parametric finite element method(PFEM) for the variational problem (3.3a),(3.3b),(3.3c), we first discrete the time as with the time step . Also a uniform partition of is given as with . Introduce the finite dimensional approximations to , and with and two given constants as
(3.4)  
(3.5)  
(3.6) 
where denotes all polynomials with degrees at most .
Let , , , and be the numerical approximations of the moving curve , the normal vector , the chemical potential , the angle and the curvature at time , respectively. The polygonal curve can be written as
(3.7) 
with these line segments are order in the clockwise direction. By noting the relationship between the parameter and , we can obtain the the following equation
(3.8) 
Thus, the tangential vector is a step function with discontinuities on the endpoints of each line segment . It can be calculated as
(3.9) 
Consequently, the normal vector and angle theta can be obtained directly
(3.10) 
where means a degree’s clockwise rotation of a vector. Moreover, given two scalar or vector functions , which is defined on , define the following discrete norms to approximate the integration over
(3.11) 
Here and are onesided limit and can be defined respectively as
(3.12) 
Now the fullydiscrete parametric finite element approximation can be described as: Take