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 $\triangle x = h$ and $\triangle t = \tau$. We have identified the limiting description of the interfacial velocity in the (i) sub-critical ($h \ll \tau$), (ii) critical ($h = O(\tau)$), and (iii) super-critical ($h \gg \tau$) 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).


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 "subcritical" 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 W is the double well potential W (u) = (1 − u 2 ) 2 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 u takes on the values 1 and −1 which are the minima of W . The function u also makes a smooth but rapid transition with thickness O( ) between the two domains. The key is then to understand the dynamics of this transition layer, in the limit of −→ 0. 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 1). Given an initial shape Ω 0 , and its boundary (or often called the interface) Γ 0 = ∂Ω 0 , a sequence of functions {u k } k≥0 is constructed in the following manner: for k = 0, we define then the following two steps are alternately performed (for k = 0, 1, 2, . . .), • diffusion step: • thresholding step: Note that the second step above is to mimick the fast reaction term which drives u to 1 or −1, the minima of W . The solution of the problem is captured by the sequence of subsets where u k attains the value 1, {x : u k (x) = 1}. Precisely, we define the time dependent set and interface as Then as τ −→ 0, Ω τ (t) (or Γ τ (t)) 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 Ω ⊂ R 2 be a bounded, smooth domain, and Γ = ∂Ω be its boundary. Let h > 0 be the spatial discretization step size. Define which are the indices of the lattice points inside Ω. Let again τ > 0 be the size of the time step. Given an initial set Ω 0 and its discrete version Ω h 0 , the discrete thresholding scheme produces u m,n k k≥0,(m,n)∈Z 2 as follows. Let For simplicity, we will use u to denote the discrete function u m,n : (m, n) ∈ Z 2 . Then similar to the continuous space case, the following two steps are alternately performed (for k = 0, 1, 2, . . .): • solution of semi-discrete heat equation: for (m, n) ∈ Z 2 and 0 < t ≤ τ , • thresholding step: for n, m ∈ Z, where the S h (·) is the solution operator of (6), i.e.
Note that the above space time discrete scheme involves two small parameter τ > 0 and h > 0. Assume that τ and h are related through Three major cases are possible: Case 2. γ = 1, i.e. τ = µh, where µ = const, called "critical".
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.

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 Γ t as the zero level set of a function u(x, t): The function u is thus often called the level set function. It solves an appropriate partial differential equations related to the motion law of Γ t . The main result in this approach is that in the space of uniformly continuous functions, there is a unique solution u and the set Γ t does not depend on the initial data u(·, 0) as long as it correctly captures the interior and exterior domains of Γ 0 . On the other hand, this set-up does not a priori ensure that Γ t corresponds to a manifold in any geometric sense. This can happen if Γ t has positive n-dimensional Lebesgue measure in which case Γ t is said to fatten or develop non-empty interior. It also means that the solution of the geometric evolution can be non-unique as ∂ {x : u(x, t) > 0} = ∂ {x : u(x, t) < 0}. 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: (See Fig. 1 is called the interior (exterior) of the front Γ t . 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 R n . We consider its motion described by a normal velocity function V of the following form: where ν is the unit (outward) normal of Γ. (Note that Dν is a symmetric matrix.) The above motion law is sometimes called anisotropic curvature motion. If V = div(Dn), 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 u needs to solve the following partial differential equation: where the function F is related to V in the following way: whereq = |p| −1 p. In order for the viscosity solution approach to work, the following monotonicity condition for V is crucial: The above property is translated to the function F as for all X ≤ Y and p, then F (X, p) ≥ F (Y, p).
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 Γ 0 = ∂Ω 0 , a usual choice of the initial data u(x, 0) = u 0 (x) for (13) is given by the sign distance function to Γ 0 : The map E t in (11) is independent of any uniformly continuous initial data as long as it has the same sign as d 0 : u 0 > 0 on Ω 0 and u 0 < 0 on Ω c 0 . 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.

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 V (Dn, n) = Div(Dn). Hence (13) becomes For convenience, we give the definition of viscosity solution specifically for this case.

Proposition 1. [2, Prop. 2.2] A locally bounded upper semicontinuous (usc) function u 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 Moreoever, if Dφ(x, t) = 0 and D 2 φ(x, t) = 0 and if the inequality holds for a sequence of h converging to 0, then With the above preparation, we are ready to present our result. We start by introducing Theorem 1 (Sub-critical case). Assume γ > 1. Then the functionsū(x, t) and u(x, t) are viscosity subsolutions and supersolutions of (18), respectively.
A few remarks about the consequence of the above statement.
1. Let u be the solution of (18) with the initial data u 0 given for example by (17). Then the sets Γ t , D + t , and D − t produced by the map E t (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 u k 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 Γ t . See the next item of remark.
2. As we are dealing with discontinuous initial data and functions u k , k ≥ 0, so in general the limit (as h, τ −→ 0) can be non-unique. In fact, we can infer using [2, Thm. 1.1, i.e. with no fattening phenomena, and in which case the boundary of the zero level set for u k converges to t>0 Γ t × {t} in the sense of Hausdorff distance.

Main Results: Critical and Super-Critical Cases
We next describe the results in the critical case, i.e. τ = µh for fixed µ > 0. To concentrate on the key ideas of our approach, we assume that locally near the origin, the boundary of the initial set Ω 0 is represented by a graph. More specifically, we assume that for some c 0 > 0, where f (x) is a C 2 , even-function satisfying The quantity κ represents the curvature of ∂Ω 0 at (0, 0). (See Fig. 2.) For simplicity of presentation of the results in this section, we consider an equivalent initializationũ m,n 0 := 1 2 + 1 2 u m,n 0 and the threshholding stepũ m,n k+1 := 1 2 + 1 2 u m,n k , whereũ m,n 0 andũ m,n k+1 are given by (5) and (7) respectively. In other words, we replace the "−1 & 1" thresholding scheme with a "0 & 1" one. This new scheme will still be denoted with u m,n 0 and u m,n k+1 . Consider the function w obtained from (6) for k = 0. Let n 0 ∈ Z be such that For simplicity, we call n 0 the "discrete normal velocity" of the boundary point (0, 0) ∈ ∂Ω 0 . The true physical normal velocity V is defined as: We will show that V exists and is still given by a function of the curvature κ. Precisely, In particular, if w 0,−n 0 (τ ) = 1 2 , then Note that discrete velocity n 0 is related to µ and κ implicitly. Nevertheless, it is straightforward to see from (34) that if κ 1 < κ 2 , we have n 0 (µ, κ 1 ) ≤ n 0 (µ, κ 2 ). It can also be seen that if κ is small enough, then n 0 = 0, i.e. the interface is pinned. Using monotonicity, the above result can be extended to the case κ = 0 giving n 0 = 0.
We next show that the result in the critical case is consistent with the sub-critical case when µ → ∞: Theorem 3. Let n 0 , µ 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.
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 µκ ≤ 0.8218.
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 (0, 0) ∈ ∂Ω forms an angle θ with the x-axis (measure in the counterclockwise sense). Without loss of generality, we assume 0 < θ ≤ π 4 . We consider the case that tan θ = p q for some positive integers p ≤ q. For simplicity, we assume that locally in the neighborhood of (0, 0), the boundary ∂Ω is given by ( . This way, the curvature 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 (np, nq), n ∈ Z, and thus bypasses a number of grid points in between. To make the motion in the normal direction precise, denote S to be the following strip: For every (s, j) ∈ S, let d(s, j) := |jq−sp| √ p 2 +q 2 be the distance from (s, j) to the line y = p q x. Next, we reorder the elements in the set S with respect to d as follows: For example, if p = 1, the points in S will be ordered as In this setting, we have the following result: Theorem 4 (Anisotropic curvature motion, critical case). Assume p, q, κ and µ are given, and the set S be ordered as in (36). Fix an integer n 0 ≥ 0. Then, for sufficiently small t, w sn 0 ,jn 0 (τ ) ≤ 1 2 holds if and only if In particular, if we assume that w sn 0 ,jn 0 (τ ) = 1 2 then In this case, d n 0 h is the normal displacement at time τ , thus dn 0 µ 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 u 0 is an L ∞ function. The solution of the above is given by In the above I |n| (α) 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 u m (t) aŝ v(ξ, t) : Using (39), we conclude that Taking the inverse Fourier transform ofv(ξ, t), we recover u m (t): 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:

Decay properties.
The following result is used to control the Green's function outside a fixed macroscopic domain.
Proof. To simplify the notation, we omit the integer part brackets, denoting 3µ Shifting the summation m = k − 3µ h and using the elementary inequality Using the lower bound for the factorial (Stirling approximation [1]), Thus by the assumption τ ≤ µh, we have the desired statement
We start with observing that it follows from Stirling approximation ( [1]) that for all ν ≥ 1: ν ν e ν Γ(ν + 1) where 0 ≤ C ν ≤ C 0 and c 0 is an absolute constant. Consequently, We now proceed with expanding each term: which can be written as Using the expansion . Using Taylor expansion, it becomes: and, therefore, Combining the four expansions above, (56) reads as 3 Proof of Theorem 1 -Sub-Critical Case.
We are going to prove that the functionū, defined in (27), is a viscosity subsolution of the equation (18). The fact that (26) is a supersolution of the same equation can be established analogously. The strategy of proof follows [2] closely. Essentially we will demonstrate how to reduce the discrete computation to [2,Prop. 4.1] which is the key consistency result needed.
Step I. First let ϕ(x, t) be a smooth test function satisfying Indeed, it follows from the coercivity condition (58) that max Z 2 ×N (u m,n k − ϕ(mh, nh; kτ )) is attained at some m = m h , n = n h an k = k τ , and up to a subsequence, there exist x = (x 1 ,x 2 ) andt such that hm h →x 1 , hn h →x 2 and τ k τ →t. This point (x,t) must be the global maximum ofū − ϕ. In addition, (60) holds, for otherwise it will contradict (59) and the definition ofū. Using the same reasoning in [2, p. 490], we conclude that the above inequality leads to where sign * is the upper semi-continuous envelope of the sign function. The above is equivalent to In terms of discrete kernels, the above inequality reads as Step II. In this step, we will show how (61) leads to (20') and (22'). We first express the left hand side of (61) as an integral. For this, we write Claim. For any G ⊂ R 2 , f : G → R and δ > 0 denote RS δ [f, G] to be the Riemann sum for f in G with step δ, i.e.  and G ⊂ R 2 we have the following error estimate between the Riemann sum and the integral: for some C > 0 which does not depend on G and δ. The estimate (66) follows from a 1D inequality which holds uniformly for all −∞ ≤ a < b ≤ +∞.
Now we proceed to analyze (64). Due to (52), we have For i = 1, 2, by (51): (69) Applying (66) to (68) and making use of (69), inequality (64) becomes with Roughly speakingQ h is a discretized set. In order to apply [2, Prop. 4.1] we need to obtain an estimate (70) in whichQ h is replaced with its continuum analog with well-controlled error. Consider Since ϕ is a smooth function which satisfies the growth condition (58), ∂Q h is bounded.
. This way, for any ( This enables us to say that which, in turn, implies that with C independent on h. After rescaling, we have Combining (70) with (72) and (73), we have Next we consider two cases. If Dφ((x 1 ,x 2 );t) = 0, we apply the result of [2,Prop. 4 which yields the desired result.
It remains to show that in the case when Dφ((x 1 ,x 2 );t) = 0, D 2 φ((x 1 ,x 2 );t) = 0, and with M h satisfying (71) for sufficiently small h, we have ∂φ ∂t ((x 1 ,x 2 );t) ≤ 0. The condition (74), strictly speaking, is different from the one in Proposition 4.1 [2], yet enables us to apply the technique developed in [2]. In particular, three cases are possible: In Case I, arguing as in [2], we deduce that for any c > 0. However, since passing to the limit in (75) as h → 0, using (74), we have ∂φ ∂t (x 1 , x 2 ; t) − c ≤ 0 for all c > 0, which yields the desired result.
The conclusion in the Cases II and III follows from Proposition 4.1 [2] without any change.
This section provides the proofs for the statements related to the critical case.

Proof of Theorem 2
The results will follow from the following statements.
(i) (Lower bound) If w 0,−n 0 (τ ) ≥ 1 2 we have We refer to page 9 for the notation concerning the set Ω. Assume n 0 ≥ 0 is such that w 0,−n 0 (τ ) ≥ 1 2 . By Lemma 1, Having shifted the summation indices, we get Similarly, if w 0,−n 0 (τ ) ≤ 1 2 , we have We briefly outline the idea of the proof. The relations (78) and (79) implicitly define n 0 in terms of h. Our goal is to obtain the description of n 0 , independent of h. To this end, we are going to establish the asymptotic (in h) expansion of the sum in the right hand side of (78) and (79): where Φ(n 0 , µ, κ) := and the desired result follows.
It is convenient to group the indices in the summation (80) into the following sets: With the above, we have Since we have Changing the order of summation in (84) using the notation (85), we have if w 0,−n 0 (τ ) ≥ 1 2 , then Using the same reasoning, we get that if w 0,−n 0 (τ ) ≤ 1 2 , we have Roughly speaking, the main technical difficulty in establishing the asymptotic behavior of the terms in (87) is caused by the term since it involves an infinite sum of the products of Bessel functions and Riemann sums. To proceed, we are going to estimate this term from below using finitely many such products, thus obtaining the lower bound (76). On the other hand, we are going to estimate this term from above, using a priori bounds for Bessel functions, and get the upper bound (77). Once we verify that those bounds match, we will obtain the asymptotic expansion of all terms in (87). Similarly, and Combining (91), (92), (93) and (94), the inequality (88) implies and since N ≥ 1 was chosen arbitrarily, letting N −→ ∞, we arrive at: Thus, (i) follows.

Proof of Theorem 3: the limiting case µ → ∞.
Recall that n 0 satisfies We first start with the following claim: as µ → ∞, up to a subsequence, we have lim µ→∞ n 0 µκ = a, with 0 < a < ∞.
First, suppose a = ∞, i.e. n 0 µ. Then, for sufficiently large µ, we have 2k µκ ≥ 1 for all k ≥ n 0 + 1. Hence, we can estimate the right hand side of (102) in the following way: It remains to show that a = 1. Suppose a > 1. Omitting the integer part notation, assume that n 0 = aµκ, then Consequently, (106) holds, which yields a contradiction with (105) if a > 1. The case a < 1 is excluded analogously.
On the other hand, the properties (44) and (45) of the heat kernel yield that (s,j)∈I θ 2 ∪I θ Evolution of a circle.
Evolution of a dumpbell.
We investigate in detail the motion of a circle. The figures below show how the radius of a circle changes for several values of the parameter µ between 0.3 and 1, as well as for various values of h between 0.05 and 0.1. The black circled curve is the exact mean curvature motion.
Evolution of the circle radius with µ.
Evolution of the circle radius with h.
The next figures illustrate the anisotropic behavior described in Theorem 4. We start with a finger-shaped domain, tilted at angle 45 • and verify that it moves with constant velocity.