text:
# Discrete models of force chain networks

## 1. Introduction

Cohesionless granular materials exhibit a range of behavior that has fascinated physicists, applied mathematicians, and engineers from Coulomb, in the late $`18^{th}`$ century, to Schaeffer, in the late $`20^{th}`$ and early $`21^{st}`$. It is therefore somewhat surprising that perhaps the simplest question one can ask about a granular material remains unanswered even now. Given a box full of sand, if a marble is placed on top of the sand, how is the weight of the marble distributed on the bottom of the box? The problem is not just quantitative; we do not even know the qualitative form of the pattern of vertical force on the bottom generated by the marble. We do not know, for example, whether the vertical force is largest directly underneath the marble or has a maximum on a ring. 111A review of the many and varied approaches that have been taken to the problem of determining stress distributions in noncohesive granular materials is beyond the scope of the present paper. The reader is referred to References \[1-13\] and the citations therein for background material.

Experiments suggest that qualitative features of the response may depend on the way in which the sand was put into the box or on the geometric and frictional properties of the grains.. Experiments on 2D layers with depths up to approximately 12 grain diameters show both two-peak and single-peak responses, depending on the geometry of the grains and/or degree of disorder in the packing.  Experiments on 3D layers have probed the response in layers as deep as 300 grain diameters, and single peaks directly beneath the applied load were observed.  A theoretical framework that allows a unified interpretation of these results would be of great interest.

Let us formulate a thought experiment that focuses on the central features governing the stress pattern at the bottom of the sandbox. The fundamental problem has to do with how stress is transmitted through the bulk of the material, so it is natural to eliminate the side walls from consideration by imagining a slab of sand that is infinite in horizontal extent. It is also tempting to eliminate the effect of the bottom boundary by imagining it to be at infinity, where all stresses are infinitesimal, and calculating the stress as a function of depth far from the bottom boundary. We will see, however, that in the context of an interesting class of models the bottom boundary condition plays an essential role.

We can also eliminate from the problem the complicating effects of gravity. We consider a slab of sand between two plates that are being pushed together. The density of the sand is assumed to be high enough that the material rigidly resists the force applied by the plates. The material is then subjected to a localized force created by poking a needle through the top plate and applying a specified force to it. Note that both plates and the confining pressure are essential in this setup. Without them, some grains at the surfaces of the slab would likely be subject to outward forces that could not be balanced, so the material would not support any force at all on the needle.

A third simplification is to assume that the bulk granular material is isotropic. That is, we assume that the geometry of the grain packing does not distinguish any particular direction in space. While the local environment of any particular grain is clearly anisotropic, we assume that this anisotropy, as measured for example by the fabric tensor which characterizes the degree to which there are preferred directions for “bonds” joining the centers of grains in contact, vanishes when coarse-grained over sufficiently large volumes. This is a strong assumption, as construction history of a real granular material may well distinguish favored directions for contacts. Nevertheless, there are highly nontrivial structures to be uncovered even in isotropic models and it is surely appropriate to understand them before attempting to include the complications and additional parameters associated with intrinsic anisotropies.

Finally, for simplicity, we can work with a two-dimensional system. The top of the slab is a line defined to be at $`z=0`$ and the bottom is defined to be $`z=d`$. The slab is infinite in the horizontal $`x`$ direction. We are interested in determining the response to a localized force applied at $`x=0`$ and $`z=0`$.

Treating the sand as an ordinary (linear) elastic medium, a standard treatment  indicates that $`\sigma _{zz}(x,z)`$ should have a single peak centered on $`x=0`$ with a width that grows linearly with depth.222 In sufficiently strongly anisotropic elastic materials, it is possible to have two peaks, both of which have widths that grow linearly with depth. This treatment assumes, however, that the stress field is governed by equations derived using the continuity of a well defined displacement field and that there is a unique stress associated with each strain configuration. Neither of these is true for a granular material consisting of rigid particles that support frictional forces, so a new approach is needed.

The problem of deriving macroscopic stress equations from known microscopic or grain scale physics has proven quite difficult. An indication of just how difficult this might be can be found in almost any experiment or numerical simulation that generates images of the intensity distribution of stresses on scales of several grain diameters. In such images, one sees immediately that the stress is concentrated on filamentary structures called force chains that support compressive stress. These chains appear to be relatively straight on scales of up to 10 grain diameters or so and give the visual impression of splitting and fusing at a variety of angles. The presence of such structures suggests that passage from the grain scale to the macroscopic stress will involve two distinct steps: we must understand how grain scale physics favors the formation of chains, and also how the interactions among chains determine the macroscopic stress field. While these two tasks must ultimately be facets of a unified theory, it may be useful to approach them separately.

Models of directed force chain networks (DFCNs) have recently been proposed to bridge the gap between the scale of individual chains and the macroscopic stress, leaving open the issue of how grain scale physics promotes chain formation. (This type of model is also referred to as a Y  Y -model.). As detailed below, a “Boltzmann equation” governing the densities of chains with specified strengths and orientations can be obtained, assuming only that whether a chain splits at a given point is determined only by the local environment around that point and that environments that lead to splitting are homogeneously distributed throughout the material. An essential feature of this equation is that the unstressed solution is unstable; perturbations of it lead to divergent responses. Thus in order to calculate physically meaningful response functions, it is necessary to perturb around a pre-stressed state. The only case for which response functions have been calculated is the special one in which all force chains support the same magnitude of stress and lie on a 6-fold symmetric set of vectors. In that case, the response consists of two peaks whose centers diverge linearly with depth but whose widths grow like the square root of depth. In other words, the only solved case of response functions for DFCNs suggests that stress propagates along characteristic directions determined by the homogeneously pre-stressed state, in marked contrast to expectations from standard elasticity theory. This propagation of stress is a feature of noisy hyperbolic equations and is therefore called a “hyperbolic response.”

It is worth noting here that hyperbolic response is known to occur in models of isostatic packings of frictionless grains  and that simple constitutive relations leading to hyperbolic stress equilibrium equations have been shown to explain nontrivial features of experiments on sandpiles and granular columns.  One must also note, however, that direct measurements of the response seem to indicate that single-peaked response is associated with strongly disordered packings, and hyperbolic response occurs only in ordered systems or at shallow depths.

The question immediately arises as to whether the hyperbolic response of the 6-fold DFCN model is an artifact of the confinement to a discrete set of directions, or perhaps of the unique property that all chains in the 6-fold DFCN have the same strength. To investigate the robustness of the 6-fold DFCN results, one would like to find other analytically tractable cases. Here we address the first step in this direction – the identification of models for which stable solutions can be found for homogeneous macroscopic stresses.

In this paper we present homogeneous solutions to the DFCN Boltzmann equation for special cases in which chains lie along a discrete set of directions with 8-fold, or more generally 4N-fold, symmetry in two dimensions. Even the homogeneous solutions – the pre-stressed states about which we might contemplate calculating response functions – have nontrivial features. The actual calculation of the response functions these solutions is beyond our present scope.

This paper also provides corrections to the presentation fo the general Boltzmann equation of Ref. . The corrections do not affect any of the calculations on discrete models reported in that paper, but do clarify some conceptual points crucial for future work on continuum models or models permitting multiple splitting angles.

## 2. The DFCN model and Boltzmann equation

### 2.1. Definitions

A DFCN is a collection of line segments with associated compressive forces assigned such that force balance is achieved at every vertex. Figure 1 shows a simple example that illustrates several features. First, each segment is assigned a strength, or force intensity, that is indicated in the figure by the thickness of the line. Every chain is assumed to be under compression, so the force exerted on a vertex due to a given chain is directed along the chain toward the vertex, and the vector sum of the forces at any given vertex must vanish.

Second, each segment is assigned a direction, indicated by the arrow, which distinguishes its “beginning” from its “end.” The physical distinction lies in the local grain configurations at either end of the chain. Certain configurations do not allow propagation of a chain through them in certain orientations and therefore give rise to force chain splitting. The relevant configurations of this type for the network shown in the figure are indicated with large circles. Other local configurations would not necessarily require a splitting event, but do permit two chains to fuse when they intersect. The relevant configurations of this type are indicated with small circles.

The directions of chains in the entire network are determined by the specification of the applied boundary forces. One must think, for now, of the force being applied at point A as coming from a needle being poked through the top plate and held at constant force. This ensures that a chain beginning at A must be present. At point B, a local configuration exists that requires the chain to end. To balance the force at B, two chains must be initiated. Hence the arrows on these to chains must point away from B.

By tracing the chains initiated by applied forces at the boundary and taking into account the type of local environments at the chain endpoints, each chain direction can be uniquely assigned, with the exception of rare cases in which fusions happen to occur precisely at places that might be mistaken for splittings. When all directions are assigned, one is likely to find some chains that end, rather than begin, on a boundary. The boundary forces required to balance these forces, shown as dashed arrows in Figure 1 must be interpreted as a response to the applied forces.

Finally, there is the possibility that a single chain can appear to have inconsistent requirements on its direction, as occurs for the chain passing through point E in Figure 1. Here two splitting events have initiated chains along the same direction that meet head on. Formally, this corresponds to a fusion event in which the outgoing chain has zero strength, which thus appears as the annihilation of two chains, and the circle marking this fusion can be imagined to lie anywhere along the chain. The probability of such configurations is determined by the grain size, or, more precisely, by the ratio of the grain size to the typical separation between force chains of opposite directions. In the following, we will take this probability to be zero. (But keep in mind that a nonzero probability can be invoked to cut off the divergence in the density of weak force chains that arises in the 8-fold model.)

Figure 2 further illustrates the difference

between chain splitting and chain fusion in the context of a hexagonal array of grains. In panel (a), two forces are applied at the top boundary and the forces propagate along chains, crossing at the central grain. In panel (b), a similar situation is shown, but this time one of the chains is assumed to be specified by fixing its position at the bottom boundary instead of the top. In this case, a fusion of the two chains at the central grain is permitted (though not required) and the resulting configuration may be different. Panel (c) shows a local configuration (a missing disk) that requires the splitting of an incoming chain.

Thus each chain is characterized by an intensity, $`f`$, and a direction $`\widehat{n}`$, or, equivalently, a vector $`\stackrel{}{f}`$. Force balance at a vertex requires that $`\stackrel{}{f}`$ for incoming chains equals $`\stackrel{}{f}`$ for outgoing chains. Positive values of $`f`$ correspond to compressive stress along a chain. negative values to tensile stress. To model noncohesive materials, which do not support tensile stress, we take $`f`$ to be always positive definite. Thus there is no ambiguity introduced by using $`\stackrel{}{f}`$ to denote the pair $`(f,\widehat{n})`$. When convenient, we will use the angular variable $`\theta `$ to indicate the direction of $`\widehat{n}`$. We will use the term “$`\stackrel{}{f}`$-chain” to refer to a chain with the given strength and direction.

### 2.2. The Boltzmann equation

The following discussion supersedes the treatment given in Ref. . In that paper, mistaken reasoning was used in writing the continuum equation. The mistake does not affect the bulk of that paper, since the forms obtained there for networks of chains confined to a discrete set of directions were the same as those derived below. Nevertheless, the corrections made here may be important for future work on continuum models and are necessary for conceptual clarity. Note especially that the definition of $`P`$ has different dimensions here from that used in Ref. . The theory is developed here explicitly for the case of two dimensions. Generalization to higher dimensions is straightforward.

Let $`P(f,\theta ,\stackrel{}{r})`$ represent the density (in an ensemble average) of force chains of intensity $`f`$ and direction $`\theta `$ passing through the spatial point $`\stackrel{}{r}`$. In other words,

$$_f^{f+\delta f}f^{}𝑑f^{}_\theta ^{\theta +\delta \theta }𝑑\theta ^{}|\widehat{n}\widehat{u}|P(f^{},\theta ^{},\stackrel{}{r})$$
(1)
is the number of chains with intensity between $`f`$ and $`f+\delta f`$ that cross a unit length line segment passing through $`\stackrel{}{r}`$ and perpendicular to $`\widehat{u}`$. With this definition, $`P(\stackrel{}{f})`$ is defined with respect to a uniform measure in the 2D space of possible forces. Given $`P(\stackrel{}{f},\stackrel{}{r})`$, the local macroscopic stress tensor is given by

$$\sigma _{\alpha \beta }=_0^{\mathrm{}}f𝑑f_\pi ^\pi 𝑑\theta n_\alpha n_\beta fP(f,\theta ,\stackrel{}{r}).$$
(2)
We wish to construct an equation that describes the ensemble average of the spatial variations in the chain densities for a system subject to specified boundary conditions. To do so, we consider the probability that a chain $`\stackrel{}{f}`$ will exist at a point $`\stackrel{}{r}+ϵ\widehat{n}`$, assuming we know $`P(f,\widehat{n},\stackrel{}{r})`$.

A given material will be characterized by a two scalar parameters, $`\lambda `$ and $`Y`$, and two angular functions, $`\varphi _s`$ and $`\varphi _f`$.

* $`\lambda `$ is the average distance in any specified direction between the point where a chain begins and the nearest point that will cause it to split; i.e., the mean length chains would have if there were no fusions.
* $`\varphi _s(\stackrel{}{f}_1|\stackrel{}{f}_2,\stackrel{}{f}_3)`$ is probability that a given splitting of a $`\stackrel{}{f}_1`$-chain results in chains with strengths and directions given by $`\stackrel{}{f}_2`$ and $`\stackrel{}{f}_3`$, normalized to unity in the sense defined below.
* $`\varphi _f(\stackrel{}{f}_1|\stackrel{}{f}_2,\stackrel{}{f}_3)`$ is relative probability that a $`\stackrel{}{f}_2`$-chain and a $`\stackrel{}{f}_3`$-chain will fuse when they intersect and thereby form a $`\stackrel{}{f}_1`$-chain, normalized to unity in the sense defined below.
* $`Y`$ is an overall efficiency with which two intersecting chains will fuse. That is, the total probability that two intersecting chains will fuse is $`Y\varphi _f`$

$`\varphi _s`$ and $`\varphi _f`$ must both include delta functions that enforce force balance at each intersection and Heaviside functions that guarantee all forces are compressive.

We neglect, for now, the finite size of individual grains. That is, we assume that $`\lambda `$ is large compared to a grain diameter and treat the granular packing underlying the force chain network as a continuous medium. The equation governing the spatial variation of $`P`$ as follows (for notational convenience, we have dropped the $`\stackrel{}{r}`$ from the argument of all of the $`P`$’s):

$`(\widehat{n})P(\stackrel{}{f})`$ $`=`$ $`{\displaystyle \frac{1}{\lambda }}{\displaystyle f^{}𝑑f^{}𝑑\theta ^{}f^{\prime \prime }𝑑f^{\prime \prime }𝑑\theta ^{\prime \prime }\left|\mathrm{sin}(\theta ^{}\theta ^{\prime \prime })\right|\varphi _s(\stackrel{}{f}|\stackrel{}{f}^{},\stackrel{}{f}^{\prime \prime })P(\stackrel{}{f})}`$ (3)
$`+{\displaystyle \frac{2}{\lambda }}{\displaystyle f^{}𝑑f^{}𝑑\theta ^{}f^{\prime \prime }𝑑f^{\prime \prime }𝑑\theta ^{\prime \prime }\left|\mathrm{sin}(\theta \theta ^{\prime \prime })\right|\varphi _s(\stackrel{}{f}^{}|\stackrel{}{f},\stackrel{}{f}^{\prime \prime })P(\stackrel{}{f}^{})}`$
$``$ $`2Y{\displaystyle f^{}𝑑f^{}𝑑\theta ^{}f^{\prime \prime }𝑑f^{\prime \prime }𝑑\theta ^{\prime \prime }\left|\mathrm{sin}(\theta \theta ^{\prime \prime })\right|\varphi _f(\stackrel{}{f}^{}|\stackrel{}{f},\stackrel{}{f}^{\prime \prime })P(\stackrel{}{f})P(\stackrel{}{f}^{\prime \prime })}`$
$`+`$ $`Y{\displaystyle f^{}𝑑f^{}𝑑\theta ^{}f^{\prime \prime }𝑑f^{\prime \prime }𝑑\theta ^{\prime \prime }\left|\mathrm{sin}(\theta ^{}\theta ^{\prime \prime })\right|\varphi _f(\stackrel{}{f}|\stackrel{}{f}^{},\stackrel{}{f}^{\prime \prime })P(\stackrel{}{f}^{})P(\stackrel{}{f}^{\prime \prime })}`$
Each term on the right hand side represents a type of event that can alter the density $`P(\stackrel{}{f})`$. We discuss each in turn:

* The first term represents the loss of $`\stackrel{}{f}`$-chains due to splitting. The $`\mathrm{sin}`$ factor is exhibited explicitly for future convenience. Here $`P(\stackrel{}{f})`$ may be taken out of the integral. To ensure that $`\lambda `$ is the mean distance between a chain will propagate before splitting (in the absence of interactions with other chains), we must therefore have

$$f^{}𝑑f^{}𝑑\theta ^{}f^{\prime \prime }𝑑f^{\prime \prime }𝑑\theta ^{\prime \prime }\left|\mathrm{sin}(\theta ^{}\theta ^{\prime \prime })\right|\varphi _s(\stackrel{}{f}|\stackrel{}{f}^{},\stackrel{}{f}^{\prime \prime })=1.$$
(4)
Note that $`\varphi _s`$ has dimensions of $`1/force^4`$.
* The second term represents the creation of $`\stackrel{}{f}`$-chains due to the splitting of chains in other directions. The splitting function here must be the same as in the first term, but with arguments exchanged, since the rate at which chains chains appear due to splitting must match the rate at which the parent chains split. The factor of $`2`$ counts the identical integral arising from the exchange of the prime and double-prime labels.
* The third term, which is nonlinear in $`P`$, represents the loss of $`\stackrel{}{f}`$-chains due to their fusions with $`\stackrel{}{f}^{\prime \prime }`$-chains (or $`\stackrel{}{f}^{}`$-chains, hence the factor of $`2`$). The quantity $`\left|\mathrm{sin}(\theta \theta ^{\prime \prime })\right|P(\stackrel{}{f})P(\stackrel{}{f}^{\prime \prime })`$ is the density of intersections of $`\stackrel{}{f}`$ and $`\stackrel{}{f}^{\prime \prime }`$-chains, and $`\varphi _f`$ is specifies the probability of fusion when two such chains meet. Here again $`P(\stackrel{}{f})`$ can be taken out of the integral. We choose to normalize $`\varphi _f`$ such that the remaining integral has a maximum value of unity over the set of $`P(\stackrel{}{f}^{\prime \prime })`$’s of the form $`P(\stackrel{}{f}^{\prime \prime })=\delta (\stackrel{}{f}^{\prime \prime }\stackrel{}{f}_0)`$, for all values of $`\stackrel{}{f}_0`$. In other words, we normalize $`\varphi _f`$ according to the type of intersection most likely to produce a fusion if all chain densities were identical. In equation form:

$$\mathrm{max}\left[f^{}𝑑f^{}𝑑\theta ^{}f^{\prime \prime }𝑑f^{\prime \prime }𝑑\theta ^{\prime \prime }\varphi _f(\stackrel{}{f}^{}|\stackrel{}{f},\stackrel{}{f}^{\prime \prime })\delta (\stackrel{}{f}^{\prime \prime }\stackrel{}{f}_0)\right]=1,$$
(5)
or, equivalently,

$$\mathrm{max}\left[f^{}𝑑f^{}𝑑\theta ^{}\varphi _f(\stackrel{}{f}^{}|\stackrel{}{f},\stackrel{}{f}_0)\right]=1,$$
(6)
where the maximization is over all possible values of $`\stackrel{}{f}`$ and $`\stackrel{}{f}_0`$. The parameter $`Y`$ then provides an absolute measure of the fusion efficiency. $`\varphi _f`$ has dimensions $`1/force^2`$, in contrast to $`\varphi _s`$.
* The fourth term represents the creation of $`\stackrel{}{f}`$-chains due to the fusion of chains in other directions. For consistency, the fusion function must be the same as in the third term.

For materials composed of perfectly rigid grains, a rescaling of all of the forces in a given DFCN yields another perfectly acceptable network. Hence, Eq. (3) must be invariant under a rescaling of all the force intensities. This is verified by straightforward dimension counting (unlike the form of the Boltzmann equation suggested in Ref. ).

Consider further the general forms of $`\varphi _s`$ and $`\varphi _f`$. For an isotropic material, $`\varphi `$ can depend only on the differences between angles. The general form is

$$\varphi _a(\stackrel{}{f}|\stackrel{}{f}^{},\stackrel{}{f}^{\prime \prime })=F_a(f,f^{},f^{\prime \prime })\delta (\stackrel{}{f}^{}+\stackrel{}{f}^{\prime \prime }\stackrel{}{f})\mathrm{\Theta }(f)\mathrm{\Theta }(f^{})\mathrm{\Theta }(f^{\prime \prime })\psi _a(\theta ^{}\theta ,\theta ^{\prime \prime }\theta ),$$
(7)
where $`F_a`$ is a combination of the force intensities that provides the correct dimensions for $`\varphi _a`$, and $`\psi _a`$ is a symmetric function of its arguments. (One might imagine functions of the ratios of force intensities multiplying the arguments of $`\psi _a`$, but these will always reduce to functions of the angles alone due to the $`\delta `$ function.)

In the case of fusions, dimension counting in Eq. (4) implies $`F_f`$ is dimensionless and hence equal to unity, up to a constant that can be absorbed in $`Y`$. In the case of splitting, on the other hand, dimension counting in Eq. (4) implies $`F_s`$ has dimension $`1/force^2`$. The relations between the $`f`$’s enforced by the $`\delta `$ function guarantee that any combination of $`f`$, $`f^{}`$, and $`f^{\prime \prime }`$ with the right dimensions is equally valid; the differences between them can simply be absorbed into $`\psi _s`$. The natural choice is $`F_s=1/(f^{}f^{\prime \prime })`$. With this choice a constant $`\psi _s=1/\pi ^2`$ corresponds to an equal probability for every possible splitting configuration.

It is clear that $`\lambda `$ can be scaled to unity without loss of generality by choice of the unit of length. Though $`Y`$ is a dimensionless parameter with physical significance, from a mathematical point of view it plays a trivial role in Eq. (3). To solve for $`P`$ for any given value of $`Y`$, we solve the case $`Y=1`$, then simply divide all $`P`$’s by $`Y`$. From here on, we take $`\lambda =1`$ and $`Y=1`$.

### 2.3. Specialization to discrete directions

We now make two simplifying assumptions to arrive at a set of ordinary differential equations that permits analytical solution. First, we assume that $`\varphi _s`$ and $`\varphi _f`$ are nonzero only for vertices of the forms shown in Figure 3. This means that all chains lie in the discrete set of directions $`\theta _n=(n\frac{1}{2})\pi /(2N)`$ (measured from the positive $`z`$ direction, which is downward), and that the strong force at a given vertex is always related to two weak ones by a factor of $`\xi =2\mathrm{cos}(\alpha )`$. Defining $`f_mf_0\xi ^m`$, the continuous function $`P(f,\theta )`$ becomes a set of discrete functions $`P(f_m,\theta _n)`$, which will be denoted $`P_{n,m}\delta (\stackrel{}{f}f_m\widehat{n}_n)`$. In this notation, it is always assumed that $`n`$ is taken modulo $`4N`$. Note that $`P_{n,m}`$ is simply a number per unit length; the dimensions of force are taken care of by the two-dimensional delta function. We will treat the case $`\alpha =\pi /(2N)`$ explicitly here. Generalization to any $`\alpha `$ that is an integer multiple of this is straightforward.

This first assumption corresponds to the following forms for the angular parts of the splitting and fusion functions:

$`\psi _s(\phi ^{},\phi ^{\prime \prime })`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left[\delta (\phi ^{}\alpha )\delta (\phi ^{\prime \prime }+\alpha )+\delta (\phi ^{}+\alpha )\delta (\phi ^{\prime \prime }\alpha )\right].`$ (8)
$`\psi _f(\phi ^{},\phi ^{\prime \prime })`$ $`=`$ $`\{\begin{array}{cc}1\hfill & \mathrm{if}\phi ^{}=\phi ^{\prime \prime }=\alpha \mathrm{or}\phi ^{}=\phi ^{\prime \prime }=\alpha \hfill \\ 0\hfill & \mathrm{otherwise}.\hfill \end{array}`$ (11)
The difference in character between $`\varphi _s`$ and $`\varphi _f`$ stems from the fact that $`\varphi _f`$ appears in integrals that are quadratic in $`P`$. When $`P`$ is a sum of two-dimensional delta functions, the fusion integrals have two more delta function factors than the splitting integrals. The same can be seen in the normalization conditions of Equations (4) and Equations (6).

Second, we assume that the boundary conditions of interest are uniform across the top and bottom of the slab. Translational symmetry in the $`x`$-direction then dictates that all horizontal gradients of $`P_{n,m}`$ vanish.

Inserting these assumptions into Eq. (3) and integrating both sides over a small volume of force space in the vicinity of $`f_m\widehat{n}_n`$, we obtain the following ordinary differential equations for the $`P_{n,m}`$’s:

$`(\mathrm{cos}\theta _n)_zP_{n,m}`$ $`=`$ $`P_{n,m}+P_{n+1,m+1}+P_{n1,m+1}`$ (12)
$`+P_{n+1,m1}P_{n1,m1}P_{n,m}P_{n2,m}P_{n,m}P_{n+2,m}.`$
The result is straightforward, but one must be careful to account for all the trigonometric factors arising from integrating over the delta functions, including the integral on both sides that is required to isolate $`P_{n,m}`$.

The different terms on the right hand side of Eq. (12) correspond to the processes described above that can create or destroy chains. In the present case, the first term represents the decay of $`P_{n,m}`$ along the direction $`\theta _n`$ due to the splitting of an existing $`(n,m)`$ chain. The next two (linear) terms account for the addition to $`P_{n,m}`$ due to the splitting of other chains in adjacent orientations. Note that if a chain is to split and create a contribution to $`P_{n,m}`$, which must have strength $`f_0\xi ^m`$, it must begin with strength $`f_0\xi ^{m+1}`$. The first nonlinear term accounts for events in which two chains of strength $`f_0\xi ^{m1}`$ fuse, and the last two nonlinear terms to events in which chains contributing to $`P_{n,m}`$ are deflected due to fusions with other chains.

In Reference , Eq. (12) was studied in detail for the case of 6-fold symmetry. In that case we have $`\xi =1`$, so the hierarchy of equations indexed by $`m`$ collapses. In addition, $`\mathrm{cos}\theta _2=\mathrm{cos}\theta _5=0`$, so the entire system reduces to four ODEs and two algebraic relations. If we assume mirror symmetry about the vertical axis, these are reduced to two equation that can be solved completely.

In order to observe the variations of $`P`$ with force strength, we must consider symmetries other than 6-fold. To avoid technical complications arising from the presence of strictly horizontal chains having $`\mathrm{cos}\theta =0`$, we consider only the cases of $`4N`$-fold symmetry for positive integers $`N`$. The simplest case is that of 8-fold symmetry, and example of which is illustrated in Figure 1.

Note that the discreteness of the possible orientations of force chains does not imply any discreteness in the possible positions of the splitting or fusion vertices. Note also that the system remains isotropic in the sense that Eq. 12, including all coefficients, is identical for all $`n`$, except for the geometric factor on the left-hand side that results from restricting attention to variations in $`z`$ and not $`x`$.

## 3. Homogeneous networks for symmetric splittings and fusions

General solutions of Eq. (12) are not yet known. We study here the important special case of homogeneous solutions; i.e., the fixed points, for which all gradients vanish. These solutions are both nontrivial and crucial for understanding the response function. Consider the trajectory of the vector $`P_{n,m}`$ as $`z`$ is varied. As detailed in Reference  for the 6-fold case, trajectories that do not pass very close to a stable (or marginally stable) fixed point inevitably lead to divergences at $`z>>\lambda `$ that are reflected in unphysical negative values of some $`P_{n,m}`$. Trajectories that avoid these divergences do so by staying in the vicinity of a fixed point over most of the range of $`z`$. As also noted in Reference , the trivial fixed point $`P_{n,m}=0`$ is unstable and hence cannot serve as a reference state for a linear response theory. Thus the response function we ultimately seek will be dominated by the DFCN structure corresponding to some nontrivial fixed point.

### 3.1. Isotropic solutions

It is instructive to begin with a search for isotropic solutions; i.e., solutions for which $`P_{n,m}=P_m`$ for all $`n`$. The fixed point equation derived from Eq. (12) reads

$$P_m+2P_{m+1}+P_{m1}^22P_m^2=0.$$
(13)
To simplify this recursion relation, define

$$d_m=P_mP_{m1}^2.$$
(14)
Eq. (13) says $`2d_{m+1}d_m=0`$, which implies

$$d_m=d_02^m.$$
(15)
For $`d_00`$, Eq. (14) gives

$$P_{m+1}=P_m^2+d_02^m.$$
(16)
For $`d_0<0`$, this recursion relation either produces a negative numbers or a divergence at large $`m`$, neither of which is physically acceptable. If $`P_0`$ is sufficiently small, the repeated squaring drives $`P_m`$ toward zero at least as fast as $`p^{2^m}`$ with $`p>1`$. The negative contribution from the $`d_0`$ term can only speed up this decay as long as $`P_m`$ is positive. But this means that $`P_m`$ decays faster than $`2^m`$, so the $`d_0`$ term eventually makes $`P_m`$ negative. Now if $`P_0`$ is large enough to avoid $`P_m`$ decaying toward $`0`$, repetitive squaring causes it to diverge like $`p^{2^m}`$ with $`p>1`$; the $`d_0`$ term becomes irrelevant at large $`m`$.

For $`d_0>0`$, on the other hand, problems arise for large, negative $`m`$. We write $`P_{m1}=\sqrt{P_md_02^m}`$ and see that complex values will be generated unless $`P_m`$ is always greater than $`d_02^m`$. But this is impossible because $`P_{m1}`$ is strictly less than $`P_m`$.

Thus we are left with $`d_0=0`$ as the only physically relevant case. When $`d_0=0`$, Eq. (12) is satisfied in a special way: both $`P_m+P_{m1}^2=0`$ and $`2P_{m+1}2P_m^2=0`$ are satisfied simultaneously. These two relations are actually identical up to a shift in $`m`$, and admit a one-parameter family of solutions

$$P_m=p^{2^m},$$
(17)
where we must have $`p>1`$ to avoid divergence at large $`m`$. Despite its apparent simplicity, it is instructive to examine this solution in a bit more detail. We consider first the behavior at large negative $`m`$, then the behavior at large positive $`m`$.

Note that $`P_m1`$ for large negative $`m`$. This means that the densities of all chains of intensity $`f_0\xi ^m`$ for large $`m`$ are the same. In other words, the total density of force chains, $`_mP_m`$ diverges. Though this may appear troublesome, it does not yield a divergence in the pressure: $`_mP_mf_m`$ remains perfectly finite. The divergence in chain density is an artifact of the restriction to a fusion function $`\varphi _f`$ that does not allow weak chains to fuse with stronger ones. Removal of this artifact requires working with a continuous distribution of force chain directions and intensities and it beyond the scope of this work.

For large positive $`m`$, we can compute the probability distribution for contact forces, a quantity that has been shown in a variety of experiments and numerical simulations to decay exponentially for large forces. For the homogeneous, isotropic solution, we have

$$𝒫(f_m)=P_m/(f_mf_{m1}),$$
(18)
where $`f_m=f_0\xi ^m`$ and $`P_m=p^{2^m}`$. The denominator in this expression is just the spacing between points with successive $`m`$’s, the required conversion factor from the discrete density $`P_m`$ to the density per unit force. This can be rewritten as

$$𝒫(f)p^{(f/f_0)^\beta }/f,\mathrm{with}\beta =\frac{\mathrm{log}2}{\mathrm{log}\xi }.$$
(19)
Recall that $`p>1`$ and $`\xi <2`$, so $`\beta >1`$ and $`𝒫`$ decays faster than exponentially for large $`f`$. For the 8-fold case, we have $`\beta =2`$ and hence a Gaussian decay. But smaller splitting angles give $`\xi `$ closer to 2, and for a splitting angle as large as $`2\alpha =60^{}`$, we have $`\beta 1.26`$. In this case, the full distribution (including the $`1/f`$ term) can look surprisingly like a simple exponential over several decades, as shown in Figure 4.

Experience with the 6-fold model suggests that the isotropic fixed points of Eq. (12) are non-generic. In the 6-fold model, linearization of the theory in the vicinity of a fixed point generally leads to hyperbolic response functions with peaks that are sharper for more strongly anisotropic fixed points. The isotropic case is pathological in that the peak width is divergent and additional cancellations conspire to produce single peak, but one that obeys different scaling laws from the predictions of elliptic models. It is therefore important to consider solutions of Eq. (12) that are not constrained to be isotropic.

### 3.2. Anisotropic solutions

To find fixed points of Eq. (12), we begin with the following ansatz :

$$P_{n,m}=p^{q_m+x_n\zeta ^m},$$
(20)
where $`p`$ is a real parameter and $`q_m`$, $`x_n`$, and $`\zeta `$ are to be determined. (This guess combines features of the isotropic solution above and the 6-fold solutions discussed in Reference .) Substituting into Eq. (12), we have

$`0`$ $`=`$ $`p^{q_m+x_n\zeta ^m}+p^{q_{m+1}+x_{n1}\zeta ^{m+1}}+p^{q_{m+1}+x_{n+1}\zeta ^{m+1}}`$ (21)
$`+p^{2q_{m1}+(x_{n+1}+x_{n1})\zeta ^{m1}}p^{2q_m+(x_n+x_{n2})\zeta ^m}p^{2q_m+(x_n+x_{n+2})\zeta ^m}.`$
Taking a cue from the isotropic case, we consider the possibility that the cancellation required by this equation occurs term by term in the following way:333Unlike in the isotropic case, we do not prove here that all physically plausible solutions must have this form.

$`p^{q_m+x_n\zeta ^m}`$ $`=`$ $`p^{2q_{m1}+(x_{n+1}+x_{n1})\zeta ^{m1}};`$ (22)
$`p^{q_{m+1}+x_{n1}\zeta ^{m+1}}`$ $`=`$ $`p^{2q_m+(x_n+x_{n2})\zeta ^m};`$ (23)
$`p^{q_{m+1}+x_{n+1}\zeta ^{m+1}}`$ $`=`$ $`p^{2q_m+(x_n+x_{n+2})\zeta ^m}.`$ (24)
Note that all three of these equations are identical up to shifts of $`n`$ and $`m`$. Since the $`q`$ terms are independent of $`n`$, we obtain two recursion relations:

$`q_{m+1}`$ $`=`$ $`2q_m;`$ (25)
$`x_n\zeta `$ $`=`$ $`x_{n1}+x_{n+1}.`$ (26)
Eq. (25) implies

$$q_m=q_02^m,$$
(28)
where the negative sign defines a convention for the sign of $`q_0`$. Eq. (26) is an eigenvalue equation for a $`4N\times 4N`$ matrix $`\stackrel{}{M}`$. The eigenvalues $`\zeta _j`$ and eigenvectors $`y_n^{(j)}`$ of $`\stackrel{}{M}`$ are easily found:

$$\zeta _j=2\mathrm{cos}\left(\frac{j\pi }{2N}\right);y_n^{(j)}=\{\begin{array}{cc}\mathrm{cos}\left(\frac{jn\pi }{2N}\right)\hfill & j=0,1,\mathrm{}2N\hfill \\ \mathrm{sin}\left(\frac{jn\pi }{2N}\right)\hfill & j=2N+1,\mathrm{}4N1\hfill \end{array}$$
(29)
Thus for each of the $`4N`$ eigenvalues and eigenvectors, we have what appears to be a three-parameter family of fixed point solutions of Eq. (12) of the form

$$P_{n,m}=p^{q_02^m+x_0y_n^{(j)}(\zeta _j)^m},$$
(30)
the free parameters being $`p`$, $`q_0`$, and $`x_0`$. Recall that the force intensity associated with the chains whose density is $`P_{n,m}`$ is $`f_m=f_0\xi ^m`$ and note that $`\xi =\zeta _1`$. Note also that for $`j=0`$, we have $`\zeta =2`$, so that the $`x_0`$ and $`q_0`$ terms can be combined. This leaves us with precisely the same isotropic solution discussed above.

Eq. (30) requires further examination to determine the range of different structures it represents. First, each three-parameter family is really only a one-parameter family of distinct physical networks. There are ways of rescaling the parameters that have no effect on the physical system begin described:

$`(q_0,x_0)`$ $``$ $`(2^\nu q_0,\zeta ^\nu x_0)\mathrm{and}`$ (31)
$`(p,q_0,x_0)`$ $``$ $`(p^{q_0},1,x_0/q_0).`$ (32)
The first of these corresponds to a redefinition of the force scale, $`f_0f_0\xi ^\nu `$, which cannot affect the physics. Strictly speaking, the system is exactly invariant under this transformation only if $`\nu `$ is an integer. Arbitrary values of $`\nu `$, however, merely shift the positions of the discrete set of $`P_{n,m}`$’s along a single smooth curve for each $`n`$, with integral $`\nu `$’s being the values that shift the entire set into itself. For our purposes, the differences between solutions related by non-integer $`\nu `$ are unimportant details. (See Figure 5.) The second transformation leads to a mathematically identical solution. (The case $`q_0=0`$ is not physically relevant, as it leads to divergence of the total stress from large $`m`$ force chains.) By performing the first transformation with $`\nu =\mathrm{log}(x_0/q_0)/\mathrm{log}(2/\zeta )`$, followed by the second, we can adjust $`q_0`$ and $`x_0`$ both to unity.

Second, note that we cannot take arbitrary linear combinations of the solutions to Eq. (26), since they would not satisfy the nonlinear Eq. (12). We can, however, take linear combinations of eigenvectors that have degenerate eigenvalues. Since $`\zeta _{4Nj}=\zeta _j`$ for $`1j2N1`$ we can generate solutions from arbitrary linear combinations of the degenerate pair of eigenvectors $`\stackrel{}{y}^{(j)}`$ and $`\stackrel{}{y}^{(4Nj)}`$. Thus $`j=0`$ and $`j=2N`$ each provide a single continuous family of solutions parameterized by $`p`$, while all other $`j`$’s come in pairs that provide families parameterized by $`p`$ and by the relative weights of the two $`\stackrel{}{y}^{(j)}`$’s.

Consider now the family spanned by $`\stackrel{}{y}^{(j)}`$ and $`\stackrel{}{y}^{(4Nj)}`$, for some specific $`j`$. By forming linear combinations, we can construct basis vectors $`\stackrel{}{e}_j`$ that are symmetric or antisymmetric under reflection through a vertical axis. The symmetric combination will give solutions for which $`P_{n,m}=P_{4N+1n,m}`$, as might be expected if the average confining force on the plates is purely normal. The antisymmetric combination will give solutions that break this symmetry, though they do not yield any negative values of $`P_{n,m}`$. These might be important for describing systems subjected to shear as well as compression. The explicit forms of the symmetric combinations (including the trivially symmetric isotropic case) can be expressed as

$$e_n^{(j)}=\{\begin{array}{cc}1\hfill & j=0\hfill \\ \sqrt{2}\mathrm{cos}(j\theta _n)\hfill & j=1,\mathrm{},2N1,\hfill \end{array}$$
(33)
where the normalization is $`\stackrel{}{e}_j=4N`$ for all $`j`$.

The antisymmetric combinations, normalized to $`4N`$, are given by

$$e_n^{(j)}=\{\begin{array}{cc}(1)^n\hfill & j=2N\hfill \\ \sqrt{2}\mathrm{sin}(j\theta _n)\hfill & j=2N+1,\mathrm{},4N1,\hfill \end{array}$$
(34)
In the following, we focus on the symmetric cases only.

The $`\stackrel{}{e}^{(0)}`$ solution is the isotropic case discussed above. The $`\stackrel{}{e}^{(1)}`$ solution describes a stress configuration with a dipole-like anisotropy. At each generation $`m`$, the angular variation $`P_{n,m}`$ has a single maximum (for $`p>1`$) at $`n=1`$ and a minimum $`n=2N`$. There are more downward-pointing than upward-pointing chains for each undirected orientation. Moreover, because $`\zeta _1`$ is positive, the decay of $`P_{n,m}`$ with $`m`$ is monotonic both at large $`m`$, where it decays to $`0`$, and at large negative $`m`$, where it decays to $`1`$. Figure 5 shows the form of $`P_{n,m}`$ for this case. Note that the term “dipole-like” applies here to the $`P`$’s. Since the stress field does not distinguish between force chains in opposite directions, the anisotropy in the stress is “quadrupole-like” – stronger in the vertical directions than in the horizontal.

A striking feature of the solutions is the occurrence of peaks for some $`n`$, which become more pronounced for larger $`p`$. Peaks are evident for all $`n`$ when the stress, $`f_mP_{n,m}`$, is plotted rather than the chain density. (See Figure 6.) Differentiating Eq. (30) with $`q_0=x_0=1`$, we find that the peak in $`P_{n,m}`$ occurs at

$$m=\frac{\mathrm{log}y_n^{(j)}+\mathrm{log}\mathrm{log}|\zeta _n|\mathrm{log}\mathrm{log}2}{\mathrm{log}2\mathrm{log}|\zeta _n|},$$
(35)
which is independent of $`p`$. (The peak in $`f_mP_{n,m}`$ does depend weakly on $`p`$.) The force scale represented by this peak must have its origin in the boundary conditions, since a linear rescaling of all forces in the bulk has no effect on the DFCN structure. Recall that part of the process required to collapse the solutions onto the one-parameter family indexed by $`p`$ was to rescale $`f_0`$, which corresponds to a rescaling of the strengths of the force chains injected at the boundaries. There are subtleties lurking here, however, as will be discussed further in Section 3.3.

A full calculation of the single contact force distribution $`𝒫(f)`$ is slightly more complicated than for the isotropic case, as it requires a sum over the different chain orientations. It is not hard to extract the dominant term at large $`f`$, however. For the symmetric $`j=1`$ case with $`q_0=x_0=1`$, Eq. (30) can be written as

$$P_{n,m}=p^{f_m^\beta +e_n^{(1)}f_m^\gamma },$$
(36)
where $`\beta =\mathrm{log}2/\mathrm{log}\xi `$ as above, and $`\gamma =\mathrm{log}\zeta _1/\mathrm{log}\xi =1`$. Since $`\beta >1`$, the dominant term at large $`f`$ will be $`P_{1,m}=p^{f_m^\beta +f}`$.

It is important to note that the anisotropy of the solution is not a consequence of any material anisotropy. The bulk properties of the material we are describing are perfectly isotropic. The boundary conditions on the system will generally not be isotropic, however, so there is no reason to expect that the pre-stressed state about which the response function should ultimately be calculated should have isotropic stress. The generic situation described by the discrete DFCN theory presented here is one in which an intrinsically isotropic material is subjected to loading forces that generate anisotropic distributions of stress chains. Because different degrees of anisotropy correspond to different fixed points of the full nonlinear theory and hence to different linearizations, the pre-stressed state may have a profound influence on the response function. This is in contrast to an ordinary elastic medium, in which the response function is not affected by an anisotropic background stress.

The solutions for higher $`j`$ have features that seem a bit unlikely from a physical standpoint, though not obviously unacceptable. For $`2j<N`$, we have higher orders of anisotropy, with $`P_{n,m}`$ exhibiting additional maxima and minima as a function of $`n`$ at each $`m`$, but still with monotonic decay at large $`m`$. For $`N<j2N`$, the eigenvalues are negative, which causes $`\zeta ^m`$ and hence $`P_{n,m}`$ to oscillate as $`m`$ varies. Whether or not these solutions are relevant hinges on the question of which boundary conditions are associated with generic physical systems.

### 3.3. Remark on boundary conditions

Given several families of homogeneous solutions, it is natural to ask which solution will be selected for a given set of boundary conditions. From a mathematical perspective, this is a straightforward question. As described in detail for the 6-fold case in Reference , the only way to specify boundary conditions on Eq. (12) that are sure to be self-consistent and not produce negative chain densities is to specify the densities of only the ingoing chains at each boundary. In general, the densities specified will not match perfectly any of the homogeneous solutions $`P_{n,m}`$ for all $`m`$ and $`n`$, even if horizontal translational symmetry is assumed. Based on experience with the 6-fold case, we might expect the trajectory in the space of $`P_{n,m}`$ as the slab is traversed to quickly approach a fixed point and stay near it over a substantial portion of the slab, then make a rapid transition to the vicinity of another fixed point, and finally move off to a point consistent with the bottom boundary condition. The details of the transition regions and selection of the fixed point remain to be worked out in the present case.

From a physical perspective, the situation is not so clear. Consider the following two possible choices of boundary conditions for the top surface. In both cases, assume that $`P_{n,m}=0`$ for all ingoing $`(n,m)`$ except the ones specified as follows: (BC1) $`P_{0,5}=P_{1,5}=1`$; and (BC2) $`P_{0,1}=P_{1,1}=\xi ^5`$. Assume also that the boundary conditions on the bottom slab are simple reflections of these through the horizontal. Thus we have boundary conditions with the same total stress stipulated, but in one case it is injected via a sparse set of strong forces, and in the other via a denser set of weaker forces. These conditions will not necessarily lead to the same fixed points in the bulk. Since the chain densities are measured with respect to the intrinsic length scale $`\lambda `$, BC1 and BC2 are not simply related by a scale transformation. We therefore expect them to lead to fixed points corresponding to different $`p`$’s, which implies that they will require different rescalings of $`f_0`$. This in turn indicates that the force scales fixed by the two different boundary conditions differ, in spite of the fact that they correspond to specifying the same overall force.

The resolution of this apparent paradox lies in the fact that the boundary conditions do not specify the densities of the outgoing chains. Ultimately, the force scale must be determined by the $`zz`$ component of the stress, which must be constant throughout the slab; i.e., $`\sigma _{zz}`$ must determine the prescribed rescaling of $`f_0`$. But the densities of outgoing chains produced by BC1 and BC2 will differ, which will cause $`\sigma _{zz}`$ to differ in the two cases. In order to determine $`\sigma _{zz}`$, which is the total pressure applied to the top surface, we must first determine the full solution, or at least which fixed point the trajectory approaches.

The converse situation, in which the stress is specified but not the individual chain densities, is more familiar, but just as subtle. One might expect it to be sufficient to specify some components of the stress tensor and its spatial derivatives on the top boundary and others on the bottom, but the DFCN equations require independent specifications of densities of chains of many directions and strengths. In an experiment in which a granular material is put under compression between to flat plates, it is not at all clear how the chains densities at the boundaries are determined. Both the distribution over directions $`n`$ and strengths $`m`$ are relevant. An important topic of future research will be to develop a method of assigning directions to force chains observed in experimental images or numerical models, and to learn from them how the boundary conditions appear to be imposed.

## 4. Conclusion

In this paper we have exhibited for the first time solutions of the Boltzmann equation for the chain densities in a model that permits the interaction of chains of different intensities. These solutions have highly nontrivial structure that raises a number of interesting questions. Most importantly, they permit calculations of both stresses and single particle force distributions in a unified way. Refinement of the model is clearly necessary: we would like to avoid divergences in the density of weak chains; and we would like to generalize the splitting and fusion functions, which might bring the single contact force distribution $`𝒫(f)`$ into closer agreement with experiment. Nevertheless, there is much to be gained from further study of the models presented here, or minor modifications of them.

The availability of closed form analytic solutions suggests that a complete theory of the stress configurations and response functions for these particular models is within reach. In addition, it appears possible that the linearized theory in the vicinity of the fixed points presented here could lead to a new derivation of constitutive relations for closing the equilibrium equations for the macroscopic stresses. Finally, the solutions presented here may provide a basis for extensions to similar discrete models in which more than one type of vertex structure is allowed. All of these topics are currently under study. The author looks forward to continued collaboration with David Schaeffer on these topics, as he is certain to offer valuable insights and pose questions that lead to fruitful calculations.

## Acknowledgments

I thank J.-P. Bouchaud, P. Claudin, M. Otto, and D. Schaeffer for extensive discussions and the group of E. Clément and the LMDH at Jussieu in Paris for its hospitality. I am especially grateful to P. Claudin for his suggestions for improving a draft of this paper. This work was supported by the NSF Grant PHY-98-70028.

Received December 2002.


format: text

text:
# Evolution of Curves and Surfaces by Mean Curvature*

## 1. Introduction

Traditionally, differential geometry has been the study of curved spaces or shapes in which, for the most part, time did not play a role. In the last few decades, on the other hand, geometers have made great strides in understanding shapes that evolve in time. There are many processes by which a curve or surface can evolve, but among them one is arguably the most natural: the mean curvature flow. This article describes the flow, some of the discoveries that have been made about it, and some unresolved questions.

## 2. Curve-shortening flow

The simplest case is that of curves in the plane. Here the flow is usually called the “curvature flow” or the “curve-shortening flow”. Consider a smooth simple closed curve in the plane, and let each point move with a velocity equal to the curvature vector at that point. What happens to the curve?

The evolution has several basic properties. First, it makes the curve smoother. Consider a portion of a bumpy curve as in figure 1(a). The portions that stick up move down and the portions that stick down move up, so the curve becomes smoother or less bumpy as in figure 1(b). The partial differential equation for the motion is a parabolic or heat-type equation, and such smoothing is a general feature of solutions to such equations. Thus, for example, even if the initial curve is only $`C^2`$, as it starts moving it immediately becomes $`C^{\mathrm{}}`$ and indeed real analytic.

1(a) leave a big chunk of space here 1(b)

Figure 1: Smoothing

However, there is an important caveat: the smoothing may only be for a short time. If the curve is $`C^2`$ at time $`0`$, it will be real-analytic for times $`t`$ in some interval $`(0,ϵ)`$. But because the equation of motion is nonlinear, the general theory of parabolic equations does not preclude later singularities. And indeed, as we shall see, any curve must eventually become singular under the curvature flow.

The simplest closed curve is of course a circle. The flow clearly preserves the symmetry, so in this case it is easy to solve the equation of motion explicitly. One finds that a circle of radius of radius $`r`$ at time $`0`$ shrinks to a circle of radius $`\sqrt{r^22t}`$ at time $`t`$, so that time $`t=r^2/2`$ the circle has collapsed to a point and thereby become singular.

The second fundamental property is that arclength decreases. The proof is as follows. For any evolution of curves,

$$\frac{d}{dt}\text{(length)}=𝕜𝕧𝑑s,$$

where $`𝕜`$ is the curvature vector, $`𝕧`$ is the velocity, and $`ds`$ is arclength. For the curvature flow, $`𝕧=𝕜`$, so the right hand side of the equation is clearly negative. Indeed, the proof shows that this flow is, in a sense, the gradient flow for the arclength functional. Thus, roughly speaking, the curve evolves so as to reduce its arclength as rapidly as possible. This explains the name “curve-shortening flow”, though many other flows also reduce arc-length.

The third property is that the flow is collision-free: two initially disjoint curves must remain disjoint. The idea of the proof is as follows. Suppose that two initially disjoint curves, one inside the other, eventually collide. At the first time $`T`$ of contact, they must touch tangentially. At the point of tangency, the curvature of the inner curve is greater than or equal to the curvature of the outer curve. Suppose for simplicity that strict inequality holds. Then (at the point of tangency) the inner curve is moving inward faster than the outer curve is. But then at a slighly earlier time $`Tϵ`$, the curves would have to cross each other. But that contradicts the choice of $`T`$ (as first time of contact), proving that contact can never occur.

This collision avoidance is a special case of the maximum principle for parabolic differential equations. The maximum principle also implies in the same way that an initially embedded curve must remain embedded.

The fourth fundamental property is that every curve $`\mathrm{\Gamma }`$ has a finite lifespan. To undertand why, consider a circle $`C`$ that contains $`\mathrm{\Gamma }`$ in its interior. Let $`\mathrm{\Gamma }`$ and $`C`$ both evolve. The circle collapses in a finite time. Since the curves remain disjoint, $`\mathrm{\Gamma }`$ must disappear before the circle collapses.

There is another nice way to see that a curve must become singular in a finite time:

Theorem 1. If $`A(t)`$ is the area enclosed by the curve at time $`t`$, then $`A^{}(t)=2\pi `$ until the curve becomes singular. Thus a singularity must develop within time $`A(0)/2\pi `$.

Proof. For any evolution,

$$A^{}(t)=_{\mathrm{\Gamma }(t)}𝕧𝕟𝑑s$$

where $`\mathrm{\Gamma }(t)`$ is the curve at time $`t`$, $`𝕧`$ is the velocity, $`𝕟`$ is the outward unit normal, and $`ds`$ is arclength. For the curvature flow, $`𝕧=𝕜`$, so

$$A^{}(t)=_\mathrm{\Gamma }𝕜𝕟𝑑s$$

which equals $`2\pi `$ by the Gauss-Bonnet theorem.∎

The first deep theorem about curvature flow was proved by Mike Gage and Richard Hamilton in 1986 \[GH\]:

Theorem 2. Under the curvature flow, a convex curve remains convex and shrinks to a point. Furthermore, it becomes asymptotically circular: if the evolving curve is dilated to keep the enclosed area constant, then the rescaled curve converges to a circle.

This theorem is often summarized by stating that convex curves shrink to round points.

The proof is too involved to describe here, but I will point out that the result is not at all obvious. Consider for example a long thin ellipse, with the major axis horizontal. The curvature is greater at the ends than at the top and bottom, so intuitively it should become rounder. But the ends are much farther from the center than the top and bottom are, so it is not clear that they all reach the center at the same time. Thus it is not obvious that the curve collapses to a point rather than a segment.

Indeed, there are natural flows that have all the above-mentioned basic properties of curvature flow but for which the Gage-Hamilton theorem fails. Consider for example a curve moving in the direction of the curvature vector but with speed equal to the cube root of the curvature. Under this flow, any ellipse remains an ellipse of the same eccentricity and thus does not become circular. For this flow, Ben Andrews \[A1\] has proved that any convex curve shrinks to an elliptical point. (See also \[AST,  SaT\].) For other flows (e.g. if “cube root” is replaced by $`r`$th root for any $`r>3`$), a convex curve must shrink to a point but in a very degenerate way: if the evolving curve is dilated to keep the enclosed area constant, then length of the rescaled curve tends to infinity \[A3\]. More generally, Andrews has studied the existence and nonexistence of asymptotic shapes for convex curves for rather general classes of flows \[A2,  A3\].

Shortly after Gage and Hamilton proved their theorem, Matt Grayson proved what is still perhaps the most beautiful theorem in the subject:

Theorem 3. \[G1\] Under the curvature flow, embedded curves become convex and thus (by the Gage-Hamilton theorem) eventually shrink to round points.

Again, the proof is too complicated to describe here, but let me indicate why the result is very surprising. Consider the annular region between a two concentric circles of radii and $`r=1`$ and $`R=2`$. Form a curve in this annular region by spiraling inward $`n`$ times, and then back out $`n`$ times to make a closed embedded curve. Figure 2 shows such a curve with $`n=3/2`$, but think of $`n`$ being very large, say $`10^{100}`$. Recall that the curve exists for a time at most $`R^2/2=2`$. By Grayson’s theorem, the curve manages, amazingly, to unwind itself and become convex in this limited time. Incidentally, notice that initially, except for two very small portions, the curve is not even moving fast: its curvature is no more than that of the inner circle.

Figure 2: Spiral

As a corollary to Grayson’s theorem, one gets an exact formula for the lifespan of any curve. Recall that the area enclosed by a curve decreases with constant speed $`2\pi `$ as long as the curve is smooth. By Grayson’s theorem, the curve remains smooth until its area becomes $`0`$. Thus the lifespan of any embedded curve must be exactly equal to the initial area divided by $`2\pi `$.

Grayson later generalized his theorem by proving that a closed curve moving on a compact surface by curvature flow must either collapse to a round point in a finite time or else converge to a simple closed geodesic as $`t\mathrm{}`$ \[G2\].

## 3. Mean curvature flow for surfaces

We now leave curves and consider instead moving surfaces. Recall that at each point of an $`n`$-dimensional hypersurface in $`^{n+1}`$, there are $`n`$-principal curvatures $`\kappa _1,\kappa _2,\mathrm{},\kappa _n`$ (given a choice $`𝕟`$ of unit normal). Their sum $`h`$ is the scalar mean curvature, and the product $`H=h𝕟`$ of the scalar mean curvature and the unit normal is the mean curvature vector. The mean curvature vector does not depend on choice of normal since replacing $`𝕟`$ by $`𝕟`$ also changes the sign of the scalar mean curvature. In the mean curvature flow, a hypersurface evolves so that its velocity at each point is equal to the mean curvature vector at that point.

The basic properties of curvature flow also hold for mean curvature flow:

The analog of the Gage-Hamilon theorem also holds, as Gerhard Huisken \[H1\] proved:

Theorem 4. For $`n2`$, an $`n`$-dimensional compact convex surface in $`^{n+1}`$ must shrink to a round point.

Oddly enough, Huisken’s proof does not apply to the case of curves ($`n=1`$) considered by Gage and Hamilton. Huisken’s proof shows that the asymptotic shape is totally umbilic: at each point $`x`$, the principal curvatures are all equal (though a priori they may vary from point to point). For $`n2`$, the only totally umbilic surfaces are spheres, but for $`n=1`$, the condition is vacuous.

The analog of Grayson’s theorem, however, is false for surfaces. Consider for example two spheres joined by a long thin tube. The spheres and the tube both shrink, but the mean curvature along the tube is much higher than on the spheres, so the middle of the tube collapses down to a point, forming a singularity. The surface then separates into two components, which eventually become convex and collapse to round points.

Thus, unlike a curve, a surface can develop singularities before it shrinks away. This raises various questions:

In the rest of this article I will describe some partial answers to these questions.

A great many results about mean curvature flow have been proved using only techniques of classical differential geometry and partial differential equations. However, most of the proofs are valid only until the time that singularities first occur. Once singularities form, the equation for the flow does not even make sense classically, so analyzing the flow after a singularity seems to require other techniques.

Fortunately, using non-classical techniques, namely the geometric measure theory of varifolds and/or the theory of viscosity or level-set solutions, one can define notions of weak solutions for mean curvature flow and one can prove existence of a solution (with a given initial surface) up until the time that the surface disappears.

The definitions are somewhat involved and will not be given here. The different definitions are equivalent to each other (and to the classical definition) until singularities form, but are not completely equivalent in general. For the purposes of this article, the reader may simply accept that there is a good way to define mean curvature flow of possibly singular surfaces and to prove existence theorems. The notion of mean curvature flow most appropriate to this article is Ilmanen’s “enhanced Brakke flow of varifolds” \[I\].

## 4. Non-uniqueness or fattening

If a surface is initially smooth, classical partial differential equations imply that there is a unique solution of the evolution equation until singularities form. However, once a singularity forms, the classical uniqueness theorems do not apply. In the early 1990’s various researchers, including De Giorgi, Evans and Spruck, and Chen, Giga, and Goto, asked whether uniqueness could in fact break down after singularity formation. They already knew that uniqueness did fail for certain initially singular surfaces, but they did not know whether an initially regular surface could later develop singularities that would result in non-uniqueness.

A technical aside: the above-named people did not phrase the question in terms of uniqueness but rather in terms of “fattening”. They were all using a level set or viscosity formulation of mean curvature flow, in which solutions are almost by definition unique. But non-uniqueness of the enhanced varifold solutions corresponds to fattening of the viscosity solution in the following sense. If a single initial surface $`M`$ gives rise to different enhanced varifold solutions $`M_t^1,M_t^2,\mathrm{},M_t^k`$, then the viscosity solution “surface” at time $`t`$ will consist of the various $`M_t^i`$’s together with all the points in between. Thus if $`k>1`$, the viscosity surface will in fact have an interior. Since the surface was initially infinitely thin, in developing an interior it has thereby “fattened”.

Recently Tom Ilmanen and I settled this question \[IW\]:

Theorem 5. There is a compact smooth embedded surface in $`^3`$ for which uniqueness of (enhanced varifold solutions of) mean curvature evolution fails. Equivalently, the viscosity (or level set) solution fattens.

The idea of the proof is as follows. Consider a solid torus of revolution about the z-axis centered at the origin, a ball centered at the origin that is disjoint from the torus, and $`n`$ radial segments in the $`xy`$-plane joining the ball to the torus. Call their union $`W`$. Now consider a nested one-parameter family of smooth surfaces $`M^ϵ`$ ($`0<ϵ<1`$) as follows. When $`ϵ`$ is small, the surface should be a smoothed version of the set of points at distance $`ϵ`$ from $`W`$. This $`M^ϵ`$ looks like a wheel with $`n`$ spokes. The portion of the $`xy`$-plane that is not contained in $`M^ϵ`$ has $`n`$ simply-connected components, which we regard as holes between the spokes of the wheel. As $`ϵ`$ increases, the spokes get thicker and the holes between the spokes get smaller. When $`ϵ`$ is close to $`1`$, the holes should be very small, and near each hole the surface should resemble a thin vertical tube.

Now let $`M^ϵ`$ flow by mean curvature. If $`ϵ`$ is small, the spokes are very thin and will quickly pinch off, separating the surface into a sphere and a torus. If $`ϵ`$ is large, the holes between the spokes are very small and will quickly pinch off, so the surface becomes (topologically) a sphere. By a continuity argument, there is an intermediate $`ϵ`$ such that both pinches occur simultaneously.

For this particular $`ϵ`$, we claim that the simultaneous pinching immediately results in non-uniqueness, at least if $`n`$ is sufficiently large. Indeed, at the moment of simultaneous pinching, the surface will resemble a sideways figure 8 curve revolved around the $`z`$-axis as indicated in figure 3(a). Of course the surface will not be fully rotationally symmetric about the axis, but it will have $`n`$-fold rotational symmetry, and here I will be imprecise and proceed as though it were rotationally symmetric.

3(a)

3(b) leave a largish chunk of space here 3(c)

Figure 3: Non-uniqueness

There is necessarily one evolution in which the surface then becomes topologically a sphere as in figure 3(b). If the angle $`\theta `$ in figure 6 is sufficiently small, there is also another evolution, in which the surface detaches itself from the $`z`$-axis and thereby becomes a torus as in figure 3(c).

One can show that as $`n\mathrm{}`$, the angle $`\theta `$ tends to $`0`$. Thus if $`n`$ is large enough, the angle will be very small and both evolutions will occur.

The proof unfortunately does not give any bound on how large an $`n`$ is required. Numerical evidence \[AIT\] suggests that $`n=4`$ suffices. The case $`n=2`$ seems to be borderline and the case $`n=3`$ has not been investigated numerically.

However, the argument completely breaks down for $`n<2`$. Indeed, I would conjecture that if the initial surface is a smooth embedded sphere or torus, then uniqueness must hold.

It is desirable to know natural conditions on the initial surface that guarantee uniqueness. As just mentioned, genus $`1`$ may be such a condition. Mean convexity (described below) and star-shapedness are known to guarantee uniqueness. The latter is interesting because the surface will typically cease to be star-shaped after a finite time, but its initially starry shape continues to ensure uniqueness \[So\].

Fortunately, uniqueness is known to hold generically in a rather strong sense. If a family of hypersurfaces foliate an open set in $`^{n+1}`$, then uniqueness will hold for all except countably many of the leaves. Of course any smooth embedded surface is a leaf of such a foliation, so by perturbing the surface slightly, we get a surface for which uniqueness holds.

## 5. The size of singular sets

For general initial surfaces, our knowledge about singular sets is very limited. Concerning the size of the singular set, Tom Ilmanen \[I\], building on earlier deep work of Ken Brakke \[BK\], proved the following theorem. (See also \[ES IV\].)

Theorem 6. For almost every initial hypersurface $`M_0`$ and for almost every time $`t`$, the surface $`M_t`$ is smooth almost everywhere.

This theorem reminds me of Kurt Friedrichs, who used to say that he did not like measure theory because when you do measure theory, you have to say “almost everywhere” almost everywhere.

Aside from objections Friedrichs might have had, the theorem is unsatisfactory in that a much stronger statement should be true. But it is a tremendous acheivement and it is the best result to date for general initial surfaces.

However, for some classes of initial surfaces, we now have a much better understanding of singularities. In particular, this is the case when the initial surface is mean-convex. The rest of this article is about such surfaces. For simplicity of language, only two-dimensional surfaces in $`^3`$ will be discussed, but the results all have analogs for $`n`$-dimensional surfaces in $`^{n+1}`$ or, more generally, in $`(n+1)`$-dimensional riemannian manifolds.

Consider a compact surface $`M`$ embedded in $`^3`$ and bounding a region $`\mathrm{\Omega }`$. The surface is said to be “mean-convex” if the mean curvature vector at each point is a nonnegative multiple of the inward unit normal (that is, the normal that points into $`\mathrm{\Omega }`$.) This is equivalent to saying that under the mean curvature flow, $`M`$ immediately moves into $`\mathrm{\Omega }`$. Mean convexity is a very natural condition for mean curvature flow:

Mean convexity, although a strong condition, does not preclude interesting singularity formation. For example, one can connect two spheres by a thin tube as described earlier in such a way that the resulting surface is mean convex. Thus neck pinch singularities do occur for some mean convex surfaces.

Theorem 7. \[W1\] A mean convex surface evolving by mean curvature flow in $`^3`$ must be completely smooth (i.e., with no singularities) at almost all times, and at no time can the singular set be more than $`1`$ dimensional.

This theorem is in some ways optimal. For example, consider a torus of revolution bounding a region $`\mathrm{\Omega }`$. If the torus is thin enough, it will be mean convex. Because the symmetry is preserved and because the surface always remains in $`\mathrm{\Omega }`$, it can only collapse to a circle. Thus at the time of collapse, the singular set is one-dimensional.

However, in other ways the result is probably not optimal. In particular, the result should hold without the mean convexity hypothesis, and singularities should occur at only finitely many times. Indeed, I would conjecture that at each time, the surface can only have finitely many singularities unless one or more connected components have collapsed to curves. That is, the surface should consist of finitely many connected components, each of which either is a curve or has only finitely many singularities.

## 6. Nature of mean-convex singularities

Recall that when a mean convex surface evolves, it starts moving inward. Because mean convexity is preserved, it must continue to move inward. Consequently, the surface at any time lies strictly inside the region bounded by the surface at any previous time. Since the motion is continuous and since the surface collapses in a finite time, this implies that region $`\mathrm{\Omega }`$ bounded by $`M_0`$ is the disjoint union of the $`M_t`$’s for $`t>0`$. It is convenient and suggestive to speak of the $`M_t`$’s forming a foliation of $`\mathrm{\Omega }`$, although it is not quite a foliation in the usual sense because some of the leaves are singular. Figure 4 shows the foliation when the initial surface is two spheres joined by a thin tube. (The entire foliation is rotationally symmetric about an axis, so suffices to show the intersection of the foliation with a plane containing that axis.)

Figure 4: Foliation by evolving mean-convex surfaces

Theorem 8. \[W2\] Consider a mean convex surface $`M_t`$ ($`t0`$) in $`^3`$ evolving by mean curvature flow. Let $`p`$ be any singular point in the region $`\mathrm{\Omega }`$ bounded by the initial surface. If we dilate about $`p`$ by a factor $`\lambda `$ and then then let $`\lambda \mathrm{}`$, the dilated foliation must converge subsequentially to a foliation of $`^3`$ consisting of

Let us call such a subsequential limit a tangent foliation at $`p`$. The tangent foliation consists of parallel planes if and only if $`p`$ is a regular point (i.e., $`p`$ has a neighborhood $`U`$ such that the $`M_tU`$ smoothly foliate $`U`$.) A tangent foliation of concentric spheres corresponds to $`M_t`$ (or a component of $`M_t`$) becoming convex and collapsing as in Huisken’s theorem (theorem 4) to the round point $`p`$.

In figure 4, the tangent foliation at the “neck pinch” point $`A`$ is a foliation by coaxial cylinders. The tangent foliations at points $`B`$ and $`C`$ are by concentric spheres. All other points are regular points and thus give rise to tangent foliations consisting of parallel planes.

Incidentally, in cases rm (1) and rm (2), the tangent foliation is unique. That is, we have convergence and not just subsequential convergence in the statement of theorem 8. However, this is not known in case rm (3). If one tangent foliation at $`p`$ consists of cylinders, then so does any other tangent folation at $`p`$ \[S\]. But it is conceivable that different sequences of $`\lambda `$’s tending to infinity could give rise to cylindrical foliations with different axes of rotational symmetries. Whether this can actually happen is a major unsolved problem, exactly analogous to the long-standing uniqueness of tangent cone problem is minimal surface theory.

Although the tangent foliation at a singular point carries much information about the singularity, there are features that it misses. For example, consider the neck pinch in figure 4, located at the point A. At a time just after the neck pinch, the two points $`D`$ and $`E`$ on the surface that are nearest to $`A`$ have very large mean curvature and are therefore moving away from $`A`$ very rapidly. However, such behavior cannot be seen in tangent foliations: the tangent foliation at any point near $`A`$ consists of parallel planes, and the tangent foliation at $`A`$ consists of coaxial cylinders.

To capture behavior such as the rapid motion away from $`A`$, rather than dilating about a fixed point as in theorem 8, one needs to track a moving point.

Theorem 9. Consider a mean convex surface $`M_t`$ ($`t0`$) in $`^3`$ evolving by mean curvature flow. Let $`p_i`$ be a sequence of points converging to a point $`p`$ in the region bounded by $`M_0`$, and let $`\lambda _i`$ be a sequence of numbers tending to infinity. Translate the $`M_t`$’s by $`p_i`$ and then dilate by $`\lambda _i`$. Then the resulting sequence of foliations must converge subsequentially to a foliation of $`^3`$ by one of the following:

The convergence is locally smooth away from the limit foliation’s singular set (a point in case rm (1), a line in case rm (2), and the empty set in the other two cases.)

A foliation obtained in this way is called a blow-up foliation at $`p`$. Of course if all the $`p_i`$’s are equal to $`p`$, we get a tangent foliation at $`p`$.

It follows from the smooth convergence that a blow-up foliation is invariant under mean-curvature flow. That is, if we let a leaf flow for a time $`t`$, the result will still be a leaf. Consequently, except for the case rm (3) of parallel planes (which do not move under the flow), given any two leaves, one will flow to the other in finite time. Thus we can index the leaves as $`M_t^{}`$ in such a way that when $`M_t^{}`$ flows for a time $`s`$, it becomes $`M_{t+s}^{}`$. In the cases rm (1) and rm (2) of compact leaves and cylindrical leaves, the indexing interval may be taken to be $`(\mathrm{},0]`$. In case rm (3), the indexing interval is $`(\mathrm{},\mathrm{})`$. Thus, except in the case of parallel planes, the blow-up foliation corresponds to a semi-eternal or eternal flow of convex sets that sweep out all of $`^3`$.

As pointed out earlier, blow-up foliations rm (1), rm (2), and rm (3) already occur as tangent foliations. (If rm (2) occurs as a tangent foliation, then the compact convex sets must all be spheres, but in general blow-up foliations, other compact convex sets might conceivably occur.)

Thus the new case is rm (4). To see how that case arises, consider in figure 4 a sequence of points $`E_i`$ on the axis of rotational symmetry that converge to the neck pinch point $`A`$. Let $`h_i`$ be the mean curvature at $`E_i`$ of the leaf of the foliation that passes through $`E_i`$.

Now if we translate the foliation by $`E_i`$ and then dilate by $`\lambda _i=h_i`$, then the dilated foliations converge to a blow-up foliation of $`^3`$ by convex non-compact surfaces. Each leaf qualitatively resembles a rotationally symmetric paraboloid $`y=x^2+z^2`$. Furthermore, the leaves are all translates of each other. In other words, if we let one of the leaves evolve by mean curvature flow, then it simply translates with constant speed.

Incidentally, the same points $`E_i`$ with different choices of $`\lambda _i`$’s can give rise to different blow-up foliations. For if the dilation factors $`\lambda _i\mathrm{}`$ quickly compared to $`h_i`$ (that is, if $`\lambda _i/h_i\mathrm{}`$), then the resulting blow-up foliation consists of parallel planes. If $`\lambda _i\mathrm{}`$ slowly compared to $`h_i`$ (so that $`\lambda _i/h_i0`$), then the resulting foliation consists of coaxial cylinders.

So what is the “right” choice of $`\lambda _i`$? In a way, it depends on what one wants to see. But this example does illustrate a general principle:

Theorem 10. \[W2\] Let $`p_iM_{t(i)}`$ be a sequence of regular points converging to a singular point $`p`$. Translate $`M_{t(i)}`$ by $`p_i`$ and the dilate by $`h_i`$ (the mean curvature of $`M_{t(i)}`$ at $`p_i`$) to get a new surface $`M_i^{}`$. Then a subsequence of the $`M_i^{}`$ will converge smoothly on bounded subset of $`^3`$ to a smooth strictly convex surface $`M^{}`$.

Of course $`M^{}`$ is one leaf of the corresponding blow-up foliation.

Theorems 8, 9, and 10 give a rather precise picture of the singular behavior, but they raise some problems that have not yet been answered:

The strongest conjecture for rm (2) is that a blow-up foliation can only consist of planes, spheres, cylinders, or the (unique) rotationally symmetric translating surfaces.

Many more eternal and semi-eternal evolutions of convex sets are known to exist. For instance, given any three positive numbers $`a`$, $`b`$, and $`c`$, there is a semi-eternal evolution of compact convex sets, each of which is symmetric about the coordinate planes and one of which cuts off segments of lengths $`a`$, $`b`$, and $`c`$ from the $`x`$, $`y`$, and $`z`$ axes, respectively. (This can be proved by a slight modification of the proof given for example 3 in the “conclusions” section of \[W2\].) The case $`a=b=c`$ is that of concentric spheres, which of course do occur as a blow-up foliation. Whether the other cases occur as blow-up foliations is not known.

A very interesting open question in this connection is: must every eternal evolution of convex sets consist of leaves that move by translating? Tom Ilmanen has recently shown that there is a one parameter family of surfaces that evolve by translation. At one extreme is the rotationally symmetric one, which does occur in blow-up foliations. At the other extreme is the Cartesian product of a certain curve with $``$:

$$\{(x,y,z):y=\mathrm{ln}\mathrm{cos}x,1<x<1\}.$$

This case does not occur in blow-up foliations (\[W2\].)

## 7. Further reading

Three distinct approaches have been very fruitful in investigating mean curvature: geometric measure theory, classical PDE, and the theory of level-set or viscosity solutions. These were pioneered in \[BK\]; \[H1\] and \[GH\]; and \[ES\] and \[CGG\] (see also \[OS\]), respectively. Surveys emphasizing the classical PDE approach may be found in \[E1\] and \[H3\]. A very readable and rather thorough introduction to the classical approach, including some new results (as well as some discussion of geometric measure theory), may be found in \[E2\]. An introduction to the geometric measure theory and viscosity approaches is included in \[I\]. See \[G\] for a more extensive introduction to the level set approach.

Theorems 7, 8, 9, and 10 about mean convex surfaces are from my papers \[W1\] and \[W2\]. These papers rely strongly on earlier work, for instance on Brakke’s regularity theorem and on Huisken’s monotonicity formula. Huisken proved theorem 8 much earlier under a hypothesis about the rate at which the curvature blows up \[H2\]. Huisken and Sinestrari \[HS 1, 2\] independently proved results very similar to theorems 8, 9, and 10, but only up to the first occurrence of singularities.

Much of the current interest in curvature flows stems from Hamilton’s spectacular work on the Ricci flow. For survey articles about Ricci flow, see \[CC\] and \[Ha\]. For discussions of some other interesting geometric flows, see the articles by Andews \[A4\] and by Bray \[BH\] in these Proceedings.

## References