A compact Eulerian representation of axisymmetric inviscid vortex sheet dynamics

A classical problem in fluid mechanics is the motion of an axisymmetric vortex sheet evolving under the action of surface tension, surrounded by an inviscid fluid. Lagrangian descriptions of these dynamics are well-known, involving complex nonlocal expressions for the radial and longitudinal velocities in terms of elliptic integrals. Here we use these prior results to arrive at a remarkably compact and exact Eulerian evolution equation for the sheet radius $r(z,t)$ in an explicit flux form associated with the conservation of enclosed volume. The flux appears as an integral involving the pairwise mutual induction formula for vortex loop pairs first derived by Helmholtz and Maxwell. We show how the well-known linear stability results for cylindrical vortex sheets in the presence of surface tension and streaming flows [A.M. Sterling and C.A. Sleicher, $J.~Fluid~Mech.$ ${\bf 68}$, 477 (1975)] can be obtained directly from this formulation. Furthermore, the inviscid limit of the empirical model of Eggers and Dupont [$J.~Fluid~Mech.$ $\textbf{262}$ 205 (1994); $SIAM~J.~Appl.~Math.$ ${\bf 60}$, 1997 (2000)], which has served as the basis for understanding singularity formation in droplet pinchoff, is derived within the present formalism as the leading order term in an asymptotic analysis for long slender axisymmetric vortex sheets, and should provide the starting point for a rigorous analysis of singularity formation.

pulled apart beyond a critical separation, rendering the film unstable. Extensive experimentation [36] has shown that the collapsing axisymmetric surface eventually breaks up through self-contact at multiple points, producing a series of satellite soap bubbles (Figure 1.1). Theoretical work on this dynamical process ranges from Maxwell's original stability analysis [26] to much more recent computational studies of fluid dynamical models [6,9] focused on the nature of the singularities.
In a broader sense, there have been two schools of thought in the study of singularity formation by moving surfaces. On the one hand, there is a very substantial body of rigorous work on the simplest geometrical law, namely motion by mean curvature [7]. However, while this law is a physically realistic description of interface motion arising from surface diffusion, it cannot be applied to the motion of soap films because it cannot account for the dynamics of the surrounding fluid and the conservation laws that follow. There are modifications of mean curvature flow that conserve volume enclosed by the surface, but they are not faithful representations of the dynamics of the surrounding fluid. On the other hand, there is the large body of more phenomenological work on singularity formation in fluid mechanics [13], where simplified PDE models have been developed to address self-similar dynamics near singularities. These models are physically realistic but have generally lacked rigor.
The two examples where the gap between these two approaches have been bridged is in the context of interface motion in two dimensions, where a systematic procedure to derive PDEs for the evolution of asymptotically thin fluid layers has been developed from exact boundary-integral formulations [18][19][20]34]. These analyses have put on a solid foundation empirical models [8,11,15] derived within lubrication or long-wave theory. Given these results, the question then arises of whether there is a comparable physical setup in three dimensions in which a PDE can be derived for the motion of a surface surrounded by an incompressible fluid. The simplest example of this would naturally be an axisymmetric surface surrounded by an inviscid fluid. If there is an affirmative answer to this question, such a PDE should have an explicit flux form.
In the inviscid limit, the moving interface can be represented by a vortex sheet with surface tension. The predominant approach to this problem has been through a Lagrangian formulation [3,10], which is the most appropriate for computational studies [2,24,29,30]. However, such an approach does not readily lend itself to the development of asymptotic models appropriate to thin necks, as would be relevant in the neighborhood of singularities. In contrast, an Eulerian formulation would not only be amenable to asymptotic analysis, but would also be subject more easily to rigorous studies.
Here, as a first step, and starting from the Lagrangian formulation, we derive an exact Eulerian dynamics for axisymmetric vortex sheets with surface tension, and show that it has an explicit flux form. Naturally, because the problem is deeply nonlocal due to the Biot-Savart interactions between distant elements of the vortex sheet, the flux is an integral over the entire sheet whose kernel is precisely the FIGURE 1.1. Collapse of a catenoid. Snapshots tens of milliseconds apart after the distance between the two supporting rings is increased just beyond the stability threshold. Imaging methods are described in [17]. mutual induction between two coaxial loops as derived by Helmholtz for fluids [22] and by Maxwell for electrical currents [27]. As the mutual induction is, by definition, the flux through one loop due to the circulation in the other, this is a very intuitive result. We show that known stability results for cylindrical vortex sheets can be recovered by direct calculation. While nonlocal, the mutual induction between two coaxial loops of the same radius is sharply peaked as the distance between the loops vanishes, and this feature suggests a natural asymptotic analysis to reduce systematically the dynamics to a local PDE. In this limit we recover at first order the inviscid version [14] of the empirical model first proposed by Eggers and Dupont [15] for this problem.

Flux Equation
We consider an infinite three-dimensional inviscid fluid of density within which is a vortex sheet that is axisymmetric about the´-axis, and whose timedependent radius is R.´; t /. In this inviscid limit there is in general a discontinuity in the tangential fluid velocity across the sheet, and this jump defines the "true vortex sheet strength" x [23,33]. In the same way that we use a minimal surface as an idealization of a static soap film, here the vortex sheet is the idealized representation of a moving film endowed with surface tension . As is well-known [3], the problem of the self-induced motion of the sheet can be expressed, in the so-called Lagrangian frame labeled by˛, as the coupled dynamics of the sheet coordinates r.˛; t / and´.˛; t /, and its "unnormalized vortex sheet strength" .˛; t /, which evolves as where Ä is the mean curvature of the surface, where˛0 is a Lagrangian label, r 0 Á r.˛0; t /, 0 Á .˛0; t /, and Performing the integrations in Â 0 yields the following expressions for the velocities in terms of K and E, the complete elliptic integrals of the first and second kind, respectively, with argument k 2 D 4rr 0 =OE.´ ´0/ 2 C .r C r 0 / 2 : To obtain a Eulerian equation of motion from these Lagrangian velocities we make the change of variables from˛to´D´.˛; t /. Setting R.´.˛; t /; t / D r.˛; t / yields [34] @R @t D r t ´t @R @´; In making the change of variables in equations (2.5a) and (2.5b), it is natural to consider instead the quantity Á D =´˛, so that d˛0 0 D d´0 Á 0 .
For this system, conservation of fluid volume should be expressible as a flux form involving the cross-sectional area R 2 , that is,

@F @´
for some function F , which we now seek. The fact that 2R@R=@t D @R 2 =@t implies that 2R.r t ´t @R=@´/ D @F=@´. However, to infer F by direct substitution of (2.5a) and (2.5b) into this expression is cumbersome and nontrivial. A much better way to find F is to make use of equations (2.3) and (2.4). After some algebra and integrations by parts, we find The integral over the variable Â 0 can be calculated in terms of elliptic integrals, with the result where M.k/ is Maxwell's function first derived for the mutual induction of a pair of circular current loops of radii r and r 0 at locations´and´0, That this remarkably compact result has not been previously obtained in this context may be a consequence of the natural emphasis on the Lagrangian formulation and its computational applications. The results in (2.9) and (2.10) are intuitive in that the mutual induction M is, by definition, the flux passing through one loop due to another, and so F is simply the sum of all the individual fluxes through the loop at position´. Moreover, in the case of a vortex loop, the integrand in F was identified by Helmholtz [22] as the stream function of the flow. Now, if we rewrite the pair (2.5) as´t D w and r t D v, then the evolution of the Eulerian sheet strength Á, using equation (2.6b), takes the flux form with the curvature now re-expressed in terms for R.´; t / as where s´D .1 C R 2 / 1=2 is the arclength metric factor. Equations (2.7), (2.8), and (2.11) comprise the flux-form Eulerian dynamics of an axisymmetric vortex sheet. The Eulerian sheet strength relates to the true vortex sheet strength as Á D s´x . That both and Á evolve through a flux form in their respective frames reflects the fact that their total integrals, in˛or´respectively, give the total conserved circulation of the sheet.
In closing this section, we demonstrate the connection between the functions F and the stream function ‰. Let W .r;´/ and V .r;´/ be the axial and radial fluid velocities, respectively, within the sheet. The evolution of R.´; t / can be expressed as @R @t D U j rDR V j rDR @R @´( cf. equation (2.6a)). Multiplying this equation by R and integrating the incompressibility relationship .1=r/.rV / r C W´D 0, one can easily show Introducing the Stokes stream function ‰, which we identify as the integral on the RHS of (2.13), with rV D @‰=@´and rW D @‰=@r, we obtain (2.14) @ @t thus confirming the connection between the functions F and ‰.

Hamiltonian Structure
In this section we discuss some issues regarding a possible Hamiltonian formulation of the present system for which there is a vortex sheet strength with nontrivial dynamical evolution. It is useful to contrast this case with that of a system of discrete vortex rings and its continuous limit [28], for despite some fundamental distinctions involving conservation laws, there are common mathematical structures involved. As first shown by Dyson [12], the Lagrangian dynamics of a discrete set of coaxial vortex rings with centers Z i on the´-axis, with radii R i and circulation i , can be written as is the mutual inductance between the two loops i and j . The first term on the RHS of (3.1a) is the approximate self-induced velocity of ring i , where a i is the core radius, which serves as a cutoff for the localized induction approximation.
The evolution equation (2.9) can now be seen as the continuum limit of (3.1b), expressed in an Eulerian form. A direct calculation shows a less obvious result, namely that the axial velocity´t in (2.5a) is the equivalent of (3.1a) without the self-induction term. In fact, By analogy to the discrete case, in which one can introduce an energy that is quadratic in the circulations i , namely, and from which one obtains the equations of motion one can take the continuum limit to obtain [31] (3.5) which, generalizing to functional derivatives, yields This is of the same form as the discrete dynamics (3.4), but with the important distinction that the vortex sheet strength itself depends on time (due to surface tension), whereas the individual circulations i in the discrete case do not. Moreover, in the discrete case it is possible to rewrite the dynamics so the left-hand sides are total time derivatives, rendering them truly Hamiltonian. In contrast, the dynamics (3.6) do not obviously have this feature. The fact that (2.1) can be thought of as a nonholonomic constraint may offer a path to obtain an "almost Hamiltonian" dynamics [16].

Stability Analysis
In this section we show how the vortex sheet dynamics in the Eulerian form reproduces known stability results [1,37] for capillary jets, both with and without a streaming velocity within the fluid enclosed by the sheet. Note that our assumption that the fluids inside and outside the sheet have the same density precludes recovering the original stability result of Rayleigh [35], which assumed vacuum outside.
We first consider the case with a quiescent fluid on both sides of a vortex sheet of radius x R and linearize the equations of motion for small perturbations in Á and R of the form (4.1) Á D y Áe i q´Cˇt and R D x R C y e i q´Cˇt : At this order the mean curvature has the simplified form Ä ' R´´ R 1 . The resulting vortex sheet evolution equation is where P D q x R. Since, in the absence of background fluid motion, Á is first order in the perturbation, the remaining factors on the RHS of (2.9) are those corresponding to a cylinder. Let Then, the linearization of (2.9) is After calculating the integral (see the Appendix), we substitute y from (4.4) into (4.2) to obtain, in agreement with previous results [1,37], the growth rate where I n and K n are the n th -order modified Bessel functions of the first and second kind, respectively. Demonstrating that the flux form PDE reproduces the stability results in the presence of streaming flows, obtained by Alterman [1] and by Sterling and Sleicher [37], follows the same procedure as the calculation above, but requires a more delicate analysis. This complexity is related to the limiting procedure of vanishing sheet thickness ı and vanishing viscosity used to arrive at the evolution equation for the vortex sheet strength (2.6). While it is not possible to determine a unique value of the tangential sheet velocity´t , the usual Lagrangian one´t D U=2 [3] would not be compatible with the new Eulerian system of reference. In fact, the Ampère integral is calculated in a frame for which the velocities at infinity are zero, and this is not true for Eulerian frames where the vortex sheet is moving in the presence of a streaming flow; hence the need to "reframe" the results.
To make sure that all the velocities involved in the calculation are measured with respect to the same laboratory frame of reference, we went back to the original boundary conditions when the vortex sheet has a nonzero infinitesimal thickness. From there we deduce what should be the corresponding tangential velocities´t in each one of the auxiliary equations for the vortex sheet evolution so that they would match the boundary conditions corresponding to the external flows on either side of the interface.
For concreteness, consider the situation in which the streaming velocities inside and outside of the cylindrical vortex sheet are U and 0, respectively. In this case the standard Lagrangian choice cannot satisfy the boundary conditions on either side of the sheet in the Eulerian frame. It is therefore necessary to reinstate the appropriate boundary conditions, which would be equivalent to taking the limit ı ! 0 at the end of the calculation.
The calculation proceeds by systematic perturbation of the terms within the flux integral. We expand separately the two, obtaining where k x R is defined in (4.3). Using and collecting terms, we obtain (see the Appendix) (4.9) Substituting into the equations of motion and retaining only first-order terms, we find (4.10)ˇy After some laborious calculation of these nontrivial integrals we obtaiň y D iP I 1 .P /K 1 .P /y Á i Á 0 2 x R P OE1 C P .I 0 .P /K 1 .P / I 1 .P /K 0 .P // y : (4.11) Using the Bessel function identity P .I 0 K 1 C I 1 K 0 / D 1, we further reduce this expression to (4.12)

ÂˇC
i P Á 0 x R K 1 .P /I 0 .P / Ã y D iP I 1 .P /K 1 .P /y Á: This equation is the evolution for perturbations to the mean radius, i.e., the mean location of the vortex sheet: R D .R C C R /=2/, where R C and R are the outer and inner radius of the infinitesimally thin sheet, respectively (R C R D ı ! 0). To reinstate the proper boundary conditions, we note that (4.12) should be split into two parts, the contribution from R C and the one from R . To do so we rewrite (4.12) as In this equation all terms involving I 0 are connected to the inner region, that is R , while those involving I 1 correspond to the outer one (R C ) yielding two equations, one for each region: (4.14) P ÂˇC i P Á 0 x R Ã K 1 .P /I 0 .P / y D iP 2 y ÁK 1 .P /I 1 .P / and (4.15) Pˇy K 0 .P /I 1 .P / D iP 2 y ÁK 1 .P /I 1 .P /: A similar expansion of the vortex sheet evolution equation (2.1) yields for the inner region (4.16) and for the outer one Substituting (4.16) and (4.17) into (4.14) and (4.15), respectively, and adding the contributions to obtain the total expansion for which coincides with previous results [1,37] for the case where the inner fluid velocity is U D Á 0 and the outer fluid is stationary.

Derivation of a Local PDE
To obtain an approximate PDE that describes the dynamics of slender necks, it is necessary to find a suitable approximation to M.k/ that would make it possible to calculate the integral in the flux equation in a controlled manner. With this purpose in mind, it is very useful to note [21] that M.k/ D Q 1=2 . /, with Q 1=2 . / the associated Legendre function of index 1 2 with variable D .2=k 2 / 1, which obeys the differential equation This identity makes it possible to obtain a uniform approximation for M.k/ by matching the inner and outer solutions ( ! 1 and ! 1, respectively). In the inner region the limiting behavior is while in the outer one is for ! 1; where A and B are constants to be adjusted by matching the inner and outer solutions, with the final result in terms of k given by The fit of M.k/ provided by the uniform approximation M u .k/ is remarkable; M.k/ M u .k/ ' 0 only in a small neighborhood of k 2 D 1 2 ; everywhere else it is almost impossible to distinguish M u .k/ from the exact solution (see Figure 5.1).
We can now approximate M.k/ by M u .k/ in the evolution equation (2.9) for R 2 .t /. In order to express M u .k/ in a form that allows calculation, we use the expression for k in terms of´,´0, R, and R 0 . In particular, by using the denominator inside the logarithm in (5.4) is expandable in Taylor series in the limit´!´0 because the ratio .R 0 R/=.´0 ´/ ! r´. In this limit, R ! R 0 , and where x is the same as defined in (4.3). The contribution to the flux F arising from the two terms in the uniform approximation (5.4) involves I 1 and I 2 , where, by the symmetry of the integrand, Performing the integrations we obtain Collecting terms and substituting r D z R in the LHS of the flux equation (2.9), we obtain (5.9) @ z R 2 @t D @ @´.
z R 2 /: While this result was obtained in a systematic way that would be useful in cases where it is necessary or desirable to obtain the next terms of the expansion, the procedure is not very intuitive. In fact, the origin of (5.9) can be more easily understood by recognizing that the first-order approximation corresponds to a cylindrical vortex sheet, where R D R 0 and Á D Á 0 . Then, starting as before from (2.9), we may write where we have used a variant of the usual variable x, namely x D .´0 ´/=2R. Using the relationship M.k/ D Q 1=2 . / and the limit P ! 0 in result (A.11) of the Appendix, we find that the value of the integral is =2. Thus (5.10) yields the same result (5.9) that was obtained for the first-order term by employing the expansion method, which makes clear its geometric interpretation.

Discussion
The present work has shown how the familiar axisymmetric vortex sheet dynamics can be recast exactly in a Eulerian form that clearly displays the inherent flux form of the motion. Previously established linear stability results are easily recovered with this formulation. Most significantly, the Eulerian dynamics is wellsuited to the development of a systematic reduction to a local approximation when the aspect ratio of the system is suitably small. This controlled approximation, which at leading order recovers well-known empirical results [15], should lend itself to more rigorous studies both of the form and validity of the governing PDEs as well as the singularities they produce. As the underlying model includes the dynamics of an incompressible fluid surrounding the vortex sheet, it constitutes a more physically relevant starting point for understanding singularity formation by moving soap films than that provided by mean curvature flows.
In light of the large body of existing work on axisymmetric vortex sheets, there are many possible extensions of the present work, including the introduction of swirl [4] or helical symmetry [3]. Greater challenges would be to incorporate weak viscous effects [5,25] and to complete the Hamiltonization of the inviscid Eulerian case in the presence of surface tension.

Appendix: Details of Stability Analysis
Here we show some intermediate steps in the stability analysis. We begin from the first equation of (4.7), where k R is as defined in (4.3); using the relationships between the elliptic functions and its derivatives, we find The corresponding first-order expansions of the relevant quantities are which yields (4.8), and where we have used the relationships dE=d k D .E K/=k and kdK=d k D E=.1 k 2 / K [21], and have substituted k D k R . Introducing the convenient notation we note that, after straightforward algebra, the quantity M 0 C M 1 , which will be needed later, reduces to k R .E K/.
Utilizing the second of the equations in (4.7) and collecting all the contributions, we obtain Using again the identity M.k R / D Q 1=2 . /, with argument D 2=k R 1, and substituting into (A.6), we obtain the first-order result Substituting (A.7) into the equations of motion and retaining only first-order terms, we find The first integral on the RHS does not contribute to the equation of motion, for 0 D 0 identically when the fluids on both sides of the vortex sheet are stationary, and the integral is a constant when there is a streaming velocity in either one of the two regions. The other integrals involve the Lagrangian parameter˛, and to be able to obtain the result in the Eulerian frame as a function of the coordinate´, we note that ds D p r 2 C´2 d˛and ds D q 1 C r 2 d´. Using these identities allows us to rewrite (A.9) d´D s 1 C r 2 r 2 C´2 d˛' .1 C O. 2 //d˛; making it possible to approximate d˛0 ' d´0 to the required order. Replacing the expressions for the perturbations (4.1) into (A.8), rewriting all exponentials in terms of´0 ´, and expressing them in terms of sines and cosines makes it clear that the only nonzero contributions to the integrals come from the cosine terms, which are the only ones that produce even integrands. Then, defining P D qR and performing the usual change of variables x D .´0 ´/=2R, we find k R .x/ D 1=.1Cx 2 / and .x/ D 1C2x 2 , yielding the first one of the two stability equations (A.10)ˇy The integrals have been calculated with the help of [21,32]. In particular, after a further change of variable, y D p 2x, the first integral acquires its canonical form 7.162.6 [21] with D 1 2 and a D p 2P ; thus (A.11) Z 1 0 Q 1=2 .1 C y 2 / cos. p 2P y/dy D 2 K 1 .P /I 1 .P /: The second integral can be divided in two contributions, To calculate I 1 we use the identity 3.2(2) in [32] for p D 1 and a D x: where N 0 and J 0 are the Bessel functions. With this replacement we obtain I 1 D =2. The second integral can be solved using identity 6.3 (5), also in [32], for a D x and b D 1: This can be reduced to the standard form 6.565.1 (with a D 2P and b D˛) found in [21] by noting that I 2 can also be written as 4P 2 C t 2 dt D 2 P OEI 1 .P /K 0 .P / I 0 .P /K 1 .P /: (A. 16) Finally, after collecting all the contributions and replacing them into (A.10), we find the first-order equatioň y D iP y K 1 .P /I 1 .P / i P 0 y 2R OE1 C P OEK 1 .P /I 0 .P / K 0 .P /I 1 .P /: Making use of standard Bessel function identities, it is possible to further reduce this expression to (A.18) ÂˇC i P 0 R K 1 .P /I 0 .P / Ã y D iPK 1 .P /I 1 .P /y :