Supersonic Shear Instabilities in Astrophysical Boundary Layers
SSupersonic Shear Instabilities in Astrophysical Boundary Layers
Mikhail A. Belyaev & Roman R. Rafikov , ABSTRACT
Disk accretion onto weakly magnetized astrophysical objects often proceeds via aboundary layer that forms near the object’s surface, in which the rotation speed of theaccreted gas changes rapidly. Here we study the initial stages of formation for such aboundary layer around a white dwarf or a young star by examining the hydrodynamicalshear instabilities that may initiate mixing and momentum transport between the twofluids of different densities moving supersonically with respect to each other. We findthat an initially laminar boundary layer is unstable to two different kinds of instabilities.One is an instability of a supersonic vortex sheet (implying a discontinuous initial profileof the angular speed of the gas) in the presence of gravity, which we find to have a growthrate of order (but less than) the orbital frequency. The other is a sonic instabilityof a finite width, supersonic shear layer, which is similar to the Papaloizou-Pringleinstability. It has a growth rate proportional to the shear inside the transition layer,which is of order the orbital frequency times the ratio of stellar radius to the boundarylayer thickness. For a boundary layer that is thin compared to the radius of the star,the shear rate is much larger than the orbital frequency. Thus, we conclude that sonicinstabilities play a dominant role in the initial stages of nonmagnetic boundary layerformation and give rise to very fast mixing between disk gas and stellar fluid in thesupersonic regime.
Subject headings: accretion, accretion disks – hydrodynamics — waves – instabilities
1. Introduction
Accretion onto astrophysical objects possessing a material surface (as opposed to accretiononto black holes) always involves a non-trivial interaction of the incoming gas with the object’souter layers. Examples of such objects include white dwarfs in cataclysmic variables (CVs), youngstars 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 objectand how it mixes with the previously accreted material. Department of Astrophysical Sciences, Princeton University, Ivy Lane, Princeton, NJ 08540;[email protected] Sloan Fellow a r X i v : . [ a s t r o - ph . S R ] D ec B ∗ , crit = 5 × (cid:18) β. (cid:19) − / (cid:18) M ∗ . M (cid:12) (cid:19) / (cid:32) ˙ M − M (cid:12) yr − (cid:33) / (cid:18) R ∗ × km (cid:19) − / G . (1)Here, we have used parameters for the mass and radius that are typical for CVs in outburst(Bergeron et al. 1992; Hachisu 1986; Patterson 1984). The parameter β is a dimensionless factorof order unity that depends on the model adopted for the disruption of the disk by the stellarmagnetic field (Ghosh & Lamb 1978).For B ∗ < B ∗ , crit , the accretion disk will extend all the way to the surface of the star, for whichthere is good observational evidence in neutron star low mass X-ray binaries (Gilfanov et al. 2003)and in dwarf nova systems in outburst (Wheatley et al. 2003). Since the rotation rate of the star,Ω ∗ , is slower than the Keplerian rotation rate at the stellar surface, Ω K ( R ∗ ), a boundary layer(BL) will exist inside of which d Ω /dR > An interesting aspect of the BL problem that has received less attention in the past is the issueof BL initiation , i.e. the initial stage of BL formation which must occur when the accreted materialfirst touches the surface of the star. In this case the physical setup is going to be quite differentfrom the steady-state BL primarily because of the much larger velocity shear, justifying the studyof BL initiation in the limit of an essentially discontinuous rotation profile at the interface betweenthe star and accreting material. 3 –Even though the initiation stage is just a transient phase in the BL evolution, one still needsto understand it to get a full picture of the BL phenomenon. Also, one need not think that the BLinitiation is a unique event for every accreting object — it may, in fact, be repetitive. The mostcommon situation in which this BL initiation is recurrent involves the dramatic increase of the massaccretion rate (due to some sort of outburst triggered by an instability) through the disk which ismagnetically disrupted outside the star under normal circumstances. According to Equation (1),a sudden increase of ˙ 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 pointat 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, Livio& Pringle (1992) have argued that the observed lag in the rise time of the UV emission relative tothe optical in a CV system transitioning to outburst can be explained if the disk is magneticallydisrupted in quiescence but not in outburst. Thus, the hottest, innermost part of the disk isevacuated 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 stellarsurface, due to the high accretion rate - ˙ M ∼ − M (cid:12) yr − (Kenyon et al. 1988). In these systems,the BL can puff up to become of order the stellar radius and the distinction between the boundarylayer, the disk, and the star becomes blurred (Popham et al. 1993).The goal of our present work is to make a first step towards understanding the formation ofthe BL. As in the case of a well-developed BL, the major issue for this initial stage lies in details ofthe angular momentum transport and mixing. However, because of the enormous shear present atthe star-disk interface at this stage, it is highly likely that purely hydrodynamical shear instabilitieswould 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 identifya particular variety of the shear instability that most efficiently initiates mixing between the twofluids, 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 insta-bilities: stratification, rotation, magnetic fields, turbulence, radiative transfer, and the supersonicnature of the flow. The latter aspect of the problem is very important and is inevitable duringthe BL initiation phase, when disk material rotating at the Keplerian speed comes into contactwith the more slowly spinning stellar surface. The differential azimuthal velocity between the twointeracting flows is then bound to be a significant fraction of the Keplerian velocity at R ∗ andshould highly exceed the sound speed both in the disk and in the outer layers of the star. 4 – The study of shear instabilities in compressible fluids is of fundamental physical significance andhas received much attention in the past. Landau (1944), Hatanaka (1949), Pai (1954), Miles (1958),and Gerwin (1968) have all studied the problem in the vortex sheet approximation, where one halfplane of compressible fluid moves at constant velocity over another. When the two fluids move ata low Mach number relative to each other, one finds that infinitesimal perturbations are governedby the classical Kelvin-Helmholtz (KH) dispersion relation for two incompressible fluids. However,Miles (1958) showed that above a critical Mach number, the vortex sheet becomes marginallystable to two-dimensional disturbances along the direction of the flow. This is a surprising resultconsidering that one might have expected increasing the shear to lead to an increased growth ratefor the instability rather than stabilization. Understanding the implications of this result for themomentum and mass transport in the astrophysical BLs is one of the goals of our study.Later, Blumen et al. (1975), Ray (1982), Choudhury & Lovelace (1984) and Glatzel (1988)studied the stability of a fluid with a continuous and monotonically varying velocity profile. Theyfound that unlike the vortex sheet, continuous velocity profiles with a supersonic velocity differenceacross them were unstable even at high Mach number. Glatzel (1988) showed that the instabilitywas similar to the Papaloizou-Pringle instability which operates in hydrodynamical disks with radialboundaries (Papaloizou & Pringle 1984). 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 ∝ /L , where L is the thickness of the shear layer. Thus, the thinner the shear layer, the faster theinstability operates! This is exactly the opposite of what one might have expected from the vortexsheet stability criterion and we will provide an explanation for this in § § § §
4, we introduce gravity as a small parameter and perform a perturbation expansion to studywhat influence this has on the stability of the vortex sheet. Finally, in §
2. Formalism
Here we discuss the formalism that we use for analytic calculations throughout the paper. Weadopt cylindrical coordinates ( (cid:36), θ, z ) and for simplicity assume that all equilibrium quantities areindependent of z . Thus, for our purposes gravity is given by − g ( (cid:36) ) ˆ (cid:36) ( g ( (cid:36) ) > d Ω /dz (cid:54) = 0, but it does contain the necessary ingredients for studying shearinstabilities. Given our assumptions, the equation of hydrostatic equilibrium reads¯ g ( (cid:36) ) ≡ g ( (cid:36) ) − (cid:36) Ω ( (cid:36) ) = − ρ dPd(cid:36) , (2)where ¯ g ( (cid:36) ) is the effective gravity. The equilibrium density, pressure, and sound speed are givenby ρ ( (cid:36) ), P ( (cid:36) ), and s ( (cid:36) ) respectively.Denoting the ( (cid:36), θ, z ) velocities by ( u, v, w ), using v = Ω (cid:36) , and assuming adiabaticity, theEuler equations in cylindrical coordinates are − ∂ρ∂t = 1 (cid:36) ∂∂(cid:36) ( ρ(cid:36)u ) + 1 (cid:36) ∂∂θ ( ρv ) + ∂∂z ( ρw ) (3) DuDt = − ρ ∂P∂(cid:36) + v (cid:36) − g (4) DvDt = − ρ(cid:36) ∂P∂θ − uv(cid:36) (5) DwDt = − ρ ∂P∂z (6) D ( P ρ − γ ) Dt = 0 (7) DDt ≡ ∂∂t + u ∂∂(cid:36) + v(cid:36) ∂∂θ + w ∂∂z (8)On top of the zeroth order equilibrium state, we consider infinitesimal two-dimensional pertur-bations of the form δf ( (cid:36) ) exp[ i ( mθ + k z z − ωt )], where the δ denotes an Eulerian perturbation. 6 –Starting from Equations (3)-(7), the linearized first order equations are given by i ¯ ωδρ = 1 (cid:36) ( (cid:36)ρδu ) (cid:48) + imρ(cid:36) δv + ik z ρδw (9) i ¯ ωδu = 1 ρ δP (cid:48) + ¯ gρ δρ − δv (10) i ¯ ωδv = im(cid:36)ρ δP + 2 Bδu (11) i ¯ ωδw = ik z ρ δP (12) i ¯ ωδρ = − C L ρδu + i ¯ ωs δP, (13)where ¯ ω = ω − m Ω is the phase speed in the locally corotating frame, B = (Ω + d ( (cid:36) Ω) /d(cid:36) ) / C L = γ − d ln P/d(cid:36) − d ln ρ/d(cid:36) is the Ledoux discriminant, and the primesdenote differentiation with respect to (cid:36) . Note that in Equations (9)-(13), we have made Cowling’sapproximation by ignoring the first order perturbation to the gravitational potential (Cowling1941). 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 δP and δu : i (cid:18) ¯ ω s − k z − m (cid:36) (cid:19) δP = ρ (cid:18) C L ¯ ω + 2 Bm(cid:36) (cid:19) δu + ¯ ω(cid:36) ( (cid:36)ρδu ) (cid:48) (14) iρ (¯ ω − ¯ gC L − κ ) δu = ¯ ωδP (cid:48) + ¯ ω ¯ gs δP − m(cid:36) δP. (15)Here we have introduced the epicyclic frequency κ = 4 B Ω.One generally expects the BL width δ BL to be small compared to the stellar radius, δ BL (cid:28) R ∗ .Under this assumption, we show in Appendix A that Equations (14)-(15) can be simplified to thefollowing set of equations: i (cid:18) ¯ ω s − k z − m R ∗ (cid:19) δP = ρ (cid:18) C L ¯ ω + mR ∗ S (cid:19) δu + ¯ ω ( ρδu ) (cid:48) (16) iρ (¯ ω − ¯ gC L ) δu = ¯ ωδP (cid:48) + ¯ ω ¯ gs δP. (17)The form of Equations (16) and (17) is the same as for a plane-parallel stratified shear flow (Alexakiset al. 2002). This suggests that the effects of Coriolis force and curvature are unimportant in aradially thin BL if shear instabilities provide the turbulence, and leads us to redefine the cylindrical( (cid:36), θ, z ) coordinate system into a Cartesian ( x, y, z ) system. The perturbed quantities now havethe form δf ( x ) exp[ i ( k y y + k z z − ωt )], where k y ≡ m/R ∗ . From this it immediately follows that¯ ω = ω − k y V y ( x ), where V y ( x ) = R ∗ Ω( x ) is the velocity profile. From here on, we will ignore therotation terms, and treat the flow in the BL as a plane parallel shear flow.From Squire’s theorem for plane parallel shear flows (Fejer & Miles 1963), any three-dimensionalperturbation is mathematically equivalent to a two-dimensional perturbation upon making the 7 –transformations k y → k y / cos θ , k z →
0, and V y → V y cos θ , where cos θ = k y / (cid:113) k y + k z . Thus, itis sufficient to only consider two-dimensional perturbations ( k z = 0) in the stability problem, andfor the rest of this paper we take perturbations in the form δf ( x ) exp[ i ( k y y − ωt )] . (18)Assuming two-dimensional perturbations, Equations (16) and (17) become i (cid:18) ¯ ω s − k y (cid:19) δP = ρ (cid:0) C L ¯ ω + k y V (cid:48) y (cid:1) δu + ¯ ω ( ρδu ) (cid:48) (19) iρ (¯ ω − ¯ gC L ) δu = ¯ ωδP (cid:48) + ¯ ω ¯ gs δP, (20)where the primes now denote differentiation with respect to x . Equations (20) and (19) can beused to obtain the generalized Rayleigh equation (Appendix B.1), which using a notation similarto Alexakis et al. (2002) reads δφ (cid:48)(cid:48) + (cid:32) k x − g ˜ ρ k s + k g ˜ W − ˜ W (cid:48)(cid:48) ˜ W (cid:33) δφ = 0 . (21)Here, δφ = − δu √ ρ/k x , (22)is a modified stream function, and for simplicity we have dropped the bar over ¯ g , so now g denotesthe effective gravity. We have also defined the quantities:˜ W = k y W (cid:112) ˜ ρ/ik x (23) W = V y − c where c = ω/k is the phase speed. (24) k x = k y ( W /s −
1) is the square of the x-component of the wavevector (25)in the absence of shear or stratification. (26)˜ ρ = ρf , where ρ is the density. (27) f = exp (cid:18)(cid:90) x k g ( ζ ) dζ (cid:19) . (28) k g = g/s is a measure of the inverse of the local scale height. (29) k s = ρ (cid:48) /ρ is the inverse stratification length scale. (30)Note that the Ledoux discriminant is given by C L = − ( k g + k s ).The boundary conditions on δφ that must be satisfied at a discontinuity in the density or thevelocity are that the upper and lower fluids stay in contact and that the pressure perturbation iscontinuous across the interface. We show explicitly in Appendix B.2 that these conditions can be 8 –formulated as δφ + ˜ W + = δφ − ˜ W − (31)˜ W + δφ (cid:48) + − ˜ W (cid:48) + δφ + − g + ˜ ρ + δφ + ˜ W + = ˜ W − δφ (cid:48)− − ˜ W (cid:48)− δφ − − g − ˜ ρ − δφ − ˜ W − . (32)The +/ − signs denote evaluation of a quantity in the upper/lower fluid at the location of theinterface.
3. The Vortex Sheet without Gravity
Equations (21)-(32) 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: V y ( x ) = (cid:40) ¯ V y , x > − ¯ V y , x < , (33)where ¯ V y is a constant. This is known as a vortex sheet approximation. It represents the simplestpossible description of the velocity variation between the two limiting values by essentially ignoringthe details of the transition. Subsequently in §
5, we explore a more realistic model of the velocityvariation, in which the transition occurs in a region of finite radial width. We point out that thevortex sheet approximation is valid for k y δ BL (cid:28)
1. We show in Appendix C.3 that in this limit thedispersion relation for a finite width layer of constant shear reduces to the vortex sheet dispersionrelation.We assume in this section that g = 0, so there is no gravity, and that ρ and s are constantabove and below the interface, but can be discontinuous across it. This case has already beenconsidered by other authors in the past including Landau (1944), Hatanaka (1949), Pai (1954),Miles (1958), and Gerwin (1968). However, the results for the case without gravity will often bereferenced later in the paper and serve as a verification of the formalism we developed in § δφ (cid:48)(cid:48)± + k x, ± δφ ± = 0 , (34)(35)where just as in §
2, the +/ − signs denote the upper/lower fluids respectively. Since k x, ± is constantfor V y ( x ) given by Equation (33), we have that δφ ± ∝ e − ik x, ± x . (36)In general, k x, ± is complex, and the sign is determined by applying the appropriate boundaryconditions. We discuss the boundary conditions shortly, which will also make clear the reason 9 –for the negative sign in the exponential of Equation (36). Plugging the expressions for δφ ± intoEquation (32) gives ˜ W + k x, + δφ + = ˜ W − k x, − δφ − , (37)where we have used δφ (cid:48)± = ik x, ± δφ ± , ˜ ρ ± = ρ ± , and ˜ W (cid:48)± = 0. Using Equation (31) to substitute for δφ − in terms of δφ + , and substituting for ˜ W in terms of W and ρ we have k y k x, + ρ + W = k y k x, − ρ − W − . (38)Introducing the density ratio (cid:15) = ρ + /ρ − , the Mach number in the upper fluid M = ¯ V y /s + , and thephase speed normalized by the sound speed in the upper fluid ϕ = c/s + , we have (cid:15) ( M − ϕ ) (cid:115) ( M + ϕ ) s s − − M + ϕ ) (cid:112) ( M − ϕ ) − . (39)By the definition of the sound speed, s = γP/ρ , the condition of pressure balance everywherethroughout the flow requires that γ − ρ + s = γ − − ρ − s − , where γ + and γ − are the adiabatic indicesabove and below the interface. Assuming γ + = γ − , we have from pressure balance that ( s − /s + ) = ρ + /ρ − = (cid:15) . Thus, the dispersion relation becomes: (cid:15) ( M − ϕ ) (cid:112) ( M + ϕ ) (cid:15) − − M + ϕ ) (cid:112) ( M − ϕ ) − . (40)Miles (1958) has studied Equation (40) and found that the stability criterion, i.e. that ϕ is purelyreal, is given by M > M crit = 12 (1 + (cid:15) / ) / . (41)This shows the surprising result that infinitesimal disturbances are stabilized at high Mach number.However, Fejer & Miles (1963) pointed out that due to Squire’s theorem, it is always possible tochoose an angle θ for the wavevector with respect to the flow velocity such that the projected Machnumber M cos θ is smaller than M crit . Thus, even at high Mach number, the vortex sheet withoutgravity is still unstable to three dimensional disturbances, which are almost perpendicular to thedirection of the flow; these unstable oblique modes resemble classical KH modes. However, in theastrophysical context, a thin disk has a scale height s/ Ω (cid:28) R ∗ , and the wavelength of the obliquemodes will be small relative to the disk scale height only for very small wavelengths λ (cid:28) s/ Ω. Ifthe BL itself has a thickness δ BL (cid:38) s/ Ω, the vortex sheet approximation for the oblique modes isinvalid, since either the modes don’t fit into a disk scale height, or the condition λ (cid:29) δ BL is notsatisfied. 10 – M (cid:29) M (cid:29) M crit , where M crit was defined inEquation (41). Squaring Equation (40), one obtains a sixth order polynomial in ϕ :( (cid:15) − ( M + ϕ ) (cid:15) )( M − ϕ ) − (1 − ( M − ϕ ) )( M + ϕ ) = 0 . (42)This polynomial has two easy to find analytic factors (Miles 1958) ϕ = (cid:40) − M (cid:32) − (cid:15) / (cid:15) / (cid:33) , − M (cid:32) (cid:15) / − (cid:15) / (cid:33)(cid:41) . (43)We now take the limit M (cid:29) M crit , in which case to terms of order O ( M − ), the other four solutionsto the polynomial in Equation (42) are ϕ = (cid:26) M + 1 + (cid:15) M + 1) , M − − (cid:15) M − , − M + (cid:15) / + (cid:15) / M − (cid:15) / ) , − M − (cid:15) / − (cid:15) / M + (cid:15) / ) (cid:41) . (44)It is easy to see that the six roots of the polynomial (42) correspond to sound waves. Startingfrom Equation (25) and rearranging terms we obtain W k y = s ( k x + k y ) . (45)As long as k x is real, this is the dispersion relation for a sound wave, since W is the square ofthe phase velocity in the frame comoving with the fluid. Plugging in the six roots from Equations(43) and (44) into Equation (25), it is straightforward to verify that k x, ± are indeed real for allof them. Each of the four roots in Equation (44) has a further simple interpretation. The firsttwo 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 propagatealmost parallel to the interface in the + y and − y directions in the lower fluid. The two roots inEquation (43) are more difficult to interpret, but the first of these corresponds to a standing wavewhen the two fluids have equal density (i.e. (cid:15) = 1).Since each of the six roots for M (cid:29) M crit yields a real value for k x, ± , the solutions do notdamp away from the interface, and we need to apply radiation boundary conditions at x = ±∞ .The proper procedure is to demand that all waves are outgoing in each fluid in a frame whichis subsonic with respect to the fluid (Miles 1957b), and it is convenient to work in the comovingframe of each fluid. The dimensionless phase velocity in the frame comoving with the upper fluidis ϕ + ,CF = ϕ − M , and the analogous expression for the lower fluid is ϕ − ,CF = ϕ + M . 11 –In order to satisfy Equation (38), k x, + and k x, − must have the same sign, so we must haveeither δφ ± ∝ e i ( | k x, ± | x − ωt ) or δφ ± ∝ e i ( −| k x, ± | x − ωt ) . Moreover, since k x, + and k x, − have the samesign, it is clear that ϕ + ,CF and ϕ − ,CF must have the opposite sign to yield outgoing waves in thecomoving frames of each of the two fluids. Only three of the six roots found above satisfy thiscondition: ϕ l = − M + (cid:15) / + (cid:15) / M − (cid:15) / ) (46) ϕ m = − M (cid:32) − (cid:15) / (cid:15) / (cid:33) (47) ϕ u = M − − (cid:15) M − . (48)We will refer to these three roots as the lower, middle, and upper branches, respectively. Further-more, it is straightforward to check that the solutions which yield outgoing waves in both the upperand lower fluids have δφ ± ∝ e i ( −| k x, ± | x − ωt ) , so k x, ± are positive given our definition (36). We now relax the assumption of M (cid:29) (cid:15) = 1 (Figure 1a) and (cid:15) = .
01 (Figure 1b). At the critical Mach numbergiven by Equation (41), the upper and lower branches merge together in the real plane and bifurcatein the complex plate. These bifurcated solutions turn into the two incompressible KH modes for M (cid:28) M crit . Unlike, the lower and upper branches, the middle branch has no incompressible analogand ceases to be a viable physical solution below a critical Mach number2 M m = 1 + (cid:15) . (49)The reason for this is that below M = M m , k x, ± switches from real to imaginary and the boundaryconditions for the middle branch at x = ±∞ can no longer be satisfied.
4. The Isothermal Vortex Sheet with Gravity
We now go beyond the simplifying assumption of no gravity used in § g is non-zero and constant. We shall shortly assume an isothermal equation of state, but fornow we consider the more general polytropic equation of state of the form P = K ± ρ n in each ofthe two fluids. For any equation of state of this form, we have k g = gs = − ∂P∂ρ dρdr γP = − nγ d ln ρdr = − nγ k s , (50) 12 – - - M j (cid:144) M G = Ε = (a) - - M j (cid:144) M G = Ε = .01 (b)
Fig. 1.— The dispersion relation ϕ/M as a function M for (cid:15) = 1 and (cid:15) = .
01. The solid anddashed curves give the real and imaginary components respectively. The blue, black, and redcurves correspond to the lower, middle, and upper branches respectively. The red and blue curveshave been slightly offset vertically so they do not overlap.which means k s + k g = (1 − γ/n ) k g . If n = γ , corresponding to the case of an adiabatic atmosphere,the second term in parentheses in Equation (21) drops out, but then the third term becomes difficultto treat analytically. However, for an isothermal atmosphere ( n = 1) the sound speed is constant,and both the second and third terms in Equation (21) reduce to a tractable form. Since we are onlyinterested in the gross, qualitative properties of the flow, we assume both fluids are isothermallystratified, since this assumption significantly simplifies the analytical treatment.For n = 1, the second term in Equation (21) becomes − g ˜ ρ k s + k g ˜ W = gk g (1 − γ ) k x k y W (51)We now again assume a vortex sheet velocity profile (Equation (33)), in which case for n = 1, thethird term in Equation (21) is given by˜ W (cid:48)(cid:48) ˜ W = ( √ ˜ ρ ) (cid:48)(cid:48) √ ˜ ρ = (cid:18) − γ k g (cid:19) (52)Putting Equations (51) and (52) into Equation (21), we have δφ (cid:48)(cid:48) + ˜ k x δφ = 0 , ˜ k x ≡ (cid:34) k x (cid:32) − ( γ − (cid:18) k g k y (cid:19) (cid:16) sW (cid:17) (cid:33) − (cid:18) − γ k g (cid:19) (cid:35) . (53) 13 –The perturbation is then given as δφ ± ∝ e − i ˜ k x, ± x . (54)We now comment on the terms present in Equation (53). In the absence of gravity, k g = 0,and we simply have ˜ k x = k x . If k g (cid:54) = 0, then the second term on the right hand side of Equation(53) can be written in a more familiar form as( γ − (cid:18) k g k y (cid:19) (cid:16) sW (cid:17) = (cid:18) Nω − k y V y (cid:19) , (55)where N = ( γ − g s (56)is the Brunt-V¨ais¨al¨a frequency. Thus, the second term on the right hand side of Equation (53)provides a lower frequency cutoff for sound waves at the Brunt-V¨ais¨al¨a frequency. The last termin Equation (53) arises because we have not made the short wavelength approximation k x (cid:29) k g .Now, we use the boundary conditions (31) and (32) to determine the dispersion relation.Substituting δφ (cid:48)± = − i ˜ k x, ± δφ ± , and ˜ W (cid:48)± = (2 − γ ) k g, ± ˜ W ± / − (cid:20) i ˜ k x, + + (cid:18) − γ gs (cid:19)(cid:21) ˜ W + δφ + − g ˜ ρ + δφ + ˜ W + = − (cid:20) i ˜ k x, − + (cid:18) − γ gs − (cid:19)(cid:21) ˜ W − δφ − − g ˜ ρ − δφ − ˜ W − . (57)Next, using the second boundary condition (31) to substitute for δφ − in terms of δφ + , and usingthe fact that ˜ ρ ± = ρ ± at the interface, we have − (cid:20) i ˜ k x, + + (cid:18) − γ gs (cid:19)(cid:21) ˜ W − gρ + = − (cid:20) i ˜ k x, − + (cid:18) − γ gs − (cid:19)(cid:21) ˜ W − − gρ − . (58)We now introduce the dimensionless gravity parameter G = gk y s , (59)which is closely related to the ratio of the wavelength to the pressure scale height, h s . Thus, G ∼ k y h s, + ∼
1, and
G(cid:15) − ∼ k y h s, − ∼
1. Using (cid:15) = ρ + /ρ − = s − /s , and performingsome algebra, we have (cid:34) i ˜ k x, + k y + (cid:18) − γ G (cid:19)(cid:35) (cid:18) k y k x, + (cid:19) (cid:18) W + s + (cid:19) (cid:15) + G (1 − (cid:15) ) = (cid:34) i ˜ k x, − k y + (cid:18) − γ G(cid:15) − (cid:19)(cid:35) (cid:18) k y k x, − (cid:19) (cid:18) W − s + (cid:19) . (60)Setting G = 0 it is clear that we recover the vortex sheet dispersion relation in the absence ofgravity (Equation (38)).We now check Equation (60) by showing that it reproduces the well-known incompressible KHdispersion relations in the limit M (cid:28) M crit , before going on to treat the case of the supersonicvortex sheet with gravity. 14 – We assume that ¯
V /s ± (cid:28)
1, which immediately implies M (cid:28) M crit , and we also assume G (cid:28) G(cid:15) − (cid:28)
1, which means that the wavelength of the perturbation is much smaller than the scaleheight in both the upper and lower fluids ( k g, ± /k y (cid:28) ϕ (cid:28) ϕ(cid:15) − / (cid:28) k x = ± ik y , and given our definition for δφ ± in Equation (54), we must chose the minussign in the upper fluid and the plus sign in the lower fluid to give vanishing solutions at x = ±∞ .Equation (60) then yields − (cid:15) (cid:18) ωk y − ¯ V y (cid:19) (cid:115) − ( γ − G ( M − ϕ ) + gk y (1 − (cid:15) ) = (cid:18) ωk y + ¯ V y (cid:19) (cid:115) − ( γ − G (cid:15) − ( M + ϕ ) . (61)Next, we note that ( γ − G / ( M − ϕ ) = N / ( ω − k y ¯ V y ) and that ( γ − G (cid:15) − / ( M + ϕ ) = N − / ( ω + k y ¯ V y ) , where N is the Brunt-V¨ais¨al¨a frequency for an isothermal atmosphere and wasgiven in Equation (56). We expect | ω ± k y ¯ V y | (cid:38) √ gk , which is the characteristic frequency ofsurface gravity waves. Since N ± ∼ (cid:112) gk g, ± , and we have already assumed k g, ± /k y (cid:28)
1, it followsthat N ± (cid:28) | ω ∓ k y ¯ V y | . Consequently, Equation (61) reduces to − (cid:15) (cid:18) ωk y − ¯ V y (cid:19) + gk y (1 − (cid:15) ) = (cid:18) ωk y + ¯ V y (cid:19) , (62)which is the well known KH dispersion relation for an incompressible fluid in the presence of gravity. As mentioned before in § G (cid:54) = 0.To answer this question, we will use our general dispersion relation (60) in which we willadditionally assume M (cid:29) M crit , since this assumption significantly simplifies the algebra. Sincewe have already obtained solutions for the case G = 0 and M (cid:29) M crit ( § G →
0. In our analysis, we will assumethat the density ratio (cid:15) ≤
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 careshould be taken when determining which solutions are physical and which are not. Miles (1957b) hasshown that for G = 0, the three supersonic KH modes can be understood in terms of sound wavesemitted from the interface between the two fluids. This interpretation is useful as well for the casewith gravity and leads to a couple of insights. First, the amplitude of sound waves is not constant 15 –as they propagate through a stratified medium. Rather, to conserve energy, waves propagatingupward (to lower densities) must increase in amplitude, and those propagating downward (to higherdensities) must decrease in amplitude. Second, since the sound waves are emitted from the interface,if ω has an imaginary component, then the amplitude of the emitted waves changes with time. Thus,if (cid:61) [ ω ] >
0, and there is an instability, then the amplitude of the waves will decay with distancefrom the interface, since the waves emitted in the past had lower amplitude. Conversely if (cid:61) [ ω ] < x → ±∞ . For thelimit G →
0, this means that the physical solutions are still the lower, middle, and upper branchesbut now modified by the presence of a weak gravitational field.We note here that although the dispersion relation (60) is valid for all values of G and not justfor small G , the simple picture of outgoing sound waves in the upper and lower fluids is only validfor G →
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 ω = N below whichpropagation is not possible. Given our definitions of k y and ˜ k x , the frequency of a sound wave in theframe comoving with the fluid is ω ∼ ( k y + ˜ k x ) s . Substituting N ∼ ω and using the definitionof N (Equation (56)) yields G < ( γ − − (1 + ( | ˜ k x | /k y ) ) for sound waves to propagate. If thiscondition is not fulfilled, then the picture of outgoing sound waves is invalid, and the boundaryconditions need to be formulated in a different way, which is beyond the scope of the present work. We begin by considering how the lower wave is modified in the limit G → M (cid:29) M crit .If G = 0 exactly, then ˜ k x = k x , and in the limit G → k x = k x (cid:0) O ( G ) (cid:1) . Keepingterms only to first order in G , Equation (60) becomes (cid:20) i k x, + k y + (cid:18) − γ G (cid:19)(cid:21) (cid:18) k y k x, + (cid:19) (cid:18) W + s + (cid:19) (cid:15) + G (1 − (cid:15) ) = (cid:20) i k x, − k y + (cid:18) − γ G(cid:15) − (cid:19)(cid:21) (cid:18) k y k x, − (cid:19) (cid:18) W − s + (cid:19) . (63)Writing this out explicitly in terms of M and ϕ gives i ( M − ϕ ) (cid:15) (cid:112) ( M − ϕ ) − − i ( M + ϕ ) (cid:112) ( M + ϕ ) (cid:15) − − G (1 − (cid:15) ) = 2 − γ G (cid:18) ( M + ϕ ) (cid:15) − ( M + ϕ ) (cid:15) − − − ( M − ϕ ) (cid:15) ( M − ϕ ) − (cid:19) . (64)Next, we assume that gravity only weakly affects the dispersion relation and make the pertur-bative expansion ϕ = ϕ + ϕ , | ϕ | / | ϕ | (cid:28) , (65) 16 –where ϕ is the solution for G = 0 and M (cid:29) M crit .For the lower branch, we can use Equation (46) for ϕ , and in the limit M (cid:29) M crit we have ϕ ≈ − M + (cid:15) / (66) k x, + k y = (cid:112) ( M − ϕ ) − ≈ M − (cid:15) / (67) k x, − k y = (cid:112) ( M + ϕ ) (cid:15) − − ≈ M − (cid:15) / . (68)Defining M l ≡ M − (cid:15) / , (69)Equation (64) becomes i ( M l − ϕ ) (cid:15) (cid:112) ( M l − ϕ ) − − i ( (cid:15) / + ϕ ) (cid:112) ( (cid:15) / + ϕ ) (cid:15) − − G (1 − (cid:15) ) = 2 − γ G (cid:32) ( (cid:15) / + ϕ ) (cid:15) − ( (cid:15) / + ϕ ) (cid:15) − − − ( M l − ϕ ) (cid:15) ( M l − ϕ ) − (cid:33) . (70)It will turn out (Equation (72)) that ϕ is proportional to G , so working to first order in G isequivalent to working to first order in ϕ . Equation (70) then simplifies to i ( M l − ϕ ) (cid:15) − iM l (cid:15) (1 + 2 (cid:15) − / ϕ )1 + M l (cid:15) − / ϕ + G (1 − (cid:15) ) = 2 − γ G (cid:32) M l (cid:15) − / ϕ M l (cid:15) − / ϕ − (cid:15) (cid:33) . (71)Assuming | M l (cid:15) − / ϕ | (cid:28)
1, and keeping terms to leading order in M l , we can solve for ϕ in termsof G . ϕ ,l ≈ − − γ GM l (cid:15) / i. (72)We immediately see two things from Equation (72). First, ϕ is purely imaginary, and second,if γ < ϕ is negative and the perturbation damps, whereas if γ > ϕ is positive and theperturbation grows. For realistic equations of state, γ ≤ /
3, so the lower wave always damps.
We can find an approximate solution for ϕ in the limit G → M (cid:29) M crit for the middleand upper branches in much the same manner as for the lower branch. The first order correctionfor the middle branch is ϕ ,m = − γ G (1 − (cid:15) / ) (cid:15) / i, (73)and for the upper branch is ϕ ,u ≈ − γ G(cid:15) / M − i. (74) 17 –Just as in the case of the lower branch, ϕ is purely imaginary for both the middle and upperbranches. We see that ϕ is negative for the middle branch if (cid:15) <
1. However, for the upper branchif γ < ϕ is positive and if γ > ϕ is negative. This is opposite from the lower branch meaningthat for any γ (cid:54) = 2 in the limit M (cid:29)
1, one of the two branches is always unstable and the otherone is damped. For a realistic equation of state, γ ≤ / so the upper branch is the unstable one ,and the lower branch damps. We verify Equations (72), (73), and (74) numerically by solving the fully general dispersionrelation (60) and comparing (cid:61) ( ϕ ) with our analytical estimate for the parameters M = 5, (cid:15) = . γ = 5 /
3. We plot (cid:61) ( ϕ ) vs. G in for both our analytical solutions (solid lines) and the onesobtained numerically (dashed lines) in Figure 2. The analytical solution converges to the numericalone in the limit G → - - - - - G I m @ j D Lower Wave (a) - - - - - G I m @ j D Middle Wave (b) G I m @ j D Upper Wave (c)
Fig. 2.— Plots of (cid:61) ( ϕ ) vs. G for the lower (a), middle (b), and upper (c) branches respectively,using M = 5, (cid:15) = . γ = 5 /
3. The solid curves show the numerical solution obtained by solvingEquation (60) and the dashed curves correspond to the approximate solutions from Equations (72),(73), and (74).
5. Sonic Instabilities in a Finite Width Layer of Constant Shear
The calculations presented in § continuously within a narrow shear layer. Unlike the vortex sheet, the finite width shear layer without gravity isknown to be unstable at high Mach number (Glatzel 1988; Choudhury & Lovelace 1984; Ray 1982),and in this case the growth rate of the instability scales inversely with the width of the shear layerand in proportion to the shear, S ∼ Ω K R ∗ /δ BL .In the following, we extend some of the findings of Glatzel (1988) to study sonic instabilities 18 –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 V ( x ) = ¯ V , x > δ BL , ¯ V x/δ BL , − δ BL ≤ x ≤ δ BL , − ¯ V , x < − δ BL (75)and the density profile ρ ( x ) = (cid:40) ρ + , x > − δ BL ρ − , x < − δ BL (76)Pressure equilibrium again requires that ρ + s = ρ − s − , which sets the sound speed everywhere inthe flow, and as before we have (cid:15) = ρ + /ρ − and M = ¯ V /s + .Although we only consider a linearly varying velocity profile in the shear layer, Ray (1982)and Choudhury & Lovelace (1984) have found that different shear profiles are qualitatively similar.Thus, we consider a constant shear to be representative of more general shear profiles. We now study the dispersion relation of the finite width shear layer. In applying the dispersionrelation to the initiation of the BL, we are most interested in the growth rate of the fastest growingmode for M (cid:29)
1. Glatzel (1988) has already obtained the dispersion relation for the case (cid:15) = 1,and using his techniques, we derive the dispersion relation for the case of arbitrary (cid:15) in AppendixC. We also show in Appendix C.3 that the dispersion relation for a finite width shear layer reducesto the dispersion relation for a vortex sheet in the limit k y → δ BL .In Figure 3, we plot the dispersion relation (C20) as a function of wavenumber for the pa-rameters M = 5 and (cid:15) = 1, (cid:15) = . (cid:15) = . (cid:15) = .
01, and (cid:15) = 0. The (cid:15) = 0 case is equivalent tohaving a hard reflecting boundary at x = − δ BL . The modes depicted are the ones that convergeto the upper and lower branches from § k y δ BL (cid:28) (cid:15) <
1, theupper branch always has a larger growth rate than the lower branch, and both the upper and lowerbranches have the same growth rate for (cid:15) = 1. In the (cid:15) = 1 case, we can identify the upper andlower branches as the n ± = 0 decoupled modes in § δ − BL , which implies that the thinnerthe shear layer, the faster the instability proceeds. The key to resolving this apparent controversy 19 –is to consider a mode having k y δ BL (cid:28)
1. Its growth rate is diminished if we decrease δ BL whilekeeping k y constant for M > M crit and becomes vanishingly small if we take the limit δ BL → k y constant,while decreasing δ BL we move leftward along the tail. This means (cid:61) [ ω ] becomes smaller, since ϕ = ω/k y s + is directly proportional to ω for constant k y even as we decrease δ BL . We next pointout that the curve in Figure 3c remains unchanged in shape or amplitude as we decrease δ BL ,keeping k y δ BL constant. Consequently, the value of k y for which the maximum in (cid:61) [ ω ] occurs k y, max ∝ δ − BL and likewise max[ (cid:61) [ ω ]] ∝ δ − BL . Thus, as δ BL is decreased, the instability shifts toshorter wavelengths and becomes more rapid, while at the same time, modes having k y δ BL (cid:28) 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 (cid:29) M crit masks the instability. To independently verify our dispersion relation (C20), we ran direct hydrodynamical simula-tions of shear layers using the Godunov code Athena (Stone et al. 2008) and compared the growthrate of the fastest growing mode predicted by our dispersion relation to that obtained in the sim-ulations. We initialized a flow along the y -direction using the setup described by Equations (75),and for all of the simulations we used γ = 5 / δ BL = 1 for the half width of the shear layer,and periodic boundary conditions in the y -direction; Table 1 summarizes the simulation specificparameters. To seed the instability, we initialized random perturbations to v x of magnitude 10 − 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 4 shows ρv x / (cid:15) = 0and (cid:15) = 1. Label M (cid:15) x -range y -range N x × N y (BC- x
1, BC- x × × × × ×
512 (reflecting, reflecting)Table 1: Parameters for the simulations with Athena. M – Mach number, (cid:15) – density ratio, x, y -range – box size in x, y -direction, N x × N y – number of cells in x, y -direction, (BC- x
1, BC- x
2) –lower and upper boundary conditions in x direction. 20 – Ε = H a L - - - k y ∆ BL R e H j L Ε = H b L k y ∆ BL I m H j L Ε = H c L k y ∆ BL I m H Ω L * ∆ B L (cid:144) s + Ε = .25 H d L - - - k y ∆ BL R e H j L Ε = .25 H e L k y ∆ BL I m H j L Ε = .25 H f L k y ∆ BL I m H Ω L * ∆ B L (cid:144) s + Ε = .1 H g L - - - k y ∆ BL R e H j L Ε = .1 H h L k y ∆ BL I m H j L Ε = .1 H i L k y ∆ BL I m H Ω L * ∆ B L (cid:144) s + Ε = .01 H j L - - - k y ∆ BL R e H j L Ε = .01 H k L k y ∆ BL I m H j L Ε = .01 H l L k y ∆ BL I m H Ω L * ∆ B L (cid:144) s + Ε = H m L - - - k y ∆ BL R e H j L Ε = H n L k y ∆ BL I m H j L Ε = H o L k y ∆ BL I m H Ω L * ∆ B L (cid:144) s + Fig. 3.— Dispersion relations for M = 5 and a finite width shear layer. The top row is for (cid:15) = 1,the second row for (cid:15) = .
25, the third row for (cid:15) = .
1, the fourth row for (cid:15) = .
01, and the bottomrow for (cid:15) = 0. The first column shows the real part of ϕ , the second the imaginary part of ϕ andthe third the imaginary part of ω . The solid line corresponds to the upper branch and the dashedline 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 ( § ABCD0 200 400 600 80010 - - - - - K E - y (a) Fig. 4.— Comparison of the analytically derived growth rates (given by the slopes of the dashedlines) to the ones obtained using Athena (solid lines). The curves A, B, C, D correspond to (cid:15) = 1, . .
1, and 0, respectively (see text for further simulation details) and have been offset verticallyfor clarity.
Panels, (a), (b), and (c) of Figure 5 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 arevery different, illustrating the different instability mechanisms which operate for the M (cid:29) M crit vs. M (cid:28) M crit cases and also for (cid:15) = 1 vs. (cid:15) = 0. The case M (cid:28) M crit , in panel (b) has beenextensively studied (e.g. (Vallis 2006)), so we will not consider it here..For the M (cid:29) M crit case with (cid:15) = 1, the instability is caused by a radiation mechanism . Thismechanism has already been discussed by Glatzel (1988), so we only mention it here briefly. Eachmode of the dispersion relation can be associated with a pseudo-energy that is conserved. If amode has a negative pseudo-energy, then radiation of energy away from the boundary layer regionwill cause the pseudo-energy to become even more negative, amplifying the mode and leading toinstability. The radiation mechanism is responsible for the smooth, broad hump in panels (b) and(c) of Figure 3.For the case of M (cid:29) M crit and (cid:15) = 0, the instability mechanism is quite different. Rather thana broad hump, panels (n) and (o) of Figure 3 exhibit sharp localized peaks. The physical cause ofthese peaks is due to over-reflection of modes that become trapped between the lower boundaryand the critical layer. Thus, we shall call this the over-reflection mechanism .The over-reflection mechanism is discussed in Narayan et al. (1987), who considered a rotatingshear layer adjacent to a reflecting wall. They explained the instability in the (cid:15) = 0 case in terms 22 – (a) (b)(c) Fig. 5.— Panels (a), (b), and (c) show v x during the linear stage of the instability for simulationsA, E, & D respectively and have been rotated by 90 degrees.of the leaking of action density past the corotation radius, which is the analog of the critical layerfor a rotating system. Like the pseudo-energy of Glatzel (1988), the action density is a conservedquantity, and instability occurs because a wavemode that is trapped between the critical layer andthe wall has a negative action density, whereas positive action leaks out past the critical layer. As aresult, 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 3. Thus, we prefer to regard them as separate mechanisms,but point out that in both cases, instability is ultimately caused by radiation emitted from theboundary layer region.We now elucidate some of the properties of the modes that are trapped between the criticallayer and the lower reflecting boundary using a simple model. Consider a sonic mode that is trappedbetween 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 23 –stationary. Let us work in this frame, since it makes the explanations more clear. In order for thewavefronts to be stationary, the equation for a wavefront is dydx = ± (cid:112) M ( x ) − . (77)Here, M ( x ) is the Mach number as a function of x , and the positive sign is for waves propagatingin the − x direction, while the negative sign is for waves propagating in the + x direction. Equation(77) is obtained from the fact that the wavefront propagates at the speed of sound, and in order forit to be stationary, the angle that the wavefront forms with the y -direction obeys sin( θ ) = 1 / M ( x ).Consider now the upper branch from § x c = δ BL ( M − /M .At the critical layer, we set dy/dx = 0, which is just to say that an upward traveling sound waveis reflected there. It immediately follows that inside the region of shear M ( x ) = − M (cid:18) − xδ BL (cid:19) , − δ BL ≤ x ≤ δ BL . (78)(79)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 (75).We now consider the fate of a wavepacket that is trapped between the lower boundary andthe critical layer. Figure 6 shows a set of wavefronts derived by integrating Equation (77) withred segments corresponding to propagation in the + x direction and blue segments correspondingto propagation in the − x direction. Consider a localized wavepacket at point A in Figure 6a. Thewavevector of the wavepacket is oriented perpendicular to the wavefront, and as the wavepacketpropagates towards point B, its wavevector is rotated by the shear. When the wavepacket reachesthe critical layer at point C, it is reflected back towards point D, and the reflected wavepackethas a higher amplitude than the incident one (see Narayan et al. (1987) for the case of a rotatingshear layer). After reflecting off the critical layer, the wavepacket propagates downward to pointD, its wavevector continuing to be rotated in the same sense as before due to the shear. Finallythe wavevector reaches the perfectly reflecting boundary at point E, whereupon the x -componentof its wavevector is reflected and the cycle begins anew. It is thus clear that the repeating cycleA-E will lead to exponential amplification of the wavepacket, due to over-reflection at the criticallayer.The distance in y between points A and E is the wavelength of the longest wavelength mode, λ max ,y , that can be trapped between the lower boundary and the critical layer. A trapped mode canhave a wavelength shorter than λ max ,y , as long as its wavelength satisfies the relation λ y = λ max ,y /n where n is an integer. Thus, the relation for the wavenumber of the n -th trapped mode in the generalcase is k n,y = 2 πnλ max,y . (80) 24 – ABCDEn =
10 5 10 15 20 25 30 35 - - X (a) n =
20 5 10 15 20 25 30 35 - - X (b) Fig. 6.— The n = 1 and n = 2 trapped modes for M = 5 and δ BL = 1. The lower black linedenotes the solid lower boundary, the upper black line denotes the top of the shear layer, and thedashed 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.We indicate using arrows the values of k n,y for the first six trapped modes in panel (o) of Figure 3,and it is clear that they agree well with the locations of the peaks in the dispersion relation, whichlends support for the over-reflection argument. We also plot the n = 2 case in Figure 6b and pointout the similarity between the analytical curves in Figure 6 and the shape of the wavefronts fromsimulation D in Figure 5c.We have discussed the radiation mechanism for (cid:15) = 1 and the over-reflection mechanism for (cid:15) = 0. We now discuss what happens in the more general case of 0 < (cid:15) <
1. As we can see fromFigure 3, both mechanisms operate simultaneously for 0 < (cid:15) <
1. In this case, there is partialreflection at the lower boundary, and as (cid:15) goes from one to zero, less and less of the radiation canescape 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 Figure3 as (cid:15) goes to zero. For (cid:15) = 0, the bump disappears entirely, and the radiation mechanism nolonger operates. On the other hand, as (cid:15) goes from one to zero the reflection mechanism becomesstronger and stronger, since more and more of the energy is reflected at the lower boundary withtotal reflection at (cid:15) = 0. Interestingly, there are still some small-scale wiggles in the dispersionrelation, even for (cid:15) = 1. This is because even if the density is everywhere uniform, radiationcan still partially reflect off the discontinuity in the velocity derivative at x = −
6. Discussion and Conclusions
We have studied supersonic shear instabilities that could drive the turbulence in the BLs ofstars for which the disk is undisrupted by a magnetic field. Our study is aimed mainly at identifyingthe instabilities that lead to the formation of the BL when the disk just touches the surface of thestar. The main result of our work is the identification of two types of instabilities that could operatein 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. Althoughthe vortex sheet is stable above a critical Mach number, the addition of a small amount of gravitydestabilizes it. We have found that the eigenfrequencies of the dispersion relation in the limit G → G . This has theeffect of making the upper branch unstable. We now consider whether the instability of the upperbranch in the limit G → ω ∼ Ω K (cid:15) / i, (81)where Ω 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 (cid:15) (cid:46)
1. It thenfollows from Equation (81) that the characteristic growth rate of the instability is (cid:46) Ω K and isindependent of the wavenumber.If the BL is radially thin, as would be expected during the initiation phase, then the growthrate given by Equation (81) is small compared to the shear rate S ∼ Ω K R ∗ /δ BL (cid:29) Ω K . Consideringthat S is the characteristic growth rate for shear instabilities, if there are other mechanisms forinstability 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 widthshear layer at high Mach number is proportional to S . The sonic instability is similar in natureto the Papaloizou-Pringle instability in that both are global instabilities and cannot be derivedfrom a local analysis. There are two destabilizing mechanisms for the sonic instability. The firstcorresponds to emission of radiation and gives instability over a broad range of wavenumbers. Thesecond corresponds to over-reflection of trapped modes and results in sharply peaked, disconnectedregions of instability in k -space. Because sonic instabilities operate on a much faster timescale thanthe gravity mechanism for a vortex sheet, this makes them an appealing candidate for the initialstages of boundary layer formation, when one might expect large shears to be present.We mention that Alexakis et al. (2004) have considered the Miles instability (Miles 1957a) inthe context of boundary layer formation. Specifically, they invoked the Miles instability to generatemixing of WD and stellar material and explain the enrichment in heavy elements observed in novaexplosions. The Miles instability was proposed as a mechanism for generating waves over waterat low wind speeds and operates through a resonant interaction between the wind and the waterwave. Due to this interaction, a component of the pressure perturbation is created that is in phase 26 –with the velocity of the air-water interface, and much like pumping a swing, swings up the interfaceto large amplitudes. However, the Miles instability was initially introduced as a way of explainingthe formation of waves on water at weak wind speeds, for which the air-water interface is stableto the KH instability. During boundary layer formation, however, large shears are generated, sowe are not in the weak wind regime, and the Miles instability is likely to be swamped by sonicinstabilities.There are two main astrophysical implications of our findings. First, we demonstrate that theinitiation of the BL is likely to take only very short amount time after the first material from thedisk arrives to the stellar surface. This is because the sonic instabilities that we explored in thispaper have an extremely short growth rate, which is very weakly dependent on the density contrastbetween the disk and stellar material. Thus, mixing of the two fluids starts almost immediatelyafter they come into contact.Second, given the efficiency with which the sonic instability operates, it is likely that it mayplay important role also in more developed phases of the BL evolution. As long as the BL pos-sesses some effective ”boundaries” (e.g. sharp changes in the velocity of density behavior) and thegas flow within it is supersonic the purely hydrodynamic sonic instabilities are going to operatein it potentially providing means for continuing mixing and angular momentum transport insidethe layer. Future numerical calculations capable of following the nonlinear development of sonicinstabilities should be able to address this issue.We are grateful to Jim Stone for useful discussions. The financial support for this work isprovided by the Sloan Foundation and NASA grant NNX08AH87G.
REFERENCES
Abramowitz, M., & Stegun, I. A. 1972, Handbook of Mathematical Functions, New York: Dover,1972,Alexakis, A., Young, Y., & Rosner, R. 2002, Phys. Rev. E, 65, 026313Alexakis, A., Calder, A. C., Dursi, L. J., et al. 2004, Physics of Fluids, 16, 3256Armitage, P. J. 2002, MNRAS, 330, 895Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214Balsara, D. S., Fisker, J. L., Godon, P., & Sion, E. M. 2009, ApJ, 702, 1536Bergeron, P., Saffer, R. A., & Liebert, J. 1992, ApJ, 394, 228 This statement is not true for the supersonic KH instability with gravity investigated in §
4, see Equation (74).
27 –Blumen, W., Drazin, P. G., & Billings, D. F. 1975, Journal of Fluid Mechanics, 71, 305Chang, I.-D., & Russell, P. E. 1965, Physics of Fluids, 8, 1018Chimonas, G. 1970, Journal of Fluid Mechanics, 43, 833Choudhury, S. R., & Lovelace, R. V. E. 1984, ApJ, 283, 331Cowling, T. G. 1941, MNRAS, 101, 367Crawford, J. A., & Kraft, R. P. 1956, ApJ, 123, 44Fejer, J. A., & Miles, J. W. 1963, Journal of Fluid Mechanics, 15, 335Fisker, J. L., & Balsara, D. S. 2005, ApJ, 635, L69Fujimoto, M. Y. 1987, A&A, 176, 53Fujimoto, M. Y. 1988, A&A, 198, 163Gerwin, R. A. 1968, Reviews of Modern Physics, 40, 652Ghosh, P., & Lamb, F. K. 1978, ApJ, 223, L83Gilfanov, M., Revnivtsev, M., & Molkov, S. 2003, A&A, 410, 217Glatzel, W. 1988, MNRAS, 231, 795Godon, P. 1995, MNRAS, 277, 157Hachisu, I. 1986, ApJS, 62, 461Hanawa, T. 1987, A&A, 179, 383Hatanaka, H. 1949, J. Soc. Sci. Culture, 2, 3Howard, L. N. 1961, Journal of Fluid Mechanics, 10, 509Inogamov, N. A., & Sunyaev, R. A. 1999, Astronomy Letters, 25, 269Kenyon, S. J., Hartmann, L., & Hewett, R. 1988, ApJ, 325, 231Kippenhahn, R., & Thomas, H.-C. 1978, A&A, 63, 265Kley, W., & Hensler, G. 1987, A&A, 172, 124Kley, W. 1989, A&A, 208, 98Kley, W. 1989, A&A, 222, 141Klu´zniak, W. 1987, Ph.D. Thesis, 28 –Koldoba, A. V., Romanova, M. M., Ustyugova, G. V., & Lovelace, R. V. E. 2002, ApJ, 576, L53Landau, L. 1944, C.R. Acad. Sci. U.S.S.R., 44, 139Livio, M., & Pringle, J. E. 1992, MNRAS, 259, 23PLynden-Bell, D., & Pringle, J. E. 1974, MNRAS, 168, 603Miles, J. W. 1957, Journal of Fluid Mechanics, 3, 185Miles, J. W. 1957, Acoustical Society of America Journal, 29, 226Miles, J. W. 1958, Journal of Fluid Mechanics, 4, 538Miles, J. W. 1961, Journal of Fluid Mechanics, 10, 496Miles, J. W. 1965, Physics of Fluids, 8, 1754Narayan, R. 1992, ApJ, 394, 261Narayan, R., Goldreich, P., & Goodman, J. 1987, MNRAS, 228, 1Pai, S. I. 1954, J. Aero. Sci. 21, 325Papaloizou, J. C. B., & Pringle, J. E. 1984, MNRAS, 208, 721Patterson, J. 1984, ApJS, 54, 443Piro, A. L., & Bildsten, L. 2004, ApJ, 610, 977Piro, A. L., & Bildsten, L. 2007, ApJ, 663, 1252Popham, R., & Narayan, R. 1992, ApJ, 394, 255Popham, R., Narayan, R., Hartmann, L., & Kenyon, S. 1993, ApJ, 415, L127Popham, R., & Narayan, R. 1995, ApJ, 442, 337Pringle, J. E. 1977, MNRAS, 178, 195Pringle, J. E. 1981, ARA&A, 19, 137Pringle, J. E., & King, A. R. 2007, Astrophysical flows by J. E. Pringle, A.R. King; CambridgeUniversity Press, 2007.Ray, T. P. 1982, MNRAS, 198, 617Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337Shakura, N. I., & Sunyaev, R. A. 1988, Advances in Space Research, 8, 13 29 –Shtemler, Y. M., Mond, M., Cherniavskii, V., Golbraikh, E., & Nissim, Y. 2008, Physics of Fluids,20, 094106Spruit, H. C. 1999, A&A, 349, 18Stone, J. M., Gardiner, T. A., Teuben, P., Hawley, J. F., & Simon, J. B. 2008, ApJS, 178, 137Sung, C.-H. 1974, A&A, 33, 99Tassoul, J.-L. 2000, Stellar rotation / Jean-Louis Tassoul. Cambridge ; New York : CambridgeUniversity Press, 2000. (Cambridge astrophysics series ; 36),Vallis, G. K. 2006, Atmospheric and Oceanic Fluid Dynamics, by Geoffrey K. Vallis, pp. 770. Cam-bridge University Press, November 2006. ISBN-10: 0521849691. ISBN-13: 9780521849692,Wheatley, P. J., Mauche, C. W., & Mattei, J. A. 2003, MNRAS, 345, 49
A. Reduction of Equations
Defining S ≡ (cid:36)d Ω /d(cid:36) to be the shear rate, and assuming δ BL (cid:28) R ∗ and Ω ∗ (cid:28) Ω K we havethat 2 B ≈ S ∼ Ω K R ∗ /δ BL (cid:29) Ω K . We next note that | ¯ ω | ≥ ( (cid:61) [¯ ω ]) with exact equality at thelocation of a critical layer, where (cid:60) [¯ ω ] = 0. According to Chimonas (1970), the growth rate of thefastest growing mode in a plane parallel stratified shear flow is of the order (cid:61) [ ω ] ∼ S − N , (A1)where N is the Brunt-V¨ais¨al¨a frequency (Equation (56)). If we assume that shear instabilitiesprovide the turbulence in the BL, then we can estimate that | ¯ ω | (cid:38) S . Using this estimate onthe left hand side of Equation (15), we find that the first term, ¯ ω , dominates the third term, κ ≈ S by a factor of R ∗ /δ BL (cid:29)
1, so we can ignore the third term. Next, we comparethe second to last term ¯ ω ¯ g/s to the last term 2Ω m/(cid:36) on the right hand side of Equation (15).Defining h s ≡ s / ¯ g , dividing the second to last term by the last term, and ignoring constantsof order unity we have ¯ ωR ∗ /h s m Ω (cid:38) SR ∗ /h s m Ω ∼ R ∗ /h s δ BL m . Next, we make the additionalreasonable assumption that the fastest growing mode has m/R ∗ ∼ h − s , since h s is of order thescale height and sets the natural length scale in the problem. Continuing our line of reasoning, wethen have R ∗ /h s δ BL m ∼ R ∗ /δ BL (cid:29)
1. This means that the last term in equation (15) is negligiblecompared to the second to last term and can be ignored. Finally, we compare the first term, C L ¯ ω ,and the second term, 2 Bm/(cid:36) , on the right hand side of Equation (14). Dividing the first term bythe second term, using 2 B ≈ S , | ¯ ω | (cid:38) S and m/R ∗ ∼ h − s , and noting that C L ∼ h − s , we have This preprint was prepared with the AAS L A TEX macros v5.2.
30 – C L ¯ ωR ∗ /mS (cid:38) h s S/h s S ∼
1. Thus, both terms can potentially be of comparable magnitude andboth must be retained.Keeping only the dominant terms for δ BL (cid:28) R ∗ , using (cid:36) ≈ R ∗ , and assuming that ρ (cid:48) /ρ (cid:29) R − ∗ and δu (cid:48) /δu (cid:29) R − ∗ , Equations (14) and (15) finally reduce to Equations (16) and (17). B. Generalized Rayleigh Equation and Boundary ConditionsB.1. Derivation of the Generalized Rayleigh Equation
We start with Equations (19) and (20). Using the definitions from §
2, Equation (19) can bewritten as δu (cid:48) + C ( x ) δu = C ( x ) (B1) C ( x ) = − ( k g + W (cid:48) /W ) (B2) C ( x ) = i (1 − W /s ) k y δP/ρW, (B3)and similarly Equation (20) can be written as δP (cid:48) + C ( x ) δP = C ( x ) (B4) C ( x ) = k g (B5) C ( x ) = − iρ ( k y W + g ( k g + k s )) δu/k y W. (B6)We can use the integrating factor method of solving first order differential equations to make achange of variables and get rid of the δu term in Equation (B1) and the δP term in Equation (B4).Making these changes of variable, we have δq (cid:48) = (1 − W /s ) k y δPρW f , δq ≡ ( iW f ) − δu (B7) δ ˜ P (cid:48) = − if ρ ( k y W + g ( k g + k s )) δuk y W , δ ˜ P ≡ f δP. (B8)We can combine Equations (B7) and (B8) into a single equation for δq : (cid:18) ˜ ρW δq (cid:48) − W /s (cid:19) (cid:48) = ˜ ρ ( k y W + g ( k g + k s )) δq. (B9)This is the same basic form as Chimonas (1970). Since Equation (B9) is a second order differentialequation for δq , we can convert it to standard form (get rid of the first order term) with the changeof variable δφ = ˜ W δq/k y . (B10)Making this change of variable yields Equation (21) which is the generalized Rayleigh equation(Alexakis et al. 2002). 31 – B.2. Derivation of the Boundary Conditions at the Interface
We present here a physically motivated derivation of the interfacial boundary conditions, Equa-tions (31) and (32). The two classical boundary conditions for the vortex sheet that should besatisfied at the interface are (1) that the two fluids should stay in contact at the interface, and(2) that the force exerted on the lower fluid by the upper fluid is equal and opposite to the forceexerted on the upper fluid by the lower fluid. These can be expressed as δξ + = δξ − (B11)∆ P + = ∆ P − , (B12)where δξ is the displacement of the interface in the x -direction, and ∆ denotes the Lagrangian dif-ferential. However, we must formulate both of these boundary conditions in terms of the generalizedstreamfunction perturbation δφ .We begin with the condition δξ + = δξ − . We start with the relation DδξDt = − i ( ω − k y V y ) δξ = δu, (B13)where D/Dt , denotes the Lagrangian derivative. From Equations (B11) and (B13) we have δu + ω − k y V y = δu − ω + k y V y . (B14)Using the definitions of δφ , ˜ W , and ˜ ρ given in §
2, and noting that ˜ ρ = ρ at the interface ( x = 0),we immediately get Equation (31): δφ + ˜ W + = δφ − ˜ W − . (B15)For the second boundary condition, we begin by writing∆ P = δP − gρδξ. (B16)Using Equation (B7) we can substitute for δP in terms of δq , which gives at the interface˜ W k y δq (cid:48) + − g + ρ + δξ + = ˜ W − k y δq (cid:48)− − g − ρ − δξ − . (B17)Using Equations (B13), (B7), and (B10) to substitute for δξ , δq , and δu in terms of δφ , we arriveat Equation (32):˜ W + δφ (cid:48) + − ˜ W (cid:48) + δφ + − g + ˜ ρ + δφ + ˜ W + = ˜ W − δφ (cid:48)− − ˜ W (cid:48)− δφ − − g − ˜ ρ − δφ − ˜ W − . (B18) 32 – C. Finite Width Shear Layer
We derive here the dispersion relation for a finite width shear layer in the absence of gravityand for a constant shear. Our analysis extends the work of Glatzel (1988) to arbitrary densityratios above and below the shear layer. We will use his notation in this section and consider thepressure perturbation δP rather than the generalized stream function δφ . Where appropriate, wedescribe how to transform the results back into the notation used in the body of the text. C.1. Setup of the Problem
Consider the velocity profile V ( x ) = − , x < − x, − ≤ x ≤ , x > , (C1)and the density profile ρ ( x ) = ρ − , x < − ρ , − ≤ x ≤ ρ + , x > . (C2)The adiabatic index and equilibrium pressure are everywhere constant so we have ρ − s − = ρ + s = ρ s . The perturbations are assumed to be of the form δP = δf ( x ) exp[ i ( k y y − ωt )]. We define M = 1 /s to be the Mach number at x = 1 inside the shear layer. We also define (cid:15) − = ρ /ρ − (C3) (cid:15) + = ρ /ρ + (C4)and adopt the following quantities from Glatzel (1988)¯ σ ≡ W = − ¯ ωk y (C5) Q ≡ ¯ σ − / δP (C6) ζ = ik y M ¯ σ (C7) χ = i k y M (C8) µ = 34 . (C9)Glatzel (1988) has shown that Q satisfies Whittaker’s equation d Qdζ + (cid:18) −
14 + χζ + 1 / − µ ζ (cid:19) Q = 0 (C10) 33 –and inside the shear layer has the solution Q = c M χ,µ ( ζ ) + c M χ, − µ ( ζ ) , (C11)where c and c are constants and M χ,µ ( ζ ) is a Whittaker function. Outside the shear layer thevelocity is constant, and the perturbation can be written as Q ∝ (cid:40) exp (cid:2) ± k y (1 − (cid:15) − − M ¯ σ ) / x (cid:3) , x < − (cid:2) ± k y (1 − (cid:15) − M ¯ σ ) / x (cid:3) , x > § C.2. Dispersion Relation
To derive the dispersion relation, we must apply the contact and pressure continuity boundaryconditions at x = ± σ are everywhere continuous,the pressure continuity boundary conditions at x = ± Q in = Q out (C13)where the subscripts “in” and “out” denote evaluation of a quantity inside or outside the shearlayer, respectively. Using the expression (Glatzel 1988) δu = ikρ ¯ σ δP (cid:48) , (C14)the contact boundary conditions at x = ± δu in = δu out (C15) ρ out δP (cid:48) in = ρ in δP (cid:48) out (C16)14 ζ + 1 Q in dQ in dζ = ± ρ in /ρ out (cid:18) ρ out ρ in − χζ (cid:19) / . (C17)In deriving the expression (C17), we have made use of relations (C5-C8), and also the results (C12)and (C13). The sign in Equation (C17) should be chosen on the basis of outgoing waves, andthe lower sign should be chosen at x = 1 and the upper sign at x = −
1. We also point out that ρ in /ρ out = (cid:15) ± at x = ± x = ±
1, and set c = 1, whichsets the normalization. Using the property of the Whittaker function (Abramowitz & Stegun 1972)that ζ dM χ,µ dζ = (cid:18) ζ − χ (cid:19) M χ,µ + (cid:18)
12 + χ + µ (cid:19) M χ +1 ,µ , (C18) 34 –we can derive the dispersion relation (cid:20) − χ + 2 ζ + + 2 (cid:15) + ζ + (cid:16) (cid:15) − − χζ + (cid:17) / (cid:21) M χ, ( ζ + ) + (4 χ + 5) M χ +1 , ( ζ + ) (cid:20) − χ + 2 ζ + + 2 (cid:15) + ζ + (cid:16) (cid:15) − − χζ + (cid:17) / (cid:21) M χ, − ( ζ + ) + (4 χ − M χ +1 , − ( ζ + ) = (cid:20) − χ + 2 ζ − − (cid:15) − ζ − (cid:16) (cid:15) − − − χζ − (cid:17) / (cid:21) M χ, ( ζ − ) + (4 χ + 5) M χ +1 , ( ζ − ) (cid:20) − χ + 2 ζ − − (cid:15) − ζ − (cid:16) (cid:15) − − − χζ − (cid:17) / (cid:21) M χ, − ( ζ − ) + (4 χ − M χ +1 , − ( ζ − ) . (C19)Using the relation between the confluent hypergeometric function and the Whittaker function(Glatzel 1988), Equation (C19) can be written in terms of the confluent hypergeometric functionas¯ σ (cid:20) − χ + 2 ζ + + 2 (cid:15) + ζ + (cid:16) (cid:15) − − χζ + (cid:17) / (cid:21) F ( − χ, , ζ + ) + (4 χ + 5) F ( − χ, , ζ + ) (cid:20) − χ + 2 ζ + + 2 (cid:15) + ζ + (cid:16) (cid:15) − − χζ + (cid:17) / (cid:21) F ( − − χ, − , ζ + ) + (4 χ − F ( − − χ, − , ζ + ) =¯ σ − (cid:20) − χ + 2 ζ − − (cid:15) − ζ − (cid:16) (cid:15) − − − χζ − (cid:17) / (cid:21) F ( − χ, , ζ − ) + (4 χ + 5) F ( − χ, , ζ − ) (cid:20) − χ + 2 ζ − − (cid:15) − ζ − (cid:16) (cid:15) − − − χζ − (cid:17) / (cid:21) F ( − − χ, − , ζ − ) + (4 χ − F ( − − χ, − , ζ − ) . (C20)Equation (C20) is a generalization of the dispersion relations considered by Glatzel (1988).For instance, taking (cid:15) ± → δu ( x = ±
1) = 0, taking (cid:15) ± → ∞ we recover his dispersion relation (5.24) for a vacuumboundary condition δP ( x = ±
1) = 0, and taking (cid:15) ± = 1 we recover his dispersion relation (5.38)for an infinite fluid with uniform density everywhere. C.3. The Vortex Sheet Limit
We can show that Equation (C20) reduces to the dispersion relation for the vortex sheet,(Equation (40)), by taking the limit k y →
0. In this limit, ζ → χ →
0, but ¯ σ , the phasespeed in the comoving frame, is finite, and4 χζ = ( M ¯ σ ) − (C21)is also finite. We also note the property of the hypergeometric function that in the limit ζ → F ( l, m, ζ ) = 1 + lm ζ + O ( ζ ) . (C22) 35 –Using these relations, Equation (C20) can be reduced to the form¯ σ (cid:15) − ( M ¯ σ − (cid:15) − − − / = ¯ σ − (cid:15) + ( M ¯ σ (cid:15) − − / (C23)Equation (40) can then be obtained by setting (cid:15) + = 1, (cid:15) − = (cid:15) , and using the relation M = 1 /s ++