Critical layers and protoplanetary disk turbulence

3.4 MB·

by  Orkan M. Umurhan

Draft version July 12, 2016
PreprinttypesetusingLATEXstyleemulateapjv.01/23/15
CRITICAL LAYERS AND PROTOPLANETARY DISK TURBULENCE
Orkan M. Umurhan1, Karim Shariff and Jeffrey N. Cuzzi
NASAAmesResearchCenter,MoffettField,CA94053,U.S.A
(Dated:)
Draft version July 12, 2016
ABSTRACT
A linear analysis of the zombie vortex instability is performed in a stratified shearing sheet setting for three model
barotropic shear flows. The linear analysis is done by utilizing a Green’s function formulation to resolve the critical
6 layersoftheassociatednormal-modeproblem. TheinstabilityistheresultofaresonantinteractionbetweenaRossby
1 wave and a gravity wave which we refer to as Z-modes. The associated critical layer is the location where the Doppler
0 shifted frequency of a distant Rossby wave equals the local Brunt-V¨ais¨al¨a frequency. The minimum required Rossby
2 number for instability, Ro=0.2, is confirmed for parameter values reported in the literature. It is also found that the
shearlayersupportstheinstabilityinthelimitwherestratificationvanishes. Thezombievortexinstabilityisexamined
l
u inajetmodel,findingthattheinstabilitycanoccurforRo=0.05. NonlinearvorticityforcingduetounstableZ-modes
J isshowntoresultinthecreationofajetflowatthecriticallayeremergingastheresultofthecompetitionbetweenthe

9 vertical lifting of perturbation radial vorticity and the radial transport of perturbation vertical vorticity. We find that
thepictureofthisinstabilityleadingtoaformofnonlinearlydrivenself-replicatingpatternofcreationanddestruction

] is warranted: a parent jet spawns a growing child jet at associated critical layers. A mature child jet creates a next
n
generation of child jets at associated critical layers while simultaneously destroying its parent jet via the Rossby wave
y
instability.
d
Keywords: hydrodynamics, instabilities, protoplanetary disks, turbulence, waves
-
u
l
.f 1. INTRODUCTION nant with distant buoyant critical layers. The resonant
s interaction at the critical layer is unstable and results
c The last few years have enjoyed a revival in the study
in the generation of jet flows similarly oriented with the
i of non-magnetic routes to turbulence in protoplanetary
s original tubular field. These “child” jets are found in
disks. Mostly with the aid of high resolution numerical
y
planes parallel to the mid-plane containing the original
h investigations, several new instability mechanisms and
“parent” tubular surplus field. We note here that this
p processeshavebeenidentifiedandshowntoleadtosome
classofprocesses,i.e.,ofgeneratingamean-flowthrough
[ kind of sustained turbulent activity, including (but not
limited to), the baroclinic instability (Klahr & Boden- the action of a buoyant critical layer, have long been
3 heimer, 2003, Petersen et al. 2007a,2007b, Lesur & Pa- considered relevant to geophysical flow phenomena – for
v example, as a way of transferring wave momentum into
paloizou, 2010, Lyra & Klahr, 2011, Klahr & Hubbard,
2 mean-flowinoceanswithoutgeneratingtoomuchmixing
2014, Lyra 2014 ) the Vertical Shear instability (Gol-
8 (Booker & Bretherton 1967).
dreich & Schubert, 1967, Fricke, 1968, Brandenberg and
3 M13 show that the original parent vorticity field used
Urpin1998,Urpin2003,Nelsonetal. 2013,Stoll&Kley
0 toshowcasethisprocessremains stablyconstituteddur-
2014, VSI hereafter), and recently the self-replicating
0 ing the growth and maturation of the children. Typi-
“zombie vortex instability” (Marcus et al. 2013, Mar-
1. cus et al. 2015, ZVI hereafter). This study is aimed at cally, the vorticity of the child jet grows in magnitude
0 developing deeper intuition of the linear instability as- until it triggers a secondary roll-up type of instability
6 sociated with the self-reproducing ZVI 2 examined in and then gets destroyed. Before these first generation
1 Marcus et al. (2013, M13 hereafter). of jets are destroyed, however, they have enough time
: to spawn a second generation of jets through the same
v M13 uncover the ZVI via numerical simulations in a
resonant critical layer interaction that brought the first
i stably stratified rotating Cartesian box model known as
X generation into existence. The numerical experiments in
the shearing box (Goldreich & Lynden-Bell 1965). The
M13 illustrate a self-reproducing process enveloping the
r model is composed of a vertical component of uniform
a domain, essentially “crystallizing” it by leaving behind
gravity and a basic, vertically uniform, Keplerian (az-
a lattice pattern of dynamically active regions of con-
imuthal) flow field varying linearly with respect to the
centrated vertical vorticity undergoing cycles of creation
(nominal)radialboxcoordinate. M13considerthephysi-
and destruction. The final result resembles a sustained
calresponseofanazimuthally-alignedandrelativelythin
turbulent state.
“tube”ofsurplusverticalvorticityinitiatedinthecenter
If this self-replicating process holds up to independent
ofthebox. Thetubularvortexrespondstoperturbations
scrutiny, its consequences are profound as far as planet
by radiating internal gravity waves which become reso-
formationisconcernedbecausesucharesultingflowstate
hasimportantimplicationswithregardstodustaccumu-
1AlsoatSETIInstitute,189BernardoWay,MountainView,
CA94043,U.S.A.Email: [email protected] lation in protoplanetary disk, magnetic-Dead-Zones of
2 An earlier version of this manuscript referred to the process protoplanetarydisks(Turner,etal.,2014). Forexample,
as the “zombie mode instability” or “ZMI”. We have changed its ascenariooftenquotedisoneinwhichsteadyanticylonic
nomenclature to be consistent with other studies found in recent
vortices can attract particles over time. But if this in-
literatureonthesubject.
2 Umurhan et al.
stability is active, while creating vortices it also destroys 3. The primary mechanism that spawns the jet flow
them and thereby disrupting the envisioned process of is the vertical tilting of perturbation radial vortic-
steady particle concentration. The consequence of this ity although the radial transport of perturbation
would be especially felt during the early stages of cold vertical vorticity can also be important in magni-
disk evolution. Equally important is that triggering this tude depending upon the type of shear flow under
dynamical process requires a minimum source of pertur- consideration.
bationvorticity(seealsoMarcusetal. 2015). Identifica-
tion of this Z-mode self-replicating process in numerical 4. Model jet flows support both the ZVI and the fa-
experiments necessitates high resolution studies because miliarRossbywaveinstability(Lovelaceetal.1999;
the dynamical activity is concentrated in a narrow re- Li et al. 2001, 2000; Umurhan 2010; Meheut et al.
gion centered on the critical layers wherein the region’s 2010, 2012; Les & Lin 2015) (RWI hereafter) the
size depends linearly on the growth rate. A similar de- latter of which induces nonlinear roll-up of the an-
mandonresolutionisneededfortheVSI,butforentirely ticyclonic part of the jet.
different reasons (Nelson et al. 2013; Stoll & Kley 2014;
5. When the amplitude of the jet’s vorticity exceeds
Barker & Latter 2015; Umurhan et al. 2016; Richard et
a minimum value, the jet gets destroyed by the
al. 2016).
RWI.IntheKeplerianflowframe,theanti-cyclonic
In this study we are interested in shedding some light side of the jet experiences destructive roll-up into
onthemechanismofthelinearinstabilityassociatedwith coherent vortices.
this phenomenon and, furthermore, developing some in-
6. Dependinguponthebaseflowfieldconsidered,the
sight as to how the the critical layer instability leads to
analysis indicates that the self-replicating ZVI can
the creation of jets. The linear stability analysis of this
be active for surplus vorticity fields whose Rossby
relatively simplified shearing box setting is a notoriously
numbers are as small as 0.05 – a figure which is
difficult problem even in the case where there is only a
smallerthanpreviouslyanticipatedbyaboutafac-
purely Keplerian flow with no other additional vorticity
tor of four.
surplus (e.g., Dubrulle et al. 2005). A normal mode
reduction of this stripped-down linearized system pro- This study is organized as follows: Section 2 presents
duces a second order differential equation with irregular the shearing box equations assumed for this analysis, in-
singularpointsandmakesgeneratingsolutions,bothan- cluding the scalings leading to them and other assump-
alyticallyandnumerically,extremelydelicateandsubtle. tions and simplifications. The exact nonlinear form of
The challenges are more compounded if one attempts to theverticalvorticityevolutionequationisalsoshown. In
linearly study the same system with the addition of the Section 3 the linearized equations of motion are derived,
surplus vertical vorticity field considered by M13 which assuming a (generalized) vertically uniform surplus ver-
varies both radially and vertically. In this setting, the tical vorticity field. After stating our assumptions, it is
linear stability problem becomes inseparable in both the shownhowtheequationsofmotioncanbewritteneither
radial and vertical coordinates, thereby amplifying the as a set of coupled integral equations or as a second or-
challenges faced by the analyst.3 der differential equation in the radial coordinate of the
Instead, we consider the linear normal mode response shearing box. The system is further decomposed into a
in the stably stratified shearing box containing a verti- potential vorticity formulation which identifies the main
cally uniform surplus vertical vorticity field – either as wave mode involved in the ZVI later explored. Section 4
a vorticity step located at the origin, or a vertically uni- presentsthemodelsurplusvorticityfieldstobetestedin
formshearlayer,oraverticallyuniformdeviationjetflow our analysis: the vorticity step, the shear layer and the
(see Figure 1). By being vertically uniform, these flows asymmetric jet. These three profiles are considered be-
are often referred to as “barotropic”. These test models causethevorticitystepshowcasesthecriticallayermech-
arechosenbecauseofanalyticalandnumericaltractabil- anisminitspuresttheoreticalformwhiletheshearlayer
ity. While these model flow profiles are different than is similar to the flow considered in M13 while the jet
what was considered in M13, we find that the linear dy- flow represents the flow to emerge from the critical layer
namicalpropertiesoftheseflowprofilesarequalitatively instability itself.
similar to the processes involved in the self-reproducing Section 5details the relevant length andtime scales in
mechanism reported in M13. We find: the problem and identifies the independent parameters
of the problem depending upon which profile is consid-
1. TheZVIischaracterizedprimarilyasthenearres-
ered. Section 6 discusses the general solution method
onance between the Doppler shifted frequency of
implemented. Because of the aforementioned difficulties
the Rossby wave associated with the distant sur-
inherenttosolutionsofdifferentialoperatorswithexplic-
plusvorticityfieldandtheBrunt-V¨ais¨al¨afrequency
itly appearing, non-removable, irregular singular points,
at the buoyant critical layer. Without the Rossby
we opt for solving the coupled integro-differential equa-
wave the buoyant critical layer is not excited.
tions developed in Section 3.
Section 7 concentrates on the normal mode results for
2. The nonlinear outcome of the linear instability is the three model shear flows considered. We detail the
to drive into existence a jet-like vertical vorticity propertiesoftheZombiemodeinstability(hereintermed
field at the critical layer(s). Z-modes)inthesimplevortexstepprofile. Wealsoshow
howtheasymmetricjetprofilealsosupportsboththeZ-
3 This is also a feature of the VSI; see the discussion in Barker modeinstabilityandtheRWI(thelatterdetailedinAp-
&Latter(2015);Umurhanetal.(2016). pendixB).Forthelatter,weanalyticallydevelopgrowth
Critical layers and disk turbulence 3
rates and conditions for marginality. Section 8 examines meansthatitisbuoyantlystableandthemodelisstably
the nonlinear forcing implied by an unstable Z-mode. stratified against buoyant instabilities.
We explicitly demonstrate how it drives into existence Inderivingtheboxequations,thereisalwaysacompo-
a jet flow at the location of the critical layer. Section nentrepresentingtheradialgradientofthegravitational
9 summarizes our main results and interprets the self- potentialofthecentralstar,andthisappearsintheform
replicatingdynamicalprocessofjetcreation/destruction oftheexpression3Ω2xfoundontherighthandsideofthe
0
in terms of the physical results garnered by the exami- x-momentumequation. Intheabsenceoftimedependent
nation peformed in this study. dynamics, the balance of this term with the correspond-
ingCoriolisterm−2Ω v admitsthebasicKeplerianflow
0
2. EQUATIONS state vk, i.e.
Following M13, we consider dynamics in the so-called 3
shearing box which represents dynamics in the frame of −2Ω0vk =3Ω20x, −→ vk =−2Ω0x.
a Cartesian “box”-section of an accretion disk which ro-
tates around the central star at a distance of R with a However, in this study, we shall consider arbitrary plane
0
rotationrateΩ andcorrespondingrotationvectorgiven parallel shear flows v0(x) containing both vk plus some
by ˆzΩ , , see G0oldreich & Lynden-Bell (1965), which is deviation/surplus shear flow V , the latter of which is
aligned0 with the disk vertical direction. We take the similarly scaled by the local sou0nd speed. Such a devia-
equivalent radial direction to be along xˆ and the disk’s tion flow is a steady solution of this system provided the
azimuthaldirectionvectorasyˆ. Lengthscalesinthebox pressure field shows some radial variation in its steady
areinunitsofthelocalverticalscaleheightH . Thebox state as well:
0
equations are considered valid because the ratio H /R v (x)=v +V , (6)
0 0 0 k 0
is generally small for cold accretion disks. We make the
whereV =(∂ Π )/(2Ω )isthedeparturefromthebasic
additionalassumptionofaBoussinesqfluidand,thereby, 0 x 0 0
Keplerian flow v . We detail in Section 4 the two types
suppressingacousticmodes. Theequationsofmotionare k
of V we analyze herein. We often refer to these flow
therefore 0
fields as barotropic flow profiles since we assume there is
(∂ +u∂ +v∂ +w∂ )u−2Ω v=−∂ Π+3Ω2x,(1) no vertical variation of V .
t x y z 0 x 0 0
(∂ +u∂ +v∂ +w∂ )v+2Ω u=−∂ Π, (2) Although this is primarily a linear stability analysis,
t x y z 0 y
(∂ +u∂ +v∂ +w∂ )w=−∂ Π+gΘ, (3) thenonlinearverticalvorticityequationisusefulinhelp-
t x y z z
ing us to interpret the consequences of the results we
(∂ +u∂ +v∂ +w∂ )Θ+βwΘ=0, (4)
t x y z report here in this work. We define the vertical compo-
∂xu+∂yv+∂zw=0. (5) nent of the fluid vorticity as
The quantities u,v,w denote the components of velocity ζ ≡∂ v−∂ u.
x y
in the radial, azimuthal and vertical directions respec-
tively. In the shearing box, the u,v and w velocity com- Inordertoobtainanevolutionequationforthisquantity
ponents are scaled by the local sound speed and this is weoperateonequation(2)by∂xandsubtractfromitthe
possible only because the box equations are for dynam- result of operating on equation (1) by ∂y. Rearranging
ics in a frame moving with the exceedingly faster Kep- the result and making use of (5) we find
lerian velocity. To facilitate comparison with the results (cid:16) (cid:17)
of M13, we make the similar assumption that is made ∂t+u∂x+v∂y+w∂z ζ =
inthatstudywhereintheverticalcomponentofthedisk (cid:104) (cid:105)
gravitational field g is taken to be a constant. (2Ω0 +ζ)∂zw− (∂xw)·(∂zv)−(∂yw)·(∂zu) . (7)
The quantity Θ can be interchangeably interpreted
as either an entropy or a buoyancy; for this analysis What this equation says is that the vertical vorticity is
we consider it as the latter, in which case Θ is a non- advected (in a Lagrangian sense) by the flow field and,
dimensional quantity by definition. Equation (3) says had the right hand side of the above equation been zero,
the entropy/buoyancy field is passively advected by the the local value of ζ would be preserved along the way.
total flow with no sources or sinks meaning to say that However, at any given position, ζ can change because of
theperturbationsareadiabatic. Implicitly,wearesaying the two effects appearing on the right hand side of equa-
that the actual cooling times of the protoplanetary disk tion (7). The first of these, arising from the expression
gasismuchlongerthanthegrowthratescalculatedlater (2Ω +ζ)∂ w, is the familiar effect of vertical stretching
0 z
on. 4 Because we are in the Boussinesq approximation, wherein the total vorticity contained in a moving fluid
the pressure field p is rewritten in terms of the quantity elementcanspinupordowndependinguponwhetheror
Π≡p/ρ , as the density ρ is taken to be a constant in notthereisconcurrentverticalstretchingintheflowfield
0 0
the mean state. The quantity β, which is in units of in- itself. This is the usual effect known from the study of
verse length, represents the vertical gradient of the disk Taylor columns. The remaining terms can be rewritten
buoyancy (itself a dimensionless quantity) and as with g to represent the well-known vortex tilting effects,
above, we assume it is constant. As taken here, β > 0 (cid:104) (cid:105)
=− (∂ w)·(∂ v)−(∂ w)·(∂ u) ,
x z y z
4 This is in contrast with the conditions giving rise to the VSI (cid:16) (cid:17) (cid:16) (cid:17)
which requires cooling times to be very short compared to both = ∂ w−∂ v ∂ w+ ∂ u−∂ w ∂ w,
y z x z x y
thelocalrotationtimeofthediskandtheassociatedgrowthrates.
SeerecentdiscussioninLin&Youdin(2015). =ζ ∂ w+ζ ∂ w, (8)
x x y y
4 Umurhan et al.
where the x and y directed vorticities are defined by noting here that the vertical velocity is equated with
ζ ≡∂ w−∂ v andζ ≡∂ u−∂ w respectively. Forex- the horizontal divergence based on the above form. The
x y z y z x
ample, the expression ζ ∂ w describes the rate in which equationsofmotionmaynowbeformallyreducedbyone
x x
the x-directed component of the vorticity is turned into order in time derivatives to get the following
vertical vorticity due to the shear along the x-direction (cid:16) (cid:17)
of the vertical velocity field w. A similar interpretation ∂ +iαv ζ(cid:48)=−(2Ω +v )D(cid:48)−v u(cid:48), (17)
t 0 0 0x 0xx
holds for the term ζ ∂ w. 5 We find below that when
y y (cid:16) (cid:17)
the so-called zombie modes appear, they lead to nonlin- ∂ +iαv D(cid:48)=−m2Π(cid:48)−θ(cid:48), (18)
t 0
eargenerationofverticalvorticitythroughthecombined
(cid:16) (cid:17)
actionofthetwovortextiltingtermsshowninexpression ∂ +iαv θ(cid:48)=βgD(cid:48), (19)
(8). Wefindthatbetweenthetwo,thevortextiltinggen- t 0
erationofverticalvorticityisusuallybytheζx∂xw term. where we have defined the quantity θ(cid:48) ≡ imΘ(cid:48). Note
that throughout this study we assume stable stratifica-
3. LINEARIZATION tion which means that the Brunt-V¨ais¨al¨a frequency N
B
We linearize around a plane-parallel shear state by in- is real, i.e. N2 ≡gβ >0. Asfar asnomenclature iscon-
B
troducing the form cerned, we mostly dispense with using the symbol NB
for the Brunt-V¨ais¨al¨a frequency and, instead, retain for
(cid:16) (cid:17)T
u,v,w,Π,Θ (cid:55)→(0,v (x),0,Π (x),0)T+ its designation the expression gβ.
0 0 The above third order system (in time) is supple-
(cid:16) (cid:17)T mented by the diagnostic condition relating the pertur-
u(cid:48)(x,t),v(cid:48)(x,t),w(cid:48)(x,t),Π(cid:48)(x,t),Θ(cid:48)(x,t) ei(αy+mz)+c.c.,
bation pressure to the other quantities appearing, and is
given as the solution of
where the T superscript means transpose, v (x) =
0 (cid:16) (cid:17)
a−r3bxitΩra0r/y2+baVro0t(rxo)p,iicnswhheaicrhpVro0(fixle) oisftohueracfhooreomsinengt.ioTnehde ∂x2−α2−m2 Π(cid:48) =2Ω0ζ(cid:48)+θ(cid:48)−2iαv0xu(cid:48). (20)
aboveansatzinsertedintotheequationsofmotionyields We observe that the vertical vorticity field is driven by
thefollowingpartialdifferentialequationsforthepertur- verticalstretching(D(cid:48))andradialadvectionofthemean
bation quantities, vorticity gradient (v u(cid:48)).
0xx
(cid:16) (cid:17) In the form as developed here, the streamfunction
∂t+iαv0 u(cid:48)−2Ω0v(cid:48)=−∂xΠ(cid:48), (9) and velocity potential solutions are written in terms of
(cid:16) (cid:17) Green’s functions, i.e.
∂ +iαv v(cid:48)+(2Ω +v )u(cid:48)=−iαΠ(cid:48), (10)
t 0 0 0x (cid:90)
(cid:16) (cid:17) ψ(cid:48) = G (x,x(cid:48))ζ(cid:48)(x(cid:48))dx(cid:48),
∂ +iαv w(cid:48)=−imΠ(cid:48)+gΘ(cid:48), (11) ψ
t 0 (cid:90)
(cid:16)∂t+iαv0(cid:17)Θ(cid:48)=−βw(cid:48), (12) φ(cid:48) = Gφ(x,x(cid:48))D(cid:48)(x(cid:48))dx(cid:48), (21)
∂ u(cid:48)+iαv(cid:48)+imw(cid:48)=0. (13) whereG(x,x(cid:48))istheappropriateGreen’sfunctionassoci-
x
atedwiththetwo-dimensionalLaplaceoperator∂2−α2,
The subsequent analysis exploits the incompressible na- x
subject to appropriate boundary conditions (see more be-
ture of the disturbances. This is done by reducing the
low),
linearizedequationsofmotionintovorticity/dilatational (cid:16) (cid:17)
form,i.e.,bydefining(respectively)theverticalvorticity ∂2−α2 G =δ(x−x(cid:48)), (22)
x φ,ψ
perturbation and horizontal velocity divergence accord-
ing to whereδ(x)istheDiracdeltafunction. Similarly,thecor-
responding solution of the pressure fluctuations is given
ζ(cid:48) ≡∂ v(cid:48)−iαu(cid:48), D(cid:48) ≡∂ u(cid:48)+iαv(cid:48). (14)
x x by
Thehorizontalvelocityfieldsmaybewritteninvelocity- (cid:90) (cid:110) (cid:111)
potential/streamfunction form, i.e., Π(cid:48) = G (x,x(cid:48)) 2Ω ζ(cid:48)(x(cid:48))+θ(cid:48)(x(cid:48))−2iαv (x(cid:48))u(cid:48)(x(cid:48)) dx(cid:48).
Π 0 0x
u(cid:48) =−iαψ(cid:48)+∂ φ(cid:48), v(cid:48) =∂ ψ(cid:48)+iαφ(cid:48) (15) (23)
x x
The Green’s function G is the solution of
where ψ(cid:48) and φ are the streamfunction and velocity- Π
(cid:16) (cid:17)
potentials,respectively. Theaboveformulationautomat- ∂2−k2 G =δ(x−x(cid:48)), k2 ≡α2+m2. (24)
x Π
ically satisfies the incompressibility equation (14) pro-
vided The system of equations (17–20) together with their
associated diagnostic relationships (15) and (16) is an
ζ(cid:48) =(∂2−α2)ψ(cid:48), D(cid:48) ≡−imw(cid:48) =(∂2−α2)φ(cid:48), (16)
x x integro-differential system which must be solved sub-
ject to boundary conditions in the radial direction. The
5 Thisisawell-knownfeatureinincompressiblerotatingflows
Green’s function strategy adopted here has been used in
where,ingeneralvectorform,thevorticityequationreads
otherdiskstudies(e.g.,Dubrulle&Knobloch,1992). For
dω this study we report upon solutions in which all quanti-
=(ω·∇)u,
dt ties exponentially decay to zero as x → ±∞. We note
where ω = 2Ω ˆz+ζ with ζ ≡ ∇×u with u the vector velocity that for the unstable localized modes which are the sub-
0
field. ject of this study, we have checked and verified that the
Critical layers and disk turbulence 5
resultsreportedhereareinsensitivetowhetherornotthe in which the following normal mode ansatz has been as-
perturbations are periodic (on scale 2L) or if there are sumed
no-normalflowboundaryconditionsimposedatx=±L. Π(cid:48) =Πˆ(x)e−iωt, (32)
The above statement becomes robust so long as the hor-
andsimilarlyfortheotherquantitiesappearing. Here,ω
izontal domain is large enough with, typically speaking,
is the unknown complex normal mode frequency. Inser-
L≥π. ThelocationoftruncationscaleLdoesnotinter-
tionofthisformintothefundamentalperturbationequa-
fere with the existence and/or expression of an unstable
tions (9–14) followed by some manipulation shows that
mode.
the following relationships hold between various quanti-
ties:
3.1. An alternative form in terms of a perturbation
potential vorticity gβm2 iσ iσm2
θˆ= Πˆ, Dˆ = θˆ= Πˆ,
It is worth noting that the processed linearized equa- σ2−gβ gβ σ2−gβ
tions(17–19)mayberecastinsteadintermsofapertur- 2Ω αΠˆ +σ∂ Πˆ
bation vertical potential vorticity, Ξ(cid:48), defined by uˆ=i 0 x , (33)
ω2−σ2
ε
(cid:18) (cid:19)
2Ω +v
Ξ(cid:48) ≡ζ(cid:48)+ 0 0x θ(cid:48). (25) as well as
βg
ω2
2Ω ζˆ=− εθˆ−2Ω v iuˆ/σ, (34)
In many shear flow applications, following the behavior 0 gβ 0 0xx
of the potential vorticity is an extremely useful diagnos-
tic (Hoskins et al. 1985). The perturbation equations whe(cid:0)re we have de(cid:1)fined σ(x) ≡ αv0(x)−ω and ωε2(x) ≡
maynowberewritteninthefollowingalternativeformu- 2Ω 2Ω +v (x) - the latter of these is the local disk
0 0 0x
lation: epicyclic frequency. Rewriting equation (20) in terms of
(cid:16) (cid:17) the normal mode ansatz means that
∂t+iαv0 Ξ(cid:48)=−v0xxu(cid:48), (26) (cid:16)∂2−α2−m2(cid:17)Πˆ =2Ω ζˆ+θˆ−2iαv uˆ, (35)
(cid:16) (cid:17) x 0 0x
∂ +iαv ζ(cid:48)=−(2Ω +v )D(cid:48)−v u(cid:48), (27)
t 0 0 0x 0xx and making use of the above relationships and some ad-
(cid:16) (cid:17) gβ (cid:16) (cid:17) ditional reduction shows that
∂ +iαv D(cid:48)=−m2Π(cid:48)+ ζ(cid:48)−Ξ(cid:48) ,(28)
t 0 2Ω +v (cid:16) (cid:17) ω2−gβ
0 0x ∂2−k2 Πˆ = ε m2Πˆ
with x gβ−σ2
(cid:16)∂x2−k2(cid:17)Π(cid:48) =(cid:20)2Ω0(22ΩΩ0 ++vv0x)−gβ(cid:21)ζ(cid:48) +σ2−1 ωε2 (cid:18)2αv0x − 2Ω0σv0xx(cid:19)(cid:16)2Ω0αΠˆ +σ∂xΠˆ(cid:17)(3.6)
0 0x
gβ This system is subject to the same boundary conditions
−2iαv u(cid:48)+ Ξ(cid:48).(29) as outlined above. 6
0x 2Ω +v
0 0x InspectionoftheODE(36)showsthatthereexiststhe
We note immediately that in the event that v is con- possibility of the system supporting several critical lay-
0x
stant (e.g. pure Keplerian flow with V =0), the system ers(Drazin2002;Drazin&Reid2004). Locationswhere
0
materially conserves the perturbation potential vorticity thedenominatorsvanisharecandidateirregularsingular
withthebasicconstantshearstate. Inslightlymoregen- points. Of interest to us are the points associated with
eral terms, if V (x) = Ω x, where Ω = constant (units expressions multiplying the pressure Πˆ (and, not any of
0 x x
of inverse time), then its derivative expressions). These points become irregu-
larifthedenominatorsortheassociatedexpressionscan
(cid:16) (cid:17) 3
∂ +iαqx Ξ(cid:48) =0, q =− Ω +Ω , (30) pass through zero linearly with respect to variations in
t 2 0 x x. Several candidate points are identified: The first of
these is the classical one associated plane parallel shear
which means that a potential vorticity perturbation is
flows and occurs at points x in which Re(σ)=0, i.e.
advected by the composite background shear flow field pp
= qx. In this way the system becomes relatively trans-
αv (x )−Re(ω)=0.
parent. In the event we consider a flow field such that 0 pp
v = 0, together with Ξ(cid:48) = 0 initially, then the linear In classical plane parallel shear flow problems, such crit-
re0xsxponse is of pure inertial gravity waves. ical layers activate only if both v0xx (cid:54)=0 and if viscosity
is included (Drazin & Reid, 1981). The normal modes
3.2. As a single second order differential equation and associated with these critical layers are called Tollmien-
critical layers Schlichting waves (TS waves) and they are unstable for
wide values of Reynolds numbers. TS-waves cease to
Previoustreatments(e.g.,Dubrulleetal. 2005)ofthis
system as a normal mode problem have instead turned 6 Had we considered no-normal flow boundary conditions at
equations (9 – 14) into a single second order differential x = ±L, then the imposition of an impenetrable flow boundary
equation. Indeed,combiningtheseequationsintoasingle conditionatthesetwolocationswouldbethesameasimposing
oneforthenormal-modepressureperturbationresultsin 2Ω0αΠˆ+σ∂xΠˆ =0,
(cid:0)∂x2−k2(cid:1)Πˆ =2Ω0ζˆ+θˆ−2iαv0xuˆ, (31) atx=±Lprovidedσ2−ωε2(cid:54)=0attheboundaries.
6 Umurhan et al.
be normal modes in the exactly inviscid problem. In Drazin & Reid, 1981) in which constant opposite veloc-
monotonic shear flows (flows in which the shear velocity itylayerssandwichauniformshearlayerofthickness2∆
strictlyincreasesordecreases)theretendstobeonlyone centered at x=0. Thus we have
such critical point x . (cid:26) (cid:20) (cid:18) (cid:19)(cid:21) (cid:20) (cid:18) (cid:19)(cid:21)(cid:27)
pp ω (cid:15) x−∆ x+∆
Withtheinclusionofbuoyancyunderstablestratifica- V = a ln cosh −ln cosh (39.)
0 2 (cid:15) (cid:15)
tion, in which gβ >0, two more critical layers emerge at
points x± which we henceforth refer to as buoyant crit-
bg
ical layers to distinguish them from the others. These
The asymmetric jet. In the absence of a Keplerian shear
locations are associated with the first term√on the right
thetriplejumpprofiledescribedbelowwillphysicallyre-
hand side of equation (36), where Re(σ)± gβ =0, i.e.
semble that of a jet. However, taken in aggregate with
αv (cid:0)x±(cid:1)±(cid:112)gβ−Re(ω)=0. (37) the background Keplerian shear, the composite flow de-
0
bg scribes a shear with a weak jet-like undulation atop of
The instability identified in Section 7 pertains to the ac- it. Nevertheless, we refer to this as “jet” in so far as the
tion of these buoyant critical layers. An inspection of deviation flow V resembles one. The flow V involves
0 0
equation(36)showsthatthesecriticallayersareintrinsi- three steps in the mean vorticity profile located at posi-
cally a three-dimensional phenomenon as the expression tions −∆−,0 and ∆+. The vorticity of V in both the
0
m2Πˆ/(gβ−σ2) requires vertical perturbations (m(cid:54)=0), regions x > ∆+ and x < −∆− are zero while the vor-
azimuthal perturbations (α (cid:54)= 0) and a radial variation ticity is given (approximately) to be ωa− for the region
in the shear flow v (x). −∆− < x < 0 and ωa+ for the region 0 < −∆+ < x.
Finally, equation0 (36) also admits the possibility of Thus we have
two more critical layers associated with those points
where the denominator of the second term on the right ω− (cid:40) (cid:20) (cid:18)x+∆−(cid:19)(cid:21) (cid:104) (cid:16)x(cid:17)(cid:105)(cid:41)
hand side of equation (36) is equal to zero, that is, V = a (cid:15) ln cosh −ln cosh
Re(σ)±ω (x) = 0. These points are often referred to 0 2 (cid:15) (cid:15)
ε
in the astrophysical literature as the corotation points ω+ (cid:40) (cid:20) (cid:18)x−∆+(cid:19)(cid:21) (cid:104) (cid:16)x(cid:17)(cid:105)(cid:41)
or Lindblad resonances of a disk (e.g., Papaloizou &
− a (cid:15) ln cosh −ln cosh (40.)
Pringle, 1984). While at first glance it may seem that 2 (cid:15) (cid:15)
the corotation points ought to be important to the dy-
namicsinvestigatedhere,wefindinourresultsthatthey Most of the variation of this mean velocity field is con-
in fact play little role. The main reason appears to be fined to within −∆− < x < ∆+ centered x = 0 and
because the expressions found in the numerator associ- once one has moved sufficiently far from this region the
ated with this term, i.e. 2Ω αΠˆ + σ∂ Πˆ, always gets flow returns to being largely Keplerian. For ω+ = ω−
0 x a a
nearly as small as the value of the denominator in this and ∆+ =∆−, in the limit (cid:15)→0 the resulting symmet-
region. This would cease to be the case if, for example, ric profile recovers the so-called triangular jet (Drazin,
self-gravitational physics were included in the analysis. 2002).
While further elucidation to clarify the inactivity of the
Lindblad resonances are surely in order, because these Wereferencethenatureofthevelocityprofileinterms
layers play no role in our results these matters are not of deviations from the Keplerian state according to the
considered forthwith. followingconvention: ifdV /dx<0thenwesaythatthe
0
profile (or the part under consideration) is anticyclonic
4. CHOICEOFBAROTROPICVELOCITYPROFILES, with respect to the Keplerian shear, while if dV /dx>0
SOMENOMENCLATURE then we say that profile (or part under consid0eration)
We consider three barotropic shear profiles V (x) that is similarly cyclonic. Often times, we will consider the
are departures from pure Keplerian flow (Figur0e 1): Rossby number, defined by
1 dV
The vorticity step. This single vorticity jump profile Ro≡ 0,
2Ω dx
(Figure 1a) , and its derivatives, are given by 0
which is the same definition as used by M13 and Marcus
(cid:104) (cid:16)x(cid:17)(cid:105)
V =ω (cid:15)ln cosh . (38) et al. (2015). We use Ro to quantify the change of the
0 a (cid:15) vorticityinpartorinthewholeofaprofilebeingtested.
The model profile is governed by two parameters, ω For example, the jump in the vorticity across the single
(unitsofinversetime)and(cid:15)(unitsoflength). Thelattear stepvorticitydefectisω whiletheeffectiveRossbynum-
a
quantitycontrolshowsharplythesheartransitionoccurs ber characterizing deviations of the shear flow about the
(aroundx)whiletheformerdictatestheshearprofilefor Keplerian state is Ro = O(|ω |/2Ω ). Generally speak-
a 0
values of |x|(cid:29)(cid:15), i.e. V ≈ω |x| and V ≈ω sgn(x) as ing, jumps in the vorticity relate to the Ro in the way
x→±∞. The useful fe0atureaof this typ0xe of maodel (and indicated and we think of them in terms of this equiva-
for the others described hereafter) is that as (cid:15) → 0 the lence hereafter.
profile resembles a piecewise linear velocity field. In this We perform a compreshensive analysis of these three
sense, as (cid:15)→0 V ≈2ω δ(x). profiles for the following reasons: The vorticity step
0xx a
showcases the essential mechanism of the critical layer
Theshearlayer. Thedoublejumpinvorticityisafacsim- mechanism. We consider the shear layer because it is
ile of the classic Rayleigh shear profile (Rayleigh, 1880, the barotropic analog of the flow profile considered in
Critical layers and disk turbulence 7
(a)SingleJump (Vortex Defect) (b)DoubleJump (ShearLayer) (c)TripleJump (AsymmetricJet)
−1 −1.45 −1
∼2∆
−1.2 −1.5 −1.2
dv−1.4 ∼ωa −1.55 −1.4 ∼∆+ ∼ωa+
0
dx −1.6
−1.6 ∼−ωa −1.65 ∼ωa −1.6 ∼∆−
−1.8 ∼ω−
−1.8 a
−1.7
−2
−1 0 1 −1.75 −2
−1 0 1 −1 0 1
x
Figure 1. Thetotalshearprofile(dv /dx=−3Ω /2+dV /dx)ofthetwobarotropicvelocityfieldsdetailedinSection4: (a)The“vortex
0 0 0
step”, (b) the “shear layer”, and (c) the “asymmetric jet”. The shear is scaled in units of Ω while the radial length x is scaled to the
0
boxscale L. Notethatall shearfieldsshownincludethe backgroundKeplerian state(=−3Ω /2)whichisindicatedforreferencebythe
0
dashed horizontal line in panel (a). The deviation from Keplerian shear profile shown in panel (b) is anti-cyclonic while in panel (c) it is
anti-cyclonicintheregion−∆−<x<0providedω−<0anditiscyclonicfor0<x<∆+ providedω+>0.
a a
M13. The analysis of the asymmetric jet is included be- profile problem with the addition of ∆/(cid:15), describing the
cause jets are flows that emerge from the critical layer ratio width of the shear layer to the size of its transi-
instability. tionzone. Theresultsforthejetprofilearedescribedby
seven parameters: α(cid:15), m(cid:15), gβ/Ω2, ∆±/(cid:15) and ω±/Ω .
5. RELEVANTLENGTHANDTIMESCALES,AND Wehavecheckedandverifiedth0atthenumericaalre0sults
DEPENDENCEOFRESULTSONPROBLEM
wedevelop(anddescribehereafter),reproducesinvariant
PARAMETERS
solutions for invariant values of the aforementioned pa-
The linearized equations and the flow fields we have
rameters. Forexample,forthevorticity-stepproblemwe
adopted to study are presented in dimensional units
verifythattheeigenvalueswefinddependstrictlyonm(cid:15),
mainly in order to facilitate comparison between these
α(cid:15)aswellasgβ/Ω2 andω /Ω only. Avariationin(cid:15)ac-
theoretical results and numerical results reported in the 0 a 0
companiedbyadjustmentsinαandmthatleavem(cid:15)and
literature (namely, M13). The shearing box equations
α(cid:15) invariant leave the results invariant up to numerical
are expressed in terms of length scales proportional to
accuracy of the computational algorithm.
thelocalscaleheightofthedisk(e.g.,Umurhan&Regev,
2004). In the absence of any other superimposed flow
fieldstherearenoothernaturallengthscalesintheprob-
lem. Given the aforementioned flow fields described in
the previous section, there are now two length scales are
introduced to the problem, namely, ∆ (or ∆±) and (cid:15) re-
spectively describing the width(s) of the shear layer (jet
flow) and the length scale of their corresponding transi-
tion zones.
Similarly, the shearing box equations are generally ex-
pressedintimescalesproportionaltotheinverserotation In more realistic disk models, the Brunt-V¨ais¨al¨a fre-
rate of the disk, Ω−1. The plain Coriolis effects (in the quency is a function of vertical position. Indeed, the
absence of any V )0are expressed in equation (36) as the vertical component of gravity has the approximate lo-
square of the epic0yclic frequency ω2, which is in units of cal dependence g ∼ Ω2z. Buoyancy, which is a non-
0
ε dimensional quantity, is best gauged by the vertical en-
Ω2. Intermsoftheshearingboxequationsutilizedinthis
0 √ tropy profile which, in turn, is strongly dependent on
work,theBrunt-V¨ais¨al¨a frequency( gβ)isanindepen- the global disk model under consideration. Provided the
dent time scale appearing in the system. The addition entropystructureisstablystratifiedandsymmetricwith
of the flow fields V , discussed in Section 4, introduces respecttothediskmidplane,wesupposethatthevertical
0
additionaltimescales: 1/ω associatedwiththevorticity gradient of buoyancy β ∼ z/(HH ), in which H is the
a β β
step flow and shear layer models or, the two timescales vertical variation scale of the (non-dimensionalized) en-
1/ω± for the jet flow model. tropywhileH istheusualpressurescale-height. Suppos-
a
Had we chosen to do so, the governing equations ingthemodeldependentH ∼O(H),thenthesquareof
β
(like, for example, equation 36) would have been non- theBrunt-V¨ais¨al¨a frequencytakesontherangeofvalues
dimensionalized according to (cid:15) (spatial scales) and Ω−1 given by the relationship: gβ ≈ Ω2(z/H)2 . Consider-
0 0
(temporal scales). The resulting equations would then ingthattheshearingboxequationsareformallyvalidto
transparently exhibit their dependence upon the non- a few scale heights above the disk midplane, and since
dimensional parameters: α(cid:15), m(cid:15), gβ/Ω2 and ω /Ω . For ourmodelformulationassumesconstantvaluesofgβ,we
thevorticity-stepmodelprofilecalculati0on,thesaefo0urpa- make sure to consider values of 0 < gβ ≤ 6Ω2, in which
0
rameters completely characterize the system’s solutions. the upper bound corresponding to about 2.5 disk scale
The results of shear layer problem depend upon five pa- heights while the lower bound corresponds to locations
rameters, i.e., the four parameters of the vorticity-step in the vicinity of the disk midplane.
8 Umurhan et al.
6. SOLUTIONMETHOD we often use this method to determine the approximate
solutiononacoarsegridandthenwerefinethissameso-
Wesolvethecoupledintegralequations(17–19)assum-
lution through Stage 2: which interpolates the solution
ing normal mode perturbations of the form
ontoafinergrid(either2or4times)andthensolvesthe
(cid:32) ζ(cid:48)(x,t) (cid:33)  ζˆ(x)  discretized system matrix operator system using a stan-
D(cid:48)(x,t) =Dˆ(x)e−iωt+c.c. (41) dard Newton-Raphson-Kantorovich solution technique.
θ(cid:48)(x,t) θˆ(x) The benefit of this approach is that Stage 1 produces
all of the normal modes permitted by the system and
together with the corresponding Green’s function solu- Stage 2 helps to identify which of the numerically de-
tions for the diagnostic relationships φ,ψ(cid:48) and Π(cid:48) found termined solutions are spurious and which are robust.
in equations (21) and (24). As noted earlier, we seek Spurious solutions are identified as those candidate nor-
solutions that show exponential decay as |x| (cid:29) 1 which mal mode solutions whose eigenvalues show no conver-
necessarilyprecludessimplewavenormalmodesolutions gencewithincreasedresolution. Z-modesarefoundtobe
that have no attenuation in the far limit. particularly tricky to obtain reliably, generally requiring
The system is a series of three coupled integral equa- anywhere from 700-1500 grid points of radial resolution
tions. Thevariables,ζ,Dandθ,arediscretizedoneither ondomainsrangingfromπ to2π (alsoseesection7.1.4).
a uniform or Gaussian grid x of N points, but for the
j 7. LINEARTHEORYRESULTS
purposes of this report we quote the results developed
using a uniform grid only. 7 Thus each variable is nu- We have scanned for solutions in the event that the
basic flow is a pure Keplerian velocity field, i.e. for
merically represented as a column vector corresponding
toitsvaluesonthegrid,i.e. ζˆ(cid:55)→ζ ,Dˆ (cid:55)→D andθˆ(cid:55)→θ . V(x) = 0. We find no converged continuous normal
j j j modes. The reason appears to us clear: the structure
The Green’s functions are turned into matrix operators
of the alternative formulation of the equations, i.e. (26–
so that, for instance, the stream function (ψˆ (cid:55)→ ψj) is 28) together with Ξ(cid:48) = 0 indicates that the system sup-
written as a matrix operation relating it to the vertical
ports shear modified inertial-gravity waves. No normal
vorticity, that is to say
modesareexpectedsincetheradialextentofthesystem
dx is infinite and, as such, is unable to support a “global”
ψj =G(jψn)⊗ζn, G(jψn) ≡−2αe−α|xj−xn|, supported mode. Normal modes are potentially possible
only if other boundary conditions are adopted (not done
wheredxisthegridspacingandthesymbol⊗isthema- here). Ofcourse,asaninitialvalueproblemwithagiven
trix multiplication operation. The corresponding deriva- initial disturbance, this system would respond by shed-
tive of dψˆ/dx(cid:55)→(cid:0)dψ(cid:1) is written as ding inertia-gravity waves which are non-normal mode
j solutions of (26–28). These waves would propagate out
(cid:16)dψ(cid:17) =dG(ψ)⊗ζ , to x→±∞, but they donot qualify as normalmodes in
jn n this case.
j
in which Normal modes do exist for V (x) (cid:54)= 0. In this case,
0
there are two kinds of modes supported which we hence-
dx
dG(ψ) ≡ sgn(x −x )e−α|xj−xn|. forth refer to as “Rossby modes”(R-modes) and “zom-
jn 2 j n bie modes” (Z-modes). The R-modes are the three-
Asimilarsetofmatrixoperationsaredefinedandimple- dimensionalcontinuationoftheclassicaltwo-dimensional
mented for the potential function φˆ(cid:55)→ φ and the pres- shearmodesexaminedintheliterature(e.g. theRayleigh
j shear layer, the RWI,the triangular jet, etc.). R-modes
sure field Πˆ (cid:55)→Π . In this construction, the exponential
j areunstablewhentwoormoreRossbywaves(sometimes
decay of solutions as |x|→∞ is ensured. The complete
known as “Rossby edge waves”), each being associated
setofequations(17–19)isconvertedintoasinglematrix
with local extrema in the radial vorticity gradient of the
form
basic shear flow, become resonantly phase-locked due
∂V
=M⊗V (42) to their mutual interaction, i.e., the counter propagat-
∂t ing Rossby wave mechanism (Baines & Mitsudera 1994;
with Heifetz etal. 1999;Umurhan 2010). Thereare nounsta-
bleR-modesinthesinglevorticitystepflowfieldbecause
(cid:16) (cid:17)T
V= ζ ,···,ζ ,···,ζ ,D ,···,D ,θ ,···,θ . it can support only a single Rossby wave precluding the
1 j N 1 N 1 N
possibility of resonant wave-wave interaction.
(43)
Ontheother hand, theZ-modes aredifferentfrom the
Assumption of the normal mode form in equation (41)
R-modes in that instability in these modes involves the
turns the above into a single matrix problem to deter-
resonantinteractionbetweenasingleRossbywaveanda
mine the unknown eigenvalues −iω. M is constructed
buoyant critical layer(s) nearby. In the following subsec-
following the method described in Umurhan (2010).
tionsweexaminethepropertiesofZ-modesforthethree
We then go through two stages to obtain a solution.
model shear flows.
Stage1 usesstandardmatrixinversionmethodstoestab-
lish both eigenvalues and the corresponding eigenfunc- 7.1. Z-modes in the vorticity step model
tions. Becausethemethodiscomputationallyexpensive,
7.1.1. Structure of the basic carrier wave
7 All solutions obtained and reported herein are equivalently Thevelocityshearfieldgiveninequation(38)cansup-
obtainedusingeitherdiscretization. portasinglelocalizedRossbyedgewavethatpropagates
Critical layers and disk turbulence 9
(a) perturbation buoyancy: imΘˆ(x)ei(αx+mz)
3
δx∼0.12
singlejump,open
2 ω =−0.42
a
gβ =1
imΘˆ(x)1 ε =0.02
α =1.2 m =3.3333
0
ω =0.036i + 0.402
−0.4 −0.2 0 0.2 0.4 0.6 0.8 1 1.2
x
(b) perturbation vorticity: ζˆ(x)ei(αx+mz)
6
4
ˆ
ζ(x)
2
0
−0.4 −0.2 0 0.2 0.4 0.6 0.8 1 1.2
x
(c) perturbation dilatation: Dˆ(x)ei(αx+mz)
3 δx∼0.12
real
2
imaginary
1
ˆ
D(x)
0
−1

−0.4 −0.2 0 0.2 0.4 0.6 0.8 1 1.2
x
Figure 2. ArrayofZ-modeeigenmodesforthevorticitystepproblem: gβ=1.0Ω2,ω =−0.42Ω ,α=1.2andm=3.33with(cid:15)=1/30.
0 a 0
ThissampleflowfieldhasanequivalentdeviationflowfieldwithRo=0.21. Thepanelsshowtheperturbationfieldsfor(a)thebuoyancy
(b)thevorticity,(c)thedilatation. Thebasicstatevorticitygradientappearsimprintedintheperturbationvorticityfieldaroundx=0.
The induced critical layers appearing at x± ≈−0.45,0.62 are most prominent in the buoyancy and dilatation fields. Numerical details:
bg
lineargridusedisN =3064withsolutionsshownbyopencircleswithafittedcurveconnectingthem. Solutionscalculatedonthedomain
−π<x<π,howeveronlytheactivepartsareshowninthefigure.
along the azimuthal direction. The core of the wave is after a little manipulation becomes
localized around x = 0 as this is the location where the
radialgradientofthebasicstateshearprofileisthegreat- (cid:18)1(cid:19)(cid:32) ω (cid:33)
x ≈ a . (44)
est. In the limit where m(cid:15) → 0, together with m (cid:28) α, (cid:0) (cid:1)
pp α −3Ω /2+ω sgn x
thisdisturbancecanbethoughtofasaverticallyuniform 0 a pp
azimuthally propagating edge wave and it was demon-
Similarly, the critical layers associ√ated with gravity ef-
stratedinUmurhan(2010)tohavearealfrequencygiven fects, i.e. those points αv (x±)± gβ−Re(ω) = 0 are
by ω =ω in the effective limit α(cid:15)→0 (α fixed, (cid:15)→0). 0 bg
a given by
The frequency response of the unstable modes reported
hereareverynearlyequaltothisvalue,i.e. Re(ω)≈ωa– (cid:18) (cid:19) √ 
andthisisespeciallytrueformodeswithvaluesofm≤2 1 ω ± gβ
(see Figure 5). We consider this Rossby edge wave to be x±bg ≈ α −3Ω /2a+ω sgn(cid:16)x±(cid:17). (45)
the basic carrier mode of the instability associated with 0 a bg
Z-modes.
Theaboveareroughguides–inactualitythepositionof
With this insight, we can predict the approximate lo-
the critical layers will differ from the above approximate
cations of the various critical layers in the limit where
form when m (cid:29) 1, since ω of the Rossby wave is no
both m is relatively small and (cid:15) (cid:28) 1.8 For the classical
longerexpectedtobegivenbyω . Nonetheless,whenthe
shear critical layer we have αv (x )−Re(ω)=0, which a
0 pp correct value of Re(ω) appropriate for the given Rossby
8 Wealsonotethatinthissenseweconsider(cid:15)tobesufficiently waveisinput,weobtainthecorrectcriticallayerposition
small if its value is below any numerically resolvable grid scale as expected.
length. In Figure 2 we show a fairly typical result involving
10 Umurhan et al.
(a) perturbation x-velocity uˆ(x)ei(αx+mz) thestronglocalizednatureofthedisturbanceswherethe
0.05 singlejump,open amplitude of the perturbations decays to zero rapidly as
uˆ(x) 0 (cid:116)a =−0.42 |x|(cid:29)(cid:15) (at least, several 100 times (cid:15)). It is generally for
g(cid:96) =1
−0.05 (cid:161) =0.02 this reason we are confident that these responses (and
−0.1 (cid:95) =1.2 m =3.3333 overall physical effects) are independent of the assumed
(cid:116) =0.036i+0.402
−1.5 −1 −0.5 0 x 0.5 1 1.5 2 radialboundaryconditionssolongastheboundariesare
(b) perturbation y-velocity: vˆ(x)ei(αx+mz) set far enough away. We have scanned the response of
0.05 this system using both impenetrable walls and periodic
vˆ(x) 0 conditions in the radial direction and we find that the
−0.05 response is essentially the same so long as the bound-
−0.1
aries themselves are set far enough apart. In principle
−1.5 −1 −0.5 0 x 0.5 1 1.5 2
(c) perturbation w-velocity: wˆ(x)ei(αx+mz) this means L being sufficiently greater than π, but in
0.6 real practice we find that even slightly larger than L = π is
0.4 imaginary good enough. Observable signatures of the critical lay-
0.2 ers at x = x± are also evident in these fields, however,
wˆ(x) 0 bg
it is most prominent in the vertical velocity. This has
−1.5 −1 −0.5 0 x 0.5 1 1.5 2
(d) perturbation pressure: Πˆ(x)ei(αx+mz) profound implications for nonlinear driving (Section 8).
0.01 The perturbation potential vorticity, is displayed in
Πˆ(x) 0 Figure 4. Consistent with our analysis of Section 3.1,
−0.01
−0.02 rimeaalginary we see that Ξ(cid:48) has power purely in the region of x = 0
−0.03 and no discernible power in the regions around the criti-

−1.5 −1 −0.5 0 x 0.5 1 1.5 2 cal points x±. This is rationalized because the potential
bg
vorticity is (linearly) forced only by the advection of the
Figure 3. Likethepreviousfigurebutshowingthefields(a)u(cid:48), meanvorticitygradientwhichhasappreciablepoweronly
(b) v(cid:48),(b) w(cid:48) and (d) Π(cid:48). All quantities show substantial drop- around x=0. For us, this observation lends further cre-
offasxapproachestheleftandrightboundariesemphasizingthe dence to our claim that the main perturbation structure
strong localized nature of these disturbances. We emphasize that is that of a Rossby wave and that it is this wave that
the solutions are calculated on the domain |x| < π but that for
triggers the critical layers at x±.
purposes of clarity only the inner portions are displayed since all
bg
solutionsshowexponentialdecaywithincreased|x|.
7.1.2. Critical layer scaling
perturbationpotentialvorticity: Ξˆ(x)ei(αx+mz)
6 The sizes of the critical layers scale approximately as
(cid:112)
5 ∼ Im(ω)/Ω /m. This fact can be seen from an ele-
0
mentary boundary layer scaling analysis of Eq. (36) in
4
Ξˆ(x) real the vicinity of x = xbg . [For an analogous critical layer
3 imaginary analysisseethethoroughexaminationofabarotropicjet
instability performed by Balmforth and Piccolo (2001)].
2 The critical layer arises from the first term on the RHS
of equation Eq. (36), and when the layer is activated,
1
this term balances the second x derivative of Πˆ:
−10. 5 −1 −0.5 0 0.5 1 1.5
x
ω2−gβ
∂x2Πˆ ≈ (cid:16)√ ε(cid:17)(cid:16)√ (cid:17)m2Πˆ.
Figure 4. Potentialvorticityfortheproblemrelatingtotheprevi- gβ−αv (x)+ω gβ+αv (x)−ω
oustwofigures. ThepotentialvorticityΞ(cid:48)isbasedonitsdefinition 0 0
inequation(25). Mostpowerinthisquantityiscontainedmainly (46)
in the region surrounding x = 0 corresponding to the maximum in which we have explicitly factored the denominator of
ofthemeanvorticitygradient. Thereisnodiscerniblepowernear
thetermontheRHSoftheaboveexpression. Foragiven
thecriticalpointsx± -asexpected.
bg solution value ω, at the critical layer we have according
to equation (37),
unstable modes of the system. The three panels show
the quantities θ(cid:48),ζ(cid:48) and D(cid:48). The vorticity shows strong (cid:112) (cid:104) (cid:105)
power near x = 0. This arises from the radial trans- gβ± αv (x )−Re(ω) =0,
0 bg
port of the mean vorticity gradient term u(cid:48) · v and
0xx
indicates that the perturbation is structurally that of a A first order Taylor series expansion of the term αv (x)
Rossby wave centered at that position. According to near x=x produces 0
equation (45) and given the parameters of the system, bg
thepositionsofthetwobuoyantgravitywavecriticallay- (cid:12)
ers are x− ≈ −0.448 and x+ ≈ 0.616. The importance αv (x)≈αv (x )+αdv0(cid:12)(cid:12) (cid:0)x−x (cid:1)
bg bg 0 0 bg dx(cid:12) bg
of these locations are self-evident in both the buoyancy x
bg
and dilatation fields, demonstrating amplified power in
the vicinity of those critical layers. Considering the assumption |Im(ω)| (cid:28) |Re(ω)|| (noting
The corresponding perturbation velocity fields and that all of the growth rates determined in this study are
pressure fluctuations are shown in Figure 3, indicating usuallyafactorof10smallerthanΩ ),puttingallofthe
0