\documentstyle[12pt]{article}

\def\Re{\mbox{\rm Re}}
\def\Im{\mbox{\rm Im}}
\def\Ag{\mbox{\rm Ag}}
\def\pist#1#2{\noindent\hangindent 2em\hangafter1\hbox to 2em{#1\hfil~~}#2}
\def\be{\begin{equation}}
\def\ee{\end{equation}}
\def\onehalf{1 \over 2}

\begin{document}
\bibliographystyle{plain}
\title{ A Study of Singularity Formation in the Kelvin--Helmholtz
      Instability with Surface Tension}

\author{Michael Siegel\\ \footnotesize{Mathematics Department, The Ohio State
University, Columbus, OH43210}}

\date{}
\maketitle
\begin{abstract}
We consider the possibility of singularity formation 
in vortex sheets acted on by surface tension.
In the evolution of a nearly flat vortex sheet with viscosity and surface tension
neglected, Kelvin-Helmholtz instability has been shown to lead to the development
of  curvature singularities on the sheet. This is not unexpected, as the linearized
equations of motion for this flow exhibit short wavelength instability: infinitesimal
perturbations grow at an infinite rate.  Since the inclusion of surface tension
stabilizes (on linear theory) the high wavenumber modes, it has
been presumed that it also suppresses the formation of singularities.
We examine the possibility of singularity formation
by generalizing
Moore's approximation for vortex sheets to include the presence of surface tension.
For the approximate theory, traveling wave solutions with
imaginary wave speed are found which develop curvature  singularities
on the interface  in finite time.
Computations of the traveling
wave solution are provided to reveal other features, and comparisons are
made between these solutions and previously computed solutions of the 
full equations.   
\end{abstract}

\section{Introduction}

Recent studies of the evolution of vortex sheets
in two-dimensional, incompressible flows
have shed much light on their behavior.  The attention
which vortex sheets have received is in part due to their importance as
models of  free
%high Reynolds number shear layers.
shear layers in high Reynolds number flows.  Examples of vortex sheet
flows  include
those in which the vorticity is produced by a velocity
discontinuity present initially in the flow, those in which
the vorticity is generated baroclinically along a curve of density
discontinuity, and those in which the vorticity is produced
by a jump in viscosity across the interface.  In each of these
examples,  an initial value problem can be prescribed in  which
the motion is subject to a well known instability:
Kelvin-Helmholtz instability in the first case,
Rayleigh-Taylor instability in the second, and Hele-Shaw instability
in the the third.

When physical regularizations such as surface tension are absent,
the form of these instabilities is such that the shortest wavelength
components of a small perturbation grow the fastest \cite{jo:sho}.
Unfortunately,
such behavior is disastrous from the standpoint of numerical computation,
as mesh refinement leads to fast growing short wavelength modes
which quickly swamp a calculation.  If this pattern of growth
persists into the non-linear regime, then the motion is ill-posed;
small deviations in the initial condition will produce
significant changes in the ensuing motion.
In addition, singularities can form on the interface in finite
time.

Singularity formation is one of the most interesting phenomenon for such 
free surface flows, and has received a lot of attention in the 
literature over the past 15 years (Caflisch \cite{ca:mat} includes a brief review).
These singularities, if they occur, are important to mathematical theory
and numerical computation, while their physical significance depends on
the particular problem.  
For example, in  flows subject to the Kelvin-Helmholtz
and Rayleigh-Taylor instabilities, the type of singularity may
be related to the type of the subsequent roll-up of the
 sheet.  This, in turn, determines the rate of
mixing in the flow.

Recently, a new picture of singularity formation  
in vortex sheets has emerged by analytically
extending the equations for the location of the interface into the
complex domain.  Specifically, let $(x(\alpha,t),y(\alpha,t))$
be a parametric representation of the interface at time
$t$, and consider $\alpha$ to be a complex valued variable.
This continuation into the complex plane is a mathematical 
procedure with no physical meaning; the physical interface 
corresponds to the restriction of $x$ and $y$ to real $\alpha$.
Nevertheless, singularities in $x$ and $y$ which may be
present in the complex plane have a long range influence on
the form of the physical interface.  In addition, these singularities 
can move towards the real $\alpha$ axis.  When one reaches
this axis, it actually appears on the interface.

For periodic disturbances to a flat vortex sheet
with non-zero mean vorticity, 
asymptotic
studies by Moore \cite{mo:num}, \cite{mo:spo} and by Caflisch, Orellana, and Siegel \cite{cos}
show how singularities in the complex $\alpha$ plane
move towards and reach the real $\alpha$ axis in finite time,
upon which a curvature singularity forms on the interface.
These studies are based on a simple approximate
theory for the evolution of vortex
sheets using using a system of differential equations in the
parameter $\alpha$ (here $\alpha$ has a physical interpretation
as the circulation along the sheet).  A major advantage of the
approximate description is that it is local in $\alpha$,
as opposed to the exact description, which is given by a non-local
integro-differential equation.  
Numerical calculations for  the full equations  by
Meiron, Baker and Orszag \cite{me:ana}, Krasny \cite{kr:sin}, and Shelley \cite{sh:stu}
verify most of the predictions of the asymptotic studies.

A similar set of approximate equations were
derived for the Rayleigh-Taylor instability in \cite{cos}.
Baker, Caflisch and Siegel \cite{ba:sin} obtained a family of exact
solutions to the approximate equations for which singularities develop on
the fluid interface.  
The resulting predictions for the time, type
and location of singularity formation were directly verified
by numerical computations of the full equations, for
a wide range of Atwood numbers (a measure of the density difference 
 of the fluids) \cite{ba:sin}.   Numerical
 evidence for singularities in the Rayleigh-Taylor instability
had been previously supplied by Pugh \cite{Pu:the} for the Boussinesq limit.
Tanveer \cite{ta:sin} provided evidence that the single fluid case does
not exhibit singularity formation.
A related investigation of singularities for Boussinesq bubbles
was performed by Cowley and Pugh \cite{cp:bous}.  


The large body of numerical and analytical work cited above
provides strong evidence that, for  vortex sheet flows in
the absence of regularization, singularities can develop in finite time.
In order to prevent singularities from reaching the real
axis, it is necessary to introduce some form of regularization. 
 Numerical computations
have been used to examine various regularizations (see, for example, 
\cite{kr:des}, \cite{tr:num}, and \cite{ba:vor}).  These computations have been able to proceed beyond 
the singularity time and show the formation of rolled-up vortex sheets.
However, many of these regularizations do not correspond to a real physical
effect (an exception is given in \cite{ba:vor}).

One of the most prominent forces present in real interfacial flows 
is that of surface tension.
This force  stabilizes the high wavenumber modes in each of
the Kelvin-Helmholtz, Rayleigh-Taylor, and Hele-Shaw instabilities.
Since surface tension removes the ill-posedness in the linear
theory, it has been suggested that it does also in the
non-linear case \cite{bi:hel}.  However, we know of no published analytical study
which has  demonstrated this in  the instabilities mentioned above.
A recent numerical treatment by Hou, Lowengrub and Shelley \cite{ho:rem}
suggests that curvature singularities form precisely at the instant
that the interface self-intersects.  Other numerical studies are \cite{pu:nst},
\cite{ra:non},  and \cite{ba:sta}.

In this paper, we address the issue of singularity
formation in the Kelvin-Helmholtz instability with
surface tension.  This is done by generalizing Moore's approximation
to include the effects of surface tension.
As in \cite{ba:sin}, exact traveling wave solutions to the approximate equations
are found which develop singularities in finite time.   From these
solutions, predictions are made concerning the type, time, and location 
of singularity formation.
The predictions of the type of singularity are
consistent with the results of Hou, Lowengrub, and Shelley (M. Shelley, private communication).

Our approximation is formally valid for small deformations of the interface
and for singularities that are far from the real axis in the
complex plane.  In the present study, the theory is applied for 
large deformations and for singularities that are near the real axis, i.e., 
beyond the range of its formal validity.  Nevertheless, the success of this
approach in the previous study \cite{ba:sin}, in which
the theory is also applied in a regime beyond its formal validity, motivates
the present study.

The exact equations of motion for a vortex sheet acted on by
gravity and surface tension are presented in \S2.
In \S3, the system of localized equations is derived,
and in \S4 traveling wave solutions to this system are obtained.
Computations of these solutions are performed in \S5.  In \S6,
we discuss several discrepancies between the traveling wave solutions
and previously computed solutions to the full 
equations.  Concluding remarks are presented in \S8.


 


\section{Governing Equations}

Consider the two-dimensional motion of two immiscible, incompressible, 
inviscid fluids separated by a sharp interface.  The subscripts $1$ and $2$ are 
used to refer to the lower and upper fluids, respectively.  We assume that the 
two fluids have constant densities $\rho_1$ and $\rho_2$.  The interfacial 
shape is assumed to be periodic in the horizontal direction with period $2\pi$.  
It is further assumed that there is no circulation in the interior of each fluid, 
so that any vorticity as confined to be on the interface.
Thus the interface is a vortex sheet.

It is convenient to use complex coordinates to specify the position of the 
interface.  We let $z(\alpha ,t)=x(\alpha ,t)+iy(\alpha ,t)$ denote the position 
of the interface in the complex plane in terms of the time $t$ and a parameter 
$\alpha$ that labels a particular marker particle on the sheet.

Define the velocities $\vec u_1$ and $\vec u_2$ as the limiting values of the 
fluid velocities as the interface is approached from below and above, 
respectively.  Conservation of mass requires that they have equal normal 
components, i.e.,
\[ \vec u_1\cdot\vec n = \vec u_2\cdot \vec n\,.\]
In the absence of viscosity, however, the tangential components can be 
discontinuous.  
The vorticity at a point $z(p)$ on the interface is due to a jump in the tangential 
velocity of the fluids near $z(p)$.  We define the vortex sheet strength 
$\widetilde{\gamma}$ by
\[ \widetilde{\gamma}=(\vec u_1-\vec u_2)\cdot\widehat{s}\]
where $\widehat s$ is a unit vector tangent to the sheet at $z(p)$, pointing 
to the right relative to the direction from lower to upper fluid.  It is also 
convenient to define the unnormalized vortex sheet strength $\gamma$ by
\[\gamma =\widetilde{\gamma}\vert z_\alpha\vert\,.\]

Although the normal component of the interfacial velocity $\vec U$ is required to 
satisfy $\vec U\cdot\vec n = \vec u_1\cdot \vec n=\vec u_2\cdot\vec n$, there is 
freedom in defining the tangential component.  Following Baker et.al. \cite{ba:gen}, we 
choose $\vec U$ to be a weighted average of $\vec u_1$ and $\vec u_2$, i.e.,
\[\vec U={1 +\beta \over 2}\vec u_1+{1 -\beta \over 2}\vec u_2\]
where $\beta\in [-1,1]$.  Note that by choosing $\beta =-1$ the Lagrangian motion 
follows the upper fluid, whereas for $\beta =0$, this motion satisfies $\vec U =
{1 \over 2}(\vec u_1+\vec u_2)$, the average velocity.  In numerical simulations, 
this definition of the sheet velocity allows one to obtain some control over the 
positioning of the Lagrangian markers.  We also find that a particular choice of 
$\beta$ simppifies the analysis of \S 4.

Given the position and strength of the vortex sheet, its complex velocity 
$\partial\overline{z}/\partial t$ is prescribed by the following
formula  derived by Baker et. al.  \cite{ba:gen},
\begin{eqnarray}
\overline{z}_t(\alpha ,t) &= &{1 \over 2\pi i}PV\int^\infty_{-\infty}{\gamma (
\alpha ',t) \over z(\alpha ,t)-z(\alpha ',t)}d\alpha ' +{\beta\gamma (\alpha ,t) 
\over2z_\alpha (\alpha ,t)} \label{eq:1}
\end{eqnarray}
where the integral is in the Cauchy principal value sense.  Subscripts $\alpha$ 
and $t$ refer to differentiation, and an overbar denotes complex conjugate.  The 
integral term in (\ref{eq:1}) was first
derived by Birkhoff and is known as the Birkhoff-Rott integral; it
gives the mean of 
the two fluid velocities as a point $z(\alpha ,t)$ on the interface is approached 
from either side.  The second term on the right hand side of (\ref{eq:1}) results from 
our definition of $\vec U$ as a weighted average of $\vec u_1$ and $\vec u_2$.

In order to complete the description of the motion, an evolution equation for 
$\gamma$ is required.  For a vortex sheet which is acted upon by gravitational 
and surface tension forces, this equation is given by
\begin{eqnarray}
\gamma_t(\alpha ,t) &= &2A\left[Re\left\{{\partial\overline{q} \over 
\partial t}z_\alpha\right\}-{1 \over 2}\beta\gamma Re\left\{{q_\alpha 
\over z_\alpha}\right\} +{1 \over 8}{\partial \over \partial\alpha}~{\gamma^2 \over 
z_\alpha\overline{z}_\alpha}+g Im \{z_\alpha\}\right]\nonumber\\
&+ &{1 \over 2}\beta {\partial \over \partial \alpha}~{\gamma^2 \over z_\alpha 
\overline{z}_\alpha}-{2\tau \over (\rho_1+\rho_1)}~{\partial \over \partial\alpha}
\left({1 \over R}\right) \label{eq:2}
\end{eqnarray}
in which $Re$ and $Im$ denote the real and imaginary parts.  
Here $A$ is the Atwood ratio, defined as
\[A = {\rho_2-\rho_1 \over \rho_1+\rho_2}\,,\]
$g$ is the acceleration due to gravity, $R$ is the radius of curvature at the 
point $\alpha$ on the sheet (defined as positive when fluid $1$ protrudes into 
fluid $2$), $\tau$ is the surface tension coefficient and 
$\overline{q}$ is the Birkhoff-Rott integral
\begin{eqnarray}
\overline{q} &= &{1 \over 2\pi i}PV\int^\infty_{-\infty}{\gamma (\alpha ',t) \over 
z(\alpha ,t)-z(\alpha ',t)} d\alpha '\,. \label{eq:3}
\end{eqnarray}
Equation (\ref{eq:2}) is obtained by subtracting the tangential components of 
Euler's 
momentum equation across the interface.  The difference between pressures on 
either side of the interface is expressed in terms of interface quantities using 
Laplaces equation $p_1-p_2=\tau /R$.  A detailed derivation of equation 
(\ref{eq:2}) 
without the surface tension term is given by Baker, Meiron and Orszag \cite{ba:gen}.  
Derivations wich include the surface tension term are provided by Tryggvason \cite{tr:num} 
and Pullin \cite{pu:nst}.

We comment briefly on the physical meaning of each of the terms in the vorticity 
evolution equation (\ref{eq:2}).  The term
\begin{eqnarray}
\Bigg[ Re\left\{{\partial \overline{q} \over \partial t} z_\alpha
\right\} &+ &{1 \over 8}~{\partial \over \partial\alpha}~{\gamma^2 \over z_\alpha 
\overline{z}_{\alpha}}\Bigg] \label{eq:4}
\end{eqnarray}
describes the rate of change of the (unnormalized) vortex sheet strength following 
marker particles due to the differential tangential accelerations of the fluids 
above and below the interface \cite{tr:num}.
If $\beta
\not= 0$, the terms $\displaystyle {1 \over 2}\beta {\partial \over \partial
\alpha}~{\gamma^2 \over z_\alpha\overline{z}_\alpha}$ and $\displaystyle A\beta
\gamma Re\left\{{q_\alpha \over z_\alpha}\right\}$ are corrections to 
this rate due to the use of a weighted average for the interface velocity.  
The terms $2AgIm\{z_\alpha\}$ and $\displaystyle {2\tau \over (\rho_1+\rho_2)}~
{\partial \over \partial\alpha}\left({1 \over R}\right)$ describe the rate of 
change of vorticity due to the action of gravity 
and surface tension, respectively.  We note that the surface tension term $
\displaystyle{{2
\tau \over (\rho_1+\rho_2)}~{\partial \over \partial\alpha}\left({1 \over R}
\right)}$ acts so as to oppose any interfacial growth or contortion \cite{ra:non}.  
This contrasts the situation in the break up of a 3-D cylindrical jet, in which 
the effect of surface tension is to induce the break up if the disturbance 
wavelength is larger than the jet circumference.

A simple solution to equations (\ref{eq:1},\ref{eq:2}) is
\begin{eqnarray}
z(\alpha ,t) &= &\alpha +{\beta\gamma_0 \over 2}t\,,~~\qquad~~\gamma (\alpha ,t)
=\gamma_0\,. \label{eq:5}
\end{eqnarray}
This solution corresponds to a flat interface with uniform vortex sheet\break 
strength $\gamma_0$, with the marker particles translating to the left with 
speed $\displaystyle {\beta \gamma_0 \over 2}$.  This would be an equilibrium 
solution in an Eulerian frame of reference.

The linearized version of equations (\ref{eq:1},\ref{eq:2}) about the flat 
interface (\ref{eq:5}) permits solutions of the form
\begin{eqnarray*}
z(\alpha ,t) &= &\alpha +{\beta\gamma_0 \over 2} t+Ae^{ik\alpha\pm
\sigma (k )t}\\
\gamma (\alpha ,t) &= &\gamma_0+Be^{ik\alpha\pm\sigma (k )t}
\end{eqnarray*}
where
\begin{eqnarray}
\sigma (k) &= &\frac{ik(A+B)\gamma_0\pm [(1-A^2)\gamma^2_0k^2-4\tau \vert k
 \vert^3+4\Ag \vert k \vert
]^{1/2}}{2} \label{eq:6}
\end{eqnarray}
and, since we have assumed $2\pi$ periodic disturbances, $k$ is an  
integer.  There is an additional branch to the dispersion relation, $\sigma =0$, 
which is usually termed a labeling mode and is an artifact of the Lagrangian 
approach (see \cite{ba:sin}).

The dispersion relation has several features which are relevant to this study:

\vspace{.2in}
\pist{1.}{The choice $\beta =-A$ in (\ref{eq:6}) removes the imaginary part of 
$\sigma$; i.e., it removes the translational component of $\sigma$.}

\vspace{.12in}
\pist{2.}{For $\gamma_0\not= 0$ and $\tau =0$, the real part of $\sigma$, which 
is the growth (or decay) rate of the linear modes, is proportional to $\vert 
k\vert$ for large $k$.  This growth rate is an indication of Kelvin-Helmholtz 
instability. For $ g \neq 0$ and $
\gamma_0=\tau=0$, the real part of $\sigma$ is proportional to 
$\displaystyle \vert k \vert^{\onehalf}$ for large $k$. This growth rate is an indication
of Rayleigh-Taylor instability.} 

\vspace{.12in}
\pist{3.}{For $\displaystyle \gamma_0\not= 0$ and $ \displaystyle 0<\tau <{(1-A^2)\gamma^2_0 
\over 4}+Ag$, there are a finite but non-zero number of unstable modes.  
Above  a cut-off wavenumber $k_c$, given by
\[k_c={(1-A^2)\gamma^2_0+[(1-A^2)^2\gamma^4_0+32\tau Ag ]^{1/2} \over 
8\tau},\]
each $\sigma (k)$  is pure imaginary.  Thus 
interfacial surface tension has a stabilizing influence.}

\vspace{.12in}
\pist{4.}{For $\gamma_0\not= 0$ and $\displaystyle \tau\ge {(1-A^2)\gamma^2_0 
\over 4} +Ag$ every mode $\sigma (k)$ is pure imaginary, indicating complete 
suppression of the linearized instability for large surface tension.}

\vspace{.2in}
Birkhoff \cite{bi:hel} was the first to argue that the unbounded growth rate in the 
linearized dispersion relation for vortex sheets with $\gamma_0\not= 0$, $\tau 
=0$ allows the possibility of singularity development.  Since then, analytic 
studies by Moore \cite{mo:num}, \cite{mo:spo},
Caflish and Orellana \cite{ca:ill} and Duchon and Robert \cite{du:glo}, 
and numerical studies by Meiron, Baker and Orszag \cite{me:ana}, Krasny \cite{kr:sin}, and 
Shelley \cite{sh:stu} have all supported the conclusion that singularities form in 
finite time.  If $\gamma_0=0$, Birkhoff's argument for singularity development based 
on linear analysis fails.  Nevertheless, studies by Pugh \cite{Pu:the} and Baker, 
Caflisch and Siegel \cite{ba:sin} in the classical Rayleigh Taylor instability (i.e., 
for $\gamma_0=0$, $\tau =0$ and $g\not= 0$) conclude that singularities do 
develop 
in finite time.  With singularity development in the unregularized problem 
generally well established, we turn our attention to the question of singularity 
formation in the problem regularized (on linear theory) by surface tension.  
Under the action of 
an imposed shear and surface tension, the vortex sheet is subject to two 
competing tendencies:  the tendency of the imposed shear to advect vorticity 
along the sheet toward an accumulation point leading to the build up of curvature, 
and the tendency of surface tension forces to avoid the build up of curvature.  
We develop a simple approximate theory which suggests that,
despite the presence of surface tension, 
curvature singularities form.

\section{The Localized Approximation}

We shall derive a simple approximate theory for the evolution of vortex sheets 
acted on by surface tension and gravity.  The approximate theory was first 
developed by Moore \cite{mo:spo} in order to analyze singularity formation in the 
Kelvin-Helmholtz problem.  It was later extended by Caflisch, Orellana and Siegel 
\cite{cos} to other interfacial flows, including the Rayleigh-Taylor instability,
and by Caflisch \cite{ca:sfo} to smooth flows.  For 
our purposes, we need only generalize the derivation given in \cite{cos} 
to include the surface tension term 
and the existence of a non-zero mean sheet vorticity.  However, for completeness, 
we shall repeat the main points of the theory.

The approximation method is most easily described by extending the parameter 
$\alpha$ to be complex.  Define an operation $\ast$ by
\[ z^\ast (\alpha )=\overline{z}(\overline{\alpha})\,.\]
Note that $z^\ast$ is the unique analytic function which equals $\overline{z} 
(\alpha )$ for $\alpha$ real.  With this definition, the evolution 
equations (\ref{eq:1},\ref{eq:2}) can be analytically extended to complex $\alpha$ as
\[\partial_t z^\ast = q^\ast +{\beta\gamma \over 2z_\alpha}\]
\begin{eqnarray}
\partial_t\gamma &= &A\Bigg[ q^\ast_tz_\alpha +q_tz^\ast_\alpha -{1 \over 2}
\beta\gamma\left({q_\alpha \over z_\alpha}+{q^\ast_\alpha \over z^\ast_\alpha}
\right)\nonumber\\
&~ &-ig(z_\alpha -z^\ast_\alpha )\Bigg] +{A+2\beta \over 4}~{\partial \over 
\partial\alpha}\left({\gamma^2 \over z_\alpha z^\ast_\alpha}\right)\label{eq:7}\\
&~ &-{2\tau \over (\rho_1+\rho_2)}~{\partial \over \partial\alpha}\Im
{z^\ast_{\alpha\alpha}z_\alpha \over (z_\alpha z^\ast_\alpha )^{3/2}}\nonumber
\end{eqnarray}
where the complex velocity $q^\ast$ is given by
\begin{eqnarray}
q^\ast (\alpha ,t) &= &{1 \over 2\pi i} PV\int^\infty_{-\infty}{\gamma (\alpha +
\xi ,t) \over z(\alpha ,t)-z(\alpha +\xi ,t)}\,d\xi\,.\label{eq:8}
\end{eqnarray}

Knowledge of the behavior of the kernel in the above principal value 
integral (i.e., the positions and types of singularities, and behavior at infinity) 
can enable the integral to be evaluated by applying the residue theorem over a 
suitably deformed contour.  This motivates a decomposition of the kernel into a 
sum of terms with known analytic behavior in the complex plane, plus a remainder term.  
To implement this kernel decomposition, we first decompose $\gamma$ and $z$ as
\begin{eqnarray}
z(\alpha ,t) &= &\alpha +s_+(\alpha ,t)+s_-(\alpha ,t)+s_0(t)\nonumber\\
\gamma (\alpha ,t) &= &\gamma_+(\alpha ,t)+\gamma_-(\alpha ,t)+\gamma_0
(t)\label{eq:9}
\end{eqnarray}
in which $s_+$, $\gamma_+$ contain the positive wavenumber components of $s$, 
$\gamma$ respectively and $s_-$, $\gamma_-$ contain the negative wavenumbers, 
i.e.,
\begin{eqnarray*}
s_+(\alpha ,t) &= &\sum_{k>0}\widehat{s}_k(t)e^{ik\alpha}\\
s_-(\alpha ,t) &= &\sum_{k<0}\widehat{s}_k(t)e^{ik\alpha}\\
\gamma_+(\alpha ,t) &= &\sum_{k>0}\widehat{\gamma}_k(t)e^{ik\alpha}\\
\gamma_-(\alpha ,t) &= &\sum_{k<0}\widehat{\gamma}_k(t)e^{ik\alpha}\,.
\end{eqnarray*}
If $s$ and $\gamma$ are analytic in the strip $\vert\Im~\alpha\vert <\rho$, then 
$s_+$ and $\gamma_+$ are analytic in the upper half plane $\Im~\alpha >-\rho$ and 
are called {\it {upper analytic}}.  Similarly, $s_-$ and $\gamma_-$ are analytic 
in the lower half plane $\Im~\alpha <\rho$ and are called {\it {lower analytic}}.  
The constant term in $z$, denoted by $s_0(t)$, plays no role in the resulting 
equations, since it does not effect the values of the Fourier coefficients 
$\widehat{s}_k(t)$, $\widehat{\gamma}_k(t)$ for $k\not= 0$.  
The constant term $\gamma_0$ in the representation of $\gamma$ is the mean vortex 
sheet strength.  Denote
\begin{eqnarray}
z_{\alpha+}=1+s_{\alpha +}~~&\qquad &~~ z_{\alpha -}=1+s_{\alpha -}\,.\label{eq:11}
\end{eqnarray}
Note that

\vspace{.2in}
\pist{1.}{The $\ast$-operator switches + components to minus components, i.e.}
\[(s^\ast )_+=(s_-)^\ast\,,~~~~(s^\ast )_-=(s_+)^\ast\,.\]

\pist{2.}{Since $\gamma$ is real when $\alpha$ is real, $\widehat{\gamma}_{-k}
=\overline{\widehat{\gamma}}_k$ so that}
\[(\gamma_+)^\ast =\gamma_-\,.\]

The analytically extended equations (~~~) can be written as
\begin{eqnarray*}
z^\ast_t &= &B_1[z_\alpha ,z^\ast_\alpha ,\gamma,\gamma^\ast ]\\
\gamma_t &= &B_2[z_\alpha ,z^\ast_\alpha ,\gamma ,\gamma^\ast ]\,.
\end{eqnarray*}
For any function as operation $A[s]$,
such as the kernel in the
principal value integral (\ref{eq:8}) or any other 
terms in (\ref{eq:7}), we can write
\[A[s] = A[s_+]+A[s_-]+E[s_+,s_-]\,.\]
This is just a definition of the term $E$.  The approximation takes the form
\[ A[s]\sim A[s_+]+A[s_-]\,.\]
If $s$ has a constant term, $s_0$, the corresponding approximation is
\[A[s]\sim A[s_++s_0]+A[s_-+s_0]-A[s_0]\,.\]
The neglected term $E$ contains only cross products of upper and lower analytic 
functions such as $s_+s_-$.  Pure products such as $s^2_+$ or $s^2_-$ are 
contained in $A[s_+]$ and $A[s_-]$.  Thus the approximation amounts to neglecting 
all cross products of upper and lower analytic functions.  It is shown in \cite{cos} 
that the size of a typical neglected term of the form $s_+s_-$ satisfies
\[\Vert s_+s_-\Vert\le e^{-2\rho}\Vert s_+\Vert~\Vert s_-\Vert\]
where $\Vert s\Vert =\sup\limits_{\vert Im~\alpha\vert\le\rho}\vert s(\alpha )
\vert$.  For $\rho >>1$, this is negligible.

The approximation is based on the following physical argument:  Non-linear 
interactions cause wavenumbers to add so that energy flows back and forth 
among wavenumbers.  However, the process of singularity development primarily 
involves flow of energy from low wavenumbers to high wavenumbers.  Our 
approximation keeps this flow of energy to higher wavenumbers, but ignores 
all flow of energy from high to low wavenumbers, i.e., all interactions between 
positive wavenumbers and negative wavenumbers are neglected.  Mathematically, 
this is the same as keeping interactions of $+$ terms with $+$ terms and $-$ 
terms with $-$ terms, but neglecting interactions of $+$ with $-$.

In the present situation, we apply the approximation for both $s$ and $\gamma$.  
Therefore, (\ref{eq:7}) becomes
\begin{eqnarray}
z^\ast_t &= &B_1 (z_{\alpha +,}(z^\ast )_{\alpha +,}\gamma_{+,}(\gamma^\ast )_+) \nonumber\\
&~ &+B_1(z_{\alpha -,}(z^\ast )_{\alpha -,}\gamma_{-,}(\gamma^\ast )_-) \nonumber\\ 
\gamma_t &= &B_2(z_{\alpha+,}(z^\ast)_{\alpha +,}\gamma_{+,}(\gamma^\ast )_+) \nonumber\\
&~ &+B_2(z_{\alpha -,}(z^\ast )_{\alpha -,}\gamma_{-,}(\gamma^\ast )_-)\,.
\label{eq:11.5}
\end{eqnarray}
For example, in the second equation of (\ref{eq:7}), the surface tension term is 
approximated as
\begin{eqnarray*}
{2\tau \over (\rho_1+\rho_2)} {\partial \over \partial\alpha}\Im{z^\ast_{\alpha
\alpha}z_\alpha \over (z_\alpha z^\ast_\alpha )^{3/2}}
&\sim &{\tau \over i(
\rho_1+\rho_2)}~{\partial \over \partial\alpha}\Bigg[ {s^\ast_{\alpha\alpha -}
s_{\alpha +}+s^\ast_{\alpha\alpha -}-s_{\alpha\alpha +}s^\ast_{\alpha -}-
s_{\alpha\alpha +} \over (1+s_{\alpha +}+s^\ast_{\alpha -}+s_{\alpha +}
s^\ast_{\alpha -})^{3/2}}\\
&~ &+{s^\ast_{\alpha\alpha +}s_{\alpha -}+s^\ast_{\alpha\alpha +}-s_{\alpha
\alpha -}s^\ast_{\alpha +}-s_{\alpha\alpha -} \over (1+s_{\alpha -}+s^\ast_{\alpha 
+}+s_{\alpha -}s^\ast_{\alpha +})^{3/2}}\Bigg]
\end{eqnarray*}
where we have used (\ref{eq:11}).  Each of the other terms 
in (\ref{eq:7}) is similarly approximated. 

An advantage of this approximation is that now the Birkhoff-Rott integral $q^\ast
[s,\gamma ]$ can be evaluated by integration along a suitably deformed contour, 
resulting in a purely local term coming from a residue.  By using this 
evaluation, then projecting equations (\ref{eq:11.5}) onto upper and lower 
analytic function space, we obtain the following 
system:
\begin{eqnarray}
{\partial s^\ast_- \over \partial t} &= &{\beta -1 \over 2}\left({\gamma_0+
\gamma_+ \over 1+s_{\alpha +}}\right)\\
{\partial s_+ \over \partial t} &= &{\beta +1 \over 2}\left({\gamma_0+\gamma_+ 
\over 1+s^\ast_{\alpha -}}\right)\\
{\partial\gamma_+ \over \partial t} &= &{A+\beta \over 2}~{\partial \over 
\partial\alpha}\left({(\gamma_0+\gamma_+)^2 \over (1+s_{\alpha +})(1+s^\ast_{\alpha 
-})}\right) -i\Ag (s_{\alpha +}-s^\ast_{\alpha -})\\
&~ &+{i\tau \over (\rho_1+\rho_2)}~{\partial \over \partial\alpha}\left[
{s^\ast_{\alpha\alpha -}s_{\alpha +}+s^\ast_{\alpha\alpha -}-s_{\alpha\alpha +}
s^\ast_{\alpha -}s_{\alpha\alpha +} \over (1+s_{\alpha +}+s^\ast_{\alpha -}+
s_{\alpha +}s^\ast_{\alpha -})^{3/2}}\right]\,.  \label{eq:12}
\end{eqnarray}
In obtaining the above, we have used the relation $(\gamma^\ast )_-=\gamma_+$, 
and have applied the $\ast$ operator to go from an equation for $\displaystyle 
{\partial s^\ast_+ \over \partial t}$ to one for $\displaystyle {\partial s_+ 
\over \partial t}$.  For details of the localized approximation procedure and
its application in similar problems, see \cite{cos}.

It is convenient to write (\ref{eq:12}) as a quasi-linear system.  Denote
\begin{eqnarray}
\phi &= &1+s_{\alpha +}~\quad~\psi =1+s^\ast_{\alpha -}~\quad~\omega =
\gamma_+ \label{eq:12.1}
\end{eqnarray}
and differentiate to first two equations of (\ref{eq:12}) to obtain
\begin{eqnarray*}
\phi_t &= &{\beta +1 \over 2}~{\partial \over \partial\alpha}\left({\omega +
\gamma_0 \over \psi}\right)\\
\psi_t &= &{\beta -1 \over 2}~{\partial \over \partial\alpha}\left({\omega +
\gamma_0 \over \phi}\right)
\end{eqnarray*}
\begin{eqnarray}
\omega_t &= &{A+\beta \over 2}~{\partial \over \partial\alpha}\left({(\omega +
\gamma_0)^2 \over \phi\psi}\right) -i\Ag (\phi -\psi )\nonumber\\
&~ &+{i\tau \over (\rho_1+\rho_2)}~{\partial \over \partial\alpha}\left[
{\psi_\alpha\phi -\phi_\alpha\psi \over (\phi\psi )^{3/2}}\right]\,. \label{eq:13}
\end{eqnarray}

The localized equations are much simpler than the full 
equations (\ref{eq:1},\ref{eq:2}), mainly in that they are purely local 
differential equations rather than non-local integro-differential 
equations.  However, as described in \cite{ba:sin}, there is a defect in the 
approximate system concerning the role of $\beta$.  For the full equations, 
solutions obtained for different values of $\beta$ are equivalent.  The 
parameterization of the solution in $\alpha$ depends on the particular value of 
$\beta$, but the resulting 
curve does not.  This is not the case for the approximate system,
since different values of $\beta$ lead to different approximations.   
The reason for this is simple.
Different choices of $\beta$ correspond to different 
parameterizations by $\alpha$ for the motion.  The approximation involves 
Fourier transformation to make the decompositions into upper and lower analytic 
functions.  Since changes of variables are not easily expressed through the 
Fourier transform, there is no simple relation between solutions of the 
approximate system for different choices of $\beta$.  Nevertheless, using an 
argument analogous to that given in Caflisch and Orellana \cite{ca:lon} or Baker, Cowley, 
Tanveer and Page \cite{co:asy}, one can establish that solutions to the approximate system for 
various $\beta$ are asymptotically close to the exact solution when the 
deformations of the interface are small and the singularities are far from the 
real axis.  Consequently, the nature of the singularities are predicted reliably 
by the approximate system independently of $\beta$.

We next consider special solutions of the localized equations (\ref{eq:13}) which take 
the form of traveling wave solutions in the complex plane.  These solutions will 
be shown to have singularities which move with constant speed toward the  real 
$\alpha$ axis.  The singularities reach this axis at a time $t=t_c$, upon which 
a curvature singularity appears on the fluid interface.

\section{Traveling Wave Solutions}

We seek solutions to the localized equations (\ref{eq:13}) that are $2\pi$ periodic in 
space.  Note that these equations contain only upper analytic variables, namely 
$\phi$, $\psi$ and $\omega$.  The only other coefficients in the system are 
constants, i.e., there are no coefficients or terms which are known functions of 
$\alpha$ and $t$.  An unusual property of such systems (or analogous lower 
analytic systems) is that they often allow exact solutions which take the form of 
periodic traveling waves in the complex plane.  Furthermore, the speed of a 
given traveling wave solution is independent of the wave amplitude and identical 
to the speed given by a linear analysis.  This contrasts the behavior of 
traveling wave solutions to general non-linear PDE's, which have speeds that 
typically depend on the amplitude of the wave.   

The existence of such traveling wave solutions in upper analytic PDE's can be 
illustrated by considering a single quasi-linear first order PDE in $\alpha$ 
and $t$ of the form
\begin{eqnarray}
f_{t+}+L[f_+]f_{\alpha +} &= &R[f_+]\label{eq:14}
\end{eqnarray}
where $f_+$ is an upper analytic variable.  It is assumed that $L[f_+]$ and $R[
f_+]$ admit a Taylor's series expansion in $f_+$ about some neighborhood of 
$f_+=0$, and that the coefficients of the Taylor's expansion are constants 
independent of $\alpha$ and $t$.  Furthermore, it is assumed that $R[0]=0$.  
We then look for solutions of the form
\begin{eqnarray}
f_+ &= &a_1e^\xi +a_2e^{2\xi}+\cdots~~\mbox{\rm with}~~\xi =i(\alpha +\sigma t)
\,,\label{eq:15}
\end{eqnarray}
where the $a_i$'s are constants.  
This is just a Fourier expansion of a periodic traveling wave with speed 
$\sigma$.  The absence of a term $a_0$ is justified by the assumption $R[0]=0$.  
Substituting this form into (\ref{eq:14}) and equating coefficients to like powers of 
$e^\xi$ leads to the following hierarchy of equations
\begin{eqnarray}
a_1F_1(\sigma ) &= &0\nonumber\\
F_2(\sigma ;a_1,a_2) &= &0\nonumber\\
~\quad\vdots\quad~ &~ &~\nonumber\\
F_j(\sigma ;a_1,a_2,\dots ,a_j) &= &0\nonumber\\
~\quad\vdots\quad~ &~ &~ \label{eq:16}
\end{eqnarray}

The $F_i$'s are functions of the indicated variables.  For the first equation to 
have non-trivial solutions with $a_1\not= 0$, $\sigma$ must satisfy $F_1 (\sigma 
)=0$.  Thus $\sigma$ is determined in the first (linear) order of the hierarchy.  
The remaining equations successively determine $a_j$ as a function of $a_{j-1},
\dots ,a_1$; $\sigma$.  The fact that (\ref{eq:14}) is a purely upper analytic equation 
and allows an outflow of energy to high wavenumbers while avoiding inward 
flow is crucial in producing the structure of this hierarchy, and results in the 
solution (\ref{eq:15},\ref{eq:16}).  Note that the traveling
 wave solution to (\ref{eq:14}) will be 
uniquely determined if we require that $f_+\sim a_1e^\xi$ as $\xi\to -\infty$ 
for a given $a_1$.  We also remark that the existence of these special solutions 
makes the localized approximation a useful procedure, even for operators which 
are already local in form (see for example \cite{ca:sfo}).

Traveling wave solutions with non-zero imaginary component of speed  are 
particularly significant with regard to singularity formation.  For these solutions, 
singularities travel in the complex plane toward the real $\alpha$ line with 
speed $\Im (\sigma )$.  A singularity $\alpha_s=x_s+iy_s$, with imaginary component 
$y_s$ at $t=0$, will hit the real $\alpha$ line at $\displaystyle t={y_s \over 
\Im (\sigma )}$.  Consequently a singularity in the traveling wave solution will 
occur on the interface in finite time, if the initial data has a singularity at 
any complex value of $\alpha$.

Motivated by the above remarks, we look for solutions to the localized equations 
of the form $\phi =\phi (\xi )$, $\psi =\psi (\xi )$ and $\omega =\omega (\xi )$ 
with $\xi =i(\alpha +\sigma t)$.  We require that $\phi\sim 1+\widehat{\phi}_1
e^\xi$, $\psi\sim 1+\widehat{\psi}_1e^\xi$ and $\omega\sim\widehat{\omega}_1 
e^\xi$ as $\xi\to -\infty$, and expect that $\sigma$ will be given by
\begin{equation}
\sigma = {(A+\beta )\gamma_0 \over 2}\pm i\left[{(1-A^2) \over 4}\gamma^2_0 
-\tau +Ag\right]^{1/2}  \label{eq:15.5}
\end{equation}
which is identical to the speed obtained from a linear analysis (this must be 
verified at the first order of the hierarchy described earlier).  Substitution 
of these forms into (\ref{eq:13}) gives the following system of ODE's:
\begin{eqnarray}
\sigma\phi_\xi &= &{\beta +1 \over 2}~{\partial \over \partial\xi}\left({\omega 
+\gamma_0 \over \psi}\right)\nonumber\\
\sigma\psi_\xi &= &{\beta -1 \over 2}~{\partial \over \partial\xi}\left({\omega +
\gamma_0 \over \phi}\right)\nonumber\\
\sigma\omega_\xi &= &{A+\beta \over 2}~{\partial \over \partial\xi}\left({(\omega 
+\gamma_0)^2 \over \phi\psi} \right) -\Ag (\phi -\psi )\nonumber\\
&~ &-{\tau \over (\rho_1+\rho_2)}~{\partial \over \partial\xi}\left[{\psi_\xi\phi 
-\phi_\xi \psi \over (\phi\psi )^{3/2}}\right]\,.\label{eq:17}
\end{eqnarray}
The first two equations can be integrated directly, yielding
\begin{eqnarray}
\phi &= &\phi (0) +{\beta + 1 \over 2\sigma}\left({\omega +\gamma_0 \over \psi}
-{\omega (0)+\gamma_0 \over \psi (0)}\right)\nonumber\\
\psi &= &\psi (0) +{\beta -1 \over 2\sigma}\left({\omega +\gamma_0 \over \phi} -
{\omega (0)+\gamma_0 \over \phi (0)}\right)\,.\label{eq:18}
\end{eqnarray}

By the definitions (\ref{eq:12.1}), the constant part in the Fourier expansions for both 
$\phi$ and $\psi$ is equal to one, whereas there is no constant part in the 
Fourier expansion of $\omega$.  Consequently, the right hand sides of (\ref{eq:18}) 
must have a mean value of $1$, and this implies that the following relationships 
must hold:
\begin{eqnarray}
\phi (0)+{\beta +1 \over 2\sigma}\left(\gamma_0-{\omega (0)+\gamma_0 \over 
\psi (0)}\right) &= &1\nonumber\\
\psi (0) +{\beta -1 \over 2\sigma}\left(\gamma_0 -{\omega (0)+\gamma_0 \over 
\phi (0)}\right) &= &1\,.\label{eq:19}
\end{eqnarray}
Define
\begin{eqnarray}
\Phi_0 &= &1 - {\beta +1 \over 2\sigma}\gamma_0\nonumber\\
\Psi_0 &= &1 - {\beta -1 \over 2\sigma}\gamma_0\,.\label{eq:20}
\end{eqnarray}
Substitution of the relations (\ref{eq:19}) and (\ref{eq:20}) into (\ref{eq:18}) then 
gives
\begin{eqnarray}
\phi &= &\Phi_0+{1 +\beta \over 2\sigma}~{\omega +\gamma_0 \over \psi}\nonumber\\
\psi &= &\Psi_0-{1 -\beta \over 2\sigma}~{\omega +\gamma_0 \over \phi}
\,.\label{eq:21}
\end{eqnarray}

It is possible to express the third equation in (\ref{eq:17}) as an ODE with a single 
dependent variable.  A convenient choice for the dependent variable is 
$\displaystyle \Omega=
{\omega +\gamma_0 \over \sigma\Phi_0\Psi_0}$.  We use equations (\ref{eq:21}) to 
express $\phi$ and $\psi$ in terms of $\Omega$:
\begin{eqnarray}
2\phi &= &\Phi_0(1+\Omega\pm\sqrt{1+2\beta\Omega +\Omega^2})\nonumber\\
2\psi &= &\Psi_0(1-\Omega\pm\sqrt{1+2\beta\Omega +\Omega^2})\,.\label{eq:22}
\end{eqnarray}
Here the square root is defined using the principal
branch cut, i.e. $(1+2\beta\Omega +\Omega^2)^{1/2}=\sqrt{M}e^{i\theta /2}$, where
$\theta\in [-\pi ,\pi ]$ is the argument of $1+2\beta\Omega +\Omega^2$, and
$M$ is its amplitude. 
% We remark that for $\beta =\pm 1$, this definition of
%branch leads to $(1+2\beta\Omega +\Omega^2)^{1/2}=-1-\beta\Omega$, so that the
%choice of the $-$ sign in (~~~), (~~~) is consistent with the requirement that
%$\phi$ and $\psi$ have mean $1$.
This choice of branch makes sense only if the
curve ${\cal C}$ traced out by $1+2\beta\Omega +\Omega^2$ as $\alpha$ goes from
$0$ to $2\pi$ does not intersect the negative real axis.  This is the case for
$\Omega '=\Omega -\Gamma_0$ small, where 
\be
\Gamma_0=\frac{\gamma_0}{ \sigma\Phi_0\Psi_0}
\label{eq:22.5}
\ee
is the mean of $\Omega$.
  The situation is more complicated when
${\cal C}$ intersects the negative real axis; this will be discussed momentarily.
The sign in front of the square root is selected by the requirement that the 
constant parts in the Fourier expansions of $\phi$ and $\psi$ equal one.  This 
requirement can lead to selection of either the $+$ or $-$ branch, depending on 
the values of the parameters $g$, $\gamma_0$ and $\tau$.  The branch of the 
square root which is to be selected in a specific case will be pointed out when 
that case is discussed in detail.  For now, we leave the $\pm$ in order to indicate
this indeterminacy.

Upon substituting (\ref{eq:22}) into the third equation of (\ref{eq:17}) and 
manipulating some, we obtain
\begin{eqnarray}
{\partial \over \partial\xi} &\Bigg\{& {(A+\beta )\sigma^2 \Omega^2 \over
v(\Omega)^2} -\sigma^2 \Omega
+\left({2 \over \Phi^3_0\Psi^3_0}\right)^{1 \over 2}
{\tau \over  (\rho_1+\rho_2)}
\frac{u(\Omega)}  {v(\Omega)}
\Omega_\xi\Bigg\}\nonumber\\
&= &{Ag \over 2}\left[\left({1 \over \Psi_0}-{1 \over \Phi_0}\right)(1+u(\Omega))
+\Omega\left({1 \over \Psi_0}+{1 \over \Phi_0}
\right)\right]\,.\label{eq:23}
\end{eqnarray}
where
\begin{equation}
u(\Omega)=\pm(1+2\beta\Omega+\Omega^2)^{\onehalf},~
v(\Omega)=(1+\beta\Omega \pm u(\Omega))^{\onehalf}.
\label{eq:23.5}
\end{equation}
Several situations lead to a considerable simplification of the above equation.  
As we have pointed out earlier, the choice $\beta =-A$ selects a convenient 
parameterization, and we see here that this choice eliminates the first term 
in (\ref{eq:23}).  The case $\tau =0$, $\gamma_0=0$, also leads to a marked 
simplification of this equation.  This case corresponds to an interface subjected 
to a pure Rayleigh-Taylor instability and is treated in detail in \cite{ba:sin}.  Since 
this is a particularly simple situation, yet previews some of the features of the 
more difficult case $\tau\not= 0$, we briefly review the results of \cite{ba:sin} here.  
For $\gamma_0=0$, $\Phi_0$ and $\Psi_0$ are both equal to $1$, and it is necessary 
to choose the $+$ branch in (\ref{eq:22},\ref{eq:23.5}) to satisfy the requirement 
that $\phi$ and $\psi$ have a constant Fourier component equal to $1$.  If $\beta$ 
is chosen to be $-A$, then equation (\ref{eq:23}) reduces to
\[\Omega_\xi = -\left({\Ag \over \sigma^2}\right)\Omega\,.\]
The requirement that $\Omega$ is upper analytic and $2\pi$ periodic in $\alpha$ 
forces $\displaystyle {\Ag \over \sigma^2}=-n$ for integer $n$.  A unique solution 
can be chosen by 
specifying the behavior of $\Omega$ as $\xi\to -\infty$.  If we select $\Omega
\sim\Omega (0)e^\xi$ as $\xi\to -\infty$, then $n=1$ and the resulting solution 
is
\begin{equation}
\Omega =\Omega (0)e^\xi~~\qquad~~\sigma =\pm i(\Ag )^{1/2}
\label{eq:24}
\end{equation}
where $\Omega (0)$ is a complex amplitude.  Note that the speed $\sigma$ is as 
expected from our linear analysis.

The traveling wave solutions correspond to a special choice of initial condition 
for the sheet position and sheet strength.  For these initial conditions, the 
nature of the singularity is easily deduced from the explicit solutions (\ref{eq:22}), 
(\ref{eq:24}).  In particular, for $\Omega (0)$ small, the functions $\phi$ and $\psi$ 
have square root singularities in the lower half complex $\alpha$ plane which 
move vertically upwards
% along the line
%        \be
%\alpha = i(Ag)^\onehalf ~ t-i \ln \left ( -A \pm i [1-A^2]^\onehalf
%\right ) +i \ln \|Omega(0)| \nonumber
%        \ee
 with constant speed.   These singularities reach the 
real $\alpha$ axis at positions
        \be
\alpha_c = \pi \pm \tan^{-1} \left(\frac{[1-A^2]^{\onehalf}}{A} \right),  \nonumber
        \ee
and at time
        \be
t_c= \frac{1}{(Ag)^{\onehalf}} \ln \frac{1}{|\Omega(0)|} \nonumber
        \ee
at which point they are physically realized.  The values for the time and 
position of singularity formation obtained from the traveling wave solution were 
directly verified by numerical computation of the full equations, with excellent 
agreement between the theoretical preductions and numerical results for values 
of Atwood number between $0$ and approximately $.9$.  Further details, including 
traveling wave solutions for other values of $\beta$, are given in \cite{ba:sin}.

The predictions in the $\tau =0$, $\gamma_0=0$ case agree well with the results 
of numerical calculations despite the fact that the approximate theory is applied 
for large deformations and for singularities that are near to the real axis, 
which is beyond its range of formal validity.  The agreement is suggestive and 
motivates us to use the localized approximation theory to examine singularity 
formation in the $\tau\not= 0$ case.

\subsection{Singularity Formation for $\tau \neq 0$}

Our main objective is to consider the case $g=0$, $\gamma_0\not= 0$, $\tau\not= 
0$; this is the simplest case in which $\tau\not= 0$ and which still exhibits 
unstable modes in the linear analysis (for sufficiently small $\tau$).  These 
parameter values correspond to a fluid interface which is subjected to a shear 
flow and, therefore, susceptible to the Kelvin-Helmholtz instability, but which 
is also acted on by the stabilizing influence of surface tension.  For $g=0$ 
the term on the right hand side of (\ref{eq:23}) is zero and consequently this equation 
can be directly integrated with respect to $\xi$.  After performing this 
integration and solving for $\Omega_\xi$, we obtain
\be
\Omega_\xi = \kappa u(\Omega) v(\Omega)
 \left\{ \sigma^2\Omega -{(A+\beta )\sigma^2\Omega^2 \over v(\Omega)^2} 
 +C\right\}
\label{eq:24.9}
\ee
where $\displaystyle \kappa={(\rho_1+\rho_2) \over \sqrt{2}\tau} (\Phi_0\Phi_0)^{3/2}$ and $C$ is 
a constant of integration.  Since $\Omega_\xi$ has zero mean (i.e., its constant 
Fourier mode is zero), the constant $C$ is chosen so that the sum contained 
within the  brackets in (\ref{eq:24.9}) has zero mean.  This gives
\be
C = -\sigma^2\Gamma_0+(A+\beta )\left[{\sigma^2\Gamma^2_0 \over (1+\beta
\Gamma_0\pm\sqrt{1+2\beta\Gamma_0+\Gamma^2_0}}\right]
\label{eq:25.5}
\ee
using definition (\ref{eq:22.5}).

We now examine the possibility of singularities in the traveling wave solution 
$\phi (\xi )$, $\psi (\xi )$, $\Omega (\xi )$.  The differential equation for 
$\Omega (\xi )$ is simplest when $\beta =-A$, so we  restrict ourselves 
to this case.  Analysis of (\ref{eq:24.9}) for $\beta\not= -A$ leads to similar conclusions
as for $\beta=-A$. 
With the simplification, (\ref{eq:24.9}) becomes
\be
\Omega_\xi=F(\Omega)     \label{eq:25}
\ee
where
\be
F(\Omega)=\kappa u(\Omega) v(\Omega) \{ \sigma^2 \Omega + C \}  \label{25.1}
\ee

We also require that $|\beta|=|A| \neq 1$,  i.e., the interface separates two fluids of
non-zero density. For certain values of 
$A$ we have been able to integrate (\ref{eq:25}) and obtain the solution for implicit 
form.  However, the implicit equation is complicated and it is simpler and 
perhaps more illuminating to analyze singularity formation without explicitly 
integrating.  Therefore, we proceed with this indirect approach.

The requirement that 
$\phi$ and $\psi$ have mean $1$ now makes it necessary to choose the $(-)$ branch in 
(\ref{eq:22},\ref{eq:23.5}).   
Since $\Omega$ is $2\pi$ periodic, upper 
analytic and has mean $\Gamma_0$, we can write $\Omega =\Gamma_0+
\widehat{\Omega}_1e^\xi+O(e^{2\xi})$ for constant $\widehat{\Omega}_1$ and $\xi
\to -\infty$.  Substituting this into equation (\ref{eq:25}) and equating coefficients 
of the leading order term $e^\xi$, gives
\[\widehat{\Omega}_1=\sigma^2 D\widehat{\Omega}_1\]
where $D=-\kappa(1+2\beta\Gamma_0+\Gamma^2_0)^{1/2}[1+\beta\gamma_0-(1+2\beta\Gamma_0
+\Gamma_0^2)^{1/2}]^{1/2}$ is the mean of the term $u(\Omega) v(\Omega)$ 
on the right hand side of (\ref{eq:25}).  The relation $\sigma^2 D=1$ determines 
$\sigma$, and a straightforward but lengthy computation verifies that the result 
is given by (\ref{eq:15.5}), as expected.  The solution to (\ref{eq:25}) is then uniquely 
determined by specifying $\displaystyle \widehat{\Omega}_1 = \lim\limits_{\xi\to -\infty}
{\Omega -\Gamma_0 \over e^\xi}$.  Each coefficient $\widehat{\Omega}_n$ for $n>1$ 
is determined from $\widehat{\Omega}_1,\dots ,\widehat{\Omega}_{n-1}$ and 
$\sigma$ as described earlier.

Singularities can occur in $\Omega (\xi )$ when $F(\Omega )=0$ or when $F(\Omega 
)=\infty$.  Corresponding  singularities
in $\phi$ and $\psi$  occur when $\Omega (\xi )$ 
is singular or when $1+2\beta\Omega +\Omega^2=0$ (this is also a situation in 
which $F(\Omega )=0$).  We now examine each of these possibilities:

\vspace{.2in}
\pist{1.}{$F(\Omega )=0$:  There are three situations which lead to $F(\Omega )
=0$.  One is  $u(\Omega)=0$, a second is $v(\Omega)
=0$, and a third 
is given by $\Omega -\Gamma_0=0$.
It can be shown that each of these situations does not 
produce a finite time singularity in $\Omega (\xi )$.  
In the first situation,  the function $\Omega(\xi)$ is
analytic in a 
neighborhood of a point $\xi_p$ where $u(\Omega)=0$.
To see this,  rewrite (\ref{eq:25}) as a differential equation in the function $u=(1+2\beta\Omega 
+\Omega^2)^{1/2}$.  This results in an equation of the form $u_\xi =q(u)$, where 
$q(u)$ is given by
\[q(u)=-[\kappa\sigma^2(\Omega (u)-\Gamma_0)][(\beta +\Omega (u))(1+\beta\Omega (u)
-u)^{1/2}]\,.\]
In the above, $\Omega (u)=-\beta +\sqrt{\beta^2-1+u^2}$; 
the branch of the square root is not important for our purposes.  The function 
$q(u)$ is analytic in a neighborhood of $u=0$.  Furthermore, $q(u=0)\not= 0$ for 
$\vert\beta\vert\not= 1$, since the relation $1+2\beta\Omega +\Omega^2=0$ 
guarantees that $\Omega (u=0)\not= -\beta$ and $1+\beta\Omega\not= 0$, and we 
assume that $\Omega (u=0)\not=\Gamma_0$ (the case $\Omega (u)=\Gamma_0$ is 
treated separately).  Consequently, $\displaystyle \xi =\int{1 \over q(u)}du=h(u)$ is an 
analytic function of $u$ in a neighborhood of $u=0$.
 Since $\displaystyle {dh \over du} (u=0)
\not= 0$, the function $\xi =h(u)$ can be locally inverted (for $u$ near $0$) to 
obtain $u=h^{-1}(\xi )$, which is analytic in a neighborhood of $\xi =h(0)$.  
%When $\vert\beta\vert =1$, the relation $1+2\beta\Omega +\Omega^2=0$ implies that 
%$1_\beta\Omega =0$, and the above arguments fail.  However, writing (~~~) as a 
%differential equation in the function $v=(1+\beta\Omega )^{1/2}$ leads to an 
%equation which behaves like $v_\xi\sim cv^2$ for constant $c$ and $v$ near zero.  
%This implies that $v=0$, or equivalently $1+2\beta\Omega +\Omega^2=0$, only for 
%infinite values of $\xi$.
Although $1+2\beta\Omega +\Omega^2=0$ can occur in finite time, the
function $\Omega(\xi)$ remains analytic in a neighborhood of the point where
this occurs.}

\vspace{.1in}
\pist{}{The second of the three situations does not occur
in finite time.  This is demonstrated by 
writing (\ref{eq:25}) as a differential equation in the variable $v=(1+\beta\Omega +
u(\Omega))^{\onehalf}$,  which leads to an equation  
that
behaves like 
$v_\xi \sim cv$  for constant $c$ and $v$ near zero.
This implies that 
 $v =0$, or equivalently $1+\beta\Omega +
u(\Omega)=0$, only for infinite values of $\xi$.} 

\vspace{.1in}
\pist{}{Finally, we consider the possibility that $\Omega -\Gamma_0=0$.  Writing (\ref{eq:25}) 
as a differential equation in the function $\Omega '=\Omega -\Gamma_0$ leads to 
an equation which has the form $\Omega '_\xi =c_1\Omega '$ for $\Omega '$ near 
$0$,  where $c_1$ is a constant.  Again we conclude that $\Omega '=0$ only for 
infinite values of $\xi$.}

\vspace{.1in}
\pist{}{Singularities can still at points where $F(\Omega )=\infty$.  We consider this 
possibility next.}

\vspace{.2in}
\pist{2.}{$F(\Omega )=\infty$:  This can only occur if $\Omega\to\infty$.  For 
$\Omega$ near $\infty$,
\[ F(\Omega )\sim d~{\Omega^{5/2}}\]
where $d$ is a constant which depends on the particular choice of branch in (\ref{eq:25}) 
but which is not zero.  Thus equation (\ref{eq:25}) allows singularities of the form $\Omega 
\sim d_1(\xi -\xi_0)^{-2/3}$, and leads to singular behavior 
in $z_\xi$ given by $z_\xi\sim d_2(\xi -\xi_0)^{-2/3}$,  where $d_1$
and $d_2$ are constants.}

\vspace{.2in}

It is necessary to verify that the singular behavior in $\Omega$ actually 
corresponds to a geometric singularity in the interfacial curve and does not 
just indicate a singularity in the parameterization used to represent the curve.  
This is most easily done by computing the sheet curvature $\displaystyle {\cal K}=\Im\left[
{z_\alpha z^\ast_{\alpha\alpha} \over (z_\alpha z^\ast_\alpha)^{3/2}}\right]$ 
near a point $\xi_0$ where $\Omega (\xi_0)=\infty$.  Near such a point we find 
that ${\cal K}\sim b~ \Omega^{1/2}$, which indicates that the curvature has a 
singularity of the form ${\cal K}\sim b_1(\xi -\xi_0)^{-1/3}$, for constant 
$b,~b_1$.  The unnormalized vortex sheet strength $\gamma =\sigma\Phi_0\Psi_0(\Omega 
+\Omega^\ast )$ blows up at $\xi_0$; however the true vortex sheet strength
$\displaystyle \widetilde{\gamma} = {\gamma \over \vert z_\alpha\vert}$ remains bounded.

In summary, for $A\not= 1$ and $\beta =-A$ singularities arise in the traveling 
wave solution at positions in the finite
complex plane  where $\Omega =\infty$.  The singularities correspond to 
infinite curvature at the interface.  Since $\Omega (\xi )$ is an upper analytic 
function, these singularities are located in the lower half of the complex 
$\alpha$ plane.  The theory predicts that they move toward the real  
axis with constant speed $\sigma$, as long as
$\tau < (1-A^2)\gamma^2/4$ (required for $\sigma$ to have
a non-zero imaginary part).  When they reach the real axis, the 
singularities appear on the sheet.

We remark that our results are consistent with those obtained by Baker, Cowley, 
Tanveer and Page \cite{co:asy}.  They employ a local analysis similar to that which led 
to the system of PDE's (\ref{eq:13}), and show that curvature singularities in the 
complex plane of the form ${\cal K}\sim c(\xi -\xi_0)^{-1/3}$ are allowed by 
their localized equation.

The information which we now have concerning the analytic behavior of $\Omega 
(\xi )$ near a point  where $u(\Omega)=\sqrt{1+2\beta\Omega +\Omega^2}=0$ forces us 
to be more precise in defining a branch for the square root.  Let the curve 
which is traced by $1+2\beta\Omega +\Omega^2$ as $\alpha$ ranges from $0$ to 
$2\pi$ be denoted by ${\cal C}$.  For $\Omega '=\Omega -\Gamma_0$ small in 
the range $\alpha=0$ to $2\pi$, i.e. for $\widehat{\Omega}_1{e}^{\sigma t}$ 
small, the curve ${\cal C}$ will not cross the negative real axis.   


As $\widehat{\Omega}_1{e}^{\sigma t}$ (and therefore $\Omega '$) increases, 
the curve ${\cal C}$ eventually goes through zero; however, we have seen that 
$\sqrt{1+2\beta\Omega +\Omega^2}$ is not singular when this occurs.  For larger 
values of $\widehat{\Omega}_1{e}^{\sigma t}$ the curve ${\cal C}$ crosses the 
negative real axis, but remains an analytic function of $\xi$  (see Figure 1).  
We define a branch of the square root which maintains analyticity of $(1+2\beta 
\Omega +\Omega^2)^{1/2}$ as a function of $\xi$ when ${\cal C}$ crosses the 
negative real axis, in the following manner.  Define $1+2\beta\Omega +\Omega^2=M
e^{i\theta}$ where $M$ is the amplitude and $\theta$ is the angle of $1+2\beta
\Omega +\Omega^2$ with respect to the $x$-axis.  The angle $\theta =0$ is 
defined so as to correspond to the right most crossing of the curve ${\cal C}$ 
with the $x$-axis, and the other values of $\theta$ are fixed by the requirement 
that $\theta$ be a continuous function of the arclength of ${\cal C}$.  For 
example, in Figure 1, point $A$ corresponds to $\theta =0$, point $B$ to an 
angle $\theta$ just greater than $\pi$, point $C$ to $\theta =2\pi$, and 
point $D$ to an angle just greater than $3\pi$.  The square root is then defined 
by $(1+2\beta\Omega +\Omega^2)^{1/2}=\sqrt{M}e^{i\theta /2}$. In the case
$\gamma_0 \neq 0,~\tau \neq 0,~g=0$,  we keep the minus 
sign in front of the square roots in (\ref{eq:23.5}) so that our branch 
definition will satisfy the requirement that $\phi$ and $\psi$ have mean $1$ 
when $\widehat{\Omega}_1{e}^{\sigma t}$ is small.  For larger values of 
$\widehat{\Omega}_1{e}^{\sigma t}$, such as in Figure 1, our branch definition 
is equivalent to using the principal branch with a minus sign in front for the 
larger circular lobe, and a plus sign in front for the smaller lobe.  Note that 
this definition of branch gives continuous dependence of the solution $\Omega 
(\xi )$ on $\widehat{\Omega}_1{e}^{\sigma t}$.  We similarly define a branch 
for the term $v(\Omega)=(1+\beta\Omega +u(\Omega ))^{1/2}$. 
It was shown earlier that the trace of $1+\beta\Omega 
+u(\Omega )$ does not pass through zero as $\widehat{\Omega}_1{e}^{\sigma t}$ 
is increased; however, it can cross the negative real axis.   
In this case, the square root definition given above maintains continuity of the 
term $(1+\beta\Omega +u(\Omega ))^{1/2}$.

The location  $\xi_s$ of a  singularity 
 in the complex plane can be determined by integration
of (\ref{eq:25}) with the requirement that $\Omega(\xi) \sim \widehat{\Omega}_1 e^{\xi}$
as $\xi \rightarrow \infty$.   It is easy to show (see \cite{ba:sin}) that 
\be 
\xi_s = \lim_{L\rightarrow\infty} \int^L_0 \left(\frac{1}{F(\Omega\prime)} -
\frac{1}{\Omega^\prime} \right) d\Omega^\prime
 + \ln L  -\ln\widehat{\Omega}_1.  \label{eq:25.2}
        \ee
This expression defines two singularity positions, one for each sign
of (principal) branch of $F$.
If the sum of the first two terms in (\ref{eq:25.2}) is denoted by $p$, then the
trajectory of the singularity $\xi_p$ is given by
\be
\alpha = - i(p - \ln \widehat{\Omega}_1 ) -\sigma t\; . \nonumber
        \ee
Thus, for a given value of surface tension $\tau$,
the distance of the singularity from the real line and the time it takes
to reach it depend logarithmically  on the size of the perturbation,
for small perturbations.  Although it is easy
to determine $p$ numerically, we instead solve (\ref{eq:25}) for $\Omega(\xi)$ so
that the sheet location may be obtained in addition to the time of singularity formation.


\section{Numerical Solutions}
\subsection{Numerical Method}

We have obtained numerical solutions to equation (\ref{eq:25}) in the case $\beta =A=0$ 
for various values of the surface tension coefficient $\tau$.  These solutions can 
be used to obtain the shape of the interface,
location of singularity formation,  and various other quantities of 
interest.
In the computations, we   set $\gamma_0=1$ and $\rho_1+\rho_2=1$;
together with a relevant length scale, such as the
length of the computational box $L=2 \pi$, this
choice makes our variables effectively dimensionless.

 Input data for the computed solution consists of the surface tension 
coefficient $\tau$ and either the Fourier coefficient $\displaystyle 
\widehat{\Omega}_1=\lim\limits_{\xi\to -\infty}{\Omega -\Gamma_0 \over e^\xi}$ 
or the initial value $\Omega (\xi_0)$ of the function at the point $\xi_0=i
\alpha_0$.  Equation (\ref{eq:25}) is then used to determine the higher Fourier coefficients 
$\widehat{\Omega}_2,\widehat{\Omega}_3,\dots$ in the former case or integrated 
with respect to $\alpha$ to obtain the full solution in the latter case.  Since 
$\Omega$ is of traveling wave form, its values at some later time $t$ can be 
found by multiplying each Fourier coefficient $\widehat{\Omega}_k$ with 
$e^{k\sigma t}$.  In practice, however, this procedure causes the small round 
off errors which exist in the high Fourier modes to be greatly magnified.  As a 
result, we follow the second approach above by giving the value of $\Omega (
\alpha =0)$ and integrating equation (\ref{eq:25}) with respect to $\alpha$.  
The solution at a later time $t$ is then found by changing  
$\Omega (0)$ to a different value  $\Omega_1(0)$
and again integrating.  The time $t$ of the solution 
passing through $\alpha =0$, $\Omega_1 (0)$ is given by
\[t={1 \over \sigma}\ln {\widetilde{\Omega}_1 \over \widetilde{\Omega}}\,,\]
where $\widetilde{\Omega}_1$, $\widetilde{\Omega}$ are the Fourier coefficients 
of $e^{i\alpha}$ in the solutions with initial values $\Omega_1(0)$, $\Omega
(0)$, respectively.

Equation (\ref{eq:25}) was integrated with respect to $\alpha $ using the SLATEC 
differential equation solver DDEABM.  The routine employs an
Adams-Moulton predictor-corrector method with adaptive step size.  The step-size 
is chosen so that the single step error does not exceed a user specified error 
per step.   In each plot which is presented, a sufficient
number of points are used to resolve all features. 




Once the function $\Omega (\xi )$ has been determined, the other two variables 
in the localized approximation equations, $\phi$ and $\psi$, can be determined 
from (\ref{eq:22}).  Integrating the relation $z_\alpha =\phi +\psi^\ast -1$ then gives 
the position of the interface.
The $x$ constant of integration is determined
by the requirement that $x(0)=0$.  It is also required
that the physical interface $ y(x)$ has no mean height; this fixes the
$y$ constant of integration.    The unnormalized vortex sheet strength is 
determined from $\gamma=\gamma_0+2Re~\omega$ with $\omega =\sigma\Phi_0\Psi_0
(\Omega -\Gamma_0)$.

\subsection{Computational Results}

\renewcommand{\thefootnote}{\fnsymbol{footnote}}
\addtocounter{footnote}{1}
The solution $\Omega (\xi )$ to (\ref{eq:25}) is 
plotted in Figure 2 at four different times, 
for $\tau =.10$
The curves shown are 
obtained by plotting $Re~(\Omega -\Gamma_0)$ on
the $x$-axis and $Im~ (\Omega -\Gamma_0)$ on the $y$-axis, for $\alpha$ in the
interval $[0,2\pi ]$.  
The corresponding interface positions and vortex sheet 
strengths are plotted in Figures 3 and 4.  For this value of $\tau$ the linearized 
theory gives a cut-off wave number $k_c=2.45$.  Consequently, in the discrete 
system there are two unstable modes corresponding to $k=1$ and $k=2$, with the 
remaining modes neutrally stable.  This case is intermediate between two extremes, 
namely that for which $k_c\to\infty$ as $\tau\to 0$ and that for which $k_c=1$ 
and $\tau =1/4$.  These two extremes will be discussed subsequently.

In  Figure 2a, the initial value of $\Omega$ is chosen as 
$\Omega (0)-\Gamma_0 =-.001$.  By choosing the value of $\Omega -\Gamma_0$ to 
be a negative real number when $\alpha =0$, we obtain an interfacial curve 
which satisfies $z(\alpha =0)=0$ and which has odd symmetry with respect to 
$\alpha =0$.   
The value of $\Omega (0)-\Gamma_0$  in Figure 2a is small enough so that 
the motion is in the linear regime.  Consequently, the solution is nearly 
circular and is closely approximated by $\Omega (\xi )=\Omega (0)e^\xi$, the 
solution to the linearized version of (\ref{eq:25}).  As $\Omega (0)$  decreases, i.e., $\vert\Omega 
(0)\vert$ increases (this is equivalent to letting the solution $\Omega (\xi)$ 
with initial data given in Figure 2a evolve forward in time) the 
motion enters the non-linear regime and the curve $(Re~ \Omega(\alpha),~Im ~\Omega(\alpha) )$ deforms 
from its previously circular form.  In Figure 2b, we see that a dimple forms 
on the curve as it tends to avoid the point $1+\Omega^2=0$, which corresponds 
to a zero of the right hand side of equation (\ref{eq:25}). 
The other
zeros of $F(\Omega)$ are farther away and do not visibly affect the shape of the curve.
This zero and its consequent 
dimple lead to the formation of two 
relative vorticity maxima (see Figure 4b).
This contrasts the zero surface tension case, in which a single
maxima forms at $\alpha=0$.  Physically,  
surface tension produces a restoring force which prevents  
accumulation of vorticity at the center of the interface, where the curvature 
gradients are large.  Instead, two relative vorticity maxima form on either side 
of the centerpoint, where the curvature gradient is smaller \cite{ra:non}.

As $\Omega (0)$ is further decreased, the curve $(Re~\Omega (\alpha ),~Im~\Omega (
\alpha ))$ eventually passes through the
point where $1+\Omega^2=0$ (Figure 2c).  
At this point $\Omega_\alpha =0$ and the curve $(Re~\Omega (\alpha ),Im~\Omega 
(\alpha ))$ achieves a zero angled cusp.  However, as demonstrated earlier, this 
point does not correspond to a singularity in any of the functions $\Omega (\xi
)$, $\phi (\xi )$ or $\psi (\xi )$, and the evolution may be smoothly continued 
beyond it.
The profile of the interface at this 
time indicates that the wave has just begun overturning.

Further  decreasing 
$\Omega (0)$ produces a curve $(Re~\Omega (\alpha ),~Im~\Omega (\alpha ))$ like 
that shown in Figure 2d.  At this time the interface, shown in Figure 3d,  has 
already developed the appearence of corners  at the turning points 
in the ${\cal S}$-shaped curve.  
The sheet strength (Figure 4d) and curvature (Figure 5d) peak sharply at values 
of $\alpha$ corresponding to the location of the `corners'. 
Between these peaks, the sheet strength is negative.
Shortly after this time, 
the two intersecting lobes of the curve $(Re 
~\Omega (\alpha ),Im~\Omega (\alpha ))$ stretch
out to infinity on the right.   
At this point, singularities appear on the 
interface and the unnormalized vortex sheet strength and curvature both 
blow-up.

Plots of the curvature corresponding to the interfacial
profiles of Figure 3 are presented in Figure 5.
Note that at intermediate times in the
motion, a localized `sine' like variation in the
curvature develops and grows.  
At late times, two extremely sharp peaks of curvature 
shoot up from the edges of this localized perturbation,
in the same locations as the vorticity peaks.
The  initial evolution of these eventual maxima are more closely
detailed in Figure 6.  Interestingly, the original
perturbation shrinks as the sharp peaks grow.  For other values of $\tau$
the maxima  are located in different positions (this is detailed below)
and the perturbed region stretches or compresses accordingly;
otherwise the curvature resembles that in Figure 5 and
will not be shown.


Many aspects of the behavior remain the same when the surface tension is 
increased, subject to the constraint that there is at least one unstable mode.  
Figures 7 and 8 depict the traveling wave solution for $\tau =.219$ and $k_c=1.1$, 
at stages of the evolution analogous
to those shown in Figures 3 and 4. 
Despite the similarities with the previous figures, there are 
some important qualitative differences.   
As shown in Figure 8, the two vorticity maxima 
are further separated in $\alpha$.  The 
interfacial profile assumes more of an
extended ${\cal S}$ shape, and the amplitude of 
the deformation at the critical time is significantly larger than for the 
previous value of $\tau$.

When surface tension is decreased, the opposite trends are observed.  For example, 
consider $\tau =.0079$, for which the cut-off wave number is $k_c=31.6$.  Figures 
9 and 10 contain the pertinent plots for this value of $\tau$.   
In this case, the two vorticity and curvature maxima are 
closer in $\alpha$ than in the $\tau =.10$ case.
Furthermore, the 
interfacial profile assumes a `triangular' shape, and the amplitude at the 
critical time is greatly decreased from the previous cases.  These trends 
continue as $\tau$ is further reduced:  the separation of the vorticity and 
curvature maxima further decreases (see Figure 11), 
the sheet profile approaches a limiting 
`triangular' shape at the critical time, and the size of the deformation at this 
time decreases.  Figure 11  further illustrates the dependence of singularity position
on surface tension.
When 
$\tau =0$, the only traveling wave solution is the flat sheet solution (\ref{eq:5}).
The form of the traveling wave solution at the critical time for decreasing
values of $\tau$ suggests an approach (in the weak sense) to the flat sheet solution
as $\tau \rightarrow 0$.


\section{Discussion}

An interesting feature of the traveling wave solutions is that,
for each
value of $\tau$, singularities form before there is any sign of interface
roll-up.   This is inconsistent with the results of  \cite{ra:non}
and \cite{ho:rem}, which show rolled-up interfaces in the presence
of surface tension.\footnote{
Our results are consistent with those of Pullen \cite{pu:nst}, in
which the interface develops irregularities prior to roll-up similar to those
produced by a singularity; however, the results of Baker and Nachbin \cite{ba:sta}
strongly suggest that the irregularities in Pullin's computations  are numerical in nature.}
There
are several possible explanations for this discrepency:
 
\vspace{.2in}
\pist{1.}{Interactions between positive and
negative wavenumber modes, absent in the approximate equations, are necessary
to produce roll-up.}
 
\vspace{.12in}
\pist{2.}{The difference in initial data, i.e., traveling wave initial data
versus linearized eigenfunction initial data in \cite{ra:non} and 
$\tau=0$ linearized eigenfunction initial data in \cite{ho:rem}, are enough to 
cause the discrepency.  We note that the traveling wave initial data
and the linearized eigenfunction initial data are identical to $O(\varepsilon^2)$,
where $\varepsilon$ is the height of the initial disturbance; the traveling 
wave and $\tau=0$ linearized eigenfunction data differ by $O(\varepsilon)$.}

\vspace{.12in}
\pist{3.}{Smoothing operations remove
a singularity which occurs prior to roll-up and enable the computations to
continue into the roll-up regime.} 

\vspace{.2in}
The careful computations of Hou et. al. \cite{ho:rem}
show roll-up occuring prior to the formation
of singularities  and therefore render
1 and 2  more plausible explanations of our lack of roll-up.
The discrepency may be related to the  late time vorticity distribution  
for  the traveling wave solution, which
contains two very narrow peaks  with a region of oppositely signed vorticity between them.
This distribution is consistent with an interfacial profile that stretches
but does not roll-up.
The vorticity exhibits this form regardless of the 
value of surface tension, as long as there
is at least one unstable mode.
However, this particular distribution may not
be generic for other initial data.

  For $\tau \stackrel{>}{\sim} .04$ ,  the computations in
\cite{ra:non} also show a doubly peaked vorticity distribution,  although
$\gamma$ remains positive for all $\alpha$.  In this situation, they find that the sheet
does not roll up.   On the other hand, for smaller $\tau$,
they find the vorticity accumulates in a single peak
at the center of the curve, leading to roll up.
It is possible that
the method of \cite{ra:non} does not differentiate between
two closely spaced vorticity peaks, since their method
involves insertion and redistribution
of vortex
sheet markers, using linear interpolation to position and define the sheet strength of
the redistributed markers.  This is an effective smoothing operation.

In contrast,
 Beale and Lowengrub \cite{be:wel}
find that the vorticity distribution always forms two peaks
and that the separation of these peaks decreases with decreasing $\tau$.
Their results  suggest that these two features are generic.
They use  the $\tau=0$ linearized eigenfunctions for initial data, and
consider values of surface tension for which their code runs
a reasonable length of time.  The presence of two vorticty peaks
was taken as evidence of the existence of two singularities
in the complex plane which approach each other as $\tau$ decreases.
Although their vorticty distribution strongly resembles
that for the traveling wave  solution, their code stops before the vorticity between the
peaks goes negative.  The computations  of Hou et. al. \cite{ho:rem}
extend the solution for this initial data to much longer
times.  However, their data does not show
the vorticity going negative until just before a late
time singularity, at which point the (rolled-up) interface and
the vorticity distribution have little resemblance to the traveling
wave solution. Nevertheless, the exponent of curvature blow-up is
the same as that predicted by the traveling wave solution.  (need to talk
more with Mike or John about this). 



\section{Conclusion}

The analysis and computations presented above  show that,
in the approximate theory, singularities develop on the fluid interface in 
finite time.  The behavior of the sheet position $z$, the vorticity
strength $\gamma$, and the curvature $\kappa$ 
in the vicinity of  a singularity $\zeta_s$ is given by
$z_\xi\sim a_1(\xi-\xi_s)^{-\frac{2}{3}},~ \gamma\sim a_2(\xi-\xi_s)^{-\frac{2}{3}}$,
 and $\kappa\sim a_2(\xi-\xi_s)^{-\frac{1}{3}}$, where $a_1$--$a_3$ are constants.
The theory is applied outside of its range of formal
validity, so the finite time  occurence of these singularities in the full equations must be 
accepted as merely plausible. 
Numerical computations 
using the traveling wave initial conditions are being performed 
by G. Baker and A. Nachbin in order to verify the
asymptotic predictions.


\clearpage
\bibliography{ref}

\pagebreak

\begin{center}
\large{\bf List of Figure Captions}
\end{center}

\begin{enumerate}
\item
The form of the curve ${\cal C}$, after it has passed through the origin.

\item
Phase plots of the solution $\Omega(\alpha)-\Gamma_0$ to 
(\ref{eq:25}), for $\alpha \in [0,2\pi]$.  Surface tension
is set at $\tau=.10$.  The plots (a)--(d) correspond to initial
values $\Omega(0)-\Gamma_0=$ -.001, -.15, -.34, -.569 and times
$t=$ 0., 11.02, 12.25, and  13.15.  The marker ($\times$) denotes
the position where $1+\Omega^2=0$.

\item
The interface profiles  $(x(\alpha),y(\alpha))$ corresponding
to the solutions in Figure 2.  Note that the $t=0$ curve
is plotted on a different scale.

\item
The vortex sheet strength $\gamma$ versus $\alpha$,
corresponding to the solutions in Figure 2.  The plots are on
different scales.
In (d), the sheet strength peaks at $\gamma(\pm.87)=173$ (not shown);
the strength at $\alpha=0$ is $\gamma(0)=-.14$.

\item
The curvature $\kappa$ versus  $\alpha$, corresponding to the 
solutions in Figure 2.  The plots are on different scales.
In (d), the curvature peaks at $\kappa(\pm.87)=\pm26.5$ (not shown).

\item
The initial evolution of curvature  spikes.
The plots correspond to initial values $\Omega(0)-\Gamma_0=$
-.54, -.55, -.56 and times $t=13.05, 13.08, 13.12$.

\item
The interface profiles $(x(\alpha), y(\alpha))$ for 
$\tau=.219$.   The curves (a)--(d) correspond to initial values 
$\Omega(0)-\Gamma_0=$ -.001, -.4, -.47, -.5402, and times
$t=$ 0., 35.7, 37.1, 38.5.  The $t=0$ curve is plotted
on a different scale.

\item
The vortex sheet strength $\gamma$ versus $\alpha$,
corresponding to the profiles in Figure 7.  The plots are on
different scales.
In (d), the sheet strength peaks at $\gamma(\pm2.31)=303$ (not shown);
the strength at $\alpha=0$ is $\gamma(0)=-.08$.

\item
The interface profiles $(x(\alpha), y(\alpha))$ for
$\tau=.0079$.   The curves (a)--(d) correspond to initial values
$\Omega(0)-\Gamma_0=$ -.00001, -.001, -.175, -.5718  and times
$t=$ 0., 6.7, 8.1, 8.2.  The $t=0$ curve is plotted
on a different scale.

\item 
The vortex sheet strength $\gamma$ versus $\alpha$, 
corresponding to the profiles in Figure 9.  The plots are on 
different scales. 
In (d), the sheet strength peaks at $\gamma(\pm.06)=145$ (not shown); 
the strength at $\alpha=0$ is $\gamma(0)=-.14$.

\item
The distance from the origin  of the singularity on the interface,
measured in terms of the arclength $s$ along the interface,
and plotted as a function of $\tau$.   For the range of $\tau$ shown,
the dependence is approximately linear.

\end{enumerate}


\end{document}
