1 Introduction↩︎

Surface diffusion is a widespread process involving the movement of adatoms, molecules and atomic clusters at solid material interfaces [1]. Due to different surface lattice orientation, an anisotropic evolution process is generated for a solid material, which is called anisotropic surface diffusion in the literature. Surface diffusion with an anisotropic surface energy plays an important role as a crucial mechanism and/or kinetics in various fields such as epitaxial growth [2], [3], surface phase formation [4], heterogeneous catalysis [5], and other pertinent fields within surface and materials science [6]–[8]. In fact, broader and consequential applications of surface diffusion have been discovered in materials science and solid-state physics, notably in areas such as the crystal growth of nanomaterials [9], [10] and solid-state dewetting [4], [11]–[17].

As shown in Fig. 1, let \(\Gamma\mathrel{\vcenter{:}}=\Gamma(t)\) be a closed curve in two dimensions (2D) associated with a given anisotropic surface energy \(\hat{\gamma}(\theta)>0\), where \(\theta\in2\pi\mathbb{T}\mathrel{\vcenter{:}}=\mathbb{R}/2\pi\mathbb{Z}\) represents the angle between the vertical axis and unit outward normal vector \(\boldsymbol{n}=\boldsymbol{n}(\theta)\mathrel{\vcenter{:}}=(-\sin\theta,\cos\theta)^T\). It should be noted that the anisotropy can also be viewed as a function \(\gamma(\boldsymbol{n})\) of the normal vector \(\boldsymbol{n}\) [13], [18], [19]. While \(\gamma(\boldsymbol{n}) = \gamma(-\sin\theta, \cos\theta)\) is equivalent to \(\hat{\gamma}(\theta)\) by the one-to-one correspondence \(\boldsymbol{n}(\theta)=(-\sin\theta, \cos\theta)^T\), the \(\hat{\gamma}(\theta)\) formulation is often more convenient and straightforward in 2D.

Suppose \(\Gamma\) is represented by \(\boldsymbol{X}\mathrel{\vcenter{:}}=\boldsymbol{X}(s,t)=(x(s,t),y(s,t))^T\), where \(s\) denotes the arc-length parameter, and \(t\) represents the time. The motion of \(\Gamma\) under anisotropic surface diffusion is governed by the following geometric flow [20], [21]: \[\label{eqn:surface32diffusion} \partial_t\boldsymbol{X}=\partial_{ss}\mu\boldsymbol{n},\tag{1}\] where \(\mu\) is the weighted curvature (or chemical potential) defined as \[\label{eqn:weighted32curvature} \mu\mathrel{\vcenter{:}}=[\hat{\gamma}(\theta)+\hat{\gamma}^{\prime\prime}(\theta)]\kappa\tag{2}\] with \(\kappa\mathrel{\vcenter{:}}=-(\partial_{ss}\boldsymbol{X})\cdot\boldsymbol{n}\) being the curvature.

The anisotropic surface diffusion 1 is a fourth-order and highly nonlinear geometric flow, which possesses two major geometric properties, i.e., the area conservation and the energy dissipation. Let \(A_c(t)\) be the area of the region \(\Omega(t)\) enclosed by \(\Gamma(t)\), and \(W_c(t)\) be the total surface free energy, which are defined as \[A_c(t)\mathrel{\vcenter{:}}=\int_{\Omega(t)}1\,\mathrm{d}x\mathrm{d}y=\int_{\Gamma(t)}y(s,t)\partial_sx(s,t)\,\mathrm{d}s,\qquad W_c(t)\mathrel{\vcenter{:}}=\int_{\Gamma(t)}\hat{\gamma}(\theta)\,\mathrm{d}s.\] One can prove that [7], [16], [22] \[\frac{\mathrm{d}}{\mathrm{d}t}A_c(t)=0,\qquad\frac{\mathrm{d}}{\mathrm{d}t}W_c(t)=-\int_{\Gamma(t)}|\partial_s\mu|^2\,\mathrm{d}s\leq 0,\] which immediately implies the anisotropic surface diffusion 1 –2 satisfies the area conservation and energy dissipation, i.e., \[A_c(t)\equiv A_c(0),\qquad W_c(t)\leq W_c(t_1)\leq W_c(0),\qquad\forall t\geq t_1\geq0.\]

When \(\hat{\gamma}(\theta)\equiv 1,\forall\theta\in 2\pi\mathbb{T}\), the weighted curvature \(\mu\) reduces to \(\kappa\), and it is referred to as isotropic surface energy. If \(\hat{\gamma}(\theta)\) is not a constant function and \(\hat{\gamma}(\theta)+\hat{\gamma}^{\prime\prime}(\theta)>0\) for all \(\theta\in2\pi\mathbb{T}\), we classify the surface energy as weakly anisotropic, otherwise, it is termed strongly anisotropic. Typical anisotropic surface energies \(\hat{\gamma}(\theta)\) which are widely employed in materials science include:

1. the \(m\)-fold anisotropic surface energy [22] \[\label{eqn:k-fold32anisotropy} \hat{\gamma}(\theta)=1+\beta\cos m(\theta-\theta_0),\qquad\theta\in 2\pi\mathbb{T},\tag{3}\] where \(m=2,3,4,6,|\beta|<1\) are dimensionless strength constants, \(\theta_0\in 2\pi\mathbb{T}\) is a constant. Note that \(\hat{\gamma}(\theta)\) is weakly anisotropic when \(|\beta|<\frac{1}{m^2-1}\), and strongly anisotropic otherwise.

2. the ellipsoidal anisotropic surface energy [23] \[\label{eqn:ellipsoidal32anisotropy} \hat{\gamma}(\theta)=\sqrt{a+b\cos^2\theta},\qquad\theta\in2\pi\mathbb{T},\tag{4}\] where \(a,b\) are two dimensionless constants satisfying \(a>0\) and \(a+b>0\).

3. the Riemannian-like metric (also called BGN) anisotropic surface energy [24], [25] \[\label{eqn:BNG32anisotropy} \hat{\gamma}(\theta)=\sum_{l=1}^L\sqrt{\boldsymbol{n}(\theta)^TG_l\boldsymbol{n}(\theta)},\qquad\boldsymbol{n}=(-\sin\theta,\cos\theta)^T,\qquad\theta\in 2\pi\mathbb{T},\tag{5}\] where \(L\geq 1\), \(G_l\in\mathbb{R}^{2\times 2}\) positive definite \(\forall 1\leq l\leq L\). When \(L=1,G_1=\text{diag}(a,a+b)\) in 5 , the Riemannian-like metric anisotropy 5 collapses to the ellipsoidal anisotropic surface energy 4 .

4. the piecewise anisotropic surface energy [26] \[\label{eqn:piecewise32anisotropy} \hat{\gamma}(\theta)=\sqrt{\left(\frac{5}{2}+\frac{3}{2}\text{sgn}(n_1)\right)n_1^2+n_2^2},\tag{6}\] with \(\boldsymbol{n}=(n_1,n_2)^T\mathrel{\vcenter{:}}=(-\sin\theta,\cos\theta)^T,\,\,\theta\in 2\pi\mathbb{T}.\)

Many numerical methods have been proposed for simulating the evolution of curves under surface diffusion, including the phase-field method [15], [27]–[29], the discontinuous Galerkin method [30], the marker particle method [31], [32] and the parametric finite element method (PFEM) [33]–[36]. Among these methods, the energy-stable PFEM (ES-PFEM) by Barrett, Garcke, and Nürnberg [33], also referred to as BGN’s method, achieves the best performance in terms of mesh quality and unconditional energy-stability in the isotropic case. The ES-PFEM has been further extended to other geometric flows, such as the solid-state deweting problem [12], [22], [37], demonstrating its adaptability and robustness. Furthermore, Bao and Zhao have recently developed a structure-preserving PFEM (SP-PFEM) [38]–[40], which can preserve the enclosed mass at the fully-discretized level while also maintaining the unconditional energy stability. Extending these PFEMs to anisotropic surface diffusion is a major focus of recent research in surface diffusion. While BGN successfully adapted their methods to a specific Riemannian metric form [24], [25], designing a SP-PFEM for anisotropic surface diffusion with arbitrary anisotropies remains challenging.

Lately, based on the \(\hat{\gamma}(\theta)\) formulation, Li and Bao introduced a surface energy matrix \(G(\theta)\) and extended the ES-PFEM from the isotropic cases to the anisotropic cases [41]. Due to the absence of a stabilizing term, their method requires a certain restrictive condition to ensure the energy stability. Subsequently, Bao, Jiang, and Li incorporated a stabilizing function within the \(\gamma(\boldsymbol{n})\) formulation. They constructed a symmetric surface energy matrix \(\boldsymbol{Z}_k(\boldsymbol{n})\) and proposed a symmetrized SP-PFEM for the anisotropic surface diffusion in [42], [43]. The symmetrized SP-PFEM with the stabilizing function works effectively for symmetric surface energy distributions (i.e., \(\gamma(\boldsymbol{n})=\gamma(-\boldsymbol{n})\)) to maintain the geometric properties. However, there are different anisotropic surface energies \(\hat{\gamma}(\theta)\) which are not symmetrically distributed or do not satisfy the specific condition, such as the \(3\)-fold anisotropic surface energy 3. Very recently, based on the \(\gamma(\boldsymbol{n})\) formulation, Bao and Li introduced a novel surface energy matrix and established a comprehensive analytical framework to demonstrate the unconditional energy stability of the proposed SP-PFEM [44], [45]. Based on this framework, they successfully reduced the requirement for the anisotropy to \(3\gamma(\boldsymbol{n})>\gamma(-\boldsymbol{n}),\,\,\forall\boldsymbol{n}\in\mathbb{S}^1\).

However, the critical situation \(3\gamma(\boldsymbol{n}^*)=\gamma(-\boldsymbol{n}^*),\,\,\boldsymbol{n}^*\in\mathbb{S}^1\) was not addressed in their analytical framework. This is because their framework relies on an estimate of the gradient of \(\gamma(\boldsymbol{n})\). The challenge comes from \(\gamma(\boldsymbol{n})\) being a function defined on the unit sphere \(\mathbb{S}^1\), which makes its gradient complicated to analyze. In contrast, the \(\hat{\gamma}(\theta)\) formulation, defined on \(2\pi\mathbb{T}\), possesses a simpler derivative and thus allows for a better analysis of the critical situation \(3\hat{\gamma}(\theta^*)=\hat{\gamma}(\theta^*-\pi),\,\,\theta^*\in2\pi\mathbb{T}\). Therefore, inspired by their analytical framework, we adopt the \(\hat{\gamma}(\theta)\) formulation to further investigate the critical situation.

The main objective of this paper is to propose a structure-preserving stabilized parametric finite element method (SPFEM) for simulating surface diffusion 1 with the surface energy \(\hat{\gamma}(\theta)\) under very mild conditions as

(i) \(3\hat{\gamma}(\theta)\geq\hat{\gamma}(\theta-\pi),\qquad\forall\theta\in2\pi\mathbb{T},\)

(ii) \(\hat{\gamma}^\prime(\theta^*)=0\), when \(3\hat{\gamma}(\theta^*)=\hat{\gamma}(\theta^*-\pi),\qquad\theta^*\in 2\pi\mathbb{T}\).

Compared to the \(\gamma(\boldsymbol{n})\) formulation, the \(\hat{\gamma}(\theta)\) formulation has the following advantages:

(i) it is more intuitive and has a simpler form in practical applications;

(ii) it enables a reduction in the regularity requirement for the anisotropy from \(C^2\) to globally \(C^1\) and piecewise \(C^2\);

(iii) it allows for a more convenient discussion of critical situations as \(3\hat{\gamma}(\theta^*)=\hat{\gamma}(\theta^*-\pi)\) for some \(\theta^*\in 2\pi\mathbb{T}\) in 2D.

The remainder of this paper is organized as follows: In section 2, we introduce a stabilized surface energy matrix \(\hat{\boldsymbol{G}}_k(\theta)\) and propose a new conservative formulation for anisotropic surface diffusion. In section 3, we present a novel weak formulation based on the conservative form, introduce its spatial semi-discretization, and propose a full discretization by SPFEM. In section 4, we analyze the structure-preserving properties of the proposed scheme, i.e., area conservation and unconditional energy stability, and develop a comprehensive framework for energy stability. It starts from defining a minimal stabilizing function \(k_0(\theta)\), then we obtain the main result through a local energy estimate under the assumption that \(k_0(\theta)\) is well-defined. The existence of \(k_0(\theta)\) is further demonstrated in section 5. In section 6, we extend the SPFEM for simulating solid-state dewetting of thin films under anisotropic surface diffusion and contact line migration. Extensive numerical results are provided in section 7 to demonstrate the efficiency, accuracy and structure-preserving properties of the proposed SPFEM. Finally, we conclude in section 8.

2 A conservative formulation↩︎

In this section, we propose a novel formulation with stabilization for 1 and derive a conservative form by introducing a stabilized surface energy matrix.

Applying the geometric identity \(\kappa\boldsymbol{n}=-\partial_{ss}\boldsymbol{X}\) [46], the anisotropic surface diffusion equations 1 –2 can be reformulated into \[\tag{7} \begin{align} &\boldsymbol{n}\cdot\partial_t\boldsymbol{X}-\partial_{ss}\mu=0,\qquad 0<s<L(t),\qquad t>0,\tag{8}\\ &\mu\boldsymbol{n}+\left[\hat{\gamma}(\theta)+\hat{\gamma}^{\prime\prime}(\theta)\right]\partial_{ss}\boldsymbol{X}=\boldsymbol{0}.\tag{9} \end{align}\] where \(L(t)=\int_{\Gamma(t)}\,\mathrm{d}s\) is the length of \(\Gamma(t)\).

For a vector \(\boldsymbol{v}=(v_1,v_2)^T\in\mathbb{R}^2\), we denote \(\boldsymbol{v}^\perp\in\mathbb{R}^2\) as its perpendicular vector which is the clockwise rotation of \(\boldsymbol{v}\) by \(\frac{\pi}{2}\), i.e. \[\boldsymbol{v}^\perp\mathrel{\vcenter{:}}=(v_2,-v_1)^T=-J\boldsymbol{v},\qquad \text{with}\quad J=\left(\begin{array}{cc} 0 & -1 \\ 1 & 0 \end{array}\right).\] Then the tangent vector \(\boldsymbol{\tau}\mathrel{\vcenter{:}}=\partial_s\boldsymbol{X}\), and unit normal vector \(\boldsymbol{n}\) can be written as \(\boldsymbol{n}=-\boldsymbol{\tau}^\perp\). And the tangent vector can also be given by \(\boldsymbol{\tau}=\boldsymbol{n}^{\perp}\).

Applying ?? , the governing geometric PDE 7 for anisotropic surface diffusion can be reformulated as the following conservative form \[\tag{17} \begin{align} &\boldsymbol{n}\cdot\partial_t\boldsymbol{X}-\partial_{ss}\mu=0,\tag{18}\\ &\mu\boldsymbol{n}+\partial_s\left(\hat{\boldsymbol{G}}_k(\theta)\partial_s\boldsymbol{X}\right)=0.\tag{19} \end{align}\]

3 A SPFEM for anisotropic surface diffusion↩︎

In this section, we first develop a novel weak formulation based on the conservative form 17 and present the spatial semi-discretization of this weak formulation. After that, a structure-preserving SPFEM is proposed by adapting the implicit-explicit Euler method in time, which preserves area conservation and energy dissipation at the discrete level.

3.1 Weak formulation↩︎

In order to derive a weak formulation of equation 17 , we introduce a time-independent variable \(\rho\) which parameterizes \(\Gamma(t)\) over a fixed domain \(\rho\in\mathbb{I}=[0,1]\) as \[\Gamma(t)\mathrel{\vcenter{:}}=\boldsymbol{X}(\rho,t)=(x(\rho,t),y(\rho,t))^T\colon\mathbb{I}\times\mathbb{R}^+\to\mathbb{R}^2.\] The arc-length parameter \(s\) can thus be computed by \(s=\int_0^\rho|\partial_\rho\boldsymbol{X}(q,t)|\,\mathrm{d}q\). (We do not distinguish \(\boldsymbol{X}(\rho,t)\) and \(\boldsymbol{X}(s,t)\) if there’s no misunderstanding.)

Introduce the following functional space with respect to the evolution curve \(\Gamma(t)\) as \[L^2(\mathbb{I})\mathrel{\vcenter{:}}=\left\{u\colon\mathbb{I}\to\mathbb{R}\mid\int_{\Gamma(t)}|u(s)|^2\,\mathrm{d}s=\int_{\mathbb{I}}|u(s(\rho,t))|^2\partial_\rho s\,\mathrm{d}s<+\infty\right\},\] equipped with the \(L^2\)-inner product \[\Bigl(u,v\Bigr)_{\Gamma(t)}\mathrel{\vcenter{:}}=\int_{\Gamma(t)}u(s)v(s)\,\mathrm{d}s=\int_{\Gamma(t)}u(s(\rho,t))v(s(\rho,t))\partial_\rho s\,\mathrm{d}\rho,\] for any scalar (or vector) functions. The Sobolev spaces are defined as \[\begin{align} H^1(\mathbb{I})&\mathrel{\vcenter{:}}=\left\{u\colon\mathbb{I}\to\mathbb{R}\mid u\in L^2(\mathbb{I}),\,\,\text{and}\,\,\partial_\rho u\in L^2(\mathbb{I})\right\},\\ H^1_p(\mathbb{I})&\mathrel{\vcenter{:}}=\left\{u\in H^1(\mathbb{I})\mid u(0)=u(1)\right\}. \end{align}\] Extensions of above definitions to the functions in \([L^2(\mathbb{I})]^2,[H^1(\mathbb{I})]^2\) and \([H^1_p(\mathbb{I})]^2\) are straightforward.

By multiplying the equation 18 by a test function \(\varphi\in H^1_p(\mathbb{I})\), integrating over \(\Gamma(t)\), and applying integration by parts, we obtain \[\label{variation14432sd} \Bigl(\boldsymbol{n}\cdot \partial_t\boldsymbol{X}, \varphi\Bigr)_{\Gamma(t)}+\Bigl(\partial_s\mu, \partial_s\varphi\Bigr)_{\Gamma(t)}=0.\tag{20}\]

Similarly, by taking the dot product of equation 19 with a test function \(\boldsymbol{\omega}=(\omega_1,\omega_2)^T\in[H^1_p(\mathbb{I})]^2\) and integrating by parts, we have \[\label{variation24432sd} \begin{align} 0&=\Bigl(\mu\boldsymbol{n}+\partial_s\left(\hat{\boldsymbol{G}}_k(\theta)\partial_s\boldsymbol{X}\right),\boldsymbol{\omega}\Bigr)_{\Gamma(t)}\\ &=\Bigl(\mu\boldsymbol{n},\boldsymbol{\omega}\Bigr)_{\Gamma(t)}+\Bigl(\partial_s(\hat{\boldsymbol{G}}_k(\theta)\partial_s\boldsymbol{X}),\boldsymbol{\omega}\Bigr)_{\Gamma(t)}\\ &=\Bigl(\mu\boldsymbol{n},\boldsymbol{\omega}\Bigr)_{\Gamma(t)}-\Bigl(\hat{\boldsymbol{G}}_k(\theta)\partial_s\boldsymbol{X},\partial_s\boldsymbol{\omega}\Bigr)_{\Gamma(t)} \end{align}\tag{21}\]

Combining 20 and 21 , we propose a new weak formulation for 17 as follows: Given an initial closed curve \(\Gamma(0)\mathrel{\vcenter{:}}=\boldsymbol{X}(\cdot,0)=\boldsymbol{X}_0\in[H^1_p(\mathbb{I})]^2\), find the solution \((\boldsymbol{X}(\cdot,t)=(x(\cdot,t),y(\cdot,t))^T,\mu(\cdot,t))\in[H^1_p(\mathbb{I})]^2\times H^1_p(\mathbb{I})\), such that: \[\tag{22} \begin{align} \begin{aligned} &\Bigl(\boldsymbol{n}\cdot \partial_t\boldsymbol{X}, \varphi\Bigr)_{\Gamma(t)}+\Bigl(\partial_s\mu, \partial_s\varphi\Bigr)_{\Gamma(t)}=0,\qquad \forall\varphi\in H^1_p(\mathbb{I}),\tag{23} \end{aligned}\\ \begin{align} &\Bigl(\mu\boldsymbol{n},\boldsymbol{\omega}\Bigr)_{\Gamma(t)}-\Bigl(\hat{\boldsymbol{G}}_k(\theta)\partial_s\boldsymbol{X},\partial_s\boldsymbol{\omega}\Bigr)_{\Gamma(t)}=0,\qquad\forall\boldsymbol{\omega}\in[H^1_p(\mathbb{I})]^2.\tag{24} \end{align} \end{align}\]

It can be demonstrated that the weak formulation 22 maintains two geometric properties, namely, area conservation and energy dissipation.

To prove the above theorem, we first introduce the following transport lemma:

Now the proof of Proposition 1 is ready to be presented:

3.2 A semi-discretization in space↩︎

To obtain the spatial discretization, let \(N>2\) be a positive integer and \(h=1/N\) be the mesh size, grid points \(\rho_j=jh\), sub-intervals \(I_j=[\rho_{j-1},\rho_j]\) for \(j=1,2,\dots,N\) and the uniform partition \(\mathbb{I}=[0,1]=\cup_{j=1}^NI_j\). The closed curve \(\Gamma(t)=\boldsymbol{X}(\cdot,t)\) is approximated by the polygonal curve \(\Gamma^h(t)=\boldsymbol{X}^h(\cdot,t)=\left(x^h(\cdot,t),y^h(\cdot,t)\right)^T\) satisfying \(\boldsymbol{X}^h(\rho_j,0)=\boldsymbol{X}(\rho_j,0)\).

The polygon \(\Gamma^h(t)\) is composed of ordered line segments \(\{\boldsymbol{h}_j(t)\}_{j=1}^N\), i.e. \[\Gamma^h(t)=\bigcup_{j=1}^N\boldsymbol{h}_j(t) \qquad\text{with}\quad\boldsymbol{h}_j(t)=(h_{j,x},h_{j,y})^T\mathrel{\vcenter{:}}=\boldsymbol{X}^h(\rho_j,t)-\boldsymbol{X}^h(\rho_{j-1}, t).\] And we always assume that \(\displaystyle h_{\min}(t)=\min_{1\leq j\leq N}|\boldsymbol{h}_j(t)|>0\) for all \(t>0\).

By using \(\boldsymbol{h}_j\), the discrete geometric quantities such as the unit tangential vector \(\boldsymbol{\tau}^h\), the outward unit normal vector \(\boldsymbol{n}^h\) and the inclination angle \(\theta^h\) can be computed on each segment as: \[\boldsymbol{\tau}^h|_{I_j}=\frac{\boldsymbol{h}_j}{|\boldsymbol{h}_j|}\mathrel{\vcenter{:}}=\boldsymbol{\tau}_j^h,\qquad\boldsymbol{n}^h|_{I_i}=-(\boldsymbol{\tau}_j^h)^\perp=-\frac{(\boldsymbol{h}_j)^\perp}{|\boldsymbol{h}_j|}\mathrel{\vcenter{:}}=\boldsymbol{n}_j^h;\] and \[\theta^h|_{I_j}\mathrel{\vcenter{:}}=\theta_j^h,\qquad\text{satisfying}\quad\cos\theta^h_j=\frac{h_{j,x}}{|\boldsymbol{h}_j|},\quad\sin\theta_j^h=\frac{h_{j,y}}{|\boldsymbol{h}_j|}.\]

We introduce the finite element subspaces \[\begin{align} &\mathbb{K}^h\mathrel{\vcenter{:}}=\left\{u^h\in C(\mathbb{I})\mid u^h|_{I_j}\in\mathcal{P}^1(I_j),\,\,\forall j=1,2,\dots,N\right\}\subseteq H^1(\mathbb{I}),\\ &\mathbb{K}_p^h\mathrel{\vcenter{:}}=\{u^h\in\mathbb{K}^h\mid u^h(0)=u^h(1)\},\qquad\mathbb{X}_p^h\mathrel{\vcenter{:}}=[H_p^1(\mathbb{I})]^2, \end{align}\] where \(\mathcal{P}^1(I_j)\) is the set of polynomials defined on \(I_j\) of degree \(\leq 1\). For \(u,v\in\mathbb{K}^h\), the mass-lumped inner product \(\Bigl(\cdot,\cdot\Bigr)^h_{\Gamma^h(t)}\) with respect to \(\Gamma^h(t)\) is defined as \[\label{eqn:mass32lumped32inner32product} \Bigl(u, v\Bigr)^h_{\Gamma^h(t)}\mathrel{\vcenter{:}}=\frac{1}{2}\sum_{j=1}^N |\boldsymbol{h}_j(t)|\,\left((u\cdot v)(\rho_{j-1}^+)+(u\cdot v)(\rho_j^-)\right).\tag{26}\] where \(u(\rho_j^\pm)=\lim\limits_{\rho\to \rho_j^\pm} u(\rho)\). And the discretized differential operator \(\partial_s\) for \(f\in\mathbb{K}^h\) is defined as \[\partial_s f|_{I_j}\mathrel{\vcenter{:}}=\frac{f(\rho_j)-f(\rho_{j-1})}{|\boldsymbol{h}_j|}.\] The above definitions also hold true for vector-valued functions.

We now propose the spatial semi-discretization for 22 as follows: Let \(\Gamma_0^h\mathrel{\vcenter{:}}= \boldsymbol{X}^h(\cdot,0)\in\mathbb{X}^h_p,\mu(\cdot)\in\mathbb{K}^h_p\) be the approximations of \(\Gamma_0\mathrel{\vcenter{:}}= \boldsymbol{X}_0(\cdot),\mu_0(\cdot)\), respectively, for \(t>0\), find the solution \(\left(\boldsymbol{X}^h(\cdot,t),\mu^h(\cdot)\right)\in\mathbb{X}_p^h\times\mathbb{K}_p^h\) such that \[\label{spatial32semi-discretization32surface32diffusion} \begin{align} &\Bigl(\boldsymbol{n}^h\cdot \partial_t\boldsymbol{X}^h, \varphi^h\Bigr)_{\Gamma^h(t)}^h+\Bigl(\partial_s\mu^h, \partial_s\varphi^h\Bigr)_{\Gamma^h(t)}^h=0,\qquad \forall\varphi^h\in \mathbb{K}_p^h,\\ &\Bigl(\mu^h\boldsymbol{n}^h,\boldsymbol{\omega}^h\Bigr)_{\Gamma^h(t)}^h-\Bigl(\hat{\boldsymbol{G}}_k(\theta^h)\partial_s\boldsymbol{X}^h,\partial_s\boldsymbol{\omega}^h\Bigr)_{\Gamma^h(t)}^h=0,\qquad\forall\boldsymbol{\omega}^h\in\mathbb{X}_p^h, \end{align}\tag{27}\] where \[\hat{\boldsymbol{G}}_k(\theta^h)=\left(\begin{array}{cc} \hat{\gamma}(\theta^h) & -\hat{\gamma}^\prime(\theta^h)\\ \hat{\gamma}^\prime(\theta^h) & \hat{\gamma}(\theta^h) \end{array}\right)+k(\theta^h)\left(\begin{array}{cc} \sin^2\theta^h & -\cos\theta^h\sin\theta^h\\ -\cos\theta^h\sin\theta^h & \cos^2\theta^h \end{array}\right).\]

Denote the enclosed area and the free energy of the polygonal curve \(\Gamma^h(t)\) as \(A^h_c(t)\) and \(W^h_c(t)\), respectively, which are given by \[\label{discrete32area} \begin{align} &A^h_c(t)\mathrel{\vcenter{:}}=\frac{1}{2}\sum_{j=1}^N(x_j^h(t)-x_{j-1}^h(t))(y_j^h(t)+y_{j-1}^h(t)),\\ &W^h_c(t)\mathrel{\vcenter{:}}=\sum_{j=1}^N|\boldsymbol{h}_j(t)|\hat{\gamma}(\theta^h_j), \end{align}\tag{28}\] where \(x_j^h(t)\mathrel{\vcenter{:}}= x^h(\rho_j, t), y_j^h(t)\mathrel{\vcenter{:}}= y^h(\rho_j, t),\,\, \forall 0\leq j\leq N\).

Following similar steps in Proposition 1, it can be proved that the two geometric properties for the semi-discretization 27 still preserves:

3.3 A structure-preserving SPFEM discretization↩︎

Let \(\tau\) be the uniform time step. Denoting the approximation of \(\Gamma(t)=\boldsymbol{X}(\cdot,t)\) at \(t_m=m\tau, m=0,1,\dots,\) as \(\Gamma^m=\boldsymbol{X}^m(\cdot)=\cup_{j=1}^N\boldsymbol{h}_j^m\) where \(\boldsymbol{h}_j^m\mathrel{\vcenter{:}}=\boldsymbol{X}^m(\rho_j)-\boldsymbol{X}^m(\rho_{j-1})\). Then the definitions of the mass lumped inner product \((\cdot,\cdot)^h_{\Gamma^m}\), the unit tangential vector \(\boldsymbol{\tau}^m\), the unit outward normal vector \(\boldsymbol{n}^m\), and the inclination angle \(\theta^m\) with respect to \(\Gamma^m\) can be given in a similar approach.

Following the ideas in [38], [40], [42], [48] to design an SP-PFEM for surface diffusion, we utilize the explicit-implicit Euler method in time. The derived fully-implicit structure-preserving discretization of SPFEM for the anisotropic surface diffusion 7 is expressed as follows:

Suppose the initial approximation \(\Gamma^0(\cdot)\in \mathbb{X}^h\) is given by \(\boldsymbol{X}^0(\rho_j)=\boldsymbol{X}_0(\rho_j), \forall 0\leq j\leq N\). For any \(m=0, 1, 2, \ldots\), find the solution \((\boldsymbol{X}^{m+1}(\cdot)=(x^{m+1}(\cdot),y^{m+1}(\cdot))^T, \mu^{m+1}(\cdot))\in \mathbb{X}^h_p\times \mathbb{K}^h_p\) such that \[\label{eqn:sp-pfem32surface32diffusion} \begin{align} &\Bigl(\boldsymbol{n}^{m+\frac{1}{2}}\cdot \partial_t\boldsymbol{X}^{m+1}, \varphi^h\Bigr)_{\Gamma^m}^h+\Bigl(\partial_s\mu^{m+1}, \partial_s\varphi^h\Bigr)_{\Gamma^m}^h=0,\qquad \forall\varphi^h\in \mathbb{K}_p^h,\\ &\Bigl(\mu^{m+1}\boldsymbol{n}^{m+\frac{1}{2}},\boldsymbol{\omega}^h\Bigr)_{\Gamma^m}^h-\Bigl(\hat{\boldsymbol{G}}_k(\theta^m)\partial_s\boldsymbol{X}^{m+1},\partial_s\boldsymbol{\omega}^h\Bigr)_{\Gamma^m}^h=0,\qquad\forall\boldsymbol{\omega}^h\in\mathbb{X}_p^h, \end{align}\tag{29}\] where \[\hat{\boldsymbol{G}}_k(\theta^m)=\left(\begin{array}{cc} \hat{\gamma}(\theta^m) & -\hat{\gamma}^\prime(\theta^m)\\ \hat{\gamma}^\prime(\theta^m) & \hat{\gamma}(\theta^m) \end{array}\right)+k(\theta^m)\left(\begin{array}{cc} \sin^2\theta^m & -\cos\theta^m\sin\theta^m\\ -\cos\theta^m\sin\theta^m & \cos^2\theta^m \end{array}\right),\] and \[\label{dis32surface32energy32matrix} \boldsymbol{n}^{m+\frac{1}{2}}\mathrel{\vcenter{:}}= -\frac{1}{2}\frac{1}{|\partial_\rho \boldsymbol{X}^m|}(\partial_\rho \boldsymbol{X}^m+\partial_\rho\boldsymbol{X}^{m+1})^\perp.\tag{30}\]

3.4 Main results↩︎

Denote the enclosed area and the free energy of the polygon \(\Gamma^m\) as \(A^m_c\) and \(W^m_c\), respectively, which are given by \[\tag{31} \begin{align} \tag{32} A^m_c&\mathrel{\vcenter{:}}=\frac{1}{2}\sum_{j=1}^N\left(x^m(\rho_j)-x^m(\rho_{j-1})\right)\left(y^m(\rho_j) +y^m(\rho_{j-1})\right),\\ \tag{33} W^m_c&\mathrel{\vcenter{:}}=\sum_{j=1}^N|\boldsymbol{h}_j^m|\hat{\gamma}(\theta^m_j). \end{align}\]

In practical applications, it’s common to encounter situations where \(\hat{\gamma}(\theta)\) lacks high regularity. In the following sections, we always assume that \(\hat{\gamma}(\theta)\) is globally \(C^1\) and piecewise \(C^2\) on \(2\pi\mathbb{T}\).

We thus introduce the following energy stable conditions on \(\hat{\gamma}(\theta)\):

The main result of this work is the following structure-preserving property of the SPFEM 29 :

The proof of area conservation part is analogous to [38], we omit here for brevity. And the energy dissipation will be established in the next section.

4 Unconditionally energy stability of the SPFEM↩︎

The key point in proving energy dissipation of 29 is to establish an energy estimate akin to \[\left(\hat{\boldsymbol{G}}_k(\theta^m)\partial_s\boldsymbol{X}^{m+1},\partial_s(\boldsymbol{X}^{m+1}-\boldsymbol{X}^m)\right)_{\Gamma^m}^h\geq W_c^{m+1}-W_c^m,\] for controlling the energy difference between two subsequent time steps with the surface energy matrix \(\hat{\boldsymbol{G}}_k\). To achieve desired inequality, we need a local version of the estimate, which is formulated by the following lemma:

4.1 The minimal stabilizing function and its properties↩︎

To prove the local energy estimate ?? , the following two auxiliary functions are introduced as: \[\label{eqn:auxiliary32functions} \begin{align} &P_\alpha(\phi,\theta)\mathrel{\vcenter{:}}= 2\sqrt{\hat{\gamma}^2(\theta)+\alpha\hat{\gamma}(\theta)\sin^2\phi},\,\,\,\,\,\quad\quad\quad\qquad\forall\phi\in2\pi\mathbb{T},\\ &Q(\phi,\theta)\mathrel{\vcenter{:}}=\hat{\gamma}(\theta-\phi)+\hat{\gamma}(\theta)\cos\phi+\hat{\gamma}^\prime(\theta)\sin\phi,\qquad\forall\phi\in2\pi\mathbb{T}. \end{align}\tag{34}\] With the help of \(P_\alpha,Q\), we present the definition of minimal stabilizing function as follows: \[\label{eqn:def32minimal32stabilizing32function} k_0(\theta)\mathrel{\vcenter{:}}=\inf\left\{\alpha\geq 0\mid P_\alpha(\phi,\theta)-Q(\phi,\theta)\geq 0,\,\,\forall\phi\in 2\pi\mathbb{T}\right\}.\tag{35}\]

The following theorem guarantees the existence of \(k_0(\theta)\):

The proof of Theorem 3 will be presented in Section 5.
Once the \(\hat{\gamma}(\theta)\) is given, the minimal stabilizing function \(k_0(\theta)\) is determined, inducing a mapping from \(\hat{\gamma}(\theta)\) to \(k_0(\theta)\). Moreover, this mapping is sublinear, i.e., it is positive homogeneity and subadditivity.

The proof is similar to [42], we omit here for brevity.

4.2 Proof of the local energy estimate↩︎

4.3 Proof of the main result↩︎

Leveraging the local energy estimate ?? as outlined in Lemma 2, we can now prove the unconditional energy stability aspect of the main result, as stated in Theorem 2.

5 Existence of the minimal stabilizing function↩︎

In this section, we present a proof of the existence of the minimal stabilizing function corresponding to \(\hat{\gamma}(\theta)\) that satisfies the energy stable condition in ?? .

Recall the definition 34 of \(P_\alpha,Q\), denote \[\label{eqn:F61P942-Q942} \begin{align} F_\alpha(\phi,\theta)&\mathrel{\vcenter{:}}= P_\alpha^2(\phi,\theta)-Q^2(\phi,\theta)\\ &=4\hat{\gamma}(\theta)\left(\hat{\gamma}(\theta)+\alpha\sin^2\phi\right)\\ &\quad-\left(\hat{\gamma}(\theta-\phi)+\hat{\gamma}(\theta)\cos\phi+\hat{\gamma}^\prime(\theta)\sin\phi\right)^2. \end{align}\tag{40}\] Suppose \(C>0\) is a positive number such that \(\frac{1}{C}\leq\hat{\gamma}(\theta)\leq C\) and \(|\hat{\gamma}^\prime(\theta)|,|\hat{\gamma}^{\prime\prime}(\theta)|\leq C,\,\,\forall\theta\in2\pi\mathbb{T}\).

Before formally commencing our proof, we first need the following two technical lemmas:

 Similar to derivations of @eq:eqn32taylor32gamma , according to the
 mean value theorem, we know that there exists positive number
 $0<\varepsilon<\frac{\pi}{2}$ and a neighborhood
 $U_\pi=\{\phi\mid|\phi-\pi|<\varepsilon\}$. Such that for any
 $\phi\in U_\pi$, there exists a $\xi$ between $\theta^*-\phi$ and
 $\theta^*-\pi$ satisfying
 $$\label{eqn:taylor32of32gamma2} \begin{align} \hat{\gamma}(\theta^*-\phi)&=\hat{\gamma}(\theta^*-\pi)-\hat{\gamma}^\prime(\theta^*-\pi)(\phi-\pi)+\frac{\hat{\gamma}^{\prime\prime}(\xi)}{2}(\phi-\pi)^2\\ &=\hat{\gamma}(\theta^*-\pi)+\frac{\hat{\gamma}^{\prime\prime}(\xi)}{2}(\phi-\pi)^2. \end{align}$$ {#eq:eqn:taylor32of32gamma2} 

 Again, by the mean value theorem, we know that there exists a
 $\xi_1$ between $\pi,\phi$,
 $$\label{eqn:taylor32of32cos2} \cos\phi=-1-\frac{\cos\xi_1}{6}(\phi-\pi)^2.$$ {#eq:eqn:taylor32of32cos2} 
 By substituting @eq:eqn:taylor32of32gamma2 and
 @eq:eqn:taylor32of32cos2 into
 $2\hat{\gamma}(\theta^*)-Q(\phi,\theta^*)$, we have
 $$\begin{align} 2\hat{\gamma}(\theta^*)-Q(\phi,\theta^*)&=2\hat{\gamma}(\theta^*)-\hat{\gamma}^\prime(\theta^*-\pi)-\frac{\hat{\gamma}^{\prime\prime}(\xi)}{2}(\phi-\pi)^2\\ &\quad +\hat{\gamma}(\theta^*)+\frac{\cos\xi_1}{6}\hat{\gamma}(\theta)(\phi-\pi)^2+\hat{\gamma}(\theta^*)\sin\phi\\ &=\frac{1}{6}\left(\hat{\gamma}(\theta^*)\cos\xi_1-3\hat{\gamma}^{\prime\prime}(\xi)\right)(\phi-\pi)^2. \end{align}$$
 Thus,
 $$\begin{align} |2\hat{\gamma}(\theta^*)-Q(\phi,\theta^*)|&\leq \frac{1}{6}\left(|\hat{\gamma}(\theta^*)\cos\xi_1|+3|\hat{\gamma}^{\prime\prime}(\xi)|\right)|\phi-\pi|^2\\ &\leq \frac{2C}{3}|\phi-\pi|^2. \end{align}$$

 Given that $|2\hat{\gamma}(\theta^*)+Q(\phi,\theta^*)|\leq 5C$, we
 can deduce that
 $$\begin{align} |4\hat{\gamma}^2(\theta^*)-Q^2(\phi,\theta^*)|&\leq |2\hat{\gamma}(\theta^*)+Q(\phi,\theta^*)||2\hat{\gamma}(\theta^*)-Q(\phi,\theta^*)|\\ &\leq \frac{10C^2}{3}|\phi-\pi|^2. \end{align}$$
 Noting that for any $\phi\in U_{\pi}$, we have
 $|\sin\phi|>\frac{2}{\pi}|\phi-\pi|$ since $\varepsilon\in(0,\pi)$.
 Therefore,
 $$\begin{align} F_\alpha(\phi,\theta^*)&=4\hat{\gamma}(\theta^*)\alpha\sin^2\phi+4\hat{\gamma}^2(\theta^*)-Q^2(\phi,\theta^*)\\ &\geq \frac{16}{C\pi^2}\alpha|\phi-\pi|^2-\frac{10C^2}{3}|\phi-\pi|^2\\ &\geq\left[\frac{16}{C\pi^2}\alpha-\frac{10C^2}{3}\right]|\phi-\pi|^2,\,\,\forall\phi\in U_{\pi}. \end{align}$$
 Therefore, there exists a positive constant
 $k_{U_\pi}\mathrel{\vcenter{:}}=\frac{5\pi^2C^3}{24}$, such that for
 $\alpha>k_{U_\pi}$, we have
 $F_{\alpha}(\phi,\theta)\geq 0,\forall\phi\in U_\pi$.

With the help of above two lemmas, Theorem 3 can now be proven:

6 Extension to solid-state dewetting↩︎

In this section, we extend the conservative formulation 17 and its SPFEM 29 for a closed curve under anisotropic surface diffusion to solid-state dewetting in materials science [12], [22], [37].

6.1 Sharp interface model and a stabilized formulation↩︎

As shown in Fig. 2, the solid-state dewetting problem in 2D is described as evolution of an open curve \(\Gamma(t)=\boldsymbol{X}(s,t)=(x(s,t),y(s,t))^T\,\,(0\leq s\leq L)\) under anisotropic surface diffusion and contact line migration. Here, \(s\) and \(t\) are the arc-length parameter and time, respectively, and \(L\mathrel{\vcenter{:}}= L(t)=|\Gamma(t)|\) represents the total length of \(\Gamma(t)\). As it was derived in the literature [12], [22], [37], a dimensionless sharp-interface model for simulating solid-state dewetting of thin films with weakly anisotropic surface energy can be formulated as: \(\boldsymbol{X}(s,t)\) satisfying anisotropic surface diffusion 1 –2 with following boundary conditions

1. contact point condition \[\label{bdc32contact32points} y(0,t)=0,\qquad y(L,t)=0,\qquad t\geq 0;\tag{43}\]

2. relaxed contact angle condition \[\label{bdc32relaxed32contact32angle} \frac{\mathrm{d}x_c^l(t)}{\mathrm{d}t}=\eta f(\theta^l_d;\sigma),\qquad \frac{\mathrm{d}x_c^r(t)}{\mathrm{d}t}=-\eta f(\theta^r_d;\sigma),\qquad t\geq 0;\tag{44}\]

3. zero-mass flux condition \[\label{bdc32zero-mass32flux} \partial_s\mu(0,t)=0,\qquad \partial_s\mu(L,t)=0,\qquad t\geq 0;\tag{45}\]

where \(x_c^l(t)=x(0,t)\leq x_c^r(t)=x(L,t)\) are the left and right contact points, \(\theta_d^l\mathrel{\vcenter{:}}=\theta_d^l(t),\theta_d^r\mathrel{\vcenter{:}}= \theta_d^r(t)\) are the contact angles at each contact points, respectively. \(f(\theta;\sigma)\) is defined as \[f(\theta;\sigma)\mathrel{\vcenter{:}}=\hat{\gamma}(\theta)\cos\theta-\hat{\gamma}^\prime(\theta)\sin\theta-\sigma,\qquad\theta\in 2\pi\mathbb{T}.\] with \(\sigma\mathrel{\vcenter{:}}=\cos\theta_Y=\frac{\gamma_{VS}-\gamma_{FS}}{\gamma_{FV}}\) and \(\theta_Y\) be the isotropic Young contact angle. \(0<\eta<+\infty\) denotes the contact line mobility [12], [22], [37].

Similarly to 17 , a stabilized formulation for solid-state dewetting can be derived as: \[\tag{46} \begin{align} &\boldsymbol{n}\cdot\partial_t\boldsymbol{X}-\partial_{ss}\mu=0,\qquad 0<s<L(t),\qquad t>0,\tag{47}\\ &\mu\boldsymbol{n}+\partial_s\Bigl(\hat{\boldsymbol{G}}_k(\theta)\partial_{s}\boldsymbol{X}\Bigr)=\boldsymbol{0},\tag{48} \end{align}\] with boundary conditions 43 –45 , where \(\boldsymbol{\hat{G}}_k(\theta)\) is the stabilizing surface energy matrix given in ?? .

6.2 A stabilized weak formulation↩︎

Suppose \(\Gamma(t)\) is parameterized by a time-independent variable \(\rho\) over a fixed domain \(\rho\in\mathbb{I}=[0,1]\), i.e. \[\Gamma(t)\mathrel{\vcenter{:}}=\boldsymbol{X}(\rho,t)=(x(\rho,t),y(\rho,t))^T\colon\mathbb{I}\times\mathbb{R}^+\to\mathbb{R}^2,\] thus the arc-length parameter \(s\) can be given as \(s=\int_0^\rho|\partial_q\boldsymbol{X}(q,t)|\,\mathrm{d}q\). (Again, we do not distinguish \(\boldsymbol{X}(s,t)\) and \(\boldsymbol{X}(\rho,t)\).)

Introducing the following functional spaces as \[H^1_0(\mathbb{I})\mathrel{\vcenter{:}}=\left\{u\in H^1(\mathbb{I})\mid u(0)=u(1)=0\right\},\qquad\mathbb{X}\mathrel{\vcenter{:}}= H^1(\mathbb{I})\times H^1_0(\mathbb{I}).\] Following similar derivations as in Section 3, we present the new variational formulation for 46 with boundary conditions 43 –45 as follows: Given an initial open curve \(\Gamma(0)=\boldsymbol{X}(\cdot,0)=\boldsymbol{X}_0\in\mathbb{X}\), find the solution \((\boldsymbol{X}(\cdot,t),\mu(\cdot,t))\in\mathbb{X}\times H^1(\mathbb{I})\) such that \[\label{new32variational32formulation32solid-state32dewetting} \begin{align} \begin{aligned} &\Bigl(\boldsymbol{n}\cdot \partial_t\boldsymbol{X}, \varphi\Bigr)_{\Gamma(t)}+\Bigl(\partial_s\mu, \partial_s\varphi\Bigr)_{\Gamma(t)}=0,\qquad \forall\varphi\in H^1(\mathbb{I}), \end{aligned}\\ \begin{align} &\Bigl(\mu\boldsymbol{n},\boldsymbol{\omega}\Bigr)_{\Gamma(t)}-\Bigl(\hat{\boldsymbol{G}}_k(\theta)\partial_s\boldsymbol{X},\partial_s\boldsymbol{\omega}\Bigr)_{\Gamma(t)}-\frac{1}{\eta}\left[\frac{\mathrm{d}x_c^l(t)}{\mathrm{d}t}\omega_1(0)+\frac{\mathrm{d}x_c^r(t)}{\mathrm{d}t}\omega_1(1)\right]\\ &\qquad\qquad\,\,\, +\sigma\left[\omega_1(1)-\omega_1(0)\right]=0,\qquad\forall\boldsymbol{\omega}=(\omega_1,\omega_2)^T\in\mathbb{X}. \end{align} \end{align}\tag{49}\]

The area \(A_o(t)\) enclosed by \(\Gamma(t)\) and the substrate and the total energy \(W_o(t)\) are defined as \[A_o(t)\mathrel{\vcenter{:}}= \int_{\Gamma(t)}y(s,t)\partial_sx(s,t)\,\mathrm{d}s,\quad W_o(t)\mathrel{\vcenter{:}}= \int_{\Gamma(t)}\hat{\gamma}(\theta)\,\mathrm{d}s-\sigma\left(x_c^r(t)-x_c^l(t)\right).\] Similar to the closed curve case, the weak formulation 49 still preserves area conservation and energy dissipation.

The proof is similar to Proposition 1, we omit the details here for brevity.

6.3 A SPFEM and its proprieties↩︎

Introducing the finite element spaces \[\mathbb{K}_0^h\mathrel{\vcenter{:}}= \left\{u^h\in\mathbb{K}^h\mid u^h(0)=u^h(1)=0\right\},\qquad\mathbb{X}^h\mathrel{\vcenter{:}}= \mathbb{K}^h\times\mathbb{K}^h_0.\] The semi-discretization in space of 49 and its area conservation and energy dissipation properties are similar to Section 3, we omit the details here for brevity.

Let \(\tau\) be the uniform time step, denotes the approximation of \(\Gamma(t)=\boldsymbol{X}(\cdot,t)\) at \(t_m=m\tau,m=0,1,\dots\) as \(\Gamma^m=\boldsymbol{X}^m(\cdot)=\cup_{j=1}^N\boldsymbol{h}_j^m\) where \(\boldsymbol{h}_j^m\mathrel{\vcenter{:}}= \boldsymbol{X}^m(\rho_j)-\boldsymbol{X}^m(\rho_{j-1})\). Then the structure-preserving discretization of SPFEM for solid-state dewetting 46 can be stated as: Suppose the initial approximation \(\Gamma^0(\cdot)\in \mathbb{X}^h\) is given by \(\boldsymbol{X}^0(\rho_j)=\boldsymbol{X}_0(\rho_j), \forall 0\leq j\leq N\). For any \(m=0, 1, 2, \ldots\), find the solution \((\boldsymbol{X}^{m+1}(\cdot)=(x^{m+1}(\cdot),y^{m+1}(\cdot))^T, \mu^{m+1}(\cdot))\in \mathbb{X}^h\times \mathbb{K}^h\), such that \[\label{eqn:sp-pfem32solid-state32dewetting} \begin{align} \begin{aligned} &\Bigl(\boldsymbol{n}^{m+\frac{1}{2}}\cdot\frac{\boldsymbol{X}^{m+1}-\boldsymbol{X}^m}{\tau}, \varphi^{h}\Bigr)_{\Gamma^m}^h+\Bigl(\partial_s\mu^{m+1}, \partial_s\varphi^h\Bigr)_{\Gamma^m}^h=0,\qquad \forall\varphi^h\in \mathbb{K}^h, \end{aligned}\\ \begin{align} &\Bigl(\mu^{m+1}\boldsymbol{n}^{m+\frac{1}{2}},\boldsymbol{\omega}^h\Bigr)_{\Gamma^m}^h-\Bigl(\hat{\boldsymbol{G}}_k(\theta^m)\partial_s\boldsymbol{X}^m,\partial_s\boldsymbol{\omega}^h\Bigr)_{\Gamma^m}^h-\frac{1}{\eta}\left[\frac{x_l^{m+1}-x_l^m}{\tau}\omega_1^h(0)+\frac{x_r^{m+1}-x_r^m}{\tau}\omega_1^h(1)\right] \\ &\qquad\qquad\qquad\,\,\,+\sigma\left[\omega_1^h(1)-\omega_1^h(0)\right]=0,\qquad\forall\boldsymbol{\omega}^h=(\omega_1^h,\omega_2^h)\in\mathbb{X}^h, \end{align} \end{align}\tag{50}\] satisfying \(x_l^{m+1}=x^{m+1}(0)\leq x_r^{m+1}=x^{m+1}(1)\). The definition of \(\boldsymbol{G}_k(\theta^m)\) is similar to 30 .

Denote the area \(A_o^m\) enclosed by the open polygonal curve \(\Gamma(t)\) and the substrate and the total surface energy \(W_o^m\) as \[\begin{align} A_o^m&\mathrel{\vcenter{:}}= \frac{1}{2}\sum_{j=1}^N\left(x^m(\rho_j)-x^m(\rho_{j-1})\right)\left(y^m(\rho_j)+y^m(\rho_{j-1})\right),\\ W_o^m&\mathrel{\vcenter{:}}= \sum_{j=1}^N|\boldsymbol{h}_j^m|\hat{\gamma}(\theta^m_j)-\sigma\left(x_r^m-x_l^m\right). \end{align}\] Then for the SPFEM 50 , we have following result on its structure preservation.

7 Numerical results↩︎

In this section, we report numerical experiments for the proposed SPFEM 29 and 50 for time evolutions of closed curves and open curves, respectively. Extensive results are provided to illustrate their efficiency, accuracy, area conservation and unconditional energy stability.

In the numerical tests, following typical anisotropic surface energies are considered in the simulations:

• Case I: \(\hat{\gamma}(\theta)=1+\beta\cos 3\theta\) with \(|\beta|<1\). It is weakly anisotropic when \(|\beta|<\frac{1}{8}\) and strongly anisotropic otherwise;
  

• Case II: \(\hat{\gamma}(\theta)=\sqrt{1+b\cos^2\theta}\) with \(b>-1\).

To compute the minimal stabilizing function \(k_0(\theta)\), we solve the optimization problem 35 for 20 uniformly distributed points in \([-\pi,\pi]\) and do linear interpolation for the intermediate points.

To verify the quadratic convergence rate in space and linear convergence rate in time, the time step \(\tau\) is always chosen as \(\tau=16h^2\) except it is stated otherwise. The manifold distance [41], [49] \[M(\Gamma_1,\Gamma_2)\mathrel{\vcenter{:}}= 2|\Omega_1\cup\Omega_2|-|\Omega_1|-|\Omega_2|,\] is employed to measure the distance between two closed curves \(\Gamma_1,\Gamma_2\), where \(\Omega_1,\Omega_2\) are the regions enclosed by \(\Gamma_1,\Gamma_2\) and \(|\Omega|\) represents the area of \(|\Omega|\). Suppose \(\Gamma^m\) is the numerical approximation of \(\Gamma^h(t=t_m\mathrel{\vcenter{:}}= m\tau)\), thus the numerical error is defined as \[e^h(t)\Big|_{t=t_m}\mathrel{\vcenter{:}}= M(\Gamma^m,\Gamma(t=t_m)).\] In the Newton’s iteration, the tolerance value is set to be \(\text{tol}=10^{-12}\).

To test the mesh quality, the energy stability and area conservation numerically, we introduce the following indicators: the weighted mesh ratio \[R_{\gamma}^h(t)\mathrel{\vcenter{:}}= \frac{\max\limits_{1\leq j\leq N}\hat{\gamma}(\theta_j)|\boldsymbol{h}_j|}{\min\limits_{1\leq j\leq N}\hat{\gamma}(\theta_j)|\boldsymbol{h}_j|},\] the normalized area loss and the normalized energy for closed curves: \[\left.\frac{\Delta A^h_c(t)}{A^h_c(0)}\right|_{t=t_m}\mathrel{\vcenter{:}}= \frac{A^m_c-A^0_c}{A_c^0},\qquad\left.\frac{W^h_c(t)}{W^h_c(0)}\right|_{t=t_m}\mathrel{\vcenter{:}}= \frac{W^m_c}{W^0_c},\] and for open curves: \[\left.\frac{\Delta A^h_o(t)}{A^h_o(0)}\right|_{t=t_m}\mathrel{\vcenter{:}}= \frac{A^m_o-A^0_o}{A_o^0},\qquad\left.\frac{W^h_o(t)}{W^h_o(0)}\right|_{t=t_m}\mathrel{\vcenter{:}}= \frac{W^m_o}{W^0_o}.\]

In the following numerical tests, the initial shapes are chosen as a complete and a half ellipse with major axis \(4\) and minor axis \(1\) for closed curves and open curves, respectively, unless stated otherwise. The exact solution \(\Gamma(t)\) is approximated by choosing \(k(\theta)=k_0(\theta)\) with a small mesh size \(h_e=2^{-8}\) and a time step \(\tau_e=2^{-12}\) in 29 . For solid-state dewetting problems, we always choose the contact line mobility \(\eta=100\).

7.1 Results for closed curves↩︎

Fig. 3 plots the convergence rates of the proposed SPFEM 29 for: (a) the \(3\)-fold anisotropy \(\hat{\gamma}(\theta)=1+\beta\cos 3\theta\) with different anisotropic strengths \(\beta\) under a fixed time \(t=0.5\); (b) the ellipsoidal anisotropy \(\hat{\gamma}(\theta)=\sqrt{1-0.8\cos^2\theta}\) at different times. It clearly demonstrates that the second-order spatial convergence remains consistent regardless of anisotropies and computational times, suggesting a high level of robustness in the convergence rate.

Fig. 4 exhibits that the weighted mesh ratio \(R_{\gamma}^h\) converge to constants as \(t\to+\infty\). This suggests an asymptotic quasi-uniform mesh distribution of the proposed SPFEM 29 .

The time evolutions of the normalized area loss \(\frac{\Delta A^h_c(t)}{A^h_c(0)}\), the number of the Newton’s iteration with \(h=2^{-7},\tau=2^{-10}\) are given in Fig. 5. And the normalized energy \(\frac{W^h_c(t)}{W^h_c(t)}\) with different \(h\) are summarized in Fig. 6.

The observation from Fig. 5–Fig. 6 reveals that:

1. The normalized area loss is at \(10^{-16}\), aligns closely with the order of the round-off error (cf. Fig. 5). This observation affirms the practical preservation of area in simulations.
   

2. The numbers of the Newton’s iteration are initially \(3\) or \(4\), and quickly descend to \(2\) (cf. Fig. 5). This discovery indicates that the proposed SPFEM 29 can be solved with high efficiency, requiring only a few iterations.
   

3. The normalized energy is monotonically decreasing when \(\hat{\gamma}(\theta)\) satisfies the energy stable conditions ?? in Definition 1 (cf. Fig. 6). Results in Fig. 6 (a) shows that the proposed SPFEM 29 still preserves good energy stability properties when \(\beta\) takes its maximum value of \(1/2\) in Remark 5. And results in Fig. 6 (b) indicate that, unlike the ES-PFEM in [41], the SPFEM remains unconditionally energy stable when \(b<-1/2\) as stated in Remark 6.

7.2 Results for open curves in solid-state dewetting↩︎

Fig. 7 plots the computation errors of the proposed SPFEM 50 for: (a) the \(3\)-fold anisotropy \(\hat{\gamma}(\theta)=1+\beta\cos 3\theta\) with different anisotropic strengths \(\beta\) under a fixed time \(t=0.5\); (b) the ellipsoidal anisotropy \(\hat{\gamma}(\theta)=\sqrt{1+2\cos^2\theta}\) at different times. The results verify the quadratic convergence rate for the proposed SPFEM 50 .

In Fig. 8, the weighted mesh ratios \(R_\gamma^h\) tend to constants as \(t\to+\infty\), showing that the SPFEM 50 still possesses the asymptotic quasi-uniform distribution.

Time evolutions of the normalized area loss \(\frac{\Delta A^h_o(t)}{A^h_o(0)}\), the number of the Newton’s iteration with \(h=2^{-7},\tau=2^{-10}\) are presented in Fig. 9. And the normalized energy \(\frac{W^h_o(t)}{W^h_o(t)}\) with different \(h\) are illustrated in Fig. 10.

It can be observed from Fig. 6–Fig. 9 that:

1. The normalized area loss is about \(10^{-16}\) at the same order of the round-off error (cf. Fig. 9), verifying that the area is conserved up to the machine precision.
   

2. The numbers of the Newton’s iteration are initially \(3\) and finally \(2\) (cf. Fig. 9). This finding suggests that, despite the fully-implicit nature of the proposed SPFEM 50 , it can be solved very efficiently.
   

3. The normalized energy is monotonically decreasing when \(\hat{\gamma}(\theta)\) satisfying ?? (cf. Fig. 10). In contrast to the ES-PFEM in [41], the proposed SPFEM 50 still guarantees the energy dissipation when \(\beta<\frac{1}{10}\) in Case I and \(b>1\) in Case II, as asserted by Remark 5 and Remark 6.

7.3 Application for morphological evolutions↩︎

Finally we apply the proposed SPFEMs 29 and 50 to simulate the morphological evolutions under the anisotropic surface diffusion. Results for both closed curves and open curves in solid-state dewetting problems are provided.

The morphological evolutions from the initial shapes to their numerical equilibriums are presented in Fig. 11–Fig. 13. For closed curve cases, the initial shape is an ellipse with major axis \(4\) and minor axis \(1\), while for open curve cases, it is an open \(4\times 1\) rectangle.

Fig. 11 plots the morphological evolutions of an ellipse with major axis \(4\) and minor axis \(1\) under anisotropic surface diffusion with four different surface energies: (a) anisotropy in Case I with \(\beta=\frac{1}{2}\), which attends the maximum value in Remark 5; (b) anisotropy in Case II with \(b=-0.8\); (c) the \(4\)-fold anisotropy \(\hat{\gamma}(\theta)=1+\frac{1}{16}\cos 4\theta\) [22]; and (d) \(\hat{\gamma}(\theta)=\sqrt{\left(\frac{5}{2}+\frac{3}{2}\text{sgn}(n_1)\right)n_1^2+n_2^2}\) with \(\boldsymbol{n}=(n_1,n_2)^T=(-\sin\theta,\cos\theta)^T\) [26].

Results in Fig. 11 (b) and Fig. 11 (c) show that, compared to the ES-PFEM in [41], the proposed SPFEM 29 demonstrates a better performance over a broader range of parameters during evolutions. Fig. 11 (d) indicates that the SPFEM 29 also works well for a globally \(C^1\) and piecewise \(C^2\) anisotropy.

And Fig. 12 illustrates the morphological evolutions under anisotropic surface diffusion with different surface energies for a four-fold star shape curve \[\left\{\begin{array}{l} x=(1+0.5\cos 4\theta)\cos\theta, \\ y=(1+0.5\cos 4\theta)\sin\theta, \end{array}\right.\theta\in 2\pi\mathbb{T},\] and an ellipse (with major axis \(4\) and minor axis \(1\)) rotated counterclockwise by \(\frac{\pi}{10}\).

In Fig. 13, we display the morphological evolutions from an open \(4\times 1\) rectangular curve to their equilibriums shapes with different surface energies: (a) anisotropy in Case I with \(\beta=\frac{1}{9}\); (b) anisotropy in Case II with \(b=2\); (c) the \(4\)-fold anisotropy \(\hat{\gamma}(\theta)=1+\frac{1}{16}\cos 4\theta\); (d) \(\hat{\gamma}(\theta)=\sqrt{\left(\frac{5}{2}+\frac{3}{2}\text{sgn}(n_1)\right)n_1^2+n_2^2}\) with \(\boldsymbol{n}=(n_1,n_2)^T=(-\sin\theta,\cos\theta)^T\).

Similar to the closed curve cases, the SPFEM 50 extends the choices in surface energies for simulating solid-state dewetting (cf. Fig. 13 (a) – Fig. 13 (c)). And Fig. 13 (d) illustrates that our method also performs effectively for \(\hat{\gamma}(\theta)\) with lower regularity.

8 Conclusions↩︎

We propose a structure-preserving stabilized parametric finite element method (SPFEM) for the anisotropic surface diffusion. This method is subject to mild conditions on \(\hat{\gamma}(\theta)\), and works effectively for closed curves and open curves with contact line migration in solid-state dewetting. By introducing a new stabilized surface energy matrix, we obtain a conservative form and its weak formulation for anisotropic surface diffusion. Based on this weak formulation, a novel SPFEM is proposed by utilizing the PFEM for spatial discretization and the implicit-explicit Euler method for temporal discretization. To analyze the unconditional energy stability, we extend the framework proposed by Bao and Li to the \(\hat{\gamma}(\theta)\) formulation. This approach starts by defining the minimal stabilizing function, proving its existence, results in a local energy estimate, and subsequently establishes unconditional energy stability. Due to the very mild requirements on the surface energy \(\hat{\gamma}(\theta)\), the methods are able to simulate over a broader range of anisotropies for both closed curves and open curves. Moreover, the SPFEMs are applicable for the globally \(C^1\) and piecewise \(C^2\) anisotropy as well, which is a capability not possessed by other PFEMs.

CRediT authorship contribution statement↩︎

Yulin Zhang: Conceptualization, Methodology, Software, Investigation, Writing - Original Draft.
Yifei Li: Project administration, Supervision, Writing - Review & Editing.
Wenjun Ying: Supervision, Writing - Review & Editing.

Declaration of competing interest↩︎

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgement↩︎

We would like to specially thank Professor Weizhu Bao for his valuable suggestions and comments. The work of Zhang was partially supported by the State Scholarship Fund (No.202306230346) by the Chinese Scholar Council (CSC) and the Zhiyuan Honors Program for Graduate Students of Shanghai Jiao Tong University (No.021071910051). The work of Li was funded by the Ministry of Education of Singapore under its AcRF Tier 2 funding MOE-T2EP20122-0002 (A-8000962-00-00). Part of the work was done when the authors were visiting the Institute of Mathematical Science at the National University of Singapore in 2024.

↩︎

References↩︎