## 1 Introduction

Motion by surface diffusion, where the normal velocity of a hypersurface in Euclidean space is given by the surface Laplacian of the mean curvature, plays an important role in various applications in material science. An important example is solid state dewetting. See, for example, [Tho12] for a general introduction of this subject. While various direct numerical approaches have been developed for this fourth order geometric evolution law – see, e.g. [BMN05, HV05, BJWZ17, BGN19] – for many applications, diffuse interface approximations are considered [WLK05, RRV06, YCG06, TLVW09, LLRV09, BN09, JBTS12, SBB15, BSB16, NBS17, SKSV17]. 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 fourth-order 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 [CENC96, RRV06, GSK08, DMW17].

The now-conventional wisdom is that, when the Cahn-Hilliard equation is paired with a polynomial free energy, for the aforementioned sharp interface convergence result to hold, it is required that the degeneracy in the so-called mobility function needs to be of sufficiently high order [Voi16]. As recently shown [DD14a, LMS15, LMS16], occasionally used second order degenerate mobility functions – see e.g., [BKB00, WKL06, TLVW09, JBTS12] – do actually not converge to surface diffusion, but contain additional bulk diffusion contributions. This is not the case for the originally proposed phase field approximation of surface diffusion [CENC96], which uses a double-obstacle potential (in the deep-quench limit) instead of the double-well polynomial potential. In a series of papers [DD12, DD14b, DD16b, DD16a], the authors conduct careful studies, both computational and theoretical, examining the effects on coarsening rates, asymptotic limits, and other solution properties, when using a range of diffusional mobility types in the Cahn-Hilliard equation, including one-sided and two-sided degeneracies of various orders. Though these papers did not solely target Cahn-Hilliard models for surface diffusion, they do confirm the findings in [DD14a, LMS15, LMS16].

In the diffuse interface model for surface diffusion proposed in [RRV06], an additional degeneracy is introduced, following similar ideas as used for the thin film limit in classical phase field models for solidification [KR98]. We refer to this model as a doubly degenerate Cahn-Hilliard (DDCH) equation. This second degeneracy does not alter the asymptotic limit [RRV06], but actually leads to more accurate surface diffusion approximations, see, e.g., [BWSV19]. In fact several simulations for realistic applications in materials science consider this additional degeneracy [ABM16, SBI15, SBV16, SBB17, SBVM17, NBS17, GBWK19, ABS19]. However, the model of [RRV06] also has a drawback. 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. Furthermore, it excludes variational derivations in more complex multi-physics applications involving surface diffusion.

This paper is organized as follows. In Section 2, we reintroduce the DDCH model of [RRV06], generalizing it in a trivial way and recalling some asymptotic convergence results. In Section 3, we introduce a new DDCH variational model for surface diffusion, and we analyze some of its properties. We make connections in Section 3 between the new DDCH variational model and (i) the non-variational DDCH model of [RRV06] and (ii) a deGennes diffuse interface model from another modeling context. We provide a simple, but powerful, numerical integration scheme in Section 4, and we compare the numerical results of the models in Section 5

with various parameter choices. The matched asymptotic analysis for the (new) variational and the (well-used) non-variational surface diffusion models are provided in Appendices

A and B, respectively.## 2 A Non-Variational Diffuse Interface Model

In this section we reintroduce a non-variational DDCH model that approximates sharp interface motion by surface diffusion. Suppose that is a bounded subset of . Let be an order/phase field parameter. Consider the following Cahn-Hilliard-type model proposed in [RRV06]:

(1) | ||||

(2) |

where is a small parameter (relative to the domain size) that is related to the thickness of the diffuse interface. denotes the mobility function, which is defined as

and is the quartic, symmetric double well potential function

Its minima, 0 and 1, represent the pure phase states. Let us assume that, for simplicity, , the phase field, and , called the chemical potential, satisfy periodic boundary conditions with respect to . Clearly, mass is conserved in the system: .

To keep the model as general as possible, let us assume that – which we shall call the *diffusion restriction function* – is defined as

The authors of [RRV06] assumed that and , and they showed using matched asymptotic analysis that, as , solutions converge formally to those of a sharp interface model of surface diffusion. Let define a hypersurface that coincides with the 0.5 level set of at some time . We say that evolves by surface diffusion if

(3) |

where is the normal velocity at the sharp interface , is the surface Laplacian, and is the mean curvature of . We show in Appendix B, by a simple calculation, that for the generalized model, upon taking

(4) |

where is the usual (Bernoulli) Gamma function, the limiting law as is again motion by surface diffusion (3). The normalization in (4) is called the *diffusion normalization*. (Note that , which recovers the result in [RRV06].)

The model (1) – (2) is a type of doubly degenerate Cahn-Hilliard (DDCH) equation, because degeneracies are associated with both the mobility and the restriction functions. We refer to this model (1) – (2) as the non-variational DDCH model (or, Model NV, for short), because it is not arrived at through a variational principal. Indeed, it is not clear whether there is a related energy that is dissipated along the solution trajectories of (1) – (2). In other words, it seems that the model is not free energy dissipative, or thermodynamically consistent. Model NV becomes degenerate as approaches the pure state values and . Indeed, both the mobility and the restriction function are degenerate when . There is ample numerical evidence to suggest that there is an important benefit owing to these degeneracies, namely, solutions to Model NV remain in the physically relevant region , as long as the initial data have this feature. This, to our knowledge, has yet to be verified theoretically.

Incidentally, this property – , for all – is often referred to as a *positivity preserving* property. To see why, note that for a typical binary model, like the traditional Cahn-Hilliard equation, the concentrations of the species and are given by and . Therefore we have the positivity of the concentrations, and , iff . We use this terminology herein.

We point out that one can regularize Model NV so that it is defined for all values of , such that the asymptotic limit should remain unchanged. This may be important for numerical implementation, since most numerical schemes cannot guarantee that solution remain positive, even when this property is known to hold for the PDE. For example, one can reintroduce the mobility as

To make the restriction function defined, continuous, positive, and differentiable on all of , one can regularize as follows:

Thus, setting we obtain the model above.

Upon choosing , , and , one obtains a more familiar (singly) degenerate Cahn-Hilliard (DCH) equation:

(5) | ||||

(6) |

Matched asymptotic analysis shows that solutions of this model also formally converges to those of a sharp interface model of surface diffusion, as . However, quite importantly, numerical results suggest such diffuse interface solutions converge more slowly to the sharp interface model compared with those of Model NV, (1) – (2). See, for example, [BWSV19]. The computational results in Section 5 support this point further. However, model (5) – (6) is free energy dissipative, that is, thermodynamically consistent. Formally, solutions to this model dissipate the free energy

at the rate

## 3 A Variational Diffuse Interface Model

Is there an energy dissipative DDCH model that is related to Model NV? The answer is a qualified *yes*. In this section, we show that there is a variational model such that Model NV is its approximation, in certain cases. To this end, consider the free energy

(7) |

where is the same quartic potential as before. Here we assume that is a singular function of the form

We call the *energy restriction* function. If necessary, can be regularized so that it is defined, continuous, and differentiable for all :

Thus, setting , we obtain the function above.

We are interested in the energy dissipative flow

(8) | ||||

(9) |

where is the chemical potential and is the variational derivative with respect to . Here, for simplicity, we have assumed natural or periodic boundary conditions on . System (8) – (9) is another type of doubly degenerate Cahn-Hilliard (DDCH) equation, which we refer to as the *variational DDCH model* (or, Model V, for short).

The derivative is, of course, singular,

which, it seems, help to keep solutions positive. The regularized version is non-singular,

though its values at 0 and 1 become increasing large as . There are some open analysis questions related to the validity of the model when the regularization vanishes (). For example, do solutions have the positivity preserving property? Do weak solutions exist and are they regular, in some sense? For what values of and does the model make sense? Our early simulation results – as well as some early theoretical results, not reported here – seem to support the validity of a positivity preserving property, when . In any case, formally, it is clear that the energy given in (7) is dissipated along the solution trajectories of the dynamical system (8) – (9), that is, system (8) – (9) is a free energy dissipative dynamical system.

### 3.1 Relation to Model NV

Let us now see how the variational/dissipative DDCH model (Model V) may be related to the non-variational DDCH model (Model NV). Assume that, to leading order in , the profile of the phase field function perpendicular to the diffuse interface is the hyperbolic tangent function

(10) |

where is the coordinate perpendicular to the interface, a type of interface distance function. Then, we have the asymptotic approximation, to leading order in ,

(11) |

(We will justify these approximations in the appendices, though they are somewhat standard.) Since

we can approximate the chemical potential as

In other words,

Thus, Model V is related to Model NV through approximation, if the interfacial profile is a hyperbolic tangent, and if . In the next section, we find that there is a logical choice for , given , for Model V, the so-called energy normalization.

### 3.2 Finite Energy and Energy Normalization

When, if ever, is the energy (7) of a diffuse interface solution with the hyperbolic tangent profile finite? To answer this, let us work in two space dimensions, for simplicity. If (10) is valid to leading order, then, using the approximation (11) the energy can be approximated as

where is the length of the diffuse interface, measured in the tangential direction along the level curve. (In three dimensions, the surface area of the iso-surface would appear.) For , taking , we find that the energy is finite and

This is the usual measure of energy for isotropic motion by mean curvature in two dimensions, namely, the energy is proportional to the interface length. In fact, we can generalize the last result: for any ,

We can use this calculation to pick the value of in the general case so that , for all . For Model V, let us choose

(12) |

We refer to this choice as the *energy normalization*. This is different from the diffusion normalization (4), though the two values are equal when and :

In any case, observe that,

Consequently, the energy cannot be made finite for the hyperbolic tangent profile when . This seems to imply that Model V cannot be made sensible for .

For , using the energy normalization, we show in Appendix A that solutions to Model V, (8) – (9), converge, as , to the solutions of the sharp interface model of surface diffusion (3). This fact leads us to the following conclusion: for , choosing the energy normalization in Model V and the diffusion normalization in Model NV, the two models are directly connected via the approximation outlined in the previous section. We will demonstrate this connection in the computational results of Section 4. For now, the case has other interesting properties that we will explore in the next section.

### 3.3 Relation to a deGennes Model

Suppose in Model V that we select the parameters

Under the assumption that solutions are restricted to the physically meaningful region of phase space , the energy (7) is formally equivalent to

(13) |

where is the obstacle potential

The energy density is called the *deGennes gradient energy density* [DWZZ19, LQZ17]. Here, it seems, the obstacle potential does not play an essential role, since the singular deGennes gradient energy density is expected to keep the solution within the physically meaningful range in a gradient flow setting. However, this version of the energy, , is instructive, as it can be viewed as the deep-quench limit of the following Flory-Huggins-deGennes-type model:

(14) |

In other words, is obtained in the (zero temperature) limit of . We note that is a model that is commonly used in the study of polymer or hydrogel materials [DWZZ19, LQZ17]. Furthermore, is also reminiscent of the standard Flory-Huggins free energy considered in [CENC96], before the deep-quench assumption is invoked.

One interesting feature of the deep-quench energy, , is that hyperbolic tangent solutions are one-dimensional minimizers. This implies that, for the associated conserved and dissipative dynamical system, evolving diffuse interfaces will have the hyperbolic tangent shape. This is in contrast to the deep-quench model considered in [CENC96], where diffuse interfaces have cosine-like profiles.

In future works, we will address some theoretical issues related to Model V, for example, the existence of weak solutions [DD16b]

, the validity of the positivity-preserving properties of the PDE and associated numerical schemes, et cetera. For the moment, we will explore some numerical solutions of Models V and NV, comparing their characteristics and properties, under the assumptions that the PDE solutions are sufficiently well-behaved.

## 4 Integration Schemes

To approximate solutions of Models NV and V, we employ the AMDiS finite element software package [VV07, WLPV15], which allows for adaptive mesh refinement around the evolving diffuse interface. To remove some of the stiffness, a linear implicit-explicit (IMEX) integration scheme is used to approximate solutions of Eqs. (8) – (9):

(15) |

where

account for the linearizations of the and terms around . The integer is the time step index; is the time stepsize at step ; and is the initial condition. We use a very small regularization constant, , as the current numerical scheme is not designed to preserve the expected positivity of the solution. Positivity-preserving schemes are being developed and will be reported on in later works. is an auxiliary variable, such that selects Model V, while enforces the approximation and allows for the integration of Model NV. Indeed, with the system can be rewritten as

(16) |

where , with and . We thus obtain a semi-implicit integration scheme for Eqs. (5) – (6). In the following, we use the system (15) with to integrate Model V and the system (16) to integrate Model NV. Notice that this distinction strictly holds true for . Adaptive mesh refinement is exploited with a fine spatial discretization within the diffuse interface, namely where , and a coarse spatial discretization elsewhere. The former is scaled with to ensure the same number of elements within the diffuse interface for any , while the latter is fixed. The model and numerical parameters are set as follows: , , , .

The results of the numerical integration of the DDCH equations are compared with the dynamics of the corresponding sharp-interface (SI) evolution in two space dimensions. In particular, we will consider the SI evolution of the closed curve, , corresponding to the level-set of the initial condition data for the DDCH models. This curve is discretized by a uniform set of points . The discrete evolution law in terms of the velocities directed along the outward-pointing unit normals read

(17) |

where is a finite difference approximation of the surface Laplacian (in this reduced dimension setting) and is the local curvature, computed as with the radius of curvature. The latter is approximated as the radius of the almost osculating circle, i.e., the circle passing through and . We note that several other methods have been proposed and used instead of the approximation in Eq. (17) (see, e.g, [BMN05, HV05, HV07, BGN08, WJBS15, BJWZ17]). Here, according to the simple evolution it has been sufficient to use a simple finite-difference scheme exploiting a forward Euler discretization of time, with a sufficiently small time stepsize.

## 5 Numerical Results

We numerically compare Model V, Model NV and the standard DCH model () with the sharp-interface (SI) solution for surface diffusion.

Figure 1 illustrates the evolution of an ellipse obtained by integrating Model V using (15). The initial condition is set by evaluating Eq. (10) with the signed distance to the considered ellipse. Figure 2 illustrates the evolution of the level-set by the SI approach, i.e. by Eq. (17). A few representative profiles during the evolution are shown in Figure 2(a), while the change over time of the x semi-axis is shown in Figure 2(b).

Figure 3 shows the behavior of the DDCH models for different values of and different . We recall that a linear scaling of the timestep and of the mesh size with have been considered. The approximation properties of the DDCH models to the sharp interface limit of surface diffusion are summarized. While all DDCH models converge to the SI solution, the error is lower for and compared to the standard DCH model (). For both Model V and Model NV can be considered. Convergence to the SI limit is similar in these cases, but we obtain slightly lower errors with Model NV (see Figure 3(d)). This may be ascribed to the larger number of operators to be considered for the integration of Model V. We recall that for Model V coincides with Model NV under the assumption .

The convergence rate is considered in Figure 4 using the relative deviation from the SI solution at time , as

(18) |

This quantity is shown for an early stage of the simulation focusing on the fast dynamics. For all DDCH models , and we obtain linear convergence. However, with a significantly smaller error for and . This confirms previous results [BWSV19] and extends them to Model V. In any case, using the DDCH models allows to consider -values to be at least doubled if compared with the classical DCH model, to reach the same accuracy.

The behavior of the DDCH models discussed so far holds true also when considering more complex surface profiles. The case of the evolution by surface diffusion of a profile where also the sign of the local curvature changes is illustrated in Figure 5.^{1}^{1}1The profile of Figure 5 is known as Camunian Rose, a stylized version of iron-age rock carvings found in northern Italy [Far97]. Figure 5(a) at shows the considered profile obtained with , along with the computational mesh adopted for the numerical simulations. In the same panel, four representative stages during the evolution are also shown. The shapes in Figure 5(a) correspond to the region with a color map highlighting the extension of the interface region. Figure 5(b) shows the different surface profiles, namely the level curve, obtained with different values of and different choice of for the Model NV. The increase of accuracy with increasing is qualitatively confirmed here. Figure 5(c) shows that Model V () delivers very similar results to Model NV as expected from the convergence results illustrated in Figure 4.

## 6 Discussion and Conclusions

In this paper we have introduced a new variational model, Model V, describing surface diffusion. This model has an energy, though this energy is defined via a singular restriction function that must be chosen carefully. One of the main purposes for the creation of the new model, was to connect it, via certain approximations, to a well known, and well used, non-variational model, Model NV. Both Models V and NV are shown to converge to sharp interface surface diffusion, via the formal method of matched asymptotics. This convergence is further supported by our numerical results, which suggest that the convergence rates with respect to are all order 1, that is, the convergence rates are all linear in .

Both models are doubly degenerate Cahn-Hilliard (DDCH) equations, and there are relative strengths and weaknesses associated to each. The new model, Model V, is more complicated to solve numerically, but has an associated energy, which makes its connection to other physics though additional energy terms, more seamless. Like for Model NV, the solutions to Model V, at a fixed value and value, more accurately approximate solutions of the sharp interface model, compared to the singly degenerate Cahn-Hillard (DCH) model (where, the only degeneracy is associated to the mobility). Solutions to Model NV, it seems, are slightly more accurate approximations than those of Model V for the case; and solutions obtained from Model NV with are more accurate than any of the solutions computed from Model V.

Several open questions related to both models remain, questions related to existence, uniqueness, and regularity of solutions, as well as to the positivity preserving property. In our computations – where we used the regularized versions of the equations, since our schemes could not maintain positivity – solutions to Model V were closer to being positive than solutions to Model NV for the case (data not shown). In other words, the overshoots of the values outside of the physically realistic range, , were smaller in a point-wise sense for Model V. Indeed, Model V has a singular energetic mechanism for positivity that Model NV may not possess. However, either way of including a restriction function seems to help in maintaining positivity, when compared to the standard singly degenerate Cahn-Hilliard model for surface diffusion. Our numerical experiments suggest that solutions to both models will remain positive, in the sense that reducing the regularization reduces the point-wise overshoot of the values outside of the physically realistic range, . In future works, we plan to report on energy stable and positivity preserving numerical schemes for Model V, as well as some other theoretical issues, including the existence of weak solutions.

## Acknowledgments

This research was partially funded by the EU H2020 FET-OPEN project microSPIRE (ID: 766955) and by the EU H2020 FET-OPEN project NARCISO (ID: 828890). We gratefully acknowledge the computing time granted by the John von Neumann Institute for Computing (NIC) and provided on the supercomputer JURECA at Jülich Supercomputing Centre (JSC), within the Project no. HDR06, and by the Information Services and High Performance Computing (ZIH) at the Technische Universität Dresden (TUD). SMW acknowledges generous financial support from the US National Science Foundation through grant NSF-DMS 1719854 and acknowledges some enlightening discussions with Prof. Shibin Dai on the topics of degenerate mobilities in the Cahn-Hilliard equation and with Prof. Cheng Wang on the topic of energy stable and positivity preserving methods for Cahn-Hilliard-deGennes-type models. Some of this work was performed while the third author (SMW) was supported as a Dresden Senior Fellow at TUD for three months in 2017 in the group of the second author (AV). SMW gratefully acknowledges this funding and the hospitality of AV and TUD during his visit.

## Appendix A Matched Asymptotic Analysis of Model V

In this appendix, we will prove, using the method of matched asymptotic expansions, that solutions of the variational DDCH system (8) – (9), which we refer to as Model V, converge to the solutions of the sharp interface model of surface diffusion (3). To keep the discussion simple, we consider only the case for which space is two dimensional. In this setting, the surface diffusion of a one-dimensional evolving curve is modelled by the equation

where is the arc length of , and is the curvature. Herein we assume that in the definition of the restriction function and we set the regularization to zero (). We closely follow the analyses in [CENC96, RRV06].

### a.1 Interface Coordinates and Asymptotic Expansions

We suppose that is a bounded open set, for example, a rectangular domain. Let us assume that there is a simply-connected domain that has a smooth boundary . Assume that is a smooth function with the property that

The idea is that, as ,

converges to the characteristic function

. We assume that is the solution to Model V, that is, (8) – (9) with initial data andClearly . We assume that remains smooth and there is always a neighborhood such that there is a well-defined signed distance function . There is a well-definite interior, where and corresponding to , and an exterior, where and corresponding to .

We assume for simplicity that is closed, and its length, denoted , is finite. Thus can be parameterized by its arc length,

. The parameterization is denoted by the vector function

which is an -periodic, twice continuously differentiable function that is a bijection when its domain is restricted to any half-open interval of length .

On the interface , there are well-defined unit tangent and unit normal vector functions, which we denote by and . The unit tangent vector points in the direction of increasing arc length, while the unit normal points in the direction of increasing interface distance . In fact, in a small enough neighborhood

and are local coordinates near the interface . The set is called the narrow band neighborhood of the interface, or just the narrow band, for short. In particular, for every , there are points and such that

In other words, for every , there is a unique pair

We can now define transformed dependent variables in the narrow band via composition: for all , define

To facilitate our asymptotic analysis, we introduce the stretched interface distance variable

Subject to this transformation, we define, for all , the stretched variables

We assume that the following Taylor expansions are valid

(19) | ||||

(20) | ||||

(21) | ||||

(22) | ||||

(23) | ||||

(24) |

Expansions (19), (20), (22), and (23) are called outer expansions; while expansions (21) and (24) are called inner expansions. The following matching conditions are employed to join the two types of expansions together:

(25) | |||

(26) |

We denote by the normal velocity of the and, by , the curvature of . Then, we have the following well-known relations in two dimensions:

As a consequence, using the chain rule,

Since we are parameterizing by arc-length, it follows that

and, therefore,

Equivalently,

In a similar way, we can derive the identity

Using these and related identities, we can express the standard operators in the transformed coordinate system :

where is any vector function that transforms as and . The full derivations of the equations above can be found in [RRV06]. Combining the expansions above, we can express the Laplacian operator as

Likewise,

### a.2 Outer Expansion

The evolution equations become singular in the outer region. While the energy for the pure states can be defined in reasonable way, in particular,

the gradient equations are undefined for the pure states. This is analogous to the deep-quench limit examined in [CENC96]. An expansion of solutions may not be available in the usual sense. Without more information, we will assume the following reasonable form of the solution in the outer regions:

where is a constant. This is enough for our analysis to go through, though we leave open the possibility that , as was assumed in [CENC96] for the deep-quench limit case. The structure of the outer solutions corresponds to matching conditions for the inner expansions of the form

### a.3 Inner Expansion

#### a.3.1 Equation (9) at

We will assume that in the inner region, the solution stays away from the pure states:

Let us make an expansion for Eq. (9), which is rewritten slightly as

where for simplicity, we have dropped the subscript on . Subscripts will have a different meaning in this section, as we will see. We have

(27) |

Here

and

Comments

There are no comments yet.