Double Bubbles Minimize

Joel Hass^{∗}and Roger Schlafly

August 2, 1999

Abstract The classical isoperimetric inequality in ℝ^{3} states that the surface of smallest area enclosing a given volume is a sphere. We show that the least area surface enclosing two equal volumes is a double bubble, a surface made of two pieces of round spheres separated by a flat disk, meeting along a single circle at an angle of 120^{∘}.

1 Introduction

In this paper we find the unique surface of smallest area enclosing two equal volumes. The surface is called a double bubble, and is made of two pieces of round spheres separated by a disk, meeting along a single circle at an angle of 120^{∘}. This is the form assumed by two equally sized soap bubbles which are brought together until their boundaries conglomerate to form a common wall. See Figure 1, due to John Sullivan.

Isoperimetric problems, which study maximizing the size of enclosed bounded regions whose boundary is of fixed size, are among the oldest problems in mathematics. For a broad discussion of isoperimetric problems see Osserman [33].

The two volume isoperimetric problem in ℝ^{3} was considered by the Belgian physicist J. Plateau [35], and appears in C.V. Boys’ famous book on soap bubbles. As Boys wrote,

When however the bubble is not single, say two have been blown in real contact with one another, again the bubbles must together take such a form that the total surface of the two spherical segments and of the part common to both, which I shall call the interface, is the smallest possible surface which will contain the two volumes of air and keep them separate.

We now outline briefly the recent history of the problem, with more details given later. Fred Almgren proved the existence of bubble clusters enclosing a given collection of volumes in three-space, and minimizing area among a general class of surfaces. Jean Taylor established a regularity result for these bubble

∗partially supported by the National Science Foundation.

Figure 1: A double bubble enclosing two equal volumes.

clusters, showing that they satisfy regularity properties physically observed by Plateau. Brian White suggested an argument for showing that a two region minimizer must be a surface of revolution, and Joel Foisy and Michael Hutchings filled in the details of this argument. At this point it was known that the minimizer consisted of pieces of Delaunay surfaces meeting at 120^{∘}. Hutchings went on to eliminate many remaining configurations, leaving just the standard double bubble and an additional family known as torus bubbles as the potential minimizers. In this paper we closely examine the family of torus bubbles and show they cannot be minimizers for the equal volume double bubble problem. Each torus bubble is determined by an ODE, and we can rigorously eliminate the entire collection by an examination of the solutions of these ODEs using a computational technique called interval arithmetic, described in Section 5. At the end, only the standard double bubble remains as a possible minimizer.

Extensive interest in the problem was generated in recent years by work of Frank Morgan. The planar case has been solved in [1] by methods special to two dimensions. In higher dimensions it is not obvious that a solution to the problem even exists. A fundamental paper of Almgren [3] established the existence of solutions to a great variety of geometric minimization problems, including multiple component isoperimetric problems in ℝ^{3}. Despite the fact that existence was established some time ago, our result is the first explicit example of a surface in ℝ^{3} solving a multiple region finite volume isoperimetric problem. Lawlor and Morgan have shown that the cone over the regular tetrahedron is a minimizing surface among all those with its boundary that separate the four faces of the tetrahedron, but it is still not known whether this surface is minimizing without the separation condition.

Such multiple region isoperimetric problems arise in many fields, including the growth and shape of biological cells [28], [40]. There were extensive studies of such problems by physical and biological scientists in the 19th century [24],[35]. Plateau established experimentally that a soap bubble cluster is a piecewise-smooth surface having only two types of singularities. The first type of singularity occurs when three smooth surfaces come together along a smooth triple curve at an angle of 120^{∘}. The second type of singularity occurs when six smooth surfaces and four triple curves converge at a point. The angles are then equal to those of the cone over the 1-skeleton of a regular tetrahedron. A mathematical proof that these types of singularities are the only ones possible in a minimizing bubble in ℝ^{3} was given by Taylor [39]. However there were no explicit minimizing bubbles known for any collection of volumes exhibiting either of these singularities. Thus we have found the first explicit example of a closed minimizing surface in ℝ^{3} known to exhibit some of the singularities predicted by Plateau.

The arguments presented in this paper are a mixture of geometrical analysis and of estimates of geometric quantities obtained by the use of numerical computation. We perform these calculations with strict estimates on the accuracy of the computations. Since it is still somewhat unusual in a mathematical proof to use digital computers to do calculations involving real numbers, we will say a few words about the nature of this part of the argument.

Computers are widely known in the mathematical community to have been used successfully for analysis of discrete and combinatorial problems, such as the 4-color theorem [6], [37]. Problems such as isoperimetric inequalities are qualitatively different, since they inherently involve real numbers. Real numbers are represented in digital computers in a floating point format which allows exact description of only a finite number of rationals. Simple calculations, like division, or even addition, lead to unrepresentable numbers, and so the computer must round off. For purposes of mathematical proof, the size of the round off must be accurately tracked throughout the calculation.

Methods for strictly estimating solutions of differential equations exist, but are not yet widely used in the mathematical community. The vast majority of numerical work is for approximation and simulation that does not meet the standards of a mathematical proof. However there have been some important results achieved through rigorous use of floating point numerical methods to achieve traditional mathematical proofs, notably Lanford’s work on the Feigenbaum Conjectures [27] and Fefferman and de la Llave’s work on the stability of matter [12],[15] and the work of McKay and Percival on non-existence of invariant tori [34]. See also recent work of Hales on the Kepler Conjecture [19] and Gabai, Meyerhoff and Thurston [17] on hyperbolic 3-manifolds. Computer calculations are essential to our proof that equal volume double bubbles minimize area. There are too many calculations to be done by hand.

The proof parameterizes the space of possible solutions by a two-dimensional rectangle, one dimension corresponding to an angle and the other to a mean curvature. This rectangle is subdivided into 15,016 smaller rectangles which are investigated by calculations involving a total of 51,256 numerical integrals. Every calculation is done with strict error bounds, and all results are precise mathematical statements. All operations conform to the IEEE 754 standard for computer arithmetic, a widely adopted standardized method for implementing real number (floating point) computations on computers [5]. Our methods indicate that numerical techniques are likely to play an important role in future geometrical arguments.

The main result is the following, announced in [20]:

Theorem 1 The unique surface of least area enclosing two equal volumes in ℝ^{3} is a double bubble.

Corollary 2 For any surface in ℝ^{3} enclosing two regions, each of volume v, the area a satisfies

a^{3} ≥ 243πv^{2},

with equality if and only if the surface is a double bubble enclosing two regions of volume v.

The first issue in establishing Theorem 1 is to find an appropriate category of surface in which to minimize area.

It suffices in this paper to consider piecewise smooth two-dimensional surfaces. It is sometimes useful to consider a much more general notion of surface, such as the (F, ϵ, δ) sets described in [3]. Our arguments actually imply that the double bubble minimizes in this larger class; see Theorem 2.2. Consideration of such a larger class of surfaces is needed primarily in the establishment of the existence and regularity of a minimizer, carried out in [3] and [39], and we will not need to be overly concerned with it in this paper.

Define a piecewise-smooth curve to be an embedded finite union of smooth curves, with any two either disjoint or having intersection contained in their endpoints. Define a piecewise-smooth surface to be an embedded finite union of smooth surfaces with piecewise-smooth boundary curves, with any two surfaces either disjoint or intersecting along piecewise-smooth curves contained in their boundaries. The set of points which are not in the interior of a smooth subsurface of a piecewise-smooth surface is called the singular set.

Define a bubble to be a piecewise-smooth surface satisfying:

1. Each two dimensional surface has constant mean curvature.

2. The singular set is of the type described by Plateau. It consists of smooth triple curves along which three smooth surfaces come together at an angle of 120^{∘} and isolated vertices where six smooth surfaces and four triple curves converge at a point. The angles at the point are equal to those of the cone over the 1-skeleton of a regular tetrahedron.

3. The mean curvatures around an edge in the singular set where three surfaces have common boundary sum up to zero.

The above conditions are necessary for no local perturbation of a piecewisesmooth surface to decrease the area while preserving the volume in each of its complementary regions.

A bubble enclosing regions of prescribed volumes is called a minimizing bubble if it minimizes area among all bubbles enclosing the same volumes. The regions are not necessarily connected.

Given positive constants v_{1} and v_{2}, let a(v_{1}, v_{2}) denote the infimum of the area of piecewise-smooth surfaces enclosing two regions R_{1} and R_{2} in ℝ^{3} which are closed bounded sets with disjoint interiors such that volume(R_{1}) = v_{1} and volume(R_{2}) = v_{2}. We will refer to the interior of the complement of R_{1} ∪ R_{2} as the exterior region R_{0}. Note that R_{0} may not be connected, in which case the bubble encloses some compact “empty regions”.

A double bubble enclosing volumes v_{1}, v_{2} is a surface made of three pieces of round spheres, meeting along a single circle at an angle of 120^{∘}, and enclosing two connected regions having volumes v_{1} and v_{2}. We consider the plane to be a sphere of infinite radius in this setting, allowing the interface of a double bubble to be a flat disk, as occurs when v_{1} and v_{2} are equal. See Figure 2. In this and all other figures, we assume that the axis of rotational symmetry is the x-axis. Furthermore, we will always take the generating curves for a surface of revolution to be in the upper half of the xy-plane.

Figure 2: Equal and unequal volume double bubbles in cross-section. The axis of revolution is the x-axis.

It has been conjectured since the work of Plateau that double bubbles give the most efficient shape for enclosing two given volumes.

Conjecture 3 The double bubble enclosing volumes v_{1}, v_{2} is the unique surface having area equal to a(v_{1}, v_{2}).

Theorem 1 solves this conjecture in the case that v_{1} = v_{2}.

Deep results of Almgren and Taylor, summarized in Section 2, imply that for any two positive numbers v_{1}, v_{2} there exists a minimizing bubble S(v_{1}, v_{2}) in ℝ^{3} which encloses volumes v_{1} and v_{2}. Arguments based on ideas of B. White and F. Morgan, and developed in [16],[22],[32], show that any solution must be a surface of revolution. We will need to refer to this proof, so we present a simple version for the case of two regions in ℝ^{3} in Theorem 2.6. The lack of such an argument for isoperimetric problems involving three or more regions in ℝ^{3} makes that problem more formidable.

While the reduction to a surface of revolution gives an enormous simplification in the scale of the problem, F. Morgan has pointed out several major topological obstacles to solving the double bubble conjecture. The first is that the regions R_{1} and R_{2} bounded by the minimizing bubble S(v_{1}, v_{2}) may not be connected, as illustrated in Figure 3.

Figure 3: Cross-section of a bubble with a non-connected region. The nonconnected region is shaded.

A second problem is that S(v_{1}, v_{2}) may enclose bounded regions in ℝ^{3} which do not form part of either R_{1} or R_{2}, as in Figure 4. These regions are called empty regions by Morgan.

Figure 4: Cross-section of a bubble with an empty torus region.

A third problem is that S(v_{1}, v_{2}) may enclose non-simply connected regions. Since the surface is a surface of revolution, these regions are homeomorphic to solid tori, and we call them torus components.

There are therefore numerous possible configurations which a minimizing bubble might take. A recent breakthrough due to Hutchings [22] has given restrictions on the type of surfaces that can arise in the double bubble problem, eliminating many possibilities. Some remaining possibilities are depicted in Figure 5. For the case of equal volumes, Hutchings showed that there were further constraints. Each region must be connected, leaving only two possible configurations, the double bubble and an additional class of possibilities, called torus bubbles, whose properties will be discussed in Section 4.

In Section 5 we describe the algorithm used in a series of computations which show that torus bubbles are not minimizers for the two equal volume isoperimetric problem in ℝ^{3}.

Figure 5: Cross-sections of some possible bubble configurations.

Techniques extending those developed in this paper also prove Conjecture 3 for other volume ratios, but at this time do not seem to suffice to cover all cases.

2 Existence and regularity

Almgren showed in [3] that there exists an area minimizing surface S(v_{1}, v_{2}) in ℝ^{3} among the set of surfaces enclosing volumes v_{1}, v_{2}. Here surface refers to a generalized notion defined using the methods of geometric measure theory, specifically what Almgren calls (F, ϵ, δ) sets. For our purposes it suffices that this class includes the piecewise-smooth surfaces.

Theorem 2.1 ([3]) An area minimizing surface S(v_{1}, v_{2}) exists and is a smooth surface almost everywhere.

The nature of the singularities of S(v_{1}, v_{2}) was established by Taylor.

Theorem 2.2 ([39]) S(v_{1}, v_{2}) is a piecewise-smooth surface. Its singularities consist of smooth triple curves along which three smooth surfaces come together at an angle of 120^{∘} and isolated points where four triple curves and six pieces of surface converge. At these isolated points the asymptotic cone is the cone over the 1-skeleton of a regular tetrahedron.

Thus Taylor’s work established that S(v_{1}, v_{2}), which Almgren’s Theorem established as a minimizer among a generalized class of surfaces called (F, ϵ, δ) surfaces, is a bubble, in our terminology. It follows that if we can establish the minimality of the double bubble among our class of bubbles then we will also show it minimizes among the more general class of surfaces considered by Almgren. We now establish some properties of minimizing bubbles. Lemma 2.3, Corollary 2.4 and Lemma 2.5, are standard results in variational geometry. We include brief proofs for completeness. Lemma 2.3 states that the sum of the (oriented) mean curvatures of all points crossed by a path starting and ending in the same region is zero. It implies that the sum of the (oriented) mean curvatures of all points crossed by any path depend only on its starting and ending points. These results hold even if the path crosses through many surfaces, and if it starts and ends in different components of the same region.

Lemma 2.3 Let γ be an oriented curve in ℝ^{3} intersecting a minimizing bubble B transversely at regular points such that the initial and final points of γ lie in the interior of the same region. Then the sum of the mean curvatures of all the points of γ ∩ B, oriented by γ, is zero.

Proof: Perturb the curve slightly so that each of its intersections with B becomes perpendicular. Consider a deformation of B which pushes points in ℝ^{3} near γ a uniform distance along the curve. To first order this preserves the volume of each region. The derivative of the area, to first order, is given by the sum of the mean curvatures over the points of γ ∩ B. If this sum is non-zero, a deformation can be defined which decreases area while preserving the volume of each region. □

A special case of the above lemma occurs when γ is a simple closed curve encircling a triple curve of the bubble. The lemma then implies that the sum of the mean curvatures around the triple curve adds up to zero. This local minimization condition is built into our definition of a bubble. More generally, the lemma implies that any two surfaces separating the same pair of regions has the same mean curvature.

Corollary 2.4 If S_{1} and S_{2} are surfaces in a bubble, each of which separate regions R_{1} and R_{2}, then their mean curvatures are equal.

Proof: Consider an arc α which starts in R_{1} passes through S_{1} into R_{2}, then transversely through the bubble and on through S_{2} back into R_{1}. A subarc β of α starts and ends in R_{2}. Applying Lemma 2.3 to both α and β implies the corollary. □

Lemma 2.5 A minimizing bubble B is connected.

Proof: If not, we can move a component by an isometry of ℝ^{3} until it touches a distinct component. The resulting singularity violates those allowed in Theorem 2.2. □

The following result was first observed by B. White and F. Morgan, and first written down by Foisy in [16]. Since we will need to refer to the proof, we present an argument for the case of two regions in ℝ^{3}. More detailed and general arguments can be found in the work of Hutchings [22] and Morgan [32].

Theorem 2.6 S(v_{1}, v_{2}) is a surface of revolution.

Proof: Since S(v_{1}, v_{2}) is a minimizing bubble, each of its faces separates a pair of distinct regions - either separating R_{1} from R_{2} or one of them from R_{0} (which is not connected if there are empty regions.)

Given any unit vector Z in ℝ^{3}, there is a plane P_{Z} in ℝ^{3} perpendicular to Z which bisects the total volume v_{1} + v_{2} enclosed by S(v_{1}, v_{2}). If Q_{Z} is a parallel plane which also bisects the total volume, then P_{Z} and Q_{Z} split the regions R_{1} and R_{2} in the same proportions. (In fact P_{Z} must coincide with Q_{Z}, but we do not need to show this. It is also irrelevant to this argument whether S(v_{1}, v_{2}) contains empty regions.) Let f(Z) denote the proportion of the volume of R_{1} on the side of P_{Z} to which Z points. The function f(Z) is well defined, continuous on the unit 2-sphere, 0 ≤ f(Z) ≤ 1 and f(−Z) = 1 − f(Z). Along any great circle of directions on the unit 2-sphere, the Intermediate Value Theorem implies that there are at least two points where f(Z) = 1/2 and the plane P_{Z} bisects the volume of both regions. We fix Z to be such a vector.

Consider the intersection of S(v_{1}, v_{2}) with each of the two half-spaces determined by P_{Z}. The intersection of a face of S(v_{1}, v_{2}) with the plane P_{Z} is given by a graph in P_{Z} with smooth edges and isolated vertices, or else the face is completely contained in P_{Z}, since the face is a smooth constant mean curvature surface [18]. Reflection of the smaller area piece of S(v_{1}, v_{2}) lying in one of these half-spaces gives a new surface S_{1}(v_{1}, v_{2}) which encloses regions of the same volumes. Since S_{1}(v_{1}, v_{2}) cannot have less area, it must have the same area as S(v_{1}, v_{2}). S_{1}(v_{1}, v_{2}) has the property that reflection through P_{Z} preserves the new regions, which we continue to call R_{0}, R_{1} and R_{2}. If the plane P_{Z} contained a face of S(v_{1}, v_{2}) then this face would now separate two components of the same region in S_{1}(v_{1}, v_{2}), and area could be reduced by removing that face from S_{1}(v_{1}, v_{2}), preserving the volume of each region. So P_{Z} cannot contain a face of S(v_{1}, v_{2}).

We repeat the above argument along a great circle of the unit 2-sphere consisting of the unit vectors which are perpendicular to Z, to find a plane P_{W} perpendicular to P_{Z} which bisects the volume of the regions R_{1} and R_{2} in S_{1}(v_{1}, v_{2}). This process constructs a surface S_{2}(v_{1}, v_{2}), still having the same area as S(v_{1}, v_{2}), enclosing regions of the same volumes, and for which reflection through each of P_{Z} and P_{W} preserves R_{1} and R_{2}.

Since composing reflections through two perpendicular planes gives a rotation of angle 180^{∘}, it follows that rotation of angle 180^{∘} about L = P_{Z} ∩ P_{W} preserves the regions R_{1} and R_{2}. Now consider any plane Q containing L. Note that half of the volume of R_{1} and R_{2} lies on each side of Q. As before the plane Q cannot contain a face of S_{2}(v_{1}, v_{2}). If the intersection of Q with S_{2}(v_{1}, v_{2}) is not perpendicular to Q, then replacing half of S_{2}(v_{1}, v_{2}) with its reflected image through Q gives a new minimizing surface which is not smooth, and has singularities not allowed by Theorem 2.2. Thus S_{2}(v_{1}, v_{2}) is perpendicular to each plane through L and is therefore a surface of revolution around L. S_{1}(v_{1}, v_{2}) and S_{2}(v_{1}, v_{2}) may not be identical, but they coincide on a half-space of ℝ^{3}. Moreover this half-space was chosen arbitrarily, so that S_{1}(v_{1}, v_{2}) is cut by P_{W} into two pieces, each half of some surface of revolution. The axis of each of these surfaces of revolution is L = P_{Z} ∩ P_{W}, so they coincide and S_{1}(v_{1}, v_{2}) itself is a surface of revolution. Similar reasoning shows that S(v_{1}, v_{2}) is a surface of revolution, and coincides with S_{2}(v_{1}, v_{2}) everywhere. If we make S_{1} by reflecting the other half of S across P_{Z}, then the same argument shows that the other half of S is also a surface of revolution about a line L′ in P_{Z}. The Almgren and Taylor regularity results, and in particular the unique continuation property of constant mean curvature surfaces, imply that L′ = L, so S itself is a surface of revolution about L. □

We next summarize some key results obtained by Hutchings in [22]. Recall that S(v_{1}, v_{2}) is a minimizing bubble that separates ℝ^{3} into regions R_{1} of volume v_{1} and R_{2} of volume v_{2} and that a(v_{1}, v_{2}) is the area of S(v_{1}, v_{2}).

Theorem 2.7 The function a(v_{1}, v_{2}) is strictly concave on [0, +∞) × [0, +∞).

Corollary 2.8 If v′_{1} > v_{1} then a(v′_{1}, v_{2}) > a(v_{1}, v_{2}).

Proof: As v_{1} → +∞, a(v_{1}, v_{2}) ≥ a(v_{1}) → +∞. Concavity of the function a implies the corollary. □

Corollary 2.9 S(v_{1}, v_{2}) has no empty regions.

Proof: If there is an empty region, join it it to one of the other regions by removing a face on its boundary, and apply the previous corollary. □

Deeper connectedness results also follow from concavity. In particular, Hutchings deduced that each of R_{1} and R_{2} is connected if the volumes are equal.

Theorem 2.10 S(v_{1}, v_{1}) encloses exactly two connected components.

Proof: See Hutchings [22, Theorem 4.2]. □

3 Delaunay Surfaces

In this section, we summarize the classification theory of constant mean curvature surfaces of revolution, the Delaunay surfaces, and present some properties of these surfaces that will be used in our study of bubbles.

The mean curvature vector field along a surface is defined to be the trace of the second fundamental form of the surface. The mean curvature can be obtained by taking the inner product of the mean curvature vector field with a unit normal vector field to the surface. The mean curvature vector field gives the direction of a variation which decreases the area of a surface as quickly as possible [38].

For a surface of revolution we have simple formulas for the mean curvature. Consider a surface of revolution about the x-axis in ℝ^{3}, with generating curve contained in the upper half-plane. If the generating curve is a graph, y = y(x), then the mean curvature h, equal to the sum of the principal curvatures, is given by the formulas:

k_{m} = \frac{−ÿ}{(1 + ẏ^{2})^{3/2}} ,

k_{p} = \frac{1}{y\sqrt{1 + ẏ^{2}}} ,

h = k_{m} + k_{p},

where k_{m} is the curvature at (x, y) of the generating curve in the xy-plane (sometimes called the meridian, or profile curve) and k_{p} is the normal curvature of the parallel curve. The value of k_{p} is equal to the reciprocal of the distance to the x-axis along the perpendicular to the generating curve [38]. Note that the mean curvature of a unit sphere with these formulas equals two. For surfaces of revolution generated by curves which are not graphs, we need to specify an orientation on the generating curve to fix the sign of the mean curvature. We use the convention that the sign of the mean curvature is given by the formulas above if the curve is oriented left-to-right, or equivalently if the graph’s tangent vector has a positive x-component. If its tangent vector has a negative xcomponent, the signs in the formulas for k_{m} and k_{p} should be reversed in the above formulas. For curves with vertical tangent vectors the above formulas do not apply directly, but in our applications vertical tangents are isolated points, and the limits of the above formulas will give k_{m} and k_{p} at vertical tangencies.

Surfaces of revolution having constant mean curvature were first studied by Euler, and classified by Delaunay [13]. The term Delaunay surface is used to refer to these surfaces.

Theorem 3.1 The surfaces of revolution of constant mean curvature are the plane, sphere, catenoid, cylinder, unduloid, and nodoid. The generating curves of a nodoid and an unduloid, called the nodary and the undulary, are periodic along the x-axis, and have exactly one local minimum and one local maximum in each period. The undulary is a graph over the x-axis. The nodary has one local maximum, one local minimum and two vertical tangencies in each period.

Proof: Expositions of the classification of Delaunay surfaces can be found in [14] and [23]. □

The nodoid and unduloid are not as well known as the other surfaces. An undulary can be obtained by rolling an ellipse along the x-axis and tracing the path taken by one of the foci. A nodary can be obtained by rolling a hyperbola [14]. See Figure 6.

With their generating curves oriented left to right, the mean curvature of each of the Delaunay surfaces is non-negative, with the exception of the nodoid. The sign for the mean curvature of a nodoid depends upon the choice of orientation for its generating nodary, and can be positive or negative.

If a surface of revolution minimizes area among all surfaces of revolution surrounding a given volume, then its generating curve satisfies an associated

Figure 6: Generating curves of a nodoid and an unduloid.

Euler-Lagrange equation. This second order ODE implies that the mean curvature is constant [14].

The mean curvature of the surface of revolution generated by a C^{2} graph y(x) is given by:

h = k_{m} + k_{p} = \frac{−ÿ}{(1 + ẏ^{2})^{3/2}} + \frac{1}{y\sqrt{1 + ẏ^{2}}} .  (1)

We can integrate this ODE once, after multiplying through by 2yẏ.

0 = h + \frac{ÿ}{(1 + ẏ^{2})^{3/2}} − \frac{1}{y\sqrt{1 + ẏ^{2}}}
= 2hyẏ + \frac{2yẏÿ}{(1 + ẏ^{2})^{3/2}} − \frac{2ẏ}{\sqrt{1 + ẏ^{2}}}
= \frac{df(x)}{dx}

where f(x) is the function

f(x) = hy^{2} − \frac{2y}{\sqrt{1 + ẏ^{2}}} .  (2)

Therefore f is constant along the graph. The fact that the equation for constant mean curvature surfaces of revolution has a first integral was known to Plateau [35, pages 138-139], who in turn references work of Beer. Korevaar, Kusner and Solomon called f the force of a constant mean curvature surface, and pointed out that it has a physical interpretation as the net force exerted by a soap film on a plane cutting off an end of the surface. They also showed that it could be defined for more general constant mean curvature surfaces [25], where it serves a useful analytical tool.

A convenient way of expressing the force of a constant mean curvature surface of revolution which holds even when the generating curve is not a graph is given by the formula

f = hy^{2} − 2y cos α .  (3)

Here α denotes the angle between the positive x-axis and the generating curve of the surface of revolution, which is an oriented curve. This makes sense for all values of α. In our applications, we alway orient the generating curve of a Delaunay surface, and thus fix an unambiguous sign for the mean curvature.

Lemma 3.2 Given h, x_{0}, y_{0} > 0 and ẏ_{0}, there exists a unique Delaunay curve δ through the point (x_{0}, y_{0}) having slope ẏ_{0} and generating a Delaunay surface D having mean curvature h. The slope ẏ_{0} can be infinite. The Delaunay surface is a sphere S if and only if

h = h_{S} = \frac{2}{y_{0}\sqrt{1 + ẏ_{0}^{2}}} .  (4)

If h > h_{S}, then D is a nodoid and δ lies strictly below the circle σ generating S in a neighborhood of (x_{0}, y_{0}). An arc of δ which is decreasing as it leaves (x_{0}, y_{0}) remains beneath σ as long as it remains a graph over the x-axis.

If h < h_{S}, δ lies strictly above σ in a neighborhood of (x_{0}, y_{0}). If 0 < h < h_{S} then D is an unduloid or cylinder, if h < 0 then D is a nodoid, and if h = 0 then D is a catenoid.

Proof: Solving Equation 1 for ÿ gives a second order ODE satisfied by the generating curves of Delaunay surfaces at points where ẏ is finite:

ÿ = −h(1 + ẏ^{2})^{3/2} + \frac{1 + ẏ^{2}}{y} .  (5)

Given a constant h and a choice of initial conditions y and ẏ this equation has a unique solution from standard existence and uniqueness of solutions to ODEs. These solutions can also be obtained by rolling an appropriate ellipse, hyperbola or parabola [14]. The latter approach also shows that we can find a nodoid passing through y with any mean curvature h ≠ 0 in the case where ẏ is infinite.

The mean curvature of the unique sphere centered on the x-axis and having slope ẏ_{0} at (x_{0}, y_{0}) is given by h_{S}. Assume now that σ is centered at (0, 0) on the xy-plane.

If h > h_{S} then Equation 5 implies that δ lies strictly below σ in a neighborhood of (x_{0}, y_{0}). If D is an unduloid then δ is an undulary whose graph crosses σ again in at least two more points (x_{1}, y_{1}) and (x_{2}, y_{2}). We can assume, after performing a left/right reflection if necessary, that x_{0} ≥ 0 and x_{1} > x_{0}. Denote by γ the subcurve of δ starting at (x_{0}, y_{0}) and running to (x_{1}, y_{1}). Let γ_{t} = γ+(t, 0) be a horizontal translate of γ by a distance t and let T = sup{t : γ_{t} ∩ S ≠ ∅}. Then γ_{T} is tangent to S at a point P and lies above σ in a neighborhood of P. This violates Equation 5, so D cannot be an unduloid. D cannot be a cylinder or a catenoid, so it can only be a nodoid. If a subarc γ of δ which is decreasing as it leaves (x_{0}, y_{0}) intersects σ while δ remains a graph then an identical argument shows that we can horizontally translate γ until it is tangent to σ while lying above it, which is impossible by Equation 5. (In fact it will remain beneath σ for somewhat longer than this, but this will not be relevant for us.)

If 0 < h < h_{S} then Equation 5 implies that δ lies strictly above σ in a neighborhood of (x_{0}, y_{0}). If D is a nodoid, consider the subarc β of δ decreasing from (x_{0}, y_{0}) until it reaches a vertical tangency (x_{v}, y_{v}). The arc β is a graph which may or may not cross σ at an additional point. In either case we can translate β horizontally to the left until the last time at which it intersects σ. The final intersection must be at an interior point of β, so after translation β becomes tangent to σ at an interior point while lying underneath σ, again impossible by Equation 5. Thus the surface must be an unduloid or cylinder.

If h < 0 then the surface must be a nodoid. □

Nodoids and unduloids can be distinguished by the sign of the product fh.

Corollary 3.3 For a nodoid, fh > 0. For an unduloid and a cylinder fh < 0. For a sphere, a catenoid and a plane, fh = 0.

Proof: For a sphere of radius r, we can compute f at a maximum point where

f = hr^{2} − 2r.

Since h = 2/r, we have f = 0. We also compute fh for an unduloid or cylinder at a local maximum. By Lemma 3.2 we know that 0 < h < h_{S} = 2/r. Then

f = hr^{2} − 2r < (2/r)r^{2} − 2r = 0.

So h > 0 and f < 0 implying that hf < 0.

We compute fh for a nodoid at a vertical tangency, where ẏ = +∞ and f = hy^{2}. Then fh = h^{2}y^{2} > 0. □

The next proposition deduces some useful properties of Delaunay surfaces.

Proposition 3.4 Let D be a Delaunay surface with generating curve δ, mean curvature h, maximum y-value y_{M} and minimum y-value y_{m}. Then

1. If δ is a nodary or undulary, the arc-length of one period of δ is 2π/h.

2. If δ is a nodary the period P of δ satisfies

2(y_{M} + y_{m}) ≤ P ≤ π(y_{M} + y_{m}).

3. If D is a nodoid,

h = 2/(y_{M} − y_{m}).

If δ is an undulary

h = 2/(y_{M} + y_{m}).

4. If D is a nodary, δ has non-zero curvature.

Proof:

1. In [23] the equation for the generating curve of a Delaunay surface is given in terms of arclength, and it is shown that nodaries and undularies are periodic, with the arc length of one period given by 2π/h.

2. The period of an undulary is equal to the perimeter of the ellipse that is rolled along the x-axis to generate it [14]. The inequalities follow by comparing twice the length of the major axis, 2(y_{M} +y_{m}), with the perimeter of an enclosing circle whose diameter is the length of the major axis, π(y_{M} + y_{m}).

3. If δ is a nodary, we calculate the force at y_{M} and y_{m}. At y_{M} we get

f = hy_{M}^{2} − 2y_{M}.

At y_{m} we get

f = hy_{m}^{2} + 2y_{m},

where the sign of the second term changes since the orientation has reversed. Setting these equal we get

h = \frac{2}{y_{M} − y_{m}} .

For an unduloid, the proof is identical except for one sign change.

4. At a point where δ has zero curvature, k_{m} = 0 so

h = k_{p} = \frac{1}{y\sqrt{1 + ẏ^{2}}}

Equation 4 gives the mean curvature of a sphere with slope ẏ at height y as:

h_{S} = \frac{2}{y\sqrt{1 + ẏ^{2}}} > 0.

Where ẏ is defined, this is exactly twice the value of h calculated above. For a nodary, h > h_{S} or h < 0 by Lemma 3.2, so we cannot have h_{S} = 2h and there are no points of zero curvature on a nodary. At points where the nodary is vertical, h = k_{m} since k_{p} = 0. Thus the curvature is non-zero at vertical tangencies except possibly when h = 0, which would imply that the Delaunay surface is a plane or a catenoid, rather than a nodoid as assumed. □

Important stability formulas for constant mean curvature surfaces were developed by J.L. Barbosa and M. do Carmo [8]. Stability refers to the behavior of a surface when a compact subsurface is deformed by a variation, while holding its boundary fixed. A subsurface of a complete constant mean curvature surface is called stable if there is no compactly supported normal variation which decreases its area while preserving the volumes on each side, and unstable otherwise. Unstable constant mean curvature surfaces cannot form part of a minimizing bubble.

Proposition 3.5 The smooth subsurfaces of a minimizing bubble S(v_{1}, v_{2}) are stable subsurfaces of Delaunay surfaces.

Proof: S(v_{1}, v_{2}) is a piecewise-smooth surface of revolution whose smooth subsurfaces have constant mean curvature by Theorem 2.6 and Theorem 2.2. Any constant mean curvature surface of revolution in ℝ^{3} is a Delaunay surface. If one of these subsurfaces is unstable, then there is a variation which maintains the volume of each region while decreasing the area. The surface then cannot be part of a minimizing bubble. □

We next deduce a formula for the horizontal distance between two points on a Delaunay curve whose y-coordinates are known. We require that the second point be above the first, but it can be either to the left or to the right. We also get a formula for the volume under a Delaunay surface.

Proposition 3.6 Let (x_{1}, y_{1}) and (x_{2}, y_{2}) be two points on a Delaunay curve.

Suppose that y_{1} < y_{2} and that between the points \frac{dy}{dx} ≠ 0. Then

x_{2} − x_{1} = ∫_{y_{1}}^{y_{2}} \frac{t}{\sqrt{(2y + t)(2y − t)}} dy  (6)

where t = t(y) = hy^{2} − f.

If x_{1} < x_{2} and the curve also satisfies \frac{dx}{dy} > 0 between the two points then

the volume underneath the surface of revolution generated by the curve between x = x_{1} and x = x_{2} is given by:

V (x_{1}, x_{2}) = ∫_{y_{1}}^{y_{2}} \frac{πy^{2}t}{\sqrt{(2y + t)(2y − t)}} dy  (7)

If instead we have x_{2} ≤ x_{1} in the previous case then the formula gives the negative of the volume underneath the surface of revolution generated by the curve between x = x_{1} and x = x_{2}.

Proof: We solve Equation 2 for dx/dy, using that dy/dx ≠ 0, and then apply the change of variables formula.

f = hy^{2} − \frac{2y}{\sqrt{1 + ẏ^{2}}}
\sqrt{1 + ẏ^{2}} = \frac{2y}{t}
\frac{dx}{dy} = \frac{t}{\sqrt{(2y + t)(2y − t)}}

Note that the sign of t = t(y) = 2y cos α is positive if and only if the curve is oriented to the right, i.e. if cos α > 0 in Equation 3. Hence the integrand give the correctly signed measure of displacement in the x direction.

The volume of the region enclosed by a surface of revolution generated by a graph y(x) between x_{1}, x_{2} with x_{1} ≤ x_{2} is given by

v = ∫_{x_{1}}^{x_{2}} πy^{2} dx .

Since we have already deduced a formula for dx, the volume formulas follow immediately. □

Note that if y_{1} > y_{2} then the formulas in Proposition 3.6 would have the this sign reversed. A choice in sign was made in taking the positive square root, and this choice is correct if the curve is oriented so that y is increasing. However we need to assume that x_{2} > x_{1} only when calculating the volumes, and not in Equation 6.

The integrals of Proposition 3.6 are singular at a local minimum or maximum, where dy/dx = 0, and we need to apply a change of variables to obtain a formula which holds near such points. Note that the formulas below do not require that y_{2} > y_{1}, but do require that the Delaunay curve is a graph over the x-axis, with x_{2} > x_{1}.

Proposition 3.7 Let (x_{1}, y_{1}) and (x_{2}, y_{2}) be two points on a Delaunay curve which is a graph over the x-axis with x_{2} > x_{1} and suppose there is exactly one critical point between the two points, with y-value labeled ymin or ymax depending on whether it is a local minimum or maximum. Let h, f be the mean curvature and force associated to the Delaunay surface generated by the curve. Set t = t(y(z)) = h(y(z))^{2} − f.

If ymin is a local minimum, let y = y(z) = ymin + z^{2}, z_{1} = −\sqrt{y_{1} − ymin} and z_{3} = \sqrt{y_{2} − ymin}. Then

x_{2} − x_{1} = ∫_{z_{1}}^{z_{3}} \frac{2t}{\sqrt{(2y + t)(2 − hymin − hy)}} dz  (8)

and

V (x_{1}, x_{2}) = ∫_{z_{1}}^{z_{3}} \frac{2πy^{2}t}{\sqrt{(2y + t)(2 − hymin − hy)}} dz .  (9)

If ymax is a local maximum, let y(z) = ymax − z^{2}, z_{2} = \sqrt{ymax − y_{2}} and z_{4} = −\sqrt{ymax − y_{1}}. Then

x_{2} − x_{1} = ∫_{z_{4}}^{z_{2}} \frac{2t}{\sqrt{(2y + t)(hy − 2 + hymax)}} dz  (10)

and

V (x_{1}, x_{2}) = ∫_{z_{4}}^{z_{2}} \frac{2πy^{2}t}{\sqrt{(2y + t)(hy − 2 + hymax)}} dz .  (11)

Finally, if the curve is an increasing graph over the x-axis between x_{1} and x_{2} then Equation 8 continues to hold if the sign of z_{1} is changed.

Proof: These formulas follow from applying a change of variables to y in Equation 6 and Equation 7. Suppose ymin is a local minimum. Divide the Delaunay curve into the two pieces on each side of the minimum, and apply Proposition 3.6 and the substitution y = ymin + z^{2} to each piece.

One of the terms in the square root in the denominator in Equation 6 factors as

2y − t = 2y − hy^{2} + f = (y − ymin)(2 − hymin − hy) = z^{2}(2 − hymin − hy) since f = hymin^{2} − 2ymin. A resulting factor of z in the denominator, which causes the singularity, cancels with a z in the numerator from the change of variables. The integrands are even functions of z, and y = ymin corresponds to z = 0, so the curve portions on opposite sides of the minimum can be mapped to positive and negative z values, and combined into one integral.

The maximum case is similar. Near a maximum the substitution y = ymax− z^{2} is used. Several minus signs appear in the derivation, resulting in the given expression.

The final statement in the Proposition is a direct application of the change of variables formula to Equation 6. □

4 Torus bubbles

We established in Lemma 3.5 that a minimizing bubble must be obtained by revolving a union of Delaunay curves contained in the upper half plane. A key case of such a bubble is the torus bubble, constructed as follows. Take two circular arcs of the same radius, facing each other, each with one endpoint and center on the x-axis, and connect the other endpoints with two different Delaunay curves meeting at 120 degrees. Rotating around the x-axis, we get a piecewise-smooth surface surrounding two components, one homeomorphic to a torus, which we call the torus component T, and one homeomorphic to a ball, which we call the ball component B. It is not immediately clear whether it is possible to make such a construction so that the curves meet at 120^{∘} angles and the mean curvatures sum to zero around each triple curve. Such torus bubbles do indeed exist, and we will need to show that none of them are minimizers. Figure 7, due to John Sullivan, shows a torus bubble.

Figure 7: A torus bubble in ℝ^{3}.

The reason that torus bubbles play so central a role in our argument is due to another key result of Hutchings:

Theorem 4.1 A minimizing bubble enclosing two equal volumes must be either a double bubble or a torus bubble.

Proof: Theorem 2.10 states that S(v, v) contains two connected components. The two possible configurations are then a consequence of Hutchings’ structure theorem for minimizing bubbles [22, Theorem 5.1]. □

We label the surfaces of a minimizing torus bubble as indicated in Figure 8. The component homeomorphic to a solid torus is denoted T and the component homeomorphic to a ball is denoted B. The volumes of the inner ball component and outer torus component are denoted by v_{i} and v_{o} respectively. The inner Delaunay surface on the boundary of T is denoted T_{i} and the outer Delaunay surface is denoted T_{o}. The generating curve of T_{o} is denoted by τ_{o} and that of T_{i} by τ_{i}. The angles subtended by the generating curves for S_{1} and S_{2} are denoted by θ_{1} and θ_{2}. The mean curvatures of T_{i} and T_{o} are h_{i} and h_{o}, with signs chosen so that so that h_{i} is positive if the mean curvature vector points into B and negative if it points into T, and h_{o} is positive if the mean curvature vector points into T and negative if it points towards outside the bubble.

Denote by (x_{1}, y_{1}) the initial point where the two Delaunay surfaces start, and (x_{2}, y_{2}) the point where they rejoin. To specify a sign for the mean curvature, we orient the generating curves τ_{i} and τ_{o} of T_{i} and T_{o} so that they run from (x_{1}, y_{1}) to (x_{2}, y_{2}).

Figure 8: Parameters of a torus bubble.

Torus bubbles for various ranges of θ_{1} are indicated in Figures 18, 19 and 20. We will show that θ_{1} and h_{o} parameterize the space of torus bubbles and that no choice of θ_{1} and h_{o} gives a minimizing bubble. To do so, we first present some properties of minimizing torus bubbles, mainly derived from the work of Hutchings.

Lemma 4.2 In a minimizing torus bubble, S_{1} and S_{2} have the same mean curvature.

Proof: This is an immediate consequence of Corollary 2.4. □

Lemma 4.3 In a minimizing torus bubble, h_{o} ≥ 0 and h_{i} = 2 − h_{o}.

Proof: The three mean curvatures around a triple curve sum to zero by Lemma 2.3 implying that h_{o} = 2 − h_{i}. If h_{o} < 0 then there is a deformation which pushes T_{o} outward, increasing the volume of the torus component and decreasing its area. Corollary 2.8 then gives a contradiction. □

Proposition 4.4 The regions bounded by a minimizing torus bubble are not all preserved by reflection through some plane P perpendicular to the x-axis.

Proof: If there is such a symmetry, then the bubble is invariant by reflection through the xy-plane and by reflection through a perpendicular plane P. The argument given in Theorem 2.6 shows that the bubble is a surface of revolution around the line of intersection of these two planes. But it is also a surface of revolution around the x-axis, and therefore a surface of revolution around two perpendicular axis. It must then be a union of spheres, and not a torus bubble. □

Corollary 4.5 Given a minimizing equal-volume torus bubble, we can rescale and reflect so that θ_{1} < θ_{2} and each of the two circular arcs generating S_{1} and S_{2} has curvature one.

Proof: If θ_{1} = θ_{2}, Lemma 4.4 gives a contradiction. If θ_{1} > θ_{2}, a left/right reflection interchanges the two angles. Lemma 4.2 implies that the two circular arcs generating S_{1} and S_{2} have the same curvature, and thus the same radius. We can rescale and assume without loss of generality that each has radius one (and thus mean curvature two.) □

From now on we will restrict attention to torus bubbles with θ_{1} < θ_{2} and with spherical caps of radius one.

Proposition 4.6 If a torus bubble contains a surface T_{1} which is a subsurface of a nodoid and the nodary generating T_{1} contains two interior vertical tangencies, then T_{1} is unstable and not a part of a minimizing bubble.

Proof: We will construct a volume preserving Jacobi field supported on a proper subset of T_{1}, establishing instability as in [8, Proposition 2.24].

We start with the constant unit vector field v which points in the direction of the positive y-axis. The vector field v is tangent to T_{1} along two circles which are generated by rotating the interior vertical tangencies of the nodary. Let T′_{1} denote the compact portion of T_{1} between these two circles. Define the vector field W =< v, N > N. where N is the normal vector field along T′_{1}. and extend W to be the zero vector field elsewhere on T_{1}. The integral of the function f =< v, N > along T_{1} is zero, since the involution (x, y) → (−x, −y) of T_{1} takes f to −f. So the normal variation defined by f preserves volume in the sense of [8, Proposition 2.24]. Moreover [8, Proposition 2.24] shows that fN is a Jacobi Field, since it is constructed by taking the normal part of a constant vector field, and since the nodary extends beyond the two vertical tangencies, the corresponding Jacobi Field is supported on a proper subset of T_{1}. The Morse Index Theorem implies that the surface T_{1} is unstable and that there is a smaller area surface separating the same volumes. □

The next proposition establishes that minimizing torus bubbles do not contain Delaunay curves which are longer than a full period.

Proposition 4.7 If τ_{i} or τ_{o} contain a full period, then the torus bubble is not a minimizer.

Proof: Suppose to the contrary that a minimizing torus bubble contains a subsurface generated by a curve which contains a full period. If a boundary surface of the torus component is a subsurface of a nodoid with width more than a period, then the associated nodary contains two vertical tangencies, which implies instability by Lemma 4.6. Similarly, if a boundary surface of the torus component is a subsurface of a cylinder of radius r and its width is greater than 2πr, the period of a cylinder, than the surface is shown to be unstable in [8].

If a boundary surface of the torus component is a subsurface U of an unduloid then the lemma is less obvious. Note that we cannot necessarily construct a volume preserving Jacobi field from the vector field ∂/∂x, since the Delaunay curve may not contain two maxima or two minima. We establish the Lemma with a reflection argument. Assume first that the unduloid containing a full period is the inner surface T_{i}.

An arc u of an undulary which is longer than a period contains three distinct points with the same y value, (x_{1}, y_{1}), (x_{2}, y_{1}), (x_{3}, y_{1}), with x_{1} < x_{2} < x_{3}, and with the slopes at (x_{1}, y_{1}), (x_{3}, y_{1}) each being non-zero and equal to the negative of the slope at (x_{2}, y_{1}). Reflect the part of u between the planes x = x_{1} and x = x_{3} through the plane x = (x_{1} + x_{3})/2 to get a new surface of revolution U′ generated by u′, the reflection of u. Reflection preserves both the area and the volume under U between the planes x = x_{1} and x = x_{3}. Suppose first that u′ does not intersect τ_{o}. Then the volume of the torus bubble is also preserved. The reflection creates a surface which is not smooth along the planes x = x_{1} and x = x_{3} since reflection changes the sign of the slope of the reflected surface along these planes, and this slope was not zero. Thus we have found a solution to the minimization problem which contains a surface which contradicts Theorem 2.2.

Possibly u′ does intersect the outer curve τ_{o}. See Figure 9. In that case the total volume enclosed by the union of U′ and T_{o} between x = x_{1} and x = x_{3} is increased, or at least not decreased. We know that U′ does not intersect the spherical caps of the torus bubble, since they do not meet the region between the planes x = x_{1} and x = x_{3}. Since the volume under U′ is equal to the volume of the ball region of the torus bubble B, the volume in the remaining (disconnected) region must be larger than that in the original torus region T. There is no increase in area, contradicting Corollary 2.8.

Note that have used here the fact that two distinct Delaunay curves intersect locally at finitely many points.

Assume now that the unduloid is the outer surface T_{o}. We again reflect it to get a new boundary surface U′. The total volume underneath U′ is the same as that underneath U. Possibly T_{i} intersects the reflected surface U′. In this case the total volume underneath the union of the two surfaces has increased. Divide the region underneath U′ so that the region under T_{i} is allocated entirely to B, as before, and all other regions are allocated to T. Then the total volume of B is preserved, that of T is increased, and area is the unchanged, a contradiction as before. □

We note that a somewhat similar argument was used by Athanassenas [7] in studying the stability of minimizing capillary surfaces. See also [41]. Our

Figure 9: Reflection of a piece of T_{i} could lead to extra intersections.

situation has some additional complications due to the possible intersections resulting from the reflection.

Lemma 4.8 If a torus bubble has θ_{1} < 60^{∘} and θ_{2} < 60^{∘} then it is not minimizing. In particular, θ_{2} ≥ 60^{∘} in a minimizing torus bubble.

Proof: If both θ_{1} and θ_{2} are less than 60^{∘}, the outer bubble T_{o} is generated by a nodary that becomes vertical twice. Lemma 4.6 implies it is unstable. Since θ_{2} is the larger angle, θ_{2} ≥ 60^{∘}. □

Lemma 4.9 Suppose that T is a torus bubble with θ_{2} ≥ 120^{∘}. The the torus component is completely contained inside the smaller region of the double bubble whose larger outer sphere subtends an angle equal to θ_{2}. Similarly the ball component contains in its interior the larger region of this double bubble.

Proof: Consider a double bubble which is symmetric around the x-axis and whose larger spherical cap S_{2} is on the right and subtends an angle of θ_{2} with the x-axis. Denote the smaller (leftmost) spherical cap by S_{o} and the spherical cap forming the interface by S_{i}. Then S_{o} is tangent to τ_{o} at (x_{1}, y_{1}). We will show that a torus region must be contained inside the smaller component B_{1} of this double bubble, the component bounded by S_{2} and S_{i}. Otherwise, we will show that τ_{i} and τ_{o} cannot rejoin after leaving (x_{1}, y_{1}), and thus cannot generate a torus region. See Figure 10 and Figure 8.

Given a vector in the upper half-plane, there is a unique circle perpendicular to the x-axis which is tangent to that vector. A Delaunay curve tangent to the same vector and generating a Delaunay surface with greater mean curvature than the corresponding sphere is a nodary, by Lemma 3.2. A Delaunay curve tangent to the same vector and generating a Delaunay surface with smaller, but still positive, mean curvature is an undulary or horizontal line.

If the mean curvature of T_{o} is equal to that of S_{o}, then the mean curvature of T_{i} is equal to that of S_{i} by Lemma 4.3. It follows from Lemma 3.2 that T_{o}, T_{i} and S_{o}, S_{i} coincide, and that there is no torus component.

Suppose now that the mean curvature of T_{o} is smaller than that of S_{o}. Then the mean curvature of T_{i} is smaller than that of S_{i} by Lemma 4.3. Since h_{o} ≥ 0, also by Lemma 4.3, it follows that T_{o} cannot be a nodoid, and therefore τ_{o} is a graph over the x-axis. In particular, τ_{o} does not recross the spherical cap S_{o}. The generating curve τ_{i} for T_{i} leaves the point (x_{1}, y_{1}) on the positive (right) side of S_{i}, and lies above S_{i}. Since S_{2} forms part of the torus bubble, τ_{i} must remeet τ_{o} without crossing S_{2} or itself if it is to generate a torus component with τ_{o}. If it is a graph over the x-axis then it cannot remeet τ_{o} at all. If it is not a graph, then it is a nodary which generates a nodoid with negative mean curvature. Such a curve crosses S_{2} before it can remeet τ_{o}.

Finally, suppose that the mean curvature of T_{o} is greater than that of S_{o}. Then Lemma 4.3 implies that the mean curvature of T_{i} is greater than that of S_{i}, and Lemma 3.2 implies that both T_{o} and T_{i} are subsurfaces of nodoids, and that near (x_{1}, y_{1}) T_{o} lies below S_{o} and T_{i} lies below S_{i}. See Figure 11. The nodary τ_{o} generating T_{o} either crosses itself before leaving B_{1} or leaves B_{1} through S_{i}. Similarly, the nodary τ_{i} generating T_{i} either crosses itself before leaving B_{1} or leaves B_{1} through S_{o}. The curves τ_{o} and τ_{i} cannot self-intersect before remeeting, since they generate a torus component. Therefore either τ_{o} and τ_{i} are contained in B_{1} or they remeet for the first time after τ_{o} leaves B_{1} through S_{i} and τ_{i} leaves B_{1} through S_{o}. The latter case is impossible, since τ_{i} must cross τ_{o} before it can leave B_{1} through S_{i}. It follows that the torus component is contained inside B_{1}.

The spherical caps S_{i} in the double bubble and S_{1} in the torus bubble are disjoint, since S_{1} has boundary circle lying to the left of S_{i}, and each is convex in opposite directions. So no part of the torus bubble intersects the interior of the larger, rightmost component of the corresponding double bubble, which must then be contained in the ball component of the torus bubble. This proves the last assertion of the Lemma. □

Figure 10: If h_{o} is too small, a torus component cannot exist.

Corollary 4.10 In a minimizing equal-volume torus bubble, θ_{2} < 120^{∘}.

Proof: If θ_{2} > 120^{∘} then Lemma 4.9 implies that the volume of the torus component is strictly less than the volume of the smaller component of the double bubble whose spherical caps subtend angles θ_{2} and (240−θ_{2})^{∘}. The other Figure 11: The torus component is trapped inside a component of an associated double bubble if θ_{2} > 120^{∘}.

component of the torus bubble contains in its interior the larger component of this double bubble, which in turn contains more than half the total volume of the double bubble. It follows that the torus component has strictly less volume than the ball component of the torus bubble, contradicting the assumption of equal volume. □

Recall that θ_{2} is the larger of the two angles subtended by the spherical caps, so the Corollary implies that both of θ_{1} and θ_{2} are less than 120^{∘}.

Lemma 4.11 In a minimizing equal-volume torus bubble, h_{o} > 0.

Proof: We know from Lemma 4.3 that h_{o} ≥ 0. If h_{o} = 0 then τ_{o} is a catenary. Since y_{2} ≥ y_{1} the catenary is given by a strictly increasing graph for {x ≥ x_{2}}, implying in particular that θ_{2} > 120^{∘} which violates Lemma 4.9. □

Knowledge of h_{o}, θ_{1}, θ_{2} can be used to give a fairly accurate qualitative picture of the torus bubble.

Proposition 4.12 A minimizing equal-volume torus bubble has the following properties:

1. τ_{i} has negative slope at (x_{1}, y_{1}) if and only if θ_{1} > 30^{∘}. Its angle with the positive x-axis is 30^{∘} − θ_{1}.

2. τ_{i} has positive slope at (x_{2}, y_{2}) if and only if θ_{2} > 30^{∘}. Its angle with the positive x-axis is θ_{2} − 30^{∘}.

3. τ_{o} has positive slope at (x_{1}, y_{1}) if and only if θ_{1} > 60^{∘}. Its angle with the positive x-axis is 150^{∘} − θ_{1}.

4. τ_{o} has angle at (x_{2}, y_{2}) with the positive x-axis equal to θ_{2} − 150^{∘}.

5. τ_{i} is a graph. It has a unique local minimum if and only if θ_{1} ≥ 30^{∘}. It never has a local maximum. If T_{i} is a nodoid then h_{i} < 0.

6. τ_{o} always has a unique local maximum. It has a vertical tangent on the left if and only if θ_{1} ≤ 60^{∘}. It never has a vertical tangent on the right. It never has a local minimum.

7. If θ_{1} ≥ 30^{∘} the local minimum ymin of τ_{i} has value

ymin = \frac{−f_{i}}{1 + \sqrt{1 + f_{i}h_{i}}}

where f_{i} is the force associated to T_{i}.

8. The local maximum ymax of τ_{o} has value

ymax = \frac{1 + \sqrt{1 + f_{o}h_{o}}}{h_{o}}

where f_{o} is the force associated to T_{o}.

Proof: The first four assertions follow from the fact that the surfaces forming a torus bubble meet at 120^{∘} angles along a triple curve.

If τ_{i} is not a graph then θ_{2} > 120^{∘}, violating Corollary 4.10. If θ_{1} ≥ 30^{∘} then τ_{i} has non-positive slope at (x_{1}, y_{1}). Since θ_{2} > θ_{1} the curve must have a local minimum in this case.

Suppose τ_{i} has a local maximum. Since θ_{2} > 60^{∘}, τ_{i} has positive slope at (x_{2}, y_{2}) and must have a local minimum as well. However θ_{2} > θ_{1} so τ_{i} must be longer than a period, contradicting Proposition 4.7. So it has no local maximum.

If h_{i} ≥ 0 and T_{i} is a nodoid, then we can compare T_{i} and T_{o} to a double bubble whose spherical interface S_{i} coincides with S_{1}, and apply Lemma 3.2 as in Lemma 4.9. T_{i} is trapped inside the component tB_{2} to the right of the interface and T_{o} is trapped inside the component B_{1} to the left of the interface, as in Figure 12. In particular, h_{o} < 0 in this case, a contradiction. So if T_{i} is a nodoid then h_{i} < 0.

Figure 12: A double bubble traps T_{i} and T_{o} if T_{i} is a nodoid with h_{i} > 0.

If θ_{1} < 30^{∘} then τ_{i} has positive slope at (x_{1}, y_{1}). If it has a minimum then it must pass through a maximum first, so it has no minimum.

τ_{o} is always oriented upwards at (x_{1}, y_{1}) and downwards at (x_{2}, y_{2}), so it always has a local maximum. If it has a local minimum then it would contain two local maxima, violating Proposition 4.7. If θ_{1} < 60^{∘} then the curve starts out to the left, therefore it has a vertical tangent at which it changes direction from left to right. Lemma 4.6 implies that there cannot have a second vertical tangent, so it has no vertical tangent on the right, as claimed. If θ_{1} ≥ 60^{∘} then τ_{o} initially is going right. If it has a vertical tangent on the right then the angle θ_{2} must be less than 60^{∘}, contradicting the assumption that θ_{1} ≤ θ_{2}. This gives the conclusion on when τ_{o} has vertical tangents.

To calculate the local minimum ymin of τ_{i}, we use Equation 2 of Section 3 for the force of T_{i},

f_{i} = h_{i}y^{2} − 2y cos α .

At a minimum, α = 0 and cos α = 1, so

f_{i} = h_{i}y^{2} − 2y.

If h ≠ 0 this quadratic expression gives two roots for y:

y = \frac{1 ± \sqrt{1 + h_{i}f_{i}}}{h_{i}} .

The value of the minimum is given by

\frac{1 − \sqrt{1 + h_{i}f_{i}}}{h_{i}}  (12)

whatever the sign of h_{i}. If h_{i}f_{i} > 0, as Corollary 3.3 implies happens at a local minimum for a nodoid, then the other root is negative and meaningless. If h_{i}f_{i} < 0, as Corollary 3.3 implies happens at a local minimum for an unduloid, then the other root corresponds to the local maximum of the unduloid. A double root occurs only in the case of a cylinder. Equation 12 becomes linear when h_{i} = 0. An equivalent expression for the minimum, which holds also when h_{i} = 0, is

ymin = \frac{−f_{i}}{1 + \sqrt{1 + f_{i}h_{i}}} .  (13)

For a local maximum, a similar analysis gives that ymax is one of the two roots

ymax = \frac{1 ± \sqrt{1 + f_{o}h_{o}}}{h_{o}} .  (14)

Since by Corollary 4.19 T_{o} is always a nodoid, h_{o} > 0 and f_{o} > 0, there is a unique positive root, and

ymax = \frac{1 + \sqrt{1 + f_{o}h_{o}}}{h_{o}} .

□

Proposition 4.13 Given θ_{1} and h_{o}, there is at most one corresponding minimizing equal-volume torus bubble.

Proof: Given θ_{1} and h_{o}, h_{i} = 2 − h_{i} is also determined. Thus S_{1}, T_{o} and T_{i} are uniquely determined. The spherical cap S_{2} is uniquely determined by the second intersection point (x_{2}, y_{2}) of τ_{o} and τ_{i}, assuming that this second point exists. If τ_{o} and τ_{i} do not meet in a second point then there is no torus bubble corresponding to θ_{1} and h_{o}. In general, τ_{o} and τ_{i} may intersect in many points, and it is necessary to show that only one of these can possibly lead to a torus bubble.

If τ_{o} and τ_{i} are graphs, then it is clear that their first intersection point with x coordinate larger than x_{1} is the unique intersection point which can give a torus bubble. In general, τ_{i} but not τ_{o} is a graph. However Proposition 4.12 shows that τ_{o} has no vertical tangents for x > x_{1}, so that both curves are graphs in this region and the same argument applies. □

Figure 13: Some eliminated 0 for cross-1 of torus bubbles.

Proposition 4.14 If a torus bubble has θ_{1} + θ_{2} > 180^{∘}, then it is not a minimizer.

Proof: Under this assumption there is a value of y realized by three distinct points on the circular arcs generating the spherical caps. Since θ_{2} ≥ θ_{1} we have in particular that θ_{2} > 90^{∘}. We know from Lemma 4.9 that θ_{2} < 120^{∘}. So 90^{∘} < θ_{2} < 120^{∘}, 60^{∘} < θ_{1} ≤ 120^{∘}. Then Lemma 4.12 implies that the curves τ_{i} and τ_{o} are both graphs. If θ_{1} + θ_{2} > 180^{∘} then there are three points on S_{1} ∪ S_{2} with y-value equal to y_{2}, (a_{1}, y_{2}), (a_{2}, y_{2}), (a_{3}, y_{2}), with a_{1} < a_{2} < a_{3} and a_{2} = x_{2}. Reflect the part of the bubble between the planes x = a_{1} and x = a_{2} through the plane x = (x_{1} +x_{2})/2 to get a new piecewise-smooth surface with the same area enclosing the same volumes, as in the proof of Lemma 4.7, but without needing to worry here about introducing self-intersections. The new piecewise-smooth surface has three surfaces meeting at an angle not equal to 120^{∘} and so is not minimizing, a contradiction. □

Lemma 4.15 For any constant b, the function f(c) = b(1 − c^{2}) + c\sqrt{3(1 − c^{2})} is at most 2-to-1 on (−1, 1).

Proof: f′(c) = −2bc − \sqrt{3}g(c) where

g(c) = (2c^{2} − 1)/\sqrt{1 − c^{2}}.

Since g′′(c) = 3(1 − c^{2})^{−5/2} > 0, g is convex on (−1, 1).

A line intersects a convex curve at most twice, so f′ has at most two roots on (−1, 1). Thus f has at most one local maximum and one local minimum on (−1, 1). Combined with f(0) = f(1) = 0, this completes the proof. □

We calculate in the next proposition the value of the forces f_{o} of T_{o} and f_{i} of T_{i}.

Proposition 4.16 For the torus bubble determined by θ_{1} and h_{o}, the forces f_{o} of T_{o} and f_{i} of T_{i} are given by:

f_{o} = (h_{o} − 1)(1 − c_{1}^{2}) + c_{1}\sqrt{3(1 − c_{1}^{2})}.

and

f_{i} = (h_{i} − 1)(1 − c_{1}^{2}) − c_{1}\sqrt{3(1 − c_{1}^{2})} = −f_{o} ,

where c_{1} = cos θ_{1}.

Proof: We first calculate the force f_{o} of T_{o}. At (x_{1}, y_{1}), y_{1} = sin θ_{1} and α = 5π/6 − θ_{1} is the angle between the positive x-axis and the graph of τ_{o},

f_{o} = h_{o}y_{1}^{2} − 2y_{1} cos α
= h_{o} sin^{2} θ_{1} − 2 sin θ_{1} cos(5π/6 − θ_{1})
= h_{o}(1 − cos^{2} θ_{1}) − 2 sin θ_{1}((−\sqrt{3}/2) cos θ_{1} − (1/2) sin (−θ_{1}))
= h_{o}(1 − c_{1}^{2}) + c_{1}\sqrt{3(1 − c_{1}^{2})} − (1 − c_{1}^{2})
= (h_{o} − 1)(1 − c_{1}^{2}) + c_{1}\sqrt{3(1 − c_{1}^{2})} .

Similarly,

f_{i} = h_{i}y_{1}^{2} − 2y_{1} cos α
= (h_{i} − 1)(1 − c_{1}^{2}) − c_{1}\sqrt{3(1 − c_{1}^{2})} .

Since h_{o} − 1 = 1 − h_{i} = −(h_{i} − 1) by Lemma 4.3, f_{o} = −f_{i}. □

Note: The fact that f_{o} = −f_{i} is an example of a general “balancing principle” for constant mean curvature surfaces, due to R. Kusner [25].

Lemma 4.17 For the torus bubble determined by θ_{1} and h_{o}, the angle θ_{2} is one of the (at most) two solutions of the equation

f(cos θ_{1}) = f(cos θ_{2}),

where

f(c) = (h_{o} − 1)(1 − c^{2}) + c\sqrt{3(1 − c^{2})}.

Proof: We calculate the force f_{o} of T_{o} at both (x_{1}, y_{1}) and (x_{2}, y_{2}).

At (x_{1}, y_{1})

f_{o} = (h_{o} − 1)(1 − c_{1}^{2}) + c_{1}\sqrt{3(1 − c_{1}^{2})},

where y_{1} = sin θ_{1} and c_{1} = cos θ_{1}. A similar calculation at θ_{2} shows

f_{o} = (h_{o} − 1)(1 − c_{2}^{2}) + c_{1}\sqrt{3(1 − c_{2}^{2})}.

Setting the two forces equal and applying Lemma 4.15 gives that there are at most two possible values for c_{2} = cos θ_{2}. Since 0 < θ_{2} < 180^{∘}, where the cosine function is monotonically decreasing, this implies that there are at most two possible values for θ_{2}. □

One solution occurs when the angles are equal. An example of such a symmetric torus bubble, depicted in Figure 7, occurs with θ_{1} = θ_{2} = 90^{∘}, and h_{o} approximately 1.9848. However these symmetric solutions never lead to a minimizing torus bubble by Corollary 4.5.

Recall that we remarked after Corollary 4.5 that we can without loss of generality restrict attention to non-symmetric torus bubbles. This means that we are working in the case where θ_{1} < θ_{2}. For each θ_{1}, h_{o}, there is at most one value of θ_{2} which can potentially lead to a stable torus bubble.

Proposition 4.18 If a minimizing torus bubble has 60^{∘} ≤ θ_{1} ≤ 120^{∘} then h_{o} > 1 − \sqrt{3} cot θ_{1} and h_{i} < 1 + \sqrt{3} cot θ_{1}. The same statement holds with θ_{2} replacing θ_{1}.

Proof: We again apply Lemma 3.2 to compare the torus bubble with a double bubble whose smaller spherical cap S_{o} subtends an angle of θ_{1}, as in the proof 1 Lemma 4.9. See Figure 14. If h_{i} is larger than the mean curvature of S_{i}, then T_{i} is a nodoid contained inside the smaller component of the double bubble, as shown in Lemma 3.2. See Figure 14. It follows that τ_{o} cannot remeet τ_{i} and the two curves cannot generate a torus component. □

Figure 14: The mean curvature of T_{o} cannot be too small if τ_{o} and τ_{i} are to meet again.

Corollary 4.19 T_{o} is a nodoid.

Proof: If 60^{∘} ≤ θ_{1} ≤ 120^{∘} then Proposition 4.18 shows that the mean curvature of T_{o} is greater than that of the sphere which is tangent to it at (x_{1}, y_{1}), so Lemma 3.2 implies that T_{o} is a nodoid. If θ_{1} < 60^{∘} then τ_{o} starts out heading to the left. It must become vertical, hence T_{o} must be a nodoid. □

Lemma 4.20 Let N be a subcurve of a nodary generating a nodoid of mean curvature h > 0, running between a vertical tangency on the left at (x_{v}, y_{v}) and a vertical tangency on the right at (x′_{v}, y′_{v}), with a unique local maximum at (x_{M}, y_{M}). Then a circle C of curvature h which is tangent to N at (x_{v}, y_{v}) lies underneath N. Moreover x_{v} − x_{m} is greater than the radius 1/h of the circle.

Proof: In a deleted neighborhood of (x_{v}, y_{v}), the curvature of N is less than that of C, so that C locally lies beneath N. If C crosses N between x_{v} and x′_{v}, let C′ be the subcurve of C starting at (x_{v}, y_{v}) and running to the first point of intersection of C with N. Let C′_{t} = C′ + (0, t) be a vertical translate of C′ by a distance of t along the y-axis. Let T = sup{t : C′_{t} ∩ N ≠ ∅}. Thus C′_{T} is tangent to N at a point P and lies above N near P. The curvature of N at P is given by k_{m} = h − k_{p} < h. Since the curvature of C′_{T} is greater than that of N, it cannot lie above it, and it follows that x_{M} − x_{v} is larger than the radius of C, x_{M} − x_{v} > 1/h. □

The next proposition rules out the existence of minimizing torus bubbles with h_{o} and θ_{1} both close to 0. The numerical calculations that we will apply are badly behaved in this region, so we use this geometric argument to exclude it.

Proposition 4.21 If a torus bubble has h_{o} ≤ 0.2 and θ_{1} ≤ 5.7^{∘}, then it is not a minimizer.

Proof: We will show that before τ_{i} and τ_{o} rejoin at (x_{2}, y_{2}), either τ_{o} has two vertical tangencies, or τ_{i} traverses more than a full period. Note that our assumption implies that h_{i} ≥ 1.8

Proposition 3.4 implies that τ_{i} is a graph and its period is no longer than 2π/h_{i} < 3.5. The bubble is not a minimizer if x_{2} − x_{1} > 3.5 by Proposition 4.7. The initial angle of τ_{i} from the positive x-axis is 30^{∘} − θ_{1}, so it cannot turn (anti-clockwise) through an angle of more than 60^{∘} +θ_{1}. When they meet again τ_{i} and τ_{o} intersect at an angle of 120^{∘}. It follows that τ_{o} does not intersect τ_{i} before it turns (clockwise) through an angle greater than 180^{∘} − θ_{1}. Therefore the two curves cannot intersect before the nodary τ_{o} reaches its maximum, which happens when τ_{o} has turned through an angle of 150^{∘} − θ_{1}. Let (x_{M}, y_{M}) be the point where τ_{o} reaches its maximum and (x_{v}, y_{v}) the point where τ_{o} first goes vertical. Then x_{2} − x_{1} > x_{M} − x_{1}.

Now h_{o} = k_{p} + k_{m} where k_{m} is the curvature of the generating curve in the xy plane and k_{p} is the principal curvature due to rotation around the axis. At (x_{v}, y_{v}), k_{p} = 0 and h_{o} = k_{m} ≤ .2 by hypothesis.

Lemma 4.20 implies that x_{M} − x_{v} > 1/h_{o} ≥ 5.

The assumptions h_{o} ≤ .2 and θ_{1} ≤ 5.7^{∘} give an upper bound for the force f_{o} of T_{o}. Calculating at height y_{1} = sin θ_{1} < θ_{1} ≤ 5.7π/180 < .1 gives

|f_{o}| = |y_{1}^{2}h_{o} − y_{1} cos (30^{∘} + θ_{1})|
< |y_{1}^{2}h_{o} − y_{1}\sqrt{3}/2|
< (.1)^{2}h_{o} + .1\sqrt{3}/2
< .09.

We bound the value of y_{v} from above by solving f_{o} = y_{v}^{2}h_{o}.

y_{v}^{2} = f_{o}/h_{o} < .09/h_{o}

so

y_{v} < \sqrt{.09/h_{o}}.

We next bound x_{1} − x_{v} from above by noting that τ_{o} is convex between x_{v} and x_{1} by Lemma 3.4. This implies that (x_{v}, y_{v}) lies to the right of the point (x, y_{v}) at which the line of slope τ̇_{o}(x_{1}) has height y_{v}, so that

|x_{v} − x_{1}| < |x − x_{1}| = |(y_{v} − y_{1})/τ̇_{o}(x_{1})| < \sqrt{3}(y_{v} − y_{1})
< \sqrt{3}y_{v} < \sqrt{3}\sqrt{.09/h_{o}} = \sqrt{.27/h_{o}}.

This implies that

x_{M} − x_{1} = (x_{M} − x_{v}) − (x_{1} − x_{v}) > 1/h_{o} − \sqrt{.27/h_{o}}.

For 0 < h_{o} ≤ .2, differentiation shows that the function 1/h_{o} − \sqrt{.27/h_{o}} is decreasing and its minimum value occurs at h_{o} = .2, so that x_{M} − x_{1} > 1/(0.2) − \sqrt{.27/(0.2)} > 3.8. Thus τ_{o} travels a distance of at least 3.8 to the right before hitting τ_{i}. But the period of τ_{i} is less than 3.5. We have shown that either τ_{o} changes direction and has two vertical tangencies or τ_{i} contains a full period. In the first case it is not a minimizer by Lemma 4.8, and in the second case it is not a minimizer by Proposition 4.7. □

The next lemma will be used to give bounds on the size of a circle enclosing a loop of a nodary.

Lemma 4.22 If an embedded curve γ in the positive quadrant of the plane is tangent to the y-axis at the origin and has curvature k > 1/r > 0, then it does not meet C(r), the circle of radius r centered at (0, r), at any point other than the origin.

Proof: If it does cross C(r), denote by γ′ the subcurve of γ starting at the origin and running to the first point of intersection with C(r), as in Figure 15. The subcurve γ′ is a graph over the x-axis since it is embedded. Let γ′_{t} = γ′+(0, t) be a vertical translate of γ′ by a vertical distance of t. Let T = sup{t : γ′_{t} ∩ C(r) ≠ ∅}. Then γ′_{T} meets C(r) at a point P in its interior, but not the interior of C(r), and thus is tangent to C(r) at P. Since the curvature of γ′_{T} is greater than that of C(r), this is a contradiction. □

Figure 15: The curvature of the curve γ cannot be greater than 1/r.

Lemma 4.23 If τ_{o} has curvature greater than k, and τ_{i} has curvature less than −k, then their union is contained in the intersection of two circles of radius r = 1/k meeting at 120^{∘}.

Proof: Construct circles C_{1} and C_{2} of radius 1/k, meeting at 120^{∘} at (x_{1}, y_{1}), with C_{1} tangent to τ_{o} and C_{2} tangent to τ_{i} as in Figure 16. Let D denote the line segment joining the two intersection points of C_{1} and C_{2}. By Lemma 4.22 we know that each of τ_{o} and τ_{i} cannot leave the intersection of these two circles without first crossing D. Since each is convex in C_{1} ∩ C_{2} it follows that they remeet before either one leaves C_{1} ∩ C_{2}. It follows that the intersection of the torus component with the upper half-plane is contained in C_{1} ∩ C_{2}. □

Figure 16: The torus component generating curves are trapped inside two circles meeting at 120^{∘}.

Lemma 4.24 If a torus bubble has h_{i} ≤ −k < 0 and the intersection of the torus component with the upper half-plane has minimum y-value equal to y_{0} with y_{0} > 1/k,then the torus component contains volume

v ≤ 2.5π(y_{0} + \frac{\sqrt{3}}{2} \frac{y_{0}}{ky_{0} − 1})(\frac{y_{0}}{ky_{0} − 1})^{2} .

Proof: In a torus bubble with h_{i} ≤ −k < 0 we have that |h_{o}| = −h_{i} + 2 ≥ k. For a Delaunay curve passing through a point (x, y) in the upper half-plane,

|k_{m}| ≥ |h| − |k_{p}|
= |h| − \frac{1}{y\sqrt{1 + ẏ^{2}}}
≥ |h| − \frac{1}{y} .

For the curves τ_{i} and τ_{o} we have y ≥ y_{0}, so that

|k_{m}| ≥ k − \frac{1}{y_{0}} .  (15)

So the curves τ_{o} and τ_{i} each have curvature with absolute value greater than

\frac{ky_{0} − 1}{y_{0}}.

By Lemma 4.23 the intersection of the torus component with the upper halfplane is contained in the intersection w of two circles of radius r = y_{0}/(ky_{0} − 1) meeting at an angle of 120^{∘}. The area of w is

a = (2π/3 − \sqrt{3}/2)r^{2} < 1.25(\frac{y_{0}}{ky_{0} − 1})^{2}.

The diameter of w is

d = \sqrt{3}\frac{y_{0}}{ky_{0} − 1}.

Since w is centrally symmetric and has minimal y-value smaller or equal to y_{0}, the center of mass of w has y-coordinate at most y_{0} + d/2. Pappus’ Theorem for the volume of a solid of revolution obtained by rotating this area around the x-axis gives that

v ≤ 2π(y_{0} + \frac{d}{2})1.25(\frac{y_{0}}{ky_{0} − 1})^{2},
≤ 2.5π(y_{0} + \frac{\sqrt{3}}{2} \frac{y_{0}}{ky_{0} − 1})(\frac{y_{0}}{ky_{0} − 1})^{2}.

□

We apply these volume estimates to reduce the range of possible mean curvatures in a minimizing equal-volume torus bubble.

Proposition 4.25 In a minimizing equal-volume torus bubble, h_{o} ≤ 10.

Proof: The hypothesis h_{o} ≤ 10 is equivalent to h_{i} ≥ −8. Suppose to the contrary that h_{i} < −8 in a minimizing equal-volume torus bubble.

Proposition 4.12 states that τ_{i} is a graph, so that the volume v_{i} of the ball component B is strictly larger than the volume of the region under S_{2}. Lemma 4.9 and Lemma 4.8 imply that 60^{∘} < θ_{2} < 120^{∘}. It follows that v_{i} > 5π/24, the volume under a unit radius spherical cap subtending an angle of 60^{∘}. We now estimate the volume v_{o} of the torus component. The curves τ_{i} and τ_{o} component contain the point (x_{2}, y_{2}) where \sqrt{3}/2 ≤ y_{2} ≤ 1. Their maximum and minimum y-values differ by less than 2/|h_{i}| ≤ 0.25 by Proposition 3.4, and therefore both curves are contained in \{y ≥ \sqrt{3}/2 − 0.25\}. The expression \frac{y_{0}}{ky_{0} − 1} is decreasing with respect to y_{0} when y_{0} > 1/k, and decreasing with

respect to k, and also 0.6 < \sqrt{3}/2 − 0.25. We apply Lemma 4.24 with k ≥ 8 and y_{0} ≥ 0.6, giving that

v_{o} < 2.5π(y_{0} + \frac{\sqrt{3}}{2} \frac{y_{0}}{ky_{0} − 1})(\frac{y_{0}}{ky_{0} − 1})^{2}
< 2.5π(1 + \frac{\sqrt{3}}{2} \frac{0.6}{4.8 − 1})(\frac{0.6}{4.8 − 1})^{2}
< 0.08π .

Thus v_{o} < 0.08π < 5π/24 < v_{i}, contradicting the assumption that the torus bubble encloses equal volumes. □

Our algorithm will need the following rather technical result about the slope of a generating curve.

Lemma 4.26 Suppose that τ_{i} satisfies h_{i} ≤ 2 and 0 ≤ θ_{1} ≤ 4^{∘}. Let α be the angle made by τ_{i} with the positive x-axis and let α_{1} be the value of this angle at (x_{1}, y_{1}). Then α ≥ α_{1} for as long as τ_{i} remains below height {y = 1/2}. In particular, α ≥ 26^{∘} below this height.

Proof: First note that if h_{i} ≤ 0 then τ_{i} is a nodary or catenary, and the angle made by τ_{i} with the positive x-axis is increasing as one moves along τ_{i}, which implies the conclusion of the Lemma. So we can assume that 0 < h_{i} ≤ 2.

Let y be any point on τ_{i}. Then

h_{i}y^{2} − 2y cos α − f_{i} = 0 .

A given value of cos α is realized by at most two y values solving this quadratic equation. Suppose that as the τ_{i} goes from a minimum to a maximum, the angle 26^{∘} ≤ α_{1} = 30^{∘} −θ_{1} ≤ 30^{∘} occurs twice. The two y-values which share the same α_{1} both solve the above equation. Thus they sum to 2 cos α_{1}/h_{i} ≥ 2(.8)/2 ≥ .8. Since the smaller root has y-coordinate sin(4^{∘}) < .1 the larger root is greater than .7 > 1/2. The angle α is increasing near the smaller root. Thus τ_{i} subtends an angle of at least α_{1} with the x-axis as long as τ_{i} remains below height 1/2. Since α_{1} = 30^{∘} − θ_{1} ≥ 26^{∘}, the final conclusion holds. □

5 Computation

We have shown that to solve Conjecture 3 in the equal volume case it suffices to show that certain torus bubbles are less efficient at enclosing two equal volumes then a double bubble. In this section we will show how these torus bubbles are ruled out by a computation. The algorithm for this computation is given below.

The basic idea is to consider a domain of torus bubbles corresponding to a product of small intervals in each of θ_{1} and h_{o}, and to calculate as much as possible about the geometry of the corresponding torus bubbles. Various calculations are then applied which rule these torus bubbles out as potential minimizers. The accuracy of these calculations depends on the size of the domain rectangle we start with. The computational scheme will succeed if these can be chosen small enough to get sufficient accuracy, yet large enough that a reasonable number of them cover all the possibilities.

The computation was performed using double precision floating point numbers. The fundamental data type used is the IEEE 754 32-bit real number, see [5]. This is a binary representation with a 23-bit mantissa (plus an implied leading bit), an 8-bit binary exponent, and a sign bit. The (finitely many) real numbers that can be represented in this standard with no error are called representable. The IEEE standard specifies that the add, subtract, multiply, divide, and square root operations be performed as if done exactly and rounded to the nearby representable number according to the rounding mode in effect. There are three rounding modes that can be chosen: up, down, or nearest. The values +∞ and −∞ are representable and behave in a specified way. Floating point exceptions are masked, but flags are sticky and available for clearing and inspection. Most computers in use today implement the IEEE standard. We denote the result of rounding a real number x down to a representable number smaller or equal to x by \underline{x}, and similarly denote by \overline{x} the result of rounding x up to a representable number greater or equal to x.

Combined with the methods of interval arithmetic, see Moore [29] and Alefeld-Herzberger [2], the IEEE standard allows numerical calculation with exact bounds on accuracy. Interval arithmetic is a method by which a real valued function on the reals can be extended to an interval valued function of intervals. An interval is formed from two representable reals. Mathematically, it represents the closed interval between the two reals. The add, subtract, multiply, divide, and square root functions are extended to intervals by the IEEE operations on reals along with directed rounding. We only use the round-to-nearest mode in calculating averages, (see the procedure avgwt below), as these are only used to divide intervals and do not need to have bounds calculated.

Generally, we say that an interval X is an extension of a real number x if x ∈ X and an interval valued function of intervals F is an extension of f if x ∈ X ⇒ f(x) ∈ F (X).

In this section, real numbers are denoted by lower case, and intervals with representable real endpoints by upper case. The reals embed into the set of intervals by mapping a real to the smallest interval containing it. The lower and upper bounds to an interval are denoted with lower and upper bars, so X = [\underline{X}, \overline{X}]. Note that when X is the image of a real this agrees with the previous usage of \underline{X} and \overline{X}. Arithmetic operations on representable reals and intervals are interpreted according to IEEE and interval rules, not by the usual mathematical definitions.

As an example, the sum of two intervals A = [\underline{A}, \overline{A}] and B = [\underline{B}, \overline{B}] is given by

A + B = [\underline{\underline{A} + \underline{B}}, \overline{\overline{A} + \overline{B}}],

where the left endpoints are rounded down when added, and the right endpoints are rounded up when added.

Other operations on reals are extended to intervals as in [29] and [2]. Relations are interpreted positively, so for example when X and Y are intervals, X < Y means that for any x ∈ X and any y ∈ Y, x < y, and X ≠ Y means that for any x ∈ X and any y ∈ Y, x ≠ y.

We have equivalent expressions X < Y ⇔ \overline{X} < \underline{Y}, and X ≠ Y ⇔ X ∩ Y = ∅. The union of intervals X and Y is somewhat nonstandard; it is the smallest interval containing both sets, and is denoted X\overline{∪}Y. The interval valued functions Absolute Value(X), Max(X, Y), Min(X, Y), and the intersection of two intervals X ∩Y are defined in the standard way, without rounding. There is no IEEE standard for transcendental functions, so we designed our program to avoid all calls to trigonometric functions.

The operations of the IEEE floating point standard can lead to undefined operations. For example the quotient 0/0 and the product 0 × +∞ result in an output of NAN (Not A Number), and are signaled by the presence of an exception flag. Operations of interval arithmetic can also lead to questionable operations, for example when the quotient of two intervals A/B is calculated, and B contain 0. Our implementation of the division operator on intervals returns [−∞, +∞] if B contains 0. The algorithm draws no conclusions when such a division occurs, but rather calls for a subdivision of the input into intervals of smaller size, where the operations are repeated with greater accuracy. Exceptions such as overflow and underflow, sometimes an issue in computer assisted proofs, are not an issue in our algorithm. In any case they do not occur.

Bounds obtained from repeated applications of interval arithmetic are potentially far from sharp, especially when using wide input intervals. If just a single value needs to be computed rigorously, one can work with intervals whose width is comparable to the roundoff error of a single operation. For our purposes these types of intervals are far too thin, and would make the number of calculations we need impractical. Our intervals are “fat” in the terminology of Fefferman [15], meaning that they are sized by the scope of the problem we are solving rather than the size of the computational rounding. Often the intervals are wide enough to make a perfectly good formula look like nonsense. Sometimes it is possible to narrow the intervals under consideration because we have knowledge about what the possible legitimate values are that could arise during a computation. For example, consider a formula involving square roots. Assuming we’re not using imaginary numbers, the validity of the formula presupposes that the argument to the square root is nonnegative. However, when we pass to an interval extension of the formula, the interval argument to the square root will often include negative values, even though these values cannot arise from the problem we are interested in. The square root of a negative number would normally trigger an exception, but the nature of the calculation allows us to define the interval square root function Sqrt to discard any negative portion of an interval argument. The justification for this apparent sleight-of-hand rests on the validity of the original formula over the reals. If our theorems tell us that the quantity whose square root we are taking is non-negative, we are justified in truncating the interval to exclude negatives.

The interval function which always returns [−∞, +∞] can be used to represent any function, but not very usefully. We use this interval to return the value of a division by an interval containing zero. An interval [a, b] is empty if a > b. All of our functions take nonempty interval arguments and return nonempty interval results except for the intersection operation, which returns an interval which can be empty.

The value of the nonrepresentable constant \sqrt{3} in our algorithm is expanded to a narrow interval. Other constants used are representable. Thus we use the expression 5H_{o} ≤ 1 rather than H_{o} ≤ 0.2, which involves the non-representable constant 0.2. Representable constants such as 5 behave identically in interval arithmetic to intervals where both endpoints have the same value, such as [5, 5], and operations involving both intervals and representable constants are treated as if the constant was first converted to a one-point interval.

The verify statement is used to assure that a particular condition holds. If the condition fails, then the entire program is stopped (aborted). If no such failure occurs, the presence of such a verify statement constitutes a proof that the condition holds for the given inputs. One use of the verify statement is in the procedure avgwt defined below, which gives a weighted average of two intervals. To apply this we want to be sure that, as we expect in the cases where it is applied, the first interval is strictly smaller than the second. We check that this is so with a verify statement. Running our program results in no violations of any statement which is tested with a verify statement.

The program examines the set of all torus bubbles to see if any can be a minimizer. A range of hypothetical torus bubbles is specified by intervals Θ_{1} and H_{o}, according to Proposition 4.13. To avoid unnecessary use of trigonometric functions our algorithm uses an equivalent parameterization of the space of torus bubbles by intervals Y_{1} and H_{o}, where Y_{1} = sin Θ_{1}. Rather than solve the mean curvature differential equation directly to find intervals (X_{2}, Y_{2}) containing the point (x_{2}, y_{2}) where τ_{i} and τ_{o} intersect, we assume existence of a torus bubble, derive Y_{2}, and then deduce X_{2} from numerical integrals. Volume calculations for V_{i} and V_{o} are obtained from additional numerical integrals. At each stage in the calculation, we check whether the torus bubbles can be rejected based on the instability results of Section 4. If not, we test whether the volumes are equal. As it turns out all torus bubbles are rejected for these reasons, so area calculations are not necessary and are not present in the algorithm we present. Most torus bubbles are rejected either because Θ_{2} is out of range, or because the x-displacement of τ_{i} and τ_{o} differ, or because the volumes enclosed in the two regions differ.

The computation of Θ_{2} and Y_{2} is based on an analysis of the force function. The integration process used to calculate the value of X_{2} is complicated by the fact that the curves may not be graphs, so we have to choose an appropriate parameterization. The most convenient parameterization is in terms of y, because the ODEs satisfied by Delaunay curves involve only y, and because the y-coordinates of the endpoints are specified by Θ_{1} and Θ_{2}. This allows x-displacement and volume to be expressed directly as integrals in terms of y, making it unnecessary to actually generate an ODE solution which produces τ_{i} and τ_{o}.

The curves τ_{i} and τ_{o} are generally not graphs as functions of y, so the integrals we use have singularities at points where ẏ = 0. A change of variables resolves the singularity, as worked out in Proposition 3.7.

Another difficulty is that the geometry sometimes degenerates at the boundary of the regions we are examining. For example, we need to exclude torus bubbles with θ_{1} arbitrarily close to 0, but some of our formulas become singular when θ_{1} = 0, and don’t make sense there. If the interval Θ_{1} contains 0, then various intervals representing y-coordinates along the Delaunay curves will also contain 0, and the integrals give results which will not allow torus bubbles to be eliminated. We get around this problem by using crude estimates which apply in a small zone near the x-axis. Near the x-axis we have some information about the initial angle and the concavity of the Delaunay curves, and this suffices to give crude bounds to their behavior near the axis. We compute integrals only when at an appropriate distance above the x-axis. This approach suffices for all cases except for when h_{o} is also near 0. In that case we rely on Proposition 4.21.

We now describe in pseudo-code the main interval procedures used in our algorithm. As noted above, aside from the procedures described in this section our algorithm makes use only of standard arithmetic operations on intervals, namely addition, subtraction, multiplication, division, absolute value, maximum, minimum and square.

The procedure avgwt calculates a number in between two inputs, used later to subdivide intervals. It returns a representable real number which is a weighted average. For convenince, we set the rounding mode to round to nearest in this procedure. In some applications of this procedure, an interval is assigned the output of the avgwt procedure, eg in Step 7 of the procedure DivideAnd-CheckRectangle. The resulting interval then consists of a single point.

real procedure avgwt
input:  X, Y, w

      verify X < Y

z := (1 − w)\overline{X} + w\underline{Y}

      return z

The interval procedure Sqrt(X) returns an interval containing \{\sqrt{x} : x ∈ Xandx ≥ 0}. The real procedure width calculates an upper bound for an interval’s width.

The interval procedure Compare checks whether the interval X is to the left of the interval Y, returning A if yes, B if X is to the right of Y, and the union A\overline{∪}B if the intervals overlap.

interval procedure Compare

input: X, Y, A, B

      if X < Y then return A
      else if X > Y then return B
      else return A\overline{∪}B

Lemma 5.1 For any x, y, a, b, in any X, Y, A, B, respectively, if x ≤ y then a ∈ Compare(X, Y, A, B), and if x ≥ y then b ∈ Compare(X, Y, A, B).

Proof: The lemma follows from a straightforward case by case analysis. □

Integrate is the basic method for numerical integrals. It gives relatively wide intervals, but is good enough for our purposes. Its input is an interval valued function of one interval variable, and two intervals which give the limits of integration.

interval procedure Integrate
input: F, A, B

R := 0
H := (\underline{B} − \overline{A})/32

      if (\overline{H} > 0) then

            begin
            for i := 0 to 31 do

R := R + F (\overline{A} + (0\overline{∪}H) + iH)

            end

return 0\overline{∪}F (A)Width(A) + 0\overline{∪}F (B)Width(B) + RH

Proposition 5.2 For any integrable f and real numbers a, b, with a ≤ b, and interval extensions F, A, B, the integral of f from a to b is contained in Integrate(F, A, B).

Proof: If A ≤ B,

∫_{a}^{b} f dx = ∫ \overline{A} ; a f dx + ∫ \underline{B} ; \overline{A} f dx + ∫ b ; \underline{B} f dx .

The term ∫ \overline{A} ; a f dx is contained in 0\overline{∪}F (A)Width(A) and similarly ∫ b ; \underline{B} f dx is

contained in 0\overline{∪}F (B)Width(B). Note that we need to replace F (A) and F (B) by 0\overline{∪}F (A) and 0\overline{∪}F (B) to deal with the issue that a can lie anywhere in A and similarly b in B. Taking the union of F with 0 ensures that our sum contains the intervals F ([a, \overline{A}]) and F ([\underline{B}, b]).

The value of ∫ \underline{B} ; \overline{A} f dx is bounded by lower and upper Riemann sums, using 32

equally sized intervals. The term F (\overline{A}+(0\overline{∪}H)+iH) evaluates F on an interval which contains [\overline{A} + ih, \overline{A} + (i + 1)h] for any h ∈ H, and so contains both an upper and a lower bound for F on this interval. The RH term in the last line, when added to the previous terms, then gives an interval containing both lower

and upper Riemann sums for ∫_{a}^{b} f dx.

Otherwise, A ≤ B is false, and since a ≤ b, A and B must overlap, so we can decompose

∫_{a}^{b} f dx = ∫_{a}^{c} f dx + ∫_{c}^{b} f dx

where a ≤ c ≤ b and c ∈ A ∩ B. Since \underline{B} < \overline{A}, R = 0 and the RH term in the last line vanishes. The result follows from [a, c] ⊂ [a, \overline{A}], [c, b] ⊂ [\underline{B}, b] and inclusion monotonicity. □

The next several procedures define functions to be integrated.

The procedure Dx is a straightforward interval extension of the expression for dx/dy given in Equation 6 in Proposition 3.6. It will be used to find the x-displacement from (x_{1}, y_{1}) to (x_{2}, y_{2}) by applying Equation 6 to both τ_{i} and τ_{o}.

Volumes under Delaunay surfaces are computed using the volume integrand Dv, which is obtained from Dx as in Proposition 3.6. It is convenient to have Dv calculate volume divided by π. Since volumes are only compared to one another, we can avoid the unnecessary step of multiplying by the non-representable constant π. The volume element computation required that dx/dy > 0 for the volume to have a positive sign. The curve τ_{o} will sometimes start out in the negative x direction, when θ_{1} < 60^{∘}, and after passing through a vertical tangency switch to going in the positive x direction. Calculating the integral of Dv in this case will give the volume under the part of the nodoid generated by the arc going right minus the volume under the arc going left. This is just what we want to compute the volume V_{o} in the torus component in the case when such an overhang exists.

At local minima and maxima, Dx and Dv are singular. Ymin and Ymax are global variables which are accessed by Dxmin, Dxmax, Dvmin and Dvmax. So near local extrema, we use a change of variables to reparameterize the curve, as in Proposition 3.7. Procedure Dxmin calculates dx/dz near a minimum and procedure Dxmax near a maximum, where z is given in Proposition 3.7. We do an algebraic manipulation to make the denominator in the square root more computationally effective, namely we replace the term HY − 2 + HYmax in procedure Dxmax by the equivalent expression H(Y + Youtermin). Here Youtermin represents the minimum y-value of the nodoid which has a maximum at Ymax. The equivalence of the two expressions is derived using Proposition 3.4 as follows:

HY − 2 + HYmax = \frac{2}{Ymax − Youtermin} (Y + Ymax) − 2 = H(Y + Youtermin)

The volume calculation procedures Dvmin, Dvmax give the formula for the volume integrand used near a critical point. The volume calculations are applied only to graphs over the x-axis.

interval procedure Dx
input: Y, H, F

T := HY^{2} − F
return T/Sqrt((2Y + T)(2Y − T))

interval procedure Dxmin
input: Z, H, F

Y := Ymin + Z^{2}
T := HY^{2} − F
return 2T/Sqrt((2Y + T)(2 − HYmin − Y H))

interval procedure Dxmax
input: Z, H, F

Y := Ymax − Z^{2} T := HY^{2} − F

      return 2T/Sqrt((2Y + T)H(Y + Youtermin))

interval procedure Dv
input: Y, H, F

      return Y^{2} Dx(Y, H, F)

interval procedure Dvmin
input: Z, H, F

      return (Ymin + Z^{2})^{2} Dxmin(Z, H, F)

interval procedure Dvmax
input: Z, H, F

      return (Ymax − Z^{2})^{2} Dxmax(Z, H, F)

The next procedure is used for testing interval rectangles of torus bubbles to see if they are potential minimizers and rejecting these intervals if they can be shown to have a property which rules out any torus bubble contained in them. The input is a pair of intervals C_{1} and H_{o}, which give a range of possible values for cos θ_{1} and h_{o}. The interval of cosine values C_{1} is used, rather than an interval of θ_{1} values, to eliminate calls to trigonometric functions. A variety of tests are performed based on the analysis of previous sections. The various intervals do not necessarily give sharp estimates, so it is not a priori clear how many rectangles will be rejected. As it turns out, the tests are sharp enough to reject all rectangles.

Boolean procedure CheckRectangle
input: C_{1}, H_{o}

1. if 1000C_{1} ≥ 996 and 5H_{o} ≤ 1

then return REJECT

2. H_{i} := 2 − H_{o} Y_{1} := Sqrt(1 − C_{1}^{2}) F_{i} := (H_{i} − 1)Y_{1}^{2} − C_{1}Y_{1}Sqrt(3) F_{o} := −F_{i} C_{2} := [−C_{1}\overline{∪}C_{1}] ∩ [−.5\overline{∪}.5] T := (2(H_{i} − 1)F_{i} + 3)/(3 + (H_{i} − 1)^{2}) − (1 − C_{1}^{2})

if (C_{2}^{2} + T) ≠ 1 return REJECT

3. Y_{2} := Sqrt(T ∩ [0, 1])

C_{2} := C_{2} ∩ \frac{(H_{i} − 1)Y_{2} − F_{i}/Y_{2}}{Sqrt(3)}

if C_{2} is empty then return REJECT

4. if C_{1} ≤ .5 and H_{o} ≤ 1−Sqrt(3)C_{1}/Y_{1}

      then return REJECT
5.    if Width(C_{2}) > .5 then return NORESULT

6. Wends := (1 − C_{1})^{2}(2 + C_{1})/3 + (1 − C_{2})^{2}(2 + C_{2})/3 Ymin := −F_{i}/(1+ Sqrt(1 + F_{i}H_{i})) Ymax := (1+Sqrt(1 + F_{o}H_{o}))/H_{o} Y := Compare(C_{1},Sqrt(3)/2, Ymin, Y_{1})

if Y H_{i} < −1 then
      begin

R := 1/(−H_{i} − 1/Y)

W := 2.5(Y + (Sqrt(3)/2)R)R^{2}

      if W < Wends then return REJECT
      end
7.  Yleft := Sqrt(F_{o}/H_{o})
if (Ymin < Y_{2} and Yleft < Ymax) then Yleft := Max(Y_{1}, Yleft)
      else return NORESULT
Y_{4} := avgwt(Yleft, Ymax, .5)

Z_{2} := Sqrt(Ymax − Y_{2}) Z_{4} := −Sqrt(Ymax − Y_{4}) if 1000C_{1} ≥ 998 then

            begin

T := avgwt(Y_{1}, Y_{2}, 1/16)

Δ_{i} := (T − Y_{1})(33/16)+Integrate(Dx(·, H_{i}, F_{i}), T, Y_{2}) Δ_{o} := −(Yleft − Y_{1})Sqrt(3) T := avgwt(Yleft, Y_{4}, 1/16) Δ_{o} := Δ_{o}+Integrate(Dx(·, H_{o}, F_{o}), T, Y_{4})

            if Δ_{o} > Δ_{i} then return REJECT

Δ_{o} := Δ_{o}+Integrate(Dxmax(·, H_{o}, F_{o}), Z_{4}, Z_{2})

            if Δ_{o} > Δ_{i} then return REJECT
            if 1 ∈ C_{1} then return NORESULT
            end
8.    verify (C_{1} ≤ Sqrt(3)/2 or Y_{1} ≤ Y_{2})

T := Sqrt(Y_{1} − Ymin)

Z_{1} := Compare(C_{1},Sqrt(3)/2, −T, T)

Z_{3} := Sqrt(Y_{2} − Ymin)

Δ_{i} := Integrate(Dxmin(·, H_{i}, F_{i}), Z_{1}, Z_{3}) Δ_{o} :=Integrate(Dx(·, H_{o}, F_{o}), Y_{1}, Y_{4})

      if Δ_{i} < Δ_{o} then return REJECT

Δ_{o} := Δ_{o}+Integrate(Dxmax(·, H_{o}, F_{o}), Z_{4}, Z_{2})

      if Δ_{i} ≠ Δ_{o} then return REJECT
9.  Wbase := Integrate(Dvmin(·, H_{i}, F_{i}), Z_{1}, Z_{3})
      W_{i} := Wends + Wbase

W_{o} :=Integrate(Dv(·, H_{o}, F_{o}), Y_{1}, Y_{4})+

            Integrate(Dvmax(·, H_{o}, F_{o}), Z_{4}, Z_{2}) − Wbase
      if W_{i} ≠ W_{o} then return REJECT
            else return NORESULT

Theorem 5.3 If CheckRectangle(C_{1}, H_{o}) returns REJECT without causing an IEEE exception, then there is no area-minimizing equal-volume torus bubble with cos θ_{1} ∈ C_{1} and h_{o} ∈ H_{o} .

Proof: The procedure CheckRectangle is fed a pair of intervals containing a range of values for cos θ_{1} and h_{o}, these intervals being denoted by C_{1} and H_{o}. We will call such a range of input values an input rectangle. If every pair of values (cos θ_{1}, h_{o}) ∈ C_{1} × H_{o} in an input rectangle can be ruled out as a possible minimizing torus bubble by some test, then the procedure returns the value REJECT. Otherwise it returns the value NORESULT (which results in a subdivision of the input rectangle into four smaller rectangles.) The tests applied to the rectangles of values for cos θ_{1} and h_{o} are based on the results developed in this paper, as we will explain below. At certain points we need to check whether some hypothesis required in our tests are satisfied. We can do so by using a verify statement. This statement will cause the program to

Figure 17: Compared volumes.

terminate if it is not satisfied. The proof proceeds by interpreting and justifying each step of the algorithm.

Step 1 rejects torus bubbles based on Proposition 4.21. C_{1} is a range of cosine values for the angle θ_{1}, so 1000C_{1} ≥ 996 corresponds to cos θ_{1} ≥ .996, which implies that θ_{1} ≤ 5.2^{∘} and the Proposition applies.

Step 2 first defines some values used in further steps. The values of F_{o} and F_{i} are calculated using the formulas in Proposition 4.16. It then examines the possible values of y_{2} = sin θ_{2}. Corollary 4.10 implies −.5 ≤ cos θ_{2}, and Lemma 4.8 implies cos θ_{2} ≤ .5. Since we assumed that θ_{1} ≤ θ_{2}, we have cos θ_{2} ≤ cos θ_{1}. Proposition 4.14 implies that − cos θ_{1} ≤ cos θ_{2}. Thus the values for cos θ_{2} arising from possible minimizing torus bubbles lie in the intersection of the intervals [−.5, .5] and − cos θ_{1}, cos θ_{1}.

The force equation for T_{i}

f_{i}(c) = (h_{i} − 1)(1 − c^{2}) − c\sqrt{3(1 − c^{2})}

has at most two roots by Lemma 4.15, one of which is c_{1}. The other one is c_{2}. Letting y = \sqrt{1 − c^{2}} and rearranging gives

f_{i} = (h_{i} − 1)y^{2} + cy\sqrt{3}
(f_{i} − (h_{i} − 1)y^{2})^{2} = 3y^{2}(1 − y^{2}) .

This is a quadratic in y^{2}. The sum of the two roots in y^{2} is

\frac{2(h_{i} − 1)f_{i} + 3}{3 + (h_{i} − 1)^{2}}

and the known root is 1 − c_{1}^{2}, so subtraction gives the other root. Thus the interval T contains y_{2}^{2}, the square of the second root. We next add the interval C_{2}^{2}, containing cos θ_{2}^{2}, to the interval T containing sin θ_{2}^{2} to see if they contain values adding to one. If not, there is no minimizing torus bubble in the chosen input rectangle C_{1} × H_{o}, and this input rectangle is rejected.

Step 3 calculates the possible values of x_{2} = cos θ_{2} corresponding to y-values whose square is in the interval T. The only possible values for y_{2} are in [0, 1], so we intersect T with [0, 1] before taking the square root of T to get an interval Y_{2} which contains possible values of y_{2}. Since

f_{i} = (h_{i} − 1)(1 − c_{2}^{2}) − c_{2}\sqrt{3(1 − c_{2}^{2})}

by Proposition 4.16, and 1 − c_{2}^{2} = y_{2}^{2}, we can solve to get an expression for c_{2} in terms of h_{i}, f_{i} and y_{2}:

C_{2} = ((H_{i} − 1)Y_{2} − F_{i}/Y_{2})/Sqrt(3) .

Note that this expression is chosen to give us c_{2} with the correct sign. We intersect the resulting interval with the interval C_{2} computed in the previous step. If the resulting interval is empty, there is no minimizing torus bubble in the chosen input rectangle and it is rejected.

Step 4 rejects torus bubbles based on Proposition 4.18. It uses the fact that cot(θ_{1}) = c_{1}/y_{1}.

Step 5 passes on a request to subdivide the rectangle if it is too wide. It does not reject anything, so there is nothing to justify. This step is included to speed the program’s execution.

Step 6 calculates rough bounds to check if the volumes of the two regions can be equal. We calculate volume divided by π in the algorithm to avoid the needless step of multiplying by the non-representable constant π. The interval Wends = Vends/π is defined to be the volume of the solid of revolution between the spherical ends and the x-axis divided by π. See Figure 17. This solid is wholly contained in the ball component. Ymin contains the y-coordinate of the local minimum of the τ_{i}, based on solving the force equation as in Part 7 of Proposition 4.12. This minimum is realized in the torus bubble if and only if θ_{1} ≥ 30^{∘}. Similarly Ymax contains the maximum height for τ_{o}, as in Part 8 of Proposition 4.12. It is always achieved. Step 6 goes on to reject potential torus bubbles if a calculation shows that the torus component has less than half the total volume. The volume of the ball component is bounded below by the volume under the spherical caps, whose value is given by Vends. The interval Y is defined to ensure that it contains the minimal y-coordinate of the torus component, which occurs either in Ymin (if θ_{1} ≥ 30^{∘}) or in Y_{1}. An upper bound V for the torus component volume is calculated using the formula derived in Lemma 4.24.

Step 7 tries to eliminate torus bubbles with small θ_{1}, namely those with cos θ_{1} ≥ .998, a condition which forces θ_{1} to be less than 3.7^{∘} and allows Lemma 4.26 to be applied. These small angles must be treated differently because the equations describing Delaunay curves are not well behaved near θ_{1} = 0. Bounds on the x-displacement x_{2} − x_{1} are calculated using each of τ_{i} and τ_{o}. The interval Δ_{i} contains an upper bound for the width of τ_{i} and the interval Δ_{o} contains a lower bound for the width of τ_{o}. If the resulting intervals satisfy Δ_{i} < Δ_{o} then τ_{i} does not travel as far to the right as τ_{o} and the input rectangle can be rejected.

Figure 18: A double bubble configuration for 0 < θ_{1} ≤ 30.

Figure 19: A double bubble configuration for 30 ≤ θ_{1} ≤ 60.

Step 7 begins by defining Yleft to be an interval containing the y-value of the point where τ_{o} has a vertical tangency on the left as in Figure 18 and Figure 19. This value is computed by taking the square root of the nonnegative part of F_{o}/H_{o}. Corollary 3.3 guarantees that this value is positive in a torus bubble, so we are justified in truncating negative values which may have appeared in the calculation before taking the square root. The outer curve is always a nodoid by Corollary 4.19, and if extended to a vertical tangency, the y-coordinate of that point is in the interval Yleft. Replacing Yleft with Max(Yleft, Y_{1}) assures that Yleft contains the y-coordinate of a vertical tangency if there is one, and the initial point otherwise. In particular it always contains the y-coordinate of

Figure 20: A double bubble configuration for 60 ≤ θ_{1} ≤ 90.

the leftmost point.

The value Y_{4} is used as an intermediate point for chopping up the intervals of integration for τ_{o} between its vertical and horizontal tangencies. See Figures 18, 19, 20. The limits of integration Z_{2} and Z_{4} are defined for use when integrating near a maximum of τ_{o}, as in Proposition 3.7. The square roots taken in these calculations are justified in discarding the negative part of any interval, because for any particular torus bubble the quantity whose root is being taken is nonnegative.

An upper bound for the x-displacement of τ_{i} as it goes from y_{1} to y_{2} is given by Δ_{i}, and a lower bound for the x-displacement of τ_{o} as it goes from y_{1} to y_{2} is given by Δ_{o}. If Δ_{i} < Δ_{o}, then a torus bubble cannot occur in the given input rectangle.

We use the Intermediate Value Theorem to estimate Δ_{i} up to height y_{1} + \frac{y_{2} − y_{1}}{16} and Δ_{o} up to height yleft. The estimates for Δ_{i} and Δ_{o} near the x-axis are justified by Proposition 3.4 in the case of τ_{o} and by Lemma 4.26 in the case of τ_{i}, as shown in Figure 21.

For θ_{1} < 4^{∘}, τ_{i} is strictly increasing, so that the x-displacement is equal to the integral of Dx from y_{1} to y_{2}. This can be broken up into an integral from Y_{1} to T and an integral from T to Y_{2}, where T is set to be one sixteenth of the way from Y_{1} to Y_{2}. The values in T are all less than 1/2, so that Lemma 4.26 applies and the slope of τ_{i} below height T is smallest at its initial point (x_{1}, y_{1}). It follows that

dx/dy ≤ cot (30 − θ_{1}) ≤ cot 26^{∘} < 33/16

on this interval and the integral from Y_{1} to T is bounded above by (33/16)(T − Y_{1}). Then Δ_{i} = (T − Y_{1})(33/16)+Integrate(Dx(·, H_{i}, F_{i}), T, Y_{2}) gives an upper bound for the width of τ_{i} as it increases from Y_{1} to Y_{2}.

With τ_{o} we first find an upper bound on the leftward displacement between y_{1} and yleft using a similar method. In this case the slope has absolute value larger than tan 30^{∘} = 1/\sqrt{3} and (Yleft − Y_{1})Sqrt(3) gives an upper bound for the leftward displacement to height yleft. An integral is calculated which gives a lower bound for the rightwards displacement of τ_{o} between yleft and y_{2}, by Figure 21: Estimating the widths of τ_{i} from above and τ_{o} from below with θ_{1} small.

integrating dx/dy from a point strictly above yleft to y_{2}. This integral is calculated in two pieces, the first to Y_{4} and the second from Y_{4} through the local maximum to Y_{2}. A lower bound for Δ_{o} is gotten by taking the upper bound for the leftwards displacement previously calculated and subtracting it from the sum of these integrals. If Δ_{o} > Δ_{i} the input rectangle is rejected.

This step can reject torus bubbles even if the interval C_{1} contains the value 1, corresponding to θ_{1} = 0. If this step test fails to reject a range of potential torus bubbles, we check whether 1 is in the interval C_{1}. If it is, then the next steps will not help, since they are ineffective for θ_{1} = 0. So if 1 is in C_{1} we return NORESULT and send the current input rectangle back for subdivision.

Step 8 calculates intervals for the x-displacement of τ_{i} and τ_{o} in the case where the angle is larger than those treated in the previous step, namely cos θ_{1} ≤ .996. In this step Δ_{o} is an interval containing the x-displacement of the outer curve and Δ_{i} an interval containing the x-displacement of τ_{i}. In contrast to Step 7, they do not contain bounds for these displacements, but the actual widths of τ_{i} and τ_{o}. Thus we can reject the input rectangle in this step if Δ_{i} and Δ_{o} do not overlap.

The interval Δ_{o} is calculated by integrating Dx from Y_{1} to an intermediate Y_{4} and then using Dxmax to integrate from Y_{4} across a maximum to Y_{2}. If the width of Δ_{o} is larger than 20 then the procedure returns NORESULT, calling for a subdivision of the input rectangle of intervals. This is done because otherwise large amounts of computational time are wasted in inefficient calculations.

The calculation of Δ_{i} is complicated somewhat because we need to deal with both the case where θ_{1} ≥ 30^{∘} and the case where θ_{1} ≤ 30^{∘}. Only in the former case, where C_{1} ≤ Sqrt(3)/2, does τ_{i} achieves an interior minimum. In that case, an integral along τ_{i} can be calculated by integrating Dxmin from Z_{1} to Z_{3}. If an interior minimum is not achieved, we can still integrate Dxmin to calculate the x-displacement from Y_{1} to Y_{2}, but need to reverse the sign of Z_{1} as in the last statement of Proposition 3.7. Our integral formulas for x-displacement need to have that Y_{1} ≤ Y_{2} in this case, so we check with a verify statement that either this holds or there is an interior minimum. The Compare statement ensures the correct sign is given to Z_{1} unless C_{1} overlaps the interval Sqrt(3)/2. Z_{1} is - Sqrt(Y_{1} − Ymin) if τ_{i} achieves an interior minimum and Sqrt(Y_{1} − Ymin) otherwise. The appropriate sign is given to Z_{1} using a Compare statement. If C_{1} overlaps the interval Sqrt(3)/2 we can’t tell which sign is correct, so the Compare statement returns an interval containing both Z_{1} and −Z_{1}, and the integral then contains the x-displacements for both cases. If Δ_{i} and Δ_{o} are disjoint, then the curves’ final points cannot coincide and the input rectangle is discarded.

Step 9 calculates volumes divided by π, as in Figure 17. We again calculate volume divided by π in the algorithm to avoid an unnecessary multiplication of all volumes by π. V_{base} is the volume surrounded by T_{i} and W_{base} = V_{base}/π. Adding Wbase to Wends gives W_{i} = V_{i}/π). The calculation of V_{o} is slightly more complicated because the overhang on the left involves a subtraction, if there is an overhang. Such an overhang occurs when θ_{1} < 60^{∘} as shown in Figures 18 and 19. The formulas are set up so that volume is counted with a negative sign when τ_{o} is oriented to the left, and positive sign when τ_{o} is oriented to the right, as discussed in Proposition 3.6. This is just what is required to calculate the volume of the overhang. The volume Vbase is subtracted from the volume inside T_{o} to get V_{o}, the volume of the torus component. Step 9 then rejects an input rectangle if the values W_{o} = V_{o}/π and W_{i} = V_{i}/π are unequal. □

Now that we have tests that are able to reject certain ranges of torus bubbles, we apply these tests to see whether they in fact reject all the possibilities. Since we don’t know in advance how fat the ranges can be, we initially feed in a rectangle of values with 0 ≤ sin θ_{1} ≤ 1, 0 ≤ H_{o} ≤ 10. We recursively subdivide further as necessary. We monitor IEEE exception flags at this level, so that no computation is trusted if it raised an exception.

Boolean procedure DivideAndCheckRectangle
input: Y_{1}, H_{i}

C_{1} := Sqrt(1 − Y_{1}^{2})

      call CheckRectangle(C_{1}, H_{o})
      if result is REJECT then return SUCCESS
      split Y_{1} in half, into Y_{1a}, Y_{1b}

      split H_{o} in half, into H_{oa}, H_{ob}
      call DivideAndCheckRectangle(Y_{1a}, H_{oa})
      call DivideAndCheckRectangle(Y_{1a}, H_{ob})
      call DivideAndCheckRectangle(Y_{1b}, H_{oa})
      call DivideAndCheckRectangle(Y_{1b}, H_{ob})
      return SUCCESS

Main program
      begin
      clear exceptions
            begin

Y_{1} := [0, 1] H_{o} := [0, 10]

            Call DivideAndCheckRectangle(Y_{1}, H_{o})
            end
      verify no exceptions raised
      Print ``All torus bubbles rejected.''
      end

Theorem 5.4 The algorithm described in Main, if run to completion without causing an exception, shows that no area minimizing torus bubble can enclose equal volumes.

Proof: By Proposition 4.13, any minimizing torus bubble is determined by θ_{1} and h_{o}, where θ_{1} is the angle of first arc and h_{o} is the mean curvature of T_{o}. It suffices to consider 0 ≤ θ_{1} ≤ 90^{∘} by Proposition 4.14, and 0 ≤ h_{o} ≤ 10 by Proposition 4.3 and Proposition 4.25.

Main calls DivideAndCheckRectangle with a rectangle containing the parameters of all possible torus bubbles. DivideAndCheckRectangle subdivides this rectangle, depending on the success of its calculations. To avoid unnecessary use of trigonometric functions in the computation, the space of angles is parameterized by a variable y_{1} = sin θ_{1}. Thus the intervals Y_{1} cover the interval [0, 1], in 1-1 correspondence with 0 ≤ θ_{1} ≤ 90^{∘}. DivideAndCheck-Rectangle begins by using Y_{1} to compute C_{1}, an interval corresponding to values of the cosine of an interval of θ_{1} angles. It then passes the C_{1} and H_{o} intervals to CheckRectangle where they are tested as previously described.

Main can only finish if CheckRectangle rejects all rectangles. By Theorem 5.3, torus bubbles rejected by CheckRectangle are not minimizers. □

We implemented the algorithm described above in C++, though we could equally well have used any language which allows for control of rounding and follows the IEEE 754 standard, or any similar precise prescription for floating point arithmetic. We ran code implementing this algorithm on a variety of machines, various Wintel, Sun, HP, SGI, UNIX and Linux platforms, with identical results. The C++ code with embedded instructions is available at the web site [21]. The code takes about 10 seconds to run on a fast (1999) PC. The input rectangle was recursively subdivided up to ten times, resulting in a total of 15,016 rectangles examined. The calculation involved a total of 51,256 integrals.

From examining the output of the program, we can see the reasons for which various rectangle were rejected. Of the 15,016 rectangles examined, Step 1 rejected 31, Step 2 rejected 1140, Steps 3, 4 and 5 rejected none, Step 6 rejected 3664, Step 7 rejected 541, Step 8 rejected 6691 and Step 9 rejected 2949. In Figure 22 we indicated which steps were used to reject torus bubbles in which input rectangles. This chart offers some insight into why torus bubbles fail to be minimizers. However there is some arbitrariness to this, as doing the steps in a different order would have produced different data.

Figure 22: Chart showing which tests were used in each input rectangle. In some rectangles there were several subdivisions, and several tests were used. Steps 3,4 and 5 were never used to reject input rectangles.

In running the program no exceptions were raised, and the program returned the statement “All torus bubbles rejected.” As a consequence we obtain a proof of our main result.

Theorem 1 The unique surface of least area enclosing two equal volumes in ℝ^{3} is a double bubble.

Proof: Theorem 5.4 eliminates the possibility of a minimizing torus bubble. Theorem 4.1 implies that the only other possibility is a double bubble.

Acknowledgements: We are indebted to Frank Morgan for introducing us to this problem, to John Sullivan for generously providing us with computer graphics of torus and double bubbles, and to Morgan, M. Hutchings, W. Kahan and W. Rossman for helpful discussions.

References

[1] M. Alfaro, J. Brock, J. Foisy, N. Hodges and J. Zimba The standard double soap bubble in R^{2} uniquely minimizes perimeter, Pac. J. Math. 159, 47-59 (1993).

[2] G. Alefeld and J. Herzberger, Introduction to interval computations, Academic Press, 1983.

[3] F.J. Almgren, Existence and regularity almost everywhere of solutions to elliptic variational problems with constraints, Memoirs Amer. Math. Soc. 4, 165-199 (1976).

[4] F.J. Almgren and J. Taylor, The geometry of soap films and soap bubbles, Sci. Amer. 235,82-93 (1976).

[5] ANSI/IEEE Standard 754-1985 for Binary Floating-Point Arithmetic, The Institute of Electrical and Electronic Engineers, New York, 1985.

[6] K. Appel and W. Haken, The solution of the four-color-map problem, Scientific Amer. 237, 108-121, (1977).

[7] M. Athanassenas, A free boundary problem for capillary surfaces, Manuscripta Math. 76, 5-19, (1992).

[8] J. L. Barbosa and M. do Carmo, Stability of hypersurfaces with constant mean curvature, Math. Z. 185, 339-353 (1984).

[9] M. Berger and B. Gostiaux, Differential Geometry, Manifolds curves and surfaces, Springer-Verlag Graduate Texts 115, (1988).

[10] C.V. Boys, Soap Bubbles, Dover Publ. Inc. NY 1959 (first edition 1911).

[11] Carthage, Encyclopedia Britannica, 19th Edition, 1929.

[12] R. de la Llave, Computer assisted proofs of stability of matter, Computer aided proofs in analysis (Cincinnati, OH, 1989), 116-126,

[13] C. Delaunay, Sur la surface de revolution dont la courbure moyenne est constante, J. Math. Pure et App. 16, 309-321 (1841).

[14] J. Eells, The surfaces of Delaunay, Math. Intelligencer 9, 53-57 (1987).

[15] C. Fefferman, The N-body problem in quantum mechanics, Comm. Pure and Applied Math. Supplement 39, S67-S109 (1986).

[16] J. Foisy, Soap bubble clusters in R^{2} and ℝ^{3}, undergraduate thesis, Williams College (1991).

[17] D. Gabai, R. Meyerhoff and N. Thurston, Homotopy Hyperbolic 3-Manifolds are Hyperbolic, MSRI Preprint 1996-058.

[18] R. Gulliver, Regularity of minimizing surfaces of prescribed mean curvature, Ann. of Math. (2) 97, 275–305 (1973).

[19] T. Hales, The Kepler conjecture, E-print math.MG/9811078,

[20] J. Hass, M. Hutchings and R. Schlafly, The Double Bubble Conjecture, ERA-AMS 1, 98-102 (1995).

[21] http://math.ucdavis.edu/~hass/doubbub.cpp

[22] M. Hutchings, The structure of area-minimizing double bubbles, J. Geom. Anal. 7, 285–304 (1997).

[23] K. Kenmotsu, Surfaces of revolution with prescribed mean curvature, Tohoku Math. J. 32, 147-153 (1980).

[24] Lord Kelvin, William Thomson, On the Division of Space with Minimum Partitional Area, Phil. Mag. 24, 503 (1887).

[25] N. Korevaar, R. Kusner and B. Solomon, The structure of complete embedded surfaces with constant mean curvature, J. Differential Geom. 30, 465-503 (1989).

[26] G. Lawlor and F. Morgan, Paired calibrations applied to soap films, immiscible fluids, and surfaces or networks minimizing other norms. Pacific J. Math. 166, no. 1, 55-83, (1994).

[27] O. Lanford, Proof of the Feigenbaum Conjectures, Bull. Amer. Math. Soc. 6, 427-434 (1982).

[28] H.P. McKean, M. Schreiber and G. H. Weiss, Isoperimetric problem with application to the figure of cells, J. Math. Physics 6, 479-484 (1994).

[29] R. E. Moore, Methods and Applications of Interval Analysis, SIAM, 1979.

[30] J.D. Moore and R. Schlafly, On equivariant isometric embeddings, Math Z. 173, 119-129 (1980).

[31] F. Morgan, Soap bubbles in R^{2} and in surfaces, Pac J Math 165 (1994), 347-361.

[32] F. Morgan, Clusters minimizing area plus length of singular curves, Math. Ann. 299, 697-714 (1994).

[33] R. Osserman, The isoperimetric inequality, Bull. Amer. Math. Soc. 84, 1182-1238 (1978).

[34] R.S. MacKay and I.C. Percival,Converse KAM: theory and practice, Comm. Math. Phys. 98, 469–512 (1985).

[35] J. Plateau, Statique experimentale et theorique des liquides soumis auz seules forces molecularies, Gathier-Villars, Paris, 1873.

[36] H. Sagan, Introduction to the Calculus of Variations, Dover, 1992.

[37] N. Robertson, D. Sanders, P. Seymour and R. Thomas, The four-colour theorem. J. Combin. Theory Ser. B 70 , no. 1, 2-44, (1997).

[38] M. Spivak, A comprehensive introduction to differential geometry, Volumes III and IV, Second Edition, Publish or Perish Inc, 1979.

[39] J. Taylor, The structure of singularities in soap-bubble-like and soapfilm-like minimal surfaces, Ann. of Math. 103, 489-539 (1976).

[40] D’arcy Thompson, On growth and form, Cambridge Univ. Press, NY, 1959, (first edition 1917).

[41] T. I. Vogel, Stability of a liquid drop trapped between two parallel planes, SIAM J. Appl. Math 47 (1987) 516-525.

Joel Hass Department of Mathematics University of California Davis, CA 95616

http://math.ucdavis.edu/~hass

e-mail: hass@math.ucdavis.edu

Roger Schlafly Department of Computer Science University of California Santa Cruz, CA 95064

http://bbs.cruzio.com/~schlafly

email: real@ieee.org