Viscous overstability in dense planetary rings - Effect of vertical motions and dense packing

We investigate the linear axisymmetric viscous overstability in dense planetary rings with typical values of the dynamical optical depth 𝜏 ≳ 0 . 5. We develop a granular flow model which accounts for the particulate nature of a planetary ring subjected to dissipative particle collisions. The model captures the dynamical evolution of the disc’s vertical thickness, temperature, and effects related to a finite volume filling factor of the ring fluid. We compute equilibrium states of self-gravitating and non-self-gravitating rings, which compare well with existing results from kinetic models and N-Body simulations. Subsequently, we conduct a linear stability analysis of our model. We briefly discuss the different linear eigenmodes of the system and compare with existing literature by applying corresponding limiting approximations. We then focus on the viscous overstability, analysing the effect of temperature variations, radial and vertical self-gravity, and for the first time the effects of vertical motions on the instability. In addition, we perform local N-body simulations incorporating radial and vertical self-gravity. Critical values for the optical depth and the filling factor for the onset of instability resulting from our N-body simulations compare well with our model predictions under the neglect of radial self-gravity. When radial self-gravity is included, agreement with N-body simulations can be achieved by adopting enhanced values of the bulk viscous stress.


INTRODUCTION
The Voyager missions in the early 1980's, and particularly the recent Cassini (NASA/ESA) mission, revealed an intriguing amount of structure in Saturn's rings, including gaps, rings, and waves, on various length scales ranging from roughly 1000km down to about 100m.This structural richness is indicative of an underlying dynamical richness (Burns & Cuzzi 2006;Schmidt et al. 2009;Cuzzi et al. 2018) which has been exposed in the rings.Due to their proximity and consequent accessibility by spacecraft, Saturn's rings can serve as a "local laboratory" to study processes, such as accretion and disk-satellite interactions, which are also at work in protoplanetary disks (Cuzzi et al. 2010;Latter et al. 2018).Furthermore, it may be expected that future missions, such as the planned ring skimmer mission (Tiscareno et al. 2020) will provide observations of Saturn's rings with an even higher level of detail, which will help to further support the improvement of theories describing astrophysical disks.
★ E-mail: mlehmann@asiaa.sinica.edu.twHowever, the equations resulting from this formalism bear a great mathematical complexity such that they have so far not been applied to dynamical problems beyond the level of linear stability calculations.The main alternative approaches constitute either N-body simulations, or the application of hydrodynamic models.The limiting factor of N-body simulations is the entailing computational expense, especially when self-gravitational forces need to be included.Furthermore, it is often not straight forward to interpret the results from such simulations, and thereby to identify relevant physical processes at work.In order to get a handle on this, the development of simplified hydrodynamic models can be useful.Moreover, hydrodynamic simulations are typically less computationally expensive than N-body simulations, thus allowing to survey larger domains and timescales.On the other hand, hydrodynamic models neglect physical aspects that may or may not be important, depending on the problem at hand.In particular, hydrodynamic models typically assume a linear (Newtonian) relation between stress and strain in the fluid.Collisions between ring particles may not be sufficiently frequent to neglect the dynamical evolution of the stress, such that stress and strain are actually not instantaneously coupled (Latter & Ogilvie 2006), thus potentially violating the above assumption.Another issue is the presumed isotropy of the effective velocity dispersion tensor of the hydrodynamic model.When taking into account collisional contributions, the extreme thinness of Saturn's main rings implies that the vertical components of particles mutual impact velocities may be small, such that the vertical component of the collisional stress can be substantially reduced as compared with the planar components.Furthermore, hydrodynamic models require, depending on the number of equations, a number of a priori unknown parameters which need to be supplemented.
One quantity that has been found to be crucial to understand how the viscous overstability operates is the volume filling factor of particles (Latter & Ogilvie 2008;Mondino-Llermanos & Salo 2023).A proper description of variations in this quantity affords inclusion of variations in the disc thickness.Furthermore, hydrodynamic studies of the viscous overstability usually incorporate a bulk viscosity in the formulation of the viscous stresses, so as to mimic stabilizing processes not captured by a fluid description.Since bulk viscosity couples to compressional motions, and since it is known from theoretical studies (Borderies et al. 1985) and N-body simulations (Salo et al. 2001) of dense rings (with appreciable volume filling factor) that radial compression of particles results in vertical expansion, the effect of a bulk viscosity is potentially flawed if vertical motions are neglected.Another motivation for our new model is to move away from the rather vague and not wholly satisfactory viscosity prescriptions of previous work, which were usually in the form of power laws (Schmit & Tscharnuter 1995;Schmidt et al. 2001;Lehmann et al. 2017, to name a few).Connecting such models with N-body simulations (and real planetary rings) is not straightforward, though notable attempts to do so were undertaken (Salo et al. 2001;Schmidt et al. 2001) .Ideally, one would want a hydrodynamic model that is securely anchored in kinetic parameters, such as particle size, filling factor, etc., allowing direct comparison and generalisation from N-body simulations.
In this paper we make a first attempt to address the above points.To this end construct expressions for the pressure, the viscosity, the heat conductivity and the collisional cooling rate that are based on purely dimensional arguments as inspired by the study of Hwang & Hutter (1995), and which are framed in fundamental quantities such as the particle size, the particle mean free path, and the particle collision rate.This formulation facilitates a comparison with N-body simulations and kinetic models, but retains the mathematical simplicity of hydrodynamic models.To test predictions resulting from our model, we additionally perform local N-body simulations which include vertical and radial self-gravity forces.Overall we find good agreement between critical values of the volume filling factor and the optical depth for the onset of overstability resulting from both approaches.However, in order to obtain good agreement in the presence of radial self-gravity, we find it necessary to adopt values of the fluid's bulk viscosity that are significantly larger than those determined by Salo et al. (2001).
The paper is organised as follows.In §2 we briefly review the dynamics of dense planetary rings and aspects of the viscous overstability, as observed in Saturn's rings.In §3 we present our hydrodynamic model, along with the definition of the transport coefficients and all the parameters involved.In §4 we compute ring equilibrium states resulting from our model, which we compare to results form our N-body simulations.In §5 we perform a detailed axisymmetric linear analysis of the model.We consider various limiting cases and derive analytical expressions for the eigenmodes, which are compared to previous studies.In §6 we focus on the viscous overstability.
There we provide a detailed investigation of its driving and damping agents and present numerically obtained linear growth rates.In §7 we present our N-body simulations and compare critical quantities for the onset of viscous overstability between our model and simulations.In addition, we also compare linear growth rates of the viscous overstability following from both approaches.Finally, in §9 we summarise and discuss our main findings, as well as directions for future research.

DENSE PLANETARY RINGS AND THE VISCOUS OVERSTABILITY
While the kinematic equilibrium of a planetary ring is essentially the same as for any type of astrophysical disk, i.e. the balance between central gravity and centrifugal forces, the thermal equilibrium is rather distinct.That is, the latter is determined by the balance between heating due to viscous dissipation and cooling due to inelastic mutual particle collisions.

Viscosity in planetary rings
Viscosity (or angular momentum transport) in planetary rings essentially stems from either particle collisions or non-axisymmeric self-gravitational filamentary (wake) structures.
The contribution due to collisions is typically divided into two distinct parts, namely the 'local' and 'non-local' contributions  local and  nl , respectively (Schmidt et al. 2009).The local contribution is governed by random motions of particles in between mutual collisions, whereas the non-local contribution describes the angular momentum transferred between particles at contact during collisions.The latter becomes important as soon as the particle mean free path becomes comparable to the particle size (Shukhman 1984;Araki & Tremaine 1986).An expression for the local (kinematic) viscosity was derived by Goldreich & Tremaine (1978a): where  local is a dimensionless factor of order unity.On the other hand, the non-local (kinematic) viscosity scales as (Shukhman 1984;Araki & Tremaine 1986) In these expressions Ω denotes the orbital frequency,  is the particle random velocity ( velocity dispersion),  the particle diameter, and   the particle collision frequency.These expressions can also be derived from the basic expression for viscosity  ∼    2 , where  denotes the particle mean free path.In the case of the local viscosity  ∼ /  if multiple collisions occur during one orbital period, whereas  ∼ /Ω if collisions are rare, since ring particles are forced to evolve on epicycles.On the other hand, for the non-local mode  ∼ .Thus, their ratio is given by depending on the value of   /Ω.Here we defined known as the Savage and Jeffrey parameter (Savage & Jeffrey 1981;Araki 1991), which in Saturn's rings is assumed to be of order unity (e.g.Salo et al. 2018).Furthermore, in dense optically thick rings such as Saturn's A and B rings one has   ≫ Ω, such that the nonlocal angular momentum transport dominates the local contribution.However, even in optically thin regions of Saturn's rings, such as the C ring or the Cassini Division, where one expects   ∼ Ω, non-local transport effects are important.
Self-gravitational transport of angular momentum in Saturn's rings is realized through the emergence of self-gravity wakes, which are non-axisymmetric trailing density enhancements that form and dissolve on the orbital time scale (Salo 1992;Daisaka et al. 2001).This mode of angular momentum transport is expected to dominate throughout the A ring and possibly in some parts of the outer B ring.

Viscous overstability
It is common belief that some of the shortest scale structures which are detectable in Saturn's rings are most likely a result of the shortscale axisymmetric viscous overstability (Thomson et al. 2007;Colwell et al. 2007;Hedman et al. 2014).The axisymmetric viscous overstability is an intrinsic, oscillatory instability of the ring flow which draws energy from the Keplerian background shear and amplifies the collective epicyclic motion of the ring particles (i.e. a density wave) through a coupling of the viscous stress-oscillations to the background orbital shear.In a hydrodynamic description this mechanism requires that the viscous stress is a steeply increasing function of the surface mass density.Schmit & Tscharnuter (1995) studied the conditions for the onset of overstability in Saturn's B-ring in terms of an isothermal hydrodynamic model, assuming a simple powerlaw parameterization for the dynamic shear viscosity (5) They adopted a value  = 1.26 which was based on results from N-body simulations by Wisdom & Tremaine (1988).They found that the condition for the viscous overstability can be formulated as  >   with with  denoting the dynamic bulk viscosity.Their model assumed  =  so that   = 1/9.On the other hand, non-gravitating N-body simulations, as well as N-body simulations including vertical selfgravity (Salo et al. 2001) revealed a substantially larger critical value   ≈ 1 for the emergence of axisymmetric overstable oscillations.In terms of hydrodynamic modeling, this can be explained by including thermal effects (Spahn et al. 2000;Schmidt et al. 2001), as well as larger values of the bulk viscosity  ∼ (2 − 4)  (Salo et al. 2001).In non-self-gravitating N-body simulations (Salo et al. 2001) of metersized particles overstability was found for optical depths  ≳ 4.However, recent time-extended simulations by Mondino-Llermanos & Salo (2023) (MS23 henceforth) revealed the onset of overstability already for  ≳ 2.5, which agrees well with the dense ring kinetic model of Latter & Ogilvie (2008) (LO08 henceforth).On the other hand, in simulations by Salo et al. (2001) where vertical self-gravity was taken into account in terms of an enhanced vertical frequency of the particles Ω  /Ω = 3.6 (the same value was originally used by Wisdom & Tremaine 1988) overstability occurred for  ≳ 1.The reason for this difference is that the onset of the viscous overstability has been found to require a critical value of the volume filling factor FF =    (7) of the particle configuration (typically its value at the ring midplane) to be exceeded (LO08; MS23).In the above expression  and   denote the volume mass density and particle bulk density, respectively.In simulations of uni-sized particles which neglect vertical self-gravity, the particles tend to spread in vertical direction, thereby reducing the filling factor.This vertical extension of the ring is significantly reduced if vertical self-gravity is included.
The N-body simulations of Salo et al. (2001), and recently MS23, including particle-particle self-gravity revealed that axisymmetric overstable oscillations can co-exist with non-axisymmetric selfgravity wakes.But if the self-gravitational perturbations are too strong (resulting in too large density contrasts), overstability is suppressed.MS23 performed an extensive survey of gravitating N-body simulations and formulated a criterion for the onset of viscous overstability in terms of a critical filling factor, as well as a critical ratio between the viscous stress due to wakes (Daisaka et al. 2001) and the non-local viscosity.Latter & Ogilvie (2006) have shown that the viscous overstability is sensitive to any non-locality of the viscous stress in time.Thus, the instability mechanism requires that the stress oscillations are sufficiently in phase with the oscillations of the (epicyclic) velocity perturbations.This affords that the relaxation time of the stress (the collisional time scale) is appreciably shorter than the orbital time scale.Similarly, Lehmann et al. (2019) found that also the perturbation due to a nearby Lindblad resonance can destroy the phase relation of the viscous stress perturbation with the epicyclic oscillation in an overstable wave, offering an explanation for the mitigation of overstable waves by large scale density waves seen in their hydrodynamic simulations.
The long-term nonlinear saturation of the viscous overstability in Saturn's rings has, thus far, mainly been investigated in terms of isothermal hydrodynamic models neglecting self-gravity (Schmidt & Salo 2003;Latter & Ogilvie 2009, 2010).An exception are the nonself-gravitating N-body simulations by (Rein & Latter 2013) and the hydrodynamic simulations of (Schmit & Tscharnuter 1999) that took into account radial self-gravity.The result of these studies is that the viscous overstability typically saturates in form of traveling waves with wavelengths ∼ 100m − 500m.The papers by Latter & Ogilvie (2009), Latter & Ogilvie (2010) and Rein & Latter (2013) showed that the saturated state in general most likely consists of counterpropagating nonlinear wave trains, separated by defect structures (so called sink and source structures).More recently, Lehmann et al. (2017) performed N-body and non-isothermal hydrodynamic simulations that incorporated radial self-gravity as well as a constant vertical self-gravity (see §3.3).The key result of this study is that radial self-gravity determines the wavelength of nonlinear saturated wave trains.That is, under the influence of radial self-gravity the wave pattern tends to evolve into a state of minimum oscillation frequency, such that the pattern's wavelength is governed by the frequency minimum of the nonlinear dispersion relation of overstable waves.

Governing equations
We apply a local fluid model to describe a small patch of a planetary ring which orbits a planet of mass   at constant distance  0 , using the shearing sheet approximation (Goldreich & Lynden-Bell 1965).We adopt Cartesian coordinates (, , ), centered at radius  0 , and rotating around the planet with the Keplerian frequency such that  points radially outward,  in the direction of orbital motion, and  is parallel to the angular velocity vector given by ì Ω ≡ Ω 0 ì   .The flow of ring particles (a granular flow) is governed by the set of vertically averaged nonlinear fluid equations and In these equations   =   +     +     and   ,   are the radial and azimuthal components of the total velocity, such that the azimuthal component   includes the linear Keplerian shear − 3 2 Ω 0 .The central planet is assumed to be spherical so that we have equality between the orbital frequency Ω(), the epicyclic frequency and the vertical frequency of test particles.The balance between central gravity and centrifugal force in equilibrium has been subtracted from the radial momentum equation to leading order in the small parameter / 0 .A bar on a quantity  indicates vertical integration: As such, the surface mass density  is related to the volumetric mass density  = (, , , ) ( §3.2) via  = .Furthermore,  is the temperature, and  is the vertical thickness parameter of the disc.The temperature  is related to the root mean squared velocity deviations  of particles from the mean flow velocity field by  ≡  2 (e.g.Shu & Stewart 1985).The remaining quantities are specified below.
Eqs. ( 9)-( 14) can be obtained from vertical integration of the three-dimensional Navier-Stokes equations and the heat flow equation, thus describing the balance laws of mass, momentum, and temperature of the ring fluid.This is shown in detail in Appendix A.
The quantities q =   ,   and  denote the planar heat flux and the pressure tensor (see §3.4).The symbols ,  and Γ stand for the dynamic shear viscosity, dynamic bulk viscosity and cooling rate ( §3.5).The latter describes the loss of random kinetic energy due to inelastic particle collisions.We assume from the outset that the planar components of the velocity   ,   and the self-gravity force  sg,x , sg,y , as well as the temperature  are independent of  (see §3.3 and Appendix A).
Equations ( 13)-( 14) govern the dynamical evolution of the disc thickness parameter , where  describes the speed of vertical 'disc breathing'.For the vertical component of velocity we assume (Papaloizou & Lin 1988) such that the velocity   is the vertical velocity of a Lagrangian fluid element following the planar flow {  ,   }, and motions are assumed to be symmetrical with respect to the midplane  = 0 (see also Stehle & Spruit 1999).The latter restriction should be acceptable for our study of the viscous overstability in dense planetary rings.The terms on the right hand side of Eq. ( 13) in turn describe effects of vertical planetary gravity, vertical self-gravity and effects related to vertical stress, respectively.Furthermore,  is a dimensionless scaling factor defined by Our approach to include vertical disc motions into an otherwise twodimensional hydrodynamic model is very similar to the one of Stehle & Spruit (1999).Recently a more generalised approach, which allows for more complex vertical motions, has been devised by Ogilvie (2018).Furthermore, the kinetic model of Shu & Stewart (1985) allows for disc bending (or corrugation) motions, but assumes the disc to be in vertical hydrostatic equilibrium.Eq. ( 13) is also closely related to Equation (19) in LO08.The difference is that we allow the disc to depart from vertical hydrostatic equilibrium, in contrast to LO08, who assumed  = 0.As a consequence, the planar momentum equations ( 10)-( 11) and the heat flow equation ( 12) contain terms arising from vertical disc motions, resulting from viscous coupling and work done by pressure.The function  (  ) appearing in the equations for   ,   and  describes viscous coupling of vertical and planar motions.This quantity, and also the vertical stress   are given in §3.4.

Vertical density stratification
In order to evaluate the various terms appearing in Eqs. ( 9)-( 14) we need to make assumptions on the vertical structure of the planetary ring.For the purpose of illustration we will consider two limiting cases.In this work we will mainly apply a Gaussian distribution where the dependence on ,  and  is implicit via  and , which evolve according to (9) and ( 13), respectively.Strictly, this distribution applies in the dilute limit (Goldreich & Tremaine 1978a;Shu & Stewart 1985), characterized by a large values of /(Ω) (cf.§2).Using (18) the midplane volume filling factor (7) yields Furthermore, for the Gaussian distribution we find  = 1.
In addition, we consider a vertically unstratified (or uniform) slab of constant volumetric mass density By using (20) the (midplane) volume filling factor becomes Furthermore, for this distribution we find  = 1/3.A convenient aspect of the uniform distribution ( 20) is that all vertical integrals can be trivially solved, allowing to provide analytical expressions for all quantities appearing in our model.A uniform profile ( 20) is expected to apply to a cool, densely packed ring in which equal-sized particles have settled to a thin layer surrounding the midplane (Shu 1984;Borderies et al. 1985).
Note that the actual thickness of a particle configuration measured by √︃  2 yields  for the Gaussian, and / √ 3 for the uniform distribution.
Indeed, N-body simulations show ( §7) that the vertical distribution of a planetary ring transitions from a Gaussian toward a more or less uniform central distribution for increasing values of the optical depth (and hence ff).

Self-gravity
Since we study self-gravitating planetary rings, we additionally need to consider the Poisson equation with the self-gravitational potential  sg .From this equation one can derive the planar components of the self-gravity force  sg,/ = − /  sg as a function of the surface density, provided that one neglects any variations of  sg, / with the vertical coordinate  (Shu 1984, see also §5).This approximation is valid as long as radial variations in the density occur on much longer length scales than the vertical disc thickness.Nevertheless, corrections to the so obtained planar self-gravity forces can be applied to account for a finite vertical disc thickness ( §5).
The vertical component of self-gravity in Eq. ( 14) in general depends on the choice of the vertical disc structure and reads for the slab and Gaussian distribution, which follow directly from Eqs. (A36)-(A37), respectively.In order to arrive at these expressions we assumed that the self-gravitational potential  sg varies much faster in vertical than in planar directions, as well as symmetry of the ring around the mid-plane.
In addition, from Equation (13) we can define an effective dimensionless frequency of vertical motions which combines vertical planetary gravity and vertical self-gravity, neglecting collisions (contributions from the pressure tensor).This quantity which is often referred to as the 'Wisdom-Tremaine' factor features prominently in studies of self-gravitating planetary rings.A constant value of Ω  = 3.6Ω 0 was originally proposed by Wisdom & Tremaine (1988) to account for the effect of vertical self-gravity in N-body simulations of Saturn's B ring.The same value was also used in N-body simulations by Salo (1995), Salo et al. (2001), Rein & Latter (2013) and Lehmann et al. (2017).Recently MS23 found that the actual vertical forces in self-gravitating simulations are closer to Ω  ≈ 2Ω 0 .Moreover, they argued that an artificially increased value of Ω  can to some extent mimmick the inclusion of particle-particle surface (or frictional) forces which can promote viscous overstability in favor of self-gravity wakes (Ballouz et al. 2017).

Constitutive relations
The hydrodynamic equations ( 9)-( 14) need to be supplemented with constitutive relations for the viscous stress tensor and the heat flux q.
We assume a vertically isothermal planetary ring, such that  is independent of .Similar to previous studies (Spahn et al. 2000;Schmidt et al. 2001; LSS2017) we parameterize the heat flux as with the dynamic heat conductivity  (see §3.5.2).
For the pressure tensor we adopt the usual Newtonian form with the scalar pressure , the dynamic shear viscosity  ( §3.5.2) and the bulk viscosity  (see below).Furthermore, describes coupling of vertical and planar disc motions.Furthermore, we define the bulk viscosity as with the dimensionless parameter  (  ), which we formally let depend on the collision frequency (defined in Eq. ( 41) below).In our calculations below we find the bulk viscosity to have a significant impact on the stability boundary for viscous overstability.That is, too small values of the bulk viscosity result in viscous overstability for very small filling factors FF ≪ 1.This is not surprising, as our hydrodynamic model does not account for the non-Newtonian behavior of the viscous stress which prevails in dilute particulate rings with FF → 0 (Latter & Ogilvie 2006).The bulk viscosity in our model is used to mimic this, and other possible stabilizing processes that may occur in N-body simulations of planetary rings, but which are not explicitly captured in hydrodynamics.We find that a constant value of  0 = 4 (the equilibrium value) leads to a reasonable overall agreement between our linear calculations and results from N-body simulations presented below.A slightly better agreement is achieved when we use the relation of Salo et al. (2001) (their Figure 17, upper right panel).We then have such that  0 drops from a value of  0 ≈ 4 with increasing collision frequency, and levels of to  0 ≈ 2. The latter asymptotic value is similar to that derived from the Enskog kinetic theory of dense systems, which is  ≈ 1.3 (Chapman & Cowling 1970).A possible explanation for the rise of  0 at smaller collision frequencies may be the increasing non-locality in time of the viscous stress in simulations, as mentioned above.The above relation is used in all calculations below that involve the bulk viscosity.However, for very small collision frequency  ,0 ≲ 0.2 we find it necessary to apply a strongly increased value of  0 ≫ 1 to suppress spurious viscous overstability.

Particle collision model
We assume a planetary ring consisting of identical, non-spinning, spherical, indestructable particles with diameter  (radius   = /2) and bulk density   .The particles undergo partially inelastic mutual binary collisions and  is their mean separation distance.While the most accurate theoretical description of such a system is likely one based on a dense gas kinetic theory (see §1), in this paper we apply a granular flow model.That is, we describe the planetary ring through the equations of hydrodynamics ( 9)-( 14), but with expressions for pressure , shear and bulk viscosity  and , the heat conductivity , and the cooling rate Γ, which to some degree take into account the particulate and collisional nature of the system.Borderies et al. (1985) applied a granular flow model to Saturn's dense rings, based on an incompressible hydrodynamic model using transport coefficients derived from a simple kinetic model of binary particle collisions derived by Haff (1983).Using this model they were able to predict the onset of viscous overstability in a sufficiently dense planetary ring with optical depth  ≳ 1.This was subsequently confirmed in N-body simulations by Mosqueira (1996).Later, Hwang & Hutter (1995) revised the model by Haff (1983) and generalised it to account for two additional physical aspects of binary particle collisions which the original model omitted, and which may be of general importance to granular flows.For one, they allowed for variations of the volumetric mass density  of the particle ensemble, which requires the particle mean separation distance  to vary.Furthermore, they accounted for a finite contact time at collisions, related the particle material properties.
Since we adopt a compressible hydrodynamic model we follow Hwang & Hutter (1995), and allow for variations in the volume mass density.However, we neglect the finite contact time at particle collisions.We start by noting that the volumetric mass density may be written as where where   denotes the maximal density of the particle configuration under consideration.We assume which is approximately the maximal volume density of uni-sized randomly packed spherical particles.Based on Eq. ( 31) we define the effective filling factor in addition to the actual filling factor given by ( 7).The motivation is that the transport coefficients defined below should diverge if the densest possible packing is achieved, i.e. if  →   <   , rather than if  →   .Thus, the effective filling factor can take values 0 < ff < 1.
If we assume  ≪  (as was done in Haff 1983) we have  =   , consistent with an incompressible model.From (31) we directly obtain the particle mean separation distance where we used (33).Thus,  depends on the local state of 'packing' of the system.More precisely, (34) describes the average distance par-ticles traverse between two successive direct1 collisions with other particles.However, even in the absence of direct collisions, the particle mean free path in a planetary ring is limited from above by the epicyclic excursion length ∼ /Ω 0 , since particles are constrained to undergo epicycles rather than stream indefinitely (Goldreich & Tremaine 1978a).We formally write the particle mean free path as where we defined the modified filling factor with the dimensionless quantity  > 0, which ensures that the mean free path , unlike , does not diverge with vanishing density  → 0.
On the other hand, for  →   we find  →  → 0. Thus, Eq. ( 35) interpolates smoothly between these limiting cases, in a similar fashion as Eq. ( 1) for the local viscosity.The quantity  follows from solving yielding where we defined the parameter  ∼ 1, which will be stipulated below.
For definition of the transport coefficients below we first need to define the collision frequency of ring particles.From the time between two collisions we can define (Hwang & Hutter 1995) such that   → 0 for ff → 0 and   → ∞ for ff → 1.However, it is to be expected that this expression holds only for a sufficiently dense ring, since at small filling factor particle motions will be affected by epicyclic oscillations, as discussed above.Indeed, for ff → 0 Eq. ( 39) does not converge to the collision frequency of a dilute ring, which is given by (Shu & Stewart 1985): Recall that FF = 0.61 ff by Eq. ( 32).We note that simply replacing  in (39) by  does not resolve this discrepancy, as then the collision frequency does not vanish for ff → 0. Thus, we seek an expression for the collision frequency which converges to (39) and ( 40) for sufficiently large and small values of ff, respectively.We find that such an expression is given by which matches reasonably well the results from N-body simulations for small and large ff, as will be shown below.Note that for ff → 0 we have tanh ff → ff, and the bracket yields unity, such that   (ff → 0) becomes proportional to (40).For ff ≲ 1 the tanh-function varies slowly and   becomes essentially equal to (39).
Furthermore, we define the dynamical optical depth of the planetary ring as denoting the particle mass, as well as which is the scaled mutual Hill radius for a pair of particles.The latter is a convenient parameter to quantify the effect of self-gravity (S18; MS23).An increasing value of  ℎ implies an increased selfgravity force compared as with tidal forces.

Pressure, viscosity and cooling rate
Following Hwang & Hutter (1995) we now define pressure via where  0 is a dimensionless constant, formally of order unity, which will be specified below.In addition to the total pressure (44) we define its local contribution (cf.§2) as: which directly follows from (44) in the limit  ≫ .The latter is formally expected to apply if ff → 0, such that (44) becomes ( 45) for ff → 0 if the constant  0 is fixed properly.Thus, we expect our expression for the pressure (44) to be adequate throughout the dilute and dense regimes, i.e. for 0 < ff ≲ 1 (or 0 < FF ≲ 0.61).Based on (44) we can define an effective velocity dispersion (cf.Schmidt et al. 2001) which effectively includes local as well non-local contributions when used in the local pressure formula (45).Thus, the quantity  eff generally takes larger values than the actual velocity dispersion .
By following again Hwang & Hutter (1995) we can define the dynamic shear viscosity via  = "mass density" × "area" "time" with the dimensionless constant  0 specified below.Note that we here consider the time between direct particle collisions /, as the latter are required to produce a viscosity, in contrast to pressure defined above.However, it turns out that Eq. ( 47) is only adequate for sufficiently large values of ff, similar to (39).In order to correctly describe the viscosity for ff → 0 we therefore combine (47) with the classical formula (1) for a dilute ring, such that where we introduced the dimensionless coefficient  1 .Note that (47) differs from the classical formula for non-local viscosity (2) only by a factor of ff − 2 (48) converges to (1) for small ff, as it should.At this point it should be noted that neither (2) nor ( 47) are rigorously derived expressions for the non-local viscosity at large filling factor.Nevertheless, based on the above, we expect our expression for the viscosity (48) to behave at least qualitatively correct across the dilute and dense regimes, which will be confirmed below.Since heat conductivity and viscosity have the same dimensions, we simply define the dynamic heat conductivity via where  0 is again a dimensionless constant of order unity.Typical values  0 ∼ 1 − 5 were derived from N-body simulations (Salo 2001) applied to Saturn's B ring.In our stability analysis below we do not find any notable influence of heat conductivity on the linear viscous overstability.However, the inclusion of a finite value  0 ≳ 0.2 is found to be necessary to avoid spurious secular instability on length scales  ≳  0 .Therefore, we simply set in all calculations presented below which involve the heat flow equation (72).Finally, the cooling rate (or energy annihilation rate) is given by (Goldreich & Tremaine 1978c;Hwang & Hutter 1995 which describes the amount of specific kinetic energy lost via collisions, characterised by the normal coefficient of restitution  (see §3.5.3).Note that we again use the time between inter-particle collisions, which are required to dissipate energy.
If we assume the uniform density distribution (20) vertical integration of above transport coefficients yields as well as If we apply the Gaussian vertical distribution (18) the vertical integrals need to be computed numerically.

Coefficient of restitution
Finally,  is the normal coefficient of restitution, describing the irreversible loss of random kinetic energy of particles subsequent to a collision, as described by (51).In Saturn's main rings, collisions are expected to be rather dissipative (Porco et al. 2008), characterised by small values of .For simplicity, and in order to facilitate comparison with the recent N-body simulations of MS23 we will mainly work with a constant  = 0.1 − 0.5.Complementary, we also consider a velocity-dependent coefficient of restitution using the piece-wise powerlaw prescription (Bridges et al. 1984) which has been used in most kinetic studies and N-body simulations of Saturn's dense rings.Here we identify   ≡  as an averaged value of the mutual normal velocity between particles, and the scale parameter   , with nominal value   ≡   = 0.077mms −1 .For illustration purposes below we will also consider the case   = 10  , which results in a significantly hotter ground state.

Remarks
Our expresions for the pressure (44) and the viscosity (47) are analogous to Eqs. ( 46) and ( 47) in Borderies et al. (1985), except that we (as Hwang & Hutter 1995) do not assume  ≪ , and we use the (corrected) mean free path expression , rather than the mean separation distance .Furthermore, using (35) the requirement  ≫ , used to derive the local pressure (45), translates to the condition  ≪ 1.
As will turn out in our computations presented below, this is realized only if  ≲ 1. Physically, the fulfilment of the condition  ≪ 1 means that the epicyclic departure length of ring particles (Eq.37) should be much larger than the particle diameter .However, as outlined in §2, in Saturn's rings this condition is not expected to be fulfilled on account of their small velocity dispersion.In our calculations below we use  = 2, but generally values 1 <  < 4 result in reasonable agreement with results from N-body simulations.

EQUILIBRIUM STATE
In this section we derive the equilibrium state of a planetary ring corresponding to stationary solutions of Eqs. ( 9)-( 13) using the constitutive relations and transport coefficients described in the previous sections.

Numerical procedure
The ring equilibrium is characterised by  0 = const.,  0 = 0,  0 = −3/2Ω 0  and  0 = 0, where  0 is a free parameter.Note that since  0 is a constant, radial self-gravity forces vanish at equilibrium.Furthermore, values for the particle bulk density   , particle diameter  and the dimensionless viscosity coefficients  0 and  1 need to be specified.A list of important parameters and scales can be found in Table 1.In what follows we will work with dimensionless quantities, such that time is scaled with 1/Ω 0 , length is scaled with , and surface density is scaled with  0 .Throughout the remainder of this paper all values of   and  0 are understood to be in units of kg/m 3 and kg/m 2 , respectively.

Energy balance -equilibrium velocity dispersion
From Equation ( 12) at equilibrium we directly find the thermal equilibrium condition which describes the balance between viscous heating and collisional cooling and where Γ 0 and  0 are given by ( 53) and ( 55) at equilibrium, respectively.For the slab profile (20) this yields a quartic equation for the velocity dispersion, whose solutions can be obtained analytically, but are unwieldy.Instead, we consider some limiting cases.In particular, if the ring is hot and dilute with  0 ≫ 1 and ff ≪ 1, the viscosity (53) becomes equal to (1), and we find For ff ≪ 1 we have   ∼  (as shown below, see also Shu & Stewart 1985) and the above equation becomes the classical Goldreich & Tremaine  cr −  relation, where the unique solution  cr separates thermally stable and unstable equilibria.On the other hand, if we consider a dense ring ( 0 ∼ 1) in the limit of large collision frequency   ≫ 1 Eq.( 48) becomes equal to (47), and the energy balance (58) for the slab distribution readily yields such that now the parameter  0 controls the amount of (non-local) viscous heating, which is balanced by collisional cooling, quantified by .Note that the equilibrium temperature  0 ≡  2 0 .

Vertical hydrostatic balance -equilibrium scaleheight
Next, from (13) at equilibrium we readily find the equilibrium disc scale height, which for the slab profile reads where we defined the dimensionless self-gravity parameter and  0 denotes the pressure (52) at equilibrium.Thus, the equilibrium scaleheight is determined by the competition between vertical planetary gravity and self-gravity acting against pressure.For a Gaussian profile (18) we find Moreover, if we model vertical self-gravity through a constant value of Ω  (cf.( 25)), we find with  given by ( 17).In absence of vertical self-gravity ( → 0, or equivalently, Ω  → 1), the equilibrium scaleheight becomes  0 = √︁  0 /.Furthermore, as stated in §3.5, in the limit  0 → 0 we require Using ( 36)-( 38) this yields the relation The ring equilibrium is then computed as follows.We first specify values for  0 ,   , ,  ,  0 and  1 which are then kept fixed.Then, for sufficiently small  0 → 0 we assume initial values for  0 and  0 .From these we readily obtain  via (38) and (66).Using ( 33) and (36) we can evaluate the vertical integrals ( 52)-( 56) and apply a multi-dimensional Newton-Raphson method to solve (58) and ( 61) concurrently.We then use the so obtained values for ,  0 and  0,loc in subsequent iterations adopting increasing values of  0 .

Results
In Figure 1 we present various equilibrium quantities as functions of the optical depth , for non-gravitating ( ℎ = 0) and self-gravitating ( ℎ = 0.5) rings, where we consider a Gaussian vertical distribution (left panels) as well as a uniform vertical distribution (right panels).We use a constant  = 0.5 and adopt a distance  0 = 10 5 km from the central planet.The particle size is  = 2m and we use  0 =  1 = 0.053.The latter values are chosen as they yield good overall agreement with equilibrium properties inferred from N-body simulations, as can clearly be seen in the left panels.Note that the significant departure of quantities Ω  and  at  ∼ 0.5 is related to the formation of a layer structure in the simulations.That is, at  = 0.5 the system is has a central peak, whereas at  = 0.75 it starts to form two layers.
The displayed values of  0 in N-body simulations represent the root mean squared distance of all particles from the mid-plane: √︁ ⟨2 ⟩.
Likewise, the hydrodynamic results correspond to the value √︃  2 (cf.§3.2).The resulting values correspond to a ring which is a few particle radii thick, and in the non-gravitating case the scaleheight slightly increases with increasing optical depth.This increase is a consequence of the increasing non-local pressure (indicated by  eff ) with increasing .In the presence of vertical self-gravity the ring is subject to a significant flattening for optical depths  ∼ 1, but eventually thickens, again on account of increasing non-local pressure.
On the other hand, the velocity dispersion  is nearly independent of vertical self-gravity and in all cases steadily decreases, approaching the value Ω  in the optically thick limit  ≫ 1.Furthermore, the "viscous slope"  is defined (as in Salo et al. 2001) via where the subscript "0" is used to indicate that the derivative is computed (numerically) using the values of  corresponding to ring equilibrium states with increasing values of  0 .As such,  ≠ ( ln / ln ) 0 , since the latter derivative is understood to be taken using fixed values of  0 and  0 , whereas the latter two parameters also vary in between ground states with different  0 , as considered in (67) 2 .That being said, we find that  exceeds unity for values of the optical depth  ∼ 1 − 2, depending on the magnitude of self-gravity.A value  ≳ 1 is expected to trigger the onset of the viscous overstability ( §2 and §5).The flattening by vertical self-gravity results in an increased filling factor.Thus, non-local effects are augmented, yielding a larger viscosity and pressure.Consequently,  approaches unity at a significantly smaller optical depth  ∼ 1 than in the non-gravitating case where  ∼ 2, such that viscous overstability should be substantially promoted.
Comparing the results for the uniform and Gaussian distributions, we see that a Gaussian results in an overall more self-consistent equilibrium state in that all quantities agree well with the simulation results.Interestingly, the viscosity, which is the key quantity for our linear analysis below takes very similar values for the two vertical distributions.It should be noted that strictly, the dilute expression for the collision frequency (40) assumes a Gaussian vertical distribution.However, even if we correct this expression for a uniform distribution the overall agreement with N-body simulations is still significantly worse than for a Gaussian distribution.
In addition, Figure 2 shows the same quantities as in Figure 1 but for a constant  = 0.1 (left panels), as well as a velocity dependent coefficient of restitution (57) (right panels).In the case with  = 0.1 we used  0 =  0 = 0.047.Here the agreement between model and simulation is not as good as in the case with  = 0.5.The equilibrium following from the simulations is significantly cooler at small optical depths.
In the right panels we consider two values of the scale parameter   =   (black curves) and   = 10  (red curves).The former system is very similar to one with a constant  = 0.5, which can seen in the panel displaying  ().On the other hand, a system with   = 10  results in a substantially hotter system, such that the velocity dispersion  ≫ 1.Note that we computed the hydrodynamic curves using   = 0.5  and   = 15  , as these values yield better agreement with the simulations.Furthermore, we used  0 =  1 = 0.053 for the case   = 0.5  and  0 = 0.019,  1 = 0.094 for the case   = 15  (i.e. a significantly enhanced local viscosity).We can clearly see that the agreement between model and simulations is significantly reduced for the hotter equilibrium state.We again Figure 1.Ring equilibrium states following from our hydrodynamic model (curves) using a coefficient of restitution  = 0.5 are compared with corresponding results from N-body simulations (open circles), described in §7.Shown are from top to bottom panels the filling factor, the scaleheight, the kinematic viscosity and the velocity dispersion (first and third columns), as well as the collision frequency, the vertical self-gravity factor, the coefficient of restitution, as well as the viscous slope (second and fourth columns).The left panels correspond to a Gaussian vertical distribution (18), whereas the right panels assume a uniform density distribution (20).We compare ring states without ( ℎ = 0) and with ( ℎ = 0.5) vertical self-gravity.In the panel displaying the collision frequency   the dashed curves represent Eq. ( 39). 1, but now assuming  = 0.1 (left), and assuming a velocity dependent coefficient of restitution  (  ) given by (57) (right).The N-body simulations adopted the scale parameters   =   (red symbols) and  = 10  (black symbols), respectively.For the computation of the hydrodynamic equilibra we used   = 0.5  (red curves) and  = 15  (black curves), respectively.suspect that the main reason for this disagreement is the inaccuracy of the classical expression (1) for the local viscosity.Nevertheless, this is issue not of great concern in the present study as below we will focus on dense cold equilibrium states.

LINEAR ANALYSIS
Now that we have established the equilibrium state of the planetary ring, we turn our attention to the linear perturbation modes it supports.We analyze the stability of short-scale axisymmetric perturbations that we add to the equilibrium state, such that the perturbed variables read with real positive radial wavenumber  and complex eigenfrequency  =   +   , such that   and   denote the growth rate and frequency, respectively, and where the subscript "0" denotes values at equilibrium.

Linearised equations
Inserting ( 68) into ( 9)-( 13) and linearising with respect to the perturbations yields the dimensionless linear equations where we use the non-dimensionalisation as stated in §4 and where describing the effect of vertical self-gravity on perturbations of the ring thickness.The thickness correction factor  () in the radial self-gravity term is provided in Appendix C. Furthermore, with the derivatives   ≡  /, etc. .Here we used the fact that the spatial dependence of the transport coefficients ( 44)-( 51) can in general be written as  =  ((),  (),  (), ), and similar expressions for ,  and Γ.
Explicit expressions for the derivatives appearing in (77)-( 78) are displayed together with the equilibrium quantities  0 ,  0 and Γ 0 for illustration in Figure 3 for the equilibrium shown in Figure 1.The plots show that deviations between the two density models (20) and ( 18) are most pronounced for variations in the scaleheight  in both limits of small and large filling factors.
Eqs. ( 69)-( 74) can be written as where  is a 6 × 6 complex matrix and ì  = ,   ,   , , ,    is a corresponding eigenvector.In the general case, we solve this eigenvalue problem with standard packages in IDL.

Isothermal, vertically static limit
In this section we neglect variations of the temperature and the disc scale height from their corresponding equilibrium values and thus also discard ( 72)-( 74).We then arrive at the standard isothermal equations that have been analyzed in numerous previous studies (e.g. Ward 1981;Schmit & Tscharnuter 1995;Latter & Ogilvie 2009).The dispersion relation, which is given by Det [] = 0 (cf.80) in this limit reads giving rise to three eigenmodes, which are the viscous instability (VI) mode and a pair of overstable modes (OS).Since the radial wavenumber  ≪ 1 we can adopt a perturbation series where the   ( = 0, 1, 2) are complex.Inserting this in (81) and expanding the latter in orders of  yields the approximate eigenfrequencies Thus, the VI is triggered if   < 0. As discussed in §4.2, this is expected to occur in a dense planetary ring.Nevertheless, radial self-gravity amplifies the instability at shorter wavelengths, as described by the term ∝  3 .On the other hand, for the onset of viscous overstability we find the condition on sufficiently long wavelengths.For a viscosity of the form (5) Eq. ( 84) yields the condition  >   , where   is given by the expression (6) derived by Schmit & Tscharnuter (1995).Inspection of Figure 3 reveals that radial self-gravity is expected to increase the growth rates of the visous overstability via the  3 -term, at least for sufficiently large ff, since then 3  >  0 .It can also readily be seen that self-gravity and pressure decrease and increase the modes oscillation frequency, respectively, at sufficiently short wavelengths.This gives rise to a frequency minimum in -space, which has been found to play a role in the nonlinear saturation process of the viscous overstability (Lehmann et al. 2017).

Vertically static limit
Next, we consider the vertically static limit, such that thickness variations are neglected (but temperature may depart from its equilibrium value).This limit has been investigated in several studies (Spahn et al. 2000;Schmidt et al. 2001;Lehmann et al. 2017).
To leading order, the VI mode is now described by  VI,hstat = −3 where for the second equality we used (58), as well as In addition, we find the thermal mode where with describing the temperature dependence of the thermal equilibrium (58), as defined in Schmidt et al. (2001).The isothermal limit can be recovered for Γ  → ∞, implying that any temperature variation is instantly extinguished.In this limit ( 85) becomes ( 82).Furthermore, the viscous overstability is described by the pair with (Schmidt et al. 2001) Thus, viscous overstability now occurs (in the limit of long wavelengths) if For the equilibrium shown in Figure 1 we find that  3 < 0 decreases for all ff such that, in agreement with previous studies, thermal effects reduce growth rates of the viscous overstability (cf.( 84)), and in addition increase the critical value of ff for its onset.On the other hand,  2 > 0 for ff ≲ 0.6, but turns negative for larger ff.Thus, for most feasible values of ff thermal effects enhance the oscillation frequency of overstable modes (see also Figure 4 below).

Isothermal limit
Finally, we include variations of the ring thickness but neglect temperature variations, such that the dispersion relation is of fifth order in the complex frequency .The corresponding expression is lengthy and will not be displayed here.To leading order in , the VI mode is now described by where the second term in the bracket encompasses effects of ring thickness variations and vertical self-gravity.Based on Figure 3 we find that vertical self-gravity has an amplifying effect on the VI, whereas the other term ∝ (   −  0 ), describing non-local pressure effects, has a damping effect.
In addition, we find the pair of vertical 'breathing' modes It can readily be seen that in the here considered limit of long wavelengths viscosity (as expected) damps these modes, and in addition, decreases their frequency.In contrast, pressure forces increase the modes' frequency, since   < 0 (Figure 3).In the dilute limit (  → 0, as well as vanishing viscosity) we find  H,iso → ± √ 2, consistent with Eq. ( 31) in Lubow (1981).
The corresponding eigenvalues describing the viscous overstability modes in the isothermal limit  OS,iso can be written in the same form as ( 90), but with  2 and  3 replaced by corresponding expressions  4 and  5 , respectively.Compared to the "thermal" expressions ( 91)-( 92), the expressions  4 and  5 are rather unwieldy and will not be displayed here.Nevertheless, the functions  4 ,  5 are plotted for the non-gravitating case  = 0, together with their "thermal" counterparts  2 ,  3 in Figure 4 .We find that compared to the effect of temperature variations, the effect of vertical disc motions (i.e.scaleheight variations) is somewhat more complicated.That is, based on the behavior of  5 the latter can in principle either reduce or increase the critical value of   (or equivalently FF) for the onset of viscous overstability, depending on the critical FF resulting from (84).Furthermore, scaleheight variations are likely to reduce the overstable modes' oscillation frequencies, unless FF  resulting from ( 84) is very small.This is opposite to the effect of temperature variations.However, at larger values of FF, scale height variations are expected to have a damping effect on overstable modes, as the function  5 becomes positive at a given ff.We find that it is the bulk viscosity which is responsible for  5 taking positive values.For  0 ≲ 7/3 we find that  5 ≲ 0 for all FF.It should be noted though that overall a bulk viscosity has a damping effect (cf.Eq. ( 93)).Moreover, vertical self-gravity mildly increases both  4 and  5 (not shown), such that  4 can attain positive values for a larger range of FF.

VISCOUS OVERSTABILITY
In this section we focus on the viscous overstability.We assume the dense ring equilibrium as shown in Figure 2 with constant  = 0.1.For these parameters (see also §5), the viscous overstability is the only instability that can occur.If not stated otherwise, presented results correspond to a Gaussian vertical density distribution.Moreover, in addition to the solution of the full system of equations ( 69)-( 74), we will consider various limiting cases as in §5.2- §5.3, which will allow us to clarify the role of different physical elements in the viscous  91) and ( 92), represented by diamond symbols.Also plotted are the coefficients  4 and  5 related to vertical motions in absence of self-gravity, as described in the text.The solid (dashed) curves apply to a Gaussian (uniform) distribution ( 18).The corresponding profiles of  2 and  3 for a uniform distribution ( 20) are identical to the Gaussian case.overstability mechanism.In this section we assume a constant bulk viscosity  0 = 4.

Driving and damping agents
In this section we neglect self-gravity, but briefly consider its effect in §6.2.In absence of self-gravity the equilibrium and linear stability properties of the planetary ring can be parameterised through a single parameter, e.g. the filling factor FF (its mid-plane value) or the optical depth  for identical particles with a given constant .This can be seen from Eqs. ( 69)-( 74).Using  = 0 the above equations depend on ff,  0 , and  0 .In absence of self-gravity the latter two quantities are fully determined by the value of , or, interchangeably, FF.The reason is that in this case the equilibrium filling factor (7) and optical (42) depth differ by a factor of the scaleheight  0 (64) which, in turn, only depends on either  or FF.This is not true in the self-gravitating case, where the equilibrium scaleheight  0 (and also the temperature  0 ) explicitly depends on the particle bulk density   or equivalentlty  ℎ (cf.Eqs. ( 61), (63) as well as Figure 1).
Figure 5 (upper panel) shows the maximal growth rate of the viscous overstability (across all wavenumbers ) for increasing midplane filling factor FF(0), assuming a Gaussian density distribution, and for the different approximations as discussed above.These are "hstat" (vertically static), "iso" (isothermal), and their combination.The label "full" corresponds to the solution of the full set of Eqs. ( 69)-( 74).The second panel shows the corresponding wavelengths of the fastest growing modes.These plots readily reveal the values of the critical filling factor FF cr for the onset of instability in the different cases, indicated by the sharp "rise" of the wavelengths of growing modes as functions of FF.The remaining panels show the pseudo-energy decomposition (Appendix B) for the different cases.
The energy decomposition confirms that the driving force of the viscous overstability is in general the increase of the viscous stress with increasing surface mass density, described by   , in agreement with Eqs. ( 84) and ( 93).This quantity is positive, increases with increasing FF (Figure 3), and is closely related to the viscous slope  defined by (67).From Eq. (B7) we see that instability then requires a phase shift  <   −    < 0 between the surface mass den- sity perturbation and the azimuthal velocity perturbation.According to Figure 6 (upper panel), where we display phase shifts between individual hydrodynamic quantities corresponding to the overstable modes resulting from the full solution as shown in Figure 5 (for a left traveling wave), we find   −    ≈ /2 in all cases, which is optimal for instability.
As explained above, the critical value of   corresponds to a unique critical filling factor ff cr in our model.The requirement of a critical filling factor to be exceeded for the onset of instability is in agreement with the kinetic model of LO08 and the N-Body simulations of MS23, and is therefore considered to be the most fundamental criterion for the viscous overstability.
The largest damping effect on overstable modes in the full model we find to be related to variations of the viscous stress with ring thickness , described by the curve labeled   .This damping can easily be understood when comparing (B7) and (B8), since variations of the surface density  and the scaleheight  are practically in phase in overstable waves (See figure 6).Thus, the increase of the viscous stress with increasing  is partially compensated by a corresponding decrease due to an increasing scaleheight.Intuitively, a damping related to ring thickness variations is expected, since a vertical expansion of ring material reduces the local volume filling factor, such that viscous overstability is expected to be weakened.This is a novel finding, as previous studies on the viscous overstability did not include the effect of vertical motions.
In agreement with our calculations in the limit of long wavelengths ( §5) we find that the inclusion of temperature variations increases FF cr , which can be seen by comparing the panels labeled "iso&hstat" and "hstat".A closer inspection of the energy-decomposition (not shown) reveals that the largest difference between the two cases lies in the term related to   , which results in a larger pressure-related damping effect compared to the isothermal limit.On the other hand, the contribution due to   is found to be subdominant.Interestingly, the inclusion of variations of the scaleheight (i.e.vertical motions) slightly decreases FF cr (seen by comparing the plots labeled "iso&hstat" and "iso").However, this reduction is not seen when comparing the cases "hstat" and "full", where we find only a very small difference.Thus, it appears that thermal effects and the effects of vertical motions are not independent.
The vertical dashed lines in the three lowest panels of Figure 5 indicate the critical values for ff cr , obtained from (84) for the case "iso&hstat", (93) for the case "hstat", and (93) with  3 replaced by  5 (see Figure 4) for the case "iso".These are in excellent agreement with the numerical results.In Table 2 we summarise critical values of different quantities for both Gaussian and slab distributions and for the different limiting cases.We find in general quite notable differences between the two density distributions, which we will, however, not further discuss here.Note also that each of the individual approximations strictly requires a different prescription of the bulk viscosity  0 , which is ignored here.
In accordance with previous studies, the longest wavelengths become unstable first, albeit with negligible growth rates.Overall, the typical wavelengths of modes (with meaningful growth rates   ≳ 10 −4 ) are  ∼ 100m, a property shared with previous hydrodynamic models (Schmit & Tscharnuter 1995;Schmidt et al. 2001;Latter & Ogilvie 2009;Lehmann et al. 2017), kinetic modeling LO08 and N-body simulations (Salo et al. 2001;Lehmann et al. 2017;MS23).On the other hand, the behaviour of the maximal growth rates with increasing filling factor is less obvious and does not seem to correlate with the corresponding wavelengths of the modes.

Vertical self-gravity
With the inclusion of vertical self-gravity, the characterisation of the equilibrium and, in principle also the perturbation quantities requires one additional parameter.Here we choose to work with the optical  depth  (replacing FF) as well as  ℎ , the latter quantifying the strength of self-gravity.Therefore, in Figure 7 we display similar quantities as in Figure 5, except using  as free parameter (the bottom frame assumes a fixed  ℎ = 0.5).Note that  0 (or ) enters the linear problem via the self-gravity parameter .The upper frame displays the growth rates as a contour plot of both  0 and  ℎ .
The critical filling factor FF  for the onset of overstability significantly reduces as 0.52 ≳ ff cr ≳ 0.37 with increasing 0 <  ℎ ≲ 0.6, such that vertical self-gravity appears to have a significant effect on overstable perturbations.The pseudo-energy decomposition reveals that vertical self-gravity affects the   0 and    -terms.The former becomes less damping, whereas the latter -in contrast to the non-selfgravitating case -becomes slightly amplifying for a range of sufficiently small filling factors.Considering Eq. (B9) and using the fact that   > 0, this implies that the relative phases between the perturbations  and   are affected by vertical self-gravity, which is indeed confirmed in Figure 6 (by comparing the upper and lower panels).Furthermore, according to Eq. (B4) the ability of vertical self-gravity to directly amplify a particular eigenmode depends on the relative phase shifts of the perturbations  and  , such that a phase shift ± is optimal, whereas a phase shift of ±/2 erases the direct effect of vertical self-gravity.Figure 6 shows that for FF ≳ FF cr the phase shift   −   ≲ −/2, but it approaches −/2 for increasing FF, such that the effect of vertical self-gravity leads to a mild reduction of FF cr , but its effect eventually vanishes at larger FF.We find (not shown) that the aforementioned indirect effect of vertical self-gravity through the    -term dominates the direct effect.It should be noted though that the main effect of vertical self-gravity on the region in ( −  ℎ )-space that is overstable enters via the altered equilibrium state (i.e. via a reduction of the equilibrium FF).

Radial self-gravity
Compared to the impact of vertical self-gravity, we find that radial self-gravity has a rather strong effect on overstable perturbations, with overall substantially larger growth rates up to an order of magnitude higher.For  ℎ ≳ 0.5 we actually find overstability for all FF ≳ 0, which is related to the "spurious branch" of overstable modes at small filling factor ( §3.4), which starts to merge with the "physical branch" of overstable modes (cf. Figure 7).We will come back to the effect of radial self-gravity on overstability in §8, where we compare results from our model with those from our N-body simulations.

N-BODY SIMULATIONS
In addition to our hydrodynamic analysis, we perform N-body simulations of a small patch of a dense planetary ring, using the lo-cal simulation method (Wisdom & Tremaine 1988).The simulation method is the same as described in Lehmann et al. (2017) (see also Salo 1992;Salo 1995, with slight modifications to the computation of vertical and radial self-gravity forces, as outlined below.Note that the azimuthal component of self-gravity is omitted.This is done so as to facilitate comparison with our axisymmetric linear hydrodynamic model for the viscous overstability without complications from non-axisymmetric self-gravity wakes. We perform a series of simulations with constant  = 0.1 (for some cases we use  = 0.5), at  0 = 100, 000km in order to determine the threshold (or critical) values of the optical depth and the volume filling factor (Mondino-Llermanos & Salo 2023) for the onset of viscous overstability, where   () describes the area of particle  cut by the plane at height .The filling factor (97) represents an average value over several orbital periods and the sum is taken over all particles.Note that prior to the onset of instability the system is uniform in planar directions.In figures below we present the mid-plane value of (97) at  = 0. We consider different strengths of self-gravity, measured by the  ℎ parameter (43).We adopt the values  ℎ = 0.61, 0.5, 0.4, 0.3, 0.2, for a range of optical depths around the threshold.For the applied distance  0 these correspond to particle bulk densities 13kgm −3 ≲   ≲ 350kgm −3 .In all simulations the radial size of the box and the number of particles are the same, i.e.   = 4776 and   = 300  , where   = 1m.Also the number of radial bins is  max = 64 in all cases.These fixed values assure that the accuracy of the statistical sampling of forces, vertical profiles, etc., is similar in all cases.The simulations are evolved for 300 orbital periods.In each of the  ℎ -series different optical depths corresponds to a different tangential size   = 25  / particle radii.Furthermore, the surface mass density  0 , is determined by the prescribed values of  and  ℎ via ( 42) and ( 43): ℎ ( 0 /100, 000 km) 3 . (98)

Calculation of self-gravity forces
The radial self-gravity calculation is made with the FFT-method as in Salo & Schmidt (2010) and Lehmann et al. (2017), except that here we correct radial self-gravity for the finite thickness of the ring in a similar fashion as in our hydrodynamic model (see §5).The simulation region of radial and tangential size   and   is divided into  max radial bins.Then, using a radial Fourier decomposition of the tangentially averaged surface density the corresponding radial force is where we now add the thickness correction factor  (  ), as given by (C8) and with We note that it is essential that this thickness correction factor is used in our simulations.That is, without the suppression of very short wavelength modes, overstability would grow much more vigorous.
In above expressions  0 is the mean surface density of the simulation region,   and   are the fractional amplitude and phase of different -components with wavenumber   (i.e.wavelength  =  x /).
In practice the Fourier decomposition is performed by tabulating the surface density in the radial bins and applying a FFT to obtain  sg, at these bins.The force at particle locations is obtained by interpolation.
Furthermore, the vertical (self)-gravity force at a given radial location is assumed to be of the form (see Eq. ( A33) with (A34) and ( 75)) which is consistent with the vertical self-gravity (25) appearing in our hydrodynamic model.The local value ()/ () is calculated via a FFT, in similar manner as used for radial forces.Note that this approximation for the vertical self-gravity force is significantly more accurate than the assumption a constant frequency Ω  ≥ Ω 0 , as it accounts for local variations of the particle configuration which arise if the system is perturbed or unstable.Below we will find that our N-body simulations using ( 102) are largely consistent with the fully self-consistent self-gravitating simulations of MS23, as long as self-gravitational wakes are expected to be weak.

Results
Figure 8 collects various quantities extracted from simulations with  ℎ = 0.2, 0.3, 0.6, and 0.7, where  ℎ increases from top to bottom rows.Different optical depths  as defined by ( 96) are represented by different colors.From left to right columns, we plot the root-meansquared radial velocity   , the time-averaged vertical profile of the volume filling factor ff() as defined by (97), the time evolution of the midplane filling factor ff(0), and the vertical self-gravity factor Ω  /Ω 0 as defined by (102).Actually the occurrence of overstability has been determined from the evolution of radial Fourier modes of the velocity field.Simulations that develop viscous overstability are characterised by a saturated value   ≳ 1m Ω 0 .The plots in the second column confirm that the vertical density distribution indeed transitions from a near-Gaussian to a near-uniform distribution for increasing  (cf.§3.2), such that for  ≳ 1 the midplane region attains a fairly flat density profile.

Threshold for viscous overstability -constant vertical self-gravity
We first consider the approximation of a constant vertical self-gravity force, described by a constant Ω  , as given by ( 25), in absence of radial self-gravity.This assumption, while artificial, provides a simple means to shift smoothly between the non-selfgravitating and strongly self-gravitating limits, as far as vertical self-gravity is concerned.Table 3 summarises critical values of the optical depth   , the filling factor FF  , and the viscous slope   , for the onset of viscous overstability, resulting from the solution of (80).We adopt the values  102)).In the first and second columns runs which did (did not) develop overstability are drawn as solid (dashed) curves.
Ω  /Ω 0 = 1, 2 for the constant vertical self-gravity, as well as the constant coefficients of restitution  = 0.1, 0.3, 0.5.Results for both the slab and Gaussian vertical distribution are shown, and are compared to corresponding results reported in MS23 for the same values of  and Ω  .For Ω  /Ω 0 = 1, corresponding to vanishing vertical self-gravity, we find our values of FF  to agree quite well with those from N-body simulations.There is only little difference between the two density distributions regarding the values of FF  , which also holds for larger values of Ω  .On the other hand, the critical values for the optical depth   found in N-body simulations are much more closely reproduced with the Gaussian distribution.However, deviations from the N-body results increase with increasing value of Ω  .This is not surprising, and is most likely related to the circumstance that the vertical particle distribution in N-body simulations using a constant Ω  /Ω 0 ≳ 2 -in contrast to the vertical distribution assumed in our model -is subject to a layering, such that the number of vertical particle layers increases with increasing optical depth in a discontinuous fashion.This results in a oscillatory behaviour of the central filling factor with increasing optical depth, as the formation of each new layer abruptly changes the midplane filling factor.

Threshold for viscous overstability -vertical and radial self-gravity
We now compare critical quantities for the onset of overstability in presence of (radially varying) vertical and radial self-gravity.Here we consider a fixed  = 0.1 in simulations and model calculations.
In Figure 9 we first plot FF cr and  cr for increasing  ℎ (i.e.increasing strength of self-gravity), resulting from our hydrodynamic model and our N-body simulations adopting only vertical self-gravity forces.Again, we find overall only mild differences between FF cr for the slab (red) and Gaussian (blue) distributions, which match reasonably well the values found from our simulations (black).On the other hand, the values of  cr resulting with the Gaussian distribution significantly better match those from N-body simulations.The open circles are results obtained using the vertically static approximation "hstat" (see §5.3), where vertical motions are ignored.These results clearly demonstrate that it is the vertical self-gravity acting on overstable waves (via Eq. 74) which lowers the critical filling factor for overstability.
If we in addition include radial self-gravity (without otherwise changing any of the model parameters), a discrepancy develops with increasing  ℎ ≳ 0.4, such that the values of FF cr resulting from the model become increasingly smaller than those resulting from the simulations.As mentioned in §6.2.2 for  ℎ ≳ 0.5 we find overstability for all FF ≳ 0 when radial self-gravity is included.Thus, the effect of radial self-gravity seen in our model (cf.§6.2) is much stronger than in our N-body simulations, where it results only in a mild reduction of FF  .Nevertheless, we find that agreement between model and simulations in presence of radial self-gravity can be attained by using values of the bulk viscosity  0 larger than (30).This is discussed in Appendix D. There we also show that the measured values of the bulk viscosity in Salo et al. (2001) are likely to be altered if radial self-gravity is included in the simulations.
Furthermore, we find that critical values for FF and  resulting from our N-body simulations including radial and vertical self-gravity are very similar to those resulting from the fully self-gravitating simulations of MS23 as long as  ℎ ≲ 0.65 (Their Figures 10 and C1).For larger values of  ℎ overstability is being suppressed by strong self-gravity wake structures, not captured by our axisymmetric selfgravity implementaton (MS23).

Linear growth rates
In the previous section we compared the threshold values of important equilibrium quantities that lead to the onset of viscous overstability in our model and simulations.Here we compare the linear growth rates of resulting overstable modes, as predicted by our model and as measured from our simulations.Figure 10 shows growth rates measured in simulations (represented by circles) without self-gravity (top panel), with only vertical self-gravity (middle panel), and with both radial and vertical self-gravity (bottom panel).All simulations adopt a constant  = 0.5.The self-gravitating cases adopt  ℎ = 0.5.Compared in each panel are different optical depths.The curves represent hydrodynamic growth rates, based on the equilibrium state shown in Figure 1 (left panel) using a Gaussian distribution.
The method to measure growth rates from our simulations is largely similar to that used in Lehmann et al. (2017).That is, we seeded sinusoidal radial velocity perturbations of all modes (with identical amplitudes   Ω 0 ) simultaneously and tracked the amplitudes of individual modes in time by applying a Fourier transform on the total radially binned velocity profile.
In the case with  ℎ = 0 we find that the hydrodynamic growth rates are substantially larger than those measured from simulations.In contrast, the wavelengths of maximal growth are substantially smaller.Similar discrepancies were observed in previous hydrodynamic studies (e.g.Salo & Schmidt 2010 (their Figure 5) and Lehmann et al. 2017 (their Figure 4)).On the other hand, the threshold for overstability as found in the simulations is quite well reproduced by our model (see Figure 9), which is reflected by the reasonable agreement for the case with  = 3.The reason for the large hydrodynamic growth rates for the cases with  ≥ 4 is that the model is very close to the densest packing with FF ∼ 0.61.The value of FF in the simulations is slightly smaller, which might partially explain the large difference in the growth rates.Appreciably better agreement is seen in the growth rates resulting with the inclusion of vertical self-gravity.One possible reason for the better agreement is that the vertical layering of particles (not captured by the hydrodynamic model) is reduced when vertical self-gravity is included (Latter & Ogilvie (2008)).However, when including radial self-gravity, deviations are again larger, which is not surprising based on our findings described in the previous section.Note that for the hydrodynamic calculations in the bottom panel we used the bulk viscosity as displayed in Figure D1.Despite the threshold values for the onset of overstability being in good agreement (Figure D1, right panels), the linear growth rates show a similar deviation from the measured values as for the non-gravitating case.

summary
We developed a fluid model to describe the dynamics of dense planetary rings, such as the Saturnian A and B rings, or the Uranian  ring, inspired by the granular flow model of Hwang & Hutter (1995).In contrast to existing hydrodynamic models, and more in line with kinetic models, our model is largely based on parameters that can be determined from observations and experiments, such as the particle size , the coefficient of restitution , and the optical depth .Our model generalises previous hydrodynamic models of planetary rings by incorporating the dynamical evolution of the disc thickness, where motions are restricted to be symmetric with respect to the midplane.This generalisation enables us to study the effect of vertical motions, different vertical density stratifications, and the effect of variations in the volume filling factor on ring dynamics.Furthermore, while our model is able to adequately capture the behavior of densely packed particulate systems, it exhibits a substantially decreased mathematical complexity as compared to dense ring kinetic models.
We applied our model to compute axisymmetric planetary ring equilibrium states, where we considered a vertically uniform distribution, as well as a Gaussian vertical distribution for the volume mass density.The computed equilibrium states compare (at least qualitatively) well with those resulting from N-body simulations for small and large optical depths.
Next, we performed an axisymmetric linear stability analysis of our model, where we reproduced the viscous instability, the thermal relaxation mode, and the pair of overstable density waves, analysed in previous studies.In addition, we find a pair of vertical breathing modes, which describe oscillations of the ring thickness.We showed that in the appropriate limits the behaviour of all these modes agree with previous studies.
We then focused on the viscous overstability.In absence of selfgravity we find the occurrence of the viscous overstability as soon as a certain critical volume filling factor of the ring material is surpassed.
This finding is in agreement with the kinetic study of LO08 and the N-body simulations of MS23.In our model we find that the largest damping effect on overstable modes results from variations in the disc thickness.Furthermore, in agreement with the hydrodynamic studies of Spahn et al. (2000) and Schmidt et al. (2001), our model predicts that temperature variations have a mitigating effect on the viscous overstability.
In addition, we performed local N-body simulations including vertical and radial self-gravity.In absence of radial self-gravity the critical values for the filling factor and the optical depth for the onset of overstability agree well between our simulations and our hydrodynamic model using a Gaussian density distribution.In particular, the critical filling factor mildly reduces with increasing particle bulk density.As shown, this reduction is caused by the vertical self-gravity acting on overstable modes, and is absent if vertical motions are neglected.On the other hand, when radial self-gravity is included, the hydrodynamic model predicts for increasing  ℎ ≳ 0.4 increasingly smaller values of the critical filling factor and optical depth than our simulations.We speculate that the strong effect of radial self-gravity in our hydrodynamic model results from a reduction of the fluid's compressibility, which in turn mitigates the damping effect of the bulk viscous stress on overstable modes.We find that agreement between model and simulations can be attained by using increased values of the bulk viscosity, depending on the value of  ℎ .We also showed that the method used to measure the bulk viscosity in Salo et al. (2001) needs to be modified when radial self-gravity is included, thus providing a possible explanation for the need of an increased bulk viscosity.Moreover, we find that our N-body simulations with vertical and radial self-gravity correctly describe the threshold for viscous overstability as found in fully self-gravitating simulations of MS23 for  ℎ ≲ 0.65, beyond which self-gravity wakes suppress the instability in the latter simulations.

Caveats and outlook
Even though our hydrodynamic model is able to reproduce well (in some parameter regimes even quantitatively) the equilibrium state of a dense planetary ring as resulting from N-body simulations, agreement is substantially reduced when perturbations to the equilibrium state are considered.This can be seen for instance when comparing linear growth rates of overstable modes.
Furthermore, despite being largely framed in kinetic parameters, such as the optical depth  and coefficient of restitution , our model still needs to rely on a number of additional parameters, such as the "epicycle-parameter" , and probably most importantly, the bulk viscosity .These quantities are strongly related to the failure of a fluid description in the dilute regime, characterised by low particle collision frequencies.Perhaps the only way of getting around the need of such parameters is the employment of a kinetic model.However, even sophisticated dense ring kinetic models (Latter & Ogilvie 2008) struggle to achieve good agreement with N-body simulations when perturbations to the equilibrium are considered.
In addition, limitations of our description of disc thickness variations likely occur when vertical motions become strong, as for instance in the vicinity of a Lindblad resonance, or in strongly nonlinear density wave trains, which should lead to a "splashing" of the ring material (Borderies et al. 1985;Salo et al. 2001) .A correct analytical description in such situations might require a fully three-dimensional analysis of the fluid equations.
Regarding the description of the viscous overstability, our hydrodynamic model may potentially be improved by the construction of better suited transport coefficients, in particular the viscosity.On the other hand, kinetic models for the viscous overstability might be further improved by adopting a more accurate description of the vertical ring structure (Araki 1991), rather than adopting a pre-specified form.Even though this would further increase the mathematical complexity of such models, it could help to pin down possible sources for disagreement between models and N-body simulations.In this regard, it may also be useful to compare hydrodynamic and kinetic models directly.
The strong impact of radial self-gravity exhibited in our linear hydrodynamic model at large particle bulk densities is puzzling, and should be addressed in future work.This may require a modification to the computation of the bulk viscosity in N-body simulations as utilised in Salo et al. (2001).Also regarding the effect of radial selfgravity a direct comparison with dense ring kinetic models could be useful.
In a future study we plan to perform nonlinear hydrodynamic simulations based on our fluid model.These simulations can be used to investigate the effect of vertical thickness variations on the longterm nonlinear saturation of the viscous overstability.Moreover, once the axisymmetric dynamics are understood, 2D non-axisymmetric shearing sheet simulations can be performed to study the interaction of the viscous overstability with self-gravity wakes, which cannot directly be accounted for in a linear model.where the terms in the second bracket arise from radial and vertical self-gravity, respectively.It is a priori not clear how the above modifications affect the values of the bulk viscosity measured in self-gravitating N-body simulations.Nevertheless, we find that we can indeed adjust the value of  0 in our model for each value of  ℎ , such that the values of FF  nearly match those from our N-body simulations in the presence of radial self-gravity.These values of  0 are shown in Figure D1 (upermost panel), and are plotted against the corresponding optical depth values for which overstability occurs in the model (using a Gaussian density distribution) for the cases of radial and vertical self-gravity, and only vertical self-gravity.The corresponding values of FF  and also   are indeed nearly identical to those from our N-body simulations, as displayed in the middle and right panels, respectively.In absence of radial self-gravity we find optimal values 0.05 ≲  0 ≲ 4. The upper values are very similar to those of relation (30) which we used throughout the paper.Interestingly, for large optical depths our linear model requires only a very small bulk viscosity to obtain the best possible agreement with the simulation results.That is, at large optical depths the effect of the bulk viscosity becomes very small, such that the use of Eq. ( 30) yields very similar results for  ≳ 1. Substantially larger values of  0 are required to attain agreement when radial self-gravity is included.This suggests that the effect of radial self-gravity in (D8) dominates over the other additional terms, but this remains to be verified in future N-body simulations.
This paper has been typeset from a T E X/L A T E X file prepared by the author.Figure 9).Open (solid) circles correspond to computations without (with) radial self-gravity.

Figure 3 .
Figure 3. equilibrium values of pressure , viscosity , and cooling rate Γ, as well as their derivatives with respect to scaleheight  and surface mass density  for increasing filling factor ff.The plots correspond to the equilibrium shown in Figure 2 with  = 0.1.Note that the -derivatives generally yield negative values.

Figure 4 .
Figure 4.The "thermal" coefficients  2 and  3 , as given by (91) and (92), represented by diamond symbols.Also plotted are the coefficients  4 and  5 related to vertical motions in absence of self-gravity, as described in the text.The solid (dashed) curves apply to a Gaussian (uniform) distribution (18).The corresponding profiles of  2 and  3 for a uniform distribution (20) are identical to the Gaussian case.

Figure 5 .
Figure5.Maximum linear growth rates and corresponding wavelengths of the viscous overstability in absence of self-gravity, for increasing filling factor (top panels).The remaining panels show the pseudo-energy decomposition (Appendix B) of the modes for the various limiting cases as explained in the text.The pseudo-energies are normalised such that the dominating component (in all presented cases   ) yields unity.Furthermore, we plot the cube root of the total energy for increased visibility.

Figure 6 .
Figure 6.Phase shifts between the different hydrodynamic quantities of fastest growing overstable modes for a left traveling wave, resulting from numerical solution of the full linear eigenvalue problem (80) for non-gravitating (upper panel) and vertically self-gravitating rings (lower panel) with  ℎ = 0.5 (  = 204kgm −3 ).

Figure 7 .
Figure 7. Contours of maximum (across all wavelengths) linear growth rates of the viscous overstability in presence of vertical self-gravity for different equilibrium optical depths  and values of  ℎ , quantifying the strength of self-gravity.In the upper panels dashed curves describe a constant filling factor FF. The lower panel shows the pseudo-energy decomposition for fixed  ℎ = 0.5, as indicated by the black dash-dotted line in the upper panel.

Figure 8 .
Figure 8. Summary of a number of conducted N-body simulations wit  = 0.1 including radial and vertical self-gravity.Different rows compare different values of  ℎ as indicated.Different colors correspond to different optical depths  (96).From left to right columns we plot the time evolution of the radial component of the root-mean-squared radial velocity   , the time-average of the vertical distribution of the filling factor FF(z) (97), the time evolution of the central (midplane) filling factor FF(0), and the vertical frequency Ω  (cf.(102)).In the first and second columns runs which did (did not) develop overstability are drawn as solid (dashed) curves.

Figure 9 .
Figure 9. Critical values of the filling factor (top) and optical depth (bottom) for the onset of viscous overstability resulting from hydrodynamic model calculations using a uniform distribution (red solid curves) and a Gaussian distribution (blue solid curves) are compared with corresponding values from N-body simulations (black dotted curves) for increasing  ℎ in the presence of vertical self-gravity.The simulations used  = 0.1 and the model calculations assumed the equilibrium shown in Figure 2 (left panels).Open circles show results using the vertically static approximation (see §5.3).

Figure 10 .
Figure 10.Comparison of linear growth rates of the viscous overstability as measured from N-body simulations (cirles) and as predicted by our fluid model (curves).From top to bottom panels we consider a non-gravitating ring, a ring subject to vertical self-gravity, and a ring subject to radial and vertical self-gravity.The results assume  = 0.5.
Strictly speaking, in order to be suitable for our model, Equation (D1) requires correction in two aspects.First of all, we need to replace     →     +     =     + effective inclusion of vertical motions in this study, in contrast to the hydrodynamic model considered inSalo et al. (2001).Furthermore, in the presence of self-gravity, further modifications occur.This can be seen as follows.The portion of the radial momentum equation relevant to this discussion reads    = • •where  denotes the unity tensor.According toLynden-Bell & Kalnajs (1972) the second term in the bracket amounts to an isotropic tension, i.e. a negative pressure.If we now define the "generalised) =  0 −  0     +

Figure D1 .
Figure D1.Values of the bulk viscosity (upper panel) that result in good agreement between our hydrodynamic model and our N-body simulations for the critical values of FF (middle panel) and  (bottom panel).The bulk viscosity is drawn against the critical values of  for which overstability occurs in the model.Note that the latter decrease with increasing  ℎ (cf.Figure9).Open (solid) circles correspond to computations without (with) radial self-gravity.

Table 1 .
Important symbols and definitions.

Table 2 .
Critical values of important model parameters for the onset of viscous overstability in linear calculations excluding self-gravity, using a constant  , and for different limiting cases, as explained in the text.

Table 3 .
Critical values of important model parameters for the onset of viscous overstability using a constant  and Ω  .The three values of  cr ,  cr and ff cr correspond to model calculations using a slab profile (S) and a Gaussian profile (G), as well as N-body simulations (NB; values from Mondino-Llermanos & Salo (2023)).