Linear dynamics of weakly viscous accretion disks: a disk analog of TollmienSchlichting waves
O. M. Umurhan^{1,2,3}  G. Shaviv^{4,5}
1  Astronomy Unit, School of Mathematical Sciences, Queen Mary
University of London, London E1 4NS, UK
2 
Department of Geophysics and Space Sciences, TelAviv University,
TelAviv, Israel
3 
Astronomy Department,
City College of San Francisco, San Francisco, CA 94112, USA
4 
Department of Physics, TechnionIsrael Institute
of
Technology, 32 000 Haifa, Israel
5 
Institute of Theoretical Astrophysics,
University of Heidelberg, 69120 Heidelberg, Germany
Received 19 February 2008 / Accepted 14 November 2008
Abstract
This paper discusses new perspectives and approaches to
the problem of disk dynamics
where, in this study, we focus on the
effects of the viscous instabilities influenced by
boundary effects.
The Boussinesq approximation of the viscous large shearing box
equations is analyzed in which the azimuthal length scale of the
disturbance is much larger than the radial and vertical scales.
We examine the stability of
a nonaxisymmetric potential vorticity mode, i.e. a PVanomaly.
in a configuration in which buoyant convection
and the stratorotational instability do not to operate.
We consider a series of boundary conditions that show
the PVanomaly to be unstable both on finite and semiinfinite radial domains.
We find these conditions lead to an instability that is the disk analog of
TollmienSchlichting waves. When the viscosity is weak, evidence of the instability
is most pronounced by the emergence of a vortex sheet at the critical layer
located away from the boundary where the instability is generated.
For some boundary conditions, a necessary criterion
for the onset of instability for vertical wavelengths
that are a sizable fraction of the layer's thickness and
when the viscosity is low is that
the appropriate Froude number of the flow be greater than one.
This instability persists if more realistic boundary conditions
are applied, although the criterion on the Froude number
is more complicated. The unstable waves studied here share qualitative
features to the instability seen in rotating Blasius boundary layers.
The implications of these results are discussed. An overall new strategy
for exploring and interpreting disk instability mechanisms
is also suggested.
Key words: instabilities  hydrodynamics  waves  stars: planetary systems: protoplanetary disks
1 Introduction
The magnetorotational instability (Balbus 2003, MRI hereinafter) is generally considered to be the leading candidate explaining the source of enhanced transport observed for disk systems. Three conditions are required for its operation: the concurrent presence of rotation and shear, a primordial (no matter how small) magnetic field, and sufficient ionization of the fluid so that the gas is in the MHD regime. Consequently, it is natural to ask the question: what happens in an accretion disk if one or more of these conditions are not satisfied?
Supposing that there are disks in which the MRI or any other MHD (dynamo) mechanism is either weakly operating or entirely absent: what else can drive activity, possibly even leading to turbulence? Attempts to answer this question include, but are not limited to, defects in the Keplerian profile (Li et al. 2000), baroclinic instabilities (Klahr & Bodenheimer 2003; Johnson & Gammie 2005; Petersen et al. 2007a,2007b), transient growth and sustained subcritical dynamics (Richard & Zahn 1999; Iounnou & Kakouris 2000, Chagelishvili et al. 2003; Tevzadze et al. 2003; Yecko 2004; Umurhan & Regev 2004; Umurhan et al. 2006; Mukhopadhyay et al. 2004; Barranco & Marcus 2005; Lesur & Longaretti 2005; Lithwick 2007; and see the experiments of Richard 2001; and Ji et al. 2006) and, unsteady wave dynamics like the PapalouizouPringle instability (Papalouizou & Pringle 1984), hereinafter ``PPI'', and the stratorotational instability (Dubrulle et al. 2005; Shalybkov & Rüdiger 2005; Umurhan 2006; Brandenburg & Rüdiger 2006), hereinafter ``SRI''. The majority of these recent investigations (excepting Yecko 2004; Afshordi et al. 2005) have focused on strictly inviscid processes. Indeed, the classical approach to such questions is to investigate processes that might lead to turbulent transport by first establishing a mechanism of linear instability from the vantage point of purely inviscid (or nearly inviscid) flow.
Since astrophysical fluids have some effective viscosity  however weak it may be  we pose the question: could a weakly viscous flow in a sheared and rotating environment undergo an intrinsically viscous type of linear instability that nonlinearly saturates with significant amplitude? A few previous studies have addressed this question. Kato (1978) demonstrated that if a fluid's viscosity is a function of the state of the fluid then a disk can experience pulsational dynamics in a way similar in quality to stellar pulsation like that in the theory of Cepheid variables. Hereafter we shall refer to this effect as the viscous pulsational instability (VPI). Latter & Ogilvie (2006) reexamined the VPI by studying how even axisymmetric fmodes, in a shearing sheet environment, create fluctuating stresses that explicitly draw energy from the shear which leads to overstability. Kleiber & Glatzel (1999) have shown that accretion tori (in this case, ones which have a constant specific angular momentum distribution) can be linearly unstable above a minimum Reynolds number. Dubrulle et al. (2005) and Shalybkov & Rüdiger (2005) also report that the growth rate of the SRI may be, under certain conditions, enhanced by viscosity.
Before proceeding we present some remarks concerning the nature of the fluid state to be studied as well as a discussion about boundary conditions and a proposed alternative way to consider their uses.
1.1 Viscous considerations
The noncommutative nature of the NavierStokes equations in the limit of with the Euler equations is a long standing fact (see the discussion in Schlichting & Gersten 2001, p. 968). For instance, viscous stresses do not necessarily vanish on the boundaries of a viscous flow when the viscosity limits to zero. In turn this implies that one may not properly take the Navier Stokes equations and naïvely substitute zero viscosity to reach the inviscid limit (Lions 1993; Joseph 2003). Aside from the generic appearance of boundary layers, other effects can appear when a weak viscosity is included into problems of study. For example, in the nearly inviscid shallowwater theory of strong shear flows (Balmforth 1999), normalmodes can emerge out of a continuous spectrum when viscosity is introduced into the dynamics.
Another example of the subtleties inherent in viscous flow is the instability associated TollmienSchlichting waves (see discussion of TS waves in Schlichting & Gersten 2001; and Schmid & Henningson 2000)  which are traveling waves appearing in wallbounded flows that neither grow nor decay in the inviscid limit and become unstable when viscosity is included in the analysis (examples include plane Pouiselleflow and Blasius boundary layers). Similar to the process discussed by Balmforth (1999), the traveling mode becomes unstable through the interaction of an inviscid normal mode and a viscous normal mode  the latter of which exists only as a member of the continuous spectrum in the inviscid limit (Baines et al. 1996). Far from being considered mathematical oddities, TS waves appear to play a prominent role in the transition to turbulence in boundarylayer flows (for a recent summary see Drazin 1984). A situation studied which closely resembles the condition of an astrophysical disk is the modal and nonmodal response of a rotating Blasius boundary layer (Yecko & Rossi 2004) in which instability is promoted when the azimuthal scale of a disturbance is longer than its vertical scale.
In a general sense, because the governing equations are of a higher order in the viscous case, a new space of possible solutions emerges which are either absent or inactive in the inviscid case. Our main query is therefore: if rotationally supported flows are (linearly) wellbehaved in the exactly inviscid limit but the viscous flow shows some type of dynamically significant behaviour  even as the inviscid limit is approached  then might it be misleading to test stability only of exactly (Re = ) inviscid flows? Perhaps the subtle nature of disks is linked to this feature. That which best summarizes this perspective is the quote attributed to the atmospheric dynamicist Eady where he is purported to have said, ``It is not the process of linearization that limits insight. It is the nature of the state we choose to linearize about'', (Bayley et al. 1988).
To prospect for an instability mechanism that might lead to sustained unsteady behavior by assuming a turbulent viscosity model a priori might seem contradictory at first. But, given the difference in behavior known to exist (in other problems) between nearly inviscid and exactly inviscid models it is therefore mandatory to clarify these differences within the context of an astrophysical disk too. A conjecture that such an investigation could address is the following. It is reasonable to suppose that disks are continuously fed with a turbulent flow field either by infall or some mass transfer processes. Could this turbulent flow field undergo a secondary transition into another dynamical state (possibly turbulence of a different stripe) due to the turbulently enhanced viscosity? Suggestions which hark on these lines of thought are found in Doering et al. (2000), Kersalé et al. (2004) and more directly in terms of secondary transitions induced by Ekmann flow as suggested in Lesur & Longaretti (2005). Overstability driven by material fluctuations in the turbulent stresses (i.e. the VPI, Kato 1978; and Latter & Ogilvie 2006) are also candidates for such secondary transitions.
1.2 Interpreting boundary conditions and their effects on dynamics
Much attention has been devoted to evaluating processes which are intrinisic in some way to the fluid  meaning to say that it is assumed that it is more valuable to study those mechanical processes which are minimally sensitive to the boundary conditions imposed and maximally ``emerging'' out of something essential about the fluid and its basic state. The MRI is an example of this as well as other more basic fluid dynamical instabilities such as the RayleighTaylor and RayleighBenard instabilities. We wonder if this approach to the question of disk turbulence may be selflimiting given that many fluid instabilities which lead to some form of turbulence in terrestrial flows are driven in large part by the boundary conditions of the system (e.g. TS waves and turbulent transition). An alternative way is to view boundary conditions as a filter for certain solutions or as a tool to classify solutions. Kersalé et al. (2004) adopt a similar philosophy by studying the linear response of an incompressible fluid in a TaylorCouette type of cylindrical flow subject to a variety of boundary conditions. Of course, the TaylorCouette setup and the boundary conditions they test are not what one would expect in a terrestrial apparatus or experiment, however, if taken as a metaphor for a disk environment then this sort of exploration allows one to test, evaluate, understand and eventually categorize the dynamical response of a fluid as a function of differing boundary conditions.
The inviscid PPI and SRI are examples of linear instabilities which come about due to the imposition of arguably artificial boundary conditions on inner and outer walls of a model disk system. On the other hand, these results may be intepreted in terms of the HayashiYoung criterion for wave instability which states that a physically separated wave pair may promote linear instability if the waves counterpropagate with respect to each other with nearly the same wavespeed and if the waves have an actionatdistance effect upon each other^{} (Hayashi & Young 1987; and see also Sakai 1989; Baines & Mitsudera 1994). Indeed Goldreich et al. (1985)^{} point out that the PPI may be viewed as a process resulting from the interaction of a pair of edgewaves mutually interacting with each other across a waveevanescent region^{}. The SRI may also be similarly rationalized (Umurhan 2008).
Thus although the counterpropagating edgewaves responsible for the PPI and SRI are understood to result from the use of unrealistic boundary conditions, it is certainly not ruled out that the general counterpropagating wave mechanism could be at work in real disks. The linear instability of disks with two or more (potential) vorticity defects (e.g. Li et al. 2000) could be interpreted as an instance of this process.
1.3 An overview of the findings in this study
From a systematic asymptotic scaling analysis we derive in Sect. 2 and Appendix A the equations appropriate to a box section of a viscous shearing accretion disk (assuming an viscosity formalism) by exploiting the smallness of the parameter which assesses the ratio of the soundspeed to the rotation speed measured at some radial point of a circumstellar disk. We refer to this model as the Viscous Large Shearing Box (VLSB) and these equations have appeared before (cf. Latter & Ogilvie 2006). We are reminded that while the velocity fluctuations in the shearing box are an order smaller than prevailing rotational (``Keplerian'') velocities, the steady accretion velocities implied by the alpha viscosity model are an order smaller than the same disk rotational velocity.
We consider the fate of a nonaxisymmetric potential vorticity disturbance (or simply ``PVanomaly'') subject to varying boundary conditions. Accordingly, in Sects. 3 and 4 the VLSB are analyzed in the limit where the perturbation's azimuthal length scale is asymptotically larger than its radial and vertical scales (i.e. the quasihydrostatic semigeostrophic limit, Umurhan 2006, QHSG for short). Additionally we assume that the vertical component of gravity and entropy gradient are constants.
In Sect. 4.1 we formulate energy integrals of this reduced system in order to better understand what can contribute to destabilizing the PVanomaly . The energy budget is characterized by a ReynoldsOrr type of equation whose sources and sinks are given by the energy which the PVanomaly can extract from the shear, receive from the boundaries or lose due to dissipation.
In the spirit of Kersalé et al. (2004) we analyze the response of the PVanomaly subject to a controlled array of boundary conditions. In discussing boundary conditions we refer to the boundary closest to central object as starside as opposed to the side furthest away from the object to which refer to as farside. We consider the dynamics as occurring on both a semiinfinite domain (farside at infinity) and on a finite domain. Below we summarize the main findings. Note that we have made sure to consider boundary conditions which filter out the SRI or PPI instabilities.
In Sect. 4.2 an asymptotic analysis is done for the limit where the scaled turbulent viscosity parameter (defined in the text as ) is small. We find instability if the Froude number of the flow exceeds 1 for modestly large vertical wavenumber. Additionally, the PVanomoly interacts with a critical layer of the flow creating a potential vorticity sheet sheet whose radial extent is the size of the vertical extent of the disk. This analysis illustrates how an inviscid edgewave phenomenon (due to the nonormal flow starside boundary condition) becomes unstable when viscosity is included. Most importantly is that the instability is driven in part by the injection of energy through the boundary.
We consider in Sect. 4.3 finite domain disturbances of the PVanomaly and let the viscosity parameter be an order 1 quantity. The fourth order normal mode problem requires us to appeal to numerical computational methods for solutions. We impose on the farside boundary that both the disturbance pressures and PVanomalies vanish. At the starside boundary we require that there be nonormal flow there (as above). The remaining starside condition takes on four possibilities: (a) the flow is rigidly coupled at the wall; (b) the perturbations are stressfree; (c) the PVanomaly is zero; (d) the PVanomaly gradient is zero. The first two of these conditions are physically realistic. The latter two offer a means to consider the the effect of energy injection (or lack thereof) through the boundaries and to compare with the analytical analysis. For rigid and stressfree boundary conditions we see clear indications of a TollmienSchlichting type of instability, similar to the instability of rotating Blasius boundary layers (Yecko & Rossi 2004) and the energy budget of the disturbances show that this process does not draw upon energy across the boundaries.
2 Viscous Large Shearing Box and its QHSG approximation
In Appendix A we consider a box section of an disk
centered about its midplane and at a distance R_{0} from the
central object. If the disk is cold, then it means that
the quantity defined by the ratio of the typical value
of the local midplane disk soundspeed, c_{s},
to the local Keplerian velocity, V_{0},
is less than 1 by some substantial amount: protoplanetary disks, for example, are believed to have an . Using now familiar scaling arguments and exploiting the smallness of we derive from the full equations of motion in cylindrical coordinates appropriate equations of motion in what we refer to as the Viscous LargeShearing Box (VLSB for short). The tactics and procedures behind this effort are the same ones employed in the derivation of the LargeShearing Box (LSB) (Goldreich & LyndenBell 1965; Umurhan & Regev 2004) however, the viscous stresses are included. We have,
in which the total
entropy is defined by
and is the ratio of the specific heats at constant pressure to the specific heat at constant volume. The vertical component of gravity is dependent on z
All primed quantities are perturbations about the basic flow. The viscous stresses are
where
(10) 
The above equations are nondimensional. Time is scaled by the local rotation time of the box. All lengths are scaled according to a length which is comparable to the disk thickness (see Appendix A). Pressures are scaled according to the product of the local midplane soundspeed and density, which is in turn based on some fiducial characteristic temperature scale. For further details see Umurhan & Regev (2004). x represents the radial (shearwise) coordinate of the SB while y is the azimuthal (streamwise) and z is the vertical coordinate (normal to the disk midplane). The velocity components, i.e. , are for the radial, azimuthal and vertical directions. It is important to keep in mind that these flow variables represent perturbations about the steady Keplerian flow. , sometimes also referred to as the Coriolis parameter, is 1 in these nondimensionalized units, meaning to say because time has been scaled according to the dimensional value of the rotation rate at R_{0}, i.e. , the local Coriolis parameter formally is equal to one. We retain this symbol in order to flag the Coriolis effects in this calculation. The local shear gradient is defined to be
(11) 
in which is the full disk rotation rate. For Keplerian disks the value of q is 3/2. The local Keplerian flow is represented here by a linear shear in the azimuthal direction, i.e. .
The steady state quantities are denoted with index b and in particular
we assume that pressure (p_{b}) and density ()
profiles satisfy the hydrostatic balance relationship
All corrections to this equations are of a higher order and ignored here. The expression responsible for the VPI can be identified as fluctuating viscosity parameter in (8), term .
We call to attention that the accretion and meridional velocities characterizing disks (Kluzniak & Kita 1999) do not appear in the VLSB equation set (see Appendix A). By definition, turbulent disks exhibit accretion velocities as they are the natural consequence of equations describing global dynamics. The shear velocities are in fact quite complex as Kluzniak & Kita (1999) showed for the particular case of an ShakuraSunyaev type of disk (Shakura & Sunyaev 1973). The radial velocities in steady state are found to be sheared in the vertical direction and, as well, there exists a vertical component to the flow with both radial and vertical dependence (the meridional flow). However the scaling arguments implemented to reach these ``shearing box'' equations, especially the relative scaling relationships between the dynamical velocities and the accretion scalings, show that the influence of the steady accretion rate appears at higher orders in the expansion procedure. In other words, dynamical perturbations on the scale of the box do not feel the effects of steady accretion and meridional flow  they only feel the effects of the steady Keplerian shear. The scaling analysis also shows that the viscosity (which is the driver of the accretion flow) does influence the dynamics at these scales and is the reason why it appears in these equations.
3 The quasihydrostatic semigeostrophic approximation of the VLSB equations
The equations of motion may be simplified for further analysis by implementing the quasihydrostatic semigeostrophic (QHSG) scaling arguments used in Umurhan (2006). The QHSG is useful in its ability to expose the essential mathematical features of the inviscidSRI (Dubrulle et al. 2005; Umurhan 2006).
We suppose that the azimuthal scales of motion are much larger
than the radial or vertical scales. We measure this with the
small parameter .
In order to maintain
asymptotic validity we
assume the following orderings
(13) 
Thus we suppose that the following operations upon dynamical quantities scale accordingly as
Then we suppose that the radial and vertical velocities are correspondingly smaller than the azimuthal velocities by this same scale, in other words
These scalings say then that
These scalings will make it easy to follow waves propagating with respect to the background Keplerian flow velocity. Therefore, the temporal dependence should also scale by the scaling appropriate to . It follows that
Furthermore we say that the density, pressure (and by implication, the entropy) fluctuations are all order 1, that is
The new issue that must be addressed here is to suggest a scaling that brings in the viscous terms at the lowest nontrivial order. To this end setting achieves this goal and we shall formally write . In sum, then, to lowest order we have the following reduced set:
where we have introduced the basic state entropy
and
its dynamically varying counterpart
which are defined by
Only the azimuthal direction stress component survives at lowest order due to this scaling argument,
Although we have invoked scaling argments leading to the above sets of equations we have not formally rewritten all of the variables to signify these assumptions as it is our desire to preserve the transparency of the subsequent presentation. Note that effect responsible for the VPI survives this scaling argument as it appears in (20) as the term .
4 Boussinesq simplification, assumptions and lLinearized dynamics
In Umurhan (2006) it was shown that the QHSG approximation of Boussinesq disk models recovers the linearized hydrodynamic behavior contained therein for concurrent small values of the azimuthal wavenumber and wavespeed. It was further demonstrated that the dynamics contained in the QHSG approximation of the LSB is faithfully represented if one considers instead the equivalent incompressible Boussinesq (Spiegel & Veronis 1960) version of QHSG approximated LSB equations. Applying this sequenced reasoning to the linearized version of (1418) gives,
where we have dropped all primes from the velocity quantities
and explicitly set
to its value of 1.
In the usual Boussinesq approximation,
density fluctuations are dynamically significant
when coupled to gravity. In these circumstances is replaced by .
The nondimensionalized temperature quantity
and
its associated steady state temperature field T_{b}(z) are characterized
by the (linearized) conservation relation,
For clarity we reexpress this thermal quantity in terms of given by ^{}. In the Boussinesq approximation fluctuating density variables influence the dynamics when coupled to gravity. Forthwith, and P_{b} are taken to be constant (set to 1) and, furthermore, is absorbed into the fluctuating pressure leading to defining the enthalpy . All unprimed velocities (i.e. u,v,w) are now understood to represent linearized disturbances. The nondimensionalized BruntV is l frequency, N, emerges in the equation for the perturbation temperature field (25) and is given by
Throughout this study Nis taken to be real (buoyantly stable). The azimuthal stress is
(27) 
The term that gives rise to the VPI appears in the above as . We proceed further by (i) operating on (23) with ; (ii) operating on (22) with ; (iii) and subtracting the results to reveal
(28) 
where the incompressibility condition was used in writing the first term on the RHS of this expression. Multiplying (25) by (2q)/N^{2}followed by operating on the result with gives
(29) 
Adding these two equations together and making use of the relationships (22) and (24) yields the following single equation for :
To showcase the elements in the above we have written in order to remind ourselves that the parameter is like the inverse of a Reynolds Number^{}. We have also defined the epicyclic Froude number
which is in general a function of the vertical coordinate zand vanishes on the symmetry axis driving the local Froude number to very large values. As we shall see, determines the character of the solutions that emerge. Umurhan (2006) demonstrates by a direct comparison of calculations that the onset of the SRI (in the QHSG limit of the LSB equations explored there) is reasonably well captured analytically when N^{2} and g are assumed to be constants. Guided by these previous results as well as similar use in a series of other studies (e.g. Tevzadze et al. 2004; Bodo et al. 2007a,2007b) we shall assume
implying that is a constant as well, and we restrict our considerations to vertically periodic solutions. We shall return to this matter in the Discussion.
In this asymptotic theory
corresponds to a the perturbed potential vorticity (Tevzadze et al. 2004; Umurhan 2006), also known as the potential vorticity anomaly (Hoskins et al. 1985) and also called vortensity in the astrophysical literature (e.g. Klahr & Bodenheimer 2003). Throughout the rest of this work we will interchangeably use the terms potential vorticity perturbation, potential vorticity disturbance and potential vorticity anomaly (i.e. PVanomaly). The PVanomaly Q relates to a vorticity pointing in the vertical direction.
We consider travelling wave normal mode solutions to (30):
the modes are assumed
to be
azimuthally periodic on scale L_{y} and vertically periodic
on scale L_{z},
(32) 
where the azimuthal wavevector k can be any positive number while the vertical wavevector is any real number. The wavespeed c can be complex: when there is growth of the wave. Furthermore we define
to be a wavenumber scaling of the viscous parameter. Thus the equation governing the structure function becomes
The Froudewavenumber is defined as . (34) will be the fundamental equation of study. Furthermore expressed in this form, the azimuthal and radial velocities are
with the vertical velocity following from evaluating the incompressibility equation
at the boundary and using (34), i.e.
(37) 
Finally, the perturbed potential vorticity is . The remainder of the boundary conditions, namely in the x direction, will be stated in the following sections according to the problem being solved.
4.1 Energy integrals
It is instructive to develop global energy integrals
as such quantities aid in developing an interpretation of
the results in the following sections. We restrict our attention
to the energetics associated with the QHSGBoussinesq model set
(2221) keeping in mind
the assumptions we made about vertical and azimuthal periodicity
and the constancy of N^{2}.
The radial conditions are left arbitrary
and they will be dealt with accordingly in each of the subsequent subsections.
We proceed by defining the perturbation
thermomechanical energy density as
by multiplying (23) by w and (25) by and integrating over a domain which is periodic in the vertical and horizontal directions, and finite in the radial direction given by 0<x<x_{1}, where x_{1} is left arbitrary, we find
in which volume integrals are
integrated on volume with the volume element
and where
the surface term is
In writing we have made use of the periodicity conditions in the vertical and horizontal directions. The bracketed term as appearing means
where is the location of the outer boundary.
The energy integral (38) is the ReynoldsOrr Equation appropriate for this QHSG system. The energy E is composed of the baroclinic thermal term () plus a kinetic energy term (w^{2}), however, the kinetic energy term contains only the azimuthal velocity contribution because the vertical and horizontal velocity contributions are small by comparison in the scaled system of equations according to the QHSG approximation of the original set (see Sect. 3 and Umurhan 2006, for details). The term represents the integrated losses and is comprised of viscous losses due to the azimuthal velocity since, by the same reasoning as above, the corresponding losses due to vertical and radial velocities (in the scaled system) are negligible. The perturbed viscous stress term, responsible for the VPI plays, a destabilizing role for these PVanomalies as it appears in the above as the the term proportional to 2qv^{2}. The total external stresses on the system is given by and is comprised of the surface integrated body pressure and the surface viscous stress  the latter of which is expressed only with the tangential stress due to w. Finally is the Reynolds stress due to the background shear state. This expression may be interpreted as accounting for the amount of energy perturbations extract from the background shear state. We note also that inspection shows that always, while the remaining terms and may be either positive or negative given the state of the perturbed flow or the boundary conditions employed.
4.2 Asymptotic theory on a semiinfinite domain
The following assumptions are made in order to proceed analytically: (i) the domain in the x direction lies between 0 and , thus we have ; (ii) we require that all quantities decay as ; (iii) there is no normalflow at x=0  this inner location is considered starside as it is closest to the central object; (iv) the viscosity is weak but finite, hence, we assume that ^{}; (v) even though there is viscosity in the problem we impose no particular stress condition at x=0 and let the fluid quantities be dictated by what emerges in the interior of the domain. This essentially means that flow stresses at the starside are allowed to adjust according to the dynamical response happening in the interior of the domain. In practice it translates to only enforcing the nonormal flow boundary condition^{}.
4.2.1 Expansions and outer solution
We assume the following expansions well aware that this
is a singular perturbation calculation because of the presence of the critical
layer (see below). The solution for the wavespeed c is assumed of the form
(41) 
and a similar series of the form
(42) 
We know aposteriori that the critical layer will generate a solution proportional to and this is why such a term appears in the above expression. Thus at we find that
(43) 
To this order the noflow boundary condition at x=0 amounts to
The solution to this equation which decays as is
where A_{0} is an arbitrary amplitude. Using this solution in boundary condition (44) amounts to selecting c_{0}, which is
This says that the wavespeed is real and positive, which means in this case that there will be a critical layer in the domain, i.e. at (see below).
At we find the equation
(47) 
Using (45) for
and dividing the equation by c_{0}x we find
the more transparent form
(48) 
In which we have defined the parameter . The delta function appearing with the arbitrary coefficient C_{1} is a formal device used to signal the presence of a critical layer in the flow (Case 1960; and more recently, Balmforth & Piccolo 2001; Balmforth et al. 2001). The coefficient will result from the critical layer analysis (see below) which follows the procedures found in similar investigations (e.g. Stewartson 1981; Balmforth & Piccolo 2001). In practice it means that we must separately develop solutions to on either side of x=x_{c}. Before doing so let us observe the reason why this series expansion fails near x_{c} by developing the solution to in the vicinity of this point. An indicial analysis shows that
(49) 
These expressions show, especially that for , that the solution begins to breakdown (i.e. break order) when the quantity xc_{0} starts to approach 0. This divergence must be controlled by considering a boundary layer calculation in and around the critical layer. Note that since the vertical vorticity is proportional to , the critical layer will appear as a vortex sheet. The formal presentation of this solution, including the region of validity and expression of the boundary conditions at this order is presented in Sect. B.1.
4.2.2 Critical layer calculation, matching, and growth rate and analysis
As we demonstrated above, the solutions begin to breakdown in the vicinity of the critical layer which are those places where xx_{c} begins to get small. We must therefore reexamine (34) in this zone and to this end we define a new inner coordinate as(50) 
According to this new coordinate (34) is reexpressed as
We intorduce a series expansion for the solution to (51) by writing
The remainder of this calculation including the matching of the inner and outer solutions and the determination of the growth rate c_{1} has been relegated to Appendix B.2. We note that the term proportional to is needed for matching purposes as the inner solution is extended out of the critical layer (for details see the full exposition in the Appendix). We find that the correction wavespeed obeys
The growth rate of this mode, i.e. Im(c_{1}), is proportional to
where we have restored the definition of in terms of Re, and where we have replaced k_{F} accordingly with . The term inside the square brackets is always greater than zero for q>0and limits to 3/2 as ^{}. Therefore there is growth when .
The instability emerges from the
inner boundary due to application of the nonormal flow condition
but it will be at the critical layer where evidence of it
appears in the form of a pronounced PVanomaly.
The radial width of this vortex zone is proportional
to
the size of the box and the amplitude will
be
times the leading order perturbation pressure field
(see the end of Appendix B.2). Restoring units and
recalling that the shearingbox has been scaled according to the
thermal scale height of the disk H , the radial
width of this vortex zone, ,
is
Further analysis of these solutions together with the ReynoldsOrr Eq. (38) and the definitions (3940) shows that to lowest order in Re^{1}
(55) 
so that
A general evaluation of the Reynolds stress term shows that it contributes first at order as well, i.e.
is the first order correction to the azimuthal/radial velocity correlations. Its details may be worked out but we leave it here in general form in order to make the following argument. The ReynoldsOrr equation evaluated for this set of boundary conditions is to lowest order in Re^{1},
plus a correction which is . In the case of instability, i.e. , it follows that . Given (54) taken together with (56) we find that under conditions of linear instability it also follows that
Let us reflect upon this for a moment: with these boundary conditions it is always the case that when . The combined action of the starside viscous stresses and the domain integrated viscous lossess still results in a net transfer of energy into the domain. Because the term within the square brackets in (57) is always positive, these conditions promote the type of corrections (at order Re^{1}) to the azimuthal and radial velocity profiles such that a positive global correlation between them emerges and, hence, resulting in when . We conclude that with these boundary conditions the instability observed is fed both by the energy entering the domain from the starside boundary as well as by the energy extracted by the shear due to the resulting order Re^{1}velocity profiles.
4.3 Finite radial domain investigations
In this section we consider the normalmode solutions
of (34) occurring on a finite radial domain
where
.
All solutions
are computed numerically using a NewtonRaphson scheme on a Chebyshev
grid of anywhere from 33 to 129 points  higher
resolution is needed for smaller values of .
All numerically generated
solutions are normalized so that
.
Because
this is a fourth order system we must specify four boundary conditions.
In all
of the following calculations two of the boundary conditions
will be that
(58) 
in other words, that the pressure fluctuation and PVanomaly are zero on the farside boundary. According to the normalmode PVanomaly, , the fixed pressure condition implies that at the far boundary . It therefore follows from (35) that the perturbed stress expression, , is zero there as well. A third condition will be that there is no normalflow on the starside boundary (as before), i.e.
(59) 
For the remaining starside boundary we shall explore four different conditions enumerated in the corresponding subsections below. The most physically plausible viscous starside condition is to set to zero either the azimuthal velocity fluctuation or the azimuthal stress fluctuation. We have also considered zero PVanomaly and zero PVanomaly gradient conditions. Although these conditions are less physically realistic, they are simpler to interpret in terms of the energy arguments developed in previous sections.
4.3.1 Rigid and stressfree starside boundary
These conditions translate to requiring (rigid) or (stressfree) at x=0. Note that the usage of the term ``stressfree'' to describe the boundary condition really refers to the perturbed part of the azimuthal flow as being stressfree. Inspection of Fig. 3 shows that there are regions in the  parameter plane in which there is instability when rigid or stressfree conditions are imposed on the starside boundary. Although the parameter range for instability is not as expansive as it is for the other lessrealistic boundary conditions explored (see both the previous section and below), the results suggest that a rough criterion for linear instability is . Furthermore, under these starside conditions we find that the ReynoldsOrr Eq. (38) simplifies tosince the zerostress and rigid conditions implies that . Thus we find that the boundary conditions are such that instability is not directly driven by the injection of energy into the domain due to the (perturbed) body stresses. Instead we interpret the instability as a consequence of the velocity profiles set up by the conditions. In other words, since is positive definite, this process experiences growth entirely due to the extraction of energy from the shear as embodied by the domain integral term . The character of the eigenfunctions are seen by inspecting the corresponding profiles for small values of in Fig. 1. In both the rigid and stressfree cases there are prominent boundary layers appearing on the starside for both the pressure and potential vorticity. The critical layer in the potential vorticity also appears here (see the inset in Figs. 1c,d) but is dwarfed by the starside boundary layer.
Figure 1: A comparison of eigenfunctions for a variety of starside boundary conditions. In all plots , and . The pressure eigenfunctions are shown on the left panel of plots while the potential vorticity eigenfunctions are shown on the right panel. The predicted growth rates are also quoted. a) Zero potential vorticity; b) zero potential vorticity gradient; c) rigid boundary; d) stressfree; and e) zeropressure fluctuation. 

Open with DEXTER 
Figure 2: Like Fig. 1 except . 

Open with DEXTER 
4.3.2 Zero PVanomaly
This boundary condition may be envisioned as the starside boundary counteracting any tendency for the development of any PVanomaly there. Although this is somewhat artificial, we present here the results of this investigation because these boundary conditions give solutions that closely resemble those obtained for the calculation on the semiinfinite domain calculation. As in the semiinfinite domain calculation, instability occurs when , and it scales as for . Inspection of the eigenfunctions in Fig. 1 for , especially the profiles for , shows that (i) the perturbed potential vorticity is strongly localized in the critical layers occurring where the real wavespeed approximately equals the background flow speed; (ii) an additional boundary layer appears at the starside boundary scaling like and; (iii) the vorticity in the critical layer follows the scaling determined in the semiinfinite domain calculation. We depict in Fig. 4a a contour plot of growth rates as a function of both the Reynolds number and the inverse of the Froude number, for fixed values of and k. The vertical axis may be understood as representing a positive increase in the wave's speed (see Eq. (46)). For the parameters depicted in Fig. 4a ( ), instability sets in for .
4.3.3 Zero radial PVanomaly gradient
Requiring no radial gradient of the PVanomaly on the starside boundary is arguably the least physically realistic but we include it here, as in the previous section, because it best reproduces the asymptotic result of the semiinfinite domain calculation. Like in the previous section, where is set to zero there, there is instability when , a critical layer emerges which also scales as for small values of . However, the boundary layer appearing near the starside boundary for the calculation of Sect. 4.3.2 vanishes here. Finally, for small values of (i.e. Re^{1}) the Reynolds Orr expression for these disturbances takes on the same leading form as (56) in Sect. 4.2; and this includes the character of . We note that growth rates here are nearly identical to the growth rates determined in Sects. 4.3.2. Figure 4b shows the landscape of instability (for ) and we see that instability also sets in when (i.e. here for ) but that it is bounded above by a more complicated function of Re. We note that the critical layer becomes harder to distinguish as viscosity (that is, Re^{1} or ) is made larger.
Figure 3: Contours on the Re plane of Im(c) for q=3/2, k=1 and . Given the value of k, according to its definition is written in terms of : a) Stressfree boundary conditions at x=0; b) Rigid boundary conditions at x=0. Shaded regions indicate growing modes. Note that according to (46) the vertical axis indicates increasing wavespeed. 

Open with DEXTER 
Figure 4: Same as Fig. 3 except: a) Zero PVanomaly at x=0; b) Zero PVanomaly gradient at x=0. 

Open with DEXTER 
5 Discussion and reflections
5.1 On the Tollmien Schlichting wave analogy
TS waves appear in flows that are are wallbounded at least on one boundary. In the classic analysis done for Blasius boundary layers (e.g. Schlichting & Gersten 2001) the instability is a solution of the 2D OrrSommerfeld equation. The mechanics leading to instability is understood to arise from the action of a purely viscous mode interacting with a nearly inviscid mode (Baines et al. 1996). The global velocity profiles set up are such that the relative phase between the horizontal and (plate) normalvelocities promotes extraction of energy from the shear which then leads to energetic growth (i.e. ). The 3D instability studied here shares some major similarities to classic TS waves: (i) the equation governing the dynamics of the potential vorticity modes (30) has similar structure to the 2D OrrSommerfeld equation; (ii) the instability emerges for both noslip and freeslip boundary conditions (but not limited to these). What stands out in our minds is that, although classical unstable TS waves come about in a wide variety of background flows excluding planeCouette flow (Baines et al. 1996; Schmid & Henningson 2000), the instability here is present for a planeCouette type of flow profile.
The results here compare qualitatively to the results of rotating Blasius boundary layers. For example, Yecko & Rossi (2004) show that three dimensional modal instability preferentially emerges in anticylonic rotating Blasius flow when the vertical wavenumber of the disturbance is large in comparison to its streamwise (azimuthal) wavenumber (e.g. see Fig. 8b in Yecko & Rossi 2004). By comparison, the asymptotic scalings we have implemented spotlights dynamics characterized by these same spatial scale disparities. Thus despite the differences in the problems investigated between these two studies (i.e. the inclusion of gravity and entropy gradients and the differing base velocity profiles) the similarities in the circumstances for instability onset suggest that such processes may be more general in environments like disks  especially near the diskstar boundary.
5.2 On the assumed constancy of g and N
To make the analysis we have exposed here analytically possible we assumed that the vertical component of gravity and the BruntVäisälä to be constant with respect to the disk vertical coordinate, (31) which permits us to assume separable normalmode solutions. On the other hand, real disks (including their small sections) are characterized by vertically varying values of g and N and this means that, in general, one cannot assume separable solutions in z and x, in particular^{}. We have checked that the results obtained in the limit where the viscosity parameter is small (i.e. ) still holds when g and N are taken to be correctly z dependent. Accordingly we have repeated the asymptotic calculation described in Sect. 4.2 where, in addition, we restricted our attention to finite vertical domains by imposing either velocity or pressure conditions on the vertical boundaries. Such disturbances will be characterized by vertical overtones labeled by an overtone wavelength  which should be thought of as being analagous to the vertical wavenumber we assumed in Sect. 4. The asymptotic calculation shows that instability sets in so long as a vertically weighted Froude number, , approximately exceeds one. The calculation is far more lengthly and does not add any new qualitative details to the one presented in this work and it is for this reason we have omitted it from the current exposition and we will expand upon it in a future study.
5.3 Relationship to the viscous pulsational instability
In an axisymmetric study of an shearing sheet section of an accretion disk of constant temperature, Latter & Ogilvie (2006) argue that the VPI (Kato 1978) is most likely to manifest itself through the destabilization of an even structured fmode. The disturbances become unstable because the viscous perturbations transfer energy from the shear into the acoustic mode through the perturbed stress T_{xy}'. Because fmodes are characteristically inertialacoustic waves, they are the likely candidates for this instability since their vertical structures are the simplest which, in turn, result in minimizing dissipative losses. We observe that T_{xy}' is proportional to the pressure fluctuation. This fluctuating stress affects the evolution of the horizontal velocity perturbation by extracting/adding energy into the disturbance. As the horizontal velocity disturbances are not in general in phase with the pressure fluctuations, especially for inertial/acoustic modes, the possibility for overstability is manifest. However, the PVanomalies examined here are distinct from inertialgravity and inertialacoustic modes (e.g. Ogilvie 1998; Tevzadze et al. 2004). The horizontal velocity perturbations of PV anomalies are proportional to the radial gradient of the perturbation pressure. Since the latter quantity is proportional to a decaying exponential (i.e. ), it follows that the phase between T_{xy}'and the energy in the PVanomaly will be radians out of phase with each other (see Sect. 4.1). Thus the fluctuating stress T_{xy}' behaves to stabilize a PVanomaly. resulting in a stabilizing relationship as we observe in Sect. 4.1. As such we understand the destabilization of the PVanomaly as as being distinct from the instability leading to the VPI.
5.4 Summary and Implications
For cold disk systems, i.e. those in which magnetic effects are not active, the prospects of identifying instability mechanisms appear to be far from exhausted. We have tried to argue that certain previously considered nonmagnetic instabilities need not be discarded as candidate mechanisms driving activity for disks. Indeed the SRI and PPI instabilities, which emerge as the interaction of edgewaves along cylinder/channel walls, could in principle operate in real disks so long as there exists, in general, interacting waves propagating separated from each other by a waveevanescent region. We have demonstrated here another possible mechanism  that the existence of disk analogs of unstable TollmienSchlichting waves could also manifest themselves in real disk systems. We have carried out the calculation within a model shearing box in which we have imposed a single boundary on one side. True disks have boundary layers separating stars from the disk which are probably far more complicated (Regev & Bertout 1995) than the model we have presented here. Nonetheless, far from being a proof, we have demonstrated in this asymptotically simplified model that such a dynamical processes is, at least, feasible. It is no stretch of the imagination to suppose that analogous unstable waves may exist near the vicinity of the stardisk boundary layer. We add a final reflection. Classical unstable TS waves emerge in flows with compliant boundaries showing that such instabilities are robust and persist even if the boundaries have a certain amount of elasticity to them (Carpenter & Garrad 1985)  although compliant walls delay the onset of instability to higher Reynolds numbers. As a stardisk boundary is probably not a rigid body transition, it would be beneficial to investigate and/or model these disk analog TS waves by considering starside boundary conditions that are appropriately compliant as well.The perspective we have adopted therefore can be broken down into two parts. The first is that (as in the recent studies of Kleiber & Glatzel 1999 and Latter & Ogilvie 2006) we have expanded the exploration of the possible destabilizing role viscosity can play. Viscosity does not always have to stabilize disturbances as there are velocity profiles, dictated by boundary conditions, wherein destabilization occurs counter to one's usual physical intuition. The second is that experimentation with boundary conditions, even within the context of the shearing box environment, followed by concerted effort toward understanding and clarifying their effects is a worthwhile endeavor given our lack of complete knowledge about the boundaries of real disk systems (a situation which is strongly contrasted by what is encountered in laboratory/terrestrial flows). If turbulent stresses in cold disks are driven by the MRI resulting in effective values of (Ogilvie 2001; Ogilvie & Proctor 2003; King et al. 2007), then TS waves like the sort here could emerge as a secondary instability. This possibility is made manifest because the TS wave instability grows faster in proportion to the value of until about a value of 0.1.
Acknowledgements
The authors are indebted to the valuable comments and suggestions of the anonymous referee. The authors would like to thank the Israeli Science Foundation for making this research possible. O.M.U. also acknowledges that this research was partly supported by BSF grant 2004087 and ISF grant 1084/06. O.M.U. also thanks the Dead Sea Regional Council and the Ein Gedi Kibbutz for their hospitality and Phil Yecko for suggestive conversations.
Appendix A: Scaling arguments leading to the viscous large shearing box equations
The derivation of the VLSB equations follows the procedure executed in Umurhan & Regev (2004). The dimensional equations of motion in cylindrical coordinates in a frame of constant rotation are, in component form, given by the followingthe equations of mass continuity and entropy
=  0,  (A.4)  
=  Q,  (A.5) 
with the operator definition
are the radial (r), azimuthal ) and vertical (z) velocities as observed in the rotating frame. The rotation rate is set to the rotation at some fiducial radius r=R_{0}. The curious notation on is meant to indicate that the azimuthal velocity observed in the laboratory frame (denoted by ) would be related to its velocity in the rotating frame by . The entropy is defined to be in which C_{V} is the specific heat at constant volume, is the ratio of specific heats, that is in which C_{P} is the specific heat at constant pressure. Q is a heat function representing some sort of nonadiabatic processes that are relevant to the disk.
The
gravitational potential
is written in the following unusual
form in order to effect some generality when it comes to the shear
it induces upon the steady state flow,
For the realistic Keplerian flow profile q = 3/2 and , where M_{*} is the mass of the central object. However, is in general a function of q and only when q=3/2 it is to be interpreted as the ``Keplerian'' rotation rate.
The viscous moments are
=  (A.7)  
=  (A.8)  
=  (A.9) 
along with the viscous stresses t_{ij},
(A.10)  
(A.11)  
(A.12)  
(A.13)  
(A.14)  
(A.15) 
in which
Because the turbulent viscosity within a disk is presumed to be driven by some sort of shear process (either the MRI or a subcritical hydrodynamic transition) the bulk viscosity will be set to zero hereinafter. According to the classic proposal of Shakura & Sunyaev (1973), the shear viscosity is parametrized as
(A.16) 
where is the steadyrotation rate induced by the generalized potential (A.6). When we consider Keplerian flows we write (and see below). The parameter is a tunable order 1 quantity.
We proceed from here onto nondimensionalization. We shall consider dynamics as taking place in a small region around the point , which we shall refer to as the box. We let measure the nondimensional size of this box, i.e.,
This motivates us to scale all the directions by R_{0} and to
define the nondimensional coordinates
(A.17) 
where it is understood that x,y,z are now taken to be order 1 quantities.
Furthermore we suppose that all density quantities are scaled by
the reference density ,
pressure quantities are supposed
similarly scaled by
where
is the dimensional scale of the sound speed of the box. In
parallel with this speed is the local rotation speed of the box
around the central object,
:
when q=3/2 this speed is sometimes referred to as
the local Keplerian speed of the disk. The fundamental
ansatz of cold thin disk theory is that the ratio of these two
quantities is small. We, in fact, identify
(A.18) 
which is the classic parameter measuring the thinness or ``coldness'' of the disk (Shakura & Sunyaev 1973; and see recently Umurhan et al. 2006, for recent exploitations of this parameter). This, in turn, is a measure of the disk's vertical scaleheight . We also note here that by equating the size of the box to the ``coldness'' of the disk means we are here looking at the viscous analog of the socalled ``Large Shearing Box'' equations formally developed in Umurhan & Regev (2004).
Thus we propose that all velocities observed in this rotating
frame are scaled by the soundspeed, meaning to say that
where u',w',v' are order 1 nondimensionalizations of the corresponding velocities. Note that in these scalings we are saying that all velocities observed in the rotating frame are order smaller than the basic rotation speed of the box V_{0}. This seemingly obvious point is emphasized because the magnitude scale of the accretion and meridional flow induced by the turbulent viscosity is order smaller than V_{0} (Shakura & Sunyaev 1973; Kluzniak & Kita 1999). It means that to the order to which matters are considered here, that is as far as the generalized ``box'' formalism is concerned (see below), the effects of accretion and meridional flow are absent.
Time and all advective derivatives are scaled
according to the local rotation time of the disk, i.e.
.
Before putting these scalings into the
governing equations we note that
Taking into account all of the nondimensionalizations, along with writing and understanding that are the nondimensionalized density and pressure quantities, we find that the equations of motion (A.1A.3) become
are as they appear in the text. To reiterate, the effects grouped in the terms consist of curvature effects, higher order corrections due to the central potential, and turbulent viscosity induced accretion/meridional flow. The equation of continuity and entropy conservation also become
To obtain the viscous large shearing box equations presented in the text, i.e. (15), we make the following identifications and assumptions
 Drop all terms and higher from (A.19A.23).
 Everywhere write in (A.19A.23) in order to eliminate the sole x term on the RHS of (A.19). The expression qx is the background shear and will be felt by perturbations.
 Write the density and pressures as being comprised of
a steady portion and a time dependent perturbed portion, i.e.
 Specifically flag all Coriolislike related effects with the symbol .
Appendix B: Semiinfinite domain calculation details
B.1 Outer solution completion
As a result of this breakdown we write the solutions to this order according to which
side of x_{c} one is on. Formally then we say
(B.1) 
for , and
(B.2) 
for . The particular solutions satisfy
(B.3) 
with solutions given by
(B.4)  
(B.5) 
where is the Exponential Integral function (Abramowitz & Stegun 1972). We define
An analysis of these solution forms for shows that
In formally writing these solutions we purposely avoid the region described by
which we will refer to as the critical layer. The size of this bounding region is such that . It will be shown that the solutions on either side of this region will asymptotically match at lowest nontrivial order to the solution emerging from the critical layer itself. Finally, the u=0 boundary condition at x=0 may be read from (36) after making use of the above relationships
B.2 Critical layer calculation and matching
The calculation will be facilitated if we consider the evolution of the residual potential vorticity quantity throughThe polynomial terms in the above expression are the first five terms of the series expansion of the lowest order outer solution, expressed in the vicinity of the critical layer. We may now rewrite (51) instead in terms of ,
(B.7) 
where is as defined in the text. At we have simply
(B.8) 
Homogeneous solutions of the operator of these involve integrals of Airy Functions associated. However, because these are functions of complex arguments these two solutions are rejected because they show exponential growth as (e.g. Bender & Orszag 1999). All homogeneous solutions involving these Airy functions are rejected henceforth. Thus we assume the following solution expansion for
The order solution will be . There is also a solution at this order proportional to but it is rejected on account of the fact that it will not match any corresponding solution from outside. The equation at order is
With . With solution
(B.10) 
In which and where
Operating on by followed by an integration by parts verifies that is the particular solution of (B.9). We note immediately that is bounded for all real values of x by since
It is instructive to note the expression
where will be related to the expressions of associated with the outer solutions (see next section). This relationship plays the role connecting the two homogeneous outer solutions. When an asymptotic evaluation of the integral shows that
We note that there exists a phase factor proportional to that arises from this operation. This shall be explicitly referred to in the next section. One final note: in the process of matching to the outer solutions the expression produces terms proportional to in the outer region. In order to have this appropriately matched (in this case, cancelled) is the reason why there is an term in the critical layer expansion (52).
We note that the PVanomaly arising from this
critical layer appears at order
when viewed in the
unstretched coordinate frame. This is because the first nontrivial
contribution arising to the potential vorticity perturbation from this zone is
=  
= 
Thus, while the size of the critical layer zone is order the magnitude of the PVanomaly will be order the leading order pressure perturbation.
B.3 Matching
We must reexpress the outer solution in a ``small'' vicinity of the critical point x_{c}. We consider first the solutions approaching from below, that is . We have
(B.14) 
while when approaching this point from above, that is as
,
we
find
(B.15) 
Now we do the same to the critical layer solution. We restore the inner
coordinate in terms
of xx_{c} and take the limit of small
(as is standard practice in boundary
layer theory, Bender & Orszag 1999) revealing
=  
(B.16) 
Identifications are made respecting powers of and xx_{c}. The matchings are straightforward since the solutions coming in from the left and from the right of x_{c} are the same:
while the remaining two terms, proportional to (xx_{c})^{3} and (xx_{c})^{4}are satisfied given the above assignments in (B.17). At we find first that
because the above relationship implies
To complete the matching we must prepare the final term . It is asymptotically correct to do this by considering the derivative of these terms up to the bounding region of of the critical layer, i.e. . Thus we consider the derivatives as one approaches this zone (and measured by the coordinates ) as
(B.19) 
that is, approaching from below and
(B.20) 
that is, approaching from above. Subtracting the two expressions gives
Now we must match this to the corresponding expression emerging from the interior of the domain. In other words we require similarly that
in which (i) the transition in orders of occurs because of the change of variables from x to ; (ii) we have identified ; and (iii) used (B.12) in writing the above expression. We noted that the leading behavior of the integral for (i.e. for the matching zone) is
Equating the RHS of (B.21) and (B.22) and making use of the asymptotic form (B.23) we see that the offending logarithmic terms cancel leaving,
The meaning of this relationship is clear  the presence of the singular layer in the flow means that the homogeneous outer region solutions must show a jump in their derivatives (in proportion to the RHS of the above expression). Another way to interpret this is to recognize that this jump corresponds to the presence of a vortex sheet at x = x_{c}.
The complex wavespeed may be now obtained from (B.18), (B.24)
and the boundary condition (B.6). We note that for there
to be a nontrivial solution to (B.18) and (B.24) the
following
relationship,
must be satisfied. This is the solvability condition. The second matter we note is that
Thus, the solvability condition (B.25) together with the above
expressions used in (B.6) combine to give
The form quoted in the text follows from restoring the definition of into the above expression.
References
 Afshordi, N., Mukhopadhyay, B., & Narayan, R. 2005, ApJ, 629, 373 [NASA ADS] [CrossRef] (In the text)
 Abramowitz, M., & Stegun, I. A. 1972, Handbook of Mathematical Functions (Dover)
 Baines, P. G., Majumdar, S., & Mitsudera, H. 1996, J. Fluid Mech., 312, 107 [NASA ADS] [CrossRef] (In the text)
 Baines, P. G., & Mitsudera, H. 1994, J. Fluid Mech., 276, 327 [NASA ADS] [CrossRef] (In the text)
 Balmforth, N. J. 1999, J. Fluid Mech., 387, 97 [NASA ADS] [CrossRef] (In the text)
 Balmforth, N. J., & Piccolo, C. 2001, J. Fluid Mech., 449, 85 [NASA ADS] [CrossRef] (In the text)
 Balmforth, N. J., Piccolo, C., & Umurhan, O. M. 2001, J. Fluid Mech., 449, 115 [NASA ADS] [CrossRef]
 Balbus, S. A. 2003, ARA&A, 41, 555 [NASA ADS] [CrossRef] (In the text)
 Barranco, J. A., & Marcus, P. S. 2005, ApJ, 623, 1157 [NASA ADS] [CrossRef] (In the text)
 Bayly, B. J., Orszag, S. A., & Herbert, T. 1988, Ann. Rev. Fluid Mech., 20, 359 [NASA ADS] [CrossRef]
 Bender, C. M., & Orszag, S. A. 1999, Advanced Mathematical Methods for Scientists and Engineers (Springer)
 Bodo, G., Tevzadze, A., Chagelishvili, G., et al. 2007a, A&A, 475, 51 [NASA ADS] [CrossRef] [EDP Sciences] (Bodo et al. A) (In the text)
 Bodo, G., Chagelishvili, G., Murante, G., et al. 2007b [arXiv:0705.3474v1]
 Carpenter, P. W., & Garrad, A. D. 1985, J. Fluid Mech., 155, 465 [NASA ADS] [CrossRef] (In the text)
 Case, K. M. 1960, Phys. of Fluids, 3, 143 [NASA ADS] [CrossRef] (In the text)
 Chagelishvili, G. D., Zahn, J.P., Tevzadze, A. G., & Lominadze, J. G. 2003, A&A, 402, 401 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Charney, J. G., & Stern, M. E. 1962, J. Atmos. Sci., 19, 159 [NASA ADS] [CrossRef]
 Davies, H. C., & Bishop, C. H. 1994 J. Atm. Sci., 51, 1930
 Doering, C. R., Spiegel, E. A., & Worthing, R. A. 2000, Phy. Fluids, 12, 1955 [NASA ADS] [CrossRef] (In the text)
 Drazin, P. G., & Reid, W. H. 1984, Hydrodynamic Stability, Cambridge
 Dubrulle, B., Marie, L., Normand, Ch., et al. 2004, A&A, 429, 1 [NASA ADS] [CrossRef]
 Frank, J., King, A. R., & Raine, D. J. 2002, Accretion Power in Astrophysics (Cambridge: Cambridge Univ. Press)
 Goldreich, P., & LyndenBell, D. 1965, MNRAS, 130, 125 [NASA ADS]
 Goldreich, P., Goodman, J., & Narayan, R. 1986, MNRAS, 221, 339 [NASA ADS]
 Hayashi, Y.Y., & Young, W. R. 1987, J. Fluid Mech., 184, 477 [NASA ADS] [CrossRef] (In the text)
 Hoskins, B. J., McIntyre, M. E., & Robertson, A. W. 1985, Quart. J. Roy. Meteor. Soc., 111, 877 [NASA ADS] [CrossRef] (In the text)
 Ioannaou, P. J., & Kakouris, A. 2001, ApJ, 550, 931 [NASA ADS] [CrossRef]
 Ji, H., Burin, M., Schartman, E., & Goodman, J. 2006, Nature, 444, 343 [NASA ADS] [CrossRef] (In the text)
 Johnson, B. M., & Gammie, C. F. 2005, ApJ, 635, 149 [NASA ADS] [CrossRef] (In the text)
 Joseph, D. D. 2003, J. Fluid. Mech., 479, 191 [NASA ADS] [CrossRef] (In the text)
 Kato, S. 1978, MNRAS, 185, 629 [NASA ADS] (In the text)
 Kersalé, E., Hughes, D. W., Ogilvie, G. I., Tobias, S. M., & Weiss, N. O. 2004, ApJ, 602, 892 [NASA ADS] [CrossRef] (In the text)
 Klahr, H. H., & Bodenheimer, P. 2003, ApJ, 582, 869 [NASA ADS] [CrossRef] (In the text)
 Kleiber, R., & Glatzel, W. 1999, MNRAS, 303, 107 [NASA ADS] [CrossRef] (In the text)
 King, A. R., Pringle, J. E., & Livio, M. 2007, MNRAS, 376, 1740 [NASA ADS] [CrossRef] (KPL07) (In the text)
 Kluzniak, W., & Kita, D. 2000, Threedimensional structure of an alpha accretion disk [arXiv:astroph/0006266]
 Latter, L. N., & Ogilvie, G. I. 2006, MNRAS, 372, 1829 [NASA ADS] [CrossRef] (In the text)
 Lesur, G., & Longaretti, P.Y. 2005, A&A, 444, 25 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Li, H., Finn, J. M., Lovelace, R. V. E., & Colgate, S. A. 2000, ApJ, 533, 1023 [NASA ADS] [CrossRef] (In the text)
 Lions, P. L. 1993, Limits incompressible et acoustique pour des fluides visqueux, compressible et isentropique, C. R. Acad. Sci. Paris Ser. I Math, 317, 1197 (In the text)
 Lithwick, Y. 2007, ApJ, 670, 789 [NASA ADS] [CrossRef] (In the text)
 Mukhopadhyay, B., Afshordi, N., & Narayan, R., 2004, ApJ, 629, 383 [NASA ADS] [CrossRef]
 Ogilvie, G. I. 1998, MNRAS, 297, 291 [NASA ADS] [CrossRef] (In the text)
 Ogilvie, G. I. 2001, MNRAS, 325, 231 [NASA ADS] [CrossRef] (In the text)
 Ogilvie, G. I., & Proctor, M. R. E. 2003, J. Fluid Mech., 476, 389 [NASA ADS] [CrossRef]
 Papaloizou, J. C. B., & Pringle, J. E. 1984, MNRAS, 208, 721 [NASA ADS]
 Petersen, M. R., Julien, K., & Stewart, G. R. 2007a, ApJ, 658, 1236 [NASA ADS] [CrossRef] (In the text)
 Petersen, M. R., Stewart, G. R., & Julien, K. 2007b, ApJ, 658, 1252 [NASA ADS] [CrossRef]
 Regev, O., & Bertout, C. 1995, MNRAS, 272, 7179 [NASA ADS] (In the text)
 Richard, D. 2001, Thèse de doctorat, Université Paris 7 (In the text)
 Richard, D., & Zahn, J.P. 1999, A&A, 347, 734 [NASA ADS] (In the text)
 Schlichting, H., & Gersten, K. 2001, Boundary Layer Theory, Springer, 94 (In the text)
 Sakai, S. 1989, J. Fluid Mech., 202, 149 [NASA ADS] [CrossRef] (In the text)
 Schmid, P. J., & Henningson, D. S. 2000, Stability and Transition in Shear Flows, Springer (In the text)
 Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 [NASA ADS]
 Shalybkov, D., & Rüdiger, G. 2005, A&A, 438, 411 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Spiegel, E. A., & Veronis, G. 1960, ApJ, 131, 442 [NASA ADS] [CrossRef] (In the text)
 Stewartson, K. 1981, IMA J. Appl. Math., 27, 133 [CrossRef] (In the text)
 Tevzadze, A. G., Chagelishvili, G. D., Zahn, J.P., Chanishvili, R. G., & Lominadze, J. G. 2003, A&A, 407, 779 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Umurhan, O. M. 2006, MNRAS, 365, 85 [NASA ADS] [CrossRef] (In the text)
 Umurhan, O. M. 2008, A&A, 489, 953 [NASA ADS] [CrossRef] [EDP Sciences]
 Umurhan, O. M., Nemirovsky, A., Regev, O., & Shaviv, G. 2006, A&A, 446, 1 [NASA ADS] [CrossRef] [EDP Sciences]
 Umurhan, O. M., & Regev, O. 2004, A&A, 427, 855 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Yecko, P. A. 2004, A&A, 425, 385 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Yecko, P. A., & Rossi, M. 2004, Phys. Fluids, 16, 2322 [NASA ADS] [CrossRef]
Footnotes
 ... other^{}
 By this it is meant to say that there exists a wave evanescent region separating the waves in question.
 ... (1985)^{}
 They show that the interior of the domain does not intrinsically support propagating waves as it behaves like a wave evanescent zone. However the imposition of the boundary conditions brings into existence waves that propagate along the radial boundaries of the domain.
 ... region^{}
 Examples of this process are well known in atmospheric flow (Charney & Stern 1965; Hoskins et al. 1985; Davies & Bishop 1994).
 ...^{}
 A momentary comparison to (18) should convince the reader that represents a perturbed entropy quantity (Dubrulle et al. 2005) making (25) a reasonable analog of the linearized form of (5).
 ... Number^{}
 However, note that the actual Reynolds number, denoted by R, is a factor of larger than Re, i.e. R = Re.
 ...^{}
 With respect to the other scalings we have called on in this work, we shall formally require that .
 ... condition^{}
 Note, however, that we show in Sect. 4.3 that the results obtained in this section, including the calculation of the growth rates, would be essentially unchanged had we imposed, instead, that the radial gradient of the potential vorticity be zero at the starside boundary.
 ...^{}
 Note, however, this limit breaksdown the QHSG approximation and is not considered.
 ... particular^{}
 In other words by assuming g and N to be constant we are able to assume solution form while if g and N are zdependent one can (at best) assume a solution in the form where is the nonseparable structure function.
All Figures
Figure 1: A comparison of eigenfunctions for a variety of starside boundary conditions. In all plots , and . The pressure eigenfunctions are shown on the left panel of plots while the potential vorticity eigenfunctions are shown on the right panel. The predicted growth rates are also quoted. a) Zero potential vorticity; b) zero potential vorticity gradient; c) rigid boundary; d) stressfree; and e) zeropressure fluctuation. 

Open with DEXTER  
In the text 
Figure 2: Like Fig. 1 except . 

Open with DEXTER  
In the text 
Figure 3: Contours on the Re plane of Im(c) for q=3/2, k=1 and . Given the value of k, according to its definition is written in terms of : a) Stressfree boundary conditions at x=0; b) Rigid boundary conditions at x=0. Shaded regions indicate growing modes. Note that according to (46) the vertical axis indicates increasing wavespeed. 

Open with DEXTER  
In the text 
Figure 4: Same as Fig. 3 except: a) Zero PVanomaly at x=0; b) Zero PVanomaly gradient at x=0. 

Open with DEXTER  
In the text 
Copyright ESO 2009