text:
\section{Introduction}

Accretion onto astrophysical objects possessing a material surface
(as opposed to accretion onto black holes) always involves a non-trivial
interaction of the incoming gas with the object's outer layers.
Examples of such objects include white dwarfs in cataclysmic
variables (CVs), young stars gaining material from a protoplanetary
disk, and accreting neutron stars. In all these systems, one must
understand how the incoming gas shares its angular momentum with
the accreting object and how it mixes with the previously accreted
material.

If the magnetic field of the central object is strong enough, it can
disrupt the disk at some distance above the surface. Inside this region,
accretion proceeds along field lines
\citep{GhoshLamb,Koldobaetal}. The
critical value of the magnetic field at the surface of the central object
for magnetic disruption is given by
\begin{eqnarray}
B_{*,\text{crit}} = 5 \times 10^4 \left(\frac{\beta}{.5} \right)^{-7/4}
\left(\frac{M_*}{.6 M_\odot}\right)^{1/4} \left(\frac{\dot{M}}{10^{-8}
  M_\odot \ \text{yr}^{-1}}\right)^{1/2}\left(\frac{R_*}{9 \times 10^3
  \ \text{km}}\right)^{-5/4} \text{G}.
\label{Bdisrupt}
\end{eqnarray}
Here, we have used parameters for the mass and radius that are typical
for CVs in outburst \citep{Bergeronetal,Hachisu,Patterson}. The
parameter $\beta$ is a dimensionless factor of order unity that depends on
the model adopted for the disruption of the disk by the stellar
magnetic field \citep{GhoshLamb}.

For $B_* < B_{*,\text{crit}}$, the accretion disk will
extend all the way to the surface of the star, for which there
is good observational evidence in neutron star low mass X-ray binaries
\citep{GilfanovRevnivtsevMolkov} and in dwarf nova systems in outburst
\citep{WheatleyMaucheMattei}. Since the rotation rate
of the star, $\Omega_*$, is slower than the Keplerian rotation rate at
the stellar surface, $\Omega_K(R_*)$, a boundary layer (BL) will exist
inside of which $d \Omega / d R > 0$ and the rotation profile of
the star attaches smoothly to that of the disk. If the accretion
rate is high enough, it is also possible for material to spread
meridionally to high latitudes forming a spreading layer (SL)
\citep{InogamovSunyaev, PiroBildsten}. For both a BL and a SL, the intense
energy release in a localized region near the stellar surface leads to
easily observable signatures such as hard spectral components,
variability of emission, and so on.

A considerable amount of effort has been previously devoted to
understanding the structure of well-developed, steady-state BLs,
which have been evolving for long enough to establish a smooth
rotation profile in their interior. An outstanding question in the
study of steady-state BLs is
identifying the mechanism of angular momentum transport in
the layer. Several potential mechanisms have been explored over the
years, among them shear instabilities \citep{KippenhahnThomas},
baroclinic instabilities \citep{Fujimoto, Hanawa}, and the
Tayler-Spruit dynamo
\citep{Spruit, PiroBildsten1}. Despite this effort, no clear answer exists at
present regarding the nature of the angular momentum transport and
mixing in well-developed BLs.

\subsection{Initiation of Boundary Layers}

An interesting aspect of the BL problem that has received less
attention in the past is the issue of BL {\it initiation}, i.e.
the initial stage of BL formation which must occur when the accreted
material first touches the surface of the star. In this case the physical
setup is going to be quite different from the steady-state BL
primarily because of the much larger velocity shear, justifying
the study of BL initiation in the limit of an
essentially discontinuous rotation profile at the interface between the
star and accreting material.

Even though the initiation stage is just a transient phase in the BL
evolution, one still needs to understand it to get a full picture of
the BL phenomenon. Also, one need not think that the BL initiation
is a unique event for every accreting object --- it may, in fact, be
repetitive. The most common situation in which this BL initiation is
recurrent involves the dramatic increase of the mass accretion
rate (due to some sort of outburst triggered by an instability)
through the disk which is magnetically disrupted outside the star
under normal circumstances. According to Equation (\ref{Bdisrupt}), a
sudden increase of $\dot M$ by several orders of magnitude as typical
for some accreting systems (e.g. dwarf novae, FU Ori outbursts, etc.)
can rapidly compress the stellar magnetosphere to the point at which
disk material starts touching the stellar surface and a BL starts
to form.

There is good observational evidence for this kind of recurrent
behavior. For instance, \citet{LivioPringle} have
argued that the observed lag in the rise time of the UV emission
relative to the optical in a CV system transitioning to outburst can
be explained if the disk is magnetically disrupted in quiescence but
not in outburst. Thus, the hottest, innermost part of the disk is
evacuated in quiescence, and UV radiation is not emitted immediately in
the transition to outburst, since the empty region requires time to
become filled.

FU Ori stars are another type of system in which the
disk can extend all the way to the stellar surface, due to the
high accretion rate - $\dot{M} \sim 10^{-4} M_\sun
\text{yr}^{-1}$ \citep{KHH}.
In these systems, the BL can puff up to become
of order the stellar radius and the distinction between the boundary
layer, the disk, and the star becomes blurred \citep{PNFuOri}.

The goal of our present work is to make a first step towards
understanding the {\it formation} of the BL. As in the case of a
well-developed BL, the major issue for this initial stage lies
in details of the angular momentum transport and mixing. However,
because of the enormous shear present at the star-disk interface
at this stage, it is highly likely that purely hydrodynamical shear
instabilities would dominate the transport of mass and momentum
rather than anything else.

For that reason, in this work we primarily focus on exploring the
properties of various shear-driven instabilities under the conditions
typical at the BL formation stage. {First, we seek to identify a
particular variety of the shear instability that most efficiently
initiates mixing between the two fluids, i.e. has the fastest growth
rate. Second, we examine the conditions needed for its operation,
such as the density contrast between fluids, initial velocity
profile, etc.}

There are several physical
ingredients which can potentially affect operation of shear
instabilities: stratification, rotation, magnetic fields, turbulence,
radiative transfer, and the supersonic nature of the flow. The latter
aspect of the problem is very important and is inevitable during
the BL initiation phase, when disk material rotating at the
Keplerian speed comes into contact with the more slowly spinning stellar
surface. The differential azimuthal velocity between the two
interacting flows is then bound to be a significant fraction of the
Keplerian velocity at $R_*$ and should highly exceed the
sound speed both in the disk and in the outer layers of the star.

\subsection{Shear Instabilities in Compressible Fluids}

The study of shear instabilities in compressible fluids is of
fundamental physical significance and has received much attention in the
past. \citet{Landau}, \citet{Hatanaka}, \citet{Pai}, \citet{MilesKH},
and \citet{Gerwin} have all studied the problem in the vortex sheet
approximation, where one half plane of compressible
fluid moves at constant velocity over another. When the two fluids
move at a low Mach
number relative to each other, one finds that infinitesimal
perturbations are governed by the
classical Kelvin-Helmholtz (KH) dispersion relation for two incompressible
fluids. However, \citet{MilesKH} showed that above a critical Mach
number, the vortex sheet
becomes marginally stable to two-dimensional disturbances along the
direction of the flow. This is a surprising result considering that
one might have expected increasing the shear to lead to
an increased growth rate for the instability rather than stabilization.
Understanding the implications of this result for the momentum and
mass transport in the astrophysical BLs is one of the goals of our
study.

Later, \citet{Blumenetal}, \citet{Ray}, \citet{ChoudhuryLovelace} and
\citet{Glatzel} studied
the stability of a fluid with a continuous and monotonically varying
velocity profile. They found that unlike the vortex sheet, continuous
velocity profiles with a supersonic velocity difference across them
were unstable even at high Mach number. \citet{Glatzel} showed that
the instability was similar to the
Papaloizou-Pringle instability which operates in hydrodynamical disks with
radial boundaries
\citep{PP}. Thus, for high Mach number flow in a compressible fluid, a
finite thickness velocity
profile exhibits fundamentally different behavior from the vortex
sheet, since it can support unstable modes. Moreover, the growth rate
of the unstable modes scales as $\propto 1/L$, where $L$ is the
thickness of the shear layer. Thus, the thinner the shear layer, the
faster the instability operates! This is exactly the opposite of what
one might have expected from the vortex sheet stability criterion and
we will provide an explanation for this in \S \ref{finite_width_SL}.

Previous studies of possible instabilities inside astrophysical BLs have
been primarily concerned with the sub-sonic regime when the flow
can be considered as almost incompressible \citep{KippenhahnThomas,
  Fujimoto1}. This
regime may indeed apply in well-developed BLs with smooth shear
profiles, even though at the present level of knowledge one can
hardly exclude the possibility of the existence of localized regions
in steady-state BLs where compressibility is still important.
Our primary goal here is to extend these studies into the
regime of highly compressible, supersonic flows and to explore the
implications for astrophysical objects. In the course
of studying supersonic shear instabilities, we will sometimes
also incorporate stratification in our calculations and explore
the role of rotation.

The paper is organized as follows. In \S \ref{goveqn} we present our
governing set of equations and describe
the formalism we use to study the stability problem.
Then in \S \ref{compress_KH}
we study the case of the compressible vortex sheet with no gravity
and show that our formalism reproduces the dispersion relation
obtained by previous authors \citep{Landau, Hatanaka, Pai, MilesKH,
Gerwin}. In \S
\ref{compress_grav}, we introduce gravity as a small parameter and
perform a perturbation expansion to study what influence this has on
the stability of the vortex sheet. Finally, in \S \ref{finite_width_SL} we
study the stability of a finite width shear layer with a linear velocity
profile. We derive new dispersion relations for this case and test
the growth rate of the fastest growing mode using the
Godunov code Athena \citep{Stoneetal}.

\section{Formalism}
\label{goveqn}

Here we discuss the formalism that we use for analytic calculations throughout
the paper. We adopt cylindrical
coordinates $(\varpi, \theta, z)$ and for simplicity assume that all
equilibrium quantities are independent of $z$. Thus, for our purposes
gravity is given by $-g(\varpi) \hat{\boldsymbol \varpi}$ ($g(\varpi) > 0$), and we ignore the vertical
stratification of the disk. Such a setup does not allow for the
baroclinic instability, which would require $d \Omega/dz \ne 0$, but it does
contain the necessary ingredients for studying shear
instabilities. Given our assumptions, the equation of hydrostatic
equilibrium reads
\begin{eqnarray}
\bar{g}(\varpi)\equiv g(\varpi) - \varpi \Omega^2(\varpi)= -\frac{1}{\rho} \frac{d P}{d \varpi},
\end{eqnarray}
where $\bar{g}(\varpi)$ is the effective gravity. The equilibrium density,
pressure, and sound speed are given by $\rho(\varpi)$, $P(\varpi)$,
and $s(\varpi)$ respectively.

 Denoting the $(\varpi, \theta, z)$ velocities by
$(u, v, w)$, using $v = \Omega \varpi$, and assuming
 adiabaticity, the Euler equations in cylindrical coordinates are
\begin{eqnarray}
\label{Eulerstart}
-\frac{\partial \rho}{\partial t} &=&\frac{1}{\varpi}
\frac{\partial}{\partial \varpi}(\rho \varpi u) +
\frac{1}{\varpi}\frac{\partial}{\partial \theta}(\rho v) +
\frac{\partial}{\partial z} (\rho w) \\
\frac{D u}{D t} &=& -\frac{1}{\rho}\frac{\partial
  P}{\partial \varpi} + \frac{v^2}{\varpi} - g \\
\frac{D v}{D t} &=& -\frac{1}{\rho \varpi} \frac{\partial
  P}{\partial \theta} - \frac{u v}{\varpi} \\
\frac{D w}{D t} &=& -\frac{1}{\rho}\frac{\partial P}{\partial z} \\
\label{Eulerend}
\frac{D \left(P\rho^{-\gamma}\right)}{D t}  &=& 0 \\
\frac{D}{D t} &\equiv& \frac{\partial}{\partial t} + u
\frac{\partial}{\partial \varpi} + \frac{v}{\varpi} \frac{\partial}{\partial
  \theta} + w \frac{\partial}{\partial z}
\end{eqnarray}
On top of the zeroth order equilibrium state, we consider infinitesimal
two-dimensional perturbations of the form $\delta f(\varpi)
\exp[i(m \theta + k_z z - \omega t)]$,
where the $\delta$ denotes an
Eulerian perturbation. Starting from Equations
(\ref{Eulerstart})-(\ref{Eulerend}), the linearized first
order equations are given by
\begin{eqnarray}
\label{rotceq}
i \bar{\omega} \delta \rho &=& \frac{1}{\varpi} (\varpi \rho \delta u)' + \frac{i m \rho}{\varpi}
\delta v + i k_z \rho \delta w \\
i \bar{\omega} \delta u &=& \frac{1}{\rho}\delta P' + \frac{\bar{g}}{\rho} \delta \rho
- 2 \Omega \delta v \\
i \bar{\omega} \delta v &=& \frac{i m}{\varpi \rho} \delta
  P + 2B \delta u\\
i \bar{\omega} \delta w &=& \frac{i k_z}{\rho} \delta P \\
\label{roteeq}
i \bar{\omega} \delta \rho &=& -C_L \rho \delta u + \frac{i \bar{\omega}}{s^2} \delta P,
\end{eqnarray}
where $\bar{\omega} = \omega - m \Omega$ is the phase speed in the locally
corotating frame, $B = (\Omega + d(\varpi \Omega)/d\varpi)/2$ is Oort's B
constant, $C_L = \gamma^{-1} d \ln P/ d \varpi - d \ln
\rho/d \varpi$ is the Ledoux discriminant, and the primes
denote differentiation with respect to
$\varpi$. Note that in Equations (\ref{rotceq})-(\ref{roteeq}), we have
made Cowling's approximation by ignoring the first order perturbation
to the gravitational potential \citep{Cowling}. This means that we are
taking the self-gravity of the BL to be negligible.

We can reduce the above set of linear equations to a pair of equations in
$\delta P$ and $\delta u$:
\begin{eqnarray}
\label{deltaPeq}
i \left(\frac{\bar{\omega}^2}{s^2} - k_z^2 - \frac{m^2}{\varpi^2}\right) \delta P &=
&\rho \left(C_L \bar{\omega} + \frac{2Bm}{\varpi} \right) \delta u
+ \frac{\bar{\omega}}{\varpi} (\varpi \rho \delta u)' \\
\label{deltaueq}
i\rho (\bar{\omega}^2 - \bar{g} C_L - \kappa^2) \delta u &=& \bar{\omega}
\delta P' + \frac{ \bar{\omega} \bar{g}}{s^2} \delta P - \frac{2 \Omega
  m}{\varpi} \delta P.
\end{eqnarray}
Here we have introduced the epicyclic frequency $\kappa^2 = 4 B \Omega$.

One generally expects the BL width $\delta_{BL}$ to be small compared to
the stellar radius, $\delta_{BL}\ll R_*$. Under this assumption, we
show in Appendix \ref{Reduction}
that Equations (\ref{deltaPeq})-(\ref{deltaueq}) can be simplified
to the following set of equations:
\begin{eqnarray}
\label{deltaPeq1}
i \left(\frac{\bar{\omega}^2}{s^2} - k_z^2 - \frac{m^2}{R_*^2}\right) \delta P &=
&\rho \left(C_L \bar{\omega} + \frac{m}{R_*} S \right) \delta u
+ \bar{\omega} (\rho \delta u)' \\
\label{deltaueq1}
i\rho (\bar{\omega}^2 - \bar{g} C_L) \delta u &=& \bar{\omega}
\delta P' + \frac{ \bar{\omega} \bar{g}}{s^2} \delta P.
\end{eqnarray}
The form of Equations (\ref{deltaPeq1}) and (\ref{deltaueq1}) is the
same as for a plane-parallel stratified shear flow \citep{Alexakis}.
This suggests that
the effects of Coriolis force and curvature are unimportant in a
radially thin BL if shear instabilities provide the turbulence, and
leads us to redefine the cylindrical $(\varpi, \theta, z)$
coordinate system into a Cartesian $(x, y, z)$ system. The perturbed
quantities now have the form $\delta f(x)
\exp[i(k_y y + k_z z - \omega t)]$, where $k_y \equiv m/R_*$. From this it
immediately follows that $\bar{\omega}  = \omega - k_y V_y(x)$, where $V_y(x) =
R_*\Omega(x)$ is the velocity profile. From here on, we will
  ignore the rotation terms, and treat the flow in the BL as a
  plane parallel shear flow.

From Squire's theorem for plane parallel shear flows \citep{Fejer},
any three-dimensional perturbation is mathematically equivalent to a
  two-dimensional perturbation upon making the transformations $k_y
  \rightarrow k_y/\cos \theta$, $k_z \rightarrow 0$, and $V_y \rightarrow
  V_y \cos \theta$, where $\cos \theta =
k_y/\sqrt{k_y^2+k_z^2}$. Thus, it is sufficient to only consider
two-dimensional perturbations ($k_z = 0$) in the stability problem,
and for the rest of this paper we take perturbations in the form
\begin{eqnarray}
\label{pertform}
\delta f(x) \exp[i(k_y y - \omega t)].
\end{eqnarray}

Assuming two-dimensional
perturbations, Equations (\ref{deltaPeq1}) and (\ref{deltaueq1}) become
\begin{eqnarray}
\label{deltaueq2}
i \left(\frac{\bar{\omega}^2}{s^2} - k_y^2\right) \delta P &=
&\rho \left(C_L \bar{\omega} + k_y V_y' \right) \delta u
+ \bar{\omega} (\rho \delta u)' \\
\label{deltaPeq2}
i\rho (\bar{\omega}^2 - \bar{g} C_L) \delta u &=& \bar{\omega}
\delta P' + \frac{ \bar{\omega} \bar{g}}{s^2} \delta P,
\end{eqnarray}
where the primes now denote differentiation with respect to
$x$. Equations (\ref{deltaPeq2}) and (\ref{deltaueq2}) can be used to
obtain the generalized Rayleigh equation (Appendix
\ref{gRapp}), which using a notation similar to \citet{Alexakis} reads
\begin{eqnarray}
\label{gen_Rayleigh}
\delta \phi'' + \left(k_x^2 - g\tilde{\rho}\frac{k_s + k_g}{\tilde{W}^2} -
  \frac{\tilde{W}''}{\tilde{W}} \right) \delta \phi &=& 0.
\end{eqnarray}
Here,
\begin{eqnarray}
\label{deltaphi_eq}
\delta \phi = -\delta u
\sqrt{\rho}/k_x,
\end{eqnarray}
is a modified stream function, and for simplicity we have dropped the bar over
$\bar{g}$, so now $g$ denotes the effective gravity. We have
also defined the quantities:
\begin{flalign}
\tilde{W} &= k_y W \sqrt{\tilde{\rho}}/ik_x \\
W &= V_y - c \text{ where }  c = \omega/k \text{ is the phase speed.}
\\
\label{kxwavevec}
k_x^2 &= k_y^2(W^2/s^2-1) \text{ is the square of the x-component of
  the wavevector} \\
& \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \text{in the absence of shear or
  stratification.} \\
\tilde{\rho} &= \rho f^2 \text{, where $\rho$ is the density.} \\
f &= \exp\left(\int_0^x k_g(\zeta) d\zeta \right). \\
k_g &= g/s^2 \text{ is a measure of the inverse of the local scale height.} \\
k_s &= \rho'/\rho \text{ is the inverse stratification length scale.}
\end{flalign}
Note that the Ledoux discriminant is given by $C_L = -(k_g + k_s)$.

The boundary conditions on $\delta \phi$ that must be satisfied at a
discontinuity in the density or the velocity are that the upper and
lower fluids stay in contact and that the pressure perturbation is
continuous across the interface. We show explicitly in Appendix
\ref{BCapp} that these conditions can be formulated as
\begin{eqnarray}
\label{BC_contact}
\frac{\delta \phi_+}{\tilde{W}_+} &=& \frac{\delta \phi_-}{\tilde{W}_-} \\
\label{BC_derivative}
\tilde{W}_+\delta \phi_+' - \tilde{W}_+'\delta \phi_+ - \frac{g_+ \tilde{\rho}_+
  \delta \phi_+}{\tilde{W}_+} &=& \tilde{W}_-\delta \phi_-' - \tilde{W}_-'\delta \phi_- - \frac{g_- \tilde{\rho}_-
  \delta \phi_-}{\tilde{W}_-}.
\end{eqnarray}
The $+$/$-$ signs denote evaluation of a
  quantity in the upper/lower fluid at the location of the
  interface.

\section{The Vortex Sheet without Gravity}
\label{compress_KH}

Equations (\ref{gen_Rayleigh})-(\ref{BC_derivative}) are fully general
and apply to an arbitrary velocity profile. In particular, we can
assume that the velocity varies discontinuously at some radius
$x=0$:
\begin{eqnarray}
\label{vsheet}
V_y(x) = \left\{
     \begin{array}{lr}
       \bar{V}_y, \ x > 0 \\
      -\bar{V}_y, \ x < 0,
     \end{array}
   \right.
\end{eqnarray}
where $\bar{V}_y$ is a constant. This is known as a
vortex sheet approximation. It represents the simplest possible
description of the velocity variation between the two limiting
values by essentially ignoring the details of the transition.
Subsequently in \S \ref{finite_width_SL}, we explore a more
realistic model of the velocity variation, in which the
transition occurs in a region of finite radial width. We
point out that the vortex sheet approximation is valid
for $k_y \delta_{BL} \ll 1$. We show in Appendix \ref{KH_conv} that in this limit
the dispersion relation for a finite width layer of constant shear
reduces to the vortex sheet dispersion relation.

We assume in this section that $g=0$, so there is no gravity,
and that $\rho$ and $s$ are constant above and below the interface, but can be
discontinuous across it. This case has already been considered by other
authors in the past including \citet{Landau}, \citet{Hatanaka},
\citet{Pai}, \citet{MilesKH}, and \citet{Gerwin}. However, the results for the case
without gravity will often be referenced later in the paper and serve as a
verification of the formalism we developed in \S \ref{goveqn}.

Using the velocity profile (\ref{vsheet}), the generalized Rayleigh
equation (\ref{gen_Rayleigh}) becomes
\begin{eqnarray}
\label{RCS}
\delta \phi_\pm'' + k_{x,\pm}^2 \delta \phi_\pm &=& 0, \\
\end{eqnarray}
where just as in \S \ref{goveqn}, the $+$/$-$ signs denote the upper/lower fluids respectively.
Since $k_{x,\pm}$ is constant for $V_y(x)$ given by Equation (\ref{vsheet}),
we have that
\begin{eqnarray}
\label{pert_eq}
\delta \phi_\pm &\propto& e^{-i k_{x,\pm} x}.
\end{eqnarray}
In general, $k_{x,\pm}$ is complex, and the sign is determined by applying
the appropriate boundary conditions. We discuss the boundary
conditions shortly, which will also make clear the reason for the negative
sign in the exponential of Equation (\ref{pert_eq}). Plugging
the expressions for $\delta \phi_\pm$ into Equation (\ref{BC_derivative})
gives
\begin{eqnarray}
\tilde{W}_+ k_{x,+} \delta \phi_+ = \tilde{W}_- k_{x,-} \delta \phi_- ,
\end{eqnarray}
where we have used $\delta \phi_\pm' = i k_{x,\pm} \delta \phi_\pm$,
$\tilde{\rho}_\pm = \rho_\pm$, and $\tilde{W}_\pm' = 0$. Using Equation (\ref{BC_contact})
to substitute for $\delta \phi_-$ in terms of $\delta \phi_+$, and
substituting for $\tilde{W}$ in terms of $W$ and $\rho$ we have
\begin{eqnarray}
\label{preprecompressdis}
\frac{k_y}{k_{x,+}}\rho_+W_+^2 = \frac{k_y}{k_{x,-}}\rho_-W_-^2.
\end{eqnarray}
Introducing the density ratio $\epsilon = \rho_+/\rho_-$, the Mach
number in the upper fluid $M =
\bar{V}_y/s_+$, and the phase speed normalized by the sound speed in
the upper fluid $\varphi = c/s_+$, we have
\begin{eqnarray}
\label{precompressdis}
\epsilon (M-\varphi)^2 \sqrt{(M+\varphi)^2 \frac{s_+^2}{s_-^2} - 1} =
(M+\varphi)^2 \sqrt{(M-\varphi)^2 - 1}.
\end{eqnarray}
By the definition of the sound speed, $s^2 = \gamma P/ \rho$, the
condition of pressure balance everywhere throughout the flow
requires that $\gamma_+^{-1}\rho_+s_+^2 = \gamma_-^{-1}\rho_-s_-^2$,
where
$\gamma_+$ and $\gamma_-$ are the adiabatic indices above and below
the interface.
Assuming $\gamma_+ = \gamma_-$, we have from pressure balance that
$(s_-/s_+)^2 = \rho_+/\rho_- = \epsilon$. Thus, the
dispersion relation becomes:
\begin{eqnarray}
\label{compressdis}
\epsilon (M-\varphi)^2 \sqrt{(M+\varphi)^2 \epsilon^{-1} - 1} =
(M+\varphi)^2 \sqrt{(M-\varphi)^2-1}.
\end{eqnarray}
\citet{MilesKH} has studied Equation (\ref{compressdis}) and found that
the stability criterion, i.e. that $\varphi$ is purely real, is
given by
\begin{eqnarray}
M > M_\text{crit} = \frac{1}{2}(1 + \epsilon^{1/3})^{3/2}.
\label{Mcritunstable}
\end{eqnarray}
This shows the surprising result that infinitesimal disturbances are
stabilized at high Mach number.

However, \citet{Fejer} pointed out
that due to Squire's theorem, it is always possible to choose an angle
$\theta$ for the wavevector with respect to the flow velocity such
that the projected Mach number $M \cos\theta$ is smaller than
$M_\text{crit}$. Thus, even at high Mach number, the vortex sheet
without gravity is still unstable to three dimensional
disturbances, which are almost perpendicular to the direction of the
flow; these unstable oblique modes resemble classical KH
modes. However, in the astrophysical context, a thin disk has a
scale height
$s/\Omega \ll R_*$, and the wavelength of the oblique modes will be small relative to
the disk scale height only for very small wavelengths $\lambda \ll
s/\Omega$. If the BL itself has a thickness $\delta_{BL} \gtrsim
s/\Omega$, the vortex sheet approximation for the oblique modes is
invalid, since either the modes don't fit into a disk
scale height, or the condition $\lambda \gg \delta_{BL}$ is not
satisfied.

\subsection{Solutions for $M \gg 1$}
\label{ratiomach}

Since we are interested in the high Mach number limit for the
initiation of the BL, we now find analytical solutions to
the dispersion relation (\ref{compressdis}) for $M \gg M_\text{crit}$,
where $M_\text{crit}$ was defined in Equation (\ref{Mcritunstable}).
Squaring Equation (\ref{compressdis}), one obtains a sixth order polynomial
in $\varphi$:
\begin{eqnarray}
\label{polynomial}
(\epsilon^2 - (M+\varphi)^2 \epsilon)(M-\varphi)^4 -
(1-(M-\varphi)^2)(M+\varphi)^4 = 0.
\end{eqnarray}
This polynomial has two easy to find analytic factors \citep{MilesKH}
\begin{eqnarray}
\label{tworoots}
\varphi = \left \{-M\left(\frac{1-\epsilon^{1/2}}{1+\epsilon^{1/2}}\right),
-M\left(\frac{1+\epsilon^{1/2}}{1-\epsilon^{1/2}}\right) \right \}.
\end{eqnarray}
We now take the limit $M \gg M_\text{crit}$, in which case to terms of order
$\mathcal{O}(M^{-4})$, the other four solutions to the polynomial in Equation
  (\ref{polynomial}) are
\begin{multline}
\label{fourroots}
\varphi = \left \{M+1+\frac{\epsilon}{2(2M+1)^2}, \ M-1-\frac{\epsilon}{2(2M-1)^2}, \right. \\
 -M + \epsilon^{1/2} + \frac{\epsilon^{1/2}}{2(2M-\epsilon^{1/2})^2}, \
\left. -M -
\epsilon^{1/2} - \frac{\epsilon^{1/2}}{2(2M+\epsilon^{1/2})^2} \right \}.
\end{multline}

It is easy to see that the six roots of the polynomial
(\ref{polynomial}) correspond to sound waves.
Starting from Equation (\ref{kxwavevec}) and rearranging terms we obtain
\begin{eqnarray}
W^2 k_y^2 = s^2(k_x^2 + k_y^2).
\end{eqnarray}
As long as $k_x$ is real, this is the dispersion relation for a sound
wave, since $W^2$ is the square of the phase velocity in the frame
comoving with the fluid. Plugging in the six roots from Equations
(\ref{tworoots}) and
(\ref{fourroots}) into Equation (\ref{kxwavevec}), it is
straightforward to verify that $k_{x,\pm}$ are indeed real for all of
them.
Each of the four roots in Equation (\ref{fourroots}) has a further
simple
interpretation. The first two correspond to sound waves that propagate
almost parallel to the interface in the $+y$ and $-y$ directions in the upper
fluid, whereas the second two
correspond to sound waves that propagate almost parallel to the interface
in the $+y$ and $-y$
directions in the lower fluid. The two roots in Equation
(\ref{tworoots}) are more difficult to
interpret, but the first of these corresponds to a standing wave when
the two fluids have equal density (i.e. $\epsilon = 1$).

Since each of the six roots for $M \gg M_\text{crit}$
yields a real value for $k_{x,\pm}$, the solutions do not damp
away from the interface, and we need to
apply radiation boundary conditions at $x=\pm \infty$. The
proper procedure is to demand
that all waves are outgoing in each fluid in a frame which is subsonic
with respect to the fluid \citep{Milessound}, and it is convenient to
work in the comoving frame of each fluid. The dimensionless phase
velocity in the frame comoving with the upper fluid is
$\varphi_{+,CF} = \varphi - M$, and the analogous expression for the
lower fluid is $\varphi_{-,CF} = \varphi + M$.

In order to satisfy Equation (\ref{preprecompressdis}), $k_{x,+}$ and
$k_{x,-}$ must have the same sign, so we
must have either $\delta \phi_\pm \propto e^{i(|k_{x,\pm}| x - \omega
  t)}$ or $\delta \phi_\pm \propto e^{i(-|k_{x,\pm}| x - \omega
  t)}$. Moreover, since $k_{x,+}$ and $k_{x,-}$ have the same sign, it
is clear that $\varphi_{+,CF}$ and  $\varphi_{-,CF}$ must have the
opposite sign to yield outgoing waves in the comoving frames of each of
the two fluids.
Only three of the six roots found above satisfy this condition:
\begin{eqnarray}
\label{lbranch}
\varphi_l &=& -M + \epsilon^{1/2} + \frac{\epsilon^{1/2}}{2(2M-\epsilon^{1/2})^2} \\
\label{mbranch}
\varphi_m &=& -M\left(\frac{1-\epsilon^{1/2}}{1+\epsilon^{1/2}}\right) \\
\label{ubranch}
\varphi_u &=& M-1-\frac{\epsilon}{2(2M-1)^2}.
\end{eqnarray}
We will refer to these three roots as the lower, middle,
and upper branches, respectively. Furthermore, it is straightforward
to check that the
solutions which yield outgoing waves in both the upper and lower
fluids have $\delta \phi_\pm \propto
e^{i(-|k_{x,\pm}| x - \omega t)}$, so $k_{x,\pm}$ are positive given
our definition (\ref{pert_eq}).

\subsection{Dispersion Relation in the General Case}

We now relax the assumption of $M \gg 1$ and numerically solve the dispersion
relation (\ref{compressdis}) at arbitrary Mach number for
$\epsilon = 1$ (Figure \ref{drfig}a) and $\epsilon = .01$ (Figure \ref{drfig}b). At the
critical Mach number given by Equation (\ref{Mcritunstable}), the
upper and lower branches merge together in the real plane and bifurcate in the
complex plate. These bifurcated solutions turn into the two
incompressible KH modes for $M \ll M_\text{crit}$.
Unlike, the lower and upper branches, the middle branch has no
incompressible analog and
ceases to be a viable physical solution below a critical Mach number
\begin{eqnarray}
2M_m = 1 + \epsilon^\frac{1}{2}.
\label{Mcritmiddle}
\end{eqnarray}
The reason for this is that below $M=M_m$, $k_{x,\pm}$ switches from real
to imaginary and the boundary conditions for the middle branch at $x =
\pm \infty$ can no longer be satisfied.

\begin{figure}[!h]
  \centering
  \subfigure[]{\includegraphics[width=.49\textwidth]{dreps1g0}}
  \subfigure[]{\includegraphics[width=.49\textwidth]{dreps01g0}}
  \caption{The dispersion relation $\varphi/M$ as a
  function $M$ for $\epsilon = 1$ and $\epsilon = .01$. The solid and
  dashed curves give the real and imaginary components respectively.
  The blue, black, and red curves
  correspond to the lower, middle, and upper branches
  respectively. The red and blue curves have been slightly offset
  vertically so they do not overlap.}
 \label{drfig}
\end{figure}

\section{The Isothermal Vortex Sheet with Gravity}
\label{compress_grav}

We now go beyond the simplifying assumption of no gravity used
in  \S \ref{compress_KH} and consider the case where $g$ is non-zero
and constant. We shall shortly assume an isothermal equation of state,
but for now we consider the more general
polytropic equation of state of the form $P = K_\pm \rho^n$ in each of
the two fluids. For any equation
of state of this form, we have
\begin{eqnarray}
k_g = \frac{g}{s^2} = -\frac{\frac{\partial P}{\partial \rho}\frac{d
    \rho}{dr}}{\gamma P} = - \frac{n}{\gamma}\frac{d \ln \rho}{dr} =
-\frac{n}{\gamma} k_s,
\end{eqnarray}
which means $k_s + k_g = (1-\gamma/n)k_g$. If $n = \gamma$,
corresponding to the case of an adiabatic atmosphere, the
second term in parentheses in Equation (\ref{gen_Rayleigh}) drops out,
but then the third term becomes difficult to treat
analytically. However, for an isothermal atmosphere ($n=1$) the sound speed
is constant, and both the second and third terms in Equation
(\ref{gen_Rayleigh}) reduce to a tractable form. Since we are only interested
in the gross, qualitative properties of the flow, we
assume both fluids are isothermally stratified, since this
assumption significantly simplifies the analytical treatment.

For $n = 1$, the
second term in Equation (\ref{gen_Rayleigh}) becomes
\begin{eqnarray}
-g\tilde{\rho}\frac{k_s + k_g}{\tilde{W}^2} &=& g k_g (1-\gamma)
\frac{k_x^2}{k_y^2 W^2}
\label{isosecondterm}
\end{eqnarray}
We now again assume a vortex sheet velocity profile (Equation
(\ref{vsheet})), in which case for $n=1$, the third term in Equation
(\ref{gen_Rayleigh})
is given by
\begin{eqnarray}
\label{isothirdterm}
\frac{\tilde{W}''}{\tilde{W}} =
\frac{(\sqrt{\tilde{\rho}})''}{\sqrt{\tilde \rho}} = \left(\frac{2-\gamma}{2}k_g\right)^2
\end{eqnarray}
Putting Equations (\ref{isosecondterm}) and (\ref{isothirdterm}) into
Equation (\ref{gen_Rayleigh}), we have
\begin{eqnarray}
&&\delta \phi'' + \tilde{k}_x^2 \delta \phi = 0, \nonumber\\
&& \tilde{k}_x^2\equiv\left[k_x^2
  \left(1-(\gamma-1) \left(\frac{k_g}{k_y}\right)^2
  \left(\frac{s}{W}\right)^2\right) - \left(\frac{2-\gamma}{2}
  k_g\right)^2 \right].
\label{Rayleigh_isothermal}
\end{eqnarray}
The perturbation is then given as
\begin{eqnarray}
\label{pertwithgrav}
\delta \phi_\pm \propto e^{-i
\tilde{k}_{x,\pm} x}.
\end{eqnarray}

We now comment on the terms present in Equation
(\ref{Rayleigh_isothermal}). In the absence of gravity, $k_g = 0$,
and we simply have $\tilde{k}_x = k_x$. If $k_g \ne 0$, then the second
term on the right hand side of Equation (\ref{Rayleigh_isothermal})
can be written in a more familiar form as
\begin{eqnarray}
(\gamma-1) \left(\frac{k_g}{k_y}\right)^2 \left(\frac{s}{W}\right)^2 =
\left(\frac{N}{\omega - k_yV_y}\right)^2,
\end{eqnarray}
where
\begin{eqnarray}
\label{eqBrunt}
N^2 = \frac{(\gamma - 1) g^2}{s^2}
\end{eqnarray}
is the Brunt-V\"{a}is\"{a}l\"{a} frequency. Thus, the second term on
the right hand side of Equation (\ref{Rayleigh_isothermal}) provides a
lower frequency cutoff for sound waves at
the Brunt-V\"{a}is\"{a}l\"{a} frequency. The last term in Equation
(\ref{Rayleigh_isothermal}) arises because we have {\it not} made the short
wavelength approximation $k_x \gg k_g$.

Now, we use the boundary conditions (\ref{BC_contact}) and
(\ref{BC_derivative}) to determine the dispersion relation.
Substituting $\delta \phi_\pm' = -i \tilde{k}_{x,\pm} \delta \phi_\pm$, and
$\tilde{W}_\pm' = (2-\gamma) k_{g,\pm} \tilde{W}_\pm/2$ into the condition
(\ref{BC_derivative}), we have
\begin{eqnarray}
-\left[i\tilde{k}_{x,+} + \left(\frac{2-\gamma}{2}
\frac{g}{s_+^2}\right)\right] \tilde{W}_+ \delta \phi_+ - \frac{g \tilde{\rho}_+ \delta \phi_+}{\tilde{W}_+}
= -\left[i\tilde{k}_{x,-} + \left(\frac{2-\gamma}{2}
\frac{g}{s_-^2}\right)\right]\tilde{W}_-\delta \phi_- - \frac{g \tilde{\rho}_- \delta \phi_-}{\tilde{W}_-}.
\end{eqnarray}
Next, using the second boundary condition (\ref{BC_contact}) to substitute for
$\delta \phi_-$ in terms of $\delta \phi_+$, and using the fact that $\tilde{\rho}_\pm
= \rho_\pm$ at the interface, we have
\begin{eqnarray}
-\left[i\tilde{k}_{x,+} + \left(\frac{2-\gamma}{2}
\frac{g}{s_+^2}\right)\right] \tilde{W}_+^2 - g \rho_+
= -\left[i\tilde{k}_{x,-} + \left(\frac{2-\gamma}{2}
\frac{g}{s_-^2}\right)\right]\tilde{W}_-^2 - g \rho_-.
\end{eqnarray}
We now introduce the dimensionless gravity parameter
\begin{eqnarray}
G = \frac{g}{k_y s_+^2},
\label{eq:grav}
\end{eqnarray}
which is closely related to the ratio
of the wavelength to the pressure scale height, $h_s$. Thus, $G \sim 1$ when
$k_y h_{s,+}
\sim 1$, and $G \epsilon^{-1} \sim 1$ when $k_y h_{s,-} \sim 1$.
Using $\epsilon =
\rho_+/\rho_- = s_-^2/s_+^2$, and performing some algebra, we have
\begin{eqnarray}
\left[i\frac{\tilde{k}_{x,+}}{k_y} + \left(\frac{2-\gamma}{2}G\right)\right]
\left(\frac{k_y}{k_{x,+}}\right)^2\left(\frac{{W}_+}{s_+}\right)^2\epsilon + G
(1-\epsilon)
= \left[i\frac{\tilde{k}_{x,-}}{k_y} + \left(\frac{2-\gamma}{2}G
\epsilon^{-1}\right)\right]\left(\frac{k_y}{k_{x,-}}\right)^2\left(\frac{{W}_-}{s_+}\right)^2.
\label{mastergravity}
\end{eqnarray}
Setting $G=0$ it is clear that we recover the vortex sheet dispersion
relation in the absence of gravity (Equation (\ref{preprecompressdis})).

We now check Equation (\ref{mastergravity}) by showing
that it reproduces the well-known incompressible KH dispersion
relations in the limit
$M \ll M_\text{crit}$, before going on to treat the case of the supersonic vortex sheet with
gravity.

\subsection{Highly Subsonic Vortex Sheet with Gravity}
\label{subsonicKH}

We assume that $\bar{V}/s_\pm \ll 1$, which immediately implies $M \ll
M_\text{crit}$, and we also assume $G \ll 1$ and $G\epsilon^{-1} \ll 1$, which
means that the wavelength of the perturbation is much smaller than the
scale height in both the upper and lower fluids ($k_{g,\pm}/k_y \ll
1$). Next, we eliminate
sound wave modes by assuming
$\varphi \ll 1$ (phase velocity much lower than sound
velocity in upper fluid) and  $\varphi \epsilon^{-1/2} \ll 1$ (phase
velocity much lower than sound velocity in lower fluid). According to
Equation (\ref{Rayleigh_isothermal}), this means $k_x = \pm i k_y$,
and given our definition for $\delta \phi_\pm$ in Equation
(\ref{pertwithgrav}), we must chose the minus sign in the upper fluid
and the plus sign in the lower fluid to give vanishing solutions at $x =
 \pm \infty$.
Equation (\ref{mastergravity}) then yields
\begin{eqnarray}
\label{subsonicdispers}
-\epsilon \left(\frac{\omega}{k_y}-\bar{V_y}\right)^2\sqrt{1-\frac{(\gamma-1)G^2}{(M-\varphi)^2}} +
\frac{g}{k_y}(1-\epsilon) = \left(\frac{\omega}{k_y} + \bar{V_y}\right)^2\sqrt{1-\frac{(\gamma-1)G^2\epsilon^{-1}}{(M+\varphi)^2}}.
\end{eqnarray}
Next, we note that $(\gamma-1)G^2/(M-\varphi)^2 =
N_+^2/(\omega-k_y\bar{V_y})^2$
and that $(\gamma-1)G^2\epsilon^{-1}/(M+\varphi)^2 =
N_-^2/(\omega+k_y\bar{V_y})^2$, where $N$ is the
Brunt-V\"{a}is\"{a}l\"{a} frequency for an isothermal atmosphere and
was given in Equation (\ref{eqBrunt}). We
expect $|\omega \pm k_y \bar{V}_y| \gtrsim \sqrt{gk}$, which is the
characteristic frequency of surface gravity waves. Since $N \pm \sim
\sqrt{g k_{g,\pm}}$, and we have already assumed $k_{g,\pm}/k_y \ll
1$, it follows that $N_\pm \ll
|\omega \mp k_y\bar{V_y}|$. Consequently, Equation
(\ref{subsonicdispers}) reduces to
\begin{eqnarray}
\label{KHwellknown}
-\epsilon \left(\frac{\omega}{k_y} - \bar{V_y}\right)^2 +
\frac{g}{k_y}(1-\epsilon) = \left(\frac{\omega}{k_y} + \bar{V_y}\right)^2,
\end{eqnarray}
which is the well known KH dispersion relation for an incompressible
fluid in the presence of gravity.

\subsection{The Weak Gravity Limit at High Mach Number}
\label{weakgravity}

As mentioned before in \S \ref{ratiomach}, \citet{MilesKH} has
demonstrated that in the absence of gravity,
the vortex sheet is stable above a critical Mach number given by
Equation (\ref{Mcritunstable}). We now address the question of whether
the system still remains stable when $G \ne 0$.

To answer this question, we will use our general dispersion relation
(\ref{mastergravity}) in which we will additionally assume $M \gg
M_\text{crit}$, since
this assumption significantly simplifies the algebra. Since we have
already obtained solutions for the case $G=0$ and $M \gg M_\text{crit}$
(\S \ref{ratiomach}), we proceed by considering gravity as a
perturbation. We consider wavelengths that are much smaller than the
pressure scale height and ask what happens in the limit
$G \rightarrow 0$. In our analysis, we will assume that the density
ratio $\epsilon \le 1$, so that the system is
stable to the Rayleigh-Taylor instability.

Because the introduction of gravity makes the upper and lower fluids
stratified, some care should be taken when determining which solutions
are physical and which are not. \citet{Milessound} has shown that for
$G=0$, the three supersonic KH modes can be understood in terms of sound
waves emitted from the interface between the two fluids. This
interpretation is useful as well for the case with gravity and
leads to a couple of insights. First, the amplitude of sound waves
is not constant as they propagate through a stratified medium. Rather,
to conserve energy, waves propagating upward (to lower densities) must
increase in amplitude, and those propagating downward (to higher densities) must
decrease in amplitude. Second, since the sound waves are emitted from the interface,
if $\omega$ has an imaginary component, then the amplitude of the
emitted waves changes with time. Thus, if $\Im[\omega] > 0$, and
there is an instability, then the amplitude of the waves will decay
with distance from the interface, since the waves emitted in the past had lower
amplitude. Conversely if $\Im[\omega] < 0$ and the perturbation is
damping in time, then the amplitude of the emitted waves will increase with
distance from the interface, since the waves emitted in the past had a
larger amplitude. Both of these effects mean that the amplitude of the
waves can blow up as we move away from the interface. Thus, we take as
physical those solutions
which yield outgoing waves (away from the interface) in both the upper
and lower fluids, even if
these solutions diverge as $x \rightarrow \pm \infty$. For the limit
$G \rightarrow 0$, this means that the physical solutions are
still the lower, middle, and upper branches but now modified by the presence
of a weak gravitational field.

We note here that although the dispersion relation
(\ref{mastergravity}) is valid for all values of $G$ and not just for
small $G$, the simple picture of outgoing sound waves in the upper and
lower fluids is only valid for $G \rightarrow 0$. From a physical
point of view, this can be attributed to the following fact. Sound waves
(p-modes) traveling in a stratified atmosphere have a frequency
cutoff at $\omega^2 = N^2$ below which propagation is not
possible. Given our definitions of $k_y$ and
$\tilde{k}_x$, the frequency of a sound wave in the frame comoving with
the fluid is $\omega^2 \sim (k_y^2 + \tilde{k}_x^2) s^2$. Substituting
$N^2 \sim \omega^2$ and using the definition of $N^2$
(Equation (\ref{eqBrunt})) yields
$G^2<(\gamma-1)^{-1}(1 + (|\tilde{k}_x|/k_y)^2)$ for sound waves to propagate.
If this condition is not fulfilled, then the
picture of outgoing sound waves is invalid, and the boundary
conditions need to be formulated in a different way, which is beyond
the scope of the present work.

\subsubsection{The Lower Branch}
\label{lbsec}

We begin by considering how the lower wave is modified in the limit $G \rightarrow
0$ and $M \gg M_\text{crit}$. If $G = 0$ exactly, then $\tilde{k}_x = k_x$, and in the limit $G
\rightarrow 0$ we have $\tilde{k}_x = k_x
\left(1 + \mathcal{O}(G^2)\right)$. Keeping terms only to first order
in $G$, Equation (\ref{mastergravity}) becomes
\begin{eqnarray}
\left[i\frac{k_{x,+}}{k_y} + \left(\frac{2-\gamma}{2}G\right)\right]
\left(\frac{k_y}{k_{x,+}}\right)^2\left(\frac{{W}_+}{s_+}\right)^2\epsilon + G
(1-\epsilon)
= \left[i\frac{k_{x,-}}{k_y} + \left(\frac{2-\gamma}{2}G
\epsilon^{-1}\right)\right]\left(\frac{k_y}{k_{x,-}}\right)^2\left(\frac{{W}_-}{s_+}\right)^2.
\label{mastergravitysimple}
\end{eqnarray}
Writing this out explicitly in terms of $M$ and $\varphi$ gives
\begin{eqnarray}
\label{l0eq}
\frac{i(M-\varphi)^2 \epsilon}{\sqrt{(M-\varphi)^2-1}} - \frac{i (M+\varphi)^2}{\sqrt{(M+\varphi)^2\epsilon^{-1}-1}} +
G(1-\epsilon) =
\frac{2-\gamma}{2}G\left(\frac{(M+\varphi)^2
  \epsilon^{-1}}{(M+\varphi)^2\epsilon^{-1}-1} - \frac{(M-\varphi)^2\epsilon}{(M-\varphi)^2-1}
\right).
\end{eqnarray}

Next, we assume that gravity only weakly affects the
dispersion relation and make the perturbative expansion
\begin{eqnarray}
\label{phiseries}
\varphi = \varphi_0 + \varphi_1, \ \
|\varphi_1|/|\varphi_0| \ll 1,
\end{eqnarray}
where $\varphi_0$ is the solution for $G=0$ and $M \gg M_\text{crit}$.

For the lower branch, we can use Equation (\ref{lbranch}) for $\varphi_0$, and in the
limit $M \gg M_\text{crit}$ we have
\begin{eqnarray}
\varphi_0 &\approx& -M + \epsilon^{1/2} \\
\label{lkpeq}
\frac{k_{x,+}}{k_y} &=& \sqrt{(M-\varphi_0)^2-1} \approx
2M-\epsilon^{1/2} \\
\frac{k_{x,-}}{k_y} &=& \sqrt{(M+\varphi_0)^2\epsilon^{-1} -1} \approx \frac{1}{2M-\epsilon^{1/2}}.
\end{eqnarray}

Defining
\begin{eqnarray}
\label{lmstar}
M_l \equiv 2M - \epsilon^{1/2},
\end{eqnarray}
Equation (\ref{l0eq}) becomes
\begin{eqnarray}
\label{l1eq}
\frac{i(M_l-\varphi_1)^2 \epsilon}{\sqrt{(M_l-\varphi_1)^2-1}} - \frac{i
  (\epsilon^{1/2}+\varphi_1)^2}{\sqrt{(\epsilon^{1/2}+\varphi_1)^2\epsilon^{-1}-1}} +
G(1-\epsilon) =
\frac{2-\gamma}{2}G\left(\frac{(\epsilon^{1/2}+\varphi_1)^2
  \epsilon^{-1}}{(\epsilon^{1/2}+\varphi_1)^2\epsilon^{-1}-1} -
  \frac{(M_l-\varphi_1)^2\epsilon}{(M_l-\varphi_1)^2-1}
\right).
\end{eqnarray}
It will turn out (Equation (\ref{lvp1eq})) that $\varphi_1$ is proportional to
  $G$, so working to first order in $G$ is equivalent to working to
  first order in $\varphi_1$. Equation (\ref{l1eq}) then simplifies to
\begin{eqnarray}
i (M_l-\varphi_1) \epsilon - \frac{i M_l \epsilon (1 +
  2\epsilon^{-1/2} \varphi_1)}{1+ M_l^2\epsilon^{-1/2} \varphi_1} +
G(1-\epsilon) = \frac{2-\gamma}{2}G\left(M_l^2\frac{1+2
  \epsilon^{-1/2} \varphi_1
  }{1 + 2 M_l^2 \epsilon^{-1/2}\varphi_1} - \epsilon \right).
\end{eqnarray}
Assuming $|M_l^2 \epsilon^{-1/2}\varphi_1| \ll 1$,
and keeping terms to leading order in $M_l$, we can solve for $\varphi_1$
  in terms of $G$.
\begin{eqnarray}
\label{lvp1eq}
\varphi_{1,l}
\approx - \frac{2-\gamma}{2}
\frac{G}{M_l\epsilon^{1/2}} i.
\end{eqnarray}
We immediately
see two things from Equation (\ref{lvp1eq}). First, $\varphi_1$ is
purely imaginary, and second, if $\gamma < 2$, $\varphi_1$ is
negative and the perturbation damps, whereas if $\gamma > 2$,
$\varphi_1$ is positive and the perturbation grows. For realistic
equations of state, $\gamma \le 5/3$, so the lower wave always damps.

\subsubsection{The Middle and Upper Branches}

We can find an approximate solution for $\varphi_1$ in the limit $G
\rightarrow 0$ and $M \gg M_\text{crit}$ for the middle and upper branches in much the same
manner as for the lower branch. The first order correction for the
middle branch is
\begin{eqnarray}
\label{vp1meq}
\varphi_{1,m} = -\frac{\gamma}{2}\frac{G(1-\epsilon^{1/2})}{\epsilon^{1/2}}i,
\end{eqnarray}
and for the upper branch is
\begin{eqnarray}
\label{uvp1}
\varphi_{1,u}
\approx \frac{2-\gamma}{2}\frac{G \epsilon^{1/2}}{2M-1}i.
\end{eqnarray}

Just as in the case of the lower branch, $\varphi_1$ is purely
imaginary for both the middle and upper branches. We see that
$\varphi_1$ is negative for the middle branch if $\epsilon < 1$. However, for the
upper branch if $\gamma < 2$, $\varphi_1$ is positive
and if $\gamma > 2$, $\varphi_1$ is negative. This is opposite from the
lower branch meaning that for any $\gamma \neq 2$ in the limit
$M \gg 1$, one of the
two branches is always unstable and the other one is damped. For a
realistic equation of state, $\gamma \le 5/3$ {\it so the upper branch is
the unstable one}, and the lower
branch damps.

\begin{comment}
One may wonder whether there are any completely new
modes introduced by the presence of gravity that have no analog in the
case $G=0$. As we show in Appendix \ref{completenessapp}, there are only three physically viable
solutions in the limit $G \rightarrow 0$, corresponding to the lower,
middle, and upper branches respectively.
\end{comment}

\subsubsection{Numerical Verification}

We verify Equations (\ref{lvp1eq}), (\ref{vp1meq}), and (\ref{uvp1})
numerically by solving the fully general
dispersion relation (\ref{mastergravity}) and comparing $\Im(\varphi)$
with our analytical estimate for the parameters $M=5$, $\epsilon = .5$,
and $\gamma = 5/3$. We plot $\Im(\varphi)$ vs. $G$ in
for both our analytical solutions (solid lines) and the ones obtained
numerically (dashed lines) in Figure \ref{modefig}. The analytical
solution converges to the numerical one in the limit $G \rightarrow 0$.

\begin{figure}[!h]
  \centering
  \subfigure[]{\includegraphics[width=.32\textwidth]{lowerwave_check}}
  \subfigure[]{\includegraphics[width=.32\textwidth]{middlewave_check}}
  \subfigure[]{\includegraphics[width=.32\textwidth]{upperwave_check}}
  \caption{Plots of $\Im(\varphi)$ vs. $G$ for the lower (a), middle
  (b), and upper (c) branches respectively, using $M=5$, $\epsilon=.5$,
  $\gamma = 5/3$. The solid curves show the
  numerical solution obtained by solving Equation
  (\ref{mastergravity}) and the dashed curves correspond to the
  approximate solutions from Equations (\ref{lvp1eq}), (\ref{vp1meq}),
  and (\ref{uvp1}).}
  \label{modefig}
\end{figure}

\section{Sonic Instabilities in a Finite Width Layer of Constant Shear}
\label{finite_width_SL}

The calculations presented in \S \ref{compress_KH},\ref{compress_grav}
were performed for the velocity profile (\ref{vsheet}) featuring
a discontinuity at some radius. We now consider a more complicated
(and more realistic at later stages of the BL evolution)
initial setup in which the velocity between the two fluids varies
{\it continuously} within a narrow shear layer. Unlike the vortex
sheet, the finite width shear layer without gravity is known
to be unstable at high Mach number
\citep{Glatzel,ChoudhuryLovelace,Ray}, and in this case the growth rate of the
instability scales inversely with the width of the shear layer and in
proportion to the shear, $S \sim \Omega_K R_*/\delta_{BL}$.

In the following, we extend some of the findings of \citet{Glatzel}
to study sonic instabilities in a finite width shear layer without
gravity and apply them to the problem of the BL initiation.
The setup we consider has the velocity profile
\begin{eqnarray}
\label{vsetup}
V(x) = \left\{
     \begin{array}{lr}
       \bar{V}, \ x > \delta_{BL}, \\
       \bar{V}x/\delta_{BL}, \ -\delta_{BL} \le x \le \delta_{BL}, \\
       -\bar{V}, \ x < -\delta_{BL}
     \end{array}
   \right.
\end{eqnarray}
and the density profile
\begin{eqnarray}
\label{rhosetup}
\rho(x) = \left\{
     \begin{array}{lr}
       \rho_+ , \ x > -\delta_{BL} \\
       \rho_-, \ x < -\delta_{BL}
     \end{array}
   \right.
\end{eqnarray}
Pressure equilibrium again requires that $\rho_+ s_+^2 = \rho_-
s_-^2$, which sets the sound speed everywhere in the flow, and as
before we have $\epsilon = \rho_+/\rho_-$ and $M = \bar{V}/s_+$.

Although we only consider a linearly varying velocity
profile in the shear layer, \citet{Ray} and
\citet{ChoudhuryLovelace} have found that different
shear profiles are qualitatively similar. Thus, we
consider a constant shear to be representative of more
general shear profiles.

\subsection{Dispersion Relation}

We now study the dispersion relation of the finite width shear
layer. In applying the dispersion relation to the initiation of the
BL, we are most interested in the
growth rate of the fastest growing mode for $M \gg 1$.
\citet{Glatzel} has already obtained the dispersion
relation for the case $\epsilon = 1$, and using his techniques, we
derive the dispersion relation for the case of arbitrary $\epsilon$ in
Appendix \ref{Glatzelapp}. We also show in Appendix \ref{KH_conv}
that the dispersion relation for a finite width shear layer reduces
to the dispersion relation for a vortex sheet in the limit $k_y
\rightarrow 0$ at constant $\delta_{BL}$.

In Figure \ref{growthfinite}, we plot the
dispersion relation (\ref{dis_hyp}) as a function of wavenumber
for the parameters $M=5$ and $\epsilon = 1$, $\epsilon = .25$, $\epsilon =
.1$, $\epsilon = .01$, and $\epsilon = 0$. The $\epsilon = 0$ case is equivalent
to having a hard reflecting boundary at $x = -\delta_{BL}$. The modes
depicted are the ones that converge to the upper and
lower branches from \S \ref{compress_KH} in the vortex sheet limit
($k_y\delta_{BL} \ll 1$). For $\epsilon
< 1$, the upper branch always has a larger growth rate than the lower
branch, and both the upper and lower branches have the same growth rate
for $\epsilon = 1$. In the $\epsilon = 1$ case, we can identify the upper and
lower branches as the $n_\pm = 0$ decoupled modes in \S 5.4 of
\citet{Glatzel}. In addition to the upper and lower waves,
\citet{Glatzel} has shown that there is an infinite spectrum of damped
modes, but we do not consider these here, since we are interested in
determining the growth rate of the fastest growing mode.

We now revisit the seemingly paradoxical
statement that the vortex sheet is stable to two dimensional
disturbances along the flow direction above a critical Mach number;
yet at the same time, the growth rate of the
fastest growing mode scales with $\delta_{BL}^{-1}$, which implies that the
thinner the shear
layer, the faster the instability proceeds. The key to resolving this
apparent controversy is to consider a mode having
$k_y\delta_{BL} \ll 1$. Its growth rate is diminished if we decrease
$\delta_{BL}$ while keeping $k_y$ constant for $M>M_\text{crit}$ and becomes
vanishingly small if we take the limit $\delta_{BL} \rightarrow 0$. This can be
seen from the dying left hand tail of the curve in
Figure \ref{growthfinite}b, and
keeping $k_y$ constant, while decreasing $\delta_{BL}$ we move
leftward along the
tail. This means $\Im[\omega]$ becomes smaller, since $\varphi = \omega/k_ys_+$
is directly proportional to $\omega$ for constant $k_y$ even as we
decrease $\delta_{BL}$.  We next point
out that the curve in Figure \ref{growthfinite}c remains unchanged in
shape or amplitude as we decrease
$\delta_{BL}$, keeping $k_y\delta_{BL}$ constant. Consequently, the value of $k_y$ for which the
maximum in $\Im[\omega]$ occurs $k_{y,\text{max}} \propto \delta_{BL}^{-1}$ and likewise
$\max[\Im[\omega]] \propto \delta_{BL}^{-1}$. Thus, as $\delta_{BL}$ is decreased, the
instability shifts to shorter wavelengths
and becomes more rapid, while at the same time, modes
having $k_y\delta_{BL} \ll 1$ are stabilized for a given value of $k_y$. Since any
real shear layer is likely to have
a nonzero width, one may therefore remark that taking
the vortex sheet limit for $M \gg M_\text{crit}$ masks the instability.

\begin{comment}
We can compare the dispersion relation for the $M \gg 1$ case to the $M \ll 1$
(incompressible) case for the setup given by Equations (\ref{vsetup}) and
(\ref{rhosetup}). For $M \ll 1$, the
dispersion relation is well known and is given by
e.g. \citep{ShtemlerMond}\footnote{Be wary of typos.} as
\begin{eqnarray}
&&\left[\frac{\omega}{k}(1+\epsilon) - \epsilon \frac{\bar{V}}{k\delta_{BL}}\right] \times
\left[\frac{\omega}{k} - \frac{\bar{V}}{k\delta_{BL}}F(2k\delta_{BL})\right] = -\epsilon
\left(\frac{\omega}{k}\right) \left(\frac{\bar{V}}{k
  L}\right)\left[G(2kL) + 1\right], \\
&&F(x) = \frac{1}{2}\exp(-2x) - \frac{1}{2} + x, \ G(x) = \exp(-2x) -
1 \\
&&F(x) = x^2 + \mathcal{O}(x^3), \ G(x) = -2x + \mathcal{O}(x^2) \text{
  for } x \ll 1.
\end{eqnarray}
Just as for $M \gg 1$, the fastest growing mode again has $k_yL
\sim 1$, $k_{y,\text{max}} \propto L^{-1}$, and $\max[\Im[\omega]] \propto
L^{-1}$. Now, taking the limit $L \rightarrow 0$, while
keeping $k_y$ constant, the instability growth rate converges to
a non-zero value, and we recover the conventional KH dispersion
relation in the vortex sheet limit (Equation (\ref{KHwellknown})).
\end{comment}

\begin{figure}[!p]
  \centering
  \subfigure{\includegraphics[width=.32\textwidth]{realphieps1}}
  \subfigure{\includegraphics[width=.32\textwidth]{imagphieps1}}
  \subfigure{\includegraphics[width=.32\textwidth]{imagomeps1}}
  \subfigure{\includegraphics[width=.32\textwidth]{realphieps025}}
  \subfigure{\includegraphics[width=.32\textwidth]{imagphieps025}}
  \subfigure{\includegraphics[width=.32\textwidth]{imagomeps025}}
  \subfigure{\includegraphics[width=.32\textwidth]{realphieps010}}
  \subfigure{\includegraphics[width=.32\textwidth]{imagphieps010}}
  \subfigure{\includegraphics[width=.32\textwidth]{imagomeps010}}
  \subfigure{\includegraphics[width=.32\textwidth]{realphieps001}}
  \subfigure{\includegraphics[width=.32\textwidth]{imagphieps001}}
  \subfigure{\includegraphics[width=.32\textwidth]{imagomeps001}}
  \subfigure{\includegraphics[width=.32\textwidth]{realphieps0}}
  \subfigure{\includegraphics[width=.32\textwidth]{imagphieps0}}
  \subfigure{\includegraphics[width=.32\textwidth]{imagomeps0}}
  \caption{Dispersion relations for $M=5$ and a finite width shear
  layer. The top row is for $\epsilon = 1$, the second row for
  $\epsilon = .25$, the third row for $\epsilon = .1$, the fourth row for $\epsilon
  = .01$, and the bottom row for $\epsilon = 0$. The first column
  shows the real part of $\varphi$, the second the imaginary part of
  $\varphi$ and the third the imaginary part of $\omega$. The solid
  line corresponds to the upper branch and the dashed line to the lower
  branch. In the top row, the upper and lower branches have the same
  growth rate. The arrows in panel (o) show the wavenumbers
  corresponding to trapped wavemodes (\S \ref{physicalintuition}).}
\label{growthfinite}
\end{figure}

\subsection{Numerical Verification}

To independently verify our dispersion relation (\ref{dis_hyp}), we
ran direct hydrodynamical simulations of shear layers using
the Godunov code Athena \citep{Stoneetal} and compared the growth rate
of the fastest growing mode predicted by our dispersion relation to that
obtained in the simulations. We initialized a flow along the $y$-direction
using the setup described by Equations (\ref{vsetup}), and for all of
the simulations we used $\gamma = 5/3$, $\delta_{BL}=1$ for the half width of
the shear layer, and periodic boundary
conditions in the $y$-direction; Table \ref{flowpartable}
summarizes the simulation specific parameters. To seed the
instability, we initialized random perturbations to $v_x$ of magnitude
$10^{-6}$ in the region $x > -1$. We found that for the $M=5$ runs, we
needed a high resolution in the $x$-direction to get converged estimates for the
growth rates. Figure \ref{athenatest} shows $\rho v_x^2/2$ averaged over box as a
function of time. The
solid lines show simulation results, and the dashed lines
show the predictions from considering the fastest growing mode for the
upper branch, using the dispersion relation (\ref{dis_hyp}). There is
good agreement between the two, especially for the limiting cases of
$\epsilon = 0$ and $\epsilon = 1$.

\begin{table}[!h]
\begin{center}
\begin{tabular}{|c|c|c|c|c|c|c|}
\hline
Label & $M$ & $\epsilon$ & $x$-range & $y$-range & $N_x \times
N_y$ & (BC-$x1$, BC-$x2$) \\
\hline
A & 5 & 1 & (-4,4) & (-16,16) & $8192 \times 8192$ & (outflow, outflow) \\
B & 5 & .25 & (-4,4) & (-16,16) & $8192 \times 8192$ & (outflow, outflow) \\
C & 5 & .1 & (-4,4) & (-16,16) & $8192 \times 8192$ & (outflow, outflow) \\
D & 5 & 0 & (-1,4) & (-16,16) & $2048 \times 2048$ & (reflecting, outflow) \\
E & .05 & 1 & (-4,4) & (-8,8) & $512 \times 512$ & (reflecting, reflecting) \\
\hline
\end{tabular}
\end{center}
\caption{Parameters for the simulations with Athena. $M$ -- Mach
  number, $\epsilon$ -- density ratio, $x,y$-range -- box size in
  $x,y$-direction, $N_x \times N_y$ -- number of cells in
  $x,y$-direction, (BC-$x1$, BC-$x2$) -- lower and upper boundary
  conditions in $x$ direction.}
\label{flowpartable}
\end{table}

\begin{figure}[!h]
  \centering
  \subfigure[]{\includegraphics[width=.6\textwidth]{athenatest}}
  \caption{Comparison of the analytically derived growth rates (given
  by the slopes of the dashed
  lines) to the ones obtained using Athena (solid lines). The curves A, B,
  C, D correspond to $\epsilon = 1$, $.25$, $.1$, and $0$, respectively
  (see text for
  further simulation details) and have been offset vertically for
  clarity.
  }
\label{athenatest}
\end{figure}

\subsection{Physical Intuition}
\label{physicalintuition}
Panels, (a), (b), and (c) of Figure \ref{eigenfig} show the spatial structure of
$v_x$ for simulations A, E \& D during
the linear stage of the instability. The
fastest growing modes in these three cases are very different, illustrating the
different instability mechanisms which operate for the $M \gg M_\text{crit}$ vs. $M
\ll M_\text{crit}$ cases and also for $\epsilon = 1$ vs. $\epsilon=0$. The
case $M \ll M_\text{crit}$, in panel (b) has been extensively studied
(e.g. \citep{Vallis}), so we will not consider it here..

\begin{figure}[!h]
  \centering
  \subfigure[]{\includegraphics[width=.49\textwidth]{M10_eigenfig}}
  \subfigure[]{\includegraphics[width=.49\textwidth]{M1_eigenfig}}
  \subfigure[]{\includegraphics[width=.75\textwidth]{M10eps0_eigenfig}}
  \caption{Panels (a), (b), and (c) show $v_x$ during the linear stage of the
  instability for simulations A, E, \& D respectively and have been
  rotated by 90 degrees.}
\label{eigenfig}
\end{figure}

For the $M \gg M_\text{crit}$ case with $\epsilon = 1$, the instability is
  caused by a {\it radiation mechanism}. This mechanism has already
  been discussed by \citet{Glatzel}, so we only mention it here
  briefly. Each mode of the dispersion relation can be associated with
  a pseudo-energy that is conserved. If a mode has a negative
  pseudo-energy, then radiation of energy away from the boundary layer
  region will cause the pseudo-energy to become even more negative,
  amplifying the mode and leading to instability. The radiation
  mechanism is responsible for the smooth, broad hump in panels
  (b) and (c) of Figure \ref{growthfinite}.

For the case of $M \gg M_\text{crit}$ and $\epsilon = 0$, the instability
  mechanism is quite different. Rather than a broad hump, panels (n)
  and (o) of Figure \ref{growthfinite} exhibit sharp localized
  peaks. The physical cause of these peaks is due to
  over-reflection of modes that become trapped between the lower
  boundary and the critical layer.
  Thus, we shall call this the {\it over-reflection mechanism}.

  The over-reflection mechanism is discussed in \citet{NGG}, who
  considered a rotating shear layer adjacent to a reflecting
  wall. They explained the instability in the $\epsilon = 0$ case in terms
  of the leaking of
  action density past the corotation radius, which is the analog of the
  critical layer for a rotating system. Like the pseudo-energy of
  \citet{Glatzel}, the action density is a conserved
  quantity, and instability occurs because a
  wavemode that is trapped between the critical layer and the wall has
  a negative action density, whereas positive action leaks out past
  the critical layer. As a result, the amplitude of the trapped
  wavemode becomes even more negative and grows with time. Put in this way, it is
  clear that the over-reflection and radiation mechanisms are
  related. However, they lead to a very different structure for the
  dispersion relation, as evidenced by panels (n) and (o) vs. panels
  (b) and (c) of Figure \ref{growthfinite}. Thus, we prefer to regard them as
  separate mechanisms, but point out that in both cases, instability is
  ultimately caused by radiation emitted from the boundary layer region.

  We now elucidate some of the properties of the modes that are trapped between
  the critical layer and the lower reflecting boundary using a simple
  model. Consider a sonic mode that is trapped between the lower boundary and the
  critical layer. If the mode has a well-defined
  phase velocity, then by performing a velocity boost one can
  transform into a frame in which the wavefronts are stationary. Let us
  work in this frame, since it makes the explanations more clear. In
  order for the wavefronts to be stationary, the equation for a wavefront is
  \begin{eqnarray}
  \label{wavefronteq}
  \frac{d y}{d x} = \pm \sqrt{\mathcal{M}(x)^2 - 1}.
  \end{eqnarray}
  Here, $\mathcal{M}(x)$ is the Mach number as a function of $x$, and
  the positive sign is for waves propagating in the $-x$ direction, while
  the negative sign is for waves propagating in the $+x$
  direction. Equation (\ref{wavefronteq}) is obtained from the fact
  that the wavefront propagates at the speed of sound, and in order
  for it to be stationary, the angle
  that the wavefront forms with the $y$-direction obeys $\sin(\theta) =
  1/\mathcal{M}(x)$.

  Consider now the upper branch from \S \ref{ratiomach}, for which the
  critical layer is at
  $x_{c} = \delta_{BL}(M-1)/M$. At the critical layer, we set $dy/dx =
  0$, which is just to say that an upward traveling sound wave is
  reflected there. It immediately follows that inside the region of shear
  \begin{eqnarray}
  \label{Mxeq}
  \mathcal{M}(x) =  -M\left(1-\frac{x}{\delta_{BL}}\right), \ -\delta_{BL} \le x \le \delta_{BL}. \\
  \end{eqnarray}
  Thus, the frame in which sound waves are stationary and reflect at
  the critical layer is boosted by $-M$ relative to the
  frame we defined in Equations (\ref{vsetup}).

  We now consider the fate of a wavepacket that is trapped between the lower
  boundary and the critical layer. Figure \ref{wavefrontfig} shows a
  set of wavefronts derived by integrating Equation
  (\ref{wavefronteq}) with red segments corresponding to propagation in
  the $+x$ direction and blue segments corresponding to propagation
  in the $-x$ direction. Consider a localized wavepacket at point A in
  Figure \ref{wavefrontfig}a. The wavevector of the wavepacket is
  oriented perpendicular to the wavefront, and as the wavepacket
  propagates towards point B, its wavevector is rotated by the shear.
  When the wavepacket reaches the critical layer at point C, it is
  reflected back towards point D, and the reflected wavepacket has a
  higher amplitude than the incident one (see \citet{NGG} for the case
  of a rotating shear layer). After
  reflecting off the critical layer, the wavepacket propagates
  downward to point D, its wavevector continuing to be
  rotated in the same sense as before due to the shear. Finally the
  wavevector reaches the perfectly
  reflecting boundary at point E, whereupon the $x$-component of its
  wavevector is reflected and the cycle begins anew. It is thus clear
  that the repeating cycle A-E will lead to exponential amplification
  of the wavepacket, due to over-reflection at the
  critical layer.

\begin{figure}[!h]
  \centering
  \subfigure[]{\includegraphics[width=.49\textwidth]{n1mode}}
  \subfigure[]{\includegraphics[width=.49\textwidth]{n2mode}}
  \caption{The $n=1$ and $n=2$ trapped modes for $M=5$ and
  $\delta_{BL} = 1$. The lower black line denotes the solid
  lower boundary, the upper black line denotes the top of the shear
  layer, and the dashed black line denotes the critical layer. The red
  segments correspond to propagation in the $+x$ direction and the
  blue segments correspond to propagation in the $-x$
  direction.}
\label{wavefrontfig}
\end{figure}

  The distance in $y$ between points A and E is the wavelength of the
  longest wavelength mode, $\lambda_{\text{max},y}$,
  that can be trapped between the lower boundary and the critical
  layer. A trapped mode can have a wavelength shorter than
  $\lambda_{\text{max},y}$, as long as its wavelength satisfies the
  relation $\lambda_y = \lambda_{\text{max},y}/n$ where $n$ is an
  integer. Thus, the relation for the wavenumber of the $n$-th trapped mode in
  the general case is
  \begin{eqnarray}
  k_{n,y} = \frac{2 \pi n}{\lambda_\text{max,y}}.
  \end{eqnarray}
  We indicate using arrows the values of $k_{n,y}$ for the first six
  trapped modes in panel (o) of Figure \ref{growthfinite}, and it is
  clear that they agree well with the locations of the peaks in the
  dispersion relation, which lends support for the
  over-reflection argument. We also plot the $n=2$ case in
  Figure \ref{wavefrontfig}b and point out the similarity between the
  analytical curves in Figure \ref{wavefrontfig} and the shape of the wavefronts from
  simulation D in Figure \ref{eigenfig}c.

  We have discussed the radiation mechanism for $\epsilon=1$ and
  the over-reflection mechanism for $\epsilon=0$. We now discuss what
  happens in the more general case of $0 < \epsilon < 1$. As we can see
  from Figure \ref{growthfinite}, {\it both} mechanisms operate
  simultaneously for $0 < \epsilon < 1$. In this case, there is partial
  reflection at the lower boundary, and as $\epsilon$ goes from one to zero, less
  and less of the radiation can escape from the boundary layer region,
  so the radiation mechanism becomes weaker and weaker. This is
  evidenced by the decreasing size of the bump in the second and third
  columns of
  Figure \ref{growthfinite} as $\epsilon$ goes to zero. For $\epsilon = 0$, the
  bump disappears entirely, and the radiation mechanism no longer
  operates. On the other hand, as $\epsilon$ goes from one to zero the reflection
  mechanism becomes stronger and stronger, since more and more of the
  energy is reflected at the lower boundary with total reflection at
  $\epsilon = 0$. Interestingly, there are still some small-scale wiggles in the
  dispersion relation, even for $\epsilon = 1$. This is because
  even if the density is everywhere uniform, radiation can still partially reflect off the
  discontinuity in the velocity derivative at $x = -1$ \citep{Glatzel}. However,
  reflections off a discontinuity in the velocity derivative are quite weak, so
  the wiggles are small, and the radiation mechanism is dominant.

\section{Discussion and Conclusions}
\label{discon}

We have studied supersonic shear instabilities that could drive the
turbulence in the BLs of stars for which the disk is undisrupted by a
magnetic field. Our study is aimed mainly at identifying the
  instabilities that lead to the formation of the BL when the disk
  just touches the surface of the star. The
main result of our work is the identification of two types of
instabilities that could operate in the BLs of such systems and had
not been previously discussed in this context.

The first is an
instability of a vortex
sheet at high Mach number caused by gravity. Although the vortex sheet
is stable above a critical Mach number, the addition of a small amount
of gravity destabilizes it. We have found that the eigenfrequencies of
the dispersion relation in the limit $G \rightarrow 0$ acquire a purely
imaginary term, which is proportional to the small parameter $G$. This
has the effect of making the
upper branch unstable. We now consider whether the instability of the
upper branch in the
limit $G \rightarrow 0$ is likely to be relevant during the initiation of
the BL, when one may expect the vortex sheet approximation to be valid.
Redimensionalizing Equation (\ref{uvp1}), we obtain
\begin{eqnarray}
\label{gfreq}
\omega_1 \sim \Omega_K \epsilon^{1/2} i,
\end{eqnarray}
where $\Omega_K$ is the Keplerian velocity at the surface of the star.
During the initiation of the BL, we
expect disk material to be less dense than stellar material, which
means that $\epsilon \lesssim 1$. It then follows from Equation
(\ref{gfreq}) that the characteristic growth rate of the instability
is $\lesssim \Omega_K$ and is independent of
the wavenumber.
\begin{comment}
The resulting
instability has a characteristic growth rate $\lesssim \Omega_K$.
  It is counterintuitive that gravity can destabilize the vortex
  sheet,
  since in the incompressible limit, $M \ll 1$, gravity stabilizes the
  KH instability, assuming the equilibrium state is stably
  stratified. Thus, we take some time to discuss this point on a physical
  level.

To simplify matters, we take the case of a stably stratified
  vortex sheet and
  consider the force per unit area across the interface between the
  upper and lower fluids. In the incompressible
  case ($M \ll 1$), the
  perturbation to the interface ($\delta \xi$) and the pressure perturbation
are in phase and of opposite sign. Therefore, there is an
upward/downward force on the peaks/troughs of the interface due to the
pressure perturbation which acts as a destabilizing influence,
and a downward/upward force on the peaks/troughs of the interface due
to gravity (buoyancy force) which acts as a stabilizing
  influence. When the upward pull
due to pressure exceeds the downward pull due to buoyancy, instability
will result.

{\bf  The situation changes at high Mach number. For the vortex sheet
  without gravity, $\bar{k}$
  becomes imaginary above a critical Mach number (Equation
  (\ref{Mcritunstable})). This causes the pressure perturbation to be
  90 degrees out
  of phase with the interface perturbation, but in phase with the
  velocity perturbation. If we then add a small amount
  of gravity, the situation resembles
  pumping a swing, since gravity still acts on the peaks/troughs of
  the interface, but
  pressure (the agent pushing the swing) acts out of phase with the
  gravity. We mention, though, that this second analogy is only an
  approximate one, since we have not considered the effects of
  compressibility in detail. However, it does capture the essence of
  why gravity is no longer a purely stabilizing influence at high Mach number.

STILL NEED BETTER EXPLANATION...

}
\end{comment}

If the BL is radially thin, as would be expected
during the initiation phase, then the growth rate given by Equation
(\ref{gfreq}) is small compared to the shear rate $S \sim \Omega_K
R_*/\delta_{BL} \gg \Omega_K$. Considering that $S$ is the characteristic
growth rate for shear instabilities, if there are other
mechanisms for instability with growth rates proportional to $S$, they
will quickly become dominant.

In particular, we have demonstrated that the growth rate of the
sonic instability of a finite width shear layer at high Mach
number is proportional to $S$. The sonic instability is similar in
nature to the Papaloizou-Pringle instability in that both are global
instabilities and cannot be derived from a local analysis. There are
two destabilizing mechanisms for the sonic instability. The first
corresponds to emission of radiation and gives instability over a broad
range of wavenumbers. The second corresponds to over-reflection of
trapped modes and results in sharply peaked, disconnected regions of instability in
$k$-space. Because sonic instabilities
operate on a much faster
timescale than the gravity mechanism for a vortex sheet, this makes
them
an appealing candidate for the initial stages of boundary
layer formation, when one might expect large shears to be
present.

We mention that \citet{AlexakisMiles} have
considered the Miles instability \citep{Mileswave} in the context of
boundary layer
formation. Specifically, they invoked the
Miles instability to generate mixing of WD and stellar
material and explain the enrichment in heavy elements observed in nova
explosions. The Miles instability was proposed as a mechanism for
generating waves over
water at low wind speeds and operates through a
resonant interaction between the wind and the water wave. Due to this
interaction, a
component of the pressure perturbation is created that is in phase with the
velocity of the air-water interface, and much like pumping a swing, swings
up the interface to large amplitudes. However, the Miles instability
was initially introduced as a way of explaining the formation of waves
on water at weak wind speeds, for which the air-water interface is
stable to the KH instability. During boundary layer formation,
however, large shears are generated, so we are not in the weak wind
regime, and the Miles instability is likely to be swamped by sonic
instabilities.

There are two main astrophysical implications of our findings. First, we
demonstrate that the initiation of the BL is likely to take only very
short amount time after the first material from the disk arrives to the
stellar surface. This is because the sonic instabilities that we explored
in this paper have an extremely short growth rate, which is very weakly
dependent on the density contrast between the disk and stellar
material.\footnote{This statement is not true for the supersonic KH
instability with gravity investigated in \S \ref{compress_grav},
see Equation (\ref{uvp1}).} Thus, mixing of the two fluids starts
almost immediately after they come into contact.

Second, given the efficiency with which the sonic instability operates, it is
likely that it may play important role also in more developed phases of
the BL evolution. As long as the BL possesses some effective "boundaries"
(e.g. sharp changes in the velocity of density behavior) and the gas flow
within it is supersonic the purely hydrodynamic sonic instabilities are
going to operate in it potentially providing means for continuing mixing
and angular momentum transport inside the layer. Future numerical
calculations capable of following the nonlinear development of sonic
instabilities should be able to address this issue.

\acknowledgements

We are grateful to Jim Stone for useful discussions.
The financial support for this work is provided by
the Sloan Foundation and NASA grant NNX08AH87G.

redpajama_set_name: RedPajamaArXiv