A Numerical Stability Analysis of Mean Curvature Flow of Noncompact Hypersurfaces with Type-II Curvature Blowup: II

Abstract

In previous work, we have presented evidence from numerical simulations that the Type-II singularities of mean curvature flow (MCF) of rotationally symmetric, complete, noncompact embedded hypersurfaces, constructed by the second and the fourth authors of this paper, are stable. In this work, we again use numerical simulations to show that MCF subject to initial perturbations that are not rotationally symmetric behaves asymptotically like it does for rotationally symmetric perturbations. In particular, if we impose sinusoidal angular dependence on the initial embeddings, we find that for perturbations near the tip, evolutions by MCF asymptotically lose their angular dependence—becoming round—and develop Type-II bowl soliton singularities. As well, if we impose sinusoidal angular dependence on the initial embeddings for perturbations sufficiently far from the tip, the angular dependence again disappears as Type-I neckpinch singularities develop. The numerical analysis carried out in this paper is an adaptation of the “overlap” method introduced in our previous work and permits angular dependence.

Keywords:

• Mean curvature flow
• Type-II singularities
• noncompact hypersurfaces
• stability analysis without rotational symmetry
• numerical methods

2010 Mathematics Subject Classification (2010):

• 53C44
• 35K59
• 65M06
• 65D18

1 Introduction

In previous work [1], we have used numerical simulations to show that Type-II singularity formations observed in mean curvature flow (MCF) of certain noncompact rotationally-symmetric embedded hypersurfaces [2–4] are stable for small perturbations near the tip of each initial embedding, as long as rotational symmetry is retained. The question then arises whether this behavior is stable for small perturbations that are not rotationally symmetric. In this paper, we present numerical simulations which indicate that Type-II singularity behaviors are also stable for such perturbations.

One says that MCF of an embedded hypersurface develops a Type-II singularity if the supremum of the product of the time to the singularity and the norm of the second fundamental form becomes infinite. If, alternatively, for MCF of a given embedded hypersurface, this product remains finite for the lifetime of the flow, then the singularity is defined to be Type-I. While Type-I singularities are expected to predominate for MCF (the same is expected to be true for Ricci flow), the work in [2–4] shows that Type-II singularities occur for MCF solutions originating from an open set of initial embeddings within the class of rotationally-symmetric noncompact hypersurfaces satisfying certain conditions. In these solutions, the maximum of the norm of the second fundamental form occurs at the tip (the left-most point on the hypersurface depicted in Figure 1), with the rate of curvature blowup at the tip consistent with the definition of a Type-II singularity. In [1], we consider numerical simulations of MCF originating from rotationally-symmetric embedded hypersurfaces of this type having “dimple” perturbations of two classes.

Figure 1: A MCF solution forming a degenerate neckpinch at spatial infinity with the curvature at the tip (on the left) blowing up at a Type-II rate [2–3].

The “Near Class” imposes a small-rotationally-symmetric dimple very close to the tip of each initial surface. The numerical simulations in [1] show that for Near Class initial data, these dimples disappear as MCF progresses, and Type-II behaviors are found to occur. For solutions originating from Near Class perturbations, as from unperturbed initial data, dilations near the tip approach “bowl solitons” (see Figure 2).

Figure 2: Numerical simulation of a MCF solution in the Near Class. Rotating each colored curve around the z-axis generates the hypersurface at the corresponding time. (The coordinates z and r are defined on page 3 below.) As the Type-II singularity develops, the dimple disappears.

The “Far Class” imposes a small rotationally-symmetric dimple in each initial surface relatively far from its tip. For this class of solutions, MCF evolutions develop “neckpinch” singularities that exhibit Type-I behaviors at the neck, with the curvatures at the tip remaining bounded (see Figure 3).

Figure 3: Numerical simulation of a MCF solution in the Far Class. Rotating each colored curve around the z-axis generates the hypersurface at the corresponding time. As the Type-I singularity develops, the dimple becomes a neckpinch.

As noted above, the main goal of this paper is to determine whether or not the results observed in the numerical simulations of [1] remain true if we remove rotational symmetry. We have chosen to explore this issue by carrying out numerical simulations of MCF originating from initial data much like that of [1], but with sinusoidal angular dependence initially imposed. As we show below in Section 4, the asymptotic behavior of these solutions is qualitatively unaffected by this angular dependence. More precisely, the angular dependence disappears as time approaches the first singular time, and the singular behaviors approach those observed in either the rotationally-symmetric Near Class or Far Class, independent of the initially angular dependence.

While the need to work in a collection of independent but intersecting coordinate patches is very familiar in geometric analysis on manifolds, proceeding in this way is less familiar in numerical simulations of geometric evolutions. As discussed in [1], the numerical simulations of MCF in that work require the use of two coordinate patches. Using the coordinate z for the direction parallel to the cylinder enveloping the embedded initial hypersurfaces, and the coordinate r for the radial direction orthogonal to z in the patch relatively close to the tip, the MCF evolutions are carried out in terms of the function $z(r,t)$. In the patch further away from the tip, this evolution is carried out in terms of the function $r(z,t)$. Our numerical simulations in [1] function with these two patches overlapping, and have required the development of novel numerical techniques that can handle this overlap.

In this work, we again focus on noncompact surfaces embedded in ${\mathbb{R}}^3$ and again work in two overlapping patches, here using $\theta$ as the angular coordinate. In the patch near the tip, a function $z(r,\theta ,t)$ represents the MCF evolution, while in the patch farther from the tip, we use the function $r(z,\theta ,t)$ instead. The numerical techniques that we use here to handle MCF evolutions in the two overlapping patches are adapted from those in [1] to allow for this angular dependence. (A numerical analysis of non-rotational perturbations of higher-dimensional embedded hypersurfaces will appear in a forthcoming work.)

In [1], we have speculated that for rotationally-symmetric initial data at the transition between the Near Class and the Far Class, MCF is likely to exhibit critical behavior, with Type-II degenerate-neckpinch singularity behavior occurring at the tip, and with curvature blowup rates consistent with such behavior. The asymptotic disappearance of angular dependence observed in the numerical simulations that we discuss in this paper suggests that the existence of asymptotically rotationally-symmetric critical behaviors of MCF solutions is likely to occur for such non-rotationally-symmetric transitional initial embeddings as well.

Because the numerical simulations carried out here (as well as in [1]) depend very much on the rigorous analysis of unperturbed MCF discussed in [2] and [3], we briefly review some results of those papers in Section 2 and there set up the equations for the non-rotationally-symmetric case. In Section 3, we write down the evolution equations for MCF of such non-rotationally-symmetric surfaces and discuss the computational techniques we use to carry out the numerical simulations done here. In Section 4, we present the results of these simulations, both for the Near Class and the Far Class of MCF originating from non-rotationally-symmetric initial data. The conclusions are summarized in Section 5.

2 MCF of unperturbed initial data

We briefly review some background analysis useful in studying the MCF solutions that are shown in [2–3] to develop Type-II singularities. These flows evolve from what we call our “unperturbed initial data”—complete rotationally-symmetric “undimpled” surfaces as shown in Figure 1.

For any point $(x_0, x_1, x_2)\in {\mathbb{R}}^3$, we write $z=x_0$ and $r=\sqrt{x_1^2+x_2^2}$. Let $\Gamma$ be a noncompact surface obtained by rotating the graph of $r(z)$, $a\le z<\infty$, around the z-axis. The function $r(z)$ is assumed to be strictly concave, so that $\Gamma$ is strictly convex and r is strictly increasing with $r(a)=0$ and $\underset{z\nearrow \infty }{\lim }r(z)=r_0$, where $r_0$ is the radius of the enveloping cylinder. The function r is assumed to be smooth except at $z=a$. The non-smoothness of r is due to the particular choice of coordinates; if the time-dependent profile function $r(z,t)$ were inverted in a suitable way, then this irregularity would be removed. We call the point where $r=0$ the tip of the hypersurface.

Consider the class $\mathcal{G}$ of complete surfaces that are rotationally symmetric, (strictly) convex,¹ smooth graphs over a ball and are asymptotic to a cylinder. Then a MCF solution starting from any surface in this class remains embedded under the flow and moves toward its open end while remaining asymptotic to a shrinking cylinder [5]. The solution disappears at spatial infinity at the same finite time T, called the vanishing time, that the cylinder collapses. This phenomenon can be viewed as a MCF singularity forming in finite time at spatial infinity. In fact, as shown in [2–3], this singularity is Type-II: for each real number $\gamma \ge 1/2$, there exist MCF solutions in the class $\mathcal{G}$ with the Type-II curvature blowup rate ${(T-t)}^{-(\gamma +1/2)}$ forming at the tip. For these solutions, the geometry near the tip is modeled by a bowl soliton, and the geometry near spatial infinity approaches a cylinder. Our analysis proceeds in two cases: $\gamma >\frac{1}{2}$ and $\gamma =\frac{1}{2}$.

If $\gamma >1/2$, we recall from [2] that it is useful to implement rescaled time and space variables defined by $$\tau =-\hspace{0.17em}\log \hspace{0.17em}(T-t),\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}y=z{(T-t)}^{\gamma -1/2},\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\varphi (y,\tau )=r(x,t){(T-t)}^{-1/2}.$$

Let $\lambda :=-1/y$. Then the MCF equation of a rotationally symmetric surface in terms of $\lambda (\varphi ,\tau )$ is (2.1) $$\begin{array}{cc}{{\partial }_{\tau }|}_{\varphi }\lambda \hfill & =\frac{{\lambda }_{\varphi \varphi }-2{\lambda }_{\varphi }^2/\lambda }{1+e^{2\gamma \tau }{\lambda }_{\varphi }^2/{\lambda }^4}+\left(\frac{1}{\varphi }-\frac{\varphi }{2}\right){\lambda }_{\varphi }+\left(\gamma -\frac{1}{2}\right)\lambda ,\hfill \end{array}$$(2.1) where the notation ${{\partial }_{\tau }|}_{\varphi }$ means taking the partial $\tau$-derivative while keeping $\varphi$ fixed, $\varphi \in (-\sqrt{2},\sqrt{2})$, and $\tau \ge {\tau }_0:=-\hspace{0.17em}\log \hspace{0.17em}(T-t_0)$ for an initial time $t_0$. Taking any real number $c>0$ (e.g., $c=1$) and letting $A:=c{2}^{\gamma -1/2}$, we define global initial data for (2.1) by smoothing the piecewise-smooth profile function (2.2) $$\begin{array}{c}{\widehat{\lambda }}_0(\varphi ):=\{\begin{array}{cc}\begin{array}{c}-A+e^{-2\gamma {\tau }_0}F(\zeta )-e^{-2\gamma {\tau }_0}F(R_1)\hfill \\ \hspace{1em}+\left[A-c{\left(2-{(R_1{e}^{-\gamma {\tau }_0})}^2\right)}^{\gamma -1/2}\right]\hfill \end{array},& 0\le \left|\zeta \right|\le R_1,\\ & \\ -c{(2-{\varphi }^2)}^{\gamma -1/2},& R_1{e}^{-\gamma {\tau }_0}\le \left|\varphi \right|<\sqrt{2},\end{array}\hfill \end{array}$$(2.2) where $\zeta =\varphi e^{\gamma \tau }$, $R_1$ is some large constant, and F is the unique solution of the ode initial value problem $$\frac{F_{\zeta \zeta }}{1+F_{\zeta }^2/A^4}+F_{\zeta }/\zeta =(\gamma -1/2)A,\hspace{1em}F(0)=F_{\zeta }(0)=0.$$

If $\gamma =1/2$, we recall from [3] that it is useful to implement rescaled variables defined by $$\begin{array}{cc}\tau \hfill & =-\hspace{0.17em}\log \hspace{0.17em}(T-t),\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}y=z+a\hspace{0.17em}\log \hspace{0.17em}(T-t),\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\varphi (y,\tau )=r(x,t){(T-t)}^{-1/2}.\hfill \end{array}$$

Let $\lambda :=-1/y$. Then the MCF equation of a rotationally symmetric surface in terms of $\lambda (\varphi ,\tau )$ is (2.3) $$\begin{array}{cc}{{\partial }_{\tau }|}_{\varphi }\lambda \hfill & =\frac{{\lambda }_{\varphi \varphi }-2{\lambda }_{\varphi }^2/\lambda }{1+e^{\tau }{\lambda }_{\varphi }^2/{\lambda }^4}+\left(\frac{1}{\varphi }-\frac{\varphi }{2}\right){\lambda }_{\varphi }-a{\lambda }^2,\hfill \end{array}$$(2.3) where $\varphi \in (-\sqrt{2},\sqrt{2})$, and $\tau \ge {\tau }_0:=-\hspace{0.17em}\log \hspace{0.17em}(T-t_0)$ for an initial time $t_0$. Taking any real number $c>0$ (e.g., $c=1$) and letting $A:=c{2}^{\gamma -1/2}$, we define global initial data for (2.3) by smoothing the piecewise-smooth profile function (2.4) $${\widehat{\lambda }}_0(\varphi ):=\{\begin{array}{cc}-A+e^{-{\tau }_0}F(\zeta )-e^{-{\tau }_0}F(R_1)& \\ \hspace{1em}+A-{\left(c-a\hspace{0.17em}\log \hspace{0.17em}(2-R_1^2{e}^{-{\tau }_0})\right)}^{-1},& 0⩽\left|\zeta \right|⩽R_1,\\ & \\ -1/\left(c-a\hspace{0.17em}\log \hspace{0.17em}(2-{\varphi }^2)\right),& R_1{e}^{-{\tau }_0/2}⩽\left|\varphi \right|<\sqrt{2},\end{array}$$(2.4) where $\zeta =\varphi e^{\gamma \tau }$, $R_1$ is some large constant, and F is the unique solution of the ode initial value problem $$\frac{F_{\zeta \zeta }}{1+F_{\zeta }^2/A^4}+F_{\zeta }/\zeta =a{A}^2,\hspace{1em}F(0)=F_{\zeta }(0)=0.$$

In both cases, the unperturbed initial surfaces are determined by the constants $(n,c,R_1,{\tau }_0)$. The initial data profile function ${\widehat{\lambda }}_0$ defined in either (2.2) or (2.4) is piecewise-smooth and continuous, which is sufficient for the numerical simulations based on the finite-differencing method that we use in this paper. In Section 3, we explain how to specify the perturbed initial data used in our simulations.

3 Numerical method

Our numerical method essentially consists of writing the mean curvature flow PDE in a standard parabolic form and then using a standard finite-difference method for parabolic equations. However, for reasons to be explained below, we find it convenient to write the equation in two different ways and to use both of them in implementing our numerical simulations.

In cylindrical coordinates, the surface is given by specifying r as a function of z and $\theta$. Under mean curvature flow, the function $r(t,z,\theta )$ satisfies the following evolution equation: (3.1) $${\partial }_{t}r=\frac{(1+r_z^2)r_{\theta \theta }+(r^2+r_{\theta }^2)r_{zz}-2r_{\theta }{r}_z{r}_{\theta z}-r^{-1}{r}_{\theta }^2}{r_{\theta }^2+r^2(1+r_z^2)}-\frac{1}{r}.$$(3.1)

Equation (3.1)(3.1) $${\partial }_{t}r=\frac{(1+r_z^2)r_{\theta \theta }+(r^2+r_{\theta }^2)r_{zz}-2r_{\theta }{r}_z{r}_{\theta z}-r^{-1}{r}_{\theta }^2}{r_{\theta }^2+r^2(1+r_z^2)}-\frac{1}{r}.$$(3.1) becomes singular where $r=0$. In our previous paper treating rotationally-symmetric surfaces, we have dealt with this issue by treating a part of the surface near $r=0$ by instead specifying z as a function of r. The generalization of this technique to the case of non-rotationally symmetric surfaces would be to specify z as a function of r and $\theta$. However, that would not be suitable, because the $(r,\theta )$ coordinate system becomes singular at $r=0$. Instead, in a neighborhood of the tip, we use a Cartesian coordinate system, and specify the surface by writing z as a function of x and y. Under mean curvature flow, the function $z(t,x,y)$ satisfies the following evolution equation: (3.2) $${\partial }_{t}z=\frac{(1+z_y^2)z_{xx}+(1+z_x^2)z_{yy}-2z_x{z}_y{z}_{xy}}{1+z_x^2+z_y^2}.$$(3.2)

The numerical evolution is done using Euler’s method: that is, given the surface at t, we use equation (3.1)(3.1) $${\partial }_{t}r=\frac{(1+r_z^2)r_{\theta \theta }+(r^2+r_{\theta }^2)r_{zz}-2r_{\theta }{r}_z{r}_{\theta z}-r^{-1}{r}_{\theta }^2}{r_{\theta }^2+r^2(1+r_z^2)}-\frac{1}{r}.$$(3.1) to evolve r to time $t+\Delta t$ by $r(t+\Delta t)=r(t)+\Delta t{\partial }_{t}r$. Here ${\partial }_{t}r$ represents the right-hand side of equation (3.1)(3.1) $${\partial }_{t}r=\frac{(1+r_z^2)r_{\theta \theta }+(r^2+r_{\theta }^2)r_{zz}-2r_{\theta }{r}_z{r}_{\theta z}-r^{-1}{r}_{\theta }^2}{r_{\theta }^2+r^2(1+r_z^2)}-\frac{1}{r}.$$(3.1) evaluated at time t by using the finite-difference approximations for spatial derivatives given below. Correspondingly, Euler’s method is applied to equation (3.2)(3.2) $${\partial }_{t}z=\frac{(1+z_y^2)z_{xx}+(1+z_x^2)z_{yy}-2z_x{z}_y{z}_{xy}}{1+z_x^2+z_y^2}.$$(3.2) to evolve z from time t to time $t+\Delta t$.

Our finite-difference approximations are standard second-order centered differences. That is, for the function $r(z,\theta )$, we let $r_k^i$ denote $r(i\Delta z,k\Delta \theta )$. Then our finite-difference approximations are (3.3) $$\begin{array}{cc}{r}_z\hfill & =\frac{r_k^{i+1}-r_k^{i-1}}{2\Delta z},\hfill \end{array}$$(3.3) (3.4) $$\begin{array}{cc}{r}_{\theta }\hfill & =\frac{r_{k+1}^i-r_{k-1}^i}{2\Delta \theta },\hfill \end{array}$$(3.4) (3.5) $$\begin{array}{cc}{r}_{zz}\hfill & =\frac{r_k^{i+1}+r_k^{i-1}-2r_k^i}{\Delta z^2},\hfill \end{array}$$(3.5) (3.6) $$\begin{array}{cc}{r}_{\theta \theta }\hfill & =\frac{r_{k+1}^i+r_{k-1}^i-2r_k^i}{\Delta {\theta }^2},\hfill \end{array}$$(3.6) (3.7) $$\begin{array}{cc}{r}_{z\theta }\hfill & =\frac{(r_{k+1}^{i+1}+r_{k-1}^{i-1})-(r_{k-1}^{i+1}+r_{k+1}^{i-1})}{4\Delta z\Delta \theta }.\hfill \end{array}$$(3.7)

Similarly, for the function $z(x,y)$, we let $z_k^i$ denote $z(i\Delta x,k\Delta y)$. Then our finite-difference approximations are (3.8) $$\begin{array}{cc}{z}_x\hfill & =\frac{z_k^{i+1}-z_k^{i-1}}{2\Delta x},\hfill \end{array}$$(3.8) (3.9) $$\begin{array}{cc}{z}_y\hfill & =\frac{z_{k+1}^i-z_{k-1}^i}{2\Delta y},\hfill \end{array}$$(3.9) (3.10) $$\begin{array}{cc}{z}_{xx}\hfill & =\frac{z_k^{i+1}+z_k^{i-1}-2z_k^i}{\Delta x^2},\hfill \end{array}$$(3.10) (3.11) $$\begin{array}{cc}{z}_{yy}\hfill & =\frac{z_{k+1}^i+z_{k-1}^i-2z_k^i}{\Delta y^2},\hfill \end{array}$$(3.11) (3.12) $$\begin{array}{cc}{z}_{xy}\hfill & =\frac{(z_{k+1}^{i+1}+z_{k-1}^{i-1})-(z_{k-1}^{i+1}+z_{k+1}^{i-1})}{4\Delta x\Delta y}.\hfill \end{array}$$(3.12)

Euler’s method is slow for parabolic equations, because numerical stability requires that, in Cartesian coordinates, the time step $\Delta t$ be smaller than a number of order 1 times the smaller of the squares of the space steps $\Delta x^2$ and $\Delta y^2$. In $(z,\theta )$ coordinates, the bounds on permissible $\Delta t$ are given by $\Delta z^2$ and ${(r\Delta \theta )}^2$. This makes the evolution quite slow when we approach a neck pinch, since $r\to 0$ there. If we were to use $(r,\theta )$ coordinates rather than $(x,y)$ coordinates, the bounds on permissible $\Delta t$ would be $\Delta r^2$ and ${(r\Delta \theta )}^2$. Since $r\to 0$ at the origin in polar coordinates, this would give rise to a very small $\Delta t$. This is another way of saying why Cartesian coordinates are preferable to polar coordinates in this case.

Our overlap method requires that we use two grids, each of which corresponds to a manifold with boundary. The boundary points of one grid correspond to interior points of the other grid, and their values are set using interpolation in the other grid. In our previous paper, one grid has a single coordinate r and a single boundary point given by $r_{\max }$ while the other grid has a single coordinate z with a single boundary point $z_{\min }$. Here the boundary of our $(z,\theta )$ grid is the circle $(z_{\min },\theta )$, while the boundary of our $(x,y)$ grid is a square of side length L with center at the origin. As with our previous method, as the surface evolves, we adjust the parameters (in this case $z_{\min }$ and L) so that each boundary remains within the interior of the other grid.

At each time step the function at the interior points of each grid are updated using (the finite difference versions of) equations (3.1)(3.1) $${\partial }_{t}r=\frac{(1+r_z^2)r_{\theta \theta }+(r^2+r_{\theta }^2)r_{zz}-2r_{\theta }{r}_z{r}_{\theta z}-r^{-1}{r}_{\theta }^2}{r_{\theta }^2+r^2(1+r_z^2)}-\frac{1}{r}.$$(3.1) and (3.2)(3.2) $${\partial }_{t}z=\frac{(1+z_y^2)z_{xx}+(1+z_x^2)z_{yy}-2z_x{z}_y{z}_{xy}}{1+z_x^2+z_y^2}.$$(3.2) . The function at the boundary points of each grid are then updated as follows: Because the grids correspond to overlappling coordinate systems, the boundary points of each grid should be thought of as contained in the interior of the other grid. The function at each boundary point is updated by linear interpolation from the (updated) nearest neighbor points of the other grid. Note that all this is done in a single pass, meaning that we actually have a simple algorithm in which all non-boundary points are updated (i.e., taken to the next time step) using the equation of motion and then all boundary points are updated by linear interpolation in the other grid.

To understand the nature of this algorithm, it is helpful to think in analogy to what is done in parallel computing for PDE. Recall that in parallel computing the grid is divided among different processors. But then to compute finite differences at all its points each processor needs some information contained in the other processors, and it must obtain that information from them. Similarly, the only thing that each of our grids uses its boundary points for is to compute finite differences at its interior points. But each of those boundary points should be thought of as living in the interior of the other grid. So each grid must query the other grid (through linear interpolation) to obtain the needed information about its boundary points.

We end this section by explaining how we choose the perturbed initial data sets for the numerical simulations. We first recall that for an unperturbed solution constructed in [2] (or [3]), respectively, its initial data set is obtained by joining a scaled bowl soliton to a cylinder at spatial infinity, and is defined precisely in the rescaled coordinates by equation (2.2)(2.2) $$\begin{array}{c}{\widehat{\lambda }}_0(\varphi ):=\{\begin{array}{cc}\begin{array}{c}-A+e^{-2\gamma {\tau }_0}F(\zeta )-e^{-2\gamma {\tau }_0}F(R_1)\hfill \\ \hspace{1em}+\left[A-c{\left(2-{(R_1{e}^{-\gamma {\tau }_0})}^2\right)}^{\gamma -1/2}\right]\hfill \end{array},& 0\le \left|\zeta \right|\le R_1,\\ & \\ -c{(2-{\varphi }^2)}^{\gamma -1/2},& R_1{e}^{-\gamma {\tau }_0}\le \left|\varphi \right|<\sqrt{2},\end{array}\hfill \end{array}$$(2.2) (or (2.4), respectively). Expressed in terms of the unscaled $(z,r)$-coordinates, the ODE for the bowl soliton takes the form (3.13) $$z_{rr}=(1+z_r^2)\left(\beta -\frac{z_r}{r}\right).$$(3.13)

For each choice of $\gamma$, we fix $c=1$ and choose ${\tau }_0$ large (in all the numerical simulations $\gamma =3/4$ and ${\tau }_0=4$) so that the matched asymptotics explained in Section 2 make sense. We then choose $\beta$ according to $\beta =c^{-1}(\gamma -1/2)2^{-(\gamma -1/2)}{e}^{-(\gamma +1/2){\tau }_0}$, and numerically integrate equation (3.13)(3.13) $$z_{rr}=(1+z_r^2)\left(\beta -\frac{z_r}{r}\right).$$(3.13) outward from $r=0$ to value $r_1$, which is obtained by writing $R_1$ (chosen to be $e^{\gamma {\tau }_0/2}$ in our numerical simulations) in the unscaled r-coordinate. For all $r>r_1$, we use the following analytic formula which follows from rewriting the second equation in (2.2) in the $(r,z)$-coordinates: (3.14) $$z(r)=z(r_1)+\frac{1}{c}\left[{\left(2e^{-{\tau }_0}-r^2\right)}^{\frac{1}{2}-\gamma }-{\left(2e^{-{\tau }_0}-r_1^2\right)}^{\frac{1}{2}-\gamma }\right].$$(3.14)

We choose the patch where we write z as a function of $(x,y)$ to be a square whose sides have length $10r_1$. In this patch, we have $r=\sqrt{x^2+y^2}$. For $r<r_1$, we specify z by the results of the numerical integration of (3.13). For $r\ge r_1$, z is given by formula (3.14).

In the patch where r is written as a function of $(z,\theta )$, we use the formula obtained by inverting expression (3.14).

For our “Near Class” of distorted surfaces, we change the surface only in the $(x,y)$ patch. We choose a number $r_m$ and distort the surface only where $r<r_m$. In this region, we choose two amplitudes $a_0$ and $a_1$ and change z by (3.15) $$z\mapsto z+a_0(1+a_{1}xy)(r^2-r_m^2)$$(3.15)

For our “Far Class” of distorted surfaces, we change the surface only in the $(z,\theta )$ patch. We choose an amplitude $a_0$, and two values of z: $z_a$ and $z_b$. We distort the surface only on the interval $(z_a,z_b)$. On this interval, we define the quantity F by (3.16) $$F\equiv \left(1+\frac{1}{4}\cos \hspace{0.17em}(n\theta )\right)\hspace{0.17em}\sin \hspace{0.17em}\left(\frac{\pi (z-z_b)}{z_a-z_b}\right).$$(3.16)

We then distort the surface by (3.17) $$r\mapsto r(1-a_0{F}^2).$$(3.17)

Here n is an integer, which we choose to be even so that the neck pinch will occur at $r=0$. An odd n would give rise to an off-center neck pinch, which results in a coordinate singularity when the shrinking surface intersects $r=0$.

4 Numerical results

We present our results primarily graphically. In the first pair of graphs—Figures 4 and 5—we represent the initial hypersurface in a neighborhood of the tip. The angular dependent dimple is represented as a simulated three-dimensional graph in Figure 4 and as a color coded “heat map” in Figure 5. The “rounding” of the mean curvature flow evolving from the initial hypersurface of Figure 5 is in a heat map at a later time in Figure 6, and then even later in Figure 7, when the hypersurface becomes almost entirely rotationally symmetric.

Figure 4: A simulated three-dimensional graph of angular dependent initial embedding for Near Class initial data in a neighborhood of the tip.

Figure 5: A colored heat map of the initial data represented in Figure 4. The colors correlate with the graph height.

Figure 6: A colored heat map of the mean curvature flow of the initial data depicted in Figure 5. It is clearly becoming rounder.

Figure 7: A colored heat map of the mean curvature flow with initial data from Figure 5 at a later time. This time, the angular dependence is gone.

We do not also include graphs of the evolution of the curvature at the tip for angular dependent Near Class data, but we do note that such graphs (similar to sub-figures (B) and (C) in Figures 2–5 in [1]) do show that, for such mean curvature flows, the numerical evidence supports the contention that Type-II singularities occur at the tip for perturbed initial data, even with angulardependence.

In Figures 8–10, we graphically represent the mean curvature flow for initial data with $\hspace{0.17em}\cos \hspace{0.17em}(2\theta )$ angular dependence. In Figure 8, each of the four colored curves represents a vertical cross-section of the mean curvature flow at the location of the evolving neck pinch for Far Class initial data at four successive times. The loss of the angular dependence is evident. In Figure 9, a similar graphical representation of vertical cross-sections of the mean curvature flow is depicted a bit to the left of the neck pinch, while in Figure 10, the same is depicted for vertical cross-sections a bit to the right of the neck pinch. In all of these figures, the roundingis evident.

Figure 8: Vertical cross-sections of the mean curvature flow for Far Class initial data at successive times, with $\hspace{0.17em}\cos \hspace{0.17em}(2\theta )$ angular dependence. These cross-sections are chosen to be the location of the developing neck pinch. The loss of the angular dependence is evident.

Figure 9: Vertical cross-sections of the MCF for the same initial data in Figure 8, but with the cross-sections obtained to the left of the neck pinch.

Figure 10: Vertical cross-sections of the MCF for the same initial data in Figure 8, with the cross-sections obtained to the right of the neck pinch.

Figures 11–13 are very similar to Figures 8–10, in their depiction of vertical cross-sections of the mean curvature flow for Far Class initial data, but in this case, the imposed angular dependence is $\hspace{0.17em}\cos \hspace{0.17em}(4\theta )$. Consequently, the initial data shows larger angular dependence. But again, we see that as the flow evolves, these cross-sections become increasingly round.

Figure 11: Vertical cross-sections similar to Figure 8, but with imposed angular dependence of $\hspace{0.17em}\cos \hspace{0.17em}(4\theta )$.

Figure 12: Vertical cross-sections similar to Figure 9, with imposed angular dependence $\hspace{0.17em}\cos \hspace{0.17em}(4\theta )$.

Figure 13: Vertical cross-sections similar to Figure 10, with imposed angular dependence $\hspace{0.17em}\cos \hspace{0.17em}(4\theta )$.

5 Conclusions

These numerical simulations indicate that mean curvature evolutions of angular dependent perturbations asymptotically approach rotationally symmetric solutions. They do not prove that the mean curvature flow for general perturbations of the initial embedded hypersurfaces considered in [2–3] always have the same Type-II singularity behavior at the tip as is found in [2–3]. However, the numerical simulations carried out here together with those in [1] provide strong evidence that this singularity behavior is indeed stable. Mathematical proof of such stability deserves exploration.

Additional information

Funding

DG thanks the National Science Foundation for support in PHY-1806219 and PHY-2102914. JI thanks the National Science Foundation for support in grant PHY-1707427. DK thanks the Simons Foundation for support in Award 635293. HW thanks the Australian Research Council for support in DE180101348.

Notes

1 Throughout this paper, ``convex'' means ``strictly convex''.