Astronomy&Astrophysicsmanuscriptno.GBMR2015-p (cid:13)cESO2016 January19,2016 Tidal inertial waves in the differentially rotating convective envelopes of low-mass stars I. Free oscillation modes M.Guenel1,C.Baruteau2,S.Mathis1,3,andM.Rieutord2 1 LaboratoireAIMParis-Saclay,CEA/DSM-CNRS-UniversitéParisDiderot,IRFU/SApCentredeSaclay,F-91191Gif-sur-Yvette Cedex,France 2 InstitutdeRechercheenAstrophysiqueetPlanétologie,ObservatoireMidi-Pyrénées,UniversitédeToulouse,14avenueEdouard Belin,31400Toulouse,France 6 3 LESIA,ObservatoiredeParis,PSLResearchUniversity,CNRS,SorbonneUniversités,UPMCUniv.Paris6,Univ.Paris-Diderot, 1 SorbonneParisCité,5placeJulesJanssen,92195Meudon,France 0 e-mail:[mathieu.guenel,stephane.mathis]@cea.fr;[clement.baruteau,michel.rieutord]@irap.omp.eu 2 Received.../accepted... n a J ABSTRACT 8 1 Context.Star-planettidalinteractionsmayresultintheexcitationofinertialwavesintheconvectiveregionofstars.Inlow-massstars, theirdissipationplaysaprominentroleinthelong-termorbitalevolutionofshort-periodplanets.Turbulentconvectioncansustain ] differentialrotationintheirenvelope,withanequatorialacceleration(asintheSun)ordeceleration,whichcanmodifythewaves’ R propagationproperties. S Aims.Weexploreinthisfirstpaperthegeneralpropagationpropertiesoffreelinearinertialwavesinadifferentiallyrotatinghomoge- . neousfluidinsideasphericalshell.Weassumethattheangularvelocitybackgroundflowdependsonthelatitudinalcoordinateonly, h closetowhatisexpectedintheexternalconvectiveenvelopeoflow-massstars. p Methods.Weusei)ananalyticalapproachintheinviscidcasetogetthedispersionrelation,fromwhichwecomputethecharacteristic o- trajectoriesalongwhichenergypropagates.Thisallowsustostudytheexistenceofattractorcyclesandinferthedifferentfamiliesof r inertialmodes;ii)high-resolutionnumericalcalculationsbasedonaspectralmethodfortheviscousproblem. t Results.Wefindthatmodesthatpropagateinthewholeshell(Dmodes)behavethesamewayaswithsolid-bodyrotation.However, s a anotherfamilyofinertialmodesexists(DTmodes),whichcanpropagateonlyinarestrictedpartoftheconvectivezone.Ourstudy [ shows that they are less common than D modes and that the characteristic rays and shear layers often focus towards a wedge – orpoint-likeattractor.Moreimportantly,wefindthatfornon-axisymmetricoscillationmodes,shearlayersmaycrossacorotation 1 resonancewithalocalaccumulationofkineticenergy.TheirdampingratescalesverydifferentlyfromwhatweobtainforstandardD v modesandweshowanexamplewhereitisindependentofviscosity(Ekmannumber)intheastrophysicalregimeinwhichitissmall. 7 1 Keywords. hydrodynamics–waves–planet-starinteractions–stars:rotation 6 4 01. Introduction Verbunt&Phinney(1995)andmorerecently,Meibom&Math- . 1 ieu (2005) who showed that the oldest systems have circular- 0The tidal interaction between a star and an orbiting companion izedorbits,thoughthespin-orbitsynchronizationisstillunclear 6results from the differential force exerted by each body on the (Meibometal.2006).Observationsofextrasolarplanetsbythe 1other.Indeed,thegravitationalpotentialcreatedbythecompan- radial-velocityandtransitsmethodshavedevelopedrapidlyover v:ion is not uniform inside the star since it is an extended body the past decade and stimulated interest in looking for signa- iand not a point-like mass. If the orbit is close enough, the star tures of tidal interactions in star-planet systems : for instance, Xexperiencesanon-negligibleperiodicalforcingwhichgenerates Pont (2009) looked for an excess rotation in stars with plan- rtidal flows. Their dissipation into heat through various friction ets (due to substantial inward migration) while Husnoo et al. a mechanisms (see Zahn 2013; Mathis & Remus 2013; Ogilvie (2012) used observed eccentricities to conclude that tidal inter- 2014)takesenergyawayfromthesystemandallowsforaredis- actions play a prominent role in the orbital evolution and sur- tributionofangularmomentum,whichleadstotheevolutionof vival of hot Jupiters (see also Lai 2012; Guillot et al. 2014). the spin and orbital parameters on secular timescales. Depend- Moreover, Jackson et al. (2008) consistently checked that the ingontheinitialdistributionofangularmomentumbetweenthe observedloweccentricitiesofclose-inplanetscanbeexplained individualspinsandtheorbit,thesystemmayevolvetowardsan by tidal evolution processes, and calibrated the required tidal equilibriumstate—wheretheorbitiscircularandthespinsof quality factors. Hansen (2010) and Penev et al. (2012) used thebodiesaresynchronizedandalignedwiththeorbit—orlead observed populations of exoplanets and tidal evolution models totheinspiralorejectionofthecompanion(Hut1980). to constrain tidal dissipation in stars, and both studies seem to agree that their results are inconsistent with the dissipation in- Observational evidence for ongoing tidal interactions was ferredfromthecircularizationofbinarystars,suggestingthata foundinclosebinarystars,forinstancebyGiuricinetal.(1984), Articlenumber,page1of17 A&Aproofs:manuscriptno.GBMR2015-p different mechanism is at play in this case. Winn et al. (2010) the effects of resonant excitation of gravito-inertial modes and andAlbrechtetal.(2012)alsoproposedthattidaldissipationin possibleresonancelocking,whichcanspeedupsignificantlythe convective regions could be responsible for the low obliquities secularevolution. of hot Jupiters orbiting around cool stars, though Mazeh et al. Unlikestablystratifiedradiativezones,stellarconvectivere- (2015)foundthatthispropertyextendstoKeplerObjectsofIn- gionsareveryslightlyunstablystratifiedandgravitymodescan- terestwithorbitalperiodsuptoatleast50days,forwhichtidal not propagate inside them. For more than thirty years, the pos- interactions should be negligible. These results have to be put sible effects of a dynamical tide in convective zones have been into perspective since many other processes are responsible for overlookedinfavourofZahn’sequilibriumtidetheory.However, shapingthearchitectureofthesystems—suchasmigrationin asimpleapproachtomodelthisdynamicaltideistosimplyig- theprotoplanetarydisk(e.g.Baruteauetal.2014),planet-planet noretheconvectivemotionsandtoconsiderlineardisturbances scattering (e.g. Chatterjee et al. 2008), Kozai oscillations (e.g. toaneutrallystratifiedfluid.Inthiscontext,thedynamicaltide Naoz et al. 2011), star-planet magnetic interactions (see Stru- consistsoflow-frequencyoscillationsprimarilyrestoredbythe garek et al. 2014), magnetic spin-down of the star (see Barker Coriolis acceleration — otherwise known as inertial waves — &Ogilvie2009;Damiani&Lanza2015),etc.—andtheglobal that propagate in a sphere (if the body is fully convective) or a pictureremainsquiteuncertain. sphericalshell(forsun-likestarswithanouterconvectiveenve- Onthetheoreticalside,thetidalresponseofastarhasbeen lope). These waves have been studied long ago in a uniformly investigated in the past, starting with the theory of the equilib- rotating,incompressibleandinviscidfluidbyThomson(1880); rium tide, which is the flow induced by the quasi-hydrostatic Poincaré (1885); Bryan (1889); Cartan (1922) (also Greenspan adjustment to the tidal potential, and its dissipation by turbu- (1968))andtheyexhibitremarkableproperties:theirfrequency lentfrictionintheconvectiveenvelopeoflow-massstars(Zahn ω˜ in the fluid frame (rotating with angular velocity Ω) is re- 1966a,b). This work was followed by the study of the radiative strictedto[−2Ω,2Ω]whiletheirspatialstructureisgovernedby damping of tidal gravity modes (Zahn 1975) in more massive ahyperbolicsecond-orderpartialdifferentialequationnamedaf- starswithanouterradiativeregion.Zahn(1977)thenestimated terPoincaré.Theypropagatealongraysthatareinclinedwitha thestrengthsofthesemechanismsinordertoderivethecharac- constant angle arcsin(ω˜/2Ω). Such hyperbolic equations along teristictimescalesforcircularization,synchronizationandalign- with boundary conditions in a closed container generally yield ment for binary stars. Since then, most works about tidal inter- an ill-posed problem whose solutions may be singular. The be- actions were based on Zahn’s approach and focused on giving haviorofinertialwavesisdrivenbytheraysproperties:inafull morerealisticdescriptionsofthemechanismsthoughttobere- sphere,theyfillthewholedomainforallfrequenciesandasetof sponsiblefortidaldissipation—especiallythedynamicaltidein globalnormalmodeswithfrequenciesthatisdensein[−2Ω,2Ω] radiativeregionsinearly-typestars.Forinstance,Savonije&Pa- hasbeenfoundinthecaseofafullsphere(Bryan1889).Thisis paloizou(1983,1984)performednumericalsimulationsofgrav- not true in the case of a spherical shell (Rieutord & Valdettaro itymodesexcitedbyaperiodictidalpotentialinanon-rotating 1997),forwhichraysoftenfocustowardslimitcyclescalledat- massive star and derived the associated timescales for the or- tractors(seeMaas&Lam1995;Rieutordetal.2001)thatexist bital evolution of massive binaries. The effects of rotation on in narrow frequency bands and depend on the size of the inner thesegravitymodesandtheassociatedtidaldissipation,ignored core.Theinclusionofviscosity—asasimplificationofturbu- byZahn,weretheninvestigatedbyRocca(1987,1989)whoin- lence—givesrisetopeculiarmodesthatpossessshearlayerslo- cludedtheeffectsoftheCoriolisacceleration,andGoldreich& calizedaroundtheaforementionedattractors(e.g.Rieutordetal. Nicholson(1989)whotookintoaccountpossibleinternaldiffer- 2001). entialrotation. The tidal force exerted by a close planetary or stellar com- After the discovery of the first exoplanetary systems in the paniononitshoststarmaythereforeexciteinertialwavesinthe mid1990’s,thestudyofstar-planettidalinteractionregainedin- external convective envelope of low-mass stars, if one or more terestandstudiesshiftedtowardslower-massstarswithradiative components of the tidal potential has a frequency that fall into coresandconvectiveenvelopes,forwhichtheacceptedtidaldis- the inertial range. The energy dissipated and the angular mo- sipationmechanismsofarwastheviscousdissipationofZahn’s mentum carried by these waves may play an important role on equilibrium tide in the convective zone. Terquem et al. (1998), theevolutionoftheorbitalarchitectureofclose-inplanetarysys- Goodman & Dickson (1998) and Witte & Savonije (2002) per- temsandontherotationoftheircomponents,yetithasonlyre- formed numerical calculations of the dynamical tide in the ra- centlystartedtobeinvestigated(seeOgilvie&Lin2004,2007; diativecoreofnon-rotatingsolar-likestarsandstudiedtheeffect Ogilvie 2009; Goodman & Lackner 2009; Rieutord & Valdet- ofresonancesinthecoreofasolar-typestar,andfoundthatthe taro 2010; Le Bars et al. 2015) in uniformly-rotating fluids. associatedtidaldissipationcouldbemoreefficientthanthevis- Thetidaldissipationinducedbyforcedinertialmodeshasbeen cous dissipation of the equilibrium tide in the convective enve- found to be a very erratic function of the excitation frequency lope.Barker&Ogilvie(2010,2011);Barker(2011)alsostudied andvariesoverseveralordersofmagnitude,thoughfrequency- thepossiblenon-linearbreakingoftidalgravitywavesnearthe averaged estimates can be obtained analytically (Ogilvie 2013) centre of a solar-type star and its influence on tidal dissipation and evaluated along stellar evolution (Mathis 2015), showing and on the survival of hot Jupiters and angular momentum de- strongvariationswithstellarmass,ageandrotation.Ithasbeen position. Progressively, the low-frequency regime — where ro- shown by Baruteau & Rieutord (2013) that differential rotation tation is likely to have the most important effects — of stellar may strongly affect the propagation properties of linear inertial oscillations regained interest with the works by Savonije et al. modesofoscillations.Theirstudywasrestrictedtoshellularand (1995); Savonije & Papaloizou (1997); Papaloizou & Savonije cylindricalrotationprofiles,butturbulentconvectioncanalsoes- (1997) who simulated the full tidal response of a uniformly- tablishvariousdifferentialrotationprofiles(seeBrun&Toomre rotatingmassivestarandfoundsignaturesoftidally-excitedin- 2002;Brownetal.2008;Mattetal.2011;Gastineetal.2014), ertial waves in the convective core, as well as gravito-inertial including conical rotation as observed in the Sun (Schou et al. wavesintheradiativeenvelope.Witte&Savonije(1999,2001) 1998;Garcíaetal.2007).Inthiswork,weexplorethepropaga- studied the tidal evolution of massive binary systems including tion and dissipation properties of inertial modes in stellar con- Articlenumber,page2of17 Gueneletal.:Tidalinertialwavesinthedifferentiallyrotatingconvectiveenvelopesoflow-massstars-I vective envelopes (for low-mass stars only) with conical differ- perioddetermineswhatwecallthe“Rossbyscale"belowwhich entialrotation,namelydifferentialrotationthatonlydependson turbulenceislittleinfluencedbyrotationandisfastcomparedto latitude. thewaveperiod.Ourmodelassumesthateddiessmallerthanthe In Sect. 2, we expose our physical model and the associ- Rossbyscalecanberepresentedbyaturbulentviscosity.Tofix atedhypothesesaswellasthedifferentialrotationprofileweuse. ideas,letusconsidertheexampleoftheSun.Theboundingfre- TheninSect.3,wecarryoutananalyticalanalysiswhichyields quency2Ω correspondstowaveswithaperiodof12.5days. max thegeneralpropagationpropertiesofinertialwavesintheinvis- FromMixing-LengthTheory,whichsaysthatthetypicalscaleof cid case. These results are compared to the viscous problem in turbulenceisthemixing-length(abouttwicethepressurescale Sect.4wherewecomputenumericallythevelocityfieldsofsuch height), we note that the largest eddies have a typical scale of inertial modes, which sometimes have very different properties 100Mm,withatypicalvelocityof50m/s.Thiswouldbethetyp- thaninthesolid-bodyrotationcase.Finally,weidentifyinSect. ical numbers for the so-called giant cells characterized by a 20 5 which of these properties are important for the forced prob- daysturnovertimescale(Mieschetal.2008).Withthesenum- lemthatwewilltreatinasubsequentpaper,andinparticularthe bersweeasilyfindthattheRossbyscaleis50Mm.Thusallthe evaluationofthetidaldissipationinducedbythesewaves. eddies of scale much smaller than 50Mm are assumed to influ- encetheinertialmodesthroughaturbulentdiffusion.Theinflu- enceofturbulentmotionswithlargerscalesisadmittedlymuch 2. Inertialmodesindifferentiallyrotating moredifficulttotakeintoaccount.However,wemaynotethatin convectiveenvelopes thecaseoftheSunthedifferentialrotationflowismuchstronger thanthelargescaleeddies(1km/scomparedto50m/s).Thus,as 2.1. Physicalmodel afirststep,neglectingthedirectinteractionoftheseeddieswith Sincewewanttostudythepropagationofinertialwavesinthe theinertialwavesisconsistentasfarasordersofmagnitudeare differentially rotating convective envelopes of low-mass-stars, concerned. wewillconsiderasimplifiedsetupwithanincompressible,vis- cousfluidinsideasphericalshellofexternalradiusRandinter- 2.2. Mathematicalformulation nalradiusηR(0<η<1)thataccountsfortheboundarybetween theradiativecoreandtheconvectiveenvelope. Inaninertialframewiththeusualsphericalcoordinates(r,θ,ϕ), AsmotivatedinSect.1,weuseaconicaldifferentialrotation andafterlinearizationoftheNavier-Stokesequationaroundthe profile—dependingonlyonthecolatitudeθ—whichreads steady state (described by the “0” subscripts), we obtain (e.g. Baruteau&Rieutord2013) Ω(θ)=Ω0(θ)/Ωref =1+εsin2θ, (1) ∂u1+Ω ∂u1+2Ω e ×u +rsinθ(u ·∇Ω )e =−∇p +ν∇2u , ∂t 0 ∂ϕ 0 z 1 1 0 ϕ 1 1 so that the dimensionless angular velocity of the background flow is 1 at the poles and 1+ε at the equator. The quantity ε (2) isaparameterthatdescribesthebehaviorofthedifferentialro- where ∂/∂t denotes the partial time-derivative, Ω is the back- tation: 0 ground angular velocity of the fluid, e = cosθe − sinθe is z r θ – ε > 0 is for solar differential rotation (equatorial accelera- the unit vector along the rotation axis, u1 is the wave velocity tion), fieldand p1 denotesthereducedpressureperturbation,whichis – ε < 0isforanti-solardifferentialrotation(equatorialdecel- the pressure perturbation divided by the reference density (ρ0). eration). Notethatonlythelasttermoftheleft-handsideisproportional totheshearofthedifferentialrotation,whiletheothertermsare The value (cid:15) ≈ 0.3 approximates the differential rotation profile formallyidenticaltothecaseofsolid-bodyrotation. intheSun’sconvectiveenvelope. Along with this equation, we use the linearized continuity We point out that the above rotation profile does not sat- equation isfy the Taylor-Proudman theorem, which states that the veloc- ∇·u =0 (3) ityfieldofarotatingincompressiblehomogeneousinviscidfluid 1 mustbeconstantalongcolumnsparalleltotherotationaxis.As sincethefluidisalsoassumedtobeincompressible. is now well-known (e.g. Brun & Toomre 2002; Brown et al. We look for velocity (u ) and reduced pressure (p ) per- 2008),thisdifferentialrotationprofileissustainedbyReynolds turbations of angular frequen1cy Ω (in the inertial fram1e) and p stressesduetoconvectiveturbulentmotions.Thismeanflowis azimuthal wavenumber m. This means they are proportional to assumed to be dynamically stable. Our analysis will show that exp(cid:16)iΩ t+imϕ(cid:17).Eqs.(2)and(3)canberecastas somemodesareunstableforsomeparametersofthemodel(see p Eq.(19)).Theseinstabilities,interestingperse,shouldbeinter- (cid:40) pretedasgivingthelimitsoftheparametersrangeofthismodel, iΩ˜pu1+2Ω0ez×u1+rsinθ(u1·∇Ω0)eϕ =−∇p1+ν∇2u1, especiallyinitsfutureusefortidallyforcedflows. ∇·u =0, 1 Thewavesthatwestudyinthefollowingarelow-frequency (4) waves: their frequency is less than 2Ω , where Ω is the max max maximum angular velocity of the fluid. Such waves propagate where Ω˜ = Ω +mΩ is the Doppler-shifted wave frequency p p 0 onaturbulent backgroundandpresumablyinteractwith eddies (that is the local wave frequency in a frame rotating with the that have similar time scale (just like the stochastically excited fluid). acoustic modes of solar-type stars (Zahn 1966b; Goldreich & WenormalizeallfrequenciesbyΩ —whichwedefineas ref Keeley 1977)). In a Kolmogorov type turbulence, the turn-over the background angular velocity at the poles i.e. Ω (θ = 0) — 0 timescaleofeddiesdecreasesas(cid:96)2/3 if(cid:96) isthescaleoftheed- and distances by R. We introduce the associated dimensionless dies. The equality between the turnover time scale and rotation quantities:theradius x = r/R,thevelocityfieldu = u /RΩ , 1 ref Articlenumber,page3of17 A&Aproofs:manuscriptno.GBMR2015-p thereducedpressurep= p /R2Ω2 ,thedifferentialrotationpro- formercase,theequationgoverningthepathsofcharacteristics 1 ref file Ω = Ω /Ω , the wave frequency ω = Ω /Ω in the in- inameridionalplanereads 0 ref p p ref ertial frame (resp. ω˜p = Ω˜p/Ωref in the fluid’s frame) and the dz 1 (cid:16) (cid:17) EkmannumberE =ν/R2Ωref.Fromnowon,weshalluseexclu- ds = 2ω˜2 Az±ξ1/2 , (9) sively the dimensionless quantities defined above. This finally p yieldsthedimensionlesssystemthatweshallsolvenumerically which shows that, contrary to the case of solid-body rotation inSect.4: (Az = 0,ξconstant),differentialrotationtendstocurvepathsof characteristicssincetheright-handsideofEq.(9)nowdepends (cid:40) iω˜pu+2Ωez×u+xsinθ(u·∇Ω)eϕ =−∇p+E∇2u, (5) both on s and z. On the other hand, no characteristics exist in ∇·u=0. theellipticdomain,thereforetherelationξ = 0defines“turning surfaces”onwhichcharacteristicsreflect(seeFriedlander1982). Note that in the numerical simulations presented in Sect. 4 Notethatiftherotationprofileissymmetricbyz→−z,thenEq. below,weactuallytakethecurlofthefirstequationinorderto (7)isalsosymmetricbyz → −z,meaningonlypositivevalues getridofthe∇pterm: of z can be considered. We detail below the case of the conical rotationprofiledefinedinEq.(1). (cid:16) (cid:17) ∇× iω˜pu+2Ωez×u+xsinθ(u·∇Ω)eϕ = E∇×∇2u, ∇·u=0. 3.2. Pathsofcharacteristics,turningsurfacesandclassesof inertialmodes (6) 3.2.1. Pathsofcharacteristics In addition to Eqs. (6), we use stress-free boundary condi- Fortherotationprofilethatweadopt,giveninEq.(1),wehave tions (u·e = 0 and e ×[σ]e = 0, where [σ] is the viscous r r r (cid:16) (cid:17) stress tensor) at the inner and outer boundaries of the spherical A =4Ω(θ) Ω(θ)+εsin2θcos2θ (10) s shell. and In the following section, we define several useful quanti- ties for the description of the propagation properties of inertial A =−4Ω(θ)εcosθsin3θ, (11) z waves, and Eq. (1) implies that they only depend on the colati- sothat tudeθandareindependentoftheaspectratioηofthespherical shell. ξ =16Ω2(θ)ε2cos2θsin6θ (cid:16) (cid:17) +4ω˜2 4Ω2(θ)+4Ω(θ)εcos2θsin2θ−ω˜2 . (12) p p 3. Inviscidanalysis:pathsofcharacteristics, We assume that ε > −1 so that the Rayleigh stability criterion existenceofturningsurfacesandcorotation (A > 0 everywhere in the shell) is always satisfied. From Eq. s resonances (9), the slope of the paths of characteristics in the hyperbolic As a first step, we study the solutions to Eqs. (5) in the invis- domainreads tceirdisltiimcsito(fEin=er0ti)a.lItwaalvloews,swuhsictohsatruedyatvheerydyunsaemfuilcstoooflcfohraruanc-- ddzs =−A(θ)cosθsin3θ+4ωξ˜41/2 (13) derstanding the solutions to the full viscous problem for small p viscosities. where ξ 4Ω2(θ)−ω˜2 3.1. ThegeneralPoincaré-likeequation = A2(θ)cos2θsin6θ+2A(θ)cos2θsin2θ+ p, 4ω˜4 ω˜2 p p FollowingBaruteau&Rieutord(2013),weadoptcylindricalco- (14) ordinates (s,ϕ,z) and we eliminate the components of the ve- locity perturbations, which yields a partial differential equation and on p only. Focusing on the second-order terms, we obtain the Ω(θ) A(θ)=2ε . (15) followingmixed-typepartialdifferentialequation: ω˜2 p ∂∂2sp2 + ωA˜2z ∂∂s2∂pz +1− ωA˜2s∂∂2zp2 +···=0, (7) Achsaeraxcpteecritsetdic,saldletpheenadboonveθqounalnyt.itiesaswellasthedynamicsof p p wherezerothandfirstordertermshavebeenomitted–thatcor- 3.2.2. Turningsurfaces respondstotheWKBJapproximation.InEq.(7), Inthesetupdescribedin3.2.1,ξ issymmetricbytheequatorial 2Ω∂(s2Ω) 2Ω∂(s2Ω) planeoftheshellsoonlyvaluesθ ∈ [0,π/2]needtobeconsid- A (s,z)= and A (s,z)= . (8) s s ∂s z s ∂z ered.Asexplainedinthepreviousparagraph,equationξ(θ) = 0 definesturningsurfacesthatseparatethesphericalshellintoan In the case of solid-body rotation (ε = 0), Eq. (7) reduces hyperbolicdomain(ξ >0)andanellipticdomain(ξ <0).Turn- to the well-known Poincaré equation for inertial waves (Cartan ingsurfacesareportionsofconesofaxiszandaperture2θ.How- 1922;Greenspan1968).Itishyperbolicinthedomainwherethe (cid:16) (cid:17) ever, their location cannot be determined analytically because discriminantξ(s,z) = A2z +4ω˜2p As−ω˜2p ispositive(whichwe ξ(θ)=0isapolynomialequationofdegree6inX =sin2θ.Still will refer to as the “hyperbolic domain”) and becomes elliptic theirexistenceleadstotwoclassesofinertialmodeswithconical inthedomainwhereξisnegative(the“ellipticdomain”).Inthe rotation: Articlenumber,page4of17 Gueneletal.:Tidalinertialwavesinthedifferentiallyrotatingconvectiveenvelopesoflow-massstars-I 1.0 1.0 0.5 0.5 (cid:182) 0.0 (cid:182) 0.0 (cid:45)0.5 (cid:45)0.5 m(cid:61)0 m(cid:61)1 (cid:45)1.0 (cid:45)1.0 (cid:45)6 (cid:45)4 (cid:45)2 0 2 4 6 (cid:45)6 (cid:45)4 (cid:45)2 0 2 4 Ω Ω p p 1.0 1.0 0.5 0.5 (cid:182) 0.0 (cid:182) 0.0 (cid:45)0.5 (cid:45)0.5 m(cid:61)2 m(cid:61)4 (cid:45)1.0 (cid:45)1.0 (cid:45)8 (cid:45)6 (cid:45)4 (cid:45)2 0 2 4 (cid:45)10 (cid:45)8 (cid:45)6 (cid:45)4 (cid:45)2 0 2 Ω Ω p p Fig.1.IllustrationofthetwokindsofinertialmodespropagatinginafluidwithconicalrotationprofileΩ(θ) = Ω (θ)/Ω = 1+εsin2θ,for 0 ref m={0,1,2,4}fromtop-lefttobottom-right.Theeigenfrequencyintheinertialframeω isshowninx-axisandthedifferentialrotationparameter p εisiny-axis.Blue:Dmodesthatexistinthewholeshell(ξ >0everywhere).White:DTmodesthatexhibitatleastoneturningsurfaceinside theshell(thesignofξ changesintheshell).Red:Noinertialmodesmayexistintheshell(ξ < 0everywhere).Orange:Themodesinthe regionbetweenthetwoorangelinesfeatureacorotationlayer(ω˜ =0)insidetheshell.Themeaningoftheblacksolid,dottedanddashedlinesis p explainedinparagraph3.2.2. 1. modes for which at least one turning surface exists within theshell,referredtoasDTmodesfollowingtheterminology introducedinBaruteau&Rieutord(2013)(Dfordifferential ωp =−m±2, (16) rotation,Tforturningsurface), andisshownbytheverticaldottedlinesinFig.1, 2. modesthatpropagateinthewholeshell,whichisentirelyhy- perbolicbecausethereisnoturningsurfacewithinit,named – ξ(θ=π/2)=0,whentheturningsurfacereachestheequator, Dmodes. whichyields The occurrence of D and DT modes is depicted in Fig. 1 ωp =(−m±2)(1+ε), (17) fortheaxisymmetriccase(m = 0)andafewnon-axisymmetric andisshownbythesolidblacklinesinFig.1. cases(m>0).WhiteareasrepresentDTmodeswhichhaveξ < 0 at least once within the shell, while blue regions represent D We point out that there may exist two turning surfaces in a modes for which ξ > 0 everywhere within the shell. The red meridional plane for a given set of parameters. For instance, in areasillustratethecasewhereξ <0everywhereintheshell,for the bottom-left panel of Fig. 1 (where m = 2), the transition whichtheshellhostsnoinertialmodes.Thecurvesthatseparate fromthecentralblueregiontothewhiteregionaroundω = 0 p DandDTmodessatisfy: showsthattwoturningsurfacesoccuratθ=π/4beforesplitting towardsthepoleandtheequator.Ourinvestigationsshowedthat – ξ(θ = 0) = 0,whichcorrespondstothecasewheretheturn- thisonlyoccursform = ±2aroundω = 0,whichwechecked p ingsurfaceisontherotationaxis,andwhichcanberecastas bysolvingnumericallytheequationξ(θ=π/4)=0foreverym, Articlenumber,page5of17 A&Aproofs:manuscriptno.GBMR2015-p whosesolutionsaredepictedbytheblackdashedlinesinFig.1. (23) Anexampleisgiveninthefollowingsubsection(seethebottom- rightpanelofFig.5). AssumingthefluidisstableagainsttheGSFinstability,wede- We point out that for D modes — for which ∀θ, ξ(θ) > 0 fine — we have |ω˜p(θ)| < 2Ω(θ), which means that for D modes, (cid:34)(cid:32) εsin2θcos2θ(cid:33) k εcosθsin3θ(cid:35)1/2 the standard criterion for propagation of inertial waves is met B= 1+ + s , (24) everywherelocally. 1+εsin2θ kz 1+εsin2θ sothatEq.(23)becomes 3.2.3. Criticallayers |k | For non-axisymmetric modes (m (cid:44) 0), ω˜p may vanish in the ω˜p =±2Ω(θ)B(cid:107)kz(cid:107). (25) shell, which corresponds to so-called corotation resonances or critical layers (Baruteau & Rieutord 2013). Since ω˜ is a func- Forsolid-bodyrotation,B = 1andEq.(25)istheusualdisper- p tionofθonly,acriticallayeristheintersectionofaconeofaxis sionrelationforinertialwaves.Consequently,thephasevelocity zwiththesphericalshell.AscanbeseenfromEq.(13),pathsof intheframerotatingwiththefluidis characteristicsmaybecomelocallyverticalatcorotationlayers. k Theirlocationisgivenby: v =±2Ω(θ)B z k, (26) p (cid:107)k(cid:107)3 (cid:114) (cid:18) ω (cid:19) sinθ= ε−1 −1− p (18) whilethegroupvelocityisgivenby m and they exist only when ωp ∈ [−m,−m(1+ε)] (if mε < 0) or v =±2Ω(θ) ks (−k e +k e ) ω ∈ [−m(1+ε),−m] (if mε > 0). Critical layers exist for the g (cid:107)k(cid:107)3 z s s z p D and DT modes located between the two orange lines in the (cid:34) 2Ω(θ)εcosθsin3θ (cid:35) panelsofFig.1. × B− (cid:107)k(cid:107)2B−1 . (27) k k s z FromEq.(25),weseethatacriticallayermayformintheshell 3.3. Dispersionrelation,phaseandgroupvelocities (ω˜ =0)inthefollowingformalcases: p Weexpandinthisparagraphthestudyofthepropagationofin- ertial waves in the inviscid limit by considering the wave dis- – |ks| → ∞atfinitekz,sothatvp = 0andvg = 0atcorotation persionrelation.Thelatterisobtainedbyassumingsolutionsto :inertialwavescannotcrosscriticallayers. Eq.(7)intheform p∝exp[i(kss+kzz)]andadoptingtheshort- – |(cid:12)kz| →(cid:12) 0 at finite ks, so that vp = 0, vg · es = 0 and wavelengthapproximation: (cid:12)(cid:12)vg·ez(cid:12)(cid:12) → ∞, indicating that inertial waves may propagate acrossacriticallayerwithlocallyverticalpathsofcharacter- (cid:34) (cid:35) ω˜2 = kz2 A − ksA , (19) istics. (cid:12) (cid:12) p (cid:107)k(cid:107)2 s kz z – B → 0(cid:12), whic(cid:12)h implies again that vp = 0 while (cid:12)(cid:12)vg·es(cid:12)(cid:12) → (cid:113) ∞ and (cid:12)(cid:12)vg·ez(cid:12)(cid:12) → ∞, which means that in this case as well where (cid:107)k(cid:107) = k2+k2. Assuming the bracket-term is positive, inertialwavesmaypropagateacrossthecorotationresonance s z which is equivalent to assuming that the background rotation withafiniteslopeforthepathsofcharacteristics,sinceB=0 profileisdynamicallystabletoinertialwaves(seethediscussion implies inSec.2.1),wecandefine: (cid:16) (cid:17) dz k 1+εsin2θ 1+cos2θ (cid:34) k (cid:35)1/2 =− s =− . (28) B˜ = A − sA (20) ds kz εcosθsin3θ s k z z Thisresultshowsthatwhenitcomestocriticallayers,con- sothatthewavedispersionbecomesω˜ =±B˜|k |/(cid:107)k(cid:107).Thephase ical and shellular rotation profiles behave similarly — in the p z and group velocities in a frame rotating with the fluid are then sensethatthethreeaforementionedcasesarepossible—while givenby: cylindrical rotation stands out as the case where characteristics of inertial waves can formally never cross critical layers (see k v =±B˜ z k (21) Baruteau&Rieutord2013).Itisimportanttokeepinmindthat p (cid:107)k(cid:107)3 this analysis only reflects the hyperbolic part of the problem while the solutions to the full viscous problem (Sect. 4) may and behaveabitdifferentlythanthecharacteristicsofthepurelyin- (cid:34) (cid:35) k A (cid:107)k(cid:107)2 viscidproblem(asshownforinstancebyFig.13ofBaruteau& vg =±(cid:107)k(cid:107)s3 (−kzes+ksez) B˜+ 2z kskzB˜−1 , (22) Rieutord2013). whichcorrespondstoEqs.(A9)and(A10)inBaruteau&Rieu- tord(2013).Notethatk·v =0asforpureinertialwavesprop- 3.4. Dynamicsofthecharacteristicsofinertialwaves: g attractorcycles,focusingpointsandLyapunov agatinginasolid-bodyrotatingfluid. exponents Applying these formulae to our conical rotation profile Ω(θ)=1+εsin2θ,thedispersionrelationisgivenby Inthecaseofsolid-bodyrotation,thepathsofcharacteristicsina meridionalplanefollowstraightlines,whichformaconstantan- (cid:34)(cid:32) (cid:33) (cid:35) k2 εsin2θcos2θ k εcosθsin3θ glewiththerotationaxis—i.e.thez-axis—,thatonlydepends ω˜2p =4Ω2(θ)(cid:107)kz(cid:107)2 1+ 1+εsin2θ + kzs 1+εsin2θ . on the (uniform) Doppler-shifted wavefrequency ω˜p = ωp +m Articlenumber,page6of17 Gueneletal.:Tidalinertialwavesinthedifferentiallyrotatingconvectiveenvelopesoflow-massstars-I 1.0 ε= 0.00 ω=-3.7500 p m= 2 η= 0.71 0.8 0.6 z 0.4 0.2 0.0 0.0 0.2 0.4 0.6 0.8 1.0 s Fig.2.Pathsofcharacteristicsforsolid-bodyrotation(ε=0).Left:Pathsofcharacteristicsdonotconvergetowardsanyperiodicorbit(m=0, ω =1.305,η=0.35andsolidbody-rotationε=0).Right:Exampleofanattractorcycleofthepathofcharacteristicsform=2,ω =−3.75, p p solid-bodyrotation(ε=0)andthesolaraspectratioη=0.71. (see the dispersion relation given in Eq. 25 for ε = 0). For a m=0 η=0.100 ε=0.300 givenaspectratiooftheshellη,certainvaluesofthefrequency 0.0 ω cause the paths of characteristics to converge towards peri- -0.5 p Λ -1.0 odicorbitscalled“attractors”,whileforsomeotherfrequencies, -1.5 paths of characteristics nearly occupy the whole shell with no -2.0 obviousperiodicpattern.Thesetwodifferentbehaviorsareillus- 0.5 1.0 1.5 2.0 trated in Fig. 2. The existence of one or several attractor(s) for ω p a given set of parameters may be investigated by studying the exponential convergence rate of two close characteristics along m=0 η=0.350 ε=0.300 multiplereflexions.ThisisquantifiedbytheLyapunovexponent 0.0 Λ(formoredetails,seeRieutordetal.2001;Baruteau&Rieu- -0.5 Λ -1.0 tord2013).TheLyapunovexponentcanbedefinedas -1.5 -2.0 Λ= lim 1 (cid:88)N ln(cid:12)(cid:12)(cid:12)(cid:12)dxk+1(cid:12)(cid:12)(cid:12)(cid:12), (29) 0.5 1.0ω 1.5 2.0 N→∞ N k=1 (cid:12) dxk (cid:12) p m=0 η=0.710 ε=0.300 where dx is the separation between two characteristics after k k 0.0 reflections.Inourcase,thisquantitycanbeevaluatedusingei- -0.5 therthereflexionsontheequator(ds ),ontherotationaxis(dz ) Λ -1.0 k k or on the inner (resp. outer) boundary of the shell dθin (resp. -1.5 k -2.0 dθout). Λ ≈ 0 corresponds to the case of “space-filling” paths k 0.5 1.0 1.5 2.0 of characteristics, whereas Λ < 0 identifies the existence of an ω attractor(seeleftandrightpanelofFig.2respectively). p Inthecaseofsolid-bodyrotation,thelocationswherechar- m=0 η=0.900 ε=0.300 acteristics bounce off the shell boundaries, the equator and the 0.0 rotationaxiscanbedeterminedanalytically,thusallowingfora -0.5 semi-analyticcalculationofΛ(Rieutordetal.2001).However, Λ -1.0 -1.5 differentialrotationrequiresthenumericalintegrationofthepath -2.0 ofcharacteristics.Forthatpurpose,wechooseastartingpointin 0.5 1.0 1.5 2.0 theshellandwefollowthepropagationbynumericallyintegrat- ω p ingEq.(13).Sincecharacteristicscannotpropagateintheelliptic domain,itisnecessarytoimposereflexionatturningsurfaces. Fig. 3. Numerical spectrum of the Lyapunov exponent for m = 0, Forillustrationpurposes,weevaluatednumericallytheLya- ε = 0.30 and ω ∈ [0.2,2.0]. From the top to bottom panels, η = punovexponentsasafunctionoftheeigenfrequencyform = 0 {0.10,0.35,0.71,0p.90}.Theredopensquarescorrespondtothemodes fordifferentvaluesofηandε,moreprecisely: showninFigs.6to8. – η={0.10,0.35,0.71,0.90}alongwithafixedsolar-likevalue ofthedifferentialrotationparameterε=0.30forFig.3, – ε = {−0.50,−0.25,0.30,0.60} along with a fixed solar-like Weused800frequenciesperunitintervaloffrequenciesand valueoftheconvectiveshellaspectratioη=0.71forFig.4. tenpairsofcharacteristicsforeachdatapointbeforeaveraging Articlenumber,page7of17 A&Aproofs:manuscriptno.GBMR2015-p m=0 η=0.710 ε=-0.500 entialrotationparameterε.Weseethattheoccurrenceofshort- 0.0 period attractors is not sensitive to ε but some features of the -0.5 Lyapunov spectra can be tracked as they are progressively al- Λ -1.0 teredandshiftedtohigherfrequencieswithincreasingε—for -1.5 -2.0 instance, the deep double-peak starting from ωp ∈ [0.3,0.5] in 0.5 1.0 1.5 2.0 thetop-leftpanelandshiftedtoωp ∈ [1.25,1.5]inthebottom- ω rightpanel. p m=0 η=0.710 ε=-0.250 Whenthefluidisdifferentiallyrotatingwithangularvelocity Ω(θ),pathsofcharacteristicsmaystillconvergetowardsattrac- 0.0 -0.5 tors — just like in the case of solid-body rotation — but with Λ -1.0 curved characteristics as illustrated in the top-left panel of Fig. -1.5 5. They also depend on the azimuthal wavenumber m through -2.0 thenon-uniformDoppler-shiftedfrequencyω˜ .Wecanalsofind 0.5 1.0 1.5 2.0 p attractorsofcharacteristicswhenaturningsurfaceexistsinthe ω p shell(top-rightpanel).Sometimes,characteristicstendtofocus m=0 η=0.710 ε=0.300 at the intersection of a turning surface with the boundaries of 0.0 the shell as illustrated by the bottom-left plot of Fig. 5, show- -0.5 ingabehaviorthatissimilartotheonefoundbyDintransetal. Λ -1.0 (1999)forgravito-inertialwaves.Finally,weemphasize thatin -1.5 rare cases, two turning surfaces can exist within the shell, af- -2.0 fectingthewaves’propagationdomainaccordingly:itiseither 0.5 1.0 1.5 2.0 restrictedtothepolarandequatorialregions(asillustratedinthe ω p bottom-right panel of Fig. 5) or restricted to mid-latitudinal re- m=0 η=0.710 ε=0.600 gions. 0.0 -0.5 Λ -1.0 -1.5 -2.0 0.5 1.0 1.5 2.0 4. Viscousproblem:shearlayers,comparisonto ω p inviscidanalysis,behaviouratcorotation resonances Fig. 4. Numerical spectrum of the Lyapunov exponent for m = 0, η = 0.71 and ω ∈ [0.2,2.0]. From the top to bottom panels, ε = p {−0.50,−0.25,0.30,0.60}. The blue solid line marks the frequency of Our study of the propagation properties of inertial waves in a thetransitionbetweenDandDTmodes.Theredopensquarescorre- differentially rotating inviscid fluid with a conical rotation pro- spondtothemodesshowninFigs.6to8. file(seesection3)showsthattwofamiliesofinertialmodesof oscillationexist:Dmodesthatmaypropagateintheentireshell, andDTmodesthataretrappedinpartoftheshellbecauseofthe theresults.TheresultsshowninFigs.3and4shouldtherefore existence of a turning surface. We find that paths of character- beconsideredratherasaqualitativeindicatorofthepresenceof isticsdependonmthroughtheDoppler-shiftedwavefrequency, short-periodattractorcyclesthanaquantitativeonebecauseour whichisnotthecaseforsolid-bodyrotation,butthemaindiffer- method probably yields substantial uncertainties on the values encebetweenm = 0andm (cid:44) 0modesisthepossibleexistence of Λ. Note also that the focusing of characteristics towards a ofcorotationresonancesinthelattercase. wedgeinthefrequencyrangeofDTmodes(seethebottom-left panel of Fig. 5 below) prevented us from computing values of InSect.4.1,webrieflypresentthenumericalmethodweused ΛfortheDTmodes.Asexpected,attractorsofvariousstrengths inordertosolvetheviscousproblemexposedinSect.2.1.Then existinsmallintervalsoffrequencieswhereΛisstrictlynegative we study in Sect. 4.2 (respectively Sect. 4.3) the general prop- whereas intervals where Λ ≈ 0 are associated with very long- ertiesofweakly-dampedviscousDmodesthatpropagateinthe periodattractorsorquasi-periodicorbitsofcharacteristics. wholesphericalshell(resp.DTmodesthatpropagateinpartof Fig.3illustrateshowtheexistenceofshort-periodattractors theshell)vianumericalsimulations.Weshowinparticularhow in the frequency range of inertial waves can be affected by the the shear layer structure can be related to the inviscid analysis size of the inner core. It seems that when the core is small, at- detailed above in Sect. 3. As explained in the introduction of tractorsofcharacteristicsareveryrarebutstillexistasillustrated this work, we will mostly focus on singular modes for which a bytheregionaroundω ≈ 1.66inthetop-leftpanelofFig.3. short-periodattractorofcharacteristicsexistssincetheymaybe p Wecheckedthatthisattractorexistsregardlessofthesizeofthe related with strong viscous dissipation (Ogilvie 2005). Finally, innercoreaslongasthelatterissmallenough,andwecomputed we examine m (cid:44) 0 modes with corotation resonances in Sect. thecorrespondingeigenmode(seeFig.8,andFig.8inBaruteau 4.4. &Rieutord(2013)).However,theotherpanelsofFig.3indicate thatmoreandmoreattractorsexistwithincreasingthevalueofη, duetoincreasinglynumerousreflectionsontheinnercore.The 4.1. Numericalmethod casewhereη = 0.90showsthatalargenumberofshortperiod attractorsexistoverthefrequencyrangeofDmodes. ThelinearizeddimensionlesssystemofEqs.(6)andtheassoci- Fig. 4 shows the evolution of the Lyapunov exponent spec- ated stress-free boundary conditions are solved using a unique trum for a given core size but for various values of the differ- decomposition of the unknown velocity field u onto vectorial Articlenumber,page8of17 Gueneletal.:Tidalinertialwavesinthedifferentiallyrotatingconvectiveenvelopesoflow-massstars-I 1.0 1.0 ε= 0.30 ε= -0.30 ω=-3.7500 ω= 1.6800 p p m= 2 m= 0 η= 0.71 η= 0.71 0.8 0.8 0.6 0.6 z z 0.4 0.4 0.2 0.2 0.0 0.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 s s 1.0 11..00 ε= 0.50 εε== --00..9900 ωp= 2.4000 ωωpp==--00..11550000 m= 0 mm== 22 η= 0.71 ηη== 00..7711 0.8 00..88 0.6 00..66 z zz 0.4 00..44 0.2 00..22 0.0 00..00 0.0 0.2 0.4 0.6 0.8 1.0 00..00 00..22 00..44 00..66 00..88 11..00 s ss Fig.5.Upperleft:ExampleofattractorcycleforaDmodewithm = 2,ω = −3.75inasolar-likeconvectiveenvelope(ε = 0.3,η = 0.71). p Upperright:SameforaDTmodewithm=0,frequencyω =1.68andanti-solarconicalrotation(ε=−0.3),thewhitedashedlineshowsthe p turningsurface.Bottomleft:Illustrationofthefocusingofthepathsofcharacteristicsattheintersectionoftheturningsurface(dashedline)and theouterboundaryoftheshell(m=0,ω =2.4andε=0.5).Bottomright:Form=2,ω =−0.126andε=−0.90,twoturningsurfacesexist p p withintheshell,allowingthepropagationofcharacteristicsneartheequatororthepolesbutnotinbetween. sphericalharmonics(Rieutord1987): 2. the right-hand side of the fourth-order differential equation satisfied by ul is a linear combination of wl±3, dwl±3/dx, (cid:88)∞ (cid:88)l (cid:110) wl±1, dwl±1/dmx, ul±2, dul±2/dx and d2ul±2/dmx2 (themse last u(x,θ,ϕ)= ulm(x)Rml (θ,ϕ)+vlm(x)Sml (θ,ϕ) thmreetermmsvanishmform=m0). m l=0 m=−l (cid:111) +wl (x)Tm(θ,ϕ) (30) Thefactthattheradialfunctionsul ,vl andwl withdifferent m l m m m with Rm = Ym(θ,ϕ)e , Sm = ∇ Ym and Tm = ∇ ×Rm where azimuthal wavenumbers m are notcoupled is aconsequence of Ym(θ,ϕl)isthel usualsrphelricalhaHrmlonicofldegreeHlandlorderm theaxisymmetryofthebackgroundflow.Thecouplingbetween l differentlatitudinaldegrees(froml−3tol+3)resultsfromthe normalizedontheunitsphere,and∇ =e ∂ +e (sinθ)−1∂ is H θ θ ϕ ϕ choiceofourspecificconicalrotationprofilethroughtheCorio- thehorizontalgradient. lisaccelerationterms. TheprojectionofthecontinuityequationonYm givesvl as l m In order to solve these equations numerically, they are dis- a simple function of ul and dul /dx. The momentum equation m m cretizedintheradialdirectionontheGauss-Lobattocollocation isthenprojectedonRm andTm,whichyieldstwolongordinary l l nodesassociatedwithChebyshevpolynomials.Eachsetofequa- differentialequationsinvolvingonlytheradialfunctionsul and m tions is truncated at order L for the spherical harmonics basis, wl (sincevl canbeeliminated): m m and at order Nr for the Chebyshev basis. This yields a finite 1. therighthandsideofthesecond-orderequationsatisfiedby eigenvalue problem of order N × L with a block banded ma- r wl is a linear combination of ul±3, dul±3/dx, ul±1, dul±1/dx trix composed of (up to) seven block bands. Then we use the m m m m m andwl±2(thistermvanishesform=0); linearsolverbasedontheincompleteArnoldi-Chebyshevalgo- m Articlenumber,page9of17 A&Aproofs:manuscriptno.GBMR2015-p rithm(seedetailsinValdettaroetal.2007)tocomputepairsof C as a function of the spherical harmonic degree : for a given l eigenvalues(equaltoiω )andeigenvectors(valuesoful ,vl and l, we take the maximum value among all the Chebyshev coef- p m m wl oneachpointoftheradialgrid),givenaninitialvalueguess ficients.Thereforewearecertainthatnumericalconvergenceis m foriω .SincethevaluesofN andLrequiredtoachievespectral achievedforthismode. p r convergencevaryalotfromonemodetoanother,especiallyfor Finally,theaxisymmetricmodeoffrequencyω ≈1.66pre- p theverydemandingsmallvaluesof E,somefigureswillbeac- sentedintheleftpanelofFig.8displaysashearlayerthatfol- companiedbythespectralcontentofthevelocityfield.Inallthe lowsanattractorofcharacteristicsthatexistsforarbitrarilysmall casespresentedinthisworkweassumedsymmetrywithrespect valuesoftheaspectratioη.Thismeansthatinconicaldifferen- totheequatorialplanewhichexplainsthatourresultsareshown tialrotation,attractorsofcharacteristicsmayexistindependently onlyforpositivevaluesofz,butanti-symmetryisalsopossible. oftheexistenceofaninnercore,whichwasalsofoundforcylin- drical and shellular differential rotation profiles by Baruteau & Rieutord (2013). We checked that this mode does exist for ar- 4.2. Axisymmetricandnon-axisymmetricDmodeswithno bitrarily small cores, which contrasts with the fact that inertial corotationlayer modesinafullsphereareregularinthecaseofsolid-bodyrota- tion(Greenspan1968;Zhangetal.2001).Ourresultshowsthat Wepresentinthissectionournumericalresultsforafewrepre- this is probably no longer the case when differential rotation is sentativeDmodeswithabackgroundconicalrotationprofileand taken into account. The case where characteristics do not con- Ekman numbers between 10−7 and 10−8, and we compare the vergetowardsanylimitcycleisshownintheright-panelofFig. structural properties of the shear layers with the inviscid anal- 8,whichdisplaysamodeoffrequencyω ≈1.23withη=0.71 ysis detailed in Sect. 3. Note that in the rest of this paper, we p and ε = −0.25. This set of parameters correspond to the open alwaysshowthekineticenergydistributionofthecomputedve- square in Fig. 4 for which Λ ≈ 0. As expected, the shear layer locity fields, but that the viscous dissipation looks qualitatively patternsvisibleherefollowthepropagationofcharacteristicsso verysimilar(e.g.Fig.13inBaruteau&Rieutord2013). that the kinetic energy is more smoothly distributed than in the Fig. 6 displays the result of one of our numerical calcula- previousattractorcases. tionsforanaxisymmetric(m = 0)Dmodewitheigenfrequency ThefourmodesshowninFigs.6to8allcorrespondtoone ω ≈ 0.91 obtained in a spherical shell of solar aspect ratio p oftheopensquaresinFigs.3and4.Thishighlightstheinterest η = 0.71withsolar-likeconicaldifferentialrotation(ε = 0.30). of the inviscid approach undertaken in Sect. 3 since the kinetic The left panel displays the spatial distribution of the mode’s energyofagivenmodeisindeedmorelocalizedalongthinshear kinetic energy in a meridional quarter-plane. The amplitude of layers when Λ < 0 and more uniformly-distributed in the shell the mode is maximum near the rotation axis which is a feature when Λ ≈ 0. As a conclusion, the inviscid analysis allows us sharedbyinertialmodesinauniformlyrotatingfluidshell.This tounderstandhowdifferentialrotationaffectsinertialmodes,al- has been demonstrated in detail in the appendix A of Rieutord though we restricted ourselves to D modes that can propagate & Valdettaro (1997) who showed that the kinetic energy along inthewholesphericalshell.WeseethatDmodesbehavequite characteristictrajectoriesgrowsass−1/2whens→0.Thestruc- similarly to inertial modes for solid-body rotation. In the next ture of this specific mode mostly consists of a shear layer fol- section we turn to DT modes, which are specific to differential lowingashort-periodattractorformedby(slightly)curvedlines, rotation. asexpectedinthecaseofdifferentialrotation(seeSect.3).This attractor is overplotted with a thick white curve, which is the predictionforthepropagationofcharacteristicsundertheshort- 4.3. Axisymmetricandnon-axisymmetricDTmodes wavelengthapproximation.Thepatternsformedbytheattractor and by the shear layers of the viscous mode are in very good In this section, we focus on modes for which a turning surface agreement.ThemodeshowninFig.6wasextractedfromase- exists in the shell. The frequency range in which these modes quence of calculations in which we followed a particular mode may exist is indicated by the white areas in Fig. 1. However, whileprogressivelydecreasingtheEkmannumberfrom10−6 to notallsetsofparametersfallinginthesewhiteareaseffectively 10−9. The damping rate of this mode is displayed as a function correspond to an eigenmode. The inviscid analysis exposed in of E in the right-panel of Fig. 6 where the dashed line clearly Sect. 3 reveals that characteristics cannot propagate in the el- showsthatitisproportionalto E1/3,whichisthescalingthatis liptic domain of the shell, i.e. where ξ < 0 or equivalently (cid:12) (cid:12) expectedintheasymptoticregimeforsolid-bodyrotation(Rieu- (cid:12)(cid:12)ω˜p(θ)(cid:12)(cid:12) > 2ω(θ).Asaconsequence—forsufficientlysmallEk- tord&Valdettaro1997),meaningthatthiskindofDmodeisnot mannumbers—DTmodesareexpectedtobetrappedinthehy- deeplyaffectedbydifferentialrotation. perbolic domain out of which characteristics cannot propagate, apropertywhichhasalwaysbeenverifiedbyournumericalcal- We show in Fig. 7 a qualitatively similar axisymmetric D mode of eigenfrequency ω ≈ 1.35 that was obtained with the culations.Wealsopresentresultsincaseswherecharacteristics p same aspect ratio but slightly lower Ekman number E = 10−8 converge towards the intersection of a turning surface with the andanti-solardifferentialrotationε = −0.25.Thistimetheam- inneroroutershellboundaries(seeFig.5). plitude of the mode is still maximum at the rotation axis but is also quite large near the inner critical latitude (where the char- 4.3.1. ExistenceofDTmodes acteristicsaretangenttotheinnercore),whichisreminiscentof the solid-body rotation case where shear layers are sometimes In order to find DT modes, as a first step we checked the dis- emittedatthecriticallatitude.TherightpanelofFig.7depicts tribution of the eigenvalues of our linear problem for a given the spectral content of u and w for this mode : the top panel setofparameters,usingaQZfactorisationmethod.Thismethod showsthemaximumChebyshevcoefficientsC asafunctionof rapidlybecomesverycostlyasthespectralresolution(andthus k theChebyshevorderk,takingthemaximumvalueamongallthe thedimensionsofthelinearproblem)increasesandthereforecan sphericalharmonicscoefficientsforagivenk;similarly,thebot- only be performed at a low spectral resolution, effectively lim- tompanelshowsthemaximumsphericalharmonicscoefficients itingustomoderatevaluesof E.Despitethislimitation,trends Articlenumber,page10of17