A sharp-interface model and its numerical approximation for solid-state dewetting with axisymmetric geometry

# A sharp-interface model and its numerical approximation for solid-state dewetting with axisymmetric geometry

## Abstract

Based on thermodynamic variation, we rigorously derive the sharp-interface model for solid-state dewetting on a flat substrate in three dimensions in the form of cylindrical symmetry. The governing equations for the model belong to fourth-order 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 sharp-interface 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 solid-state dewetting.

###### keywords:
Solid-state dewetting, axisymmetry, thermodynamic variation, surface diffusion, contact line migration, parametric finite element method (PFEM).
1\biboptions

sort&compress

## 1 Introduction

A drop of water on a leaf or a soap bubble in the air tend to form near-spherical 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 nano-scale, 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 solid-state dewetting (3). Compared to liquid dewetting, solid-state dewetting can be strongly influenced by the anisotropy of surface energy (1); (8); (7); (9); (10), and the mode of mass transport for solid-state dewetting is dominated by surface diffusion (11); (12) rather than fluid dynamics.

In recent year, the solid-state 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 solid-state 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 two-dimensional case (13) and three-dimensional 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 two-dimensional solid-state 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 fourth-order (for weak anisotropy) or sixth-order (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 re-meshing to redistribute the mesh points for the numerical stability. On the contrary, the PFEM overcomes these drawbacks by a semi-implicit 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 solid-state dewetting.

In this work, we study the three-dimensional solid-state dewetting problem by imposing the axis-symmetry 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 cross-section 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

 (s,φ)→S(s,φ):=(r(s)cosφ,r(s)sinφ,z(s)). (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 solid-state dewetting problems can be written as

 W=∬SγFV(N)dS+(γFS−γVS)A(Γo/Γi)Substrateenergy. (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

 θ=arctanzsrs,γ(θ)−γ(−θ)=0,∀θ∈[0,π],γ(θ)∈C2([0,π]). (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

 W=∬Sγ(θ)dS+(γFS−γVS)(πr2o−πr2i)Substrateenergy. (1.4)

The main objectives of this paper are as follows: (1) to derive the sharp interface-model for solid-state 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 solid-state 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 solid-state 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

The unit outer normal vector of the surface is then given by

 N=Ss×Sφ|Ss×Sφ|=(−zscosφ,−zssinφ,rs). (2.2)

We consider as the first and as the second parameter for the surface, then the first fundamental form is given by

 I=Eds2+2Fdsdφ+Gdφ2, (2.3)

with Or it can be written in the metric tensor notation as

 (gij)=(100r2),  g=r2. (2.4)

Consider a small axis-symmetric perturbation of the surface and denote the perturbed surface by :

 Sε:=(rε(s)cosφ,rε(s)sinφ,zε(s)). (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

 Γε:=Xε=X+εX(1), (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

 Wε = ∬Sεγ(θε)dSε+(γFS−γVS)(π(rεo)2−π(rεi)2)Substrateenergy, (2.7) = ∫2π0∫L0γ(θε)√gεds+(γFS−γVS)π[rε(L)2−rε(0)2], = 2π∫L0γ(θε)|Xεs|rεds+(γFS−γVS)π[rε(L)2−rε(0)2].

Denote and as the unit tangential and outer normal vector of the curve respectively, we expand the following terms at and obtain

 rε=r+r(1)ε+O(ε2), |Xεs|=1+(τ⋅X(1)s)ε+O(ε2), γ(θε)=γ(θ)+γ′(θ)(n⋅X(1)s)ε+O(ε2). (2.8)

Now in order to calculate the first variation of the total energy of the system, we can write the total energy as

 Wε=W+εW(1)+O(ε2). (2.9)

Here is the coefficient of the term, substituting Eq. (2.8) into Eq. (2.7), we get

 W(1) = 2π∫L0[γ′(θ)(n⋅X(1)s)r+γ(θ)(τ⋅X(1)s)r+γ(θ)r(1)]ds (2.10) +2π(γFS−γVS)[r(L)r(1)(L)−r(0)r(1)(0)]

Integration by parts, we obtain

 W(1) = 2π∫L0[−(γ(θ)τ+γ′(θ)n)s⋅X(1)r−(γ(θ)τ+γ′(θ)n)⋅X(1)rs (2.11) +γ(θ)r(1)]ds+2π[[γ(θ)τ+γ′(θ)n]⋅X(1)r]∣∣s=Ls=0 +2π(γFS−γVS)[r(L)r(1)(L)−r(0)r(1)(0)].

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:

 ns=κτ,n∣∣s=0=(−sinθid,cosθid),n∣∣s=L=(−sinθod,cosθod),θs=−κ, τs=−κn,τ∣∣s=0=(cosθid,sinθid),τ∣∣s=L=(cosθod,sinθod),τ⋅X(1)rs−r(1)=zsX(1)⋅n.

In addition, we requires that the two contact lines are always on the substrate after the small cylindrical perturbation, that is to say

 X(1)∣∣s=0=(r(1)(0),0),X(1)∣∣s=L=(r(1)(L),0).

Then we have

 W(1) = 2π∫L0[(γ(θ)+γ′′(θ))κ(n⋅X(1))r−(γ(θ)zs+γ′(θ)rs)(X(1)⋅n)]ds (2.12) +2πr(L)r(1)(L)[γ(θod)cosθod−γ′(θod)sinθod+(γFS−γVS)] −2πr(0)r(1)(0)[γ(θid)cosθid−γ′(θid)sinθid+(γFS−γVS)] = ∬S[(γ(θ)+γ′′(θ))κ0−γ(θ)zs+γ′(θ)rsr]n⋅X(1)dS +∫Γo[γ(θod)cosθod−γ′(θod)sinθod+(γFS−γVS)]r(1)(L)dΓ −∫Γi[γ(θid)cosθid−γ′(θid)sinθid+(γFS−γVS)]r(1)(0)dΓ

From the above equations, we immediately obtain the variation of the total energy with respect to the surface and the two contact lines as

 δWδS=(γ(θ)+γ′′(θ)κ−γ(θ)zs+γ′(θ)rsr, (2.13) δWδΓo=γ(θod)cosθod−γ′(θod)sinθod+(γFS−γVS), (2.14) δWδΓi=−[γ(θid)cosθid−γ′(θid)sinθid+(γFS−γVS)]. (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

 X(1)∣∣s=0=(0,z(1)(0)),X(1)∣∣s=L=(r(1)(L),0). (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 Gibbs-Thomson relation (23), the chemical potential can be defined as

 μ=Ω0δWδS=Ω0[(γ(θ)+γ′′(θ))κ−γ(θ)zs+γ′(θ)rsr], (2.17)

with representing the atomic volume. Subsequently the normal velocity of the surface is given by surface diffusion (12)

 vn=DsνΩ0kBTe∇2Sμ. (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 time-dependent Ginzburg-Landau kinetic equations (18); (19),

 voc=−ηδWδΓo=−η[γ(θod)cosθod−γ′(θod)sinθod+(γFS−γVS)], (2.19) vic=−ηδWδΓi=η[γ(θid)cosθid−γ′(θid)sinθid+(γFS−γVS)], (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 solid-state dewetting with anisotropic surface energy in three dimensions with cylindrical symmetry can be described in the following dimensionless form

 ∂tX=1r∂s(r∂sμ)n,00, (2.21a) μ=[γ(θ)+γ′′(θ)]κ−γ(θ)∂sz+γ′(θ)∂srr, (2.21b) κ=−∂ssX⋅n; (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

 X(s,0):=X0(s)=(r(s,0),z(s,0))=(r0(s),z0(s)),0≤s≤L0:=L(0), (2.22)

The above governing equations are together with the following boundary conditions

• contact line condition

 z(L,t)=0,{z(0,t)=0,ifr(0,t)>0,∂sz(0,t)=0,otherwise,t≥0; (2.23)
• relaxed (or dissipative) contact angle condition

 ∂tr(L,t)=−ηf(θod;σ),{∂tr(0,t)=ηf(θid;σ),ifr(0,t)>0,r(0,t)=0,otherwise,t≥0; (2.24)
• zero-mass flux condition

 ∂sμ(0,t)=0,∂sμ(L,t)=0,t≥0; (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

 f(θ;σ)=γ(θ)cosθ−γ′(θ)sinθ−σ,θ∈[−π,π],σ=γVS−γFSγ0. (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:

 V(t)=V(0),∀0≤t≤T.
###### Proof.

We know we can write as

 V(t)=∫2π0∫L(t)0rzrsdsdφ.

Take derivative with respect to for we obtain directly

 ddtV(t) = 2πddt∫10rrρzdρ=2π∫10(rtrρz+rrρzt)dρ+2π∫10rrρ,tzdρ = 2π∫10(rρzrt+rrρzt)dρ−2π∫10(rz)ρrtdρ+2π(rzrt)∣∣ρ=1ρ=0

By noting Eqs. (2.23),(2.24), the boundary terms will be cancelled, and the above equation can be reduced to

 ddtV(t) = 2π∫10(rrρzt−rzρrt)dρ=2π∫10rXt⋅(−Xρ)⊥dρ = 2π∫L(t)0rXt⋅nds=2π∫L(t)0∂s(r∂sμ)ds=0.

Here we have using the Eq. (2.21a) to replace by . Thus the mass is conserved during the time evolution. ∎

###### Proposition 2.2 (Energy dissipation).

Let and satisfies the PDEs (2.21a)-(2.21b) with boundary conditions (2.23)-(2.25), let denote the energy of the whole system, we then have:

 W(0)≤W(t1)≤W(t2),∀0≤t1≤t2.
###### Proof.

By noting the Eq. (2.12), we can easily obtain that the time derivative of the total energy as

 ddtW(t) = 2π∫L(t)0rn⋅Xt[(γ(θ)+γ′′(θ))κ−zsγ(θ)+rsγ′(θ)r]ds (2.27) +∫Γo[γ(θod)cosθod−γ′(θod)sinθod−σ)]drodtdΓ −∫Γi[γ(θid)cosθid−γ′(θid)sinθid−σ)]dridtdΓ.

Combining with Eq. (2.21a), (2.21b), we have

 ddtW(t) = Missing or unrecognized delimiter for \big = −2π∫L(t)0r(μs)2ds−2πη[ro(drodt)2+ri(dridt)2]≤0,

which immediately implies the 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 solid-state 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

 minΩ(∫2π0∫L0γ(θ)dsdϕ−σ(πr2o−πr2i)),s.t.|Ω|=const. (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

 0=δW=∬Sμn⋅X(1)dS+σ∫Γof(θod,σ)r(1)(L)dΓ−σ∫Γif(θid,σ)r(1)(0)dΓ, (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

 μ(s):=[γ(θ)+γ′′(θ)]κ−γ(θ)∂sz+γ′(θ)∂sr)r≡C,a.e.s∈[0,L], (2.30a) r(0)=0,f(θco;σ)=γ(θoc)cosθ−γ′(θoc)sinθoc−σ=0. (2.30b)

Following the Winterbottom construction (24); (20), we know the equation for the equilibrium shape for solid-state dewetting under the cylindrical symmetry can be given as follows:

 {r(θ)=γ(θ)sinθ+γ′(θ)cosθ,z(θ)=γ(θ)cosθ−γ′(θ)sinθ−σ,θ∈[0,θoc]. (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

 L2(I)={f:I→R,and∫Γ(t)|f(s)|2ds=∫I|f(s(ρ,t))|2∂ρsdρ<+∞},

equipped with the inner product for any scalar (or vector) functions . Moreover, we define the the following two special functional spaces:

 H(r)a,b(I)={ϕ∈H1(I): ϕ(0)=a;ϕ(1)=b}, (3.1) H(z)a,b(I)={ϕ∈H1(I): ϕ(1)=0;ifa>0, ϕ(0)=0}, (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 solid-state 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:

 ⟨∂tX⋅n, ϕr⟩Γ+⟨∂sϕ, ∂sμr⟩Γ=0,∀ϕ∈H1(I), (3.3a) ⟨μ, φr⟩Γ−⟨˜γ(θ)κ, φr⟩Γ+⟨b(θ)⋅τ, φ⟩Γ=0,∀φ∈H1(I), (3.3b) Missing or unrecognized delimiter for \big (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 re-formulating 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

 Vh:={u∈C(I):u∣Ij∈P1,j=1,2,…,N}⊂H1(I), (3.4) Vh,(r)a,b:={u∈Vh:u(0)=a; u(1)=b}⊂H(r)a,b(I), (3.5) Vh,(z)a,b:={u∈Vh: if a>0,u(0)=0; u(1)=0}⊂H(z)a,b(I), (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

 Γm=N⋃j=1¯hmj,{hmj}Nj=1are line segments of the curveΓm, (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

 ∂sϕ=∂ρϕ∂ρs=∂ρϕ|∂ρXm|. (3.8)

Thus, the tangential vector is a step function with discontinuities on the endpoints of each line segment . It can be calculated as

 τm∣∣hmj=∂sXm∣∣hmj=Xm(ρj)−Xm(ρj−1)|Xm(ρj)−Xm(ρj−1)|,∀1≤j≤N. (3.9)

Consequently, the normal vector and angle theta can be obtained directly

 nm∣∣hmj=−τm⊥∣∣hmj,θm∣∣hmj=arctanzm(ρj)−zm(ρj−1)rm(ρj)−rm(ρj−1). (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

 ⟨f, g⟩Γm=16N∑j=1|hmj|[f((ρj−1)+)⋅g((ρj−1)+)+4f(ρj−12)⋅g(ρj−12)+f((ρj)−)⋅g((ρj)−)]. (3.11)

Here and are one-sided limit and can be defined respectively as

 f(ρ±)=limε→0+f(ρ±ε). (3.12)

Now the fully-discrete parametric finite element approximation can be described as: Take