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


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:


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

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

  • diffusion step:

  • thresholding step:


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

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


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


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 ,

  • thresholding step: for ,


    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


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 :


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:


Figure 1.

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


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:


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:


where the function is related to in the following way:


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

is nondecreasing, i.e. for all and , then . (15)

The above property is translated to the function as

for all and , then . (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 :


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


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




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




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


respectively, (21) and


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 ,


Moreoever, if and and if the inequality


holds for a sequence of converging to , then


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



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


    In other words,

    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

    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 ,


where is a , even-function satisfying


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


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


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


In particular, if , then


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

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:


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


such that

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

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


In particular, if we assume that then


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)


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




In the above is the Modified Bessel Function

In the following we will use the following notation,


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

Using (39), we conclude that

and thus

Taking the inverse Fourier transform of , we recover :

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


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


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


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

Shifting the summation and using the elementary inequality , we have

since, by property (44),

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

Thus by the assumption , we have

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:

  • Meissel formula [13], which holds uniformly for and :



    and is a polynomial of the form


    defined through the recursive relation and

Proposition 3.

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

  • For fixed ,

  • Let . Then there exists an absolute constant and such that


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


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



so that