Density Waves and the Viscous Overstability in Saturn's Rings
DDraft version October 19, 2018
Typeset using L A TEX default style in AASTeX62
Density Waves and the Viscous Overstability in Saturn’s Rings
Marius Lehmann, J¨urgen Schmidt, and Heikki Salo Astronomy Research Unit, University of Oulu, Finland
Submitted to A&AABSTRACTThis paper addresses resonantly forced spiral density waves in a dense planetary ring which is closeto the threshold for viscous overstability. We solve numerically the hydrodynamical equations fora dense thin disk in the vicinity of an inner Lindblad resonance with a perturbing satellite. Ournumerical scheme is one-dimensional so that the spiral shape of a density wave is taken into accountthrough a suitable approximation of the advective terms arising from the fluid orbital motion. Thispaper is a first attempt to model the co-existence of resonantly forced density waves and short-scalefree overstable wavetrains as observed in Saturn’s rings, by conducting large-scale hydrodynamicalintegrations. These integrations reveal that the two wave types undergo complex interactions, nottaken into account in existing models for the damping of density waves. In particular it is foundthat, depending on the relative magnitude of both wave types, the presence of viscous overstabilitycan lead to a damping of an unstable density wave and vice versa. The damping of the short-scaleviscous overstability by a density wave is investigated further by employing a simplified model of anaxisymmetric ring perturbed by a nearby Lindblad resonance. A linear hydrodynamic stability analysisas well as local N-body simulations of this model system are performed and support the results of ourlarge-scale hydrodynamical integrations.
Keywords: planets and satellites: rings, hydrodynamics, waves, instabilities INTRODUCTIONThe Cassini mission to Saturn has revealed a vast abundance of structures in the planet’s ring system, spanning a widerange of length scales. The finest of these structures have been detected by several Cassini instruments (Colwell et al.(2007); Thomson et al. (2007); Hedman et al. (2014)) and are periodic and quasi-axisymmetric with wavelengthsof some 100 m. It is generally accepted that this periodic micro structure originates from the viscous overstabilitymechanism which has been studied so far only in terms of axisymmetric models (Schmit and Tscharnuter (1995, 1999);Spahn et al. (2000); Salo et al. (2001); Schmidt and Salo (2003); Latter and Ogilvie (2008, 2009, 2010); Rein and Latter(2013); Lehmann et al. (2017)). On much greater scales, typically 10’s to 100’s of kilometers, numerous spiral densitywaves propagate through the rings, as these are excited at radii where the orbiting ring particles are in resonance withthe gravitational perturbation of one of the moons orbiting the ring system.The process of excitation and damping of resonantly forced density waves has been thoroughly studied, mostly interms of hydrodynamic models (Goldreich and Tremaine 1978a,b, 1979; Shu 1984; Shu et al. 1985; Shu et al. 1985;Borderies et al. 1985, 1986; Lehmann et al. 2016). Throughout the literature one typically distinguishes between linear and nonlinear density waves. The former are the ring’s response to a relatively small, resonantly perturbing force inthe sense that the excited surface mass density perturbation is small compared to the equilibrium value. In this casethe governing hydrodynamic equations can be linearized and as a consequence the density wave appears sinusoidal inshape. Corresponding author: Marius Lehmannmarius.lehmann@oulu.fi Upper limits for the cant-angle determined for these structures are within 1-3 degrees. a r X i v : . [ a s t r o - ph . E P ] O c t Lehmann et al.
The studies by Shu et al. (1985), Borderies et al. (1986) (BGT86 henceforth) and Lehmann et al. (2016) (LSS2016henceforth) considered the damping behavior of nonlinear density waves in a dense planetary ring, such as Saturn’sB ring. For a nonlinear density wave the surface density perturbation is of the same order of magnitude as theequilibrium value. Within a fluid description of the ring dynamics, the damping of a density wave is governed bydifferent components of the pressure tensor. The model by Shu et al. (1985) computes the pressure tensor from thekinetic second order moment equations, using a Krook-collision term. The model predicts reasonable damping lengthsof a density wave for assumed ground state optical depths (or surface mass densities) that do not exceed a certaincritical value (which depends on the details of the collision term). For optical depths larger than this critical value,the wave damping becomes very weak so that the resulting wavetrains propagate with ever increasing amplitude andnonlinearity. That said, the model fails to describe the damping of nonlinear waves in dense ring regions with highmutual collision frequencies of the ring particles, such as the wave excited at the 2:1 inner Lindblad resonance (ILR)with the moon Janus, propagating in Saturn’s B ring. The main reason for this behavior of the model at large collisionfrequencies is most likely the neglect of nonlocal contributions to the (angular) momentum transport (Shukhman(1984); Araki and Tremaine (1986)) in their kinetic model. On the other hand, BGT86 compute the pressure tensorfrom a fluid model (Borderies et al. (1985)), as well as by using empirical formulae, which yield the correct qualitativebehavior of the pressure tensor in a dense ring with a large volume filling factor. The computed damping lengths foroptical depths relevant to Saturn’s dense rings are fairly long and the authors suspect this to be a consequence of thefluid approximation.Borderies et al. (1985) have shown that density waves are unstable in a sufficiently dense ring (such as Saturn’sB ring), whereas they are stable in dilute rings of small optical depth. Schmidt et al. (2016) pointed out that theinstability condition of spiral density waves is identical to the criterion for spontaneous viscous overstability (Schmitand Tscharnuter (1995)) in the limit of long wavelengths. In LSS2016 we derived the damping of nonlinear densitywaves from a different view point compared to the approaches by BGT86 and Shu et al. (1985), which are based onthe streamline formalism (see Longaretti and Borderies (1991)). We considered the density wave as a pattern thatforms in response to this instability, using techniques that are widely applied in the studies of pattern formation insystems outside of equilibrium (Cross and Hohenberg (1993)). Consequently, the wave damping is described in termsof a nonlinear amplitude equation. The resulting damping behavior is very similar to what is predicted by the BGT86model.While the models by BGT86 and LSS2016 can predict steady state profiles of density waves alone in an overstable ringregion (see also Stewart (2016)), they do not take into account the possible presence of additional wave structures thatcan spontaneously arise in response to the viscous overstability, independent of a perturbing satellite. A first attemptto study the presence of multiple modes in a narrow ring within the streamline formalism was due to Longaretti (1989),but further improvements are required to model the (nonlinear) interaction of different modes. The possibility of co-existence of resonant spiral density waves and short-scale near-axisymmetric periodic micro structure was discovered byanalyzing stellar occultations of Saturn’s A ring, recorded with the Cassini Visual and Infrared Mapping Spectrometer(Hedman et al. (2014)). This paper is concerned with a modeling of this co-existence and a qualitative understanding ofinteractions between a resonantly forced density wave and the short-scale waves generated by the viscous overstability.In our one-dimensional hydrodynamical scheme we need to assume that both the density wave and the short-scalewaves are non-axisymmetric with the same azimuthal periodicity. However, since the short-scale waves resulting fromspontaneous viscous overstability have wavelengths of some 100 m (implying very small cant-angles of 10 − − − degrees), their dynamical evolution is expected to be very similar to that of the extensively studied axisymmetric modes(see the aforementioned papers). Hydrodynamical integrations presented in this paper confirm this expectation.In Section 2 we outline the basic hydrodynamic model equations. Section 3 explains the numerical scheme applied toperform large-scale integrations of the hydrodynamical equations. Sections 4, 5 and 6 discuss specific terms appearingin these equations that arise from the forcing by the satellite, the advection due to orbital motion of the ring fluid, aswell as the collective self-gravity forces, respectively. Results of large-scale hydrodynamical integrations are presentedin Section 7. Here we first describe the excitation process of a density wave as it follows from our integrations.Subsequently we test our scheme against the nonlinear models by BGT86 and LSS2016 in a marginally stable ring. Inaddition, we present some illustrative examples of density waves which propagate through a ring region which containssharp radial gradients in the background surface mass density. We then consider waves that propagate in an overstablering. In order to facilitate an interpretation of the results from our large-scale integrations, we introduce a simplifiedaxisymmetric model to describe the perturbation of a ring due to a nearby ILR. We perform a linear hydrodynamicstability analysis of this model to compute linear growth rates of axisymmetric overstable waves in the perturbed ring.By employing the same model we then perform local N-body simulations of viscous overstability in a perturbed ring.Finally, Section 8 provides a discussion of the main results. HYDRODYNAMIC MODELFrom the vertically integrated isothermal balance equations for a dense planetary ring we derive the model equations(Stewart et al. (1984); Schmidt et al. (2009) ∂ t τ = − [Ω − Ω L ] ∂ θ τ − u∂ r τ − τ ∂ r u,∂ t u = − [Ω − Ω L ] ∂ θ u − u∂ r u + 2Ω v − ∂ r (cid:2) φ d + φ s (cid:3) − σ ∂ r ˆ P rr ,∂ t v = − [Ω − Ω L ] ∂ θ v − u∂ r v −
12 Ω u − r ∂ θ φ s − σ ∂ r ˆ P rθ , (1)in a cylindrical frame ( r, θ, z = 0) with origin at r = r L , rotating rigidly with angular frequency Ω L = Ω( r L ) where r L denotes the radial location of a specific inner Lindblad resonance (ILR) with a perturbing satellite andΩ = (cid:20) GM P r (cid:21) / (2)with Saturn’s mass M P = 5 . · kg and the gravitational constant G = 6 . · − . In what follows we will alsomake use of the radial distance x = r − r L (3)as well as its scaled version ˜ x = x/r L .The quantity σ is the rings’ surface mass density and τ = σ/σ with the ground state surface mass density σ . Thesymbols u , v stand for the radial and azimuthal components of the velocity on top of the orbital velocity [Ω − Ω L ] r inthe rigidly rotating frame. Furthermore, ˆ P is the pressure tensor (see below). The central planet is assumed sphericalso that Ω = κ , the latter denoting the epicyclic frequency of ring particles. The rings’ ground state which describesthe balance of central gravity and centrifugal force is subtracted from above equations and we neglect the large-scaleviscous evolution of the rings which occurs on time scales much longer than those considered in this study.We neglect curvature terms containing factors 1 /r since these scale as λ/r ∼ − compared to radial derivatives.Here λ denotes the typical radial wavelength of a spiral density wave near its related Lindblad resonance where ˜ x (cid:28) θ we retain only the advective terms arising from the Keplerianmotion, i.e. the first terms on the right hand sides of Equations (1). All other θ -derivatives scale as ( mλ ) /r comparedto radial derivatives ( m denoting the number of spiral arms of the density wave), i.e. the same as curvature terms.Poisson’s equation for a thin disk ( ∂ r + ∂ z ) φ d = 4 πGσδ ( z ) , (4)establishes a relation between the self-gravity potential φ d and the surface density σ .The viscous stress is assumed to be of Newtonian form such that in the cylindrical frame we can writeˆ P = ˆ P rr ˆ P rθ ˆ P θr ˆ P θθ = (cid:32) p − η (cid:0) + ˆ γ (cid:1) ∂ r u − η (cid:0) − Ω + ∂ r v (cid:1) − η (cid:0) − Ω + ∂ r v (cid:1) p + η (cid:0) − ˆ γ (cid:1) ∂ r u (cid:33) . (5)It is thus completely described by radial gradients of the velocities u , v , the dynamic shear viscosity η as well as theisotropic pressure p (see below). The ratio of the bulk and shear viscosity is denoted by ˆ γ , which is assumed to be Lehmann et al. constant (Schmit and Tscharnuter (1995)). The isotropic pressure and the dynamic shear viscosity take the simpleforms p = p τ p σ , (6) η = ν σ τ β +1 . (7)In this study we assume p σ = 1, i.e. the equation of state for an ideal gas. The ground state pressure can be definedin terms of an effective ground state velocity dispersion c such that (Schmidt et al. (2001)) p = σ c . (8)The ground state is characterized by σ = const., u = 0, v = 0, together with the parameters in (6) and (7).We neglect azimuthal contributions due to collective self-gravity forces. This neglect is adequate as long as theexerted satellite torque is much smaller than the unperturbed viscous angular momentum luminosity of the ring. Thatis, the (self-gravitational) angular momentum luminosity carried by the wave is negligible compared to the viscousluminosity. The linear inviscid satellite torque deposited at the resonance site reads (Goldreich and Tremaine (1979)) T s = − mπ σ D Ω L ( (cid:15) r L Ω L ) [ ∂ r φ s − m φ s ] r L (9)where (cid:15) = 2 πGσ r L D (10)and (Cuzzi et al. (1984)) D = 3 ( m −
1) Ω L . (11)The viscous angular momentum luminosity in the unperturbed disk is given by (Lynden-Bell and Pringle (1974)) L visc = 3 πν σ Ω r . In addition it should be mentioned that we are not concerned with the long-term redistribution of ring surface massdensity which occurs in response to the presence of very strong density waves (BGT86) so that we assume σ = const as mentioned before.For the sake of definiteness we will restrict to parameters corresponding to the Prometheus 7:6 ILR, located at r ∼ ,
000 km in Saturn’s A ring. We take values of the rings’ ground state shear viscosity ν and surface massdensity σ (see Table 1) that can be estimated from corresponding values obtained by Tiscareno et al. (2007) forthis ring region. The nominal values of β and γ correspond to values found in N-body simulations with an opticaldepth τ dyn = 1 [see LSS2016 (Section 3)]. Besides the nominal values we will use a range of values for β [Equation(7)] and also T s [Equation (9)], in order to explore a variety of qualitatively different scenarios for the damping ofdensity waves. The adopted value for the ground state velocity dispersion c is larger then what results from localnon-gravitating N-body simulations for optical depths relevant to this study (e.g. Salo (1991)) but corresponds roughlyto expected values for Saturn’s A ring from self-gravitating N-body simulations exhibiting gravitational wakes andassuming meter-sized particles (Daisaka et al. (2001); Salo et al. (2018)). Furthermore, the value is still small enoughto ignore pressure effects on the density waves’ dispersion relation (Section 7.1).Our hydrodynamic model exhibits spontaneous viscous overstability on finite wavelengths if the viscous parameter β exceeds a critical value. To see this, let us ignore the satellite forcing φ s for the time being. We restrict to shortradial length scales so that Ω = Ω L can be considered constant except that we use[Ω − Ω L ] ∂ θ = −
32 Ω L x∂ y in Equations (1). Our 1D numerical method to solve Equations (1) assumes that any mode which forms has m -foldazimuthal periodicity (see Section 5). Hence let us introduce non-axisymmetric oscillatory perturbations such that τuv = + τ (cid:48) u (cid:48) v (cid:48) exp (cid:26) ωt + i (cid:18) k x + 32 m − m k y Ω L t (cid:19) x + ik y y (cid:27) , (12)with complex oscillation frequency ω = ω R + i ω I and real-valued radial and azimuthal wavenumbers k x > k y = mr L , respectively. The time-dependent contribution to the radial wavenumber in (12) stems from the winding ofthe perturbations due to Keplerian shear [see Meyer-Vernet and Sicardy (1987) and Equation (46)]. Since we knowthat the linear growth and the nonlinear saturation of spontaneous viscous overstability occurs on wavelengths oftypically hundreds of meters it turns out that we can neglect the effect of winding in (12). That is, for the relevantmodes the time it takes for the winding term to become equal to k x is given by t ∼ m r L λ x ORB . With λ x = 2 π/k x ∼
100 m this yields some 100,000 orbits, which is much longer than the time scale of the nonlinearevolution of the modes (i.e. thousands of orbits, see Latter and Ogilvie (2010); Rein and Latter (2013); LSS2017).Furthermore, Poisson’s equation (4) yields the relationship φ (cid:48) d = − πGσ k x τ (cid:48) (13)for a single wavelength mode (Binney and Tremaine (1987)).In the remainder of this section we apply dimensional scalings such that time is scaled with 1 / Ω L and length is scaledwith c / Ω L . Inserting (12) and (13) into (1) and linearizing with respect to the perturbations (the primed quantities),results in the eigenvalue problem 0 = − ω + ω (cid:20) − (cid:18)
73 + γ (cid:19) k x ν + 92 im ˜ x (cid:21) + ω (cid:20) − gk x − k x − (cid:18)
43 + γ (cid:19) k x ν + i (7 + 3 γ ) k x m ˜ xν + 274 m ˜ x (cid:21) − k x (cid:0) β − gk x + k x (cid:1) ν + 12 im ˜ x (cid:0) − gk x + 3 k x + (4 + 3 γ ) k x ν (cid:1) + 34 (7 + 3 γ ) k x m ˜ x ν − im ˜ x (14)for ω = ω R + iω I . The non-dimensional distance ˜ x is defined as below Equation (3). This equation can be used toobtain the growth rate ω R ( k x ) and oscillation frequency ω I ( k x ) of a given mode k x . This procedure has been carriedout for axisymmetric modes (with m = 0) in several papers [see Lehmann et al. (2017) ( LSS2017 hereafter) andreferences therein for more details]. It can be shown that the growth rates ω R ( k x ) following from Equation (14) areindependent of m (i.e. independent of k y ) and agree with those of previous studies.In the remainder of the paper the symbol k denotes the radial wavenumber of a given mode. The threshold forviscous oscillatory overstability, i.e. a vanishing growth rate ω R ( k ) = 0, can be obtained by setting ω = iω I and solvingthe imaginary and real parts of Equation (14) for ω I and β , respectively, for a given wavenumber k . This yields thecritical frequency pair ω c ( k ) = 32 m ˜ x ± (cid:115) − gk + k + (cid:18)
43 + γ (cid:19) ν k (15) The third critical frequency is associated with the diffusive viscous instability, not considered in this paper.
Lehmann et al. and the critical value of the viscosity parameter β c ( k ) = 13 (cid:18) γ − (cid:19) − (cid:18)
43 + γ (cid:19) gk + 13 (cid:18)
43 + γ (cid:19) k + 127 (cid:0)
28 + 33 γ + 9 γ (cid:1) ν k , (16)which describes the stability boundary for viscous overstability and which is also independent of m . The frequencies(15) are Doppler-shifted by m ˜ x as compared to the frequencies of axisymmetric modes. Note that due to the fact thatEquations (1) are defined in a frame rotating with Ω L this Doppler-shift is very small as ˜ x ∼ − for all cases consideredin this paper. The Doppler-shift can therefore be neglected. These results show that linear free non-axisymmetricshort-scale modes due to spontaneous viscous overstability in our hydrodynamic model behave essentially the same asaxisymmetric modes with m = 0.The curve β c ( k ) possesses a minimum at finite wavelength if g >
0, i.e. for a non-vanishing collective self-gravityforce. This wavelength is roughly two times the Jeans-wavelength λ J = c / ( Gσ ). In the above equations we define g = πGσ Ω c , (17)denoting the inverse of the hydrodynamic Toomre-parameter (a full list of symbols is provided in Table 2). Table 1 : Hydrodynamic Parameters parameter Prometheus 7:6 (
P r c [10 − ms − ] 1.5 ν [10 − m s − ] 1 γ β σ [kg m − ] 350 r L [10 m] 1.26 T s [10 kg m ] 4.56 L visc [10 kg m ] 72003. NUMERICAL METHODSFor numerical solution of Equations (1) we apply a finite difference Flux Vector Splitting method employing aWeighted Essentially Non-Oscillatory (WENO) reconstruction of the flux vector components. The method is identicalto that used in LSS2017, apart from the reconstruction of the flux vector.We define the flux-conservative variables U = ττ uτ v so that Equations (1) can be written as ∂ t U = − ∂ r F + S (18)with the flux vector F = τ uτ u + τ c τ uv and the source term S = − [Ω − Ω L ] ∂ θ τ − [Ω − Ω L ] ∂ θ ( uτ ) + 2Ω τ v − τ ∂ r ( φ d + φ s ) + σ ∂ r ˆΠ rr − [Ω − Ω L ] ∂ θ ( vτ ) − Ω τ u + σ ∂ r ˆΠ rθ . (19)In the last expression ˆΠ = p ˆ U − ˆ P is the viscous stress tensor with ˆ U denoting the unity tensor.We solve (18) on a radial domain of size L r . The domain is discretized by defining nodes r j ( j = 1 , , . . . , n ) withconstant inter-spacing h = r j +1 − r j . We adopt periodic boundary conditions in all integrations. Since a density waveis not periodic in radial direction this requires the radial domain size L r to be large enough so that the Lindbladresonance is located sufficiently far from the inner domain boundary and that an excited density wave is fully dampedbefore reaching the outer domain boundary. The discretization of the flux derivative ∂ r F is outlined in Appendix E.The source term (19) contains radial derivatives of the stress tensor which are evaluated with central discretizations of12th order. Furthermore, the evaluation of the derivatives with respect to θ and the self-gravity force ∂ r φ d appearingin (19) will be discussed in Sections 5 and 6, respectively.Due to the presence of the satellite forcing terms in (19) it turns out that the simple time step criterion arisingfrom a one-dimensional advection-diffusion problem, which was used in LSS2017, is unnecessarily strict. This criterionreads ∆ t ≤ min (cid:18) h ν , ν Λ (cid:19) , (20)where Λ is identified with the maximal eigenvalue of the Jacobianˆ A = ∂ F ( U ) ∂ U (21)of Equations (18) for the whole grid and ˆ ν is to be identified with the maximal value of the coefficient in front of theterm containing the second radial derivative ∂ r u in (1), which isˆ ν = ν (cid:18)
43 + γ (cid:19) τ β . The three eigenvalues of (21) read Λ = u, Λ (2 / = u ± c . For most integrations presented in this paper the grid spacings h are large enough so that the second term in (20)is by far the smallest and can take values down to some 10 − ORB. We find, however, that time steps in the range∆ t = 1 − · − ORB are suitable for all presented integrations, indicating that the criterion (20) cannot be appropriate.We have checked for some integrations with strong satellite forcing that reducing the time step by a factor of 0 . e kin = 1 L r (cid:90) [ L r ] d r σ (cid:0) u + v (cid:1) . (22) Lehmann et al.
Table 2 : List of Symbols
Symbol Meaning G gravitational constant M P planet’s mass M s mass of perturbing satellite a s semimajor axis of perturbing satellite e s eccentricity of perturbing satellite T s linear inviscid satellite torque φ s satellite potential L visc viscous angular momentum luminosityΩ Kepler frequencyΩ L Kepler frequency at r = r L κ epicyclic frequencyΩ Z vertical frequency of ring particles x = r − r L radial coordinate˜ x = r − r L r L scaled radial coordinate t time ω s satellite forcing frequency in the frame rotating with Ω L ˆ ω s satellite forcing frequency in the inertial frame ω complex frequency of overstable waves σ surface mass density τ = σσ scaled surface mass density τ dyn dynamical optical depth u , v planar velocity components φ d self-gravity potential φ p planetary potential c = (cid:113) p σ effective isothermal velocity dispersion ν ground state kinematic shear viscosity p isotropic pressure η dynamic shear viscosity β viscosity parameter γ constant ratio of bulk and shear viscosityˆ P pressure tensor g = πGσ Ω c inverse ground state Toomre-parameterˆ γ phase variable of a fluid streamline q nonlinearity parameter of a fluid streamline∆ phase angle of a fluid streamline a semimajor axis of a fluid streamline e eccentricity of a fluid streamline4. SATELLITE FORCING TERMSFor simplicity, we restrict to density waves that correspond to a particular inner Lindblad resonance of first order,so that the forcing satellite orbits exterior to the considered ring portion. The wave is excited by a particular Fouriermode of the gravitational potential due to this satellite with mass M s and semi-major axis a s and reads (cf. Section In the current approximation a Lindblad resonance coincides with a mean motion resonance. φ s ( t, ˆ θ ) = − GM s a s b m / exp (cid:110) i (cid:16) m ˆ θ − ˆ ω s t (cid:17)(cid:111) , valid in an inertial frame denoted by ( r, ˆ θ ). The symbol b m / = 2 π (cid:90) π dΨ cos ( m Ψ) (cid:112) ρ − ρ cos (Ψ)is a Laplace-coefficient with ρ = ra s . In the current approximation the forcing frequency readsˆ ω s = m Ω s = ( m −
1) Ω L (23)with the satellite mean motion Ω s . Upon changing to the frame rotating with frequency Ω L , denoted by ( r, θ ), wehave ( m ˆ θ − ˆ ω s t ) → ( m [ θ + Ω L t ] − ˆ ω s t ) ≡ mθ − ω s t, yielding the forcing frequency in the rotating frame ω s = ˆ ω s − m Ω L = − Ω L , (24)where we used (23). Therefore, the radial forcing component appearing in Equations (1) reads − ∂ r φ s = GM s a s ∂ r b m / exp { i ( mθ + Ω L t ) } . (25)Similarly, the azimuthal component is given by − r ∂ θ φ s = GM s a s b m / m exp { i ( mθ + Ω L t + π/ } . (26)These terms are evaluated at r = r L . AZIMUTHAL DERIVATIVESThe persistent spiral shape of a (long) density wave is generated by the resonant interaction between the ringmaterial and the perturbing satellite potential, as well as the collective self-gravity force. Since our integrations areone-dimensional, it is not possible to describe azimuthal structures directly. Therefore we need to restrict Equations (1)to a radial cut which we choose to be θ = 0 without loss of generality. The information about the azimuthal structure ofthe density pattern (the number of spiral arms m ) is then contained solely in the terms describing azimuthal advectiondue to orbital motion. i.e. the first terms on the right hand sides of Equations (1). In what follows we refer tothese terms simply as “azimuthal derivatives“. Thus, the requirement is to prescribe proper values of the azimuthalderivatives at θ = 0 for each time step of the integration.We again adopt the cylindrical coordinate frame ( r, θ ) of Section 2 which rotates with angular velocity Ω L relativeto an inertial frame denoted by ( r, ˆ θ ) so that θ = ˆ θ − Ω L t. Lehmann et al.
If we linearize Equations (1) with respect to the variables τ , u , v and φ d , so that we restrict these to describe lineardensity waves, it is possible to solve the equations in the complex plane by splitting the solution vector Ψ into its realand imaginary parts: Ψ = Ψ R + i Ψ I = τ ( r, θ, t ) u ( r, θ, t ) v ( r, θ, t ) = τ R ( r, θ, t ) u R ( r, θ, t ) v R ( r, θ, t ) + i τ I ( r, θ, t ) u I ( r, θ, t ) v I ( r, θ, t ) . (27)An evolving m -armed linear density wave can then be described through the complex vector of state τ ( r, θ, t ) u ( r, θ, t ) v ( r, θ, t ) = A τ ( r, t ) A u ( r, t ) A v ( r, t ) exp { i ( mθ − ω s t ) } (28)with A τ ( r, t ), A u ( r, t ) and A v ( r, t ) being complex amplitudes in this notation which depend on time and the radialcoordinate. The time dependence of the amplitudes is generally much slower than the oscillatory terms and vanishesonce the integration reaches a steady state. When inserted to the linearized Eqs. (1) we obtain two sets of threeequations that are possibly coupled through the azimuthal derivatives and self-gravity, depending on the appliedimplementation (cf. Appendices D and F). Note that in practice we will exclusively use the full nonlinear equations(1). For sufficiently small amplitudes A τ ( r, t ), A u ( r, t ) , A v ( r, t ) the nonlinear terms in (1) are negligible and theequations are essentially linear.In order to describe nonlinear density waves , it is necessary to make an approximation for the azimuthal dependencyof the wave. To obtain such an approximation we assume that (28) holds also in the nonlinear case. We will discussthe validity of this assumption a bit more at the end of Section 5.1. We have found two such implementations for theazimuthal derivatives (simply referred to as Methods A and B ) that yield a stationary final state for Equations (1) inthe nonlinear regime. It turns out that the application of Method A (Section 5.1) results in nonlinear wave profilesthat agree better with existing nonlinear models and therefore the results presented in subsequent sections are basedon integrations using this method. We additionally outline
Method B (Appendix D) as we have found it to work wellin the weakly nonlinear regime. Note that for sufficiently linear waves, both methods are exact down to the numericalerror. 5.1.
Method A
One implementation of the azimuthal derivatives can be derived if one considers the vector of state of the weaklynonlinear model of LSS2016 in the first order approximation, where contributions from higher wave harmonics are1omitted. That is, we have ˜ φ d ( r, θ, t )˜ u ( r, θ, t )˜ v ( r, θ, t ) = − i D ˜Ω˜ D ˜ x + (cid:16) ˜ ω s − m [ ˜Ω − ˜Ω L ] (cid:17) i (cid:16) ˜ ω s − m [ ˜Ω − ˜Ω L ] (cid:17) ˜ D ˜ x + (cid:16) ˜ ω s − m [ ˜Ω − ˜Ω L ] (cid:17) × A ( r ) exp { i ( Φ( r ) + mθ − ω s t ) } + c.c. , (29)where ˜ φ d , ˜ u , ˜ v , as well as ˜Ω, ˜Ω L , ˜ D and ˜ ω s are understood to be scaled according to Table 1 in LSS2016 and ˜ x isscaled according to Table 2. Note that ˜ ω s is the scaled version of (24). Equation (29) corresponds to Equations (35)and (45) in LSS2016, except that (29) is written in the rotating frame and is not expanded to the lowest order in ˜ x ,although small corrections due to pressure and viscosity are neglected.From the solution of the Poisson-Equation we have (Shu (1984)) τ = i sgn( k )2 πGσ ∂ r φ d (30)where φ d is assumed to be given by the unscaled first component of (29) and k = ∂ r Φ( r ) denotes the wavenumber ofthe density wave. Note that both sides of Equation (30) are understood to be real-valued since the sgn-function takesdifferent signs for the two conjugate complex exponential modes in (29). The azimuthal derivative of τ is then givenby ∂ θ τ = − m πGσ ∂ r φ d (31)if we assume k >
0. Equation (31) shows that the azimuthal derivative of τ is directly proportional to the radialcomponent of the self-gravity force, the computation of which we discuss in Section 6. The azimuthal derivatives of u and v can be directly obtained from (29) and read ∂ θ u = − m Ω ( ω s − m [Ω − Ω L ]) D ˜ x + ( ω s − m [Ω − Ω L ]) v ≈ mv, (32) ∂ θ v = m D ˜ x + ( ω s − m [Ω − Ω L ])
2Ω ( ω s − m [Ω − Ω L ]) u ≈ − mu. (33)The approximate expressions follow if we neglect D ˜ x , which is fully justified since ˜ x ∼ − for all cases consideredhere.As mentioned before, Equation (28) can only be used as an approximation for a nonlinear density wave. Assumethat the latter is correctly described by an infinite series τ ( r, θ, t ) u ( r, θ, t ) v ( r, θ, t ) = ∞ (cid:88) l =1 A τ,l ( r, t ) A u,l ( r, t ) A v,l ( r, t ) exp { li ( mθ − ω s t ) } + c.c. , (34)where the terms with l > Lehmann et al. the errors of the azimuthal derivatives themselves. Let us for the time being assume that the surface density τ in a(steady state) density wave can be described through [Borderies et al. (1983), see also Section 7.4.3 and Appendix G] τ ( r, φ ) = 11 − q ( r ) cos ( mφ + Φ( r )) (35)in a cylindrical frame ( r, φ, z = 0) rotating with the satellite’s mean motion frequency Ω s = ˜ ω s /m [see Equation (23)].Furthermore, q is the nonlinearity parameter fulfilling 0 ≤ q < r ) is the radial phase function of the densitywave [cf. (29)]. Clearly, for q not much smaller than unity the variation of the surface density τ deviates significantlyfrom a simple harmonic. Taking the azimuthal derivative of (35) yields ∂ φ τ = − mq ( r ) sin ( mφ + Φ( r ))[1 − q ( r ) cos ( mφ + Φ( r ))] . (36)By expanding this expression to first order in q , which amounts to our linear treatment of the azimuthal derivative inEquation (31), we obtain ∂ φ τ = − mq ( r ) sin ( mφ + Φ( r )) + O (cid:0) q (cid:1) . (37)This would imply that for 0 . < q < . ∼ −
60 %.Despite the considerable error that may be induced by the approximation (31) [and (32), (33)] we will see in Section7.2 that the resulting error in the radial density wave profiles is actually small. As for the approximation (31) thereason is that in our integrations we evaluate this term by using an accurate (nonlinear) expression for the self-gravityforce ∂ r φ d (Section 6) so that the actual error resulting from (31) is much smaller than what would result from thelinearized expression (37). This can be understood by realizing that Equation (30) holds also for the higher harmonics[ l > τ and φ d (see Appendix B of LSS2016 for more details) which stems from the fact that Poisson’sequation is linear. Thus, the actual error which is then made with the approximation (31) is that the contributionsof the higher harmonics l > underestimated by factors of 1 /l , but not entirely neglected. On the otherhand, from LSS2016 (Section 4.5) follows that the approximations (32), (33) hold also for the second harmonics of thevelocity fields upto a factor 1 /l (with l = 2) and the same is expected to apply to all higher harmonics l ≥
3. Hence,the error made with (32) and (33) is also a suppression of the higher harmonics l ≥ /l .Finally, note that our approximations for the azimuthal derivatives (31), (32) and (33) imply that any mode whichforms during an integration on top of the equilibrium state will be non-axisymmetric with azimuthal periodicity m .For this reason the short-scale overstable waves which appear in our integrations (Section 7.4) are non-axisymmetricwith the same m as the resonantly forced density wave. As outlined in Section 2 it is expected that the dynamicalevolution of these modes is very similar to that of axisymmetric modes. This expectation will be confirmed in Section7.4.1. SELF-GRAVITYFor most integrations we use the same implementation of collective radial self-gravity forces as described in detailin LSS2017. The model approximates the ring material as a collection of infinite straight wires (neglect of curvature)and predicts a self-gravity force at grid point j : f dj = − Ghσ n (cid:88) i =1 ,i (cid:54) = j τ ( r i ) r j − r i | r j − r i | . (38)where we defined f d = − ∂ r φ d . This relation (38) does not include the force generated by mass contained in the bin j itself, which can be approximated through∆ f d (0) = 2 Gσ (cid:2) ∂ r τ (0) h + O (cid:0) h (cid:1)(cid:3) . The analysis in LSS2016 is restricted to second order harmonics. τ ( r i ) is periodic with period n the sum (38) can be replaced by the convolution f dj = n − (cid:88) i = − n τ i f kernj − i (39)of τ i = τ ( r i ) with the force kernel, which reads f kernj − i = − Ghσ r j − r i | r j − r i | . Equation (39) can then be solved with a FFT method. However, since the density pattern τ i of a resonantly forceddensity wave is not periodic we need to pad one half of the array τ i with zeros in order to avoid false contributionsfrom grid points outside the actual grid (Binney and Tremaine (1987)), e.g. gravitational coupling of material insidethe resonance with material at far positive distances from resonance across the boundaries.For comparison with existing models for nonlinear density waves in Section 7.2 we also perform integrations whichadopt the WKB-approximation for the radial self-gravity force. The corresponding implementation is described inAppendix F. RESULTS7.1.
Excitation of Density Waves
In this section we illustrate the excitation of a resonantly forced density wave as it results from our integrations. Weuse integrations employing the
P r L r = 450 km, ∆ t = 5 · − ORB, and h = 45 m. Furthermore, Method A for the azimuthal derivatives (Section 5.1) and the
Straight Wire self-gravity model (Section 6) were em-ployed. The initial state of each integration is the Keplerian shear flow with τ ( t = 0) = 1, u ( t = 0) = 0, v ( t = 0) = 0(Section 2) and the satellite forcing is introduced at time t = 0.It is expected that during the excitation process the envelope of a density wave evolves in radial direction with thelocal group velocity (Toomre (1969); Shu (1984)). For a linear density wave described by the perturbed surface density σ = σ + A ( r ) · exp (cid:40) i (cid:90) r k ( s ) d s (cid:41) · exp { i ( mθ − ω s t ) } ] (40)with wavenumber k , one obtains in the frame rotating with frequency Ω L the dispersion relation (Goldreich andTremaine (1978a); Shu (1984)) κ − ( ω s − m [Ω − Ω L ]) + k c − πGσ | k | = 0 . (41)Taking the derivative with respect to k on both sides and re-arranging terms yields the group velocity (Toomre (1969)) v g = d ω s d k = sgn( k ) πGσ − | k | c m [Ω − Ω L ] − ω s , (42)where ω s is given by (24) and sgn( k ) denotes the sign of k . By defining D = κ − ( ω s − m [Ω − Ω L ]) and expanding this expression about the Lindblad resonance r = r L (using the approximation κ = Ω) so that D = D ˜ x + O (cid:0) x (cid:1) with D given by (11), one obtains from (41) the wavenumber dispersion for linear density waves k = ˜ x(cid:15)r L , (43)where (cid:15) is given by (10) and where the effects of pressure, expressed through the term quadratic in k in (41), areignored.4 Lehmann et al.
An expression for the nonlinear group velocity can be obtained from the nonlinear dispersion relation of spiral densitywaves [e.g. Equation (87) of Shu et al. (1985)] in the WKB-approximation κ − ( ω s − m [Ω − Ω L ]) + I ( q ) k c − πGσ | k | H ( q ) = 0 . (44)In this expression the contributions due to pressure and self-gravity are modified and depend on the nonlinearityparameter q with 0 ≤ q < H ( q ) = 1 π ∞ (cid:90) −∞ d u sin uu I (cid:18) q sin uu (cid:19) describes the nonlinear effects of self-gravity (Shu et al. (1985)) and, similarly, I ( q ) = 2 q (cid:104)(cid:0) − q (cid:1) − / − (cid:105) , describes nonlinear pressure effects. The resulting group velocity reads v nlg = sgn( k ) πGσ (cid:104) H + 2 k qq (cid:48) ∂H∂q (cid:105) − | k | c (cid:104) I + k qq (cid:48) ∂I∂q (cid:105) m [Ω − Ω L ] − ω s , (45)where a prime stands for the derivative ∂/∂k . The integral functions H ( q ) and I ( q ) fulfill H ( q ) ≥ I ( q ) ≥ q ≥
0. Furthermore, it can be verified that all other quantities enclosed in the brackets are real-valued and positive.In the linear limit q → v nlg is identical to (42), as expected. The wavenumber of thedensity wave k and the nonlinearity parameter q depend on the radial distance from resonance and for a tightly wounddensity wave we have q ∝ k (Shu et al. (1985); BGT86, see also Section 7.4.3). For typical values of the velocitydispersion c in Saturn’s dense rings the self-gravity term in (45) will always dominate the pressure term so that thenonlinear group velocity is expected to be larger than the linear limit (42).Figures 1, 2 and 3 show stroboscopic space-time diagrams [time-resolution of Ω L / (2 π )] of integrations with scaledlinear satellite torques of ˜ T s = 10 − , ˜ T s = 1 and ˜ T s = 4, where ˜ T s = 1 corresponds to the nominal forcing strength forthe Prometheus 7:6 ILR (Table 1). In these figures the gray shading measures the value of τ so that brighter regionscorrespond to larger values of τ . Since at t = 0 the spatially constant satellite forcing is introduced and the disk ishomogeneous, initially the hydrodynamic quantities u , v and τ oscillate uniformly (with infinite wavelength). Due toKeplerian shear the pattern starts to wrap up at a constant rate. This transient behavior was derived by Meyer-Vernetand Sicardy (1987) who studied the interaction of a satellite with an initially homogeneous disk in the vicinity of aLindblad resonance and in the linear limit. They showed that the wavelength of the pattern evolves as λ p ( t ) = 4 πr L m −
1) Ω L t . (46)This result was obtained in the absence of collective forces. Meyer-Vernet and Sicardy (1987) argued that aftersufficiently long time the transient behavior vanishes and the system settles on a stationary solution which is governedby collective effects (self-gravity, pressure and viscosity). They proved this for the case of a simple friction law assuminga force f = − Q u with u = ( u, v ) in the momentum equation. In the present situation self-gravity is the dominantcollective force and the disk excites a long trailing density wave propagating outward from the ILR with group velocityapproximately given by (42) (Goldreich and Tremaine (1978a); Shu (1984)).As the wavelength of the pattern decreases with time, at a certain radial location and at a certain time the wavelengthwill fulfill the dispersion relation (44) [and also (41) if q is sufficiently small]. As soon as this is the case, the wavelengthis “locked” to this value. In the figures 1-3, the region which becomes “locked“ is enclosed by the dashed and solid bluelines. The former marks the resonance, while the latter is the predicted path of the wave front assuming it propagateswith the linear group velocity (42). All wave structures outside this region eventually damp as they become increasinglywound up. An exception are the short-scale waves generated by viscous overstability (Sections 2 and 7.4). Also plottedare radial profiles of τ at four different times during the excitation process.In Figure 1 the blue solid line describes well the propagation of the wave front, until a steady state is reached(around 8,000 orbital periods) and the wave’s amplitude profile remains stationary. We find a number of differences5when comparing the figures. First of all, with increasing torque value the wave profiles attain the typical peakyappearance of nonlinear density wavetrains in thin disks (Shu et al. (1985); BGT86; Salo et al. (2001)). Furthermore,the group velocity of the waves increasingly departs from the linear prediction (42), albeit mildly. One notes that thereremains a very slow phase-drift of the wave pattern towards the resonance, indicating an increasing phase velocity withdecreasing distance from resonance. Theoretically, at resonance the wavenumber of the density wave (43) vanishes sothat the wave’s phase velocity ω s /k diverges. It can therefore in general not be expected from a numerical method tocorrectly describe the wave pattern at the exact resonance location.Figure 4 shows for the integration with ˜ T s = 10 − (Figure 1) the average wavelength of the forming pattern, sampledwithin the radial region 50 km ≤ r − r L ≤
150 km. The agreement with Equation (46) is excellent for about 3,000orbital periods. After that deviations become notable as a steady state is reached where self-gravity prevents furthershortening of the wavelength. The closer to the resonance r = r L , the earlier a steady state is attained as the resonantdensity wave pattern emerges at the resonance and propagates outward with its local group velocity.In Saturn’s rings an initial transient pattern as seen in our integrations might be observable for density waves drivenby the co-orbital satellites Janus and Epimetheus. These satellites interchange orbits every 4 years so that theirresonance locations in the rings shift periodically by tens of kilometers. Every time a resonance location is changedthe wave excited at the preceding location continues to propagate while a new density wave is launched at the newlocation (Tiscareno et al. (2006)). L [km]02000400060008000 t i m e [ O R B ] Figure 1 : Stroboscopic space-time diagram showing the evolution of the scaled surface density τ for an integrationwith the P r T s = 10 − . Brighter regions correspond to larger τ -values.The blue solid line marks the path of a signal traveling with the linear group velocity (42) starting from resonance r = r L at time t = 0. Also shown are profiles of τ at different stages of the evolution. Due to the plot beingstroboscopic, the density wave pattern eventually becomes stationary as the oscillation with frequency ω s = − Ω L iseffectively removed. 7.2. Comparison with the Nonlinear Models of BGT1986 and LSS2016
In this section we compare results of our hydrodynamical integrations with the nonlinear models of BGT86 andLSS2016, which we refer to as the BGT and the WNL (Weakly Nonlinear) model, respectively. This section isrestricted to stable waves in the sense that β < β c ( λ ) for all wavelengths λ [cf. Equation (16)], i.e. no overstabilityoccurs in the system. All hydrodynamical integrations presented in this section were conducted with L r = 450 km,time steps ∆ t = 5 · − ORB, and spatial resolution h = 45 m and used the P r
Method A for the azimuthal derivatives (Section 5.1) and the
Straight Wire Lehmann et al. L [km]02000400060008000 t i m e [ O R B ] Figure 2 : Same as Figure 1 except that ˜ T s = 1. L [km]02000400060008000 t i m e [ O R B ] Figure 3 : Same as Figure 1 except that ˜ T s = 4.self-gravity model (Section 6). Presented BGT model wave profiles were computed using the method of BGT86 (seetheir Section IVa), using the pressure tensor (5) with (6) and (7). This model takes into account secular changes inthe background surface mass density σ that accompany the steady state density wave in order to ensure conservationof angular momentum at all radii in the steady state. These modifications are neglected in our hydrodynamicalintegrations as well as the WNL model since the latter neglect the angular momentum luminosity carried by thedensity wave. To facilitate a comparison between the three different methods the profiles of τ resulting from the BGTmodel are scaled with the modified background surface mass density σ ( r ) which will not be shown.7
100 1000 10000time [ORB]110100 λ p [ k m ] pattern wavelengthEq (46) Figure 4 : Average wavelength of the forming pattern in the course of the integration shown in Figure 1. Thewavelength is obtained by counting the number of complete wave cycles of the (sinusoidal) surface density τ withinthe radial region 50 km ≤ r − r L ≤
150 km.Figure 16 (Appendix A) displays steady state profiles of the hydrodynamic quantities τ , u , v as these result fromintegrations together with profiles obtained using the WNL model (LSS2016). The profiles in the left and right columnsresult from integrations which applied Method A and
Method B for the azimuthal derivatives, respectively. The self-gravity is computed with the
Straight Wire model. As in the previous section, the satellite forcing strengths areindicated by the fractional torque ˜ T s such that ˜ T s = 1 corresponds to the nominal forcing strength at the Prometheus7:6 ILR and results in a nonlinear density wavetrain. The value ˜ T s = 9 · − corresponds to a weakly nonlinearwave. For the latter case both methods A and B produce very similar results in good agreement with the WNL model.Inspection of the nonlinear case ˜ T s = 1 reveals significant departures at larger distances from resonance between bothmethods, and Method A produces a clearly better match with the WNL model. All integrations presented in thefollowing sections were conducted with
Method A .In Figure 17 (Appendix A) we present wave profiles along with their Morlet wavelet powers (Torrence and Compo(1998)) for the case ˜ T s = 4. Also for this strongly nonlinear wave, we observe an overall good agreement for boththe amplitude profiles and wavenumber dispersions. Note that the WNL model takes into account only the first twoharmonics of the wave [cf. Equation (34)], which is clearly visible in the wavelet power. This is also the reason why τ can take values below 0.5 (see LSS2016 for details).Finally, Figure 18 (Appendix A) compares profiles obtained from integrations with the Straight Wire self-gravity (leftpanels) and the WKB self-gravity (Appendix F) using Equation (F9) (right panels) for the cases ˜ T s = 9 · − (upperpanels) and ˜ T s = 4 (lower panels). Comparison with corresponding BGT model wave profiles shows that the WKB-approximation is fully adequate for the weakly nonlinear wave with ˜ T s = 9 · − in that it yields indistinguishableresults from the Straight Wire self-gravity. For the strongly nonlinear case ˜ T s = 4 the WKB-approximation has anotable effect. As expected, its application yields overall an even better agreement with the BGT (and WNL) model.We have verified that reducing the time step or the grid spacings by factors of 1 / Wave Propagation through Density Structures
In this section we present a few illustrative examples of hydrodynamical integrations of density waves propagatingthrough an inhomogeneous ring. We restrict to the cases of jumps in the equilibrium surface density, but in principlewe could also vary other parameters with radial distance, such as the viscosity parameter β . All integrations adoptedthe P r
Method A for the azimuthal derivatives (Section 5.1) as well as the
Straight Wire self-gravity model (Section 6). Figure 19 (Appendix B) shows space-time plots of a density wave passing a region of8
Lehmann et al. increased surface density ( τ = 3, left panel) as well as a region of decreased surface density ( τ = 0 .
5, right panel),in both cases of radial width 40 km. The jumps in the equilibrium surface density, whose locations are revealed in thespace-time plots, act like additional sources for the density wave in the sense that the wave profile can change at theselocations prior to the expected arrival time of the wave front at these locations, the latter being indicated by the solidblue line (cf. Figures 1-3). It is, however, not clear how this is affected by the assumption imposed by our azimuthalderivatives that the hydrodynamic quantities describe an m -armed pattern right from the start of the integration.Figure 20 (Appendix B) shows steady state profiles of τ along with corresponding wavelet-power spectra of densitywaves passing through regions of varying equilibrium surface density. For reference, the first row shows a density wavein a homogeneous ring. The second case, with a region of increased τ , bears some similarities with Figures 4 and 5 inHedman and Nicholson (2016), showing profiles of the Mimas 5:2 density wave in Saturn’s B ring which passes througha region of radial width ∼
60 km where the normal optical depth increases sharply from about 1.5 to values 3 −
5. Inthe region of enhanced surface density in Figure 20 the wave damping is reduced due to its decreased wavenumber.Therefore, after passing the barrier the wave amplitude is enhanced as compared with the wave in the homogeneousregion. The last case represents a situation with a narrow region of mildly decreased surface density τ = 0 .
5. In thisregion the wavenumber is enhanced, resulting in stronger wave damping.If a density wave encounters a sharp discontinuity in the background surface density, such as a sharp ring edge, it istheoretically expected that it (partially) reflects at the boundary (see Longaretti (2018) and references therein). Forthe examples presented in Figure 20 the jumps in the background surface density are not sufficiently sharp to cause anotable reflection. However, Figure 5 shows a space-time plot (left panel) of a wave with ˜ T s = 9 · − as it encountersa sharp edge near r − r L = 100 km where τ changes from 1 to 0 .
2. The plot clearly shows that the long trailing waveis partially reflected as a long leading wave, which rapidly damps as it propagates back towards the resonance r = r L .The remaining part of the incoming trailing wave is transmitted into the rarefied region and quickly attenuates as itpropagates with strongly reduced wavelength. Note that the long trailing wave has a negative phase velocity − Ω L /k in our coordinate frame while its group velocity [Equation 42] is positive since k >
0. For the reflected leading wavewhich has k < r − r L (cid:46)
100 km, where thereflected wave has a substantial amplitude, the resulting density pattern behaves as a left-traveling (negative phasevelocity) wave which additionally undergoes a standing-wave motion. The standing-wave motion rapidly diminisheswith increasing distance from the edge due to the rapid damping of the reflected wave. The blue and red dashed curvesare curves of constant phase of the incoming and reflected wave, respectively, parametrized as [cf. Equation (40)] t ( x ) = t ± ω s x (cid:90) k ( s ) d s, (47)where k ( s ) is the wavenumber of a long density wave [Equation (43)]. The initial value t is chosen such that thecurves follow the path of a density maximum of the corresponding wave. The plot in the right panel shows the surfacedensity τ evaluated along these lines of equal phase, represented by the solid and dashed curves. From the solid curvewe can estimate the amplitude ratio of the incoming and the reflected waves by measuring the variation of τ near theedge as indicated by the arrows. That is, we have A max = A I + A R ,A min = A I − A R , where the subscripts I and R stand for the incoming and the reflected wave, respectively. Hence, A R A I = A max + A min A max − A min ≈ . , (48)which means that the major fraction of the incoming wave is reflected. We note that it is not clear how the reflectionis affected by our approximation of the azimuthal derivatives (31). Note also that the (viscous) time-scale on whichthe initially imposed density jumps change in a notable manner, is much longer then the applied integration times.9
20 40 60 80 100 120r-r L (km)05101520 t - [ O R B ]
40 50 60 70 80 90r-r L [km]0.80.91.01.11.21.3 τ IncomingReflected A max A min Figure 5 : Plots illustrating the reflection of a (long) trailing density wave at a sharp boundary at r − r L ∼
100 kmwhere τ is reduced by a factor of 1 /
5. In the space time plot (left panel) the blue (red) dashed curve traces a lineof equal phase of the incoming (reflected) wave so that it follows a density maximum. The plot in the right panelshows τ evaluated along these curves. As explained in the text, one can estimate the amplitude ratio of the incomingand the reflected waves from the indicated values A max and A min of τ [Equation (48)]. Since the considered wave isweakly nonlinear it follows H ( q ) (cid:38) σ (by 0 . k from (43), which is used in (47), to obtain a better fit to the locations of equal phase in the left panel.7.4. Density Waves and Viscous Overstability
We now turn to our hydrodynamical integrations of forced spiral density waves in a model ring which is sub-jected to viscous overstability such that β > β c ( λ ) [Equation (16)] for a non-zero range of wavelengths λ . Fig-ure 6 displays the linear stability curve for the P r β [Equation (7)] in the integrations discussed in the following. These values are β = 0 . , . , . , . , .
25 and 1 .
35. Viscous overstability is expected to develop for all but the smallest of thesevalues, resulting in wavetrains which are believed to produce parts of the observed periodic micro-structure in Sat-urn’s A and B rings (Thomson et al. (2007); Colwell et al. (2007); Latter and Ogilvie (2009)). For the values β = 1 . , .
16 and 1 .
20 linear viscous overstability is restricted to a relatively narrow band of wavelengths and theforced spiral density wave is stable. In contrast, for the two largest values β = 1 .
25 and 1 .
35 all wavelengths largerthan a critical one are unstable. For these two cases the forced spiral density wave itself is unstable and it is expectedfrom existing models that it retains a finite amplitude [i.e. a finite nonlinearity parameter q ] at large distance fromresonance (see BGT86 and LSS2016 for details).However, these models do not take into account the presence of the waves which are spontaneously generated byviscous overstability and which do not depend on the resonant forcing by an external gravitational potential. In thissection we study the interplay of both types of structure in a qualitative manner. All large-scale integrations presentedin this section were conducted using a grid with L r = 450 km, h = 25 m and applied Method A for the azimuthalderivatives (Section 5.1) as well as the
Straight Wire self-gravity model (Section 6).7.4.1.
Hydrodynamical Integrations without Forcing
For reference, Figures 7 and 8 describe an integration using β = 1 .
25, without forcing by the satellite ( ˜ T s = 0).The seed for this integration consists of a small amplitude superposition of linear left and right traveling overstablemodes on all wavelengths down to about 200 m. Note that without any seed and in the absence of satellite forcing,no perturbations develop. Figure 7 (left panel) shows a profile of τ (top) after about 20,000 orbital periods, along0 Lehmann et al.
100 1000 10000 λ [m]0.81.01.21.41.61.8 β c stableoverstable Figure 6 : Linear stability curve [Equation (16)] corresponding to the
P r β [Equation (7)] that are used for the large-scale integrationsof resonantly forced density waves discussed in Section 7.4. Viscous overstability is expected to develop for all valuesof β larger than the minimal value min[ β c ] λ = 1 .
03 which occurs for λ ∼
260 m. For β (cid:38) .
03 only a narrow band ofwavelengths is unstable. For β (cid:38) .
22 all wavelengths larger than a critical one are unstable, implying instability ofthe resonantly forced density wave.with its wavelet power (bottom). The structure on wavelengths λ ∼ −
400 m represents the nonlinear saturatedstate of viscous overstability. This state consists of left- and right traveling wave patches, separated by source andsink structures (Latter and Ogilvie (2009, 2010); LSS2017). This can also be seen in the stroboscopic space-timediagram (right panel), showing the evolution of τ over 600 orbits in the saturated state within a small portion of thecomputational domain near the nominal resonance location. The green dashed lines represent the expected nonlinearphase velocity v ph = ω I /k of overstable modes (Figure 8, left panel) of wavelength λ = 300 m, with ω I and k denotingthe nonlinear oscillation frequency and wavenumber of the wave. Although the modes seen in Figure 7 are in fact non-axisymmetric with azimuthal periodicity m = 7 (see Section 5), their phase velocity [cf. Equation (15)] is practicallythe same as for axisymmetric modes ( m = 0) since we are in a frame rotating with the orbital frequency at resonanceΩ L . Note that in Section 2 the symbol ω I was used to describe the linear oscillation frequency of overstable waves.The sharp decay of the density pattern near the domain boundaries is due to the inclusion of buffer-regions where β < min[ β c ] λ , so that the condition for linear viscous overstability is not fulfilled for any wavelength within theseregions. We included such buffer-regions in all large-scale integrations presented in the following. Furthermore, Figure8 (right panel) displays for the same integration as in Figure 7 the power spectral density of τ at two different times,as well as the evolution of the kinetic energy density (the insert). At the early time (200 ORB) the overstable wavesare still in the linear growth phase and the power spectrum corresponds directly to the linear growth rates ω R ( λ )(cf. Section 2 and the curve corresponding to β = 1 .
25 with q = 0 in Figure 11). During this stage the kineticenergy density increases rapidly. At later times nonlinear effects slow down the evolution and the power spectrum at t = 20 ,
000 ORB reflects the nonlinear saturation of the overstable waves.In LSS2017 we have shown that axisymmetric viscous overstability in a self-gravitating disk evolves towards a stateof minimal nonlinear oscillation frequency ω I (Figure 8, left panel), or equivalently, towards a state of vanishingnonlinear group velocity d ω I / d k . The dashed blue lines in Figure 7 (lower left panel) and Figure 8 indicate thewavelength corresponding to this frequency minimum of the nonlinear dispersion relation (by margins ±
20 m). Thewavelet power in Figure 7 reveals that in the region r > r L the saturation wavelength of viscous overstability is veryclose to the expected value. The overstable waves in this region are responsible for the sharp peak in the powerspectrum for t = 20 ,
000 ORB at λ (cid:46)
300 m (Figure 8). On the other hand, the region r < r L contains a left travelingwave with a wavelength that gradually departs from the expectation value towards the left, measuring λ ∼
400 m at theedge of the buffer-region. We observe the presence of weak long-wavelength undulations on top of the overstable wavesin the region r > r L . These mild, persistent undulations seem to adhere to the (long) density wave dispersion relationand result from the azimuthal derivative terms in the hydrodynamic equations (Section 5). They seem to prevent the1 L (km)0100200300400500600 t - [ O R B ] Figure 7 : Left : Radial surface density ( τ ) profile and its wavelet-power at t ∼ ,
000 ORB of an hydrodynamicalintegration using the
P r β = 1 .
25 and no satellite forcing ( ˜ T s = 0). The short-wavelength structuresare due to viscous overstability. The red dashed curve represents the linear density wave dispersion relation (43), whichsome persistent small amplitude undulations, resulting from the azimuthal derivative terms [Section 5], seem to follow.The blue dashed lines indicate the expected nonlinear saturation wavelength of viscous overstability by margins ±
20 m(see the text). The initial state of the integration is a small amplitude linear combination of left and right travelinglinear overstable waves on all wavelengths down to about 200 m.
Right : Stroboscopic space-time diagram showing theevolution near t = 20 ,
000 ORB of a small radial section at the nominal resonance location. Two source structures arelocated at x ∼ x ∼
14 km, respectively, sending out traveling waves both radially inward and outward. Inbetween the sources (at x ∼ ω I /k of nonlinear overstable waves, obtained from Figure 8 (left panel), for a wavelengthof λ = 300 m. Since the space-time diagram is stroboscopic, the apparent phase velocity of the waves in this plot isreduced (in absolute value) by Ω L /k , compared to the value obtained from the curve in Figure 8.saturation wavelength of overstable waves in the region r > r L to exceed the nonlinear frequency minimum. As such,the azimuthal derivatives seem to effectively remove the artificial influence of the periodic boundary conditions on thelong-term nonlinear saturation of the viscous overstability in that they sustain mild perturbations on the wave trains,at least in the region r > r L (see Latter and Ogilvie (2010) and LSS2017 for more details). This is further illustratedin Figure 23 (Appendix C) where we compare these results with those of an integrationwithout the azimuthal derivative terms. These plots confirm our finding in Section 2 that the linear behavior of small-scale axisymmetric and non-axisymmetric overstable modes is identical in our model. Furthermore, the nonlinearsaturation behavior is very similar as well. Due to the lack of persistent perturbations in the axisymmetric ( m = 0)integration all but one of the source and sink pairs will eventually merge and disappear so that the entire box willbe filled out by a single wavetrain which originates from the left buffer-region and whose wavelength increases withincreasing distance from its origin. At the time t = 35 ,
000 ORB its wavefront, which travels with a group velocity ofseveral meters per orbital period, has reached a radial distance of r − r L ∼
60 km. Latter and Ogilvie (2010) havepointed out that this long-term behavior is actually an artifact of the applied periodic boundary conditions.Note that due to the relatively low grid-resolution the wave profiles are not fully developed since the higher harmonicsare diminished. This, and also the effect of the azimuthal derivatives should, however, not affect our qualitativediscussion of the interaction between the density wave and viscous overstability in the following.2
Lehmann et al.
200 300 400 500 600 700 800 λ [m]0.940.960.981.00 ω I [ Ω ] linearnonlinear
200 400 600 800 λ [m]0.00.20.40.60.81.0 no r m a li z ed FFT - PS D t= 200 ORBt=20000 ORB ]0.000.050.100.150.200.250.30 e k i n [ σ c ] Figure 8 : Left : Linear and nonlinear frequencies of overstable waves adopting the
P r β = 1 . Right : Power spectral density of τ from the same integration as in Figure7 during the linear growth phase ( t = 200 ORB) and in the saturated state ( t ∼ ,
000 ORB) of viscous overstability.The insert displays the evolution of the kinetic energy density (22). The blue dashed lines in both panels indicate theexpected nonlinear saturation wavelength of viscous overstability with margins ±
20 m (see the text).7.4.2.
Co-Existence of Density Waves and Viscous Overstability
In Figure 21 (Appendix C) we compare integrations with a fixed forcing strength ˜ T s = 9 · − and varying valueof the viscosity parameter β . The first integration shows the same case as considere in Section 7.2 with β < min[ β c ] λ so that no viscous overstability develops. The integrations in rows 2-4 use min[ β c ] λ < β < β c ( λ → ∞ ) ≡ β ∞ c , whilethe integration shown in the bottom panel adopts β > β ∞ c . From top to bottom the results show an increasingsaturation amplitude of the viscous overstability in the evanescent region of the density wave ( r < r L ), as well asfor large distances r (cid:29) r L from resonance, where the density wave is already strongly damped. Furthermore, as areaction on the increased value of β the amplitude of the density wave shows a mild increase as well, particularly atlarger distances. These trends are expected from existing models for the nonlinear saturation of viscous overstability(Schmidt and Salo (2003); Latter and Ogilvie (2009)), as well as the BGT and the WNL models for nonlinear densitywaves. However, the behavior seen in the integration with β = 1 .
35 is not correctly described by the latter models.Since in this case β > β ∞ c , i.e. all wavelengths should be overstable, it is expected from these models that the densitywave does not damp but retains a finite (saturation) amplitude at large distances from resonance (BGT86; LSS2016).In contrast, our integration shows a damping of the wave very similar to the cases with β (cid:46) β ∞ c . In this case theviscous overstability possesses a sufficiently large amplitude to withstand the perturbation by the density wave atall distances r − r L , albeit with strongly diminished amplitude in the region of largest density wave amplitude. Incontrast, in the first three integrations viscous overstability is fully damped for a range of distances where the densitywave amplitude takes the largest values.Figure 22 (Appendix C) shows a series of integrations with increasing forcing strength ˜ T s = 10 − − .
16 and fixedvalue β = 1 . > β ∞ c . The first wave, excited by a small torque ˜ T s = 10 − is a linear wave. The development ofthe viscous overstability is very similar to the case without forcing (Figure 7). With increasing torque, the overstablewaves become increasingly distorted by the density wave, showing many similarities to those in Figure 21. Eventually,the density wave in the bottom panel is sufficiently strong to suppress viscous overstability in the far wave zone, andthe former wave attains a finite saturation amplitude at large distance from resonance until it hits the buffer-zone near3 x = 300 km. At times t (cid:38) ,
000 ORB there remain small distortions in the wave profile. It is possible that theseresult from the approximative treatment of the azimuthal derivatives (Section 5). The saturation amplitude τ ∼ . τ ∼ .
15. It is possible that this is aconsequence of our approximation for the azimuthal derivatives.Some details of the wave patterns encountered in our integrations are illustrated in Figures 9 and 10, which describethe integration with β = 1 .
20 and ˜ T s = 9 · − (same as in Figure 21, third row). Figure 9 shows a stroboscopicspace-time plot of a section of the radial τ -profile for times t = 3 , − ,
000 ORB. During this time the densitywave front traverses the considered region (indicated by the blue solid line) and clears overstable waves past a radialdistance r − r L (cid:38)
90 km (see also Figure 21, fourth row). The density wave corresponds to the nearly vertical patternwith radially decreasing wavelength λ ∼
10 km − ω I /k of these waves, assuming a wavelength λ = 250 m (Figure 8, left panel). Note that thefrequencies drawn in Figure 8 correspond to β = 1 .
25, but the dependence of the overstable frequency on β is weak sothat the corresponding curves for β = 1 .
20 are almost identical. The wavelength of overstable waves is modulated asthey traverse the peaks and troughs of the density wave. That is, the green dashed line matches quite well the phasevelocity of the overstable waves within the density wave peaks. In the troughs the phase velocity is notably increased,which follows from the decreased wavenumber of the overstable waves in these regions. Furthermore, a profile of τ at t = 1 ,
600 ORB (as marked by the arrow) is overplotted. If we assume that the density wave at a given time canbe described through Equation (35), we can estimate q ∼ .
25 in the region where overstability is damped. The reddashed lines in Figure 9 indicate minimum and maximum values of τ resulting from (35) for q = 0 .
25. In a similarmanner it follows that values of q ∼ .
23 and q ∼ .
20 lead to a damping of viscous overstability in the cases β = 1 . β = 1 .
10 (Figure 21), respectively. The associated q -values where overstability reappears at larger distances fromresonance seem to be slightly smaller in all cases. The mitigation of viscous overstability by a density wave will bediscussed in more detail in the following section.Figure 10 shows an orbit-resolved space-time plot of τ , illustrating how overstable wavetrains are distorted in directvicinity of the Lindblad resonance in the integration displayed in Figure 21 corresponding to β = 1 .
20. In the evanescentregion r < r L we recognize a right traveling overstable wave whose phase velocity undergoes periodic perturbationson the orbital time scale. These perturbations become stronger as the wave approaches the resonance r = r L . Inthe region r > r L the overstable waves seem to be unable to travel over any notable distance as their phase velocityrapidly changes its sign. In this region the amplitude of the overstable waves becomes strongly diminished (cf. Figure21, fourth row). 7.4.3. Viscous Overstability in a Perturbed Ring: Axisymmetric Approximation
The hydrodynamical integrations presented above reveal a variety of structures resulting from interactions betweena spiral density wave and the free short-scale waves associated with spontaneous viscous overstability. We find that asufficiently strong spiral density wave completely mitigates the growth of viscous overstability. In this section we willconsider this aspect in a more simplified, axisymmetric model which, on the one hand, allows us to conduct a simplehydrodynamic stability analysis, and, on the other hand, can be investigated with local N-body simulations. In whatfollows we assume that the perturbation is due to a nearby ILR. To that end consider the axisymmetric equations ∂ t τ = − u∂ x τ − τ ∂ x u, (49) ∂ t u = − u∂ x u + (Ω L r + 2Ω L v ) − ∂ x (cid:2) φ d + φ p (cid:3) − σ ∂ x ˆ P xx , (50) ∂ t v = − u∂ x v − L u − σ ∂ x ˆ P xy , (51)in the shearing sheet approximation (Goldreich and Lynden-Bell (1965)), using a rectangular frame ( x, y = 0) rotatingwith Ω L where x is given by (3) and where, in contrast to Equations (1), Ω = Ω L is now a constant. Note that thecomponents of the pressure tensor ˆ P xx and ˆ P xy are identical to ˆ P rr and ˆ P rθ given by (5), respectively, since we hadalready neglected curvature terms in the latter expressions. Moreover, v denotes here the total azimuthal velocity4 Lehmann et al. ✻(cid:0) ✼(cid:0) ✽(cid:0) ✾(cid:0) ✶(cid:0)(cid:0) ✶✶(cid:0)r✁r▲ ✥✂✄☎(cid:0)✺(cid:0)(cid:0)✶(cid:0)(cid:0)(cid:0)✶✺(cid:0)(cid:0)t✆✝✞✞✞✟✠✡☛☞
Figure 9 : Stroboscopic space-time diagram of a 60 km-section of τ resulting from the same hydrodynamical integrationusing β = 1 .
20 and scaled torque ˜ T s = 9 · − as displayed in Figure 21 for times t = 3 , − ,
000 ORB. Theover-plotted τ -profile describes the state at time t = 1 ,
600 ORB and the red dashed lines indicate the maximum andminimum values of τ predicted by Equation (35) for q = 0 .
25. The blue solid line describes the expected locationof the density wave front based on the group velocity (42). Overstable waves at radial distances r − r L (cid:38)
90 km arefully damped once the density wave front has passed this region. At considerably larger radial distances where thedensity wave has damped substantially, overstability reappears (see Figure 21, fourth row). The green dashed lineindicates the nonlinear phase velocity ω I /k of overstable waves with λ = 250 m (Figure 8). Overstable waves existingat distances r − r L ∼ −
110 km before the density wave front arrives possess small amplitude perturbations (theslowly left-traveling narrow features). These are expected to propagate with the nonlinear group velocity d ω I / d k ofthe overstable waves (Latter and Ogilvie (2010); LSS2017), which is small since the wavelength λ of the waves is closeto the (nonlinear) frequency minimum (Figure 8).and φ p is the gravitational potential due to the planet. We now introduce the perturbed oscillatory ground state(Mosqueira (1996)) τ = 11 − q cos ( mφ + m ∆) ,u = Ω L x q sin ( mφ + m ∆)1 − q cos ( mφ + m ∆) ,v = −
32 Ω L x − (4 / q cos ( mφ + m ∆)1 − q cos ( mφ + m ∆) . (52)As shown in Appendix G Equations (52) are valid in the vicinity of a Lindblad resonance where fluid streamlines canbe described by m -lobed orbits r ( a ) = a [1 − e ( a ) cos ( mφ + m ∆)] , (53)in a cylindrical frame ( r, φ, z = 0) rotating with the satellite’s mean motion frequency Ω s = ˆ ω s /m [Equation (23)] inthe present context. The quantities a and e denote a streamline’s semi-major axis and eccentricity, respectively and5 ✵ ✶✵ ✷✵ ✸✵ ✹✵r(cid:0)r▲ ✥✁✂✄✵✺✶✵✶✺✷✵✷✺t☎✆✝✝✝✝✞✟✠✡☛ Figure 10 : Orbit-resolved space-time diagram of τ for a region near the density wave resonance resulting from theintegration shown in Figure 21 with β = 1 .
20 for times t (cid:38) ,
000 ORB. The nearly horizontal pattern which becomesincreasingly pronounced with increasing r > r L represents the density wave. The smaller-scale structures are overstablewaves which are perturbed by the satellite resonance, causing the ”wiggles” in their appearance, in contrast to thewaves displayed in Figure 7. The green dashed line is the expected phase speed ω I /k of (unperturbed) overstablewaves with λ = 300 m obtained from Figure 8. A profile of τ is drawn for the time indicated by the arrow, showinghow the amplitude of overstable waves is reduced in approaching the resonance from smaller radii. Visible as well arethe first wave-cycles of the density wave.∆ is a phase angle. Furthermore, q is the nonlinearity parameter (Borderies et al. (1983)) fulfilling q = (cid:114) d( ae )d a + mae d∆d a . (54)In the limit q → r, θ, z = 0) which rotates with the local Kepler frequency Ω L at the resonance, we have mφ = mθ − m (cid:18) ˆ ω s m − Ω L (cid:19) t = mθ + Ω L t. (55)Using (52) and (55) the Equations (49)-(51) are identically fulfilled if one assumes a consistently expanded planetarypotential φ p = − Ω L (cid:2) r L − r L x + x (cid:3) at y = 0 and if one neglects the radial dependencies of the phase angle ∆ and the nonlinearity parameter q . Notethat the terms arising from orbital advection [the azimuthal derivatives (Section 5)] are neglected in the axisymmetricequations (49)-(51). These terms would scale relative to the other terms as x/r L (cid:46) − in the present situation.Since we assume that we are close to the Lindblad resonance ( | x | →
0) we can ignore the radial variation of ∆ in thearguments of the sine and cosine functions appearing in (52). That is, in the evanescent region close to the resonance( x (cid:46)
0) one can approximate q ∼ a d e/ d a since the eccentricity increases steeply towards the resonance and the disk’s6 Lehmann et al. response to the perturbation is not wavelike, i.e. d e/ d a (cid:29) e d∆ / d a (see for instance Hahn et al. (2009)). For x (cid:38) q [Eq. (54)] arises only from the radial variation of thephase angle ∆. This is the tight-winding approximation for the disk’s response in form of a long spiral density wavepropagating outward with radial wavenumber m d∆ / d a (cid:29) /a . However, even in this region we can approximate ∆ asa constant in (52), as long as the wavenumber of the density wave is much smaller than that of the overlying periodicmicro structure that we wish to analyze. Thus, the neglect of the radial variation of ∆ restricts the applicability ofthe above model to a small region at the resonance, since for sufficiently large x > not much greater than that of the overstable waves (cf. Figures 21, 22). As for the necessary approximationof a constant q we rely on previous studies which imply that for typical length scales of overstable wavetrains (severalkilometers) q varies slowly (see for instance Figure 3 of Borderies et al. (1986), Figure 2 of Hahn et al. (2009), as wellas Longaretti and Borderies (1986) and Rappaport et al. (2009) on the Mimas 5:3 wave). Note that the tight-windingapproximation applied in Borderies et al. (1986) inevitably assumes q ( x = 0) = 0).Thus, in what follows we assume that the phase angle ∆ is a constant and we assume (without loss of generality)that the ring is in the uncompressed state at initial time t = 0. Then Equations (52) together with (55) yield τ = 11 − q sin (Ω L t ) ,u = − Ω L x q cos (Ω L t )1 − q sin (Ω L t ) ,v = −
32 Ω L x − (4 / q sin (Ω L t )1 − q sin (Ω L t ) . (56)To the ground state (56) we now add axisymmetric perturbationsΨ( x, t ) = τ (cid:48) ( x, t ) u (cid:48) ( x, t ) v (cid:48) ( x, t ) = ˆ τ ( x, t )ˆ u ( x, t )ˆ v ( x, t ) exp [ ik ( t ) x ] (57)with time-dependent wavenumber k ( t ) = k − q sin Ω L t . (58)The time dependence in (58) stems from the periodic variation of the radial width of a streamline resulting from theperturbation by the density wave. This ansatz is chosen since Equations (56) contain only the kinematic effect ofthe density wave on the considered ring region. That is, (56) describe how a single fluid streamline behaves in thepresence of a density wave whose wavelength is assumed much larger than the extent of the streamline. For a nonlinearstudy of viscous overstability (implying longer time scales) in a perturbed ring region, the dynamical evolution of astreamline due to neighboring streamlines should be considered as well, which requires a more sophisticated treatmentthan the one adopted here. The behavior of the wavenumber according to Equation (58) is also seen directly in N-bodysimulations (cf. Figure 14).We assume that the x -dependency of the quantities ˆ τ ( x, t ), ˆ u ( x, t ) and ˆ v ( x, t ) in (57) is only weak so that ∂ x ˆ τ ( x, t )ˆ u ( x, t )ˆ v ( x, t ) (cid:28) k ( t ) ˆ τ ( x, t )ˆ u ( x, t )ˆ v ( x, t ) (59)and can be ignored.We insert the resulting expressions for τ = τ + τ (cid:48) , u = u + u (cid:48) , v = v + v (cid:48) into (49)-(51) and linearize with respectto the perturbations τ (cid:48) , u (cid:48) and v (cid:48) . This procedure yields the linear system ∂ t Ψ( x, t ) = ˆ M ( x, t )Ψ( x, t ) (60)7where the radial location x is a parameter andˆ M ( x, t ) = M M M M M M M M M (61)with M = q (1 + ik ( t ) x ) J − cos t,M = − ik ( t ) J − ,M = 0 ,M = i (cid:2) g − k ( t ) (cid:0) J + ( β + 1) αν qJ − β cos t (cid:1)(cid:3) ,M = J − q cos t [1 + ik ( t ) x ] − αν k ( t ) J − β ,M = 2 ,M = 12 i ( β + 1) k ( t ) ν J − β (4 q sin t − ,M = − J,M = − k ( t ) (cid:2) k ( t ) ν J − β − iqJ − x cos t (cid:3) . (62)To arrive at (62) we also used (13) and (17) where k x has been replaced by (58). Here we apply scalings so that timeand length are scaled with 1 / Ω L and c / Ω L , respectively. We also define the dimensionless quantities J = 1 − q sin t,α = 43 + γ, (63)for notational brevity. To illustrate the procedure for obtaining (62) let us consider the linearization of the continuityequation (49). The latter can be written as ∂ t (cid:16) τ + τ (cid:48) (cid:17) = − (cid:16) u + u (cid:48) (cid:17) ∂ x (cid:16) τ + τ (cid:48) (cid:17) − (cid:16) τ + τ (cid:48) (cid:17) ∂ x (cid:16) u + u (cid:48) (cid:17) . Since the ground state quantities are an exact solution in the current approximation, we end up with ∂ t τ (cid:48) = − u ∂ x τ (cid:48) − τ (cid:48) ∂ x u − τ ∂ x u (cid:48) , where we used ∂ x τ = 0 [Equation (56)]. Using (57) and (59) yields ∂ t τ (cid:48) = − [ ik ( t ) u + ∂ x u ] τ (cid:48) − ik ( t ) τ u (cid:48) . By applying (56), (58) and (63), as well as the aforementioned scalings, we obtain M , M and M as given in (62).All other matrix components are derived in the same fashion.The aim is now to investigate whether a seeded overstable wavetrain (57) will decay or grow in amplitude byintegrating (60) over a given time range. In the case q = 0 one can assumeΨ( x, t ) q =0 = τ (cid:48) ( x, t ) u (cid:48) ( x, t ) v (cid:48) ( x, t ) = ˆ τ ˆ u ˆ v exp [ ωt + ik x ] , Lehmann et al. i.e. the solution is a traveling wave with constant growth (or decay) rate, determined by the imaginary part of ω (Schmit and Tscharnuter (1995); Schmidt et al. (2001); Latter and Ogilvie (2009)). For q >
0, the behavior is morecomplicated.We integrate the complex-valued system of equations (60) numerically with a 4th-order Runge-Kutta method on agrid of radial size L x = 2 km. As initial state we use an eigenvector of (61) in the limit q → β = β c which reads (Schmidt and Salo (2003))Ψ( x, t ) = k ( ik ν + s )2 s ( ik ν + s ) − is + k ( ν + αν s ) + αν k exp [ ik x ] , (64)where s = (cid:113) − gk + αk ν + k is the unperturbed frequency of the overstable mode at marginal stability [c.f. Equation (15)]. The initial statecorresponds to a right traveling wave.In order to obtain the growth rate of a seeded mode λ n = L x /n where n denotes the radial mode number, we writeΨ( x, t ) = A n ( t ) exp [ ik n ( t ) x ]where k n ( t ) = (2 π/λ n ) J − ( t ). The complex amplitude is then obtained by numerical solution of A n ( t ) = 1 L x L x / (cid:90) − L x / d x Ψ( x, t ) exp [ − ik n ( t ) x ] (65)for each time step. Since (64) is not an exact eigenvector of (61) for q >
0, the first orbital periods of integrations with q > A n for different radial modes n resulting with the P r β and q are drawn in Figure 11. These plots show a monotonic decrease of the growth rates with increasing nonlinearityparameter q on all wavelengths. At the same time the maxima of the curves shift towards larger wavelengths. Fromthese plots we can estimate for given β the critical values q c that yield negative growth rates on all wavelengths. In thepresence of a perturbation with q ≥ q c no axisymmetric viscous overstability is expected to develop. The so obtainedvalues of q c seem to agree quite well with those estimated from the large-scale integrations in Section 7.4.2.As an illustration, in Figure 24 (Appendix C) we show space-time diagrams of the radial velocity perturbation u (cid:48) of the mode n = 10 for the case β = 1 .
35 with different values of q . While for q = 0 and q = 0 . q = 0 . q the overstable pattern becomes less of a uniform traveling wave.The waves seen in these space-time plots show many similarities to the overstable waves encountered in our large-scaleintegrations (Figure 10).In terms of the here applied hydrodynamical model the mitigation of overstable oscillations by the satellite perturba-tion can be explained by an increasing de-synchronization of specific terms appearing in the dynamical Equation (51).That is, the viscous overstability mechanism describes a transfer of energy from the background azimuthal shear intothe epicyclic fluid motion through a coupling of the viscous stress to the Keplerian shear (Latter and Ogilvie (2006,2009)), resulting in an oscillating angular momentum flux which instigates the epicyclic oscillation if the following twoconditions are met (Latter and Ogilvie (2006, 2008)). On the one hand, the viscous stress needs to possess a sufficientlysteep dependence on the surface mass density. This condition is expressed in terms of a critical (wavelength-dependent)viscosity parameter β c ( λ ) that must be exceeded for a given wavelength λ . Figure 6 shows β c ( λ ) for an unperturbedring ( q = 0), with minimal value Min [ β c ] λ = 1 .
03. Figure 11 reveals that β c ( λ ) increases with increasing q for all λ .For instance, for q = 0 . . > Min [ β c ] λ > .
1. Furthermore, for q = 0 . β c ] λ ∼ . β =1.10
100 200 300 400 500 λ [m]-0.04-0.020.000.02 ω R [ Ω ] q=0.0q=0.1q=0.2 q=0.3q=0.4q=0.5 β =1.16
100 200 300 400 500 λ [m]-0.04-0.020.000.02 ω R [ Ω ] q=0.0q=0.1q=0.2 q=0.3q=0.4q=0.5 β =1.25
100 200 300 400 500 λ [m]-0.04-0.020.000.020.04 ω R [ Ω ] q=0.0q=0.1q=0.2 q=0.3q=0.4q=0.5 β =1.35
100 200 300 400 500 λ [m]-0.04-0.020.000.020.04 ω R [ Ω ] q=0.0q=0.1q=0.2 q=0.3q=0.4q=0.5 Figure 11 : Linear hydrodynamic growth rates of overstable modes in a perturbed ring for different values of thenonlinearity parameter q describing the satellite perturbation. The wavelengths λ correspond to the uncompressedstate of the model ring adopted at times t = lπ/ Ω L , with non-negative integer l [cf. Equation (58)]. The usedparameters are the P r β and all growth rates are scaled with Ω = Ω L . In allpanels the growth rates for q = 0, obtained by numerical solution of (14), are plotted additionally as black solid curves.perturbation v (cid:48) that describes the coupling of the viscous stress to the Keplerian shear and which can be written M τ (cid:48) [cf. (60)-(62)]. The used parameters are the same as in Figure 24 and the snapshots cover one orbital periodin equal time-intervals. Over-plotted for the same instances of time is the epicyclic term M u (cid:48) , appearing in thesame equation. For q (cid:46) . q the phase difference of these terms drifts constantly. Therefore the energy transfer into the epicyclic oscillation is tooinefficient, resulting in a damping of seeded wavetrains.To quantify the phase relation between the angular momentum flux and the epicyclic oscillation associated with anoverstable wave we define the cross correlation of the two aforementioned relevant terms CC ( t s ) = (cid:90) [ l π ] (cid:104) M ( t ) τ (cid:48) ( t ) M ( t + t s ) u (cid:48) ( t + t s ) (cid:105) d t, (66)where the shift parameter t s takes values between 0 and 2 π and the integration is performed over l orbital periods.If the two quantities M τ (cid:48) and M u (cid:48) are periodic with a constant phase shift, their cross correlation will posses a0 Lehmann et al. sharp maximum for a specific value of the time shift t s . On the other hand, if their relative phase shift varies in amore or less uniform manner during one orbital period, the cross correlation will be small for all values of t s . Thus,as a measure for the synchronization of the two oscillatory quantities we take the maximum absolute value of thecross correlation max |CC| t s . Figure 12 displays this value for the parameters as used in Figure 11. As anticipated, thecurves show a steady decrease with increasing q and exhibit a fairly sharp drop for q = 0 . − .
4. This explains whyfor q (cid:38) . β used here. To explain the different criticalvalues q c for which the growth rates become negative for the different β -values (e.g. q c ∼ . β = 1 . β . In the unperturbed case ( q = 0) the growthrates depend linearly on the factor ( β − β c ). max. cross correlation m a x | CC ( t s ) | β = . - . Figure 12 : Curves showing the maximum absolute value of the cross correlation of the quantities M τ (cid:48) and M u (cid:48) [Equation (66)] associated with an overstable mode with λ = 200 m for different values of q , and the values of β usedfor the plots in Figure 11. The arrow indicates the direction of increasing β . The monotonic decrease of these curveswith increasing q is the reason why the growth rates of overstable modes (Figure 11) become smaller with increasing q . All curves have been normalized to yield unity for q = 0.Note that in a dilute ring, better described in terms of a kinetic model than a hydrodynamic one, the aforementionedde-synchronization can occur already in absence of an external perturbation. In a kinetic model of a dilute ringof sufficiently low dynamical optical depth the viscous stress tensor components are subjected to long (collisional)relaxation time-scales. Therefore, these cannot follow the (fast) epicyclic oscillation of the ring flow on the orbital timescale, which also prevents viscous overstability (Latter and Ogilvie (2006)).We complement these findings with a series of local N-Body simulations of a perturbed ring that include aspectsof the vertical self-gravity force in terms of an enhancement of the frequency of vertical oscillations (Wisdom andTremaine (1988)) and also the effect of collective radial self-gravity. A detailed description of the simulation methodcan be found in Salo et al. (2018) and references therein. In particular, the force method for particle impacts asintroduced in Salo (1995) is used and the radial self-gravity is calculated as in Salo and Schmidt (2010). The latter isparametrized through a pre-specified ground state surface mass density (see also LSS2017). Here we apply modifiedinitial conditions and boundary conditions that account for a perturbed mean flow in the ring following the methodof Mosqueira (1996). We perform simulations with meter-sized particles in a periodic box of uncompressed radial size L x = L / (cid:112) − q and azimuthal size L y = 10 m. The number of particles is slightly less than 10,000 in all simulations.The quantity L x is chosen such that the time-averaged ground state optical depth of the system is the same for differentvalues of q , as shown below. The radial size of the simulation region changes periodically as L x ( t, q ) = L − q cos Ω L t (cid:112) − q . (67)1 n=15 A m p . [ a . u .] q=0.0 q=0.1 q=0.2 q=0.3q=0.4 q=0.5 q=0.6 n=15 ω R [ Ω ] n=20 A m p . [ a . u .] n=20 ω R [ Ω ] n=25 A m p . [ a . u .] n=25 ω R [ Ω ] Figure 13 : Determination of linear growth rates of three different seeded overstable modes (indicated by their radialmode number n ) in N-body simulations for different values of the nonlinearity parameter q , quantifying the amountof perturbation in the ground state that corresponds to a density wave. The left panels display the time evolutionof amplitudes A n [Equation (65)] with a sampling interval of one orbital period. The right panels show the resultinggrowth rates, obtained from linear fits as drawn in the left frames. The simulations used time-averaged values of theground state optical depth and surface mass density of τ = 1 . σ = 300 kg m − , respectively, as well as a verticalfrequency enhancement of Ω z / Ω = 2 to mimic vertical self-gravity.For a fixed azimuthal width and a fixed number of simulation particles the ground state dynamical optical depth isthen given by τ dyn ( t, q ) = τ (cid:112) − q − q cos Ω L t (68)where τ is the time-averaged ground state optical depth over one period 2 π/ Ω L , i.e. τ ≡ (cid:104) τ dyn ( t, q ) (cid:105) t and is indepen-dent of q . Hence, our choice of L x removes the purely geometrical increase of the time-averaged mean optical depthwith increasing q and isolates the effect of the perturbation on the evolution of viscous overstability. Note that thequantity τ is not to be confused with the scaled surface density τ used elsewhere in this paper. In what follows we use τ = 1 .
5, and L = 2 km. The ground state surface mass density σ ( t, q ) of the simulation region will vary in the sameway as the optical depth (68) and we adopt a time-averaged ground state surface mass density of σ = 300 kg m − .Furthermore, the simulations utilize the velocity dependent normal coefficient of restitution by Bridges et al. (1984),while particle spins are neglected.Figure 13 shows measurements of linear growth rates of three different seeded overstable modes in N-body simulationsusing different values of the nonlinearity parameter q = 0 − .
6. The initial state corresponds to a standing linearoverstable wave (see Equation (37) of Schmidt et al. (2001) expanded to the lowest order of the scaled wavenumber2
Lehmann et al. -1000 0 1000x [m]0.00.51.01.52.02.53.0 t i m e [ O R B ] Figure 14 : The evolution of the radial velocity perturbation u (cid:48) during the first three orbital periods of the N-bodysimulation with n = 25 and q = 0 . k ) in a phase where only the perturbation in the radial velocity u (cid:48) has a non-zero amplitude. When this radialperturbation velocity is seeded with small amplitude, then the simulation practically starts on an overstable eigenvectorof the linear hydrodynamic model. The radial modes n = 15 −
25 correspond to time-averaged wavelengths (cid:104) λ (cid:105) t =( L /n ) / (cid:112) − q = (133 m −
80 m) / (cid:112) − q [see Equation (67)]. Figure 14 illustrates the evolution the the radialvelocity perturbation for the case n = 25 with q = 0 . u (cid:48) ( x ).In accordance with the hydrodynamic growth rates (Figure 11) the measured growth rates in N-body simulations(the right panels in Figure 13) decrease with increasing magnitude of the perturbation, quantified through q . Note thatthe overstable modes considered here would be stable in the hydrodynamic model for all used values of β (Figure 11).Since our N-body simulations do not include particle-particle gravitational forces (but merely the radial componentof its mean-field approximation), the attained (equilibrium) velocity dispersion takes considerably smaller values thanwhat we adopt for our hydrodynamic integrations (Table 1). This is the main reason why viscous overstability in theN-body simulations considered here occurs on smaller wavelengths than in our hydrodynamic model (LSS2017).For the same parameters Figure 15 illustrates the nonlinear saturation of overstability in simulations that startedfrom white noise. The left column displays curves of the kinetic energy density of perturbations that have developedon top of the ground state (52). The curves are sampled at quadrature (i.e. at times t where Ω L t = π/ l π withinteger l ) where the radially averaged (mean) optical depth takes the value (cid:104) τ dyn (cid:105) x = τ (cid:112) − q so as to mask outthe orbital oscillation due to the perturbation itself. While the curves for q = 0 − . q ≥ .
3. The right column shows snapshots of the radial particle number density profile near the endof each simulation where also (cid:104) τ dyn (cid:105) x = τ (cid:112) − q . In agreement with the energy curves we find clearly developednonlinear overstable wavetrains for q = 0 − .
2, while for q ≥ . × -3 × -3 e k i n × -3 × -3 e k i n × -3 × -3 e k i n × -3 × -3 e k i n × -3 × -3 e k i n × -3 × -3 e k i n × -3 × -3 e k i n -1000 1000x [R p ]0.51.52.5 σ / < σ > x q=0.6 σ / < σ > x q=0.5 σ / < σ > x q=0.4 σ / < σ > x q=0.3 σ / < σ > x q=0.2 σ / < σ > x q=0.1 σ / < σ > x q=0.0 Figure 15 : Nonlinear evolution of viscous overstability in N-body simulations of a perturbed ring.
Left : Curves ofkinetic energy density (22) of perturbations in units σ R P Ω , excluding the ground state velocities (52). The curvesare sampled at quadrature such that the radially averaged optical depth takes the value (cid:104) τ dyn (cid:105) x = τ (cid:112) − q withtime-averaged ground state optical depth τ = 1 . Right : Snapshots of the radial surface density profile for each q taken at the same final time at quadrature where (cid:104) τ dyn (cid:105) x = τ (cid:112) − q and the radial width of the simulation regiontakes the value 2 km / (cid:112) − q . Note that the profiles are scaled with the radially averaged surface mass density (cid:104) σ (cid:105) x which oscillates about its time-average σ = 300 kg m − in the same way as (68) . Furthermore, a vertical frequencyenhancement of Ω z / Ω = 2 was used.From these results we estimate a critical value q c ∼ . DISCUSSIONWe developed a one-dimensional hydrodynamical scheme to study the excitation of a resonantly forced spiral densitywave in a dense planetary ring. Due to the restriction to one space-dimension, the advection caused by orbital shearneeds to be approximated. We constructed corresponding azimuthal derivative terms from the weakly nonlinear modelof LSS2016. Profiles of nonlinear density waves in a viscously stable ring computed with our scheme show goodagreement with those resulting from the models by BGT86 and LSS2016.We applied our scheme to investigate the damping behavior of spiral density waves in a planetary ring which issubject to viscous overstability. The results of our large-scale hydrodynamical integrations confirm the observationthat resonantly forced spiral density waves can co-exist with short-scale waves generated by the viscous overstability4
Lehmann et al. (Hedman et al. (2014)), an aspect not taken into account in existing models for the damping of density waves.Due to our approximation of the azimuthal derivative terms the free short-scale overstable modes appearing in ourhydrodynamical integrations are also non-axisymmetric with the same azimuthal periodicity m as the spiral densitywave. We have shown that the nonlinear evolution of these short-scale modes is very similar to that of the strictlyaxisymmetric short-scale overstable modes investigated in earlier studies.We find that the damping behavior of a spiral density wave can be very different from what is predicted by existingmodels, depending on its resonance strength. A sufficiently strong spiral density wave damps the short-scale viscousoverstability. Furthermore, if the density wave is sufficiently strong and it is itself overstable it behaves according tothe models by BGT86 and LSS2016 in that it retains a finite saturation amplitude at large distances from resonance.If, on the other hand, the density wave is overstable and sufficiently weak the short-scale modes dominate and dampthe density wave.It should be noted that these results are quantitatively (but most likely not qualitatively) affected by the approxima-tion of the azimuthal derivatives in our numerical scheme. That is, although we have shown that this approximationworks well if we consider a nonlinear density wave alone, or the nonlinear evolution of viscous overstability in absenceof a density wave, it cannot be ruled out that in cases where both wave types co-exist certain nonlinear terms in thehydrodynamic Equations (1) would produce spurious quasi-resonant higher-order coupling terms between both wavetypes. Such spurious terms would be quasi-resonant due to the approximation that all terms in Equations (1) areassumed to have m -fold periodicity and the fact that the wavelength of the density wave is much greater than thatof the short-scale overstable modes, at least in close vicinity of the resonance radius. It is very unlikely though thatsuch terms would dominate the many physical coupling terms. Therefore we believe that our findings are qualitativelycorrect despite the approximation of the azimuthal derivatives.We verified the damping of viscous overstability by the density wave by performing N-body simulations as wellas a linear hydrodynamic stability analysis of a simplified axisymmetric model for a ring perturbed by a nearbyILR. Our N-Body simulations, using modified initial and boundary conditions as introduced by Mosqueira (1996),confirm the formation of nonlinear overstable wavetrains if the perturbation by the ILR is not too strong, as well as acomplete suppression of overstability if the nonlinearity parameter q associated with the perturbation exceeds a certainvalue. Critical values of q which result in a damping of viscous overstability obtained from our large-scale integrationscompare well with those that follow from the linear stability analysis of the axisymmetric model. Based on our resultswe conclude that the mitigation of viscous overstability by a density wave is due to a destruction of the phase relationof the oscillating angular momentum flux and the epicyclic oscillation associated with overstable waves. Note thata quantitative match of q -values in the aforementioned comparison should not be expected though. That is, on theone hand the overstable waves found in our large-scale hydrodynamical integrations suffer from the relatively lowspatial resolution of the computational grid, which is expected to reduce the estimates of q c from these integrations.Furthermore, the applied approximation of the azimuthal derivatives could affect these values as well. On the otherhand, the neglect of the variation of the phase angle ∆ due to the density wave in the linear stability analysis is mostlikely not justified in the far wave region of a density wave. It is not clear how this affects the computed growth ratesof overstable modes and associated values of q c .Although we understand the mitigation of viscous overstability by a density wave, an explanation for the dampingof overstable density waves, such as those presented in the the last panel of Figure 21 and the first three panels of 22,remains to be sought for. One difficulty is that in this case the nonlinear interaction of the two different modes needsto be considered. It is noteworthy that the observed density waves in Saturn’s A ring associated with the 7:6 ILR andthe 10:9 ILR with the moons Atlas and Pan, respectively, seem to correspond to this case (see Hedman et al. (2014),Figure 5).Furthermore, it should be noted that due to the neglect of particle-particle self-gravity in our modeling, self-gravitational wakes (Salo (1992)) do not form. In principle their effect on the density wave profile may be describedin terms of a gravitational viscosity (Daisaka et al. (2001)). In parts of Saturn’s dense rings (particularly the A ring)it is expected that this gravitational viscosity is the dominant mode of viscosity. However, the wakes will interactwith viscous overstability in a more or less complex manner (Salo et al. (2001); Ballouz et al. (2017)) and as suchthey will indirectly affect the damping of a density wave. Moreover, in the regions of strong density waves the sizeof self-gravitational wakes is expected to be much increased. That is, the ”straw“-like structures observed in opticalCassini images of strong density waves (e.g. Tiscareno (2017)), are believed to represent self-gravity wakes of kilometerlength scales. These length scales are much greater than the typical length scale of self-gravity wakes that can form5in an unperturbed planetary ring (Salo (1992)). An increased size of self-gravity wakes in the troughs of nonlineardensity waves is also expected from N-body simulations (Salo and Schmidt (2014)) and theoretical studies (Stewart(2017)). Stewart (2017) has shown that the characteristic length scale of self-gravitational perturbations (the Toomre-wavelength) is greatly enhanced in the troughs of strongly nonlinear density waves. This result suggests that thegravitational viscosity, which scales with the square of the Toomre-wavelength (Daisaka et al. (2001)), can be greatlyenhanced in the wave region, consequently leading a stronger damping of the density wave.Our numerical scheme allows for in principle straightforward extensions, such as the inclusion of the energy equationby using the numerical method of LSS2017. Ultimately, a two-dimensional scheme should be developed to overcomethe necessity to approximate the orbital advection terms.ACKNOWLEDGMENTSWe thank the reviewer Glen Stewart for helpful and constructive comments. We are grateful to Pierre-Yves Longarettifor discussions that greatly improved the manuscript. We acknowledge support from the Academy of Finland. MLacknowledges funding from the University of Oulu Graduate School and the University of Oulu Scholarship Foundation.REFERENCES S. Araki and S. Tremaine. The dynamics of dense particledisks.
Icarus , 65:83–109, 1986.R. L. Ballouz, D. C. Richardson, and R. Morishima.Numerical Simulations of Saturn’s B Ring: GranularFriction as a Mediator between Self-gravity Wakes andViscous Overstability. AJ , 153:146–155, 2017.J. Binney and S. Tremaine. Galactic Dynamics . PrincetonUniversity Press, 1987.N. Borderies, P. Goldreich, and S. Tremaine. Perturbedparticle disks.
Icarus , 55:124–132, 1983.N. Borderies, P. Goldreich, and S. Tremaine. A granularflow model for dense planetary rings.
Icarus , 63:406–420,1985.N. Borderies, P. Goldreich, and S. Tremaine. Nonlineardensity waves in planetary rings.
Icarus , 68:522–533,1986.R. Borges, M. Carmona, B. Costa, and W. S. Don. Animproved weighted essentially non-oscillatory scheme forhyperbolic conservation laws.
Journal of ComputationalPhysics , 227:3191–3211, 2008.F.G. Bridges, A.P. Hatzes, and D.N.C. Lin. Structure,stability and evolution of Saturn’s rings.
Nature , 309:333–338, 1984.J. E. Colwell, L. W. Esposito, M. Sremˇcevi´c, G. R. Stewart,and W. E. McClintock. Self-gravity wakes and radialstructure of Saturn’s B ring.
Icarus , 190:127–144, 2007.M.C. Cross and P.C. Hohenberg. Pattern formationoutside of equilibrium.
Reviews of Modern Physics , 65(3):851, 1993.J. N. Cuzzi, J.J. Lissauer, L.W. Esposito, J.B. Holberg,E.A. , G.L. Tyler, and A. Boischot. Saturn’s rings:Properties and processes. In R. Greenberg andA. Brahic, editors,
Planetary Rings , pages 73–199. TheUniversity of Arizona Press, 1984. H. Daisaka, H. Tanaka, and S. Ida. Viscosity in a denseplanetary ring with self-gravitating particles.
Icarus ,154:296-312, 2001.P. Goldreich and D. Lynden-Bell. II. Spiral arms as shearedgravitational instabilities.
MNRAS , 130:125, 1965.P. Goldreich and S. Tremaine. The formation of the Cassinidivision in Saturn’s rings.
Icarus , 34:240–253, 1978.P. Goldreich and S. Tremaine. The excitation andevolution of density waves.
ApJ , 222:850–858, 1978.P. Goldreich and S. Tremaine. The excitation of densitywaves at the Lindblad and corotation resonances by anexternal potential.
ApJ , 233:857–871, 1979.J. M. Hahn, J. N. Spitale, and C. C. Porco. Dynamics ofthe Sharp Edges of Broad Planetary Rings.
ApJ , 699:686–710, 2009.M. M. Hedman and P. D. Nicholson. The B-ring’s surfacemass density from hidden density waves: Less than meetsthe eye?
Icarus , 279:109–124, 2016.M. M. Hedman, P. D. Nicholson, and H. Salo. ExploringOverstabilities in Saturn’s A Ring Using Two StellarOccultations. AJ , 148:15, 2014.G.-S. Jiang and C.-W. Shu. Efficient Implementation ofWeighted ENO Schemes. Journal of ComputationalPhysics , 126:202–228, 1996.H. N. Latter and G. I. Ogilvie. The linear stability ofdilute particulate rings.
Icarus , 184:498–516, 2006.H. N. Latter and G. I. Ogilvie. Dense planetary rings andthe viscous overstability.
Icarus , 195:725–751, 2008.H. N. Latter and G. I. Ogilvie. The viscous overstability,nonlinear wavetrains, and finescale structure in denseplanetary rings.
Icarus , 202:565–583, 2009.H. N. Latter and G. I. Ogilvie. Hydrodynamicalsimulations of viscous overstability in Saturn’s rings.
Icarus , 210:318–329, 2010. Lehmann et al.
M. Lehmann, J. Schmidt, and H. Salo. A Weakly NonlinearModel for the Damping of Resonantly Forced DensityWaves in Dense Planetary Rings.
ApJ , 829:75, 2016.M. Lehmann, J. Schmidt, and H. Salo. ViscousOverstability in Saturn’s Rings: Influence of CollectiveSelf-gravity.
ApJ , 851:125, 2017.M. S. Liou and C. J. Steffen. A New Flux Splitting Scheme.
Journal of Computational Physics , 107:23–39, 1993.P. Y. Longaretti and N. Borderies. Nonlinear study of themimas 5:3 density wave.
Icarus , 67:211–223, 1986.Longaretti, P.-Y. (1989). Uranian ring dynamics - ananalysis of multimode motions.
Icarus , 82:281–287.P. Y. Longaretti and N. Borderies. Streamline formalismand ring orbit determination.
Icarus , 94:165–170, 1991.P. Y. Longaretti. Theory of Narrow Rings and SharpEdges, in Planetary Ring Systems, pages 225–276, ed.M. S. Tiscareno, & C. D. Murray. Cambridge UniversityPress, 2018.D. Lynden-Bell and J. E. Pringle. The evolution of viscousdiscs and the origin of the nebular variables.
MNRAS ,168:603–637, 1974.Meyer-Vernet, N. and Sicardy, B. (1987). On the physics ofresonant disk-satellite interaction. ”Icarus” , 69:157–175.Ignacio Mosqueira. Local simulations of perturbed denseplanetary rings.
Icarus , 122:128–152, 1996.N. J. Rappaport, P.-Y. Longaretti, R. G. French, E. A.Marouf, and C. A. McGhee. A procedure to analyzenonlinear density waves in Saturn’s rings using severaloccultation profiles.
Icarus , 199:154–173, 2009.H. Rein and H. N. Latter. Large-scale N-body simulationsof the viscous overstability in Saturn’s rings.
MNRAS ,431:145–158, 2013.H. Salo. Numerical simulations of dense collisionalsystems.
Icarus , 90:254–270, 1991.H. Salo. Gravitational wakes in Saturn’s rings.
Nature ,359:619–621, 1992.H. Salo. Simulations of dense planetary rings. III.Self-gravitating identical particles.
Icarus , 117:287–312,1995.H. Salo, J. Schmidt, and F. Spahn. Viscous overstability inSaturn’s B ring: I. Direct simulations and mesurement oftransport coefficients.
Icarus , 153:295–315, 2001.H. Salo and J. Schmidt. N-body simulations of viscousinstability of planetary rings.
Icarus , 206:390–409, 2010.H. Salo and J. Schmidt. Excess noise in synthetic stellaroccultation data from N-body simulations of Saturn’srings.
European Planetary Science Congress ,9:EPSC2014–744, 2014. H. Salo, K. Ohtsuki, and M. C. Lewis.
ComputerSimulations of Planetary Rings , in Planetary RingSystems, pages 434–494, ed. M. S. Tiscareno, &C. D. Murray. Cambridge University Press, 2018.J. Schmidt and H. Salo. A weakly nonlinear model forviscous overstability in Saturn’s dense rings.
PhRvL , 90(6):061102, 2003.J. Schmidt, H. Salo, F. Spahn, and Olaf Petzschmann.Viscous overstability in Saturn’s B ring: II.Hydrodynamic theory and comparison to simulations.
Icarus , 153:316–331, 2001.J. Schmidt, K. Ohtsuki, N. Rappaport, H. Salo, andF. Spahn.
Dynamics of Saturn’s Dense Rings , pages413–458. 2009.J. Schmidt, J. E. Colwell, M. Lehmann, E. A. Marouf,H. Salo, F. Spahn, and M. S. Tiscareno. On the LinearDamping Relation for Density Waves in Saturn’s Rings.
ApJ , 824:35, 2016.U. Schmit and W.M. Tscharnuter. A fluid dynamicaltreatment of the common action of self-gravitation,collisions, and rotation in Saturn’s B-ring.
Icarus , 115:304–319, 1995.U. Schmit and W.M. Tscharnuter. On the formation of thefine–scale structure in Saturn’s B ring.
Icarus , 138:173–187, 1999.C.-W. Shu and S. Osher. Efficient Implementation ofEssentially Non-oscillatory Shock-Capturing Schemes.
Journal of Computational Physics , 77:439–471, 1988.F. H. Shu. Waves in planetary rings. In R. Greenberg andA. Brahic, editors,
Planetary Rings , pages 513–561,Tucson Arizona, 1984. Univ. of Arizona Press.F. H. Shu, L. Dones, J. J. Lissauer, C. Yuan, and J. N.Cuzzi. Nonlinear spiral density waves - viscous damping.
ApJ , 299:542–573, 1985.F.H. Shu, C. Yuan, and J.J. Lissauer. Nonlinear spiraldensity waves: an inviscid theory.
ApJ , 291:356–376,1985.I.G. Shukhman. Collisional dynamics of particles inSaturn’s rings.
SvA , 28:574–585, 1984.Spahn, F., Schmidt, J., Petzschmann, O., and Salo, H.(2000). Stability analysis of a Keplerian disk of granulargrains: influence of thermal diffusion.
Icarus ,145:657–660.G. R. Stewart, D. N. C. Lin, and P. Bodenheimer.Collision-induced transport processes in planetary rings.In Planetary Rings ed. R. Greenberg and A. Brahic,pages 447–512, Univ. of Arizona Press, Tucson Arizona,1984. G. R. Stewart. On the Extraordinary Propagation of theJanus 2:1 Density Wave: Synergy between Density Wavesand Viscous Overstability. In
AAS/Division forPlanetary Sciences Meeting , volume 48 of
AAS/Divisionfor Planetary Sciences Meeting , page 203.02, 2016.G. R. Stewart. Straw Formation and Enhanced Dampingof Strong Density Waves in Saturn’s Rings. In
AAS/Division of Dynamical Astronomy Meeting ,volume 48 of
AAS/Division of Dynamical AstronomyMeeting , page 400.02, 2017.F. S. Thomson, E. A. Marouf, G. L. Tyler, R. G. French,and N. J. Rappoport. Periodic microstructure inSaturn’s rings A and B.
GRL , 34:L24203.1–L24203.6,2007.M. S. Tiscareno, P. D. Nicholson, J. A. Burns, M. M.Hedman, and C. C. Porco. Unravelling TemporalVariability in Saturn’s Spiral Density Waves: Results andPredictions.
The Astrophysical Journal , 651:L65–L68,2006. M. S. Tiscareno, J. A. Burns, P. D. Nicholson, M. M.Hedman, and C. C. Porco. Cassini imaging of Saturn’srings II: A wavelet technique for analysis of density wavesand other radial structure in the rings.
Icarus ,189:14–34, 2007.M. S. Tiscareno and Cassini Imaging Team (2017).High-resolution imaging of Saturn’s main rings duringthe Cassini Ring-Grazing Orbits and Grand Finale. In
AAS/Division for Planetary Sciences Meeting Abstracts , volume 49 of
AAS/Division for Planetary SciencesMeeting Abstracts , page 108.02.A. Toomre. Group Velocity of Spiral Waves in GalacticDisks.
ApJ , 158:899–914, 1969.C. Torrence and G. P. Compo. A Practical Guide toWavelet Analysis.
Bulletin of the AmericanMeteorological Society , 79:61–78, 1998.J. Wisdom and S. Tremaine. Local simulations ofplanetary rings. AJ , 95:925–940, 1988. Lehmann et al.
APPENDIX A. FIGURES OF SECTION 7.2 -50 0 50 100 150 200 0.80.91.01.11.21.3 τ Hydro IntWNL Model T s = 9.00e-02 ILR 7: 6 -50 0 50 100 150 200 -20-1001020 u / c -50 0 50 100 150 200r-r L [km]-10-50510 v / c -50 0 50 100 150 200 0.80.91.01.11.21.3 τ Hydro IntWNL Model T s = 9.00e-02 ILR 7: 6 -50 0 50 100 150 200 -20-1001020 u / c -50 0 50 100 150 200r-r L [km]-10-50510 v / c -50 0 50 100 150 200 0.51.01.52.02.5 τ Hydro IntWNL Model T s = 1.00e+00 ILR 7: 6 -50 0 50 100 150 200 -80-4004080 u / c -50 0 50 100 150 200r-r L [km]-40-2002040 v / c -50 0 50 100 150 200 0.51.01.52.02.5 τ Hydro IntWNL Model T s = 1.00e+00 ILR 7: 6 -50 0 50 100 150 200 -80-4004080 u / c -50 0 50 100 150 200r-r L [km]-40-2002040 v / c Figure 16 : Comparison of state variables resulting from hydrodynamical integrations and the WNL model using the
P r T s = 9 · − (top panels) and ˜ T s = 1 (bottom panels). The integration shown in the left(right) panels applied method A ( method B ) for the azimuthal derivatives (Section 5).9 L [km]0.51.01.52.02.53.03.5 τ Hydro Int L [km]0.51.01.52.02.53.03.5 τ BGT L [km]0.51.01.52.02.53.03.5 τ WNL
Figure 17 : Comparison of profiles of τ along with their Morlet wavelet powers resulting from a hydrodynamicalintegration and the WNL and BGT models using the P r T s = 4. The dashed red lines representthe linear dispersion relation (43). Straight Wire SG (T s = 9.00e-02) L [km]0.60.81.01.21.41.6 τ BGTHydro Int
WKB SG (T s = 9.00e-02) L [km]0.60.81.01.21.41.6 τ BGTHydro Int
Straight Wire SG (T s = 4.00e+00) L [km]01234 τ BGTHydro Int
WKB SG (T s = 4.00e+00) L [km]01234 τ BGTHydro Int
Figure 18 : Comparison of profiles of τ resulting from integrations with the Straight Wire self-gravity (left column) andthe WKB self-gravity (right column) with corresponding waves of the BGT model. Note that the WKB-approximationis intrinsic to the BGT model.0
Lehmann et al. B. FIGURES OF SECTION 7.3
Figure 19 : Space-time plots of a m = 7 density wave with scaled torque ˜ T s = 9 · − passing through a region( r − r L ∼ −
100 km) of increased ( τ = 3, left frame) and decreased ( τ = 0 .
5, right frame) equilibrium surface massdensity. As in Figures 1-3 the blue solid curve indicates the expected curve of the density wave front in a homogeneousring [Equation (42)]. Note that these plots only show the density perturbation on top of the background density. L [km]0.00.51.01.52.02.53.03.5 τ L [km]0.00.51.01.52.02.53.03.5 τ L [km]0.00.51.01.52.02.53.03.5 τ Figure 20 : Comparison of profiles of τ along with their Morlet wavelet powers resulting from hydrodynamicalintegrations using the P r T s = 9 · − . From top to bottom the equilibrium surfacedensity τ is homogeneous, elevated ( τ = 3), as well as decreased ( τ = 0 .
5) within regions of radial width ∼
40 km.1 C. FIGURES OF SECTION 7.4 -100 0 100 200 300r-r L [km]0.60.81.01.21.41.61.82.0 τ β = 0.85 -100 0 100 200 300r-r L [km]0.60.81.01.21.41.61.82.0 τ β = 1.10 -100 0 100 200 300r-r L [km]0.60.81.01.21.41.61.82.0 τ β = 1.16 -100 0 100 200 300r-r L [km]0.60.81.01.21.41.61.82.0 τ β = 1.20 -100 0 100 200 300r-r L [km]0.60.81.01.21.41.61.82.0 τ β = 1.35 Figure 21 : Hydrodynamical integrations using the
P r β = 0 . − .
35 (see Figure 6) from top to bottom. The surface density profiles (left) as well as theirwavelet-powers (right) reveal co-existence of a resonantly forced density wave and short-scale viscous overstability forthe cases β = 1 . − .
35. All integrations use a scaled torque ˜ T s = 9 · − , a grid of size L r = 450 km and resolution h = 25 m. The plots correspond to times t (cid:38) ,
000 ORB.2
Lehmann et al. -100 0 100 200 300r-r L [km]0.60.81.01.21.41.61.82.0 τ T s = 1.00e-04 -100 0 100 200 300r-r L [km]0.60.81.01.21.41.61.82.0 τ T s = 4.00e-02 -100 0 100 200 300r-r L [km]0.60.81.01.21.41.61.82.0 τ T s = 9.00e-02 -100 0 100 200 300r-r L [km]0.60.81.01.21.41.61.82.0 τ T s = 1.60e-01 Figure 22 : Hydrodynamical integrations using the
P r T s = 10 − − . β = 1 .
25, a gridof size L r = 450 km and resolution h = 25 m. As in Figure 7 the blue dashed lines indicate the wavelength of vanishingnonlinear group velocity of overstable waves (by margins ±
20 m). The plots correspond to times t (cid:38) ,
000 ORB.3 ]0.00.10.20.30.4 e k i n [ σ c ] (axisymmetric)(non-axisymmetric) Figure 23 : Comparison of the nonlinear evolution of free ( T s = 0) viscous overstability in hydrodynamical integrationswith (non-axisymmetric, m = 7) and without (axisymmetric, m = 0) the azimuthal derivative terms (Section 5).Shown ae profiles of the surface density τ along with their wavelet powers for two different times ( t = 500 ORB and t = 35 ,
000 ORB ). The blue dashed lines indicate the expected nonlinear saturation wavelength of axisymmetric( m = 0) viscous overstability by margins ±
20 m (See Section 7.4.1). The red dashed curves represent the lineardensity wave dispersion relation (43). Note that in the axisymmetric case this curve has no physical meaning. Thebottom frame displays the evolution of the kinetic energy density for both integrations. The insert plot indicates thatthe linear growth phases ( t (cid:46)
200 ORB) of non-axisymmetric and axisymmetric modes are practically identical, inagreement with our considerations in Section 2. The higher saturation energy of the axisymmetric integration is dueto the slightly larger saturation wavelength (see Section 7.4.1 for explanations).4
Lehmann et al. -1.0 -0.5 0.0 0.5 1.0r-r L [km]02468101214 t i m e [ O R B ] q=0.0 -1.0 -0.5 0.0 0.5 1.0r-r L [km]02468101214 t i m e [ O R B ] q=0.1 -1.0 -0.5 0.0 0.5 1.0r-r L [km]02468101214 t i m e [ O R B ] q=0.4 Figure 24 : Space-time diagrams showing the linear evolution of the radial velocity field u (cid:48) of an initially seededtraveling wave in a ring with β = 1 .
35 perturbed by an ILR (at r = r L ). The perturbation, quantified by q , increasesfrom left to right. At initial time t = 0 the model ring is in the uncompressed state [Equation (58)] and the initialwavelength λ = 200 m. q = 0.0 -8-4048 M τ , M u t = 4.00 ORB -8-4048 M τ , M u t = 4.25 ORB -8-4048 M τ , M u t = 4.50 ORB -8-4048 M τ , M u t = 4.75 ORB -1.0 -0.5 0.0 0.5 1.0r-r L [km]-8-4048 M τ , M u t = 5.00 ORB q = 0.3 -6-3036 M τ , M u t = 4.00 ORB -6-3036 M τ , M u t = 4.25 ORB -6-3036 M τ , M u t = 4.50 ORB -6-3036 M τ , M u t = 4.75 ORB -1.0 -0.5 0.0 0.5 1.0r-r L [km]-6-3036 M τ , M u t = 5.00 ORB q = 0.4 -4-2024 M τ , M u t = 4.00 ORB -4-2024 M τ , M u t = 4.25 ORB -4-2024 M τ , M u t = 4.50 ORB -4-2024 M τ , M u t = 4.75 ORB -1.0 -0.5 0.0 0.5 1.0r-r L [km]-4-2024 M τ , M u t = 5.00 ORB Figure 25 : Snapshots of the terms M τ (cid:48) and M u (cid:48) that appear in the equation for the azimuthal velocity perturba-tion and which must be sufficiently in phase for the viscous overstability mechanism to work. The snapshots are fromintegrations with λ = 200 m, β = 1 .
35, and cover one orbital period in equal time-intervals. With increasing strengthof the satellite perturbation (quantified through the nonlinearity parameter q ) these terms become increasingly out ofphase. For q = 0 . − π occur, which explains the negative growthrates of overstable modes on all wavelengths (Figure 11, lower right panel). For clarity the two quantities have beenrescaled so as to posses equal amplitudes in all plots.5 D. METHOD B FOR THE AZIMUTHAL DERIVATIVESD.1.
Linear Waves
When restricting to linear density waves for the time being and adopting the notation (27), (28), then the azimuthalderivative of the solution vector can be written as ∂ θ Ψ | θ =0 = ∂ θ A τ ( r, t ) A u ( r, t ) A v ( r, t ) exp { i ( mθ − ω s t ) } θ =0 = − m τ I ( r, t ) u I ( r, t ) v I ( r, t ) + im τ R ( r, t ) u R ( r, t ) v R ( r, t ) , (D1)As a result the azimuthal derivatives (D1) induce a coupling between the real and imaginary parts of (1).D.2. Nonlinear Waves
For the description of nonlinear density waves (the amplitudes A τ , A u , A v are not small) for Equations (1) a splittingin real and imaginary part is not suitable. However, we can retain the description in terms of a coupled set of equationswhich can be seen as follows. First we note that performing the azimuthal derivative of the vector of state (28) isequal to a phase shift of π/ m , so that ∂ θ Ψ( r, θ, t ) = m Ψ( r, θ + π/ , t )= m [Ψ R ( r, θ + π/ , t ) + i Ψ I ( r, θ + π/ , t ) + c.c. ] (D2)This can also be expressed in terms of a time shift of P/
4, i.e. ∂ θ Ψ( r, θ, t ) = m Ψ( r, θ, t + P/ P = π Ω L .Inspection of the forcing terms (25), (26) shows that the imaginary part of the forcing function equals the real part,but with a phase shift of − π/
2. This means that we can consider the real and imaginary parts of (1) to describe thesame forced density wave, but with a phase shift of − π/ τ R ( r, θ, t ) u R ( r, θ, t ) v R ( r, θ, t ) = τ I ( r, θ + π , t ) u I ( r, θ + π , t ) v I ( r, θ + π , t ) . (D3)This observation is independent of whether the equations are linear or nonlinear. The idea is now to define two sets ofthe same nonlinear Equations (1), where one set is forced with the real parts of (25), (26) and is assumed to possessthe solution vector Ψ R ( r, θ, t ) = τ R ( r, θ, t ) u R ( r, θ, t ) v R ( r, θ, t ) . The other set is forced with the imaginary parts of (25), (26) and is assumed to possess the solution vectorΨ I ( r, θ, t ) = τ I ( r, θ, t ) u I ( r, θ, t ) v I ( r, θ, t ) . Lehmann et al.
The amplitudes A τ , A u , A v will now be affected by nonlinear terms in (1).Combining Equations (D2) and (D3) yields ∂ θ τ R ( r, θ, t ) u R ( r, θ, t ) v R ( r, θ, t ) = m τ I ( r, θ, t + π ) u I ( r, θ, t + π ) v I ( r, θ, t + π ) = − m τ I ( r, θ, t ) u I ( r, θ, t ) v I ( r, θ, t ) and ∂ θ τ I ( r, θ, t ) u I ( r, θ, t ) v I ( r, θ, t ) = m τ R ( r, θ, t ) u R ( r, θ, t ) v R ( r, θ, t ) . Note that although we retain the notation with subscripts R and I , the interpretation of the expressions denoting realand imaginary parts is only valid in the linear regime. In the nonlinear case the two quantities Ψ R and Ψ R merelydescribe the same density wave up to a constant relative phase shift. E. WENO RECONSTRUCTION OF THE FLUX-VECTORThe computation of the flux derivative ∂ r F in (18) includes a splitting of F according to the method of Liou andSteffen (1993) which was also used in LSS2017 and a WENO reconstruction of its individual components. In shortterms the reconstruction is as follows. We have (Shu and Osher (1988)) ∂ r F ( r j ) = f j +1 / − f j − / h (E4)with the numerical flux f , implicitly defined through F ( r j ) = 1 h r j +1 / (cid:90) r j − / f ( ξ )d ξ (E5)so that (E4) is exactly fulfilled. In these expressions the subscripts j ± / r j ± h . Equation (E5) can be used to obtain interpolating polynomials for f at a given location r , since the nodalvalues of the physical flux F ( r j ) are known for all j . We denote the so obtained unique 5th-order accurate polynomialapproximation for the numerical flux values at half nodes (see Section 4.1.1 of LSS2017) by ˆ f (5) j ± / = f j ± / + O (cid:0) h (cid:1) where ˆ f (5) j +1 / and ˆ f (5) j − / use the 5-point stencils [ r j − , r j − , . . . , r j +2 ] and [ r j − , r j − , . . . , r j +1 ], respectively.The starting point of the WENO reconstruction is the replacement of ˆ f (5) j ± / byˆ f j ± / = (cid:88) k =0 w k ˆ f kj ± / (E6)where ˆ f kj +1 / = (cid:88) l =0 c kl F j − k + l are the (unique) third-order accurate polynomial approximations for f j +1 / using the three 3-point stencils[ r j , r j +1 , r j +2 ], [ r j − , r j , r j +1 ] and [ r j − , r j − , r j ], respectively (for f j − / these are shifted accordingly by − w k Equation (E6) does yield ˆ f j ± / = ˆ f (5) j ± / . The key point of thedecomposition (E6) is an adequate assignment of the weights w k so that these yield the standard 5th-order accurateLagrange interpolation ˆ f (5) j ± / wherever F behaves smoothly across the entire 5-point stencil. If, however, in someregion the solution vector contains a discontinuity in one of the three sub-stencils, the corresponding weight should7diminish in order to avoid spurious oscillations of the solution vector. We use the WENO-Z weights introduced byBorges et al. (2008) which yield improved accuracy near extrema, as compared with the original WENO weights (Jiangand Shu (1996)). This improved accuracy is important as we are modeling wave systems that exhibit a wide rangeof length scales, where the shortest length scales will span only several grid points, and where the state variables cancontain sharp gradients.Since we apply a splitting of the flux F → F + + F − so that ∂ ( F + / − ) /∂ U possess only non-negative/non-positiveeigenvalues, the reconstruction outlined above applies to f + j ± / , whereas f − j ± / is reconstructed using stencils that areshifted by +1 so as to ensure a correct upwinding (cf. LSS2017). F. WKB-APPROXIMATION FOR SELF-GRAVITYFor integrations of linear density waves one can implement the self-gravity terms that arise from the solution for theself-gravity potential φ d in the WKB-approximation [cf. Equation (30)] f d = f dR + i f dI = − ∂ r φ d ( r, t ) = i πGσ τ ( r, t )= i πGσ [ τ R ( r, t ) + iτ I ( r, t )]= − πGσ τ I ( r, t ) + i πGσ τ R ( r, t ) . (F7)In this approximation, the self-gravity force at a certain grid point is governed by the value of the surface mass densityat this particular grid point only. This implementation of the self-gravity force couples the real and imaginary partsof (1).An alternative way to implement the WKB self-gravity which does not induce an additional coupling between theequations and which turns out to work also in the nonlinear regime is derived from Equations (35), (45), (52) and(53a) in LSS2016. From these relations follows that the disk potential φ d ( l ) and radial velocity u ( l ) are related through φ d ( l ) = D (cid:15)r L Ω L u ( l ) where l = 1 , r − r L r L . The exact relation (in the rotating frame) is [cf. Equation (29)] φ d ( l ) = − D (cid:15)r L ω s − m [Ω − Ω L ] u ( l ) . (F8)If we assume this relation holds for all higher harmonics l = 3 , , . . . , we can write the self-gravity force as f d = − ∂ r φ d = 2 πGσ ω s − m [Ω − Ω L ] ∂ r u. (F9)For sufficiently linear waves the implementations (F7) and (F9) yield identical results. G. DERIVATION OF THE PERTURBED GROUND STATEWe start with Equation (53) describing an m -lobed fluid streamline and the expressions for the radial and azimuthalvelocities (Borderies et al. (1983)) u = Ω ae sin m ( φ + ∆) , (G10) v = r Ω [1 + 2 e cos m ( φ + ∆)] , (G11)where the latter expressions are valid in an inertial frame ( r, φ ) with ϕ = φ + Ω s ( t − t ). Radial compression of thering matter is described by J = ∂ a r = 1 − q cos ( mφ + m ∆ + ˆ γ ) (G12)8 Lehmann et al. with q cos ˆ γ = ∂ a ( ae ) , (G13) q sin ˆ γ = mae∂ a ∆ , (G14)where q is the nonlinearity parameter [Equation (54)] and ∂ a denotes the derivative with respect to a . This results inthe scaled surface mass density τ = 1 J . (G15)The linearized velocity fields near x = x are given by u ( x ) = u ( x ) + x [ ∂ r u ] x = x , (G16) v ( x ) = v ( x ) + x [ a∂ r ( v/r )] x = x , (G17)where ∂ r = (1 /J ) ∂ a by Equation (G12). For the radial velocity u we need to compute ∂ a u = sin m ( φ + ∆) [ ae∂ a Ω + Ω e + Ω a∂ a e ] + Ω ae cos m ( φ + ∆) m∂ a ∆= sin m ( φ + ∆) (cid:20) −
32 Ω e + Ω e + Ω q cos ˆ γ (cid:21) + Ω q sin ˆ γ cos m ( φ + ∆)= − Ω e m ( φ + ∆) + Ω q [sin m ( φ + ∆) cos ˆ γ + cos m ( φ + ∆) sin ˆ γ ]= − Ω e m ( φ + ∆) + Ω q sin [ m ( φ + ∆) + ˆ γ ] . (G18)Thus, we have u = [Ω ae sin m ( φ + ∆)] x = x + Ω xJ (cid:104) − e m ( φ + ∆) + q sin [ m ( φ + ∆) + ˆ γ ] (cid:105) x = x . (G19)For the azimuthal velocity v consider ∂ a ( v/r ) = ∂ a [Ω (1 + 2 e cos m ( φ + ∆))]= − a (1 + 2 e cos m ( φ + ∆)) + 2Ω [ ∂ a e cos m ( φ + ∆) − e sin m ( φ + ∆) m∂ a ∆]= − a (1 + 2 e cos m ( φ + ∆)) + 2Ω a [ q cos ˆ γ cos m ( φ + ∆) − q sin ˆ γ sin m ( φ + ∆)]= Ω a (cid:20) − − e cos m ( φ + ∆) + 2 q cos [ m ( φ + ∆)] (cid:21) . (G20)This leads to v = [ r Ω (1 + 2 e cos m ( φ + ∆))] x = x + Ω xJ (cid:20) − − e cos m ( φ + ∆) + 2 q cos [ m ( φ + ∆) + ˆ γ ] (cid:21) x = x . (G21)Mosqueira (1996) assumes that nonlinearity (i.e. q >
0) arises solely from an eccentricity gradient, implying ˆ γ = 0[c.f. (G13), (G14)], and that the eccentricity e vanishes at x = x . These assumptions are expected to be fulfilled inthe evanescent region of the density wave, close to the Lindblad resonance, i.e. for x (cid:46)
0. In the frame rotating withΩ( x = x )) we then have τ = 11 − q cos m ( φ + ∆) , (G22) u = Ω qx sin m ( φ + ∆)1 − q cos m ( φ + ∆) , (G23) v = −
32 Ω x − q cos m ( φ + ∆)1 − q cos m ( φ + ∆) . (G24)9If, on the other hand, we are in the density wave propagation region ( x >
0) we can assume that nonlinearity arisessolely due to the variation of the phase angle ∆, such that ˆ γ ∼ π/
2. Let us rewrite (G19) and (G21) using the firstline of (G18) and the second line of (G20) such that u = [Ω ae sin m ( φ + ∆)] x = x + (cid:34) sin m ( φ + ∆) [ ∂ a Ω ae + Ω e + Ω a∂ a e ] + Ω ae cos m ( φ + ∆) m∂ a ∆ (cid:35) xJ ,v = [ r Ω (1 + 2 e cos m ( φ + ∆))] x = x + (cid:20) − e cos m ( φ + ∆)) + 2Ω a [ ∂ a e cos m ( φ + ∆) − e sin m ( φ + ∆) m∂ a ∆] (cid:21) xJ . (G25)Since we now have mae∂ a ∆ (cid:29) ∂ a ( ae ) as well as mx∂ a ∆ (cid:29) xJ dominates for both velocities and we arrive (in the frame rotating with Ω( x = x ))) at τ = 11 + q sin m ( φ + ∆) , (G26) u = Ω qx cos m ( φ + ∆)1 + q sin m ( φ + ∆) , (G27) v = −
32 Ω x q sin m ( φ + ∆)1 + q sin m ( φ + ∆) , (G28)which is identical to (G22)-(G24), up to an irrelevant constant phase shift of π/π/