Doubly Degenerate Diffuse Interface Models of Surface Diffusion

We discuss two doubly degenerate Cahn-Hilliard (DDCH) models for surface diffusion. Degeneracy is introduced in both the mobility function and a restriction function associated to the chemical potential. Our computational results suggest that the restriction functions yield more accurate approximations of surface diffusion. We consider a slight generalization of a model that has appeared before, which is non-variational, meaning there is no clear energy that is dissipated along the solution trajectories. We also introduce a new variational and, more precisely, energy dissipative model, which can be related to the generalized non-variational model. For both models we use formal matched asymptotics to show the convergence to the sharp interface limit of surface diffusion.


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 [WLK + 05, RRV06, YCG06, TLVW09, LLRV09, BN09, JBTS12, SBB + 15, BSB + 16, NBS + 17, 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 doubleobstacle 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, SBI + 15, SBV16, SBB + 17, SBVM17, NBS + 17, GBWK19, ABS + 19]. 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.

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 R d . Let u : Ω → R be an order/phase field parameter. Consider the following Cahn-Hilliard-type model proposed in [RRV06]: where ε > 0 is a small parameter (relative to the domain size) that is related to the thickness of the diffuse interface. M 0 denotes the mobility function, which is defined as and f is the quartic, symmetric double well potential function Its minima, 0 and 1, represent the pure phase states. Let us assume that, for simplicity, u, the phase field, and w, 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 G 0 -which we shall call the diffusion restriction function -is defined as The authors of [RRV06] assumed that p = 2 and η = 30, and they showed using matched asymptotic analysis that, as ε 0, 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 u at some time t. We say that Σ evolves by surface diffusion if where v 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 where Γ is the usual (Bernoulli) Gamma function, the limiting law as ε 0 is again motion by surface diffusion (3). The normalization in (4) is called the diffusion normalization. (Note that η (p = 2) = 30, 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 u approaches the pure state values u = 0 and u = 1. Indeed, both the mobility M 0 and the restriction function G 0 are degenerate when u = 0, 1. 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 0 < u < 1, as long as the initial data have this feature. This, to our knowledge, has yet to be verified theoretically.
Incidentally, this property -0 < u( · , 0) < 1 =⇒ 0 < u( · , t) < 1, for all t ≥ 0 -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 A and B are given by u A = u and u B = 1 − u. Therefore we have the positivity of the concentrations, u A > 0 and u B > 0, iff 0 < u < 1. We use this terminology herein.
We point out that one can regularize Model NV so that it is defined for all values of u, 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 R, one can regularize as follows: Thus, setting α = 0 we obtain the model above. Upon choosing p = 0, τ = 1, and α = 0, one obtains a more familiar (singly) degenerate Cahn-Hilliard (DCH) equation: Matched asymptotic analysis shows that solutions of this model also formally converges to those of a sharp interface model of surface diffusion, as ε 0. 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

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 where f is the same quartic potential as before. Here we assume that g 0 is a singular function of the form We call g 0 the energy restriction function. If necessary, g 0 can be regularized so that it is defined, continuous, and differentiable for all u: Thus, setting α = 0, we obtain the function above. We are interested in the energy dissipative flow where w = δ u F is the chemical potential and δ u F is the variational derivative with respect to u.
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 g 0 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 α 0.
There are some open analysis questions related to the validity of the model when the regularization vanishes (α = 0). For example, do solutions have the positivity preserving property? Do weak solutions exist and are they regular, in some sense? For what values of p 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 0 < p < 2. In any case, formally, it is clear that the energy F 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.

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 u perpendicular to the diffuse interface is the hyperbolic tangent function where z is the coordinate perpendicular to the interface, a type of interface distance function. Then, we have the asymptotic approximation, to leading order in ε, (We will justify these approximations in the appendices, though they are somewhat standard.) Since ∇ · [g 0 (u)∇u] = g 0 ∆u + g 0 (u)|∇u| 2 , 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 p, for Model V, the so-called energy normalization.

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 u = 0.5 level curve. (In three dimensions, the surface area of the 0.5 iso-surface would appear.) For p = 1, taking γ = 6, 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 p ∈ [0, 2), We can use this calculation to pick the value of γ in the general case so that F [u] ≈ |Σ|, for all p ∈ [0, 2). For Model V, let us choose We refer to this choice as the energy normalization. This is different from the diffusion normalization (4), though the two values are equal when p = 0 and p = 1: In any case, observe that, lim Consequently, the energy cannot be made finite for the hyperbolic tangent profile when p ≥ 2. This seems to imply that Model V cannot be made sensible for p ≥ 2.
For 0 ≤ p < 2, using the energy normalization, we show in Appendix A that solutions to Model V, (8) -(9), converge, as ε 0, to the solutions of the sharp interface model of surface diffusion (3). This fact leads us to the following conclusion: for p = 1, 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 p = 1 has other interesting properties that we will explore in the next section.
Under the assumption that solutions are restricted to the physically meaningful region of phase space 0 ≤ u ≤ 1, the energy (7) is formally equivalent to where I is the obstacle potential The energy density ε 2 g 0 (u)|∇u| 2 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, F 0 , is instructive, as it can be viewed as the deep-quench limit of the following Flory-Huggins-deGennes-type model: In other words, F 0 is obtained in the (zero temperature) limit θ 0 of F θ . We note that F θ is a model that is commonly used in the study of polymer or hydrogel materials [DWZZ19,LQZ17]. Furthermore, F θ 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, F 0 , 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.

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): where account for the linearizations of the g α (u n+1 )f (u n+1 ) and g α (u n+1 )f (u n+1 ) terms around u n . The integer n is the time step index; τ n > 0 is the time stepsize at step n; and u 0 is the initial condition. We use a very small regularization constant, α > 0, 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 χ = 1 selects Model V, while χ = 0 enforces the approximation f (u)/ε ≈ (ε/2)|∇u| 2 and allows for the integration of Model NV. Indeed, with χ = 0 the system can be rewritten as where G α (u n ) = 1/g α (u n ), with p = 1 and γ = η = 6. We thus obtain a semi-implicit integration scheme for Eqs. (5) -(6). In the following, we use the system (15) with χ = 1 to integrate Model V and the system (16) to integrate Model NV. Notice that this distinction strictly holds true for p > 0. Adaptive mesh refinement is exploited with a fine spatial discretization h in within the diffuse interface, namely where 0.05 < u < 0.095, and a coarse spatial discretization h out 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: α = 10 −4 , τ n = τ = 10 −5 ε, h in = ε/10, h out = 0.125. 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 0.5 level-set of the initial condition data for the DDCH models. This curve is discretized by a uniform set of points {x i }. The discrete evolution law in terms of the velocities directed along the outward-pointing unit normalsn read where∆ Σ is a finite difference approximation of the surface Laplacian (in this reduced dimension setting) and κ i is the local curvature, computed as 1/r i with r i the radius of curvature. The latter is approximated as the radius of the almost osculating circle, i.e., the circle passing through x i and x i±1 . 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.

Numerical Results
We numerically compare Model V, Model NV and the standard DCH model (p = 0) with the sharp-interface (SI) solution for surface diffusion.  (15)) of an ellipse with semi-axis along x axis (A x = 1.0) and along the y axis (A y = 0.5). Illustrated by means of u (top) and w (bottom), the latter shown in the region where 0.1 < u < 0.9. ε = 0.2, p = 1. Figure 1 illustrates the evolution of an ellipse obtained by integrating Model V using (15). The initial condition u 0 is set by evaluating Eq. (10) with r the signed distance to the considered ellipse. Figure 2 illustrates the evolution of the u = 0.5 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 A x (t) is shown in Figure 2 Figure 3 shows the behavior of the DDCH models for different values of p 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 p = 1 and p = 2 compared to the standard DCH model (p = 0). For p = 1 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 p = 1 Model V coincides with Model NV under the assumption (1/ε)f (u) ≈ (ε/2)|∇u| 2 .
The convergence rate is considered in Figure 4 using the relative deviation from the SI solution at timet, as This quantity is shown for an early stage of the simulation focusing on the fast dynamics. For all DDCH models p = 0, p = 1 and p = 2 we obtain linear convergence. However, with a significantly smaller error for p = 1 and p = 2. 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 Figure 5(a) at t = 0 shows the considered profile obtained with ε = 0.2, 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 u > 0.5 with a color map highlighting the extension of the interface region. Figure 5

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 p 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 p = 1 case; and solutions obtained from Model NV with p = 2 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 p = 1 case (data not shown). In other words, the overshoots of the values outside of the physically realistic range, 0 ≤ u ≤ 1, 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, 0 ≤ u ≤ 1. 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.

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 onedimensional evolving curve is modelled by the equation where s is the arc length of Σ, and κ is the curvature. Herein we assume that 0 ≤ p < 2 in the definition of the restriction function g α and we set the regularization to zero (α = 0). We closely follow the analyses in [CENC96,RRV06].

A.1 Interface Coordinates and Asymptotic Expansions
We suppose that Ω ⊂ R 2 is a bounded open set, for example, a rectangular domain. Let us assume that there is a simply-connected domain Ω 0 0 ⊂ Ω that has a smooth boundary Σ 0 := ∂Ω 0 . Assume that u 0,ε : Ω → [0, 1] is a smooth function with the property that The idea is that, as ε 0, u 0,ε converges to the characteristic function 1 Ω\Ω0 . We assume that u is the solution to Model V, that is, (8) -(9) with initial data u 0,ε and Clearly Σ 0 = Σ(0; ε). We assume that Σ(t, ε) remains smooth and there is always a neighborhood N (t; ε) ⊃ Σ(t; ε) such that there is a well-defined signed distance function r(x, y, t; ε). There is a well-definite interior, where r < 0 and corresponding to Ω 0 (t, ε), and an exterior, where r > 0 and corresponding to Ω \ Ω 0 (t, ε).
We assume for simplicity that Σ(t, ε) is closed, and its length, denoted L(t; ε), is finite. Thus Σ(t, ε) can be parameterized by its arc length, s ∈ [0, L(t; ε)]. The parameterization is denoted by the vector function which is an L(t; ε)-periodic, twice continuously differentiable function that is a bijection when its domain is restricted to any half-open interval [a, b) of length L(t; ε). On the interface Σ(t; ε), there are well-defined unit tangent and unit normal vector functions, which we denote by τ (s, t; ε) and n(s, t; ε). The unit tangent vector points in the direction of increasing arc length, while the unit normal points in the direction of increasing interface distance r. In fact, in a small enough neighborhood N ρ (t; ε) := {(x, y) ∈ Ω | |r(x, y, t; ε)| < ρ} ⊆ N (t; ε), s and r are local coordinates near the interface Σ(t; ε). The set N ρ (t; ε) is called the narrow band neighborhood of the interface, or just the narrow band, for short. In particular, for every x = (x, y) ∈ N ρ (t; ε), there are points r ∈ (−ρ, ρ) and s ∈ [0, L(t; ε)) such that (x, y) = x(r, s) = c(s, t; ε) + r(x, y, t; ε)n(s, t; ε).
To facilitate our asymptotic analysis, we introduce the stretched interface distance variable z = r ε .
We denote by v the normal velocity of the Σ(t, ε) and, by κ, the curvature of Σ(t, ε). Then, we have the following well-known relations in two dimensions: ∂ s τ = −κn and ∂ s n = κτ .
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 ∇ x u · n = ∂ rû .

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: u = 0 + εβ + O(ε 2 ) ≥ 0 and u = 1 − εβ + O(ε 2 ) ≤ 1, where β ≥ 0 is a constant. This is enough for our analysis to go through, though we leave open the possibility that β = 0, 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 This implies that M (U 0 )∂ z W 0 = C 2 (s, t).