Geometry of catenoidal soap ﬁlm collapse induced by boundary deformation

Experimental and theoretical work reported here on the collapse of catenoidal soap ﬁlms of various viscosities reveal the existence of a robust geometric feature that appears not to have been analyzed previously; prior to the ultimate pinchoff event on the central axis, which is associated with the formation of a well-studied local double-cone structure folded back on itself, the ﬁlm transiently consists of two acute-angle cones connected to the supporting rings, joined by a central quasicylindrical region. As the cylindrical region becomes unstable and pinches, the opening angle of those cones is found to be universal, independent of ﬁlm viscosity. Moreover, that same opening angle at pinching is found when the transition occurs in a hemicatenoid bounded by a surface. The approach to the conical structure is found to obey classical Keller-Miksis scaling of the minimum radius as a function of time, down to very small but ﬁnite radii. While there is a large body of work on the detailed structure of the singularities associated with ultimate pinchoff events, these large-scale features have not been addressed. Here we study these geometrical aspects of ﬁlm collapse by several distinct approaches, including a systematic analysis of the linear and weakly nonlinear dynamics in the neighborhood of the saddle node bifurcation leading to collapse, both within mean curvature ﬂow and the physically realistic Euler ﬂow associated with the incompressible dynamics of the surrounding air. These analyses are used to show how much of the geometry of collapsing catenoids is accurately captured by a few active modes triggered by boundary deformation. A separate analysis based on a mathematical sequence of shapes progressing from the critical catenoid towards the Goldschmidt solution is shown to predict accurately the cone angle at pinching. We suggest that the approach to the conical structures can be viewed as passage close to an unstable ﬁxed point of conical similarity solutions. The overall analysis provides the basis for the systematic study of more complex problems of surface instabilities triggered by deformations of the supporting boundaries.


I. INTRODUCTION
The nature of singularity formation during the collapse of a soap film with boundary has been addressed in configurations ranging from the prototypical catenoid to less familiar one-sided and twisted surfaces [1][2][3][4][5][6][7][8]. Of particular interest has been the location of the associated singularity-in the bulk of the film or on the boundary-and how that location depends on the topology of the film. An emerging conjecture [5] is that the structure of the linear stability operator of the surface (through the associated Jacobi fields that become the unstable modes) determines the ultimate location of the singularity. Because the most fundamental bulk singularity arises from the collapse of catenoidal necks, it is natural to examine the stability problem of equilibrium catenoids to understand the approach to that singularity. Despite the long history of the study of stability of minimal surfaces, primarily from an energetic point of view [9], and even though such problems arise in many different contexts in fluid mechanics [10] and soft matter [11][12][13], there have been only a few discussions of the details of the linear and weakly nonlinear problems [2,[14][15][16]. And these have either focused on simplified laws of motion such as mean curvature flow, a model that has been extensively studied in the pure mathematics literature [17][18][19][20] without accounting for flow of the surrounding air, or considered a dynamics governed only by the inertia of the film [21,22], or did not provide systematic stability calculations.
Many years ago, Chen and Steen [2] outlined the characteristic shapes adopted by soaps films as they progress from an unstable catenoid to the topological transformation at pinchoff. Some of these, from experiments described below, are shown in Fig. 1, including the quasicatenoidal initial shape (i) that steadily contracts into (ii), then a shape with an approximately cylindrical central region (iii), and finally a well-known double-cone geometry folded back on itself (v) that leads to breakup. In the course of studying this problem we discovered that there is an additional, intermediate state shown in (iv), in which the transiently formed cylindrical region pinches at its two ends, leading to conical films that connect to each of the supporting circular frames, with a unique opening angle to the cones. These same features, with essentially the same cone geometry (iv), appear if the threedimensional catenoid is split in half with a glass plate along its axis of revolution, so that there are two hemicatenoids, each with a moving contact line on the glass plate subject to viscous drag ( Fig. 2) [23]. (In this case, since the moving contact line is the junction between a soap film and a solid, it is termed a Plateau border.) These findings raise several questions: What is special about shape (iii), composed of surfaces with vanishing Gaussian curvature, in the context of surface evolution? What determines the selected cone angle at pinching? To what extent can a weakly nonlinear analysis around the unstable catenoid account for the experimental observations, and the evolution toward the double-cone structure (iv)?
As noted above, the late stages of the breakup process that begins in panel (iv) of Fig. 1 leads to the configuration shown in panel (v), and exhibits a double cone. The geometry of this FIG. 2. Collapse dynamics of a three dimensional catenoid compared to that of a hemicatenoid in contact with a glass plate. (a) Time evolution of the scaled area A (where A = (film area)/2π R 2 and R is the hoop radius) as the film evolves from the original catenoid to two discs. Points labeled (i-viii) correspond to characteristic stages in the evolution shown in the respective images. (b) As in panel (a), for the hemicatenoid collapse. Note the much longer timescale for collapse compared to panel (a). 035105-2 configuration was experimentally discovered some years ago in the study of catenoid collapse [2,24] and soon thereafter studied theoretically [25,26] using an approach based on a combination of potential flow and similarity solutions [27] for the film shape. These theoretical results were inspired not only by experiments, but also by numerical computations of surface dynamics coupled to inviscid fluid motion in which the surface was closed, without boundaries. This is an appropriate approximation for the local problem of the double cone. However, earlier stages of the collapse, such as (ii) and (iii) in Fig. 1, require accounting for the fact that the film is open and thus air can escape the domain, as seen in numerical simulations [2]. Otherwise, instead of the conical shape seen in Fig. 1 (iii), a bulging quasihemispherical shape is seen, as found in studies of capillary bridge breakup [28] with impermeable boundaries pinning the two ends of the bridge.
In Sec. II we describe experimental studies of the collapse of catenoids under boundary deformation which reveal the geometrical features described above. Motivated by these observations, Sec. III presents an in-depth analysis of the stability of catenoids under several different laws of motion. After a short review of general results in the theory of surface motion and equilibrium catenoids, we develop the framework necessary to treat both linear stability problems and weakly nonlinear ones in the neighborhood of the critical catenoid. As an example we first apply this to the case of mean curvature flow (MCF) and then proceed to the realistic model in which the surface motion is coupled to that of the surrounding air. From these results it becomes clear that the complete set of radial modes obtained from the Sturm-Liouville problems within MCF is quantitatively very close to those within Euler flow (albeit with growth rates given by very different physics) and can serve as a means to understand the dynamical evolution of catenoids in a model-independent manner. Section IV develops this idea by a numerical projection of the experimental collapsing interface shapes onto those modes. We find that with even as few as four modes the overall shape of the interface is very accurately captured, and the experimental growth rate of the unstable mode is accurately captured by the Euler dynamics. In Sec. V we use these results to explain the experimental data on the time evolution of the film area as the catenoid evolves toward collapse. An alternative approach involving the construction of a mathematical sequence whose limit is the Goldschmidt solution [29] is used to derive the value of the critical cone angle at pinching. We also present a heuristic argument for this value by utilizing the geometry of the unstable mode and connect these results to the potential flow arguments of earlier work [25], leading to the suggestion that the conical film shape prior to the ultimate collapse implies that the film evolves close to an unstable fixed point of similarity solutions. We conclude in Sec. VI with a discussion of future work. Some calculational details are given in the Appendix.

II. EXPERIMENTS
Catenoids were obtained by pulling two coaxial circular rings of radius R = 4 cm out of surfactant solutions. Initially the half separation distance d is smaller than the critical distance d c and is raised slowly until the instability is triggered.
While not precisely controlled, the typical displacement of the rings from that of the last stable catenoid is 3-5 mm. For the case of half-catenoids, a glass plate that fits snugly inside the rings, and was pre-wetted with the surfactant solution, was slowly lowered along the axis of the catenoid to separate it into two identical hemicatenoids. The solutions were obtained by dissolving tetradecyl trimethyl ammonium bromide (TTAB) in deionized water and adding various amounts of glycerol, up to 85% in mass, to increase the viscosity. The concentration of TTAB was 3 g/l for the aqueous solution without glycerol and was raised to 6 g/l for that with glycerol to enhance the stability of the soap films. The viscosity of the solutions ranged from 1 to 77 mPa s, while the surface tension remained nearly constant (between 38 and 35 mPa m [30,31]). We found that the collapse dynamics were independent of the viscosity over the range examined, which emphasizes that the phenomenon is driven primarily by a balance between surface tension and inertia of the surrounding air [2,34]. As discussed elsewhere [32], we have also performed experiments with surfactant solutions having higher interfacial moduli [33]. In these cases, we have observed that both the dynamics and the patterns are very different from those presented here. In what follows, results are presented with surfactant solutions of viscosity 1.0, 4.3, and 77 mPa s. Fluorescein was added to the solutions to render them fluorescent under illumination by arrays of cyan LEDs as described previously [3]. The collapse dynamics were recorded with a color high-speed camera (Phantom V641, Ametek) at 4 000 or 5 600 frames per second. Typical images are shown in Figs. 1 and 2; a representative video of the collapse can be found in the Supplemental Material [35].
Taking advantage of its symmetry, we isolate one quadrant of the soap film for analysis. The contour shape defined by its radial coordinate ζ (z) as a function of the axial position z is discretized with 20 or more points as needed to keep computational errors below 2%. To make contact with analytical results described later, we define a rescaled interface shape r = ζ /a c and coordinate u = z/a c via the quantity a c = 0.5524R that appears in the theory of equilibrium catenoids (see below). Likewise, we define A to be the physical area normalized by 2π R 2 , a scaling such that A ≈ 1.2 for the critical catenoid (with d = d c ), while A = 1 in the configuration of two discs spanning the frames. Figures 2(a) and 2(b) show the collapse sequence for catenoids and hemicatenoids along with the time evolution of their respective areas A. Figure 3 shows A, the respective neck radii r n (defined as the minimum radius, found initially in the midplane, u n = 0, and later at secondary minima symmetrically displaced away at |u n | > 0) and tangent angles θ at the attachment point ( Fig. 1) for experiments with different film viscosities. Figures 2(a) and 3(a) show that A relaxes towards 1 in less than 30 ms. The neck radius decreases and reaches zero when the soap film pinches off. The angle θ increases rather regularly in time from 57 • (very close to 56.5 • = arctan(R/d c ) as expected at the onset of the instability [14]). At the end of the relaxation of the area, θ is still far from the 90 • expected for the two-disc solution; this value is reached on a much larger timescale. When the soap film pinches, we measure θ c = 71.5 • .
A more detailed examination of the film shapes is displayed in Figure 4. For each profile shown in Fig. 4(a), a linear interpolation is performed to give an estimate of the average angle in the r − u plane. The quality of that linear interpolation is quantified by the standard deviation of this angle, which is plotted as a function of time in Fig. 4(b), where we see a strong local minimum with standard deviation ∼ 0.1 • , which corresponds to an average angle θ * = 68 • that is slightly less than θ c due to the residual curvature of the interface. At this time, the top and bottom of the soap film are the most accurately conical, while the center section is a pinched cylinder; we call this the "Martini-glass" configuration. Note that this configuration occurs at the same time step as the one of the pinch-off. The two values θ c and θ * serve as lower and upper estimates of the Martini glass angle, which we take to be their mean, i.e., θ M ∼ 69.8 • .
In their classic work on soap film collapse, Chen and Steen [2] noted from numerical computations that the catenoid neck radius r n shown in Fig. 3(b) appears to collapse with the classic scaling expected by dimensional analysis for surface tension driven Euler flows [27], namely, that characteristic lengths evolve in time t relative to a singularity as (γ t 2 /ρ a ) 1/3 , where γ is the surface tension and ρ a is the density of the surrounding air. Actually this regime is observed twice: first after the location of the minimum radius bifurcates from its initial value of u = 0 to two symmetrically placed points above and below the midplane (named "roll-off", with the law r n ∝ t 2/3 attributed to be an artifact of the symmetry of the setup, with a stagnation point at u = 0), then in the final self-similar regime before pinching. In between there is a short regime following r n ∝ t 2/5 . In our experiments, quantifying the final self-similar regime is below the experimental resolution and the necks seem to pinch at some values ±u p straddling u = 0. Figure 5 shows that our data (obtained for several values of the film viscosity) for the pinching neck radius r p = r(u p ) displays that same exponent of 2/3. Thus, the approximate scaling holds away from the central stagnation point. Additionally, film viscosity plays only a minor role in the approach to pinching. We conclude that despite the fact that the ultimate topological transformation takes place later, through the appearance of the folded double cone structure, the approach to the single cone state exhibits Keller-Miksis scaling down to very small (but finite) length scales.

A. Geometrical preliminaries
We start with a brief review of the standard description of a surface of revolution S spanning two parallel circular where t is time,n is the surface's unit normal, and U is determined by the underlying physics. For a surface of revolution, we haven · ρ t = ζ t / 1 + ζ 2 z , so ζ evolves as The time evolution of geometrical quantities follows directly from these definitions, as presented briefly in earlier work [36] and in detail in the Appenedix. The area functional I[ζ ], written for convenience as 2π A, where and √ G = ζ 1 + ζ 2 z , has the time evolution where dS = √ Gdz is the element of surface area, we have used Eq. (1), and H is the mean curvature, The volume V = π d −d dzζ 2 enclosed by the surface evolves in time as As shown in the Appendix, H itself evolves as where is the Gaussian curvature and the covariant Laplace-Beltrami operator is In highlighting the appearance of a conical film shape we noted that a cone has vanishing Gaussian curvature. This raises the more general question of the spatiotemporal evolution of K during the collapse process. Figure 6 shows for one particular experiment the evolution of r(u) and K (u) up to the appearance of the cone. We see that K clearly develops a singularity as the neck radius tends toward zero at u ≈ 0.35, while at the same time K is very nearly zero over the region u 0.7, which is further evidence of the formation of a conical structure.
The singularity that develops in K arises from the collapsing radius of the film; by the (t c − t ) 2/3 scaling of lengths we expect K ∼ (t c − t ) −4/3 . But this must be an integrable singularity since there is a constraint on the total curvature K = dSK. For an up-down symmetric surface of revolution, direct integration of Eq. (8) yields where θ is defined in Fig. 1, and can alternatively be expressed as θ = π/2 − φ in terms of the tangent angle of the surface at the upper boundary (cot φ = ζ z ). When the catenoid is not bounded, but extends to infinity, then θ → π/2 and the well-known result K = −4π is recovered. As the metric factor dS ∼ ζ the integral within the total curvature scales as (t c − t ) −2/3 , and when integrated over a collapsing length scale leads to a finite result. These balances of terms can be seen in the graph of the product of K and the metric factor shown in Fig. 6(c).
The data on the time evolution of θ in Fig. 3(c) shows that it has only a modest change from its initial value to that when the cone appears. From the general expression for the time evolution of K in an axisymmetric setting, one deduces by direct integration that the evolution of K is given by the tangent angle and derivative of the normal velocity at the boundary, Figure 7 shows the slow variation of K averaged over four different values of the film viscosity. K changes relatively little from its initial value to that at the pinch time.

B. Equilibrium catenoids
For a surface to be minimal, the functional derivative of its area must vanish. From Eqs. (3) and (5) we may write For a surface of revolution spanning frames of equal radii R, the minimal surface condition H = 0 is solved by for any a, and the value associated with a particular physical situation is determined by the condition R = a cosh(d/a), that the film meets the frames at z = ±d. In our experiments the frame radius R is fixed and only d can be varied, so we rescale by R, setting The boundary condition can thus be written as As shown in Figs In the subcritical region the larger solution (α + ) is the stable one of interest. The rescaled areas A ± = A ± /R 2 of the two branches of solutions are These two branches meet at D c with area A c = 1.199 . . .
[ Fig. 8(c)]. Note that the quantities A ± are the ratios of the catenoid areas to the combined area of discs spanning each individual circular frame. The latter is equal to the total area of the Goldschmidt solution [29]: two discs spanning the two circular frames, with an infinitesimal thread connecting the two. We see in Fig. 8(c) that the Goldschmidt solution has a smaller area than the stable catenoid for D > D * , where D * = 0.528. For spacing D < D c the two branches of solutions for α have the structure of a saddle node bifurcation. We define the dimensionless boundary deformation via with < 0 in this subcritical case and > 0 in the supercritical case discussed below. A Taylor expansion around the critical point of Eq. (16) yields where c = √ 2α c /D c , and we have used the condition (∂/∂α)C(α, D)| α c ,D c = 0, which holds at the critical point. Figure 8(b) shows that this is an accurate approximation to the two branches, even for = O(1).

C. Stretched catenoids
In the supercritical regime, there is no minimal surface solution of the form of Eq. (14), and any surface spanning the two frames is unstable. To make analytical progress, it is thus necessary to define a suitable initial condition and we adopt the stretched catenoid introduced previously [14], which satisfies the boundary conditions ζ (±d ) = R for any value of a, and reduces to the minimal solution when a satisfies Eq. (16) for d < d c . Since beyond the critical spacing there are no minimal solutions, we may choose a shape with any value of a as an initial condition. Close to the critical catenoid, it is natural to set a = a c , and thus, the same parameter used in the subcritical case, = (D − D c )/D c , now with > 0 can be used to parametrize the instability. This leads to the continuation of the pitchfork bifurcation in Fig. 8(b) along the path in the D − α plane shown by the dotted line. The area of the stretched film can be calculated directly from Eq. (20) [14], where we use the notation th u ≡ tanhu, sh ≡ sinh, etc. and x c = D c /α c = 1.1997 . . ., and extends smoothly beyond A c as shown in Fig. 8(c).

D. Dynamics of perturbations within MCF
A natural first step in analyzing the dynamical evolution of a catenoid towards collapse is to study the problem within the simplest model, namely, mean curvature flow (MCF). This is an overdamped temporal evolution, in which the surface area plays the role of a Lyapunov function, dissipation is purely local, and there is no global conservation law for the volume enclosed by the film. As such, it is incapable of representing the dynamics of real macroscopic soap films, whose evolution is more nearly Hamiltonian and constrained by the incompressibility of the surrounding air. Yet, the simplicity of MCF allows for a clear understanding of the perturbative calculations necessary for linear and weakly nonlinear problems, and the comparison between MCF and Euler flows is very instructive regarding the mode structure. To make this comparison we first present in this section the linear and nonlinear dynamics of perturbed catenoids within MCF. In the analysis of stability and dynamics near a catenoid there are two relevant classes of shape perturbations: radial and normal. The former is most appropriate when considering an Eulerian PDE for the surface evolution, as in MCF, while the latter is more suited to comparisons with formal results from differential geometry that utilize intrinsic coordinate systems. In the following we focus on radial perturbations and comment later on the connection between the two.
MCF is a time evolution in which local viscous drag balances the force per unit area arising from surface curvature, i.e., μU = γ H, where μ is a drag coefficient and γ is the surface tension. We use the parameter a to define the dimensionless time t = γ t/2a 2 μ and position r = ρ/a, and drop the prime to obtain the scaled MCF n · r t = 2aH.
Rescaling the coordinate and radial displacement as MCF in an Eulerian reference frame is given by [20] It follows that the rescaled area S = A/a 2 monotonically decreases in time, We note parenthetically that a dynamics akin to Eq. (24), in which the Laplace pressure balances the inertia of the film itself (rather than the surrounding air), has often been used to describe dynamics of supported catenoidal films [14,21,22]. If ρ f is the film density and is its thickness, then the equation of motion in a rescaled time variable t = γ /2ρ f a 2 t is (dropping the prime) The stability analyses for subcritical and supercritical cases differ primarily in the nature of the base state, as the minimal surfaces in the former regime do not exist in the latter. So for the sake of generality at the moment we keep the base state f arbitrary and define the dimensionless film profile as (27) where with this sign convention g is positive for an inward displacement of the surface. Then H can be expanded in a series involving operators of increasing order of nonlinearity in g as and m = 1 + f 2 u .

E. Subcritical catenoids
In the subcritical and critical regimes, we set f = cosh(u) and obtain H 0 = 0 for any ring spacing and a Sturm-Liouville operator with weight function unity. The associated eigenfunctions of L vanish on the boundary of the domain, so L is self-adjoint and its eigenfunctions V n (u) for n = 0, 1, 2, . . . form a complete orthogonal basis, and its eigenvalues σ n are real and distinct. Thus, within MCF, we may write Since r = f − g, the equation of motion for g is and thus ḃ n V n = −H. If we now restrict attention to the leading term H 1 , then we findḃ n = σ n b n , where {(V n , σ n )} satisfy the eigenvalue equation where V n vanishes on the boundary (u = ±D/α). For the critical catenoid, the well-known zero mode is where N is a normalization factor set by the condition so N = 0.8482 . . .. Before proceeding further, we comment on the alternative approach to perturbations, one based on displacements normal to the catenoid. This formulation, unlike our analysis thus far, is most appropriate to Lagrangian descriptions of the dynamics. In this approach [14], when a normal perturbation to a surface of revolution ρ(z) is applied, the displaced surface is z are the displaced vertical and radial coordinates. Taylor expanding to find ζ (z) yields the leading order result in our rescaled units g = w 1 + f 2 u , where w = ω/a and, for the catenoid base state, g = w cosh u. This implies that the neutral normal mode is W 0 = N (1 − u tanhu). Substituting this into the equation of motion g t = −H 1 yields an eigenvalue problem closely related to Eq. (33), The operatorL can be written as a 2 (∇ 2 − 2K ), which is the familiar Jacobi operator [9], where K = −1/a 2 cosh 4 u is the Gaussian curvature of the catenoid. Thus, if V n is an eigenfunction of L it follows that W n = V n cosh u is an eigenfunction ofL with the same eigenvalue σ n . A similar relationship involving radial and normal perturbations is found by examining the second order contribution A 2 to the total surface area [37], which can be written in several forms: where dS = cosh 2 u du in Eq. (37a) is the element of surface area. The result in Eq. (85) is the well-known quadratic contribution for normal perturbations [9,38]. The calculation of higher-order contributions to the area functional requires a careful analysis of the connection between radial, tangential, and normal perturbations [39][40][41][42][43] Returning to the stability analysis, Fig. 9(a) shows the first five modes V n for the critical catenoid, ordered by the number of nodes, obtained by an eigenvalue solver (bvp4c) in Matlab. For comparison, Fig. 9(b) shows the modes for a subcritical catenoid ( = −0.05), which appear as compressed versions of the critical modes. The growth rates σ n for these two cases are shown in Fig. 9(c). In the critical case, all modes are 035105-8 linearly stable (with σ n < 0), except for the zero mode, which is neutral (σ 0 = 0). In the subcritical case, all modes have negative eigenvalues and hence are linearly stable.
To understand the dependence of the growth rates σ on we make use of the fact that the subcritical modes differ only quantitatively from the critical ones, which speaks to the suitability of a regular perturbation theory. The key ingredient in this approach is the fact that the parameter a + used to define the rescaled variable u = z/a + appearing in the eigenvalue problem Eq. (33) admits an expansion in powers of √ − near the critical point [cf. Eq. (19)], where we define x = z/a c . Expanding L, σ and V within this same framework, we have Expanding L, we find L (0) (x) = L as in Eq. (30a) and Collecting terms, we find (L (0) − σ (0) n )V (0) 0 = 0 and, at O( √ − ), an inhomogeneous equation for the perturbed eigenfunction and eigenvalue, The boundary condition V n (±D/α) = 0 can likewise be expanded near the critical point using The eigenvalue correction σ (1) n is then found by orthogonalization with the null space of the operator on the left-hand side of Eq. (41). At leading order, we use the functions V (0) n and obtain

F. Nonlinear instability of the critical catenoid
We next consider the first situation in which the catenoid collapses as is increased from negative values: critical catenoids ( = 0) under finite amplitude perturbations. As noted by Chen and Steen [2] in the context of inertial dynamics, perturbation to the critical catenoid in the form of arbitrary initial values of the relevant modes produce oscillatory motion superimposed on the collapse process. However, for the special case of an initial condition involving only the neutral mode, the evolution is monotonic. The robustness of the pinching phenomenon to the initial absence of higher modes allows us to focus purely on neutral perturbations of where {p (n) k } are projection coefficients. Note that the functions H n in Eq. (28) are defined for an arbitrary perturbation g which need not satisfy the boundary conditions g(±d ) = 0 that must hold for a physical system, whereas the modes V k do vanish there. Thus, the mode decomposition exhibits a kind of Gibbs phenomenon at the domain boundaries. Table I gives the first four nonzero projection coefficients at second and third order.
Within this approximation, the dominant mode is V 0 and the nonlinear evolution of its amplitude has the forṁ where the negative signs arise from the particular choice Eq. (27) for the sign of g. Defining B = b 0 (0)/b 0 (t ) and ν = p (3) 0 /p (2) 0 , Eq. (44) can be solved implicitly as At very early times, we may neglect the contribution from the cubic operator and obtain the explicit solution Unlike the exponential growth exhibited by a linear instability, Eq. (44) shows algebraic growth. Since p (2) 0 < 0, Eq. (46) shows that if b 0 (0) > 0, then the amplitude grows in time and the system has a finite-time singularity, when the catenoid waist vanishes, at a time t p such that b 0 (t p )V 0 (0) = 1, or As p (3) 0 < 0, inclusion of the cubic term shortens the singularity time.

G. Supercritical catenoids
For catenoids with D > D c we substitute Eq. (20) into the governing MCF Eq. (24) and define and m β = β 2 + sinh 2 x. Then for 0 we have where again x = z/a c , and the operator L β is the second relation being in the standard Sturm-Liouville form. The operator L β reduces to L as β → 1 ( → 0) and m β → cosh 2 x. From Eq. (49a) the stretched catenoid Eq. (20) has nonzero mean curvature H 0 when β = 1.
As with the scheme used for the subcritical case, we expand Eqs. (49a) and (49b) to first order in , using β 1 + + · · · and m β cosh 2 x + 2 + · · · , to obtain These are expressed in terms of the scaled coordinate x ∈ [−x c (1 + ), x c (1 + )], and to obtain the full expansion in we introduce the rescaled coordinate y = x/(1 + ) to map the problem back to the domain of the critical catenoid. Further expanding in , the leading-order results are L β = [1 + 2 (th 2 y − ythy)]L (0) (y) and the linear dynamics of the perturbation is The homogeneous Sturm-Liouville problem L β V n = σ n V n can be treated perturbatively in a manner analogous to that for subcritical catenoids in Eqs. (39), except that the expansion is now linear in . Figure 10 shows the eigenfunctions, eigenvalues and eigenvalue corrections for the first modes, illustrating the accuracy of the linear shift in comparison to the exact calculations on the stretched domain. All eigenvalues σ n acquire positive contributions linear in , such that for moderate stretches of, say, = 0.1-0.2, the mode that is neutral at = 0 acquires a positive growth rate, while all other modes remain strongly damped. This is the typical situation in a pattern-forming system in a finite domain [44], and we would expect the strongly damped modes to be slaved to the weakly unstable one, whose evolution we first calculate. If we set g(x, t ) = b 0 (t )V 0 (x) in Eq. (53)  catenoid is initially stretched but otherwise not perturbed, then b 0 (0) = 0 and the unstable mode grows on the slow timescale σ (1) 0 t due to the initial nonzero mean curvature, with The amplitudes of the strongly damped higher modes n > 0 are obtained at quadratic order by considering the separate contributions from the terms in the equation of motion g t = −H 0 + L β g − H 2 , where H 2 is given in Eq. (29b). If g = b n V n , then the left-hand side isḃ n V n , while those on the right are, respectively, of order , σ (0) n b n V n , and b 2 0 p (2) n V n , where p (2) n is the projection on V n of the operator H 2 , evaluated with g = V 0 , as in Eq. (43). Ignoring the O( ) contribution, the slaving approximationḃ n = 0 implies Since the growth rates σ (0) n grow in magnitude roughly quadratically in n (cf. Fig. 10), in this slaving approximation the dominant even mode will be n = 2, with b 2 > 0 at this level of approximation.

H. Euler flow
In our experiments, with catenoids supported on rings of radius R = 4 cm, the typical timescale for collapse is ∼0.1 s, giving a characteristic speed U ∼ 40 cm/s. Using the kinematic viscosity ν a of air, the Reynolds number Re associated with the expulsion of air from the interior is ∼10 3 , justifying a treatment of the air motion as an Euler flow, particularly at the early stages of collapse. The simplest model for the motion of the soap film itself is a vortex sheet [45] coupled to that Euler flow. This dynamics, unlike mean curvature flow, is intrinsically nonlocal. While this problem has historically been studied from a Lagrangian perspective [45], a recent Eulerian reformulation [46], motivated by prior work on interface motion in lubrication theory [47][48][49][50], is particularly suitable to the present analysis. In this approach, the motion of an axisymmetric soap film is given by the coupled dynamics of the radius ζ and vortex sheet strengthη, where ρ a is the density of air, F is the flux associated with the conservation law for cross-sectional area and w is the longitudinal fluid velocity, both of which can be expressed as nonlocal integrals involving ζ andη [46]. As before, H is the mean curvature. The flux form of Eq. (56a), which embodies the conservation law associated with the air surrounding the film, distinguishes this dynamics from MCF. Using the general relationship Eq. (2) between ζ t and the normal velocity, we have and thus the simple result (assuming up-down symmetry) which illustrates that F is the flux of air past a given point. For long slender necks, F can be approximated by the local expression F =ηζ 2 , reducing the dynamics to the well-known one-dimensional model [51].
To analyze the linear instability of a stretched catenoid under this local approximation we note that to leading order the longitudinal velocity vanishes (motion starts from rest). With the rescalings of ζ and z used previously and the rescaled time t = t/τ cap , where is the capillary time, the rescaled vortex strength is η = η(2a c ρ a /γ ) −1/2 , and the linearized dynamics of the perturbation g away from the stretched catenoid f is where, up to linear order,Ĥ 0,1 = m −1/2 H 0,1 as per Eqs. (28), (49a), and (49b). This pair of first-order equations can be converted into a single second-order one by differentiating Eq. (60a) with respect to time and using Eq. (60b) to obtain whereL β = m −1/2 L β . Compared to the mode dynamics Eq. (53) within MCF, we see not only inertial aspects (through the second time derivative) but also the effect of the conservation law (through additional spatial derivatives). An example of the applicability of the local approximation is the Rayleigh-Plateau instability of a cylinder under the Euler dynamics Eq. (61), where one can compare with the exact result incorporating nonlocality [46]. With a cylindrical base function of radius R and with time rescaling t = t/(ρ a R 3 /γ ) 1/2 , the linearized dynamics is where x = z/R here. Resolved into Fourier modes of rescaled wave vector q, the growth rate is σ (q) = q 2 (1 − q 2 ) versus the exact result q 2 (1 − q 2 )I 1 (q)K 1 (q) given in terms of modified Bessel functions. The approximate result is qualitatively identical to the exact one, while being quantitatively larger by a factor ∼1.6. We first consider the nonlinear dynamics of the critical catenoid under this Euler dynamics, analogous to the MCF version of the problem described in Sec. III F. As remarked in Sec. III F, we avoid the possibility of oscillatory excitations by considering a perturbation to the critical catenoid involving only the marginal mode, and set g(t ) = b 0 (t )V 0 (u). As with MCF, since the critical catenoid has zero mean curvature and L β V 0 = 0 when β = 1, the leading order equation of motion is wherê This yields the equation of motion for b 0 , where q (2) 0 is the projection of the right-hand side of Eq. (64) onto V 0 . Equation (65) has the solution where v = (2q (2) 15 • , and we have assumed zero initial velocity for b(t ). As noted previously [2], the early time behavior is quadratic in t, FIG. 11. The linear stability problem for Euler flow for critical and supercritical catenoids. (a, b) First five eigenfunctions W n (u) of the linear operator in Eq. (61) for a the critical catenoid ( = 0) and a weakly supercritical one ( = 0.1). (c) Growth rates σ n and oscillation frequencies ω n for the first seven modes for the cases in panels (a) and (b). Note the frequency scale (ω n ) increases downward, while the growth rate scale (σ n ) increases upward and is magnified by a factor of ten. b 0 (t )/b 0 (0) = 1 + (q (2) 0 /2)t 2 + · · · , in contrast to the linear behavior Eq. (46) within MCF.
Finally, we return to the case of the stretched catenoid under Euler flow Eq. (61). A calculation completely analogous to that within MCF for the stretched catenoid leads to the Euler-dynamics modes for the critical and stretched catenoids and their growth rates, as shown in Fig. 11. As in the case of MCF, the modes of the stretched system appear simply as stretched versions of those for the critical case. Under this inertial dynamics the modes of the critical catenoid for n 1 are oscillatory, with real frequencies shown in Fig. 11(c), while the zero mode is neutral. For > 0 the oscillatory modes acquire small shifts in frequency but do not otherwise change If we again assume an unperturbed stretched catenoid, then b 0 (0) = 0, and if furthermore it is not initially moving, sȯ b 0 (0) = 0, then This result provides the sought-after form of the exponential growth of the mode amplitude.

IV. MODE DECOMPOSITION OF EXPERIMENTAL COLLAPSE DYNAMICS
The results of Secs. III E-III H, as summarized in Figs. 9-11, reveal that the form of the modes associated with the linear stability of catenoids is essentially independent of the underlying equations of motion and that it is instead the growth rates which depend on the dynamics (MCF or Euler). This quasiuniversality suggests that a mode decomposition of the experimental shapes will give model-independent insight into the actual collapse process. The dynamics of the mode growth is also a check on the validity of the Euler calculation.
Such a projection proceeds along the lines by which we defined the base and perturbed states of catenoids in Eq. (27), namely, from the physical interface shape ζ (z, t ) we define the perturbation function g(u) = f (u) − ζ /a, where f (u) = cosh(u) is the critical catenoid, with u = z/a ∈ [−d/a, d/a] and a = α c R is given by the measured hoop radius R and the critical value α c = 0.5524 . . .. As f is defined to match the hoop radius at z = d, g vanishes identically at the domain boundaries. Using, for example, the critical modes V n (u) that satisfy Eq. (33) and are shown in Fig. 9, we write and use orthonormality to obtain the b n (t ) by direct numerical integration using the V n . By the experimental z → −z symmetry, b n = 0 for odd n.
The results of this projection are shown in Fig. 12 through the function g itself and the first four nonzero projection coefficients. The central result of this analysis is the very strong dominance of the shape evolution by the mode V 0 that is neutral for the critical catenoid. Figure 12(d) shows that its amplitude grows exponentially until very close to the critical time, when the higher modes begin a rapid growth while 035105-13 remaining relatively small even at t c . The subdominant mode until just near t c is V 2 , but it is rapidly overtaken by V 4 whose amplitude b 4 (t ) < 0. With V 4 as in Fig. 9(a), we infer that splitting of the minimum of the film shape into two incipient pinch points arises from the negative value of b 4 .
A fit of the linear portion of the semilogarithmic plot in Fig. 12(d) to the form b 0 (t ) = A exp(σ exp t ) yields the experimental growth rate σ exp = 28.4 ± 0.1 s −1 . From the original rescaling of the equations of motion introduced above Eqs. (60) and the perturbative result in Eq. (68) we expect the dimensional growth rate to be where the factor of arises from the stretch of the catenoid. Using experimental parameters a c = 0.55...R, with R = 4 cm, ρ a = 10 −3 g/cm 3 , and γ = 35 dyn/cm, we have τ −1 cap ≈ 40 s −1 , and with q (1) 0 2.83 in Eq. (67) by direct numerical integration, we find σ dim 67 √ s −1 . With in the range of 0.1-0.2 for our studies as described in Sec. II, we therefore find good agreement between the Euler calculation and experiment.

A. Cone angle at pinching: Approach based on a mathematical sequence
Here we use a simple geometric approach to calculate θ * , the cone angle at pinching introduced earlier. It is based on approximating the observed shape of the evolving soap film surface by a minimizing sequence of shapes [52]. This method, which is suitable for area minimization when the solutions are continuous, differentiable functions has been extended to discontinuous solutions of functionals of the form [53] where ζ (z) and ζ z are the radial coordinate and its z-derivative, respectively, which clearly includes the area functional for a surface of revolution. That is, the minimizing sequence method can be applied to study the approach toward the Goldschmidt solution. Following Koshelev and Morozov [53], ifζ is the minimizing function of the functional J[ζ ], then any sequence {ζ ε } (where ε is a label and does not indicate differentiation) for which the limits exist, is a minimizing sequence. Because the solution need not be differentiable, it is possible to construct a sequence consisting of right circular cones and cylinders with discontinuities in derivatives.
To proceed, we note that it is possible to construct from truncated piecewise conical surfaces, which are clearly not minimal, a suitable approximation to any catenoid [54,55]. In fact, the area 2S (c) of the critical catenoid, i.e., the last stable catenoid bounded by two circles of radius R separated by a distance 2d, is equal to that of the surface constructed with two right circular cones of height d and base radius R, each of whose vertices touch, and where θ 0 is the half apex angle of the cone. S (c) is the area of the first element of a sequence of surfaces consisting of just two separate cones, whose vertices are connected by a line, and whose half apex angle becomes successively larger, i.e., θ 0 < θ 1 < θ 2 < · · · π/2, each with half-area we denote by S (c) n , which decreases as n increases. While this sequence satisfies the conditions Eq. (72), the connecting line is unphysical.
We now construct a more physical sequence such that (i) it tends to the Goldschmidt solution, (ii) the area S n of its nth element is greater than S (c) n , and (iii) it eliminates the unphysical connecting line. It consists of shapes composed of two oppositely oriented truncated cones with apex angle θ n , joined by a cylinder of radius h(θ n ): the "Martini glass" shape ( Fig. 13). The elements of the new sequence then have area where the half-length of the connecting cylinder is It is possible to construct a sequence ζ ε (z) which is an upper bound to the area of the successive surfaces approaching the Goldschmidt solution. Introducing the dimensionless λ ε = ζ /R, λ = h/R, δ = d/R, and s = z/R, and setting ε = 1/ tan θ n , we have where λ(ε) is as yet unspecified. There is an infinite number of ways to choose λ(ε). One of the simplest is to let h(θ n ) 035105-14 be the radius of the cylinder joining two annuli separated by a distance 2d with inner radius R − h(θ n ) and outer radius R, such that the combined area of the three elements (annuli and cylinder) is 2S (c) n , or S (c) n = π [R 2 − h(θ n ) 2 ] + 2π h(θ n )d. In dimensionless units, where we have chosen the root of the quadratic equation that has λ → 0 as ε → 0 (the minus sign before the outer radical); that is, the solution tends (in dimensional units) to the Goldschmidt solution To verify that this sequence is a suitable one, it is only left to show that lim ε→0 J[ζ ε (z)] = J[ζ G (z)] [53]. Direct integration yields which becomes 2π R 2 , as ε → 0, as required.
A test of the suitability of the sequence to describe the physical system, when λ(ε) is given by Eq. (78), is its prediction of the cone angle at pinching. From the Rayleigh criterion that a cylinder is unstable under the action of surface tension when its length exceeds its circumference, the cylindrical section of the Martini glass shape will be unstable when b = π h. With S n = S n /(π R 2 ), Eq. (75) yields As n grows and S n decreases from the value associated with the critical catenoid, the first nonnegative solution of Eq. (81) for ε appears when S n 1.1244, which corresponds to θ * 69.5 • , equivalent to 2θ * 139 • , within 3% of the experimental value. The entire evolution of experimental shapes, from the initial catenoid to that with the cone angle at pinching, can also be compared to the sequence. To do this, we return to the general expressions Eqs. (75) and (76), and use the pairs (S n , θ n ) obtained from experiment to calculate the radius as the solution of the resulting quadratic equation where, again, it is necessary to choose the root that takes λ → 0 when ε → 0. Implementing this procedure for a typical set of experimental results yields for each pair (S n , θ n ) a point in the ε − λ plane, as shown in Fig. 14. Overlaid on these data is a contour plot of the area function, and arrows indicating the local gradient of that function at each point. After an initial transient period, the evolution becomes progressively more aligned with the direction of steepest descent, indicating that the system evolves as a gradient flow, and providing further evidence of the accuracy of the representation of the shapes as combinations of cylinders and cones. FIG. 14. Mapping experimental data onto the sequence. Contour plot of the energy function S in the parameter space of cone angle θ and scaled radius λ. Arrows on experimental data points indicate −∇S.

B. Cone angle at pinching: Heuristic argument
An alternative approach to calculate the cone angle at pinching involves finding functional relationships between surface area, neck radius, and the tangent angle φ at the point of support. This can be done within perturbation theory around the critical catenoid by noting that MCF and Euler dynamics share the same critical mode V 0 (x) given by Eq.
Likewise, the minimum neck radius is Since the catenoid around which we are perturbing is not only a minimal surface, but also the critical one, the first nonzero correction to the area is third order in s, With f (x) = cosh x and V 0 given in Eq. (34), we obtain With the results Eqs. (83), (84), and (86) we have the three geometric quantities given in terms of κ, and thus instead of considering the time evolution of each, we may relate one to another with κ as a parameter and compare these binary relations with the corresponding experimental results (Fig. 15). The good agreement found lends further support to the idea that much of the geometry of the collapsing catenoid is a consequence of the dominance of the fundamental unstable mode V 0 .

C. Conical similarity solutions
We conclude in this section by discussing the arguments first advanced by Day, Hinch, and Lister [25] to explain the inverted double-cone cone structure that leads to the ultimate topological transformation, seen in Fig. 1(v). If the collapse dynamics obeys a similarity solution so that the velocity potential has the scaling form where P = ρ/(γ τ 2 /ρ a ) 1/3 and τ = t c − t, then the potential function obeys the full Laplace equation in two dimensions (assuming axisymmetry). The kinematic and Bernoulli boundary conditions at the interface are consistent with a conical solution at large P, with ∼ |P| 1/2 , equivalent to φ ∼ ρ 1/2 , where ρ = |ρ|. Regularity of the solution of Laplace's equation then requires φ ∼ ρ 1/2 P 1/2 (cos θ ), where P 1/2 is the Legendre function of order 1/2. Such a velocity potential produces a contribution to Bernoulli's equation |∇φ| 2 proportional to 1/ρ, which can naturally balance the Laplace contribution to the pressure proportional to the mean curvature of the surface, as H ∼ 1/ρ for a cone, provided that curvature is negative, as in the larger of the two cones seen in Fig. 1(v). For the smaller cone, whose curvature is positive, an additional source term proportional to τ /ρ can provide that balance. The functional form of Eq. (88) for the velocity potential is the same invoked by Taylor [56] to explain the conical shapes of water drops in strong electric fields. There, the requirement that the surface of the drop be an equipotential fixes the cone angle θ 0 to be that of the zero of P 1/2 . Unlike that situation, the solution outlined above provides a balance of forces that supports a conical solution, but does not, by itself, determine the angles; they emerge from the full analysis of the similarity solution. In the case of the acute angle cone studied here, we have seen that the conical structure reaches from the pinch point all the way to the supporting rings. The arguments based on the sequence of shapes in Sec. V A and the dominance of the unstable mode in Sec. V B suggest that it is fundamentally the Rayleigh instability of the central cylinder that dictates the pinch point and thus the acute angle of the large cone.
The power-law decay of the minimum neck radius shown in Fig. 5 and the blow-up of the Gaussian curvature there both suggest that the acute-angle cone that we have focused on is describable in terms of a classic similarity solution, but with two important differences. First, the pinching region has two cones with curvatures of the same sign: the large cone we have focused on, reaching the supporting frame, and the smaller one connecting to the central cylindrical region that ultimately  forms the satellite bubble. Second, the pinching process does not go to completion. Instead, after the critical time t c those two cones begin to retract, connected by a very thin thread. Thus, an important issue is the extent to which the pinching interface conforms not only to the Keller-Miksis scaling but the more general behavior of self-similarity. If self-similar, then we should expect the data to collapse when we adopt instead of (u, r) the similarity variables (χ, ξ ) with χ = u − u n r n and ξ = r r n .
(89) Figure 16 shows such a plot, with the data used in previous figures. While at a coarse scale there appears to be an evolution towards a fixed scaling function, the magnified view in Fig. 16(b) shows the failure to converge to a true timeindependent similarity form. The above suggests a picture in which the acute-angle cone formation is associated with an unstable similarity solution, as depicted in Fig. 17 in the space of the mean curvatures H 1 and H 2 of two surfaces that meet at an incipient pinch point. The initial minimal surface at H 1 = H 2 = 0 flows under the collapse dynamics toward the acute-angle fixed point along an attracting direction, only to be repelled before the true singularity, flowing then toward the stable [26] obtuse-angle conical fixed point and the ultimate topological rearrangement.

VI. CONCLUSIONS
In this paper, we have focused primarily on the collapse of three-dimensional catenoidal soap films induced by deformations of the boundary, and have shown by several methods how one can understand global geometric features of the films near pinching. Analysis of the linear and weakly nonlinear dynamics of both stable and unstable catenoids shows that the underlying modes are qualitatively the same whether the dynamics is mean curvature flow or Euler flow, even though the dynamics of the modes in the two cases are fundamentally distinct. Our study of the mode dynamics just beyond the critical catenoid reveals a familiar structure in the field of pattern formation, with a fundamental unstable mode and an infinite ladder of stable ones. This structure raises the question of whether it may be possible develop an "inertial manifold" dynamics for the unstable mode in a manner analogous to that which has been described for interface pinching singularities in Hele-Shaw flow [48,49]. The unified analysis of the role of boundary deformations on a paradigmatic problem that we have presented in this paper constitutes a framework to study more complex problems. Among them is the collapse of hemicatenoids bounded by solid surfaces, with a free contact line [57,58] where, as we discuss elsewhere [59], the motion is dominated by viscous drag in the moving contact line, and more significantly the collapse dynamics of soap films with exotic topologies [3][4][5][6][7].

ACKNOWLEDGMENTS
We thank Gareth Alexander, Markus Deserno, Keith Moffatt, Thomas Powers, and Michael Shelley for discussions, and we are grateful to Paul H. Steen for discussions prior to his untimely death. We dedicate this paper to his memory. This work was supported in part by Established Career

APPENDIX: EVOLUTION OF GEOMETRIC QUANTITIES
In this Appendix we provide a self-contained summary of the key results for the time evolution of geometric quantities relevant to the calculations in the main body of the paper. We start by assuming that we have a surface traced out by the vector x(u α , t ) through internal coordinates u α , α = 1, 2, and varying with time t. The surface has a metric tensor with components where here subscripts denote differentiation, so x α ≡ ∂x/∂u α . The x α are the unnormalized tangent vectors to the surface. They define the unit normaln as, say, or its opposite if we switch 1 ↔ 2. Here, g = det(g) is the determinant of the metric tensor, where is the Levi-Civita symbol and the summation convention holds. The inverse metric tensor, with coefficients g αβ is the usual The coefficients h αβ of the second fundamental form are defined by the normal and tangent vectors as If we denote by Tr the conventional trace of a matrix, then the mean and Gaussian curvature, H and K, of the surface can be expressed in terms of h as H = 1 2 Tr h α β , (A6) where h α β = g αμ h μβ .

Dynamics
For the dynamical problem, the surface is assumed to move purely in the normal direction, with a velocity U that is some local or nonlocal function of the surface shape, Given this law of motion, we calculate the time evolution of g, h, and hence H and K. The fundamental assumption we make is that the internal coordinates u α are independent of time.
Let us start with the normal vector. Using the relationship n · x α = 0, we differentiate with respect to time to obtainn t · x α = −n · x αt . Now, x αt = ∂ α x t = ∂ α (Un) and since n ·n α = 0, we haven t · x α = −∂ α U . To obtain n t itself, we note that it can be expanded in terms of the tangent vectors asn t = A μxμ , for some coefficients A μ , yielding A μ g μα = −∂ α U . Multiplying by x α so, recognizing the middle expression asn t again, we havê Next, we examine the components of the metric tensor. Differentiating directly we find ∂ t g αβ = x αt · x β + x α · x βt . (A11) Again noting thatn · x α = 0, we obtain ∂ t g αβ = [n α · x β + x α ·n β ]U = −2h αβ U.
From this, we can obtain the time evolution of det g, Comparing with Eq. (A6), we deduce Note the sign difference compared to Eq. (18) of Ref. [36]. We can do two checks on this result. The first is to note that the sign of the mean curvature depends on the particular parametrization chosen. For example, the normal vector of a sphere may point inward or outward depending on the order in which the two coordinates are chosen, but the mean curvature will also have opposite signs in the two cases. If we flip the normal, then H changes sign but so too will U , leaving ∂ t g unchanged.
In the particular case of mean curvature flow (MCF) of a sphere, we can carry through the entire calculation easily.
Thus, ζ ζ t = −1, which yields ζ 2 = ζ 2 0 − 2t. Choosing the origin of time to when ζ vanishes, we find This is the famous self-shrinker of MCF [60], apart from a factor of 2 that differs from that found in the mathematics literature due to the definition of H.
Returning to the derivation, we now examine the coefficients h αβ of the second fundamental form. Acting directly on the definition Eq. (A5), we obtain and one can show that n α ·n β = 2Hh αβ − Kg αβ .
For the Gaussian curvature evolution we obtain This is consistent with prior results [36], apart from some differences in notation for Christoffel symbols.

Surface of revolution
For an axisymmetric surface described by 035105-18