Doubly Degenerate Diffuse Interface Models of Anisotropic Surface Diffusion

We extend the doubly degenerate Cahn-Hilliard (DDCH) models for isotropic surface diffusion [arXiv:1909.04458], which yield more accurate approximations than classical degenerate Cahn-Hilliard (DCH) models, to the anisotropic case. We consider both weak and strong anisotropies and demonstrate the capabilities of the approach for these cases numerically. The proposed model provides a variational and energy dissipative approach for anisotropic surface diffusion, enabling large scale simulations with material-specific parameters.


INTRODUCTION
Surface diffusion is an important transport mechanism in materials science, e.g. in processes such as solid-state dewetting of semiconductors or coarsening of bulk nanoporous metals. In the isotropic setting, surface diffusion is modeled by a geometric evolution equation which relates the normal velocity of a hypersurface in Euclidean space to the surface Laplacian of the mean curvature 2 , e.g., = Δ Σ . In more realistic anisotropic settings the mean curvature is replaced by a weighted mean curvature, defined as the surface divergence of the Cahn-Hoffmann vector 3 , e.g., = ∇ Σ ⋅ , with = (̂ ) an anisotropic surface energy and̂ the outward-pointing surface normal. Even though some direct numerical approaches for these equations, in the isotropic and anisotropic setting exist -see, e.g., [4,5,6,7,8,9] -for most applications in materials science diffuse interface approximations are preferred -see, e.g., [10,11,12,13,14,15,16,17,18,19]. These diffuse interface approaches capture the motion of the interface implicitly as the evolution of an iso-surface of a phase-field function. Typically, they are fourthorder nonlinear diffusion equations of Cahn-Hilliard type, whose solutions formally converge to those of their sharp interface counterpart, as the interface thickness tends to zero 11,20,21,22 . They require a degenerate mobility function and are termed degenerate Cahn-Hilliard (DCH) equations. In the model proposed in [11], an additional degeneracy is introduced, following similar ideas as used for the thin film limit in classical phase field models for solidification 23 . We call such models doubly degenerate Cahn-Hilliard (DDCH) models. This second degeneracy does not alter the asymptotic limit 11 , but actually leads to more accurate surface diffusion approximations. See, e.g., the discussion in [24]. In fact several simulations for realistic applications in materials science consider this model -see, e.g., [25,26,27,28,29,30,31,32,33,34,35,36] -and claim that their large scale simulations would not be feasible without the additional degeneracy. The drawback of this DDCH model, which is termed RRV model in several publications 21,30,35 , is that it is non-variational. That is, there is no known free energy that is dissipated along solution trajectories. This makes it harder to prove properties of solutions and derive certain numerical stabilities. In [1] this problem was solved by introducing a new variational DDCH model for surface diffusion, which can be connected with the non-variational DDCH model of [11]. This was done for the isotropic case and will here be generalized to the anisotropic case.
The paper is organized as follows: In Sect. 2 we first review the different isotropic DDCH models, their connection, and then extend the idea to weak as well as strong anisotropies. The latter case requires an additional curvature regularisation, leading to higher-order equations; In Sect. 3 we illustrate the numerical approach which uses a semi-implicit time-integration scheme and adaptive Finite Elements for discretization in space; two-and three-dimensional numerical results are shown in Sect. 4, including comparisons between different models and with sharp-interface solutions. Different choices for surface-energy densities are reported, also matching real material properties, and applications are illustrated to show the wide applicability. In Sect. 5 we draw our conclusions.

Isotropic, variational DDCH model
Suppose that Ω ⊂ ℝ , = 2, 3, is a bounded, open set. Let us consider the free energy with = 4 2 (1 − ) 2 , = 72 a quartic, symmetric double well potential and 0 a singular function of the form 0 can be regularised so that it is defined, continuous, and differentiable for all : The dynamics is described as the −1 gradient flow by where = is the chemical potential with the variational derivative of the free energy with respect to . 0 is the degenerate mobility function. Its regularised form is and 0 is obtained setting = 0.

Γ(4−2 )
with Γ the usual (Bernoulli) Gamma-function. The proposed system, Eqs. (2) and (3), is a free energy dissipative dynamical system. With the asymptotic approximation 1 ( ) ≈ 2 |∇ | 2 -which holds when the interface has a hyperbolic tangent profile -we obtain ≈ 0 ( ) 1 ′ ( ) − 0 ( )∇ 2 , which simplifies Eq. (3) and leads to This model is strikingly similar to the isotropic version of the model in [11], but with one important caveat. The model in [11] corresponds to the choice = 2, which is not defined for Eqs. (2) -(3) 1 . Interestingly, Eqs. (5) -(6) can be defined in a reasonable way for any ∈ [0, ∞). The choice of will not be the same as for Eqs. (2) -(3), in general. However, for = 1 they are identical. See [1] for details. Numerical solutions indicate a higher accuracy for the variational DDCH model Eqs. (2) and (3) and the non-variational DDCH model Eqs. (5) and (6) if compared with the classical DCH model ( = 0 and = 1). Another advantage of the DDCH models is that deviations of from the pure phase values, 0 and 1, are smaller in a point-wise sense, when compared with the solutions of the classical DCH model. This property can further be elaborated to guarantee, in the non-regularised case = 0, that 0 ≤ ≤ 1 for DDCH approximations. So-called positivity preserving schemes of the variational DDCH model Eqs. (2) and (3), are possible and will be discussed in later publications.

Anisotropic, variational DDCH model
Instead of a constant surface-energy density = 1, most materials have anisotropic surface-energy density = (̂ ). In the phase-field context, is extended by defininĝ 0 = − ∇ |∇ | . In other words, is defined everywhere that ∇ ≠ . We consider the free energy given by In contrast to classical approaches 37 , which consider the anisotropic surface-energy density only in the gradient term of the free energy, we here follow the approach of [13], which ensures that the interface thickness is independent of orientation. As in [1], we include the energy restriction function 0 , which is singular at the pure phase values = 0, 1. The asymptotic analysis of [1] can be repeated, under the assumption that the diffuse interface has the usual hyperbolic tangent profile, to show that

Model regularisation
We note that (7) is not defined for arbitrary smooth functions . This can be an enormous problem for practical computations.
To fix this, we regularise the energy as follows. Set = ∇ . Define the regularised, extended normal and observe that this vector is now (i) only approximately of unit length and (ii) only approximately normal to the level sets of . The regularised free energy is defined as which is now well-defined for arbitrary smooth phase-field functions. We note that it is possible, and perhaps even advantageous in some settings, to use two separate regularisation parameters for̂ 0 and 0 , but we will avoid this technical discussion here and instead focus on formulation and numerics.
In the computation of the variational derivative of , denoted , we need the following calculation: , a regularised projection-like matrix. The regularised doubly degenerate energy dissipative flow now reads where = is the chemical potential. We point out that the chemical potential becomes singular and ill-defined when the regularisation is formally turned off ( = 0), in particular, when ∇ = = , = 0, or = 1.
Also for the anisotropic cases the asymptotic approximation 1 ( ) ≈ 2 |∇ | 2 can be used to simplify the equations and, in turn, their numerical integration. However, a one-to-one correspondence with the non-variational models proposed in [11] is not possible. The corresponding version, including the surface energy density as in [13], reads with = 1, = 6 or = 2, = 30, and = 0 for weak and > 0 for strong anisotropy. Also these models formally converge asymptotically to the expected anisotropic surface diffusion models as → 0 and have been frequently used in applicationssee, e.g., [17,18,34]. To be practically useful also Eqs. (15) -(17) require an appropriate regularisation 11 . In the following, we will refer to Eqs. (15) -(17) as RRV model.

NUMERICS
Our goal is to demonstrate the advantageous solution properties of the DDCH models, as well as to show that the variational DDCH models are suitable for simulating the formation of complex faceted morphologies. We will in the following consider the case = 1 and = 6 for the variational model (12) - (14). The classical DCH model can be straightforwardly considered by setting = 0 and = 1. (Thus, we can take ≡ 1). For the time-integration of the RRV model we refer to [11,13]. We employ the regularised functionŝ , , and using the value = 10 −6 . We solve the models using adaptive finite elements, implemented in the software package AMDiS 47,48 . The time-integration scheme is semi-implicit and follows the scheme employed for the isotropic case 1 , which reads, account for the linearizations of ( +1 ) ′ ( +1 ) and ′ ( +1 ) ( +1 ) around . The integer is the time step index, the time stepsize at step , and 0 the initial condition.
To employ this scheme for the anisotropic cases, (̂ ) has to be included. The scheme for the variational model (12) -(14) reads where For weak anisotropies we set = 0 and Eq. (22) is not needed. The system is discretized in space by piecewise linear finite elements. We consider simplicial meshes which are adaptively refined by bisection to ensure a resolution ranging between 5ℎ ∼ and 10ℎ ∼ , with mesh size ℎ, within the diffuse interface. The resulting linear system is solved by an iterative solver exploiting the biconjugate gradient stabilized method (BiCGstab(l)) 49 , with a block Jacobi preconditioner (bjacobi) applied to the blocks resulting from an element-wise domain decomposition, with a local sparse direct solver UMFPACK 50 .

NUMERICAL RESULTS
With this choice, (̂ ) is convex (weakly anisotropy) for values < 1 15 and non-convex (strongly anisotropy) for > 1 15 . The evolution obtained with Eqs. (9) -(10) or Eqs. (13) - (14), depending on , is illustrated in Fig. 1 . The initial setting is a 2D shape with varying positive and negative curvatures as in [1]. The final shape, the obtained equilibrium configuration, is the approximation of the corresponding Wulff shape, with facets and rounded corners for the non-convex (strongly anisotropic) free energy. For both cases the evolution can also be obtained with the RRV and DCH models. (Recall that the DCH model is obtained from the DDCH model using ≡ 1.) Visually there is no difference. However, as already discussed for the isotropic case, the variational DDCH model has better "positivity" properties. Note, however, that since we are using a regularised model, and since our numerical method is not designed for positivity, the numerical approximations of the DDCH model do not preserve the positivity property 0 < < 1 precisely. However the overshoot values are very small for the DDCH and RRV models. The DDCH model performs the best, as the overshoots are closer to the pure phase values = 0 and = 1.
Clear evidence of this aspect in the anisotropic case is reported in Fig. 2 . It shows the deviation from the nominal values of , namely 0 and 1, away from the interface through −min( ) and max( ) − 1, respectively. The variational DDCH model shows values, which are between one and two orders of magnitudes closer to the nominal values (0 in both plots) with respect to the RRV model and three to four orders of magnitude closer with respect to the DCH model. It is worth mentioning that these values depend on the specific simulation set up. However, we expect the relative difference between the models to hold true generally.  The better approximation properties of the variational DDCH and RRV models, if compared with the classical DCH model in the isotropic case, have already been shown in [1,24], respectively. These properties result in the possibility to use larger and thus lower grid resolution, to reach the same accuracy, which enables the large-scale simulations in various applications. In the anisotropic setting the differences between the variational DDCH and RRV models and the classical DCH model become even more severe and might lead to divergence of the DCH model for strong anisotropies and desired error tolerances. To demonstrate this we compare the reached equilibrium shape with the desired Wulff shape of the sharp interface problem. We closely follow convergence studies considered in [5] for a sharp interface model for strongly anisotropic surface diffusion. Fig. 3 (a) -(c) show the evolution from a circle with radius = 5 towards the equilibrium shape for different . We consider = 0.6 and (̂ ) as in Eq. (24) with = 1 2 . Focusing on the rounded upper corner of the desired Wulff shape allows to construct an asymptotic Comparison of ( ) close to the corner obtained by the variational DDCH, RRV and DCH models for two spatial discretizations for = 1.0 (f) and = 0.25 (g). The legend in (f) also holds for (g).
solution for the limiting sharp interface problem, see [51] for details. Briefly, the idea is to take the sharp corner equilibrium shape ( = 0) as the outer solution and to derive an inner solution for the equilibrium shape near the corner as an expansion in 1∕2 . Let be the arclength with = 0 at the corner of the outer solution and Θ( ) = ( ∕ 1∕2 ) the rescaled local orientation.
Expanding Θ( ) = Θ 0 ( ) + 1∕2 Θ 1 ( ) + ⋯, one obtains the lowest order term Θ 0 ( ) by inverting with = − ( ) cos( ) + ′ ( ) sin( ) and = ± 4 the corner orientation of the outer solution. The composite solution in the neighborhood of a corner with inner solution inner is then given by ( ) = inner ( ) + outer ( ) − match , with match = ± 4 , with the sign depending on the sign of . Fig. 3 (d),(e) show the solution for various in an ( , )-and ( , )-plot, respectively. The comparison of our numerical solutions of the variational DDCH, RRV and DCH models for fixed with these asymptotic solutions is shown in Fig. 3 (f),(g), for = 1 and = 0.25, respectively. Two spatial discretizations are shown corresponding to 5 (ℎ = 0.12 and 10 (ℎ = 0.06) grid points across the diffuse interface, demonstrating the good approximation properties of the double degenerate models already for moderate values of .
As a last example we demonstrate the possibility to simulate coarsening of bulk nanoporous materials, as e.g. considered with (̂ ) = 1 in [30,35], in a strongly anisotropic setting. We highlight a design project for the Digital Archive of Mathematical Models (DAMM) 55 . Within computational domains, defined by the Wulff shapes, a classical Cahn-Hilliard model with random perturbations of = 1 2 as initial condition, is solved. After spinodal decomposition the solution is used as initial condition for the anisotropic variational DDCH model, which is solved for several time steps similarly to Fig. 5 . The isosurface = 1 2 is extracted, postprocessed and 3D printed. Fig. 6 shows three examples of the resulting objects, a sphere (isotropic) and a square and a tetrahedra with rounded corners and edges.

CONCLUSIONS
In this paper we have extended a recently proposed variational DDCH model for isotropic surface diffusion 1 to the anisotropic case. We consider weak and strong anisotropies. The first case only requires the energy with a singular restriction function to be multiplied by a surface energy (̂ ), as in [13]. The second case requires an additional regularisation by a Willmore functional, again following [13]. A direct connection of the resulting models with the well known, and well used, non-variational RRV model, [11], as possible in the isotropic case, cannot be achieved for anisotropic surface energies. However, numerical comparisons show the similarities of both approaches in terms of approximation properties to the sharp interface limit. We omit the asymptotic analysis of the anisotropic variation DDCH model. Formal convergence results to anisotropic surface diffusion can be obtained by combining the analysis in [1,11,13].
The slightly more complex variational DDCH models, if compared with the non-variational RRV model, [11], not only gives energy dissipation and a mathematical foundation for numerical analysis, it also provides better positivity preserving properties. Besides the superior approximations properties, this positivity preserving property enables the wide applicability of the model in [11] for large-scale applications. Our numerical results indicate that this property is further improved in the variational DDCH models. This is beneficial for several reasons. For a numerical point of view it allows to optimize criteria based on thresholds or changes of as, e.g., for refining spatial discretization or adaptive timestepping. Also, these models are often coupled with additional equations in the bulk phases, which are characterized by the limiting values of , e.g., elasticity 11,32 or compositions 36 . In these cases, any improvement on positivity preservation, leads to increased accuracy.
However, besides these advantages, which will further foster to application of these models in materials science, as for the isotropic setting, several questions related to existence, uniqueness, and regularity of solutions, as well as to the positivity preserving property, remain open and need to be addressed. Dresden (TUD). SMW acknowledges generous financial support from the US National Science Foundation through grant NSF-DMS 1719854.