# Density Waves and the Viscous Overstability in Saturn’s Rings

###### Abstract

This paper addresses resonantly forced spiral density waves in a dense planetary ring which is close to the threshold for viscous overstability. We solve numerically the hydrodynamical equations for a dense, axisymmetric thin disk in the vicinity of an inner Lindblad resonance with a perturbing satellite. The spiral shape of a density wave is taken into account through a suitable approximation of the advective terms arising from the fluid orbital motion. This paper is a first attempt to model the co-existence of resonantly forced density waves and short-scale axisymmetric overstable wavetrains in Saturn’s rings by conducting large-scale hydrodynamical integrations. These integrations reveal that the two wave types undergo complex interactions, not taken into account in existing models for the damping of density waves. In particular it is found that, depending on the relative magnitude of both wave types, the presence of viscous overstability can lead to a damping of an unstable density wave and vice versa. The damping of viscous overstability by a density wave is investigated further by employing a simplified model of an axisymmetric ring perturbed by a nearby Lindblad resonance. A linear hydrodynamic stability analysis as well as local N-body simulations of this model system are performed and support the results of our large-scale hydrodynamical integrations.

^{†}

^{†}Corresponding author: Marius Lehmann

^{†}

^{†}marius.lehmann@oulu.fi

0000-0002-0496-3539]Marius Lehmann \move@AU\move@AF\@affiliationAstronomy Research Unit, University of Oulu, Finland

Astronomy Research Unit, University of Oulu, Finland

0000-0002-4400-042X]Heikki Salo \move@AU\move@AF\@affiliationAstronomy Research Unit, University of Oulu, Finland

## 1 Introduction

The Cassini mission to Saturn has revealed a vast abundance of structures in the planet’s ring system, spanning a wide range 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 with wavelengths of some . It is generally accepted that this periodic micro structure originates from the axisymmetric viscous overstability mechanism (Schmit and Tscharnuter (1995, 1999); 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 density waves propagate through the rings, as these are excited at radii where the orbiting ring particles are in resonance with the 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 in terms of hydrodynamic models (Goldreich and Tremaine, 1978a, b, 1979; Shu, 1984; Shu et al., 1985, 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 in the sense that the excited surface mass density perturbation is small compared to the equilibrium value. In this case the governing hydrodynamic equations can be linearized and as a consequence the density wave appears sinusoidal in shape.

The studies by Shu et al. (1985), Borderies et al. (1986) (\hyper@linkcitecite.bgt1986BGT86 henceforth) and Lehmann et al. (2016) (\hyper@linkcitecite.lehmann2016LSS2016 henceforth) considered the damping behavior of nonlinear density waves in a dense planetary ring, such as Saturn’s B ring. In this situation the surface density perturbation is of the same order of magnitude as the equilibrium value. Within a fluid description of the ring dynamics, the damping of a density wave is governed by different components of the pressure tensor. The model by Shu et al. (1985) computes the pressure tensor from the kinetic second order moment equations, using a Krook-collision term. The model predicts reasonable damping lengths of a density wave for assumed ground state optical depths (or surface mass densities) that do not exceed a certain critical 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 wave trains propagate with ever increasing amplitude and nonlinearity. That said, the model fails to describe the damping of nonlinear waves in dense ring regions with high mutual 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 collision frequencies 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, \hyper@linkcitecite.bgt1986BGT86 compute the pressure tensor from a fluid model (Borderies et al. (1985)), as well as by using empirical formulae, which yield the correct qualitative behavior of the pressure tensor in a dense ring with a large volume filling factor. The computed damping lengths for optical depths relevant to Saturn’s dense rings are fairly long and the authors suspect this to be a consequence of the fluid approximation.

Borderies et al. (1985) pointed out that density waves are unstable in a sufficiently dense ring (such as Saturn’s B ring), whereas they are stable in dilute rings of small optical depth. Schmidt et al. (2016) have shown that the instability condition of spiral density waves is identical to the criterion for spontaneous viscous overstability in the limit of long wavelengths. In \hyper@linkcitecite.lehmann2016LSS2016 we derived the damping of nonlinear density waves from a different view point compared to the approaches by \hyper@linkcitecite.bgt1986BGT86 and Shu et al. (1985). We considered the density wave as a pattern that forms in response to this instability, using techniques that are widely applied in the studies of pattern formation in systems outside of equilibrium (Cross and Hohenberg (1993)). Consequently, the wave damping is described in terms of an amplitude equation. The resulting damping behavior is very similar to what is predicted by the \hyper@linkcitecite.bgt1986BGT86 model.

While the models by \hyper@linkcitecite.bgt1986BGT86 and \hyper@linkcitecite.lehmann2016LSS2016 can predict steady state profiles of density waves in an overstable ring region, they do not take into account the possible presence of axisymmetric wave structures that arise in response to the viscous overstability, independent of the perturbing satellite potential. The possibility of co-existence of spiral density waves with smaller wavelength axisymmetric periodic micro structure was discovered by analyzing 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 of interactions between a density wave and the waves generated by the viscous overstability. In Section 2 we outline the basic hydrodynamic model equations. Section 3 explains the numerical scheme applied to perform large-scale integrations of the hydrodynamical equations. Sections 4, 5 and 6 discuss specific terms appearing in these equations that arise from the forcing by the satellite, the advection due to orbital motion of the ring fluid, as well as the collective self-gravity forces, respectively. Results of large-scale hydrodynamical integrations are presented in 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 \hyper@linkcitecite.bgt1986BGT86 and \hyper@linkcitecite.lehmann2016LSS2016 in a marginally stable ring. In addition, we present some illustrative examples of density waves which propagate through a ring region which contains sharp radial gradients in the background surface mass density. We then consider waves that propagate in an overstable ring. In order to facilitate an interpretation of the results from our large-scale integrations, we introduce a simplified axisymmetric model to describe the perturbation of a ring due to a nearby ILR. We perform a linear hydrodynamic stability 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.

## 2 Hydrodynamic Model

We apply the vertically integrated isothermal balance equations for a dense planetary ring (Stewart et al. (1984); Schmidt et al. (2009)

(1) |

in a cylindrical frame with origin at , rotating with angular frequency where denotes the radial location of a specific inner Lindblad resonance (ILR) with a perturbing satellite and

(2) |

with Saturn’s mass and the gravitational constant . In what follows we will also make use of the radial distance

(3) |

as well as its scaled version .

The quantity is the rings’ surface mass density and with the ground state surface mass density . The symbols , stand for the radial and azimuthal components of the velocity on top of the orbital velocity . Furthermore, is the pressure tensor (see below). The central planet is assumed spherical so that , the latter denoting the epicyclic frequency of ring particles. The rings’ ground state which describes the balance of central gravity and centrifugal force is subtracted from above equations and we neglect the large-scale viscous evolution of the rings which occurs on time scales much longer than those considered in this study.

We neglect curvature terms containing factors since these scale as compared to radial derivatives. Here denotes the typical radial wavelength of a spiral density wave near its related Lindblad resonance where . From all terms containing derivatives with respect to we retain only the advective terms arising from the Keplerian motion, i.e. the first terms on the right hand sides of Equations (1). All other -derivatives scale as compared to radial derivatives ( denoting the number of spiral arms of the density wave), i.e. the same as curvature terms.

Poisson’s equation for a thin disk

establishes a relation between the self-gravity potential and the surface density .

The viscous stress is assumed to be of Newtonian form such that in the cylindrical frame we can write

(4) |

It is thus completely described by radial gradients of the velocities , , the dynamic shear viscosity as well as the isotropic pressure (see below). The ratio of the bulk and shear viscosity is denoted by , which is assumed to be constant (Schmit and Tscharnuter (1995)). The isotropic pressure and the dynamic shear viscosity take the simple forms

(5) |

(6) |

In this study we assume , i.e. the equation of state for an ideal gas. The ground state pressure can be defined in terms of an effective ground state velocity dispersion such that (Schmidt et al. (2001))

(7) |

The ground state is characterized by , , , together with the parameters in (5,6).

We neglect azimuthal contributions due to collective self-gravity forces. This neglect is adequate as long as the exerted satellite torque is much smaller than the unperturbed viscous angular momentum luminosity of the ring. That is, the (self-gravitational) angular momentum luminosity carried by the wave is negligible compared to the viscous luminosity. The linear inviscid satellite torque deposited at the resonance site reads (Goldreich and Tremaine (1979))

(8) |

where

(9) |

and (Cuzzi et al. (1984))

(10) |

The viscous angular momentum luminosity in the unperturbed disk is given by (Lynden-Bell and Pringle (1974))

In addition it should be mentioned that we are not concerned with the long-term redistribution of ring surface mass density which occurs in response to the presence of very strong density waves (\hyper@linkcitecite.bgt1986BGT86) so that we assume as mentioned before.

For the sake of definiteness we will restrict to parameters corresponding to the Prometheus 7:6 ILR, located at in Saturn’s A ring. We take values of the rings’ ground state shear viscosity and surface mass density (see Table 2) that can be estimated from corresponding values obtained by Tiscareno et al. (2007) for this ring region. The nominal values of and correspond to values found in N-body simulations with an optical depth [see \hyper@linkcitecite.lehmann2016LSS2016 (Section 3)]. Besides the nominal values we will use a range of values for [Equation (6)] and also [Equation (8)], in order to explore a variety of qualitatively different scenarios for the damping of density waves. The adopted value for the ground state velocity dispersion is larger then what results from local non-gravitating N-body simulations for optical depths relevant to this study (e.g. Salo (1991)) but corresponds roughly to expected values for Saturn’s A ring from self-gravitating N-body simulations exhibiting gravitational wakes and assuming meter-sized particles (Daisaka et al. (2001); Salo et al. (2018)). Furthermore, the value is still small enough to ignore pressure effects on the density waves’ dispersion relation (Section 7.1).

Our hydrodynamic model exhibits axisymmetric viscous overstability on finite wavelengths if the viscous parameter exceeds a critical value. To see this, let us ignore the satellite forcing for the time being. Furthermore, let us restrict to short radial length scales so that can be considered a constant simply denoted by . We introduce axisymmetric oscillatory perturbations such that

(11) |

with complex oscillation frequency and real-valued wavenumber . Poisson’s equation yields the relationship

(12) |

for a single wavelength mode (Binney and Tremaine (1987)).

In the remainder of this section we apply the dimensional scalings such that time is scaled with and length is scaled with . Inserting (11) and (12) into (1) and linearizing with respect to the perturbations (the primed quantities), results in an eigenvalue problem

(13) |

for . This equation can be used to obtain the growth rate and oscillation frequency of a given mode. This procedure has been carried out in several papers [see Lehmann et al. (2017) ( \hyper@linkcitecite.lehmann2017LSS2017 hereafter) and references therein for more details] and will not be repeated here. We merely point out that the threshold for viscous oscillatory overstability, i.e. a vanishing growth rate and hence for a given wavenumber is obtained from the real part of Equation (13), which yields a critical value of the viscosity parameter

(14) |

where we define

denoting the inverse of the hydrodynamic Toomre-parameter. A full list of symbols is provided in Table 3.

tableHydrodynamic Parameters Prometheus 7:6 () [] 1.5 [] 1 4.37 0.85 350 [ 1.26 [ 4.56 [ 7200

## 3 Numerical Methods

For numerical solution of Equations (1) we apply a finite difference Flux Vector Splitting method employing a Weighted Essentially Non-Oscillatory (WENO) reconstruction of the flux vector components. The method is identical to that used in \hyper@linkcitecite.lehmann2017LSS2017, apart from the reconstruction of the flux vector.

We define the flux-conservative variables

so that Equations (1) can be written as

with the flux vector

and the source term

(15) |

In the last expression

is the viscous stress tensor with denoting the unity tensor.

If not stated otherwise we adopt periodic boundary conditions on a radial domain whose size we denote by . The domain is discretized by defining nodes () with constant inter-spacing . The source term (15) contains radial derivatives of the stress tensor which are evaluated with central discretizations of 12th order. The computation of the flux derivative includes a splitting of according to the method of Liou and Steffen (1993) and a WENO reconstruction of the individual components.

In short terms the reconstruction is as follows. We have (Shu and Osher (1988))

(16) |

with the numerical flux which is implicitly defined through

(17) |

so that (16) is exactly fulfilled. In these expressions the subscripts denote evaluations at radial locations . Equation (17) can be used to obtain interpolating polynomials for at a given location , since the nodal values of the physical flux are known for all . We denote the so obtained unique 5th-order accurate polynomial approximation for the numerical flux values at half nodes (see Section 4.1.1 of \hyper@linkcitecite.lehmann2017LSS2017) by where and use the 5-point stencils and , respectively.

The starting point of the WENO reconstruction is the replacement of by

(18) |

where

are the (unique) third-order accurate polynomial approximations for using the three 3-point stencils , and , respectively (for these are shifted accordingly by ).

For a particular choice of the “weights” Equation (18) does yield . The key point of the decomposition (18) is an adequate assignment of the weights so that these yield the standard 5th-order accurate Lagrange interpolation wherever behaves smoothly across the entire 5-point stencil. If, however, in some region the solution vector contains a discontinuity in one of the three sub-stencils, the corresponding weight should diminish in order to avoid spurious oscillations of the solution vector. We use the WENO-Z weights introduced by Borges et al. (2008) which yield improved accuracy near extrema, as compared with the original WENO weights (Jiang and Shu (1996)). This improved accuracy is important as we are modeling wave systems that exhibit a wide range of length scales, where the shortest length scales will span only several grid points, and where the state variables can contain sharp gradients.

Since we apply a splitting of the flux so that possess only non-negative/non-positive eigenvalues, the reconstruction outlined above applies to , whereas is reconstructed using stencils that are shifted by so as to ensure a correct upwinding (cf. \hyper@linkcitecite.lehmann2017LSS2017). Furthermore, the evaluation of the derivatives with respect to and the treatment of the radial self-gravity force appearing in the source term (15) will be discussed in Sections 5 and 6, respectively.

Due to the presence of the satellite forcing terms in (15) it turns out that the simple time step criterion arising from a one-dimensional advection-diffusion problem, which was used in \hyper@linkcitecite.lehmann2017LSS2017, is unnecessarily strict. This criterion reads

(19) |

where is identified with the maximal eigenvalue of the Jacobian

(20) |

of equations (16) for the whole grid and is to be identified with the maximal value of the coefficient in front of the term containing the second radial derivative in (1), which is

The three eigenvalues of (20) read

For most integrations presented in this paper the grid spacings are large enough so that the second term in (19) is by far the smallest and can take values down to some . We find, however, that time steps in the range are suitable for all presented integrations, indicating that the criterion (19) cannot be appropriate. We have checked for some integrations with strong satellite forcing that reducing the time step by a factor of does not lead to any notable changes. For later use we also define the mean kinetic energy density within the computational domain as

(21) |

tableList of Symbols (gravitational constant) (planet’s mass) (mass of perturbing satellite) (semimajor axis of perturbing satellite) (eccentricity of perturbing satellite) (linear inviscid satellite torque) (viscous angular momentum luminosity) (Kepler frequency) (Kepler frequency at ) (epicyclic frequency) (vertical frequency of ring particles) (radial coordinate) (scaled radial coordinate) (time) (satellite forcing frequency in the frame rotating with ) (satellite forcing frequency in the inertial frame) (complex frequency of overstable waves) (surface mass density) (scaled surface mass density) (dynamical optical depth) , (planar velocity components) (self-gravity potential) (satellite potential) (planetary potential) (effective isothermal velocity dispersion) (ground state kinematic shear viscosity) (isotropic pressure) (dynamic shear viscosity) (viscosity parameter) (constant ratio of bulk and shear viscosity) (pressure tensor) (inverse ground state Toomre-parameter) (phase variable of a fluid streamline) (nonlinearity parameter of a fluid streamline) (phase angle of a fluid streamline) (semimajor axis of a fluid streamline) (eccentricity of a fluid streamline)

## 4 Satellite Forcing Terms

For simplicity, we restrict to density waves that correspond to a particular inner Lindblad resonance^{3}^{3}3In the current
approximation a Lindblad
resonance coincides with a mean motion resonance. of first order, so that the forcing satellite orbits exterior to the considered ring portion. The
wave is
excited by a particular Fourier mode of the gravitational potential due to this satellite with mass and semi-major axis and reads (cf. Section
5 in \hyper@linkcitecite.lehmann2016LSS2016)

valid in an inertial frame denoted by . The symbol

is a Laplace-coefficient with

In the current approximation the forcing frequency reads

(22) |

with the satellite mean motion . Upon changing to the frame rotating with frequency , denoted by , we have

yielding the forcing frequency in the rotating frame

(23) |

where we used (22). Therefore, the radial forcing component appearing in Equations (1) reads

(24) |

Similarly, the azimuthal component is given by

(25) |

These terms are evaluated at .

## 5 Azimuthal Derivatives

The persistent spiral shape of a density wave is generated by the resonant interaction between the ring material and the perturbing satellite potential, as well as the collective self-gravity force. Since our integrations are one-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 without loss of generality. The information about the azimuthal structure of the density pattern (the number of spiral arms ) is then contained solely in the terms describing azimuthal advection due to orbital motion. In what follows we refer to these terms simply as “azimuthal derivatives“. Thus, the requirement is to prescribe proper values of the azimuthal derivatives at for each time step of the integration.

We again adopt the cylindrical coordinate frame of Section 2 which rotates with angular velocity relative to an inertial frame denoted by so that

If we linearize Equations (1) with respect to the variables , , and , so that we restrict these to describe linear density waves, it is possible to solve the equations in the complex plane by splitting the solution vector into its real and imaginary parts:

(26) |

An evolving -armed linear density wave can then be described through the complex vector of state

(27) |

with , and being complex amplitudes in this notation which depend on time and the radial coordinate. The time dependence of the amplitudes is generally much slower than the oscillatory terms and vanishes once the integration reaches a steady state. When inserted to the linearized Eqs. (1) we obtain two sets of three equations that are possibly coupled through the azimuthal derivatives and self-gravity, depending on the applied implementation (see Sections 5.2, 6). Note that in practice we will exclusively use the full nonlinear equations (1). For sufficiently small amplitudes , , the nonlinear terms in (1) are negligible and the equations are essentially linear.

In order to describe nonlinear density waves, it is necessary to make an approximation for the azimuthal dependency of the wave. To obtain such an approximation we assume that (27) holds also in the nonlinear case. We have found two different implementations for the azimuthal derivatives (simply referred to as Methods A and B, described below) that yield a stationary final state for Equations (1) in the nonlinear regime. It turns out that the application of Method A (Section 5.1) results in nonlinear wave profiles that agree better with existing nonlinear models and therefore the results presented in subsequent sections are based on integrations using this method. We additionally outline Method B since we have found that it works well in the weakly nonlinear regime. However, the reader can skip the corresponding subsection 5.2 if preferred. Note that for sufficiently linear waves, both methods are exact down to the numerical error.

### 5.1 Method A

One implementation of the azimuthal derivatives can be derived from the vector of state of the weakly nonlinear model of \hyper@linkcitecite.lehmann2016LSS2016 in the first order approximation, where contributions from higher wave harmonics are omitted. We have

(28) |

where , , and are scaled according to Table 1 in \hyper@linkcitecite.lehmann2016LSS2016 and is scaled according to Table 3. Furthermore, is given by Equation (23). This relation follows from Equations (45) and (99a) in \hyper@linkcitecite.lehmann2016LSS2016 if the corresponding expressions are written in the rotating frame and if they are not expanded to the lowest order in (unlike in \hyper@linkcitecite.lehmann2016LSS2016). Furthermore, small corrections due to pressure and viscosity are neglected.

From the solution of the Poisson-Equation we have (cf. Appendix B in \hyper@linkcitecite.lehmann2016LSS2016)

(29) |

where is assumed to be given by the unscaled first component of (28) and denotes the wavenumber of the density wave. Note that both sides of Equation (29) are understood to be real-valued since the -function takes different signs for the two complex conjugate exponential modes in (28). The azimuthal derivative is then given by

(30) |

if we assume . Equation (30) shows that the azimuthal derivative of is directly proportional to the self-gravity force, the computation of which will be discussed in Section 6.

The azimuthal derivatives of and can be directly obtained from (28) and read

As mentioned before, Equation (27) [and (28)] can only be used as an approximation for a nonlinear density wave. Assume that a nonlinear density wave is correctly described by an infinite series

(31) |

where the terms with describe the wave’s higher harmonics. It is not straightforward to quantify the error made with the approximation (27), but the contributions to Equations (1) arising from higher order azimuthal derivatives () should be small. Note that only the azimuthal derivative terms are affected by the applied approximation.

### 5.2 Method B

#### 5.2.1 Linear Waves

#### 5.2.2 Nonlinear Waves

For the description of nonlinear density waves (the amplitudes , , are not small) for Equations (1) a splitting in real and imaginary part is not suitable. However, we can retain the description in terms of a coupled set of equations which can be seen as follows. First we note that performing the azimuthal derivative of the vector of state (27) is equal to a phase shift of and a multiplication by , so that

(33) |

This can also be expressed in terms of a time shift of , i.e.

where .

Inspection of the forcing terms (24), (25) shows that the imaginary part of the forcing function equals the real part, but with a phase shift of . This means that we can consider the real and imaginary parts of (1) to describe the same forced density wave, but with a phase shift of so that

(34) |

This observation is independent of whether the equations are linear or nonlinear. The idea is now to define two sets of the same nonlinear Equations (1), where one set is forced with the real parts of (24), (25) and is assumed to possess the solution vector

The other set is forced with the imaginary parts of (24), (25) and is assumed to possess the solution vector

The amplitudes , , will now be affected by nonlinear terms in (1).

Combining Equations (33) and (34) yields

and

Note that although we retain the notation with subscripts and , the interpretation of the expressions denoting real and imaginary parts is only valid in the linear regime. In the nonlinear case the two quantities and merely describe the same density wave up to a constant relative phase shift.

## 6 Self-Gravity

### 6.1 Straight Wire Model

We use the same implementation of collective radial self-gravity forces as described in detail in \hyper@linkcitecite.lehmann2017LSS2017. 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 :

(35) |

where we defined . This relation (35) does not include the force generated by mass contained in the bin itself, which can be approximated through

If is periodic with period the sum (35) can be replaced by the convolution

(36) |

of with the force kernel, which reads

Equation (36) can then be solved with a FFT method. However, since the density pattern of a resonantly forced density wave is not periodic we need to pad one half of the array with zeros in order to avoid false contributions from grid points outside the actual grid (Binney and Tremaine (1987)), e.g. gravitational coupling of material inside the resonance with material at far positive distances from resonance across the boundaries.

### 6.2 WKB-Approximation

For integrations of linear density waves one can implement the self gravity terms that arise from the solution for the self-gravity potential in the WKB-approximation [cf. Equation (98) in \hyper@linkcitecite.lehmann2016LSS2016]

In this approximation, the self-gravity force at a certain grid point is governed by the value of the surface mass density at this particular grid point only. This implementation of the self-gravity force couples the real and imaginary parts of (1).

An alternative way to implement the WKB self-gravity which does not induce an additional coupling between the equations and which turns out to work also in the nonlinear regime is derived from Equations (35), (45), (52) and (53a) in \hyper@linkcitecite.lehmann2016LSS2016. From these relations follows that the disk potential and radial velocity are related through

where denote the first and second harmonics of these quantities. This relation holds to the lowest order in . The exact relation (in the rotating frame) is [cf. Equation (28)]

(37) |

If we assume that this relation is valid for all higher harmonics , we can write for the self-gravity force

(38) |

## 7 Results

### 7.1 Excitation of Density Waves

In this section we illustrate the excitation process of a resonantly forced density wave. We use integrations employing the -parameters (Table 2) with different values of the forcing strength to elucidate nonlinear effects. All integrations were conducted with , time steps , and spatial resolution . Furthermore, Method A for the azimuthal derivatives (Section 5.1) and the Straight Wire self-gravity model (Section 6) were employed.

It is expected that during the excitation process the envelope of a density wave evolves in radial direction with the local group velocity. For a linear density wave described by the perturbed surface density

with wavenumber , one obtains in the frame rotating with frequency the dispersion relation (Goldreich and Tremaine (1978a); Shu (1984))

(39) |

Taking the derivative with respect to on both sides and re-arranging terms yields the group velocity (Toomre (1969))

(40) |

where is given by (23). By defining

and expanding this expression about the Lindblad resonance (using the approximation ) so that with given by (10), one obtains from (39) the wavenumber dispersion for linear density waves

(41) |

where is given by (9) and where the effects of pressure, expressed through the term quadratic in in (39), are ignored.

An expression for the nonlinear group velocity can be obtained from the nonlinear dispersion relation of spiral density waves [e.g. Equation (87) of Shu et al. (1985)] in the WKB-approximation

(42) |

In this expression the contributions due to pressure and self-gravity are modified and depend on the nonlinearity parameter with (see Section 7.4.3 and references therein). The integral

describes the nonlinear effects of self-gravity (Shu et al. (1985)) and, similarly,

describes nonlinear pressure effects. The resulting group velocity reads

(43) |

The integral functions fulfill and for . For realistic values of the velocity dispersion the nonlinear group velocity (43) is expected to be always larger than the linear limit (40).

Figures 7.1, 7.1 and 7.1 show stroboscopic space-time diagrams [time-resolution of ] of integrations with scaled linear satellite torques of , and , where corresponds to the nominal forcing strength for the Prometheus 7:6 ILR (Table 2). In these figures the gray shading measures the value of so that brighter regions correspond to larger values of . As the satellite forcing is taken to be spatially constant, initially the hydrodynamic quantities , and perform uniform oscillations (with infinite wavelength). Due to the Keplerian shear the pattern starts to wrap up at a constant rate. As the wavelength of the pattern decreases with time, at a certain radial location and at a certain time the wavelength will fulfill the dispersion relation (42) [and also (39) if is sufficiently small]. As soon as this is the case, the wavelength is “locked” to this value. In the figures, the region which becomes “locked“ in this way is enclosed by the dashed and solid blue lines. The former marks the resonance, while the latter is the predicted path of the wave front assuming it propagates with the linear group velocity (40). All wave structures outside this region eventually damp as they become increasingly wound up. An exception are axisymmetric waves generated by viscous overstability (Section 7.4). Also plotted are radial profiles of at four different times during the excitation process.

In Figure 7.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 differences when comparing the figures. First of all, with increasing torque value the wave profiles attain the typical peaky appearance of nonlinear density wave trains in thin disks (Shu et al. (1985); \hyper@linkcitecite.bgt1986BGT86; Salo et al. (2001)). Furthermore, the group velocity of the waves increasingly departs from the linear prediction (40), albeit mildly. One notes that there remains a very slow phase-drift of the wave pattern towards the resonance, indicating an increasing phase velocity with decreasing distance from resonance. Theoretically, at resonance the wavenumber of the density wave (41) vanishes so that the wave’s phase velocity diverges. It can therefore in general not be expected from a numerical method to correctly describe the wave pattern at the exact resonance location.

### 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 \hyper@linkcitecite.bgt1986BGT86 and \hyper@linkcitecite.lehmann2016LSS2016, which we refer to as the BGT and the WNL (Weakly Nonlinear) model, respectively. This section is restricted to stable waves in the sense that for all wavelengths [cf. Equation (14)], i.e. no overstability occurs in the system. All hydrodynamical integrations presented in this section were conducted with , time steps , and spatial resolution and used the -parameters (Table 2). If not stated otherwise, all integrations employed Method A for the azimuthal derivatives (Section 5.1) and the Straight Wire self-gravity model (Section 6). Presented BGT model wave profiles were computed using the method of \hyper@linkcitecite.bgt1986BGT86 (see their Section IVa), using the pressure tensor (4) with (5) and (6). Note that this model takes into account secular changes in the background surface mass density that accompany the steady state density wave in order to ensure conservation of angular momentum at all radii in the steady state. These modifications are neglected in our hydrodynamical integrations as well as the WNL model since the latter neglect the angular momentum luminosity carried by the density wave. To facilitate a comparison between the three different methods the profiles of resulting from the BGT model are scaled with the modified background surface mass density which will not be shown.

Figure 7.2 displays steady state profiles of the hydrodynamic quantities , , as these result from integrations with profiles computed with the WNL model (\hyper@linkcitecite.lehmann2016LSS2016). The profiles in the left and right columns result 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 are indicated by the fractional torque such that corresponds to the nominal forcing strength at the Prometheus 7:6 ILR and results in a nonlinear density wave train. The value corresponds to a weakly nonlinear wave. 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 reveals significant departures at larger distances from resonance between both methods, and Method A produces a clearly better match with the WNL model. All integrations presented in the following sections were conducted with Method A.

In Figure 7.2 we present wave profiles along with their Morlet wavelet powers (Torrence and Compo (1998)) for the case . Also for this strongly nonlinear wave, we observe an overall good agreement for both the amplitude profiles and wavenumber dispersions. Note that the WNL model takes into account only the first two harmonics of the wave [cf. Equation (31)], which is clearly visible in the wavelet power. This is also the reason why can take values below 0.5 (see \hyper@linkcitecite.lehmann2016LSS2016 for details).

Finally, Figure 7.2 compares profiles obtained from integrations with the Straight Wire self-gravity (left panels) and the WKB self-gravity using Equation (38) (right panels) for the cases (upper panels) and (lower panels). Comparison with corresponding BGT model wave profiles shows that the WKB-approximation is fully adequate for the weakly nonlinear wave with in that it yields indistinguishable results from the Straight Wire self-gravity. For the strongly nonlinear case the WKB-approximation has a notable effect. As expected, its application yields overall an even better agreement with the BGT (and WNL) model, apart from a few wave cycles at intermediate distances from resonance. A possible reason for this discrepancy could be that relation (37) might not hold for higher harmonics . In that case the resulting error should be largest where the wave’s higher harmonics take their largest amplitudes. We have verified that reducing the time step or the grid spacings by factors of does not change the outcome of all integrations presented in this section. The remaining differences between the integrated wave profiles and the model profiles are most likely due to the approximative implementation of the azimuthal derivatives. Nevertheless, the results presented here make us confident that our numerical integrations yield qualitatively correct behavior even of strongly nonlinear density waves.

### 7.3 Wave Propagation through Density Structures

In this section we present a few illustrative examples of hydrodynamical integrations of density waves propagating through an inhomogeneous ring. We restrict to the cases of jumps in the equilibrium surface density, but in principle we could also vary other parameters with radial distance, such as the viscosity parameter . All integrations adopted the -parameters and employed Method A for the azimuthal derivatives (Section 5.1) as well as the Straight Wire self-gravity model (Section 6).

Figure 7.3 shows space-time plots of a density wave passing a region of increased surface density (, left panel) as well as a region of decreased surface density (, right panel), in both cases of radial width . The sharp radial gradients in the equilibrium surface density, whose locations are revealed in the space-time plots, act like additional sources for the density wave in the sense that the wave profile can change at these locations prior to the expected arrival time of the wave front at these locations, the latter being indicated by the solid blue line (cf. Figures 7.1-7.1). It is, however, not clear how this is affected by the assumption imposed by our azimuthal derivatives that the hydrodynamic quantities describe an -armed pattern right from the start of the integration.

Figure 7.3 shows steady state profiles of along with corresponding wavelet-power spectra of density waves passing through regions of varying equilibrium surface density. For reference, the first row shows a density wave in a homogeneous ring. The second case, with a region of increased , bears some similarities with Figures 4 and 5 in Hedman and Nicholson (2016), showing profiles of the Mimas 5:2 density wave in Saturn’s B ring which passes through a region of radial width where the normal optical depth increases sharply from about 1.5 to values . In the region of enhanced surface density in Figure 7.3 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 homogeneous region. The two last cases represent situations with narrow regions of mildly decreased surface density and , respectively. In a region of decreased the wavenumber is enhanced, resulting in stronger wave damping. In the last case the wave damps out within this region. Note that the (viscous) time-scale on which the initially imposed sharp density gradients change in a notable manner, is much longer then the applied integration times.

### 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 subjected to viscous overstability such that [Equation (14)] for a non-zero range of wavelengths . Figure 7.4 displays the linear stability curve for the -parameters (Table 2) along with the different values adopted for the viscosity parameter [Equation (6)] in the integrations discussed in the following. Viscous overstability is expected to develop for all of the used -values, at least for a range of wavelengths, resulting in axisymmetric wave-trains as they are believed to produce parts of the observed periodic micro-structure in Saturn’s A and B ring (Thomson et al. (2007); Colwell et al. (2007); Latter and Ogilvie (2009)). For the three lowest values, linear viscous overstability is restricted to a relatively narrow band of wavelengths. In contrast, for the two largest values all wavelengths larger than a critical one are unstable. For these two cases it is expected from nonlinear models that a density wave retains a finite amplitude at large distance from resonance (\hyper@linkcitecite.lehmann2016LSS2016). However, these models do not take into account the presence of the waves generated by overstability which does not depend on the resonant forcing by an external gravitational potential. In this section we study the interplay of both types of structure. All large-scale integrations presented in this section were conducted using a grid with , and used Method A for the azimuthal derivatives (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.4.1 and 7.4.1 describe an integration using , without forcing by the satellite (). Figure 7.4.1 (left panel) shows a profile of (top) after about 20,000 orbital periods, along with its wavelet power (bottom). The structure on wavelengths represents the nonlinear saturated state of viscous overstability. This state consists of left- and right traveling wave patches, separated by source and sink structures (Latter and Ogilvie (2009, 2010); \hyper@linkcitecite.lehmann2017LSS2017). This can also be seen in the stroboscopic space-time diagram (right panel), showing the evolution of over 600 orbits in the saturated state within a small portion of the computational domain near the nominal resonance location. The green dashed lines represent the expected nonlinear phase velocity of axisymmetric viscous overstability (Figure 7.4.1, left panel) for a wavelength , with and denoting the nonlinear oscillation frequency and wavenumber of the overstable waves. Note that in Section 2 the symbol