Convergence of Space-Time Discrete Threshold Dynamics to Anisotropic Motion by Mean Curvature

# Convergence of Space-Time Discrete Threshold Dynamics to Anisotropic Motion by Mean Curvature

## Abstract

We analyze the continuum limit of a thresholding algorithm for motion by mean curvature of one dimensional interfaces in various space-time discrete regimes. The algorithm can be viewed as a time-splitting scheme for the Allen-Cahn equation which is a typical model for the motion of materials phase boundaries. Our results extend the existing statements which are applicable mostly in semi-discrete (continuous in space and discrete in time) settings. The motivations of this work are twofolds: to investigate the interaction between multiple small parameters in nonlinear singularly perturbed problems, and to understand the anisotropy in curvature for interfaces in spatially discrete environments. In the current work, the small parameters are the the spatial and temporal discretization step sizes and . We have identified the limiting description of the interfacial velocity in the (i) sub-critical (), (ii) critical (), and (iii) super-critical () regimes. The first case gives the classical isotropic motion by mean curvature, while the second produces intricate pinning and de-pinning phenomena and anisotropy in the velocity function of the interface. The last case produces no motion (complete pinning).

## 1 Introduction and Main Results

The current paper addresses convergence issues related to a thresholding scheme for motion by mean curvature. The key is the analysis of the algorithm in the space-time discrete setting in which there are two small parameters - the step sizes in the spatial and temporal directions. The ultimate results depend on the relative sizes of these parameters.

The analysis of motion by mean curvature (in which the normal velocity of a moving manifold is given by its mean curvature) is an active area. Not only it is interesting in geometry in its own right, it also finds many applications in materials science and image processing. It is a prototype of a gradient flow with respect to the area functional. Due to the possibility of singularity formation and topological changes of the evolving surface, elaborate approaches need to be used. These include (i) varifold formulation, (ii) the viscosity solutions, and (iii) singularly perturbed reaction diffusion equations.

The thresholding scheme is a particularly simple algorithm to capture the key feature of (iii). It is essentially a time splitting scheme. The first step is diffusion while the second step is thresholding to mimick the fast reaction due to the nonlinear term. This is heuristically proposed in [5] and rigorously proved in [9, 2] in the continuous space and discrete time setting. See also the work [7] for an analysis of an reaction diffusion equation in which both space and time variables are discrete. However, so far all the rigorous results essentially works in the case when the interfacial structure is well-resolved. We call this the “sub-critical” regime. When this is not the case, intricate pinning and depinning of the interface can happen. This is analogous to a gradient flow in a highly wiggling or oscillatory energy landscape. The motion also demonstrates anisotropy of the normal velocity. The motivation of the current paper is to capture these phenomena quantitatively and relate them to the underlying small parameters in the algorithm.

The most relevant reaction diffusion equation for motion by mean curvature is the following Allen-Cahn equation:

 ∂u∂t=△u−1ϵ2W′(u) (1)

In the above is the double well potential and is a small parameter. The qualitative behavior of the solution is that the underlying ambient space is quickly partitioned into two domains on which takes on the values and which are the minima of . The function also makes a smooth but rapid transition with thickness between the two domains. The key is then to understand the dynamics of this transition layer, in the limit of . It is proved in various settings that the limiting motion is motion by mean curvature [6, 8, 7, 11, 10].

As the thresholding scheme is very simple to implement and describe, we embark on its analysis demonstrating the interplay between two small parameters. The scheme is a time splitting approach to solve (1) (in the regime ). Given an initial shape , and its boundary (or often called the interface) , a sequence of functions is constructed in the following manner: for , we define

 Unknown environment '%

then the following two steps are alternately performed (for ),

• diffusion step:

 ∂v∂t=△vfor0
• thresholding step:

 uk+1(x)=sign(v(x,τ)) (3)

Note that the second step above is to mimick the fast reaction term which drives to or , the minima of . The solution of the problem is captured by the sequence of subsets where attains the value , . Precisely, we define the time dependent set and interface as

 Ωτ(t)={x:uk(x)≥0,forkτ≤t<(k+1)τ}andΓτ(t)=∂Ωτ(t).

Then as , (or ) has been shown to converge to motion by mean curvature in the viscosity setting [9, 2].

Now we describe some notations and the algorithm for the space-time discrete version of the above thresholding scheme (2) and (3). Let be a bounded, smooth domain, and be its boundary. Let be the spatial discretization step size. Define

 Ωh:={(m,n)∈Z2:dist[(nh,mh),Ω]≤h} (4)

which are the indices of the lattice points inside . Let again be the size of the time step. Given an initial set and its discrete version , the discrete thresholding scheme produces as follows. Let

 um,n0=1Ωh0(m,n)−1(Ωh0)c(m,n)={1for(m,n)∈Ωh0,−1for(m,n)∈(Ωh0)c, (5)

For simplicity, we will use to denote the discrete function . Then similar to the continuous space case, the following two steps are alternately performed (for ):

• solution of semi-discrete heat equation: for and ,

 {ddtwm,n(t)=1h2[wm+1,n(t)+wm−1,n(t)+wm,n+1(t)+wm,n−1(t)−4wm,n(t)],wm,n(0)=um,nk (6)
• thresholding step: for ,

 um,nk+1:=sign[wm,n(τ)]=sign((Sh(τ)[uk])m,n), (7)

where the is the solution operator of (6), i.e. .

Note that the above space time discrete scheme involves two small parameter and . Assume that and are related through

 h=C(τ)γ, where γ>0 and C>0 (8)

Three major cases are possible:

Case 1.

, i.e. , called “subcritical”.

Case 2.

, i.e. , where , called “critical”.

Case 3.

, i.e. , called “supercritical”.

Roughly speaking, the main result in Case 1 is that the level curves of a discrete heat equation move according to the motion by mean curvature. This gives the same result as [2]. Case 2 gives a version of anisotropic curvature dependent motion which demonstrates pinning of the interface when the curvature of the interface is too small. In Case 3, there is no motion at all.

### 1.1 Curvature Dependent Motion and Viscosity Solutions

As mentioned earlier, singularities and topological changes can occur for motion by mean curvature. Different mathematical approaches are invented to define the solution for all time. Due to the presence of maximum principle, we find the viscosity solution to be the most suitable and convenient for our problem. We spend a moment to briefly describe this method which can produce a “unique” global in time solution.

The essential idea is to represent the moving interface as the zero level set of a function :

 Γt={x:u(x,t)=0} (9)

The function is thus often called the level set function. It solves an appropriate partial differential equations related to the motion law of . The main result in this approach is that in the space of uniformly continuous functions, there is a unique solution and the set does not depend on the initial data as long as it correctly captures the interior and exterior domains of . On the other hand, this set-up does not a priori ensure that corresponds to a manifold in any geometric sense. This can happen if has positive -dimensional Lebesgue measure in which case is said to fatten or develop non-empty interior. It also means that the solution of the geometric evolution can be non-unique as . Conditions preventing this from happening are discussed in [3]. On the other hand, a definition of generalized front is used so that a “unique solution” can be defined. This approach defines the interface as the following triplet of objects:

 Γt={x:u(x,t)=0},D+t={x:u(x,t)>0},D−t={x:u(x,t)<0}. (10)

Figure 1.

(See Fig. 1.) Then () is called the interior (exterior) of the front . It is shown that the map:

 Et:(Γ0,D+0,D−0)⟶(Γt,D+t,D−t) (11)

is well defined. We refer to [4] for a more detailed description.

Next we describe the equation for general curvature dependent front propagation. We follow the exposition in [12]. Given an interface (hypersurface) in . We consider its motion described by a normal velocity function of the following form:

 V=V(Dν,ν) (12)

where is the unit (outward) normal of . (Note that is a symmetric matrix.) The above motion law is sometimes called anisotropic curvature motion. If , then the motion is called (isotropic) motion by mean curvature. If we want to represent the moving interface by (9) or (10), then the function needs to solve the following partial differential equation:

 ut+F(D2u,Du)=0 (13)

where the function is related to in the following way:

 F(X,p)=−|p|V(−|p|−1(I−¯p×¯p)X(I−¯p×¯p),−¯p) (14)

where . In order for the viscosity solution approach to work, the following monotonicity condition for is crucial:

 V is nondecreasing, i.e. for all X≤Y and p, then V(X,p)≥V(Y,p). (15)

The above property is translated to the function as

 for all X≤Y and p, then F(X,p)≥F(Y,p). (16)

With the above set-up, then it can be shown that equation (13) is well-posed in the space of uniformly continuous functions. Given an initial manifold , a usual choice of the initial data for (13) is given by the sign distance function to :

 d0(x)=sdist(x,Γ0)={dist(x,Γ0),x∈Ω0,−dist(x,Γ0),x∈Ωc0, (17)

The map in (11) is independent of any uniformly continuous initial data as long as it has the same sign as : on and on .

The definition of viscosity solution of (13) is given in Section 1.2 where a general approach to prove convergence of various approximating schemes is also given. Of the stability and consistency conditions generally required for most convergence proof, for the current algorithm, the former is quite easy to satisfy by means of maximum principle. The crux of the matter is the latter condition which is the key result of our paper for the case of space time discrete thresholding scheme. Once we have this, then we can more or less quote the general convergence result.

### 1.2 Main result: Sub-Critical Case

Our most complete result is the sub-critical case for which we can prove convergence to motion by mean curvature in the viscosity sense. In this case, the velocity function is given by . Hence (13) becomes

 ∂u∂t=|∇u|div(∇u|∇u|)=Δu−(D2uDu|Du)|Du|2 (18)

For convenience, we give the definition of viscosity solution specifically for this case.

###### Definition 1.

A locally bounded upper semicontinuous (usc) function (respectively, lower semicontinuous (lsc)) function is a viscosity subsolution (respectively, supersolution) of (18), if for all , and if is a local maximum point of , then one has

 ∂ϕ∂t(x,t)−(△ϕ−(D2ϕDϕ|Dϕ)|Dϕ|2)(x,t)≤0,if Dϕ(x,t)≠0, (19)

and

 ∂ϕ∂t(x,t)−△ϕ(x,t)+λmin(D2ϕ(x,t))≤0,if Dϕ(x,t)=0, (20)

where is the least eigenvalue of . (Respectively, if for all , and if is a local minimum point of , then one has

 ∂ϕ∂t(x,t)−(△ϕ−(D2ϕDϕ|Dϕ)|Dϕ|2)(x,t)≥0,if Dϕ(x,t)≠0, (21)

and

 ∂ϕ∂t(x,t)−△ϕ(x,t)+λmax(D2ϕ(x,t))≥0,if Dϕ(x,t)=0, (22)

where is the maximum (or principle) eigenvalue of .)

On the other hand, a simpler characterization can be given.

###### Proposition 1.

[2, Prop. 2.2] A locally bounded upper semicontinuous (usc) function is a viscosity subsolution (respectively supersolution) of (18) iff if satisfies (19) and

 ∂φ∂t(x,t)≤0  if Dφ(x,t)=0  and D2φ(x,t)=0, (\ref{SubDeg}’)

respectively, (21) and

 ∂φ∂t(x,t)≥0  if Dφ(x,t)=0  and D2φ(x,t)=0, (\ref{SuperDeg}’)

The consistency proof of the thresholding scheme relies on the following result

###### Proposition 2.

[2, Prop 4.1] If is a sequence of smooth functions bounded in and converging locally in to a function and is a sequence of points converging to such that , then if ,

 liminfh1h12⎛⎜⎝12−1(4πh)N2∫{ϕh(⋅,th−h)≥0}exp(−|xh−y|24h)dy⎞⎟⎠≥12√π|Dϕ(x,t)|(∂ϕ∂t−△ϕ+(D2ϕDϕ|Dϕ)|Dϕ|2)(x,t). (23)

Moreoever, if and and if the inequality

 12−1(4πh)N2∫{ϕh(⋅,th−h)≥0}exp(−|xh−y|24h)dy≤0 (24)

holds for a sequence of converging to , then

 ∂ϕ∂t(x,t)≤0 (25)

With the above preparation, we are ready to present our result. We start by introducing

 u––(x,t)=liminfkτ→t;(mh,nh)→xum,nk (26)

and

 ¯u(x,t)=limsupkτ→t;(mh,nh)→xum,nk (27)
###### Theorem 1 (Sub-critical case).

Assume . Then the functions and are viscosity subsolutions and supersolutions of (18), respectively.

A few remarks about the consequence of the above statement.

1. Let be the solution of (18) with the initial data given for example by (17). Then the sets , , and produced by the map (11) is well-defined. As shown in [2, Thm. 1.2], the above statement implies that

 u––(x,t)=1 in D+t and ¯u(x,t)=−1 in D−t, (28)

In other words,

 ∂{x:u––(x,t)≥1},∂{x:¯u(x,t)≤−1}⊂Γt.

It is in this sense that we say the zero level set of in the limit moves according to the motion by mean curvature. Note that in general we can only say the the limiting interface is contained in but might not equal to . See the next item of remark.

2. As we are dealing with discontinuous initial data and functions , , so in general the limit (as ) can be non-unique. In fact, we can infer using [2, Thm. 1.1, Cor. 1.3] that we have a unique limit if

 ⋃t>0Γt×{t}=∂{(x,t):u(x,t)>0}=∂{(x,t):u(x,t)<0}

i.e. with no fattening phenomena, and in which case the boundary of the zero level set for converges to in the sense of Hausdorff distance.

### 1.3 Main Results: Critical and Super-Critical Cases

We next describe the results in the critical case, i.e. for fixed . To concentrate on the key ideas of our approach, we assume that locally near the origin, the boundary of the initial set is represented by a graph. More specifically, we assume that for some ,

 Missing or unrecognized delimiter for \right (29)

where is a , even-function satisfying

 f(0)=f′(0)=0,f′′(0)=−κ≤0. (30)

The quantity represents the curvature of at . (See Fig. 2.)

Figure 2.

For simplicity of presentation of the results in this section, we consider an equivalent initialization and the threshholding step , where and are given by (5) and (7) respectively. In other words, we replace the “ ” thresholding scheme with a “ ” one. This new scheme will still be denoted with and . Consider the function obtained from (6) for . Let be such that

 w0,−n0(τ)≤12,andw0,−n0−1(τ)>12 (31)

For simplicity, we call the “discrete normal velocity” of the boundary point . The true physical normal velocity is defined as:

 V=limh,τ→0n0hτ=limh,τ→0n0μ (32)

We will show that exists and is still given by a function of the curvature . Precisely,

###### Theorem 2 (Critical case).

Assume . For sufficiently small , holds if and only if

 n−1∑k=1∫√2kμκ0e−x24dx+12∫√2nμκ0e−x24dx≤12∫∞√2nμκe−x24dx+∞∑k=n+1∫∞√2kμκe−x24dx. (33)

In particular, if , then

 n0−1∑k=1∫√2kμκ0e−x24dx+12∫√2n0μκ0e−x24dx=12∫∞√2n0μκe−x24dx+∞∑k=n0+1∫∞√2kμκe−x24dx. (34)

Note that discrete velocity is related to and implicitly. Nevertheless, it is straightforward to see from (34) that if , we have . It can also be seen that if is small enough, then , i.e. the interface is pinned. Using monotonicity, the above result can be extended to the case giving .

We next show that the result in the critical case is consistent with the sub-critical case when :

###### Theorem 3.

Let , and satisfy (31). Then

 limμ→∞n0μ=κ,

i.e. we get the mean curvature motion as in the subcritical case.

The main result in the super-critical case follows easily from Theorem 2.

###### Corollary 1 (Super-critical case).

Assume . If is sufficiently small, then , i.e. the front does not move.

The above statement indicates that the front will not move if either or is sufficiently small. From numerical calculation, we find that the smallness condition is quantified by .

Finally, we obtain an extension of Theorem 2 to the case of an anisotropic curvature motion. In particular, we want to calculate the normal velocity of a boundary point if the normal line at the point forms an angle with the coordinate axis. For concreteness, let the normal line at forms an angle with the -axis (measure in the counter-clockwise sense). Without loss of generality, we assume . We consider the case that for some positive integers .

For simplicity, we assume that locally in the neighborhood of , the boundary is given by , with , where . This way, the curvature of at is . In this setting the notion of normal motion has to be defined somewhat differently from the one in Theorem 2 since the line in the normal direction intersects the grid only at the points , and thus bypasses a number of grid points in between. To make the motion in the normal direction precise, denote to be the following strip:

 S:={(s,j):0≤s≤q−1;−∞

For every , let be the distance from to the line . Next, we reorder the elements in the set with respect to as follows:

 S:={(s1,j1),(s2,j2),(s3,j3),...} (36)

such that

 0

For example, if , the points in will be ordered as

 {(1,0),(2,0),...,(q−1,0),(0,−1),(1,−1),(2,−1),...,(q−1,−1),(0,−2),(1,−2),...}.

In this setting, we have the following result:

###### Theorem 4 (Anisotropic curvature motion, critical case).

Assume , , and are given, and the set be ordered as in (36). Fix an integer . Then, for sufficiently small , holds if and only if

 12∫√2dn0μκ0e−x24dx+n0−1∑i=1∫√2diμκ0e−x24dx≤12∫∞√2dn0μκe−x24dx+∞∑i=n0+1∫∞√2diμκe−x24dx (37)

In particular, if we assume that then

 12∫√2dn0μκ0e−x24dx+n0−1∑i=1∫√2diμκ0e−x24dx=12∫∞√2dn0μκe−x24dx+∞∑i=n0+1∫∞√2diμκe−x24dx (38)

In this case, is the normal displacement at time , thus is the normal velocity.

As the proof of all of our results relies heavily on the asymptotics of discrete heat kernel, we first collect their key properties and connection to the continuum heat kernel.

## 2 Properties of discrete heat kernels

### 2.1 Derivation of discrete heat kernel and its elementary properties.

We first consider a one-dimensional analog of (6)

 {unτ(t)=1h2[un+1(t)+un−1(t)−2un(t)],n∈Z,t≥0;un(0)=u0(nh) (39)

where the initial data is an function. The solution of the above is given by

 um(t)=∞∑k=−∞Gm−k(2th2)uk0 (40)

where

 Gn(α):=12π∫π−πcos(nξ)eα(cosξ−1)dξ=e−αI|n|(α) (41)

In the above is the Modified Bessel Function

 In(α)=∞∑m=0(α2)2m+nm!(m+n)!,n≥0.

In the following we will use the following notation,

 α=2th2 (42)

In order to establish (40), define the Fourier transform of a sequence as

 ^v(ξ,t):=1√2π∞∑m=−∞e−imξum(t).

Using (39), we conclude that

 ∂^v∂t=1√2π∞∑m=−∞e−imξumt(t)=1h2√2π∞∑m=−∞(eiξ+e−iξ−2)e−imξum(t)=2h2(cosξ−1)^v

and thus

 ^v(ξ,t)=^v(ξ,0)e2th2(cosξ−1).

Taking the inverse Fourier transform of , we recover :

 Missing or unrecognized delimiter for \right

Using separation of variables, it is straightforward to verify that the solution to the original problem (6) (in two dimensions) is given as

 wm,n(t)=∞∑s,j=−∞Gs−m(2th2)Gj−n(2th2)ws,j(0) (43)

The discrete heat kernel (41) has the following elementary properties:

 ∞∑k=−∞Gk(α)=1 (44)
 Gk(α)=G−k(α),k≥1 (45)
 12G0(α)+∞∑k=1Gk(α)=12 (46)

### 2.2 Decay properties.

The following result is used to control the Green’s function outside a fixed macroscopic domain.

###### Lemma 1.

Let be the discrete heat kernel (41) with . For any fixed , suppose . Then

 ∞∑k=[3μh]Gk(α)=o(e3)μh,h→0. (47)
###### Proof.

To simplify the notation, we omit the integer part brackets, denoting in (47) with . We have

 ∞∑k=3μhGk(α)=e−2τh2∞∑k=3μh∞∑n=0(τh2)2n+kn!(n+k)!

Shifting the summation and using the elementary inequality , we have

 ∞∑k=3μhGk(α)=e−2τh2∞∑m=0∞∑n=0(τh2)2n+m+3μhn!(n+m+3μh)!≤(τh2)3μh(3μh)!e−2τh2∞∑m=0∞∑n=0(τh2)2n+mn!(n+m)!≤(τh2)3μh(3μh)!

since, by property (44),

 e−2τh2∞∑m=0∞∑n=0(τh2)2n+mn!(n+m)!≤∞∑m=−∞Gm(2τh2)=1.

Using the lower bound for the factorial (Stirling approximation [1]),

 (3μh)!≥√2π3μh(3μeh)3μh.

Thus by the assumption , we have

 ∞∑k=3μhGk(α)≤1√2π3μh(eτ3μh)3μh=o(e3)μh,

the desired statement ∎

### 2.3 Asymptotic expansions.

In what follows, we are going to use the asymptotic expansions for modified Bessel function, [1, p. 199] (see also [13]):

• The expansion for for any fixed index is given by:

 Iν(z)=ez√2πz(1−4ν3−18z+O(1z2)),asz→∞. (48)
• Meissel formula [13], which holds uniformly for and :

 Iν(νz)=ννeνΓ(ν+1)eνη(1+z2)14(1+∞∑k=1vk(ξ)νk) (49)

where

 ξ=1√1+z2,η=√1+z2+lnz1+√1+z2,

and is a polynomial of the form

 vk(ξ)=k∑s=0a(k)sξk+2s, (50)

defined through the recursive relation and

 vk+1(ξ)=12ξ2(1−ξ2)v′k(ξ)+18∫ξ0(1−5s2)vk(s)ds
###### Proposition 3.

Assume is such that as . Then the following asymptotic expansions hold.

• For fixed ,

 Gn(2τh2)=1√4πh√τ(1−h2τ4n3−116+O(h4τ2)),h→0 (51)
• Let . Then there exists an absolute constant and such that

 G√τxh(2τh2)=1√4πh√τe−x24(1+Ch,xhx√τ) (52)
###### Proof.

Expansion (51) follows from (41) and (48):

 Gn(2τh2)=e−2τh2In(2τh2) (53)

For (52), we will make use of the Meissel formula (49). Borrowing the notations there, we set

 ν=x√τh,z=2√τxh.

Then

 η=√1+z2+lnz1+√1+z2=√x2h2+4τxh+ln2√τhx+√x2h2+4τ (54)

so that

 eην=e√τh2√x2h2+4τ(2√τhx+√x2h