ON THE SELF-SIMILAR SOLUTIONS OF THE CRYSTALLINE
MEAN CURVATURE FLOW IN THREE DIMENSIONSNORBERT POŽÁR

Abstract. We present two types of self-similar shrinking solutions of positive
genus for the crystalline mean curvature flow in three dimensions analogous
to the solutions known for the standard mean curvature flow. We use themto test a numerical implementation of a level set algorithm for the crystalline
mean curvature flow in three dimensions based on the minimizing movements
scheme of A. Chambolle, Interfaces Free Bound. 6 (2004). We implement a
finite element method discretization that seems to improve the handling ofedges in three dimensions compared to the standard finite diﬀerence method
and illustrate its behavior on a few examples.

1. Crystalline mean curvature flow
The understanding of the evolution of small crystals has been a challenging problem of material science and mathematical modeling. In this regime, the evolution
seems to be governed by the surface energy, whose eﬀects are usually modeled bymean curvature terms. Due to the lattice structure of a typical crystal, the surface
energy density is anisotropic. In fact, it is postulated that the optimal shape (Wulﬀshape) is a convex polytope and such anisotropies are called crystalline. This causes
diﬃculties for the definition of an anisotropic (crystalline) mean curvature and asuitable notion of solutions of the resulting surface evolution problem, and makes
the development of an eﬃcient numerical method challenging.
The crystalline mean curvature was introduced independently by S. B. Angenentand M. E. Gurtin [5] and J. E. Taylor [57] to model the growth of small crystals, seealso [11, 37]. The surfaces of solid and liquid bodies have a surface energy, which
is usually expressed as the surface integral of a surface energy density σ : S n−1 →(0, ∞) over the boundary of a set E ⊂ Rn , representing the body,∫
F(E) :=
σ(ν) dS,∂E

where ν : ∂E → S
is the unit outer normal of E. Here n is the dimension, usually
2 or 3 in applications. For many materials, especially liquids, σ is given by thesurface tension coeﬃcient and is therefore constant on the unit sphere S n−1 . Thissurface energy is a manifestation of the fact that the atoms or molecules forming
the body have a smaller interaction energy when surrounded by the particles of then−1

Date: June 7, 2018.
2000 Mathematics Subject Classification. 53C44, 35K93, 65M99.
Key words and phrases. crystalline mean curvature flow, level set method, self-similarsolutions.
This work was partially supported by JSPS KAKENHI Grant No. 26800068 (Wakate B) and
No. 18K13440 (Wakate).
1

2

N. POŽÁR

x2

∂E

W
ν
σ<1
x1

E
Figure 1. Left: Crystal of atoms on a regular lattice. Right: The
1-level set of the associated surface energy density σ and its Wulﬀ
shape W.
same kind. Liquids, typically, do not have any preferred direction in the distributionof particles, and therefore the surface energy density is isotropic.
The situation is quite diﬀerent for crystals. Let us give a simplified illustration.
If we suppose that the atoms are distributed along the regular square lattice in two
dimensions, every atom inside the body has exactly four neighbors with which it
creates chemical bonds, see Figure 1. On the surface, however, some of these bonds
are broken where a neighbor is missing and the surface energy is proportional to
the number of the broken bonds. This number is given in terms of the taxicab (ℓ1 )length of the surface, not the usual Euclidean length. In particular, in this caseσ(ν) ∼ ∥ν∥1 := |ν1 |+|ν2 |, where ν is the macroscopic unit outer normal. This is thebasic motivation for the introduction of nondiﬀerentiable surface energy densities.
For convenience, we will assume that σ is positively one-homogeneously extended
from S n−1 to Rn asp),p ∈ Rn \ {0},σ(p) = |p|σ( |p|∑σ(0) = 0, where |p| := ( i |pi |2 )1/2 is the usual Euclidean norm.
If σ is convex on Rn and σ(p) > 0 for all p ∈ Rn \{0}, we call it an anisotropy. We
are in particular interested in anisotropies
σ that are piece-wise linear, for instance∑the ℓ1 -norm σ(p) = ∥p∥1 := i |pi |. Such anisotropies will be called crystalline
anisotropies.
The optimal shape of the crystal, that is, the shape with minimal surface energy
for a given volume, is a translation and scaling of the Wulﬀ shape(1.1)

W := {x : x · p ≤ σ(p), p ∈ Rn },see [56].
The evolution {Et }t≥0 of a body driven by the dissipation of the surface energy
then leads to a formal gradient flow
(1.2)

V = β(ν)(κσ + f )on ∂Et ,where V is the normal velocity of the surface ∂Et , f = f (x, t) is a given external
force, and β : S n−1 → (0, ∞) is a mobility. Finally, −κσ is the first variation of thesurface energy F at Et . κσ is usually called the anisotropic mean curvature of thesurface.

SOLUTIONS OF THE CRYSTALLINE MEAN CURVATURE FLOW

3

If σ ∈ C 2 (Rn \ {0}) and {p : σ(p) < 1} is strictly convex, then it is well-known[26] that the anisotropic mean curvature can be evaluated as
κσ = − div∂E (∇σ(ν)),
where div∂E is the surface divergence on ∂E.
When σ is crystalline, the situation is significantly more complicated. In particular, even if ∂E is smooth, ∇σ(ν) might be discontinuous (or not even defined) on
parts of the surface. Therefore κσ might not be defined, or might not be a function.
Instead, following for example [11], we define the subdiﬀerential of σ as
∂σ(p) := {ξ ∈ Rn : σ(p + h) − σ(p) ≥ ξ · h, h ∈ Rn },where ξ · h is the usual inner product on Rn . Note that ∂σ(p) is a nonempty
compact convex subset of Rn . We replace ∇σ(ν) by a vector field z : ∂E → Rn ,usually called a Cahn-Hoﬀman vector field, that is a selection of ∂σ(ν(x)) on ∂E,that is, z(x) ∈ ∂σ(ν(x)), x ∈ ∂E. However, now there are multiple choices of z
which potentially lead to diﬀerent values of κσ = − div∂E z. It turns out that areasonable choice is a vector field zmin that minimizes ∥− div∂E z + f ∥L2 (∂E) . The
crystalline (mean) curvature is then defined asκσ := − div∂E (zmin ).
Such a choice is motivated by the standard theory of monotone operators due toY. Kōmura and H. Brézis [16, 41]. Furthermore, since the Euler-Lagrange equationof the minimization problem is ∇(− div∂E z + f ) = 0, this choice yields κσ + f thatis constant, if possible, on flat parts, or facets, of the crystal parallel to the flat partsof the Wulﬀ shape W. Therefore facets are usually preserved during the evolution,
as expected. However, κσ might be even discontinuous on facets, and then facet
breaking or bending occurs, see [14] and Figure 14, Figure 15. This poses a serious
diﬃculty for introducing a suitable notion of solutions for this problem. Since κσis itself given as a solution of a minimization problem, it is in general diﬃcult toevaluate it, except in special circumstances. Moreover, κσ is a nonlocal quantity
on the facets of the crystal as the following example shows.
∑n1.1. Example. Consider the cubic anisotropy σ(p) = ∥p∥1 := i=1 |pi | and suppose that the initial shape is the cube centered at 0 with side-length L0 > 0, β ≡ 1.
Let us try to find {Et }t≥0 . It is not diﬃcult to see that for a cube QL = (− L2 , L2 )n ,L > 0, the vector field z(x) = σ◦x(x) is a Cahn-Hoﬀman vector field on ∂QL , whereσ ◦ (x) := sup {x · p : σ(p) ≤ 1}. Here σ ◦ is the convex polar of σ and it is the dual
norm of σ, σ ◦ (x) = ∥x∥∞ := max1≤i≤n |xi |, see [52]. In particular, σ ◦ (x) = L2 on∂QL . Therefore div∂QL z = L2 (n − 1). Since it is a constant on the facets, z actually
minimizes ∥− div z∥L2 (∂QL ) among all Cahn-Hoﬀman vector fields, and thereforeκσ = − L2 (n−1) if f ≡ const. We deduce that, if L20 (n−1) ≥ f ≡ const so that solution is shrinking, the solution of (1.2) is {Et }t≥0 , Et = QL(t) , where L(0) = L0 and√
L′ = − L2 (n − 1) + f . When f = 0, the unique solution is L(t) = L20 − 4(n − 1)t.
L20
the cube vanishes.
At the extinction time t1 := 4(n−1)
If the forcing term f is strong enough, the crystal will grow. However, it will
only stay a cube as long as the velocity of corners is less than f since β ≡ 1. If the
speed of the corners is bigger, the corners will round up, as can be easily seen by
the comparison principle. See for example [32] and references therein.

4

N. POŽÁR

Solutions of the crystalline mean curvature flow. Introducing a notion of solutionsfor (1.2) with the crystalline anisotropy have been a challenging problem. In two
dimensions, if f is constant on facets, the situation is somewhat simpler since κσis constant on facets of the crystal parallel to the facets of the Wulﬀ shape W.
Therefore if the initial shape is a polygon with edges parallel to edges of the Wulﬀ
shape, the facets will move without breaking or bending. Their evolution canbe tracked by the crystalline algorithm [57], which also yields eﬃcient numericalmethods. However, these methods cannot treat evolutions that are not strictly
faceted. For fully general situations, the level set method was successfully used to
introduce a notion of viscosity solutions to (1.2) in two dimensions by M.-H. Gigaand Y. Giga [27–29], and a numerical algorithm was developed by A. Chambolle [17]and further extended to the crystalline case by A. Oberman, S. Osher, R. Takei
and R. Tsai [48].
In three dimensions, the situation is significantly more complicated by the possible bending or breaking of facets. There is an extensive number of publications
that is beyond the scope of this paper, for instance [12–15], and [11] for an introduction to the topic and references. Recently, A. Chambolle, M. Morini and
M. Ponsiglione [19], introduced a well-posed notion of solutions for the particularvelocity law V = σ(ν)κσ . In a subsequent paper with M. Novaga [18], they generalized the theory to V = β(ν)(κσ + f ), where β is an anisotropy and f = f (x, t)
is a Lipschitz continuous function. They define solutions using the signed distancefunction to the evolving set, which is required to be a solution of a certain PDE ina sense of distributions, and prove the existence of such a solution using the minimizing movements algorithm. Independently, Y. Giga and the author introduced
a well-posed notion of viscosity solutions for the level set formulation of (1.2) inthe fully general form V = F (ν, κσ + f ) where the continuous nonlinearity F is
nondecreasing in the second variable, but with constant f and for bounded crystals [33, 34]. See also Section 2.1. It was shown in [18] that both of these notions
coincide whenever they both apply.
As for the available numerical results, published results concerning the purely
crystalline anisotropy so far seem to only treat the two dimensional evolution. However, the algorithm proposed in [17, 48] generalizes naturally to three dimensions,and can easily accommodate a general external force as explained in Section 3.
We present some of the results of this implementation below. Let us also mention
the three dimensional results of J. W. Barrett, H. Garcke and R. Nürnberg [6–8],
who develop a parametric finite element method for the anisotropic mean curvature
flow and apply it to the Stefan problem with Gibbs-Thomson law that features analmost-crystalline, but still smooth, anisotropic curvature. This method does notseem to be able to handle topological changes. For a flow with topological changes
their phase field method [9, 10] is available. See also [25] for a survey of numericalapproaches.
Self-similar solutions. The study of self-similar solutions of the classical mean curvature flow has been important for the understanding of the singularities of the flow.
In two dimensions, it is known that any simple initial closed curve will become aboundary of a convex set in a finite time and therefore the only compact embedded
self-similar solution is the circle [24, 36].
In three dimension the situation is more interesting. It is known that the sphere is
the only convex self-similar solution [38]. The first embedded nonconvex self-similar

SOLUTIONS OF THE CRYSTALLINE MEAN CURVATURE FLOW

5solution was the “shrinking doughnut” solution constructed by S. B. Angenent [4].
This solution can be used to rigorously show the neck-pinching singularity startingwith a dumbbell-shaped initial data. Self-similar solutions of higher genus werediscovered numerically by D. L. Chopp [21], but to the author’s knowledge their
existence have not been proven rigorously. The construction of the first highergenus embedded self-similar solution was done by X. H. Nguyen [44–46], but thissolution is diﬀerent from the one found by Chopp in [21].
The behavior is more complex in the crystalline case. There are convex selfsimilar solutions other than the Wulﬀ shape even in two dimensions. For example,
any axes-aligned rectangle will generate a self-similar solution of the crystalline
curvature flow with the ℓ1 anisotropy. For “non-rectangular” even anisotropies, theWulﬀ shapes in two dimensions are stable [54]. The stability of the Wulﬀ shapesolution in three dimensions was considered in [47], with more examples of nonWulﬀ convex solutions. There are also examples of nonconvex self-similar solutions
in two dimensions [40]. For a construction of self-similar solutions in a sector see[30]. The self-similar solutions of positive genus constructed below appear to be
new.
Main results and the outline. We construct the crystalline shrinking doughnut
in Section 2.2 and the crystalline sponge-like solution in Section 2.3. In Section 3
we present an implementation of a numerical algorithm for the crystalline mean
curvature flow and test its accuracy using the self-similar solutions.
2. Self-similar solutions with positive genus
It is known that the solution of V = σ(ν)κσ with initial data given by the Wulﬀ
shape W of the anisotropy σ is self-similar up to the vanishing time, see Example 1.1.
In this section we construct self-similar solutions of the anisotropic or crystalline
mean curvature flow with positive genus. One is a torus-like solution, analogous
to the solution constructed by S. B. Angenent for the isotropic mean curvature
flow [4], and the other is a sponge-like solution of genus 5, similar to the solution
for the isotropic mean curvature flow discovered numerically by D. L. Chopp [21].
Such solutions are useful for understanding possible singularities of the flow, forinstance showing that a neck-pinching occurs [4]. Furthermore, the solutions canbe constructed explicitly for certain anisotropies and may be therefore useful for
testing numerical methods, as we will do in Section 3.4.
2.1. Notion of solutions using the level set method. We need to first clarify
what we mean by a solution of the crystalline mean curvature flow (1.2). The levelset method for the mean curvature flow was introduced and developed in [20, 23, 49].
The basic idea is to introduce an auxiliary function u : Rn × [0, ∞), whose evolutionof every level set {{x ∈ Rn : u(x, t) < c}}t≥0 , c ∈ R, satisfies the velocity law (1.2).
It is easy to see [26] that in this case
V =−

ut
,
|∇u|

ν=

∇u
,
|∇u|

and

κσ = − div(∇σ(∇u)).

Therefore u formally satisfies the equation
(
)
∇u
ut
=β
(2.1)
(− div(∇σ(∇u)) + f )
−
|∇u|
|∇u|

in Rn × 0.

6

N. POŽÁR

If σ is a crystalline anisotropy, ∇σ might be discontinuous and therefore thediﬀerential operator on the right-hand side is very singular. In fact, it is a nonlocal
operator on the flat parts of the surface of the crystal parallel to the flat parts of theWulﬀ shape. Therefore this equation does not fit within the classical framework ofviscosity solutions for geometric equations [20, 23]. The extension of the viscosity
theory to (2.1) had been a challenging open problem. In one dimension, which also
covers two-dimensional crystals, the theory was developed by M.-H. Giga, Y. Giga,
P. Rybka and others [27, 29, 31]. Y. Giga and the author recently introduced a new
notion of viscosity solutions for (2.1) with f independent of the space variable thatapplies to the crystalline anisotropy [33, 34]. This notion is well-posed for bounded
crystals and stable with respect to a regularization of the anisotropy, that is, withrespect to the approximation of the crystalline curvature by smooth anisotropic
curvatures. The main idea is to restrict the space of test functions to faceted test
functions for which it is possible define the operator div[∇σ(∇u)] as the divergenceof a minimizing Cahn-Hoﬀman vector field z ∈ ∂σ(∇u) as explained above. The
solutions constructed below are viscosity solutions of the level set formulations.
2.2. The shrinking doughnut. In [4], S. B. Angenent showed the existence of adoughnut-like self-similar shrinking solution for the isotropic mean curvature flowin dimensions n ≥ 3.
An analogue of this solution can be also constructed in the anisotropic case. Let
us fix a dimension n ≥ 3. We consider a “cylindrical” anisotropyσ(p) := σ̃(p′ ) + |pn |,

p = (p′ , pn ) ∈ Rn ,where σ̃ is an even anisotropy on Rn−1 (not necessarily crystalline), i.e., we requirethat σ̃(−p′ ) = σ̃(p′ ). We define the mobility(2.2)β(p) := σ̃(p′ ) + µ|pn |,

p = (p′ , pn ) ∈ Rn ,where µ > 0 is a given constant. We will see later that we must take µ = 12 forn = 3 to get a self-similar solution. Note thatσ ◦ (x) = max(σ̃ ◦ (x′ ), |xn |),x = (x′ , xn ) ∈ Rn ,{}where σ̃ ◦ (x′ ) := sup x′ · p′ : p′ ∈ Rn−1 , σ̃(p′ ) ≤ 1 is the convex polar of σ̃, see for
instance [52] for this and other convex analysis results.
We consider the anisotropic mean curvature flow
(2.3)

V = β(ν)κσ .

Let R > r > 0, h > 0 be real parameters. Define the set
Tr,R,h := {x = (x′ , xn ) ∈ Rn : r < σ̃ ◦ (x′ ) < R, |xn | < h}.
This is a “torus” with the hole aligned with the xn -direction, see Figure 2. We will
use this set to construct a self-similar evolution{
Tr(t),R(t),h(t) , 0 ≤ t < t∗ ,Et :=∅,t ≥ t∗ ,given appropriate functions r(t), R(t) and h(t), where t∗ is the extinction time.
Let us first calculate the anisotropic curvature on the surface of Tr,R,h . We will
split the surface Γ = ∂Tr,R,h into (o) the outer surface Γo where σ̃ ◦ (x′ ) = R, (i) the
inner surface Γi where σ̃ ◦ (x′ ) = r, and (s) the side facets Γs with |xn | = h.

SOLUTIONS OF THE CRYSTALLINE MEAN CURVATURE FLOW

7

2R
2r

2h

Figure 2. Torus Tr,R,h for n = 3 and σ̃(p′ ) = |p′ |1 (left), andwith hexagonal anisotropy (right).
To compute the anisotropic curvature, we need to find a Cahn-Hoﬀman vector
field z : Γ → Rn on the surface such that z(x) ∈ ∂σ(ν(x)) for Hn−1 -a.e. x withminimal ∥divΓ z∥L2 (Γ) . It can be shown that divΓ z minimizes the L2 norm exactly
when divΓ z is constant on the flat parts of Γ, see [11].
One such z can be given explicitly as
(
xn )
(2.4)
z(x) = g(σ̃ ◦ (x′ ))x′ ,,x = (x′ , xn ) ∈ Γ,hwhereg(s) := as−n+1 + b,
with
Rn−2 rn−2 (R + r)
Rn−2 + rn−2,
b = n−1
.
n−1
n−1
R
−r
R
− rn−1
The constants a and b are chosen so that g(r)r = −1 and g(R)R = 1. Since r < R
and g(s)s is increasing in s as a < 0, we conclude that −1 < g(s)s < 1 for r < s < R.
In particular,
a=−

(2.5)

σ̃ ◦ (g(σ̃ ◦ (x′ ))x′ ) ≤ 1 ⇔ g(σ̃ ◦ (x′ ))x′ ∈ ∂ σ̃(0)for all x = (x′ , xn ) ∈ Γ.

Furthermore,g ′ (s)s + (n − 1)g(s) = (n − 1)band thus(2.6)

divx′ (g(σ̃ ◦ (x′ ))x′ ) = g ′ (σ̃ ◦ (x′ ))x′ · ∇σ̃ ◦ (x′ ) + (n − 1)g(σ̃ ◦ (x′ ))= g ′ (σ̃ ◦ (x′ ))σ̃ ◦ (x′ ) + (n − 1)g(σ̃ ◦ (x′ )) = (n − 1)bwhenever ∇σ̃ ◦ (x′ ) exists, where we used that x′ · ∇σ̃ ◦ (x′ ) = σ̃ ◦ (x′ ), see [52].
Let us check that z is indeed a Cahn-Hoﬀman vector field on Γ. A convenient
approach is to write Γ as the level set of a Lipschitz continuous function. Consider(2.7)

ψ(x) = max(r − σ̃ ◦ (x′ ), σ̃ ◦ (x′ ) − R, |xn | − h),x = (x′ , xn ) ∈ Rn .

Clearly ψ is Lipschitz, Γ = {x : ψ(x) = 0} and Tr,R,h = {x : ψ(x) < 0}. ∇ψ is anouter normal vector of Γ (with respect to Tr,R,h ) and exists Hn−1 -a.e. on Γ. Since∂σ is positively zero-homogeneous, it is enough to show that z(x) ∈ ∂σ(∇ψ(x))
Hn−1 -a.e. x ∈ Γ.

8

N. POŽÁR

Let ri Γj denote the relative interior of the surface Γj . Thus suppose that x ∈ Γsuch that ∇ψ(x) exists. We have three cases:x ∈ ri Γo : In a neighborhood of such a point, we see that ψ(y) = σ̃(y ′ ) − R,and so in this case we have ∇ψ(x) = (∇σ̃ ◦ (x′ ), 0), from which we deduce
(
)
x′
xn
z(x) =,∈ ∂ σ̃(∇σ̃ ◦ (x′ )) × [−1, 1] = ∂σ(∇ψ(x)).
σ̃ ◦ (x′ ) hx ∈ ri Γi : This can be handled as the previous case, recalling that σ̃ is even′and thus − σ̃◦x(x′ ) ∈ ∂ σ̃(−∇σ̃ ◦ (x′ )).
x ∈ ri Γs : Now we are on the top or the bottom flat(facet, ψ(y)= |yn | − h in)
the neighborhood of this point, and so ∇ψ(x) = 0, |xxnn | . Thus, recalling(2.5), we deduce
(
)
{
}
xn
xn
z(x) = g(σ̃ ◦ (x′ ))x′ ,∈ ∂ σ̃(0) ×= ∂σ(∇ψ(x)).
|xn |
|xn |
We have proved that z defined in (2.4) is a Cahn-Hoﬀman field on Γ.
Now we show that divΓ z is constant on the flat parts of Γ, and hence it minimizes∥divΓ z∥L2 (Γ) . The surface (tangential) divergence divΓ z(x) for x ∈ Γ can becomputed easily using the level set method by constructing an extension z̄ of z
away from the facet such that z̄ ∈ ∂σ(∇ψ) a.e. in a neighborhood of x ∈ Γ so thatdiv z̄(x) exists. Then we have divΓ z(x) = div z̄(x), see [26].
Following the above verification that z is a Cahn-Hoﬀman vector field, we canconstruct the extension z̄ in a neighborhood of x ∈ Γ in the following three casesin the following way:)( ′x ∈ ri Γo : Here we can simply take z̄(y) = σ̃◦y(y′ ) , yhn which yields
n−2
1
+ ,
R
h
if the divergence exists, where we followed the computation in (2.6) withg(s) = 1s .
x ∈ ri Γi : Similarly to the previous situation, we may take
)
(
yn
y′
,
z̄(y) = − ◦ ′ ,σ̃ (y ) hwhich yields
n−2
1divΓ z(x) = div z̄(x) = −+rhif the divergence exists.
)(x ∈ ri Γs : Here we may take z̄(y) = g(σ̃ ◦ (y ′ ))y ′ , |yynn | , yielding
divΓ z(x) = div z̄(x) =divΓ z(x) = div z̄(x) = (n − 1)busing (2.6).
Since divΓ z is constant on the flat parts of the surface, we deduce that it minimizes
the L2 norm among all Cahn-Hoﬀman vector fields and therefore we can defineκσ = − divΓ z.
Let us now find the normal velocity V of the surface Γt = ∂Tr(t),R(t),h(t) wherer(t), R(t) and h(t) are some C 1 functions satisfying 0 < r(t) < R(t), 0 < h(t). It is
straightforward to use the level set method. Let ψ(x, t) be the function defined in

SOLUTIONS OF THE CRYSTALLINE MEAN CURVATURE FLOW

9(2.7) with the time dependent parameters r(t), R(t), h(t) above. Let us fix t andx ∈ Γt so that ψt (x, t) and ∇ψ(x, t) are defined. We can again consider three cases
ψtand use the formula V (x, t) = − |∇ψ|(x, t) [26]: R′ (t)x ∈ ri Γo ,′ )| , |∇σ̃◦ (xr ′ (t)(2.8)
V (x, t) = − |∇σ̃◦ (x′ )| , x ∈ ri Γi , ′h (t),x ∈ ri Γs .
To relate V to κσ using the law (2.3), it is left to compute β(ν) on Γ. Using∇ψ
the level set method again, ν = |∇ψ|
whenever ∇ψ exists and is nonzero [26], andtherefore(∇σ̃ ◦ (x′ ),0)x ∈ ri Γo , |∇σ̃◦ (x′ )| ,(2.9)

ν(x) =◦′(x ),0)− (∇σ̃
|∇σ̃ ◦ (x′ )| , x ∈ ri Γi ,(0, xn ),x ∈ ri Γs .
|xn |Recalling β given in (2.2) and that σ̃(±∇σ̃ ◦ (x′ )) = 1 whenever ∇σ̃ ◦ (x′ ) exists [52],we deduce{1x ∈ ri Γo ∪ ri Γi ,◦′ ,(2.10)β(ν(x)) = |∇σ̃ (x )|µ,x ∈ ri Γs .
Relating the anisotropic curvature κσ = − divΓ z with (2.8), (2.9) and (2.10) via(2.3), we must have
n−2
1
− ,
R
h
n−2
1
r′ = −+ ,(2.11)rh
Rn−2 + rn−2h′ = −(n − 1)µ n−1.
R− rn−1For the evolution to be self-similar, we therefore need( )′( r )′h=0=.
RRR′ = −

Thus there are constants γ ∈ (0, 1) and λ > 0 such thatr = γR,h = λR.

Plugging this into (2.11) and multiplying each equation by R, we get(2.12a)(2.12b)(2.12c)1
λn−2
1=− 2 +γ
λγ(n − 1)µ 1 + γ n−2·.
=−
λ1 − γ n−1

RR′ = −(n − 2) −From (2.12a) and (2.12b), we obtain(2.13)

λ=γ.
(n − 2)(1 − γ)

10N. POŽÁR

On the other hand, (2.12a) and (2.12c) yield()
1
1 + γ n−2(2.14)
λ=(n − 1)µ−1 .
n−2
1 − γ n−1
Hence combining (2.13) and (2.14) and multiplying both sides by (n − 2)(1 − γ n−1 ),we have
1 + γ + · · · + γ n−2 = (n − 1)µ(1 + γ n−2 ).
This equation has a solution γ ∈ (0, 1) only for certain µ. Let us only considern = 3. Then we have
0 = (2µ − 1)(1 + γ),
which has a solution γ ∈ (0, 1) only if µ = 21 . In this case there are infinitely manysolutions, for instance,
γ = 12 ,λ = 1.

Then the solution of (2.3) with initial data E0 = T 21 ,1,1 is the self-similar evolution{√
1 − 4tE0 , 0 ≤ t < t∗ = 41 ,Et =∅,t∗ ≤ t.
It is possible to show that this is the unique (open) evolution given by theviscosity solution of (2.3) in the sense of [34], but this is beyond the scope of thisnote. For more details, see [28].
2.3. The sponge. In this section we construct a shrinking “sponge” in dimension
n = 3. It is a self similar solution of the crystalline mean curvature flow with cubic
anisotropy σ(p) = |p|1 with genus 5 and resembles the level 1 Menger sponge. Its
shape is given by a cube centered at the origin from which the neighborhood (inthe maximum (ℓ∞ ) norm) of all three axes has been removed, see Figure 3. It is ananalogue of the solution that was discovered numerically by D. L. Chopp [21] forthe standard isotropic mean curvature flow.
Let us thus consider the crystalline mean curvature flow
(2.15)

V = σ(ν)κσ

where σ(p) = |p|1 is the cubic (ℓ1 ) anisotropy, with initial data
(2.16)

E0 := Sr0 ,R0

for some constants 0 < r0 < R0 , where{}
Sr,R := x ∈ R3 : |x|∞ < R, |xi | > r for at least two of i = 1, 2, 3
for 0 < r < R, see Figure 3.
We claim that the open evolution {Et }t≥0 solving (2.15) is of the form
Et := Sr(t),R(t) ,where r(t), R(t) are solutions of a certain ODE system. Our goal is to find 0 <r0 < R0 so that the evolution {Et }t≥0 is self-similar.
For 0 < r < R, due to the symmetries, Sr,R has only exactly two types of facets:(o) the outer facet given by a square of side 2R with square of side 2r with the same
center removed, and (i) the inner facet given by a rectangle with sides of lengths

SOLUTIONS OF THE CRYSTALLINE MEAN CURVATURE FLOW

(a)

(b)

(c)

11

(d)

Figure 3. The shrinking sponge. (A) The initial set E0 . The
outer edge has length 2R0 while the inner edge has length 2r0 .
(B) Self-similar sponge solution with initial size 0.8 at time t =0.015 computed with domain resolution M = 256. A numericalsolution uses Chambolle’s algorithm; since the solution is unstable,
the accumulated numerical errors will cause it to eventually diverge
from the self-similar solution. (C) Evolution for Rr00 > ξ ∗ . (D)Evolution for Rr00 < ξ ∗ .
2r and R − r. The crystalline curvature on each of these facets can be expressed as
the ratio of their signed perimeter and their area, that is,κo = −

8(R + r)2=−,2
2
4R − 4r
R−r

κi =

1
2
R − 3r−=.
r
R−rr(R − r)For a detailed explanation of why these are indeed the curvatures, see Section 2.2
for an explicit Cahn-Hoﬀman field construction, or [43].
Following the reasoning in Section 2.2, the normal velocity of the outer facet of{Et } is V = R′ (t) and the inner facet is V = −r′ (t). Note also that since the facetsof Sr,R are axes-aligned, we have σ(ν) = |ν|1 = 1 a.e. on ∂Sr,R . Therefore (2.15)
simplifies for {Et } into the system of ODEs

2

′

,
 R = κo = −R−r(2.17)
t > 0,
R − 3r, r′ = −κi = −r(R − r)with initial data r(0) = r0 , R(0) = R0 .
Let us perform a more detailed analysis of the behavior of the system (2.17).
The evolution will be self-similar if
R(t)≡ const.
ξ(t) :=r(t)Rewriting the system (2.17) for ξ and r, we have
ξ 2 − 3ξ − 2′,ξ = 2r (ξ − 1)(2.18)
t > 0,
ξ−3
 r′ = −,r(ξ − 1)
with initial data ξ0 = Rr00 > 1, r(0) = r0 > 0.

12

N. POŽÁR

We setg(ξ) := ξ 2 − 3ξ − 2.
√

3+ 17
∗
The constant ξ ∗ :=is the unique2{} value > 1 such that g(ξ ) = 0. Let us also
3
introduce QR := x ∈ R : |x|∞ < R
We can classify the behavior of the solution of (2.18) according to the initialratio ξ0 = Rr00 as follows:(a) 1 < ξ0 < ξ ∗ : Since ξ ′ < 0 and r(t) ≤ R0 , there exists a time t∗ > 0such that limt→t∗ − ξ(t) = 1, that is, limt→t∗ − R(t) − r(t) = 0. The sponge
gets thinner with time and converges to the edges of the cube QR(t∗ ) , andvanishes at t∗ ; Figure 3(D).
√(b) ξ0 = ξ ∗ : The solution is ξ(t) = ξ0 = ξ ∗ , r(t) = r02 − 2ct for 0 < t < t∗ ,
∗

r2

∗
∗
0
where c = ξξ∗ −3−1 > 0 and t = 2c > 0. At t = t the sponge vanishes. The
evolution is self-similar.
(c) ξ0 > ξ ∗ : There is time t□ > 0 such that limt→t□ − (ξ(t), r(t), R(t)) =(+∞, 0, R(t□ )). In this case, the holes close up and the sponge becomesa cube at time t□ ; Figure 3(D). After t = t□ , R evolves with R′ = − R2 . The
cube vanishes at a later time t∗ > t□ .

See Figure 3 for a numerical solution using Chambolle’s algorithm discussed in
Section 3. Note that (b) is unstable since g ′ (ξ ∗ ) > 0, and therefore numerically wewill observe either (a) or (c).
Let us summarize the solution of (2.15) in terms of the values 0 < r0 < R0 :
1 < Rr00 ≤ ξ ∗ :R0∗r0 > ξ :{Sr(t),R(t) , 0 ≤ t < t∗ ,Et :=∅,t ≥ t∗ .
□Sr(t),R(t) , 0 ≤ t < t ,□Et := QR(t) ,t ≤ t < t∗ ,∅,t ≥ t∗ .

As in the shrinking doughnut case by following [28], one can show that this is
the unique (open) evolution given by the viscosity solution of (2.15) in the sense of[34].
3. The numerical algorithm
An eﬃcient method for the mean curvature flow (1.2) is based on a minimizing
movement formulation due to A. Chambolle [17], which can be eﬃciently solved bya split Bregman iteration proposed by [35, 48]. Suppose that Ω ⊂ Rn is a bounded
convex domain and that the evolving set is contained in Ω. Furthermore, we needto assume that the one-homogeneous extension of β as in (1.1) is convex.
The insight of Chambolle is to formulate the minimizing scheme of F. Almgren,
J. E. Taylor and L. Wang [2] in terms of the signed distance function, so that theSOLUTIONS OF THE CRYSTALLINE MEAN CURVATURE FLOW

13

evolving set is its level set. In [17], he proposed the time discretization by theminimization problem (f ≡ 0 in [17])
)
(
∫
1
2
(3.1)
∥v − wm ∥ +σ(Dv) dx − ⟨fm , v⟩ ,vm+1 ← arg min
2hv∈L2 (Ω)Ω
where h > 0 is a chosen time step, wm is the signed distance function of the levelset {vm < 0} at the previous time step m, induced by the metric given by themobility β and fm = f (·, mh). The minimization is performed over all v ∈ L2 (Ω),
and ∥·∥ and ⟨·, ·⟩∫ are the L2 (Ω)-norm and inner product, respectively. The total
variation energy Ω ∫σ(Dv) dx is defined as the L2 lower semicontinuous envelope ofthe functional v 7→ Ω σ(∇v) dx defined for v in the Sobolev space W 1,1 (Ω). Notethat the functional in (3.1) is the Moreau-Yosida regularization of the total variationenergy with parameter h, and the minimization problem (3.1) is equivalent to theresolvent problem for the total variation energy. In other words, it is the implicit
Euler discretization of the total variation flow [3]. Since |∇vm+1 | ≈ |∇wm | = 1
when h is small, we can deduce that vm+1h−wm ≈ V using (2.1).
The full minimizing movements algorithm to find a discrete sequence of approximations Em , m = 0, 1, 2, . . . of the evolving set {Et }t≥0 at time steps tm = mhreads: Set E0 as the initial data and then iteratively for m = 0, . . . do

(3.2)

wm ← signdistβ Em)(µ2∥v − wm ∥ + ∥σ(∇v)∥1 − ⟨fm , v⟩ ,vm+1 ← arg minv∈L2 (Ω) 2
Em+1 ← {vm+1 < 0},where µ = h1 , and f is a given source, and we write ∥σ(∇v)∥1 =∫σ(Dv) dx. NoteΩthat the minimization is equivalent to the minimization ofµ2v − (wm + µ−1 fm ) + ∥σ(∇v)∥1 .
2
The signed distance function must correspond to the anisotropy β. Recall that
we assume that β is one-homogeneously extended to Rn as in (1.1), and such a
extension is an anisotropy. In particular, β is assumed to be convex (but it needs
not to be symmetric with respect to the origin). We define the signed distance assigndistβ Em (x) := inf β ◦ (x − y) − infc β ◦ (y − x)
y∈Em
y∈Em
{
inf y∈Em β ◦ (x − y),x∈/ Em ,=◦c β (y − x),− inf y∈Em
x ∈ Em ,where β ◦ (x) := sup {x · p : β(p) ≤ 1} is the convex polar of β. If β(p) = |p|,signdistβ is just the standard signed distance function induced by the Euclidean
metric.
Let us motivate the above choice of wm . If we set wm to be the signed distanceabove, we have β(∇wm ) = 1 a.e., see [52]. Then performing (3.2), vm+1 − wm is
approximately (κ − fm )h and the free boundary advances in the normal direction
by
(
)
∇wm
(κ − fm )h= −β(κ − fm )h = β(ν)(−κ + fm )h−|∇wm ||∇wm |
yielding the correct free boundary velocity.

14

N. POŽÁR

As h → 0, the evolution will converge to a continuous evolution {Et }t≥0 , see[17]. Provided that there is no fattening, this evolution will be the unique solutionof the anisotropic mean curvature flow [18, 19, 39] and, in the crystalline case inparticular, the unique viscosity solution solution of the crystalline mean curvature
flow [33, 34].
It might seem that the minimization problem in (3.2) is rather diﬃcult for numerical computation, mainly due to the non-diﬀerentiable second term. In particular,
the standard minimization methods like conjugate gradients or Newton iteration
are poorly suited. For this reason, Chambolle proposed an iterative algorithm in[17]. More recently, it was recognized in [48] that the minimization problem canbe addressed by the so-called proximal algorithms [50], the alternative direction
method of multipliers (ADMM) or the split Bregman method [35]. To find the2minimizer v of µ2 ∥v − u∥ + ∥σ(∇v)∥1 , we choose λ > 0, set b0 = d0 = 0 and theniterate for k = 0, 1, . . .
(3.3a)(3.3b)(3.3c)µλ22∥v − u∥ + ∥dk − ∇v − bk ∥ ,
2
2
v
λ
2
dk+1 ← arg min ∥σ(d)∥1 + ∥d − ∇vk+1 − bk ∥ ,2
dbk+1 ← bk + ∇vk+1 − dk+1 ,vk+1 ← arg min

until some stopping condition is reached, typically when ∥vk+1 − vk ∥2 is suﬃciently
small. Heuristically, this scheme introduces a new gradient variable d, and thenenforces the constraint d = ∇v by a quadratic penalty. Since we are minimizing asum of convex terms over two variables, the problems can be decoupled into iterating(3.3a) and (3.3b). (3.3c) is called a Bregman iteration, and it helps to enforce theconstraint exactly. Note that when convergence is achieved, (3.3c) implies d = ∇v.
For a detailed discussion of the motivation, convergence and other properties, see[35].
The advantage of this iteration process is the simplicity of the subproblems. The
first minimization problem (3.3a) is equivalent to finding the solution v of
(3.4)(µ − λ∆)v = µu + λ div(bk − dk )in Ω,with an appropriate boundary condition, for instance Neumann. It is not necessary
to solve it accurately, so one or two Gauss-Seidel iterations are suﬃcient [35].
To find a numerical solution, we use the finite diﬀerence method (FDM) or thefinite element method (FEM) to discretize (3.2), see Section 3.1. In both cases, we
represent d, b by discrete values on the lattice (FDM) or on the elements (FEM).
This has the important consequence that the minimization problem (3.3b) for d in
the discrete case completely decouples. Then, for each node or element i the ith
component of the minimizer dk+1,i is given by the so-called shrink operator [35](3.5)

dk+1,i = shrinkσ ((∇vk+1 + bk )i , 1/λ).

Note that the shrink operator can be expressed using the orthogonal projection on
the Wulﬀ shape W of σ [48],shrinkσ (ξ, 1/λ) := (I − PW/λ )(ξ).
In typical cases of isotropic, cubic and hexagonal anisotropies, that is, when the
Wulﬀ shape W is a sphere, a cube or a hexagonal prism, respectively, the orthogonal

SOLUTIONS OF THE CRYSTALLINE MEAN CURVATURE FLOW

15

Figure 4. Rounding of edges in the computation in Figure 18;
neighborhood of the top left vertex (top row) and bottom rightvertex (bottom row) are magnified: FDM (left column) vs. FEM(right column).
projection is very simple. More general Wulﬀ shapes can be handled by the method
proposed in [48].
It is interesting to relate the Bregman iterate b to the Cahn-Hoﬀman vector fieldin the definition of the anisotropic (crystalline) mean curvature. In particular, for anode or an element i, if convergence is achieved, (3.5) yields di = shrinkσ (di +bi , 1/λ)and we have bi = PW/λ (di + bi ). From the last equality, we see that either di = 0,
but then λbi ∈ W, or di ̸= 0, but then λbi ∈ ∂W and di is a normal of ∂W at λbi .
In other words, λbi ∈ ∂σ(di ). From (3.4) we deduce that v−u
h = λ div b. Hence λb
is a discrete Cahn-Hoﬀman vector field for the resolvent problem (3.1).
3.1. Discretization: FDM vs FEM. The standard way [35, 48] to discretize(3.3) is to use the finite diﬀerence method (FDM). This is quite straightforward,
see also [17], and seems suﬃcient in two dimensions and for applications in image
processing. For instance, (3.4) is discretized using the standard central diﬀerence
scheme on a 2n + 1 point stencil. However, in three dimensions this approach
introduces unwanted artifacts, for instance rounding of edges in certain directions,
see Figure 4. This is most likely caused by the fact that the gradient d and theCahn-Hoﬀman vector field b on a given cubic cell are given only by four out of theeight nodes. This restriction of the number of degrees of freedom for the gradient
limits the ability of the discretization to account for a rapidly changing vector fieldnear an edge.
We propose a finite element method (FEM) discretization. We minimize thefunctional in (3.3a) on the space of piece-wise linear functions (P 1 elements) on atetrahedral mesh. The vector fields ∇v, d and b are then approximated by piecewise constant vector fields (P 0 elements) on this mesh. The minimization problem(3.3b) on the space of piece-wise constant (P 0 elements) completely decouples oneach element and the minimizer is given by (3.5).
The merit of this approach is the six-fold increase of the number of degrees of
freedom of d and b (for the cost of also increasing the memory and computationaltime requirements), which allows for a better resolution of edges, see Figure 4.
Additionally, it is consistent with our redistance approach in Section 3.2. However,

16

N. POŽÁR

Figure 5. A triangular anisotropy demonstrates a topologicalchange in 2D for non-even anisotropies. Left: FDM Right: FEM.
Parameters h = 10−4 , M = 64, and plot step 0.005 = 50h.
(0, 1) = ξ2,1

(1, 1) = ξ1,3 = ξ2,3
T2

T1
(0, 0) = ξ1,1 = ξ2,1

(1, 0) = ξ1,2

Figure 6. Tessellation of the unit cube in n = 2.
at some extreme cases in two dimensions the FEM performs worse than FDM, see
Figure 5 for an example.
Let us present the mesh construction in an arbitrary dimension to fix the idea.
We assume that the computational domain is a cube Ω is subdivided into M nsmaller cubes of equal size, M ∈ N is the resolution. We split each of the smaller
cubes into n! simplices in the following way. By translation and scaling, we mayn!
assume that the small cube is Q = [0, 1]n . Let {pm }m=1 be the n! permutationsof (1, . . . , n). For each m = 1, . . . , n!, we define the sequence of n + 1 vertices
ξm,1 , . . . , ξm,n+1 of Q usingξm,1 = (0, . . . , 0),{(ξm,j )i + 1(ξm,j+1 )i =(ξm,j )i

if i = (pm )j ,otherwise,i = 1, . . . , n.

Then their convex hulls,conv(ξm,1 , . . . , ξm,n+1 ),1 ≤ m ≤ n!,are simplices that form a tessellation of Q, see Figure 6.
We use the mesh that is given by the above tessellation of each of the M n cubesin Ω. The total number of elements is n!M n .
3.2. Redistance. At each time step of the algorithm (3.2), the signed distancefunction signdistβ {vm < 0} to the 0-level set of the solution of the minimization
problem has to be recomputed. (This is not strictly true, see the discussion in
Section 3.3.) It is sometimes referred to as redistance. The problem amounts to

SOLUTIONS OF THE CRYSTALLINE MEAN CURVATURE FLOW

17

solving the (anisotropic) eikonal equation |∇w|β := β ◦ (∇w) = 1 with boundary
data w = 0 on ∂E. There are various eﬃcient methods for doing this, including
an iteration scheme [55], as well as more direct algorithms like the fast marching
method [53] or the fast sweeping method [58]. We choose the fast sweeping method
due to its simplicity and eﬃciency.
The boundary ∂E is given as the 0-level set of a function v with discrete values
on a regular grid. Unfortunately, in general the 0-level set of the resulting approximation w of the signed distance function will be diﬀerent from the 0-level set ofthe original function v. The fast marching and fast sweeping methods require an
initialization step, where the distance function is assigned at grid points that are
direct neighbors of the 0-level set of v. A special care must be taken so that theinterface is not moved unnecessarily, as these eﬀects might quickly accumulate overa series of consecutive time steps. This is especially important at points where
the surface should not move, a typical case for non-convex and non-concave facets,
where this unwanted redistance eﬀect might dominate the evolution.
There are a few standard schemes to initialize the nearby values in the literature,by analyzing the intersection of the level set with the grid lines [1, 17]. However,they do not seem to reproduce the correct value of the distance function even if
the level set is flat, which is a common situation in the crystalline mean curvature flow. Furthermore, the generalization to three dimensions seems unnecessarily
complicated.
We choose a naive method that appears to be superior in our case, is very simple
to implement in an arbitrary dimension, and that computes the exact signed distance function in a neighborhood of the flat facets (away from vertices and edges).
The idea is to split the M n cubes of the uniform grid into n! simplices as explainedin the construction of the mesh for the finite element method in Section 3.1 and suppose that v is aﬃne on each of them. This is very much in the spirit of the marchingtetrahedra method [22]. Then for all elements that the 0-level set intersects, thatis, on which v changes sign, we set the initial value of the signed distance function
at each vertex of the element to be the value of v normalized by the β-norm of thegradient of the aﬃne function given by v on the element. If a given grid node isa vertex of multiple elements intersecting the 0-level set, we set the initial value
to be the minimum over all the elements. To be more explicit, suppose that vi ,wi , i = 1, . . . , N are the values at the grid nodes xi , and Tj , j = 1, . . . , K are thesimplices on which v changes sign. We initialize wi to
wi = (sign vi ) inf|vi |1≤j≤K β ◦ ( ∇v|T )jxi ∈∂Tj,

where the infimum is defined as +∞ if it is over an empty set and ∇v|Tj is the(constant) value of ∇v on the simplex Tj . This initialization method is second
order accurate, O(M −2 ), near a smooth surface, in contrast to the first order accuracy of the initialization in [1, 17]. Moreover, the unwanted artifacts caused by
the redistance seem to be reduced, for instance compare the preservation of nonconvex/non-concave facets in Figure 17, even for a modest space resolution M = 64.
After this initialization step, we perform the 2n sweeps of the fast sweeping
method [58].

18

N. POŽÁR

3.3. Notes on the implementation. Let us give a few notes on our implementation from a practical point of view.
In the discussion of the iteration (3.3), we suggested to initialize b0 = d0 = 0.
However, this is unnecessary since the iteration converges to the unique minimizer
no matter the initial guess. Since the vector fields b and d vary relatively slowly from
one time step to another, we can reuse the value of b and d from the previous time
step to start the iteration. The main consequence is that the number of iterationsof (3.3) necessary to obtain a reasonably accurate result is dramatically decreased
as the work is spread out over consecutive time steps, see Table 1.
The choice of the stopping criterion for the iteration (3.3) is important. In our
implementation we stop once ∥vk+1 − vk ∥ < εbtol for some given εbtol > 0, where∥·∥ is for example the ℓ2 norm. When εbtol is chosen too large, the facets do not
become flat as there are too few iterations to reduce the lowest frequency mode ofthe error, whose wavelength is proportional to the size of the facets. On the other
hand, if εbtol is chosen too small, the number of iterations becomes unnecessarily
high for no gain in accuracy, which is limited by the time step error O(h) andthe ability to resolve the corners and edges of the crystal with accuracy O(M −1 ).
Unfortunately, there does not seem to be an explicit way to estimate the necessary
εbtol and some experimentation is required, see Figure 8.
Another source of numerical problems is the computation of the distance function in (3.2), see the discussion in Section 3.2. However, often vm is a very good
approximation of wm close to the zero-level set and hence we can take wm = vm
and skip the distance computation for a few time steps, after which we need torecompute the distance function. This can help reducing the redistance artifacts
that sometimes appear, especially if the time step is very small and there are partsof the surface that do not move, such as non-convex/non-concave facets or curved
sections. We did not use this in the results presented here.
3.4. Numerical results. To illustrate the performance, we present a few simple
numerical results based on our implementation of the above algorithm in the Rust
programming language. The domain is always taken to be Ω = (− 12 , 12 )n andλ = µ8 in (3.3a), (3.3b), see Section 3.4.2. The stopping condition is chosen as1∥vk+1 − vk ∥ℓ2 < 10−5 for n = 2 and ∥vk+1 − vk ∥ℓ2 < 10−4 M 2 for n = 3, where∥·∥ℓ2 is the discrete ℓ2 -norm, so that the average stopping tolerance per element
scales like M −1 . This choice seems to give reasonable results across all computations
presented here. We use M = 64 in n = 2 and M = 256 for n = 3. The performancefor n = 3 is ≈ 1min/timestep for FDM and ≈ 5min/timestep for FEM on a single
core of Intel Core i7-4770K @ 3.50 GHz. The figures depicting the solution show
the zero-level set of the numerical solution vm .
3.4.1. Estimating the error of solutions. To quantify the error between the numerical and the exact solution, we consider the Hausdorﬀ distance between two surfaces
Γ1 , Γ2 ⊂ Rn in the maximum and L2 norms as()
distH,∞ (Γ1 , Γ2 ) := max max dist(x, Γ2 ), max dist(x, Γ1 ) ,
(∫

x∈Γ1

x∈Γ2

2

dist(x, Γ2 ) dH

distH,2 (Γ1 , Γ2 ) :=
Γ1

n−1

∫
2

dist(x, Γ1 ) dH

+
Γ2

n−1

) 21
.

SOLUTIONS OF THE CRYSTALLINE MEAN CURVATURE FLOW

circlerel. error

100

10 1

10 2

10 2

10 3

10 3

10 4

10 4

10 5

10 5
100

101

M

square

100

10 1

102

19

L
L2

100

101

M

102

Figure 7. Relative error of the numerical Hausdorﬀ distance fortwo concentric circles of radii 0.3 and 0.4 (left) and two concentricsquares of side-length 0.6 and 0.8 (right) as a function of the meshresolution M . The distance function in our implementation of thefast sweeping method is only first-order accurate, which influences
the computation for larger M . This is more apparent in the caseof the maximum distance. Dotted lines are M −1 and M −2 forcomparison.
Note that distH,∞ is the usual Hausdorﬀ distance. Both are defined to be ∞ ifexactly one of the surfaces is empty.
Numerically, these two distances between two level sets are computed with thediscrete distance function given by the fast sweeping method with the initialization
described in Section 3.2. The integrals and maxima are approximated by thoseof the piece-wise linear distance function over a mesh generated from the levelset function by the marching tetrahedra method [22] in two or three dimensions.
Both distances work well for smooth surfaces, but distH,∞ penalizes corners andoveremphasizes errors on a small part of the boundary. In Figure 7 we test theaccuracy of this numerical algorithm in cases of the level set of (x21 +x22 )1/2 −r (circle)
and max(|x1 |, |x2 |)−r (square) with r = 0.3, 0.4 for each. Initially, when M −1 ≈ 0.1,
the accuracy is dominated by the accuracy of the level set reconstruction using themarching tetrahedra and the distance function initialization explained in Section 3.2
and both of these are second-order accurate. Once M −1 ≪ 0.1, the first-order
discretization error in the fast sweeping algorithm dominates. Therefore we report
the errors using distH,2 , which appears to be less sensitive to this error.
3.4.2. Choice of λ. Let us address the choice of the parameter λ in proportion toµ. Theoretically, the split Bregman iteration will converge for any λ > 0. However,this parameter significantly aﬀects the speed of convergence. In the original paper[35] it was suggested to use λ = 2µ. But after somewhat extensive testing in ourcase it appears that the value λ = µ8 gives the fastest convergence and smallest error
overall over both two and three dimensional computations. In this test, we simply
ran the computation for many choices of parameters µλ and stopping conditions εbtol

20

N. POŽÁR

6
4
2
0

M = 512, h = 1 × 10 5
btol = 0.01 × 2 k

btol = 0.01 × 2 k

M = 64, h = 1 × 10 45
10
= 0.001 × 2k

15

6
4
2
0

5
10
= 0.001 × 2k

15

Figure 8. Performance for diﬀerent µλ and εbtol for the shrinking
hexagon in n = 2. The area of the discs is proportional to the logarithm of the error, while the area of the diamonds is proportional
to the logarithm of the computational time. µλ = 0.001 × 27 ≈ 18
is generally well-performing.
in ∥vk+1 − vk ∥ℓ2 < εbtol for a self-similar shrinking hexagon in two dimensions andshrinking cube in three dimensions with various M and h, and plotted the errors
and computational time. See Figure 8 for an example. Due to the computational
cost in three dimensions, the testing was much less exhaustive there.
3.4.3. Self-similar solutions. We test the method in two and three dimensions using
self-similar solutions of V = β(ν)κσ . In two dimensions, we use the shrinking Wulﬀ
shape solution with hexagonal symmetry: the Wulﬀ shape is the regular hexagon
with edge of length 1, Figure 9. The results are given for three diﬀerent time
step sizes h. For implementation reasons,this and the below tests with hexagonal√
mobility have time scaled by the factor 2 3 3 . Therefore the extinction time in n = 2√√220.42 3≈ 0.069, and 0.48 3 ≈ 0.035 in n = 3.
is not t∗ = 0.4
2 = 0.08 as expected but
4
In three dimensions, we use the shrinking Wulﬀ shape solution when the Wulﬀ shapeis the hexagonal prism with height 2 and regular hexagonal base with edges of length
1, Figure 10. Finally, we use the self-similar solutions constructed in Section 2. Dueto their instability, presence of diﬀerent types of facets and simplicity, they seem
to provide a convenient benchmark. For the crystalline shrinking doughnut, weuse cubic and hexagonal anisotropies σ̃, see Figure 11 and Figure 12 respectively.
The error of the sponge solution is presented in Figure 13. In all of these, theerror is estimated as maxm∈M distH,2 (∂E(tm ), ∂Em ) where M is a set of selected
time steps such that for all m ∈ M we have tm = 0.002z for some z ∈ Z, andthat tm are far enough from the extinction time t∗ , where all the solutions diverge
significantly due to the instability. Table 1 shows the details on the performance of
the algorithm for n = 3 and M = 64, 256.
3.4.4. Qualitative tests. Samples of evolutions in two dimensions for various anisotropies are shown in Figure 17. In this case, the stopping condition is taken as∥vk+1 − vk ∥ℓ2 < 10−4 to illustrate the eﬀect of choosing too large of a stopping
tolerance, and therefore the very long edges in the bottom right figure are notcompletely straight. The stationary parts of the boundary are very well preserved
and redistancing artifacts are not noticeable. In the rest of the figures, we show

SOLUTIONS OF THE CRYSTALLINE MEAN CURVATURE FLOW21

M = 64
FDM
FEM
Test#Breg/step time/Breg (s) #Breg/step time/Breg (s)Hex Wulﬀ
24
0.017
29
0.066
ℓ1 doughnut
35
0.017
52
0.068
Hex doughnut
38
0.017
39
0.075
Sponge
28
0.014
33
0.044
M = 256
FDM
FEM
Test
#Breg/step time/Breg (s) #Breg/step time/Breg (s)Hex Wulﬀ
42
0.93
42
2.8
ℓ1 doughnut
100
0.68
109
3.0
Hex doughnut
93
0.78
98
3.5
Sponge
133
0.69
125
3.1
Table 1. Performance of the numerical method in n = 3 with
M = 64 and M = 256 for tests in Section 3.4.3 based on the
average number of Bregman iteration per time step (#Breg/step)and average time per one Bregman iterations (time/Breg (s)) inseconds. Ran on Intel Xeon CPU E5-4650 v2 @ 2.40GHz.

qualitative tests of facet breaking, bending and topological changes in n = 3, see
Figure 14, Figure 15, Figure 16 and Figure 18.
3.4.5. Discussion. The implementation of the numerical method presented in thepaper is able to reproduce all of the tested features of the highly singular crystalline
mean curvature flow without the need to regularize the crystalline curvature: flat
faces, sharp edges and vertices, correct facet breaking and bending, and finally topological changes like neck pinching and facet-edge or edge-vertex collisions. Moreover,the method appears to be first order accurate both in space and time. We also presented a FEM discretization of the total variation energy that seems to perform
visually somewhat better than the FDM discretization for non-cubic anisotropies
for the cost of increased computational time. Our scheme for the reinitializationof the distance function seems to perform well without introducing redistancing
artifacts.
Various optimizations, such as performing the computation only in a small neighborhood of the level set to significantly reduce the computational complexity, as wellas the coupling of the curvature flow with the heat equation via the Gibbs-Thomson
relation are under investigation [51].
Acknowledgments. Parts of this paper are based on an extended abstract for atalk that the author gave at the 2017 Spring Meeting of the Mathematical Societyof Japan.

N. POŽÁR

maxdistH, 2

22h = 1.6 × 10 3
h = 4.0 × 10 4
h = 1.0 × 10 4

10 2

10 3
32

64

128
M

256

512

maxdistH, 2

Figure 9. The error of a self-similar shrinking hexagonal Wulﬀ
shape of side-length 0.4 in n = 2 with the numerical solution using
the FEM (solid) and FDM (dashed) discretization (maximum overt ∈ [0, 0.05]). Dotted line is M −1 .

FEM
FDM
M 1

10 1
10 2

distH, 2

8

16

32

M

64

128

0.015
t

0.020

0.025

256

FEM
FDM
0.01
0.00

0.000

0.005

0.010

0.030

Figure 10. The error of a self-similar shrinking hexagonal Wulﬀ
shape with a hexagonal base with edges of length 0.4 and height0.8 in n = 3 with h = 10−4 as the function of M (top, maximum
over t ∈ [0, 0.026]) and t for M = 128 (bottom).

maxdistH, 2

SOLUTIONS OF THE CRYSTALLINE MEAN CURVATURE FLOW

10 1
10 2
0.03

distH, 2

FEM
FDM
M 1

8

16

32

M

64

128

256

FEM
FDM

0.02
0.01
0.00

0.000 0.005 0.010 0.015 0.020 0.025 0.030 0.035t

maxdistH, 2

Figure 11. The error of a self-similar crystalline shrinking doughnut with ℓ1 anisotropy from Section 2.2 of side-length 0.8 in
n = 3 with h = 10−4 as the function of M (top, maximum overt ∈ [0, 0.02]) and t for M = 128 (bottom).

FEM
FDM
M 1

10 1
10 2
8

32

M

64

128

0.015
t

0.020

0.025

256

FEM
FDM

0.02
distH, 2

16

0.01
0.00

0.000

0.005

0.010

0.030

Figure 12. The error of a self-similar crystalline shrinking doughnut with hexagonal anisotropy from Section 2.2 of hexagonal basewith side-length 0.4 and height 0.8 in n = 3 with h = 10−4 as thefunction of M (top, maximum over t ∈ [0, 0.02]) and t for M = 128(bottom).

23

24

N. POŽÁR

FEM
FDM
M 1

maxdistH, 2

10 1
10 2
8

32

0.005

0.010

M

64

128

256

FEM
FDM

0.06
distH, 2

16

0.04
0.02
0.00

0.000

0.015
t

0.020

0.025

Figure 13. The error of a self-similar shrinking sponge from Section 2.3 with R0 = 0.4 in n = 3 with h = 10−4 as the function of
M (top, maximum over t ∈ [0, 0.014]) and t for M = 128 (bottom).

Figure 14. Facet breaking example in 3D [14] at three selected
times—L-shaped facets break into rectangular facets.

Figure 15. Facet bending (the top facet) in an axes-aligned initial
shape for the cubic (ℓ1 ) anisotropy [42], with a detail on the right.

SOLUTIONS OF THE CRYSTALLINE MEAN CURVATURE FLOW

Figure 16. Topological change in 3D with a hexagonal
anisotropy—a pinch-oﬀ of an initially connected dumbbell shape.
FDM discretization; some rounding of edges is apparent.

Figure 17. Flow V = κσ with various anisotropies: h = 10−4 ,M = 64, FEM discretization and plot step 0.005 = 50h. The
bounding box of the initial curve is the square [−0.4, 0.4]2 inside
the computational domain (− 12 , 12 )2 . Even at this relatively low
resolution, artifacts caused by the redistance are not apparent.

25

26

N. POŽÁR

Figure 18. Topological change in 3D with a triangularanisotropy—a pinch-oﬀ due to a collision of a facet and an edge.
FEM discretization.

References[1] D. Adalsteinsson and J. A. Sethian, The fast construction of extension velocities in level set
methods, J. Comput. Phys. 148 (1999), no. 1, 2–22, DOI 10.1006/jcph.1998.6090.
[2] F. Almgren, J. E. Taylor, and L. Wang, Curvature-driven flows: a variational approach, SIAM
J. Control Optim. 31 (1993), no. 2, 387–438, DOI 10.1137/0331020.
[3] F. Andreu-Vaillo, V. Caselles, and J. M. Mazón, Parabolic quasilinear equations minimizing
linear growth functionals, Progress in Mathematics, vol. 223, Birkhäuser Verlag, Basel, 2004.
MR2033382 (2005c:35002)[4] S. B. Angenent, Shrinking doughnuts, Nonlinear diﬀusion equations and their equilibrium
states, 3 (Gregynog, 1989), Progr. Nonlinear Diﬀerential Equations Appl., vol. 7, Birkhäuser
Boston, Boston, MA, 1992, pp. 21–38. MR1167827[5] S. Angenent and M. E. Gurtin, Multiphase thermomechanics with interfacial structure. II.
Evolution of an isothermal interface, Arch. Rational Mech. Anal. 108 (1989), no. 4, 323–391,
DOI 10.1007/BF01041068.
[6] J. W. Barrett, H. Garcke, and R. Nürnberg, A variational formulation of anisotropic geometric evolution equations in higher dimensions, Numer. Math. 109 (2008), no. 1, 1–44, DOI
10.1007/s00211-007-0135-5. MR2377611[7], Parametric approximation of surface clusters driven by isotropic and anisotropic
surface energies, Interfaces Free Bound. 12 (2010), no. 2, 187–234, DOI 10.4171/IFB/232.
MR2652017[8], Finite-element approximation of one-sided Stefan problems with anisotropic, approximately crystalline, Gibbs-Thomson law, Adv. Diﬀerential Equations 18 (2013), no. 3-4,
383–432.
[9], On the stable discretization of strongly anisotropic phase field models with applications to crystal growth, ZAMM Z. Angew. Math. Mech. 93 (2013), no. 10-11, 719–732, DOI
10.1002/zamm.201200147. MR3118772[10], Stable phase field approximations of anisotropic solidification, IMA J. Numer. Anal.
34 (2014), no. 4, 1289–1327, DOI 10.1093/imanum/drt044. MR3269427[11] G. Bellettini, An introduction to anisotropic and crystalline mean curvature flow, Hokkaido
Univ. Tech. Rep. Ser. in Math. 145 (2010), 102–162.
[12] G. Bellettini, V. Caselles, A. Chambolle, and M. Novaga, Crystalline mean curvature flow ofconvex sets, Arch. Ration. Mech. Anal. 179 (2006), no. 1, 109–152, DOI 10.1007/s00205-0050387-0.
[13] G. Bellettini and M. Novaga, Approximation and comparison for nonsmooth anisotropic
motion by mean curvature in RN , Math. Models Methods Appl. Sci. 10 (2000), no. 1, 1–10,
DOI 10.1142/S0218202500000021.
[14] G. Bellettini, M. Novaga, and M. Paolini, Facet-breaking for three-dimensional crystals evolving by mean curvature, Interfaces Free Bound. 1 (1999), no. 1, 39–55, DOI 10.4171/IFB/3.
[15], On a crystalline variational problem. I. First variation and global L∞ regularity,
Arch. Ration. Mech. Anal. 157 (2001), no. 3, 165–191, DOI 10.1007/s002050010127.
[16] H. Brézis, Monotonicity methods in Hilbert spaces and some applications to nonlinear partial diﬀerential equations, Contributions to nonlinear functional analysis (Proc. Sympos.,SOLUTIONS OF THE CRYSTALLINE MEAN CURVATURE FLOW

27

Math. Res. Center, Univ. Wisconsin, Madison, Wis., 1971), Academic Press, New York, 1971,
pp. 101–156.
[17] A. Chambolle, An algorithm for mean curvature motion, Interfaces Free Bound. 6 (2004),
no. 2, 195–218, DOI 10.4171/IFB/97.
[18] A. Chambolle, M. Morini, M. Novaga, and M. Ponsiglione, Existence and uniqueness foranisotropic and crystalline mean curvature flows, preprint, available at https://arxiv.org/
abs/1702.03094.
[19] A. Chambolle, M. Morini, and M. Ponsiglione, Existence and uniqueness for a crystalline mean curvature flow, Comm. Pure Appl. Math. 70 (2017), no. 6, 1084–1114, DOI
10.1002/cpa.21668. MR3639320[20] Y. G. Chen, Y. Giga, and S. Goto, Uniqueness and existence of viscosity solutions of generalized mean curvature flow equations, J. Diﬀerential Geom. 33 (1991), no. 3, 749–786.
[21] D. L. Chopp, Computation of self-similar solutions for mean curvature flow, Experiment.
Math. 3 (1994), no. 1, 1–15.
[22] A. Doi and A. Koide, An Eﬃcient Method of Triangulating Equi-Valued Surfaces by Using
Tetrahedral Cells, IEICE TRANSACTIONS on Information and Systems E74-D (1991), no. 1,
214–224.
[23] L. C. Evans and J. Spruck, Motion of level sets by mean curvature. I, J. Diﬀerential Geom.
33 (1991), no. 3, 635–681.
[24] M. Gage and R. S. Hamilton, The heat equation shrinking convex plane curves, J. Diﬀerential
Geom. 23 (1986), no. 1, 69–96. MR840401[25] H. Garcke, Curvature driven interface evolution, Jahresber. Dtsch. Math.-Ver. 115 (2013),
no. 2, 63–100, DOI 10.1365/s13291-013-0066-2.
[26] Y. Giga, Surface evolution equations - a level set approach, Monographs in Mathematics,
vol. 99, Birkhäuser Verlag, Basel, 2006. (earlier version: Lipschitz Lecture Notes 44, Universityof Bonn, 2002).
[27] M.-H. Giga and Y. Giga, Evolving graphs by singular weighted curvature, Arch. Rational
Mech. Anal. 141 (1998), no. 2, 117–198.
[28], Crystalline and level set flow—convergence of a crystalline algorithm for a generalanisotropic curvature flow in the plane, Free boundary problems: theory and applications, I(Chiba, 1999), GAKUTO Internat. Ser. Math. Sci. Appl., vol. 13, Gakkōtosho, Tokyo, 2000,
pp. 64–79. MR1793023[29], Generalized motion by nonlocal curvature in the plane, Arch. Ration. Mech. Anal.
159 (2001), no. 4, 295–333, DOI 10.1007/s002050100154.
[30] M.-H. Giga, Y. Giga, and H. Hontani, Self-similar expanding solutions in a sector for acrystalline flow, SIAM J. Math. Anal. 37 (2005), no. 4, 1207–1226, DOI 10.1137/040614372.
MR2192293[31] M.-H. Giga, Y. Giga, and P. Rybka, A comparison principle for singular diﬀusion equations
with spatially inhomogeneous driving force for graphs, Arch. Ration. Mech. Anal. 211 (2014),
no. 2, 419–453, DOI 10.1007/s00205-013-0676-y.
[32] Y. Giga, P. Górka, and P. Rybka, Bent rectangles as viscosity solutions over a circle, Nonlinear Anal. 125 (2015), 518–549, DOI 10.1016/j.na.2015.05.033.
[33] Y. Giga and N. Požár, A level set crystalline mean curvature flow of surfaces, Adv. Diﬀerential Equations 21 (2016), no. 7-8, 631–698.
[34], Approximation of general facets by admissible facets for anisotropic total variationenergies and its application to the crystalline mean curvature flow, Comm. Pure Appl. Math.
71 (2018), no. 7, 1461–1491, DOI doi.org/10.1002/cpa.21752.
[35] T. Goldstein and S. Osher, The split Bregman method for L1-regularized problems, SIAM J.
Imaging Sci. 2 (2009), no. 2, 323–343, DOI 10.1137/080725891.
[36] M. A. Grayson, The heat equation shrinks embedded plane curves to round points, J. Diﬀerential Geom. 26 (1987), no. 2, 285–314. MR906392[37] M. E. Gurtin, Thermomechanics of evolving phase boundaries in the plane, Oxford Mathematical Monographs, The Clarendon Press, Oxford University Press, New York, 1993.
[38] G. Huisken, Flow by mean curvature of convex surfaces into spheres, J. Diﬀerential Geom.
20 (1984), no. 1, 237–266. MR772132[39] K. Ishii, An approximation scheme for the anisotropic and nonlocal mean curvature
flow, NoDEA Nonlinear Diﬀerential Equations Appl. 21 (2014), no. 2, 219–252, DOI
10.1007/s00030-013-0244-z.

28

N. POŽÁR[40] T. Ishiwata, T. K. Ushijima, H. Yagisita, and S. Yazaki, Two examples of nonconvex selfsimilar solution curves for a crystalline curvature flow, Proc. Japan Acad. Ser. A Math. Sci.
80 (2004), no. 8, 151–154. MR2099341[41] Y. Kōmura, Nonlinear semi-groups in Hilbert space, J. Math. Soc. Japan 19 (1967), 493–507.
[42] M. Lasica, Personal communication (2016).
[43] M. Lasica, S. Moll, and P. B. Mucha, Total Variation Denoising in l1 Anisotropy, SIAM J.
Imaging Sci. 10 (2017), no. 4, 1691–1723. MR3709886[44] X. H. Nguyen, Construction of complete embedded self-similar surfaces under mean curvature
flow. I, Trans. Amer. Math. Soc. 361 (2009), no. 4, 1683–1701, DOI 10.1090/S0002-9947-0804748-X. MR2465812[45], Construction of complete embedded self-similar surfaces under mean curvature flow.
II, Adv. Diﬀerential Equations 15 (2010), no. 5-6, 503–530. MR2643233[46], Construction of complete embedded self-similar surfaces under mean curvature flow,
Part III, Duke Math. J. 163 (2014), no. 11, 2023–2056, DOI 10.1215/00127094-2795108.
MR3263027[47] M. Novaga and E. Paolini, Stability of crystalline evolutions, Math. Models Methods Appl.
Sci. 15 (2005), no. 6, 921–937, DOI 10.1142/S0218202505000571. MR2149929[48] A. Oberman, S. Osher, R. Takei, and R. Tsai, Numerical methods for anisotropic mean
curvature flow based on a discrete time variational formulation, Commun. Math. Sci. 9 (2011),
no. 3, 637–662, DOI 10.4310/CMS.2011.v9.n3.a1.
[49] S. Osher and J. A. Sethian, Fronts propagating with curvature-dependent speed: algorithms
based on Hamilton-Jacobi formulations, J. Comput. Phys. 79 (1988), no. no. 1, 12–49.
[50] N. Parikh and S. Boyd, Proximal Algorithms, Found. Trends Optim. 1 (2014), no. 3, 127–239,
DOI 10.1561/2400000003.
[51] N. Požár, A numerical level set method for the Stefan problem with a crystalline GibbsThomson law (preprint).
[52] R. T. Rockafellar, Convex analysis, Princeton Mathematical Series, No. 28, Princeton University Press, Princeton, N.J., 1970. MR0274683 (43 #445)[53] J. A. Sethian, A fast marching level set method for monotonically advancing fronts, Proc.
Nat. Acad. Sci. U.S.A. 93 (1996), no. 4, 1591–1595, DOI 10.1073/pnas.93.4.1591.
[54] A. Stancu, Asymptotic behavior of solutions to a crystalline flow, Hokkaido Math. J. 27(1998), no. 2, 303–320, DOI 10.14492/hokmj/1351001287. MR1637988[55] M. Sussman, P. Smereka, and S. Osher, A Level Set Approach for Computing Solutions to
Incompressible Two-Phase Flow, Journal of Computational Physics 114 (1994), no. 1, 146 159, DOI http://dx.doi.org/10.1006/jcph.1994.1155.
[56] J. E. Taylor, Existence and structure of solutions to a class of nonelliptic variational problems, Symposia Mathematica, Vol. XIV (Convegno di Teoria Geometrica dell’Integrazione e
Varietà Minimali, INDAM, Roma, Maggio 1973), Academic Press, London, 1974, pp. 499–508.
MR0420407[57], Constructions and conjectures in crystalline nondiﬀerential geometry, Diﬀerential
geometry, Proceedings of the Conference on Diﬀerential Geometry, Rio de Janeiro (B. Lawson
and K. Tanenblat, eds.), Pitman Monogr. Surveys Pure Appl. Math., vol. 52, Longman Sci.
Tech., Harlow, 1991, pp. 321–336, DOI 10.1111/j.1439-0388.1991.tb00191.x.
[58] H. Zhao, A fast sweeping method for eikonal equations, Math. Comp. 74 (2005), no. 250,
603–627, DOI 10.1090/S0025-5718-04-01678-3.
(N. Požár) Faculty of Mathematics and Physics, Institute of Science and Engineering,
Kanazawa University, Kakuma town, Kanazawa, Ishikawa 920-1192, Japan.
E-mail address: npozar@se.kanazawa-u.ac.jp