Mon.Not.R.Astron.Soc.000,000–000(0000) Printed4January2012 (MNLATEXstylefilev2.2) Density waves in debris discs and galactic nuclei Mir Abbas Jalali1,2(cid:63) and Scott Tremaine2† 1 Sharif University of Technology, Azadi Avenue, Tehran, Iran 2 School of Natural Sciences, Institute for Advanced Study, Einstein Drive, Princeton, NJ 08540, U.S.A. 4January2012 2 ABSTRACT 1 We study the linear perturbations of collisionless near-Keplerian discs. Such systems 0 aremodelsfordebrisdiscsaroundstarsandthestellardiscssurroundingsupermassive 2 black holes at the centres of galaxies. Using a finite-element method, we solve the n linearizedcollisionlessBoltzmannequationandPoisson’sequationforawiderangeof a disc masses and rms orbital eccentricities to obtain the eigenfrequencies and shapes J of normal modes. We find that these discs can support large-scale ‘slow’ modes, in 1 which the frequency is proportional to the disc mass. Slow modes are present for arbitrarily small disc mass so long as the self-gravity of the disc is the dominant ] sourceofapsidalprecession.Wefindthatslowmodesareoftwogeneraltypes:parent P modes and hybrid child modes, the latter arising from resonant interactions between E parentmodesandsingularvanKampenmodes.Themostprominentslowmodeshave . h azimuthal wavenumbers m = 1 and m = 2. We illustrate how slow modes in debris p discsareexcitedduringafly-byofaneighbouringstar.Manyofthenon-axisymmetric - features seen in debris discs (clumps, eccentricity, spiral waves) that are commonly o attributed to planets could instead arise from slow modes; the two hypotheses can be r t distinguished by long-term measurements of the pattern speed of the features. s a Key words: stellardynamics,methods:numerical,galaxies:kinematicsanddynam- [ ics, galaxies: nuclei, planets and satellites: formation, protoplanetary discs 2 v 1 5 5 1 INTRODUCTION ofthesefeaturescanbeproducedbyplanets(seeWyatt2009 4 forareview).Moreover,inthecaseofβ Pictoris(Lagrange Debrisdiscsareplanetesimaldiscsthataredetectedthrough . et al. 2010), and perhaps Fomalhaut (Kalas et al. 2008), 0 thermal infrared emission or scattered starlight from dust planets have been detected that may indeed be responsible 1 formed in recent planetesimal collisions. The bolometric lu- forsomeorallofthesefeatures.Nevertheless,itisimportant 1 minosity from detectable debris discs is typically (cid:38)10−5 of 1 to ask what long-lived structures could arise in debris discs thestellarluminosity,theinferreddustmassesaretypically v: (cid:46) 1M⊕, and the ages of the host stars range from 10 Myr withoInuttphliasnpeatps.er we examine the possibility that low-mass i to 10 Gyr (see Wyatt 2008 for a review). X discs can support long-lived normal modes maintained by A variety of features in debris discs have been inter- the self-gravity of the disc. Normally it is assumed that de- r preted as evidence for planets. These include structures in a brisdiscscannotsupportsuchmodesbecauseoftheirsmall the β Pictoris disc, including a warp (Heap et al. 2000), masses;alllocalizeddisturbancesaredispersedbytheKep- a system of tilted rings (Wahhaj et al. 2003) and a bright lerianshear.However,aspecialfeatureofKeplerianorbitsis clump(Telescoetal.2005);clumpsinthediscsaroundVega thateccentricorbitsdonotprecess.Thustheevolutionofec- (Wyatt2003),(cid:15)Eridani(Greavesetal.2005),ηCorvi(Wy- centricdisturbancesinadebrisdiscisgovernedbythenon- att et al. 2005), and HD 107146 (Corder et al. 2009); the Keplerianforces,howeversmallthesemaybe.Inthispaper eccentricity of the discs around HR 4796A and Fomalhaut we shall focus on the non-Keplerian forces arising from the (Telesco et al. 2000; Kalas et al. 2005); spiral structure in self-gravity of the disc. We neglect other possible perturba- thediscaroundHD141569(Clampinetal.2003);andsharp tionsforavarietyofreasons.Weignoregravitationalforces innerorouteredgesinthediscsaroundFomalhautandHD fromplanetsbecauseourprincipalgoalistounderstandthe 92945 (Kalas et al. 2005; Golimowski et al. 2011). propertiesofdiscsinthesimplestcase,whennoplanetsare Detaileddynamicalmodelshaveshownthatmostorall present. We ignore radiation pressure, even though this af- fects the dynamics of the dust that dominates the thermal (cid:63) [email protected](MAJ) infrared emission and the scattered light; our justification † [email protected](ST) is that the large planetesimals that generate the dust are (cid:13)c 0000RAS 2 M.A. Jalali and S. Tremaine unaffected by radiation pressure but we recognize that the arethepropertiesofthefrequencyspectraofnear-Keplerian distribution of (invisible) parent bodies and (visible) dust discs?(ii)Arethereanyunstablemodes?(iii)Arethereiso- is likely to be different. We ignore gas drag since old debris latedoscillatorymodesinthespectrumthatsurviveLandau discscontainlittleornogas,andsincetheplanetesimalsare damping?(iii)Whatarethedifferencesbetweenthespectra likelytobelargeenoughtobeinsensitivetodrag.Weignore of cold and warm discs? (iv) How can stable density waves collisions between planetesimals because they are likely to be excited in such discs? be rare; indeed such collisions probably drive the long-term We introduce a family of axisymmetric near-Keplerian erosion of the disc in which case the collision time cannot discs in §2 and construct their equilibrium phase-space dis- be much less than the disc age. tribution functions (DFs) in §3. We obtain the governing Debrisdiscsaredistinctfromprotoplanetarydiscs:the equations of the perturbed dynamics in §4 and explain the latterarecomprisedmostlyofgas,notdustorplanetesimals; numerical solution procedure in §5. We present the fre- they are much younger (typically less than a few Myr) and quency spectra of our discs in §6 and discuss the charac- moremassive(0.001–0.1M )thandebrisdiscs(seeWilliams teristicsofeigenmodesinwarmandcolddiscs.Wedescribe (cid:12) & Cieza 2011 for a review). Protoplanetary discs are de- how these waves can be excited by tidal forces in §7. The pletedbyvariousprocesses,includingphotoevaporation,ac- readerwhoismainlyinterestedintheapplicationofourre- cretion onto the host star, condensation of refractory ele- sultstodebrisdiscsandgalacticnucleicanfocusonFigures ments into dust grains and then planetesimals, and stellar 7 and 10 and the discussion in §8. winds. Eventually they are likely to evolve into planetesi- mal/debrisdiscs.Althoughouranalysishereisrestrictedto collisionless systems,many ofour results—in particularthe 2 THE MODEL existence of stable, slow, lopsided modes supported by the self-gravity of the disc—also apply to protoplanetary gas We introduce a simple model of annular discs around mas- discs and may explain some non-axisymmetric features of sive objects by subtracting two Toomre (1963) discs with these discs. n=1 and n=2; the resulting surface density is To summarize we treat debris discs as collisionless sys- (cid:26) (cid:27) 3M 1 1 temscomposedofparticlesinfluencedonlybythegravityof S (r)= d − , thehoststarandtheself-gravityofthedisc.Theirdynamics d 4πb2 [1+(r/b)2]3/2 [1+(r/b)2]5/2 is therefore similar to the dynamics of discs of stars orbit- 3M (r/b)2 = d , (1) ing near the supermassive black holes found in the centres 4πb2[1+(r/b)2]5/2 of most galaxies. Examples of these include the disc(s) of where M is the disc mass, b is a length scale, and r is the young stars found in the central parsec of the Milky Way d radialdistancetothecentralstar.Thepotentialcorrespond- (Genzel et al. 2010), the two discs—one of young stars at ing to the surface density S is ∼0.1pcandoneofoldstarsat∼1pc—foundatthecentre d of M31 (Bender et al. 2005), and the stellar discs that are GM 1+2(r/b)2 Φ (r)=− d , (2) inferredtoformintheouterpartsofquasaraccretiondiscs d 2b [1+(r/b)2]3/2 (Goodman 2003) with G being the gravitation constant. For a central star of The properties of the normal modes of low-mass near- massM ,thetotalpotentialgoverningthemotionofparti- Keplerian discs were investigated by Tremaine (2001; here- (cid:63) cles is after T01), who found that (i) the frequency of the mode is proportional to the ratio µ of the masses of the disc and Φ (r)=−GM(cid:63) +Φ (r). (3) central star, but the shape of the mode is independent of µ 0 r d so long as µ(cid:28)1 (hence these are called ‘slow’ modes); (ii) We define allslowmodesarestable;(iii)indiscswithrmseccentricity M e (cid:28)1allslowmodeshaveazimuthalwavenumberm=1, µ= d, R=r/b, (4) rms M i.e., they are lopsided. (cid:63) TheresultsinT01arebasedonlinearnormal-modecal- and work with the dimensionless unperturbed potential culations for discs composed of particles in circular orbits, bΦ 1 µ 1+2R2 withsoftenedself-gravityusedtomimictheeffectsoftheve- V0(R)≡ GM0 =−R − 2(1+R2)3/2, (5) locity dispersion or non-zero eccentricities of the particles. (cid:63) These calculations are supplemented by analytic results us- and density ing the WKB (short-wavelength) approximation, which as- b2S 3µ R2 sumes that the wavelengths of the normal modes are small Σ (R)≡ d = . (6) 0 M 4π(1+R2)5/2 comparedtotheradius.TheWKBresultsappeartoprovide (cid:63) ausefulguideeventhoughthisshort-wavelengthapproxima- The top panel in Figure 1 shows the radial profile of Σ /µ. 0 tionisnotrealisticforsomeofthediscmodes.Inthispaper The velocity of circular orbits, v (R), is determined c the effects of the velocity dispersion are computed directly, from andweexaminediscswitharangeofrmseccentricitiese , rms dV 1 µR2(2R2−1) fromnearlyzero(‘cold’discs)to∼0.35(‘warm’discs).Our v2(R)=R 0 = + . (7) c dR R 2 (1+R2)5/2 numerical results are derived using a finite-element method (FEM)forstudyingthelinearnormalmodesofcollisionless The second term on the right side of (7) becomes negative self-gravitatingdiscs,asdescribedinJalali(2010).Inpartic- forR2 <1/2.Thismeansthatourdiscscannotexistinthe ular,weintendtoaddressthefollowingquestions:(i)What absenceofacentralpointmass.Moreprecisely,v2 ≥0atall c (cid:13)c 0000RAS,MNRAS000,000–000 Density waves in near-Keplerian discs 3 from 0.05 2π (cid:73) dR Ω (J) J (cid:73) dR 0.04 = , φ = φ . (11) Ω (J) p (R,J) Ω (J) 2π R2p (R,J) R R R R µ 0.03 /0 In the limit µ → 0, the potential is Keplerian and we have Σ 0.02 Ω =Ω =a−3/2 withabeingtheorbitalsemi-majoraxis. R φ 0.01 For 0<µ(cid:28)1 the radial and azimuthal frequencies are no 0 longerequal,buttheirdifferenceΩpr =Ωφ−ΩRissmall,and -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 thatistheprecessionrateofthelineofapsides.TheTaylor log (R) 10 expansion of Ω begins with terms of O(µ) because Ω pr pr vanishes for Keplerian orbits. Consequently, for µ(cid:28)1, the e=0.0 precession rate is proportional to the disc mass. For nearly 0.005 e=0.4 e=0.6 circular orbits, the precession rate is given analytically by e=0.8 0 e=0.996 Ω = 3µR3/2(1−4R2) +O(µ2). (12) pr 4 (1+R2)7/2 -0.005 Insteadoftheactionsonemayusethesemi-majoraxis pr Ω a(J) and eccentricity e(J) defined by -0.01 R (J)+R (J) R (J)−R (J) a= min max , e= max min , (13) 2 R (J)+R (J) -0.015 min max whereR (J)andR (J)aretheminimumandmaximum min max -0.02 distancesofparticlesfromthecentralstar.Thesedefinitions areconsistentwiththestandardKepleriandefinitionswhen -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 the disc mass vanishes. In the bottom panel of Figure 1, log (a) 10 we have plotted the variation of Ω versus a for µ = 0.1 pr and several choices of e. It is seen that the precession rate Figure 1.Top:Surface-densityprofileofthecompositeToomre of orbits—of any eccentricity—has a positive peak within disc(eq.1).Bottom:VariationoftheprecessionrateΩpr forµ= 0.1.Forradialorbits(eccentricitye=1)Ωpr=0. the region where Σ0 is rising, and then switches sign and remains negative in the outer regions. The precession rate crosses through zero near a = 0.5 at all eccentricities. The maximumprecessionratefornearlycircularorbitsandµ(cid:28) Rifandonlyifµ≤55/2;thisisnotalimitationinpractice 1 is given by equation (12) as ω =0.05861µ, which occurs 0 since protoplanetary discs are expected to have µ(cid:28)1. atR=0.2859.In§6,weshallshowthatthepatternspeeds Werestrictourselvestorazor-thindiscssincetheverti- of stable waves are closely related to ω . 0 cal structure of thin discs should not strongly affect their large-scale response. Using the polar coordinates (R,φ) and their corresponding generalized momenta (pR,pφ), the 3 PHASE-SPACE DISTRIBUTION FUNCTION Hamiltonian function governing the motion of particles reads Particle orbits in collisionless discs are not necessarily cir- cular.Wethereforeconstructphase-spacedistributionfunc- H (p ,p ,R)≡E = p2R + p2φ +V (R). (8) tions (DFs) that enable us to distribute non-circular orbits 0 R φ 2 2R2 0 in the disc. We seek DFs of the form (Sawamura 1988; Pi- Since φ is a cyclic coordinate, its conjugate momentum p chon & Lynden-Bell 1996) φ isaconstantofmotionintheunperturbeddisc.Theorbital f (E,L)=L2K+2g (E), E =−E, (14) 0 K energyE isanotherintegralofmotion.Canonicalperturba- tiontheoriesdescribingthemotionofparticles,andtheper- where 0≤L≤Lc(E), Lc(E) is the angular momentum of a turbed collisionless Boltzmann equation, are substantially circular orbit with energy E =−E, and K is a positive in- simplified by using the action variables J = (J ,J ) and teger.ToreproducethesurfacedensitytheDFmustsatisfy R φ their conjugate angles w=(w ,w ) with the relation R φ 1 (cid:73) 1 (cid:73) Σ (R)=2(cid:90) ΨdE(cid:90) Lmax √f0(E,L) dL, Ψ=−V , (15) JR = 2π pRdR, Jφ = 2π pφdφ=pφ. (9) 0 0 0 L2max−L2 0 where L = R[2(Ψ−E)]1/2. Substituting (14) into (15) These integrals are taken along the orbits, which consist of max and performing the integral over L we find slowly precessing Kepler ellipses when µ (cid:28) 1. The unper- turbedHamiltonianH dependsonlyontheactions,notthe (cid:18) 3 1(cid:19)(cid:90) Ψ Σ (R) 0 2K+1B K+ , dEg (E)(Ψ−E)K+1 = 0 , (16) angles. The action Jφ =L is the magnitude of the angular- 2 2 0 K R2K+2 momentum vector. In the angle-action space, the equations where B(p,q) is the beta function. Taking the (K +2)th- of motion become orderderivativeofbothsidesof(16)withrespecttoΨgives w˙ =Ω(J)= ∂H0(J), J˙=0, (10) an explicit analytic form of gK, ∂J 1 dK+2 Σ (R) g (Ψ)= √ 0 . (17) and the orbital frequencies Ω(J)=(ΩR,Ωφ) are computed K π2K+1Γ(K+3/2)dΨK+2R2K+2 (cid:13)c 0000RAS,MNRAS000,000–000 4 M.A. Jalali and S. Tremaine One needs to know explicitly the function R(Ψ) before 1 doing the derivatives on the right side of (16). Since µ is 0.9 0 small in the discs we are considering, we utilize a pertur- 0.8 -1 bation method to compute R in terms of Ψ. Let us define u=1/R and rewrite (5) in the form 0.7 -2 0.6 -3 1 u(2+u2) e 0.5 -4 Ψ=u+µQ(u), Q(u)= . (18) 2(1+u2)3/2 0.4 -5 0.3 We now assume a formal series expansion for u in terms of 0.2 µ as (Bellman 1964) 0.1 u(Ψ)=u0(Ψ)+µu1(Ψ)+µ2u2(Ψ)+··· , (19) 0 -2 -1 0 1 2 andsubstitutethisinto(18).Thefunctionsu (Ψ)arerecur- log (a) j 10 sively determined by putting equal to zero the coefficients of µj (j = 0,1,2,···). The recursion begins with u = Ψ. 1 0 Up to the third-order terms, we obtain 0.9 0.8 0 u1 =−Q(u0), 0.7 -1 u2 =−u1Q(cid:48)(u0), 0.6 -2 1 e 0.5 -3 u =−u Q(cid:48)(u )− u2Q(cid:48)(cid:48)(u ), (20) 3 2 0 2 1 0 0.4 -4 whereQ(cid:48)(u)=dQ/du.Theseriesforuconvergesrapidlyso 0.3 keepingthetermsofO(µ2)isquitesufficientforcomputing 0.2 R(Ψ)=1/u(Ψ) in discs with µ≤0.1. 0.1 The functions g (E) admit negative values for K = 0 K -2 -1 0 1 2 0,1,buttheyarepositive-definiteandthereforephysicalfor log (a) 10 plausible values of µ<1 when K ≥2. We have plotted the contoursoflog10(f0/µ)using(a,e)asindependentvariables Figure 2. Contours of log10[f0(a,e)/µ] for µ = 0.1. The max- in Figure 2 for K = 5 and K = 29. The mean and rms imum surface density Σ0(R) occurs at R = 0.8165. There- eccentricity of the disc, e¯and e , are given by fore, the highest phase-space density appears in the vicinity of rms log (a)(cid:39)0.Top:K=5.Bottom:K=29. (cid:82) ef (J) d2J Γ(3)Γ(5 +K) 10 e¯= (cid:82) 0 = 2 2 +O(µ) f (J) d2J Γ(3+K) 0 e2rms =(cid:82)(cid:82)ef2f(0J(J))dd2J2J = 2K2+5 +O(µ). (21) 00..45 0 Larger values of K correspond to colder discs. For µ (cid:28) 1 ms 0.3 the mean eccentricity e¯= 0.329 for K = 5 and e¯= 0.159 er 0.2 for K = 29. When K (cid:29) 1 the DF at a given energy or 0.1 semi-major axis approaches the Schwarzschild or Rayleigh K=5 DF, 1.6 K=10 K=20 K=29 f0(e2)de2 ∝exp(−e2/e20)de2, e−02 =K+1/2, (22) σR 1.1 Inthislimitthemeanandrmseccentricityarerelatedtoe 0.6 √ 0 by e¯= πe /2, e =e . 0 rms 0 0.1 A necessary condition for stability to small-scale ax- isymmetric disturbances is that Toomre’s Q > 1; here 3 Q=σRΩR/(3.36Σ0) where σR is the radial velocity disper- Q 2 µ sion. The models in this paper with µ(cid:28)1 have Q>0.5/µ everywhere and thus are stable in this sense. The top two 1 panelsofFigure3showthermseccentricityandσR asfunc- 0 -1.5 -1 -0.5 0 0.5 1 1.5 tions of radius; for µ (cid:28) 1 these are independent of µ. The log (R) 10 bottom panel shows µQ which is also independent of µ for µ(cid:28)1.Noteinparticularthatthermseccentricityisalmost Figure3.Thermseccentricity,radialvelocitydispersion,andµQ independent of radius. (the disc/star mass ratio times Toomre’s stability parameter Q) asfunctionsofradius.Whenµ(cid:28)1allthreeplotsareindependent ofµ;thecurvesarefromnumericalmodelswithµ=0.1. 4 PERTURBED DYNAMICS We assume that gas drag, collisions, and other non- gravitationaleffectsarenegligiblesothedisccanbetreated (cid:13)c 0000RAS,MNRAS000,000–000 Density waves in near-Keplerian discs 5 as a collisionless fluid. We impose small-amplitude distur- 5 THE FINITE-ELEMENT METHOD bances to the surface density, potential, and DF: The dynamics and stability of collisionless discs are usually Σ(R,φ,t)=Σ0(R)+(cid:15)Σ1(R,φ,t), (23) studiedbyoneoftwonumericalmethods:(i)N-bodysimu- V(R,φ,t)=V (R)+(cid:15)V (R,φ,t)+(cid:15)V (R,φ,t), (24) lations(e.g.,Sellwood1987);(ii)expansionoftheperturbed 0 1 e gravitational potential in a set of basis functions, followed f(w,J,t)=f (J)+(cid:15)f (w,J,t), (25) 0 1 by the evaluation of a matrix representing the response of where (cid:15) (cid:28) 1 and Ve is an external perturbing potential, the disc to a given imposed potential (e.g., Kalnajs 1977). perhapsinducedbyabinarycompanion,anencounterwith Neitherofthesemethods,however,isidealforinvestigation a passing star, or the tidal field of the birth cluster. The of the oscillations and response of low-mass near-Keplerian perturbedsurfacedensityΣ1anditscorrespondingpotential discs, for several reasons: (i) slow oscillations are stable V1 are related through Poisson’s integral: (T01) and therefore more difficult to detect than growing (cid:90)(cid:90) Σ (R(cid:48),φ(cid:48),t)R(cid:48)dR(cid:48)dφ(cid:48) modes; (ii) slow oscillations have low frequencies, and thus V1(R,φ,t) = −G (cid:112) 1 N-body simulations must be followed for many dynamical R2+R(cid:48)2−2RR(cid:48)cos(φ−φ(cid:48)) times;(iii)low-massdiscsalsosupportshort-wavelengthfast (cid:90)(cid:90) Σ (R(cid:48),φ(cid:48),t)cos(φ−φ(cid:48))dR(cid:48)dφ(cid:48) + GR 1 ,(26) (i.e.,frequencyindependentofµ)oscillationsandthesecan- R(cid:48) not be resolved without a large set of basis functions; (iv) andweconsiderself-consistentdensityperturbationssothat weshallfindthatsomeslowoscillationshavenearlysingular (cid:90) components.Here,weadoptafinite-elementmethod(FEM) Σ1 = f1 d2v. (27) and reduce the linearized CBE to a system of ordinary dif- ferential equations that describes the temporal evolution of The second term on the right side of (26) is the indirect the disc, both the eigenfrequency spectrum of an isolated potentialperturbationthatarisesbecauseweareworkingin disc and the response of a disc to external perturbations. anon-inertialreferenceframecentredonthestar.Itisnon- We use a C FEM (all functions are continuous, but not 0 zero only for m=1 perturbations since perturbations with necessarily differentiable at boundaries between elements) m(cid:54)=1leavethecentreofmassofthediscunchanged.Fora in the configuration space. particlewithactionsJ,theradialdistanceR andexp(imφ) Inthissection,webrieflyreviewtheprinciplesofFEM can be expanded as Fourier series in the angle variables, modelling. For a general introduction see Zienkiewicz et +∞ +∞ al. (2005). Detailed descriptions of the application of an R= (cid:88) ξ(J)eilwR, eimφ =eimwφ (cid:88)η(J)eilwR. (28) FEMtocollisionlessself-gravitatingsystemscanbefoundin l l l=−∞ l=−∞ Jalali(2010)forperturbedsystemsandinJalali&Tremaine Any function of R and φ that is 2π-periodic in φ can thus (2010) for equilibrium models. beexpressedinthe(w,J)coordinates.FortheHamiltonian function that governs the motion of particles, we write 5.1 Finite ring elements in the configuration space H=H (J)+(cid:15)V (w,J,t)+(cid:15)V (w,J,t). (29) 0 1 e We split the configuration space into N ring elements. The where H0 is defined in equation (8). Therefore, the per- nth element is characterized by its nodes at Rn and Rn+1, turbed equations of motion become and by a linear interpolating vector Gn(R) defined by ∂H ∂ G =(cid:2) G G (cid:3), G =1−R¯, G =R¯, (35) w˙ = =Ω(J)+(cid:15) (V +V ), (30) n 1,n 2,n 1,n 2,n ∂J ∂J 1 e J˙=− ∂∂Hw =−(cid:15)∂∂w(V1+Ve). (31) wwehearreeRi¯nt=ere(Rste−dRonnl)y/∆inRlinneaanrdp∆erRtunrb=atRionn+s,1d−isRtunr.baSninccees of different azimuthal wavenumber m are independent. For It is obvious that the actions vary slowly in the perturbed thewavenumberm,thepotentialV andthesurfacedensity disc. Subtracting the evolutionary equations of w and w 1 φ R Σ are thus computed from gives the apsidal precession rate in the perturbed disc, 1 (cid:18) (cid:19) N w˙φ−w˙R =Ωpr(J)+(cid:15) ∂∂J − ∂J∂ (V1+Ve). (32) V1(R,φ,t)=Re(cid:88)Hn(R)Gn(R)·an(t)eimφ, (36) φ R n=1 Since Ω =O(µ), for low-mass discs (µ(cid:28)1) this equation N containsprtwo small parameters, µ and (cid:15). Σ1(R,φ,t)=Re(cid:88)Hn(R)Gn(R)·bn(t)eimφ. (37) The DF in the perturbed disc obeys the collisionless n=1 Boltzmann equation (CBE), The function H (R) is unity for R ≤ R ≤ R and zero n n n+1 df ∂f otherwise. The column vectors = +[f,H]=0, (33) dt ∂t a =(cid:2) a a (cid:3)T, b =(cid:2) b b (cid:3)T, n n1 n2 n n1 n2 where [·,·] denotes a Poisson bracket. Here we confine our- contain the nodal potentials and densities, respectively. selves to the linearized equation: According to the definition of G (R), Σ is equal to n 1 ∂f1 +[f ,H ]+[f ,V ]=−[f ,V ]. (34) Rebn1exp(imφ) at R = Rn and to Rebn2exp(imφ) at ∂t 1 0 0 1 0 e R = R . Similarly, the nodal potentials at these radii n+1 The remainder of this paper is devoted to the study of so- involve a and a . The perturbed surface density and its n1 n2 lutionsofthisequationandtheirapplicationtocollisionless correspondingpotentialarecontinuousanddifferentiablein- near-Keplerian discs. side elements and the continuity of these functions at the (cid:13)c 0000RAS,MNRAS000,000–000 6 M.A. Jalali and S. Tremaine boundaries of elements (nodes of rings) implies HereU ,U andU areconstantsquarematricesofdimen- 1 2 3 sionN ×N ,andZ(t)isanN -dimensionalcolumnvector, t t l t a =a , b =b . (38) n2 n+1,1 n2 n+1,1 which is the Galerkin projection of −[f ,V ]. 0 e This means that for a given m we have N = N +1 inde- OnecanalsoverifythattheGalerkinprojectionsof(26) t pendent nodal potentials/densities. and (27) respectively become Theangle-actionrepresentationoftheperturbedpoten- (cid:88) p(t)=C·d(t), d(t)= F(l)·z(t). (48) l tial V reads 1 l V (w,J,t)=Re (cid:88)+∞ h˜ (J,t) ei(lwR+mwφ), (39) The constant matrices C and F(l) are of dimension Nt × 1 1,l N . We combine (47) and (48) to express p(t) in terms of t l=−∞ z(t), and transform (47) to a non-homogeneous ordinary l where differential equation for z(t): l N dz(t) h˜1,l(J,t) = (cid:88)Ψl(n,J)·an(t), (40) dlt = −iU−11(l)·U2(l)·zl(t)+iU−11(l)·Zl(t) n=1 +∞ Ψ(n,J) = 1 (cid:73) H (R)G eim(φ−wφ)e−ilwR dw . (41) + (cid:88) iU−11(l)·U3(l)·C·F(l(cid:48))·zl(cid:48)(t), (49) l 2π n n R l(cid:48)=−∞ The external disturbance V can also be expressed in terms for l,l(cid:48) =0,±1,±2,···. By defining e ofangleandactionvariables.TocomputetheperturbedDF z(t)=(cid:2) ... zT zT zT zT zT ... (cid:3)T, (50) f (w,J,t),weuseFourierseriesofanglevariablesandwrite −2 −1 0 +1 +2 1 and collecting the elements of U−1(l) · Z(t) (for all l = N ∞ 1 l f (w,J,t)=Re(cid:88) (cid:88) E(n,J)·zn(t) ei(lwR+mwφ), (42) 0,±1,±2,···) in the global forcing vector F(t), the system 1 l l (49)canbecastintothestandardformoflinearevolutionary n=1l=−∞ equations: where d (cid:2) (cid:3) z(t)=−iA·z(t)+iF(t). (51) El(n,J)= El1(n,J) El2(n,J) , (43) dt isaninterpolatingvectorintheactionspace(tobespecified In the absence of external disturbances, F(t) = 0, the cor- in §5.3) and responding homogeneous equation admits a solution of the form z(t)=exp(−iωt)c that yields the linear eigensystem: zn =(cid:2) zn zn (cid:3)T, (44) l l1 l2 A·c=ωc. (52) isacolumnvectorofto-be-determinedDFswhoseelements Wefindthespectrumofω usingHessenbergtransformation should satisfy the continuity condition ofAfollowedbyQRfactorization.Theeigenvectorconjugate zn =z(n+1). (45) to a given eigenfrequency ωj is then computed using the l2 l1 singular value decomposition Equation(42)calculatesthedistributionofperturbedorbits A−ω I=VT·W·V , (53) based on their passage through ring elements in the config- j 1 2 urationspace.Ifanorbitstaysonlyinsidethenthelement, where W is a diagonal matrix whose elements are the sin- its DF becomes gular values of A−ω I, and I is the identity matrix. The j fˆn(w,J,t)=Re (cid:88)∞ El(n,J)·znl(t) ei(lwR+mwφ). (46) cisoltuhmeneigoefnVve2ctcoorrrze(sjp)oansdsoincgiatteodthweitshmωajl.lest singular value l=−∞ In general, eccentric orbits may visit more than one ring 5.3 Interpolating functions in the action space element. The summation over n in (42) takes this behavior into account. InourC0 FEManalysis,thelocalinterpolatingvectorfunc- tionsG (R)(alsoknownasshapefunctions)canreconstruct n thespatialprofileofanyoscillatorywavewhosewavelength issufficientlylargecomparedtothesizesofelements.How- 5.2 Projected evolutionary equations ever,wemustalsointerpolatef intheactionspace,which 1 Weusetheconditions(38)andassemblethenodaldensities requiresdefiningtheinterpolatingvectorsEl(n,J)(eq.42). b (t)andpotentialsa (t)intheglobalN -dimensionalvec- To do this we use arbitrary dynamic solutions of the lin- n n t torsd(t)andp(t),respectively.Similarly,zn(t)arecollected earized collisionless Boltzmann equation, which should be l in z(t). We now take the inner product of (34) with an adequate representation of the DF for the purposes of l interpolation.Intheangle-actionspace,andusingequation e−i(l(cid:48)wR+mwφ)[El(cid:48)(n(cid:48),J)]T, (39), one can show and integrate the resulting systems of equations over the +∞ angle-actionspacetoobtaintheGalerkin-weightedresidual eimφGn(R)=V˜1(n,J,w)= (cid:88) Ψl(n,J) eilwR+imwφ. (54) form of (34) as l=−∞ We assume ∂f /∂t=−iγf , substitute V˜ (n,J,w) into the dz(t) 1 1 1 U1(l)· dlt =−iU2(l)·zl(t)+iU3(l)·p(t)+iZl(t). (47) linearized CBE and solve the resulting equation to obtain (cid:13)c 0000RAS,MNRAS000,000–000 Density waves in near-Keplerian discs 7 the interpolating vector in action space (cf. eq. 42) 6 PROGRADE WAVES Our finite-element mesh is uniform in log radius, R =10−α1+α2y(n,N), (56) n 1 n−1 El(n,J)= l∂f0l/Ω∂J+R+mΩm∂−f0/γ∂JφΨl(n,J). (55) y(n,N)=2(N +1) + N +1, n=1,2,··· ,N. (57) R φ The numerical computation of the Fourier coefficients Ψ(n,J) and then the interpolating vectors E(n,J) needs l l a mesh in the (a,e) space. Such a grid is not arbitrary be- cause there must be at least one orbit that visits the nth The physical eigenfrequency ω will thus be equal to ω +γ c ring element in configuration space. We fulfill this require- with ω being the computed eigenvalue of (52). Varying γ c ment using the following two-dimensional grid tests the robustness of our numerical methods since the re- sultsshouldbeindependentofγ.Ourtestsshowthatingen- (a ,e )=[R ,y(j,M )], i j i e eral our results are insensitive to variations of γ. Neverthe- less, the choice γ (cid:38) max[Ω (J)] offers better performance, wherethegridpointsalongthea-directionexactlycoincide pr with the boundary nodes of the mesh in the configuration particularly in colder models. We originally used γ = 0, space, and there are j = 1,2,··· ,M + 1 grid points in correspondingtotheuseofstaticCBEsolutionsasinterpo- e the e-direction. A circular orbit is indeed assigned to each latingvectorsinactionspace,butwiththischoicewefound boundary of a ring element. This is particularly helpful in occasional spurious growing modes. colddiscswhereonemustinterpolatethepopulationofcir- We remark that we do not generate a finite-element cular orbits. The parameters α and α are chosen so that 1 2 meshinactionspace,fortworeasons:(i)toreducethesizeof thecomputeddiscmass4π2(cid:82) f(J)dJusingthegridpoints the Galerkin-weighted evolutionary equations; (ii) to avoid in the (a,e)-space agrees with the actual disc mass within creatingspuriousgrowingmodes.Thesecondoftheseprop- 1%. erties has a straightforward mathematical explanation: the In this paper we focus on slow modes with azimuthal numberofreachableeigenmodesintheconfigurationspaceis wavenumber m = 1. Slow modes exist with larger m, so equaltothenumberofindependentnodalvariables,whichis long as the calculation includes Fourier terms with index N+1inourFEManalysis.However,theeigenvalueproblem l=−m. In particular, we have found a number of isolated, (52)hasbeenformulatedinthephasespaceandthenumber non-singular m=2 modes; these are present only if we use ofcomputedeigenmodesisequalto(N+1)×(l −l +1) max min a fine FEM mesh, since they are more compact and have wherel <0andl >0arethelowerandupperbounds min max shorterwavelengthsthanthem=1modes.Thewavelengths in the l sums. Since the nodal densities d are related to z l of m = 2 modes shrink to zero as the disc becomes colder throughequation(48),therewillbe(N+1)×(l −l ) max min (see Appendix). This behavior is expected since the only computedeigenmodesmorethanN+1eigenmodesthatthe large-scale slow modes in cold low-mass discs have m = 1. dimension of d determines. The extra modes should there- Wefoundnounstablemodes,whichisalsoexpectedforlow- fore overlap in groups of (l −l ) members to avoid max min mass discs (T01). spuriousmodes.Thishappensinournumericalcalculations We began our calculations with N = M = 70 and e performed in §6 when the Fourier expansions over w con- R l = −1, and increased the number of Fourier terms and verge inside all ring elements, as is expected for the recon- ring elements until the eigenfrequencies of stable modes structionofV1 andf1 inthe(wr,wφ)-subspace.Generating foundfrom(52)convergedtoafractionalaccuracyof10−4. a finite-element mesh, let us say with N nodes in the two- a Typically this required computing all Fourier terms with dimensional J-space, will result in 2N ×(l −l +1) a max min −2 ≤ l ≤ 3 and a grid with N = 160 and M = 140 e modes,butassumingtheconvergenceofFourierseries,only (N = 180 and M = 140 for the models with the lowest e (N+1)groupsofthemwillcorrespondtoeigenmodesinthe rms eccentricity, corresponding to K = 29). We have also configuration space. Consequently, 2N −N +1 computed a experimented with including terms with larger values of |l| modeswillbespurious,andourcalculationsshowthatsuch butthesehadonlyasmalleffectonourresults.Takinggrid spurious modes are growing. Working with 2N = N +1 a pointsintheregionswithtinyvaluesoff(a,e)(seeFigure2) will not help because it does not necessarily guarantee the leadstolargeerrorsinthepropertiesofthecalculatedden- convergence of FEM model in the action space. sity waves because the FEM discretization errors become OnlyfewmodesoutofN+1possiblestatesinthecon- larger than the absolute magnitudes of physical quantities. figurationspace(see§6)arephysical.Therestareeithersin- We evade this difficulty by generating the FEM mesh only gular,ordonotsatisfytheboundaryconditionsasR→∞. in the annular region 0.01≤R≤100 using the parameters Note that for frequencies ω that lie between the maximum (α ,α )=(2,4) in equation (57). 1 2 and minimum of the precession frequency Ω the singular All non-singular eigenmodes with m=1 were found to pr modesmaybevanKampenmodes(restrictedtothesurface beprograde(ω>0).Wefindtwogeneraltypesofmodes:a in action space on which ω=Ω ) which would be damped parent family that is already present when only the l=−1 pr by the Landau mechanism. However, not all modes with Fouriercomponentisincludedinthecalculation,andachild frequency in this range are necessarily van Kampen modes, familythatbifurcatesfromtheparentfamilyasmorel-terms since the mode may not be produced by the orbits associ- are included. The eigenfrequencies of child modes are very atedwiththatresonance.Thusdiscretemodesmayoverlap closetothoseoftheirparentmode(typicallywithin1–2%). in frequency space with continuous modes. They emerge from resonant interactions between two ap- (cid:13)c 0000RAS,MNRAS000,000–000 8 M.A. Jalali and S. Tremaine 0.4 avoid overcrowding the diagrams. Although the maximum µ=0.025 A5 A3 A1 precessionrateω0isproportionaltoµ,thespectraofω¯agree 0.3 to within 1% in the models with µ = 0.025 and µ = 0.05. B5 B3 B1 This scaling shows that our results can be directly applied _e 0.2 C5 C3 C1 toalldiscswithmassratiosµ(cid:28)1,inparticulartothetiny mass ratios µ(cid:46)O(10−3) characteristic of debris discs. D5 D3 D1 Figure 4 shows that the modes become more closely 0.1 spaced as their frequency decreases and the minimum fre- quencyineachspectrumisanaccumulationpoint.Thisim- E5 E3 E1 µ=0.05 plies the existence of prograde waves with arbitrarily short wavelengths. There is also a nice correlation between the 0.3 F5 F3 F1 precessionrateofthemosteccentricorbitinthemodel(see _e 0.2 G5 G3 G1 Figures 1 and 2) and the lowest frequency in the spectrum. Models with highly eccentric orbits have an accumulation H5 H3 H1 point of lower frequency. Figure 4 shows that the number 0.1 of modes increases with decreasing e¯. In the limit of e¯→0, however, dispersion-supported waves (or p-modes) can not WKB exist according to WKB results (T01). The frequency spacing between modes B and B is 1 2 0.3 larger than the spacing between C and C , which in turn 1 2 is larger than the spacing between D and D (which is 1 2 _e 0.2 so small that the two points are indistinguishable in the Figure).SimilarbehaviorisseenintheF,G,andH families 0.1 in the middle panel of the Figure. In the limit e¯→ 0, the parent modes tagged with the numbers 2j +1 and 2j + 00.7 1 2 3 4 5 2 (j = 0,1,2,···) become degenerate. In the language of ω/ω 0 T01,theyformadegenerateleading/trailingpairofp-modes (seeAppendix).Thepairingprocessbeginsfrommodeswith Figure 4. Eigenfrequency spectra of stable, prograde density highest pattern speeds, for the resonant cavities of those wavesinnear-Kepleriandiscs.Onlyparentmodesareshown.The modes are fed mostly by near-circular orbits, which are the verticalaxisisthemeaneccentricitye¯andthehorizontalaxisis ω¯ = ω/ω0 where ω0 = 0.05861µ+O(µ2). These models corre- only population used in WKB analysis. The child modes spondtoDFsoftheform(14)withK=5,10,20,29.Thecalcula- of degenerate pairs also disappear because their supporting tionincludesFouriertermsl=−2,−1,0,1.Notethelogarithmic eccentric orbits disappear as e¯→0. Modes with ω¯ →1 and scale of the horizontal axis. Top: µ = 0.025 and ω0 = 0.00146. sufficiently large e¯ engage highly eccentric orbits and thus Middle: µ = 0.05 and ω0 = 0.00293. The eigenfrequencies of leadtomorecomplexdynamics.Eccentricorbitsareindeed modes D1 and D2 are very close and indistinguishable in the the backbones of discs, and when perturbed, they affect a plots.Theyareω¯D1 =4.659andω¯D2 =4.647.Similarly,wehave vast radial domain while near-circular orbits have only a ω¯H1 =4.700andω¯H2 =4.689.Notethesimilarityofthespectra local influence on developing patterns. inthetopandmiddlediagramsdespitethechangeofafactorof Wenowexaminetheshapesofthemodes.Afterfinding twointhediscmassµ;thisfeatureischaracteristicofslowmodes. ω,wecalculateitscorrespondingeigenvectorc,andusethis Mode shapes associated with the labelled frequencies have been to compute the nodal potentials p and nodal densities d plottedinFigures5,6and7.Bottom:Eigenfrequencyspectrade- rived from the WKB approximation described in the Appendix. from (48). Defining Each plotted point represents a degenerate leading/trailing pair N ofmodes. (cid:88) X(R) = Re H (R)G (R)·b , (58) n n n n=1 N (cid:88) proximatemodesthatareweaklycoupled:theparentmodes Y(R) = Im H (R)G (R)·b , (59) n n n and singular van Kampen modes. For l = 0 and l = +1 n=1 the singular components of the child modes correspond to one can compute the perturbed density patterns thecorotation(CR)andouterLindblad(OLR)resonances, respectively. The coupling between slow and van Kampen Σ (R,φ,t)=X(R)cos(mφ−ωt)−Y(R)sin(mφ−ωt), (60) 1 modesisprobablyduemostlytohighlyeccentricorbitsthat forasinglewavenumberm.Notethatb areextractedfrom areperturbedbythegravityfrombothwaveforms.Themain n the elements of d using the following formula: evidence for this is that as the mean eccentricity e¯shrinks, child modes collapse to singular modes and disappear. b =(cid:2) d d (cid:3)T, n=1,2,··· ,N. (61) n n n+1 We denote the maximum precession rate of circular orbits by ω = max[Ω (J)]; from equation (12) ω = Figure5showstheprofileofX(R)forthelabelledpar- 0 pr 0 0.05861µ+O(µ2).Wethenplotourresultsusingthenormal- entmodesofFigure4.Notonlyarethenormalizedfrequen- izedfrequencyω¯ =ω/ω .Figure4showstheeigenfrequency cies of the modes A and E identical, but also their mode 0 j j spectraofprogradem=1parentmodesforthemassratios shapes are very similar. These remarks apply to the pairs µ=0.025andµ=0.05,andforfourvaluesofthemeanec- (B ,F ), (C ,G ) and (D ,H ) as well, and demonstrate j j j j j j centricitye¯.Thefrequenciesofchildmodesarenotshownto that the waveforms are independent of µ so long as µ(cid:28)1, (cid:13)c 0000RAS,MNRAS000,000–000 Density waves in near-Keplerian discs 9 K=5 K=10 1 AE11 1 BF11 0.5 0.5 R) 0 R) 0 X( X( -0.5 -0.5 -1 -1 1 AE33 1 BF33 0.5 0.5 R) 0 R) 0 X( X( -0.5 -0.5 -1 -1 1 AE55 1 BF55 0.5 0.5 R) 0 R) 0 X( X( -0.5 -0.5 -1 -1 0.01 0.1 1 10 0.01 0.1 1 10 R R K=20 K=29 1 GC11 1 DH11 0.5 0.5 R) 0 R) 0 X( X( -0.5 -0.5 -1 -1 1 GC33 1 DH33 0.5 0.5 R) 0 R) 0 X( X( -0.5 -0.5 -1 -1 1 GC55 1 DH55 0.5 0.5 R) 0 R) 0 X( X( -0.5 -0.5 -1 -1 0.01 0.1 1 10 0.01 0.1 1 10 R R Figure 5. Perturbed density components X(R) (cf. eq. 60) for some stable modes in near-Keplerian discs of µ = 0.025 (solid lines) and µ=0.05 (dotted lines). In several panels the dotted line is not visible because it lies under the solid line. There are N =160 ring elementsintheconfigurationspaceforK=5,10and20,andN =180elementsforK=29.Filledcirclesmarkthelocationsofelement nodes.Inallpanels,themaximumof|X(R)|hasbeennormalizedtounity. as one would expect for slow modes. The figure also shows blad resonances, respectively. In low-mass discs, these reso- thatthewavelengthofoscillationsincreaseswiththepattern nances are at large radii where the surface density is small, speed ω in a given disc, and decreases as the mean eccen- so the singular component of a child mode involves only a tricity of the disc shrinks. The number of nodes increases small fraction of the mass involved in the parent mode. As as the frequency decreases. An interesting property of the the disc mass shrinks to zero the child modes merge with wavesshowingmultiplenodesisthattheirdensitypeaksare the parent mode. The reason is that the eigenfrequency of approximately equally spaced in logarithmic scales. the parent mode is proportional to the disc mass so with very small disc masses the corotation and outer Lindblad Thechildmodesarehybridmodesthatinheritthefea- resonances are at extremely large radii where the surface tures of their parents in the central regions of the disc, but density is negligible. Thus the distinction between parent have a spike at the location of singular modes that cou- and child modes is unimportant for low-mass discs such as ple to them. Figure 6 displays the parent mode D of fre- 8 debris discs. quencyω¯ =2.083(seeFigure4)anditschildrenD with 8,CR ω¯ =2.0765andD withω¯ =2.0728,whichcontainsin- Figure 7 displays shaded contour plots of the pattern 8,OLR gularvanKampenmodesatthecorotationandouterLind- ofΣ (R,φ,t)forsomemodelswithµ=0.025(modeshapes 1 (cid:13)c 0000RAS,MNRAS000,000–000 10 M.A. Jalali and S. Tremaine K=29 Lindblad and corotation resonances are at very large radii where the disc surface density is negligible. 1 Parent D8 0.5 X(R) 0 7 EXCITATION OF OSCILLATORY WAVES -0.5 Protostarsliveintheharshenvironmentsoftheirbirthclus- -1 1 Child D8 (CR) ters.SimulationsoftheOrionnebula(Scally&Clarke2001) show that about 10 per cent of stars can have encounters 0.5 closer than 100 AU within 107 years. Such encounters can X(R) 0 excite waves in planetesimal/debris discs. Encounters were -0.5 invokedasapossibleexplanationfortheasymmetriesinthe βPictorisdebrisdiscbyKalas&Jewitt(1995)andLarwood -1 1 Child D8 (OLR) &Kalas(2001)buttheseauthorstreatedthedebrisdiscas 0.5 a collection of test particles, which can give misleading re- X(R) 0 spureltcsesssiniocne.theself-gravityofthediscdominatestheapsidal -0.5 Since our goal here is only to illustrate this process we -1 confineourselvestoin-planeparabolicencounters.Consider 0.01 0.1 1 10 100 a disc particle orbiting around a star of mass M , and as- R (cid:63) sumeaperturberofmassM .Asinearliersections,wescale p alllengthssothatthedisclengthscalebisunity,anddenote Figure 6. Perturbed density components X(R) for the par- thenormalizedpositionvectorsoftheparticleandperturber ent mode D8 and two of its associated child modes D8,CR and D8,OLR, in which the parent mode is coupled to singular modes (with respect to the host star) by R and Rp, respectively. at the corotation (CR) and outer Lindblad (OLR) resonances. The equation of motion for a disc particle is Themodelparametersareµ=0.025andK=29. d2R =−∇[a ·R+V (R)+V (R,t)+V (R,R )], (62) dt2 (cid:63) 0 1 p p where a is the acceleration vector of the host star in corresponding to µ = 0.05 are similar). It is seen that the (cid:63) an inertial frame, V is the unperturbed gravitational po- wave packets are more radially compact in the colder (K = 0 tential due to M and the self-gravity of the disc, V is 29) model than warmer (K =5,10) ones. (cid:63) 1 the perturbed self-gravitational potential of the disc, and The properties of these modes can be explored using V =−(M /M )/|R −R| is the potential field of the per- the short-wavelength or WKB approximation described in p p (cid:63) p turber. The gradient ∇ is taken over the R space, and the the Appendix. The validity of this approximation requires normalizedtimetisrelatedtotheactualtimet through k>h/Rwherekisthewavenumberandhisadimensionless actual t/t =(GM /b3)1/2. numberoforderunity.Iftwoadjacentnodesofawaveareat actual (cid:63) R and R then (cid:82)R2kdR =π so the condition for validity We assume that a(cid:63) is due to the encounter; thus we 1 2 R1 ignore the cluster’s tidal field. Consequently, oftheWKBapproximationmaybewrittenπ>hlogR /R 2 1 or log10R2/R1 <1.36/h. Inspection of Figure 6 shows that a ·R=(cid:18)Mp(cid:19)Rp·R, R =|R |. (63) for h = 1 this condition is satisfied by all of the modes we (cid:63) M R3 p p (cid:63) p have computed, though not by much in some cases. For a distant encounter, R (cid:28) R , the potential V can be The bottom panel of Figure 4 shows the WKB fre- p p expanded as the following series quency spectrum. Each point corresponds to a degenerate potahirerofomf tordaeilsi,nogn.eTchoemsepomsoeddeosfaleraisdeinfgrosmpirwaalvweasvienstahnedrtehse- Vp =−(cid:18)MMp(cid:19)R1 (cid:88)∞ (cid:18)RR (cid:19)iPi[cos(φ−φp)], (64) (cid:63) p p onant cavities defined by the closed frequency contours in i=0 R ·R Figure A1. The WKB approximation correctly reproduces cos(φ−φ )= p , several striking features of the FEM frequency spectra: (i) p RpR all modes are prograde (ω > 0); (ii) both the number of where φ and φ are, respectively, the azimuths of the disc p modes and the maximum frequency grow as the mean ec- particle and perturber measured from an inertial reference centricitye¯ofthediscshrinks;(iii)thereisanaccumulation line, and P are Legendre polynomials. The effective poten- i point of modes near ω/ω0 = 1 in the FEM spectra and at tial due to flying-by perturber thus reads ω/ω =1 in the WKB spectra; (iv) there is also reasonable 0 quantitative agreement between the frequencies derived by Ve = a(cid:63)·R+Vp, thetwomethods,atleastforthediscswiththelowestmean = −(cid:18)Mp(cid:19) 1 (cid:88)∞ (cid:18) R (cid:19)iP [cos(φ−φ )], (65) eccentricity.TheWKBanalysisintheAppendixfailstofind M R R i p (cid:63) p p thechildmodes,fortworeasons:(i)itisbasedontheepicy- i=2 cle approximation, which assumes that the eccentricity is where we have dropped the i = 0 term in (64) because it smallandthusneglectsthehighlyeccentricorbitsthatcou- makes no contribution to the force, and the i=1 term has ple the slow and van Kampen modes; (ii) it is based on the been cancelled by a ·R. (cid:63) approximation that the disc mass µ → 0, and in this limit Modes having azimuthal wavenumber m = 1 can only thepatternspeedoftheslowmodegoestozerosotheouter be excited by those i ≥ 3 terms of V that produce cosφ e (cid:13)c 0000RAS,MNRAS000,000–000