# Stationary dynamics of planetary ring edge modes at Saturn's B ring edge

Pierre-Yves Longaretti

► To cite this version:

Pierre-Yves Longaretti. Stationary dynamics of planetary ring edge modes at Saturn's B ring edge. Icarus, 2024, 412, pp.115978. 10.1016/j.icarus.2024.115978. hal-04778403

HAL Id: hal-04778403

https://hal.science/hal-04778403v1

Submitted on 12 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Stationary dynamics of planetary ring edge modes at
Saturn’s B ring edge

Pierre-Yves Longaretti<sup>a,b</sup>

<sup>a</sup>Université Grenoble Alpes (UGA)

<sup>b</sup>CNRS/INSU — Institut de Planétologie et d'Astrophysique de Grenoble (IPAG) UMR 5274

**Abstract**

Recent analyses of the structure of the B ring edge have revealed the existence of a rather large number of mode-like features<fnref n="1" /> of azimuthal wavenumber ranging from *m* = 1 to *m* = 5. Several such mode-like features are associated with each *m*. The presence of so many features is a puzzle in itself.

The present work investigates whether these features can all be associated
with independent edge modes and develops to this effect a formal extension
of the still embryonic theory of edge modes as trapped waves. Based on this
formal extension, forced and free edge modes are explored in some detail. The
analysis explicitly shows how edge mode properties such as their resonance
location (for free modes) and edge amplitude (for the *m* = 2 forced mode) are
related to the ring surface density. Furthermore, the apsidal misalignment of
the forced *m* = 2 mode has been known for some time to be related to the
ring viscous dissipation in the edge region.

This investigation leads to a number of conclusions:

i/ Only one mode-like structure associated with a given *m* can be associated
with an edge mode, except possibly for *m* = 1 modes where each of the
three detected features can be interpreted as an independent edge mode.

ii/ The ring surface density at the edge is in the 200 g/cm² range under the assumption of dynamically independent edge modes, but a preliminary

Email address: pierre-yves.longaretti@univ-grenoble-alpes.fr (Pierre-Yves Longaretti)

<footnote n="1">¹In this work the expression mode-like features refers to sinusoidal in azimuth (at least in first approximation), variations of the edge shape.</footnote>
crude account of collective mode effects brings this down to a more realistic value of ~100 g/cm². This applies in the last 20 or 30 kilometers to the edge ($m \neq 1$). Further inside (one to two hundred kilometers), the surface density is somewhat more modest, ~60 g/cm² ($m = 1$).

iii/ The viscous dissipation at the edge is constrained by the observed misalignment between the mean longitude of Mimas and the forced edge mode. The implied level of dissipation at the edge is consistent with the now predominant idea that transport in the B ring, away from the edge, is dominated by self-gravity wakes.

iv/ The $m = 2$ mode forced by Mimas is never nonlinear enough by itself to produce the viscous angular flux momentum reversal required to make the edge effectively confined by the associated satellite torque. This explains the failure of Hahn et al. (2009) to explain the confinement of the edge by Mimas except for unrealistic stress tensor properties. It is argued that the presence of so many modes and mode-like features provides a way to restore in a somewhat different form the edge confinement process of Borderies et al. (1982), thereby resolving the conundrum raised by Hahn et al. (2009).

*Keywords:* Planets: ring dynamics, Dynamics: edge modes, Saturn: B ring

# 1. Introduction

The B ring edge exhibits a complex time varying shape that is apparently due to the presence of a number of sinusoidal (in azimuth and time) radial variations (Nicholson et al., 2014; French et al., 2023). The mode driven by Mimas is probably the oldest-known feature contributing to these patterns. The Mimas 2:1 resonance has long been associated with the B ring outer edge, and the dynamics of the confinement of the edge by this satellite has been elucidated in landmark papers (Borderies et al., 1982, 1984) although some puzzling questions remain (Hahn et al., 2009). However the Cassini data analysis reveals the existence of a host of other mode-like features at the B ring edge besides the one excited by Mimas.

These features are found by fitting a superposition of sinusoidal oscillations to the actual ring shape. Each of these is characterized by an amplitude, an azimuthal wavenumber $m$, a pattern speed $\Omega_p$ and a phase. The pattern speed and wavenumber are associated to a resonance location $a_r$ in the usual way, i.e., $m[\Omega(a_r) - \Omega_p] = \kappa(a_r)$. The search for such features is limited to
$a_r$ “close” to the B ring outer edge radius, i.e., within a few tens of kilometers inside the B ring outer edge boundary for $m \neq 1$ and a few hundreds kilometers for $m = 1$. In practice, a systematic scan of the expected range of pattern speeds is performed for each $m$ in order to identify possible mode-like features. A positive identification is made when a substantial reduction is obtained in the residual of the least-square fit of the model to the actual edge shape. Quite surprisingly, two mode-like features associated to $m = 2$ have been found, one forced by Mimas, as expected, and a second one with close pattern speed (in the above meaning) and amplitude (Spitale and Porco, 2010; Nicholson et al., 2014). Multiple mode-like features have also been found for $m = 1, 3, 4$ and 5. All the features associated with a given $m$ can also be viewed as slow variations of a single, fundamental mode-like feature, and this point of view is adopted in Nicholson et al. (2014) and French et al. (2023), but this is essentially conventional from a kinematic point of view. The reader is referred to these two publications for more detail on the fitting procedure.

The first question raised by these findings is whether these mode-like features are the signature of true edge modes or not. Edge modes are characterized by their number of radial nodes besides their azimuthal wavenumber, pattern speed, amplitude and phase, and can be interpreted as standing density waves trapped between the resonance radius and the ring edge. This picture has long been put forward by Peter Goldreich for all kind of standing modes, including eccentric narrow rings. The basic dynamical description of these modes has been outlined in Longaretti (2018), following an unpublished lead by Peter Goldreich.

The first purpose of the present work is to expand upon this preliminary analysis of edge modes, in an ideal setting where mode interactions are neglected, with specific focus on the B ring outer edge. Such an independent mode setting is meaningful at least in the linear limit. However, all types of prominent modal motions in major planetary ring systems are nonlinear to some extent; this is the case, e.g., of the strongest density waves in Saturn's rings, or of the largest narrow rings in the Saturnian and Uranian systems. Edge modes do not differ much in this respect, neither in the level of nonlinearity observed, nor in the fact that they share with these other features the same basic dynamics. Consequently, one may expect that nonlinear interactions between modes may change to some extent the properties of edge modes, and the ring diagnostics derived from them.

Still, as no extensive study of such edge modes has appeared in print yet
(besides the preliminary analysis outlined in Longaretti 2018), it is of some
interest to present a more detailed analysis of such modes, in particular of the
ones actually detected at ring edges — the B ring edge in the present case.
Furthermore, it is argued in the concluding section that mode interactions
will not significantly alter the results obtained in the independent mode limit.
In any case, this constitutes a necessary first step in itself, and a useful and
informative comparison point for more refined investigations.

The dynamical analysis presented in this paper is used most notably to
extract information on ring properties, in particular the ring surface density
close to the edge, and the ring effective viscosity in the edge region. This
is achieved first through a semi-quantitative analysis of both free and forced
modes which also clarifies the physics involved, and next through detailed
numerical solutions of the mode structure equations. As a side benefit, this
investigation will pinpoint the main reason why Hahn et al. (2009) were not
successful in their detailed investigation of the B ring edge confinement by
Mimas.

This work relies on the streamline formalism, as restated in terms of
epicyclic variables (Borderies and Longaretti, 1987; Longaretti and Borderies,
1991; Borderies-Rappaport and Longaretti, 1994) from its original elliptic
variable formulation (Goldreich and Tremaine, 1979a; Borderies et al., 1983,
1985, 1986), referred to collectively as the BGT papers. This epicyclic vari-
ables reformulation leads to a closer connection with observations while re-
moving a number of internal inconsistencies in the formalism, without affect-
ing its main formal characteristics and physical conclusions (see Longaretti
1992 and Longaretti 2018 for a review).

This paper is organized as follows. Section 2 recalls the formalism kine-
matics in epicyclic variables; this also constitutes an opportunity to formulate
the kinematic foundation of the trapped wave picture of edge modes in a more
precise and extensive way than in Longaretti (2018). Section 3 introduces
the relevant stationary dynamics for single modes, and the two levels of ap-
proximations used in subsequent numerical solutions of the mode structure
equations. The first level solves for the eccentricity profile of the mode while
neglecting the effect of the stress tensor; the second level solves next for the
phase profile of the mode, making use of the eccentricity profile just estab-
lished. The self-consistency of the procedure is ensured by the smallness of
the ring stress tensor compared to its self-gravity.

The remainder of the paper mostly focuses on the B ring edge context.
Section 4 develops a simplified, semi-quantitative analysis of the stationary
dynamics of a single (free or forced) mode based on a two-streamline approximation adapted to edge modes; this is explicitly applied wherever possible to the B ring edge modes. The last part of the paper is devoted to more precise numerical analyses of single mode structure and properties. Section 5 discusses single free modes, while section 6 deals with the $m = 2$ forced mode excited by Mimas. These two sections are completed by a discussion of boundary conditions (Appendix C), a point which has never been broached in the streamline formalism, due to the limited number of numerical analyses of the dynamics in this framework (with the notable exception of Borderies et al. 1986, 1989). Section 7 collects and discusses the main findings and their implications.

## 2. Kinematics

### 2.1. Epicyclic variable streamline kinematics

The streamline formalism used in this work formulates ring dynamics as a fluid problem in Lagrangian form. The original publications (Goldreich and Tremaine, 1979a; Borderies et al., 1982, 1983, 1985, 1986) — especially the earlier ones — were not quite clear on whether they dealt with a particle or fluid problem. This was first rigorously discussed in Longaretti (1992). The following presentation of the kinematics closely follows Longaretti (2018).

On the one hand, ring fluid particles must follow eccentric orbits, as the Navier-Stokes equations of fluid motion (or the moment equations of more fundamental kinetic approaches) in Lagrangian variables with all forces neglected but the planet’s (zeroth order approximation) reduces to a test particle equation of motion ($dr/dt = -\nabla\Phi_p$ where $\Phi_p$ is the potential of the planet). To lowest order in deviation from circularity, the resulting positions of test fluid particles read:

$$r = a_e[1 - \epsilon_e \cos M_e], \qquad (1)$$

$$\theta = \varpi_e + M_e + 2\epsilon_e \sin M_e = \varphi_e + 2\epsilon_e \sin M_e, \qquad (2)$$

where $a_e$, $\epsilon_e$, $\varpi_e$ and $M_e$ are the epicyclic semi-major axis, epicyclic eccentricity, epicyclic periapse angle and epicyclic mean anomaly; in Eq. (2), terms of relative order $J_2\epsilon_e$ have been neglected (the ratio $\Omega/\kappa$ has been set to 1 in the first order terms). The mean epicyclic longitude is $\varphi_e = \varpi_e + M_e$, and $\Omega$ and $\kappa$ are the epicyclic angular velocity and radial frequency at $a_e$. For unperturbed motions, $dM_e/dt = \kappa$ and $d\varphi_e/dt = \Omega = \dot{\varpi}_e + \kappa$.
The related radial and angular velocities read (to the same level of precision in $J_2\epsilon_e$)

$$ \frac{dr}{dt} = a_e \epsilon_e \Omega \sin M_e, \qquad (3) $$

$$ \frac{d\theta}{dt} = \Omega + 2\Omega\epsilon_e \cos M_e. \qquad (4) $$

In terms of the leading planetary zonal gravity harmonic coefficient $J_2$, the epicyclic angular velocity and epicyclic frequency read

$$ \Omega(a_e) = n(a_e) \times \left[ 1 + \frac{3}{4} \left( \frac{R_p}{a_e} \right)^2 J_2 \right], \qquad (5) $$

$$ \kappa(a_e) = n(a_e) \times \left[ 1 - \frac{3}{4} \left( \frac{R_p}{a_e} \right)^2 J_2 \right], \qquad (6) $$

where $M_p$ and $R_p$ are the planet mass and radius, and $n$ is the elliptic mean motion:

$$ n(a_e) = \left( \frac{GM_p}{a_e^3} \right)^{1/2} . \qquad (7) $$

More complete expressions can be found in Borderies and Longaretti (1987) and Borderies-Rappaport and Longaretti (1994).

On the other hand, fluid particle positions are assumed to be specified by an $m$-lobe shape. This assumption is motivated by the form of the deviations from circular motion observed, e.g., for an edge mode or a narrow ring:

$$ r = a \left(1 - \epsilon(a, t) \cos[m(\varphi - \Omega_p t) + m\Delta(a, t)]\right), \qquad (8) $$

$$ \theta = \varphi + 2\epsilon(a, t) \sin[m(\varphi - \Omega_p t) + m\Delta(a, t)], \qquad (9) $$

where $\varphi = \varphi_0 + \Omega t$ is the fluid particle mean epicyclic longitude, $\epsilon$ its eccentricity, $m\Delta$ the apsidal shift and $\Omega_p$ the pattern speed associated with the mode, and self-consistently determined by its dynamics. The azimuthal wavenumber is assumed to be positive ($m > 0$) consistently with the fact that one is interested by edge modes at an outer ring edge, corresponding to an inner Lindblad resonance.

The velocity field associated with Eqs. (8) and (9) is given by:

$$ u_r = a\epsilon\Omega\sin[m(\varphi - \Omega_p t) + m\Delta(a, t)], \qquad (10) $$

$$ u_{\theta} = \Omega a (1 + \epsilon \cos[m(\varphi - \Omega_p t) + m\Delta(a, t)]). \qquad (11) $$
To leading order in eccentricity, Eqs. (8) and (9) lead to

$$r = a(1 - \epsilon(a, t) \cos[m(\theta - \Omega_p t) + m\Delta(a, t)]), \quad (12)$$

which describes the shape of the quasi-stationary path followed by fluid particles in the frame rotating at the pattern speed, i.e., fluid streamlines, to a high degree of precision.

Eqs. (1), (2) and (8), (9) and the related velocity relations which follow from Eqs. (3), (4) and (10), (11) can be satisfied simultaneously only if $a = a_e$, $\epsilon = \epsilon_e$, $\varphi_e = \varphi$ and

$$M_e = m(\varphi - \Omega_p t) + m\Delta, \quad (13)$$

i.e., if the following relations hold

$$\frac{dm\Delta}{dt} = -m(\Omega - \Omega_p) + \Omega - \dot{\varpi}_e, \quad (14)$$

$$\varpi_0 = \varphi_0(1 - m) - m\Delta_0, \quad (15)$$

where $\varpi_0$ and $\varphi_0$ are the periapse angle and mean epicyclic longitude of the fluid particle at $t = 0$, and where $m\Delta_0$ is the streamline apsidal shift at the same time. The second relation expresses the constraint that all particles with identical $a$ must satisfy in order to belong to a common $m$-lobe streamline in the rotating frame, while the first is a necessary condition for the pattern to be maintained at all times. Thus both equations represent kinematic constraints for all particles to undergo the type of collective motions of interest. In this relation, $\dot{\varpi}_e$ includes the all perturbations, and most notably collective effects (self-gravity and stress tensor).

The first relation is used to obtain the perturbation equation for $m\Delta$ from the perturbation equation of the periapse angle. In these relations, the contribution of $d\varphi_e/dt$ has been neglected, as it is smaller by a factor $\epsilon$ than the leading contributions. In rings, the contribution of the self-gravity is usually the largest one and maintains this relation everywhere in the vicinity of the resonance location, in particular in stationary settings ($dm\Delta/dt = 0$).

A resonance radius $a_r$ can be associated with the pattern speed, and is implicitly given by

$$m[\Omega(a_r) - \Omega_p] = \kappa(a_r). \quad (16)$$

By definition, this radius is the resonance location of free edge modes. For forced mode, this resonance location is the usual resonance with the satellite forcing potential.
Finally, following Borderies et al. (1983) and Shu et al. (1985), let us define a complex eccentricity:

$$Z \equiv \epsilon \exp im\Delta. \qquad (17)$$

The nonlinear parameters $(q, \gamma)$ are related to $Z$ through

$$qe^{i\gamma} = a \frac{dZ}{da} e^{-im\Delta}. \qquad (18)$$

In the streamline formalism, these parameters allow us to relate the perturbed $(\sigma)$ and unperturbed $(\sigma_0)$ surface densities

$$\sigma(a, \varphi, t) = \frac{\sigma_0(a, t)}{J} = \frac{\sigma_0(a, t)}{1 - q \cos[m(\varphi - \Omega_p t) + m\Delta + \gamma]} \qquad (19)$$

This relation follows from the conservation of mass $\delta M$ between the unperturbed (at $a, \varphi$ within $\delta a, \delta \varphi$) and perturbed (at $r, \theta$ within $\delta r, \delta \theta$) elementary surface elements, to leading order in eccentricity: $\delta M = \sigma_0 a \delta a \delta \varphi = \sigma r \delta r \delta \theta = J \sigma a \delta a \delta \varphi$ ($J > 0$ as $q < 1$ due to collective effects).

## 2.2. Trapped wave picture of edge modes

When the trapped wave picture of edge modes is invoked, edge mode eccentricities are related to underlying amplitudes and phases of the long trailing and leading density waves reflected into one another at the edge and at the resonance implicitly defined by their common pattern speed. The superposition of these two waves creates a standing wave, which appears as an edge mode. The radial extent between the resonance location and the ring edge constitutes the cavity of this trapped standing wave.

In this picture, one needs to represent the radial position of fluid particles as a superposition of an outgoing long trailing wave and an ingoing long leading wave. The common (except for the sign) radial wavenumber of these ingoing and outgoing waves is specified by the dispersion relation of density waves; they also have the same amplitude because their propagation equations are identical, and because they reflect without loss of angular momentum flux into one another:

$$r = a \left[ 1 - \varepsilon \cos[m(\varphi - \Omega_p t) + m\Delta_{+}^{dw}] - \varepsilon \cos[m(\varphi - \Omega_p t) + m\Delta_{-}^{dw}] \right],$$

where $\varepsilon$ is the (common) wave eccentricity and $m\Delta_{\pm}^{dw}$ are the trailing (+) and leading (-) wave apsidal shifts (the superscript $dw$ serves as a reminder
that we are dealing here with density waves producing an edge mode, and not with the edge mode itself). These two phases are related to their common radial wavenumber ±k through

$$m\Delta_{\pm}^{dw} = \pm \int_{0}^{x_e} kdx + \phi_{\pm} \equiv \pm m\Delta_{0}^{dw} + \phi_{\pm}, \quad (20)$$

where $x = (a - a_r)/a_r$ is the relative distance to the resonance $a_r$, $x_e = (a_e - a_r)/a_r$, and $k$ is the trailing wave radial wave number, $\propto x$ (see, e.g., Goldreich and Tremaine 1979b, or Shu 1984). $\phi_{\pm}$ is the constant of integration (value of the apsidal shift) of the outgoing (+) and ingoing (-) wave and is needed to account for the effect of dissipation on edge mode structure. Indeed, the outgoing and ingoing waves are both damped, breaking the symmetry of their phases and amplitudes that is implied by their dynamical equations in the absence of dissipation; this term is absent from the preliminary analysis of Longaretti (2018).

Defining $m\Delta = (m\Delta_+ + m\Delta_-)/2$ and $\delta m\Delta = (m\Delta_+ - m\Delta_-)/2$, one has:

$$r = a \left[ 1 - \varepsilon \cos(m\Delta_0^{dw} + \delta m\Delta) \cos[m(\varphi - \Omega_p t) + m\Delta] \right]. \quad (21)$$

In this expression, $m\Delta$ is the arbitrary phase of the resulting standing wave pattern, except for the forced mode where this phase closely follows the phase of its forcing potential (see sections 4 and 6.2 on this point). In forced edge modes, and in the absence of dissipation due to ring particle collisions, this mode apsidal shift is independent of location and can be set to zero, an assumption made at the onset in Longaretti (2018) but which must be lifted for the present purpose.

Comparing Eqs. (8) and (20) leads to $\varepsilon = \varepsilon \cos(m\Delta_0 + \delta m\Delta)$. For edge modes forced by an external satellite, the linear asymptotic theory of forced density waves gives $\delta m\Delta = -\pi/4$. There is no such simple prescription for free waves, and $\delta m\Delta$ will be determined later on from numerical simulations.

Note that the identification $\varepsilon = \varepsilon \cos(m\Delta_0 + \delta m\Delta)$ implies that $\varepsilon$ is not necessarily positive. This convention makes both the amplitude $a\varepsilon$ and the phase $m\Delta$ continuous and differentiable with $a$. The sign of the amplitude can be arbitrarily chosen for a free mode, as a change of its sign amounts to a change of $m\Delta$ by $\pi$, and as the phase of free modes at $t = 0$ is arbitrary. This is not true of forced modes, where the phase is imposed by the satellite forcing potential, but even in this case, a simultaneous change of sign of the amplitude and of the phase by $\pi$ leaves the mode unchanged.
This trapped wave picture description is not strictly necessary and the kinematic relations of section 2.1 can be assumed from the onset to describe edge modes. However, this gives some intuitive understanding of the physical nature of edge modes.

## 2.3. Mode-like features

In this work, the denomination “mode” or “edge mode” refers explicitly to trapped waves. Viscous overstabilities are expected to contribute to the growth of modes (Borderies et al., 1986; Papaloizou and Lin, 1988) but also to the growth of mode librations (Longaretti and Rappaport, 1995), which are intrinsically non-stationary oscillations of a mode eccentricity and apsidal shift. Although such librations can formally be kinematically described in a way similar to modes at the edge (Nicholson et al., 2014; French et al., 2023), they are not independent modes, but features belonging to actual modes. Such features are referred to as mode-like features in this work. Other potential sources of such features besides viscous overstabilities are non-resonant mode interactions mediated through the ring self-gravity and stress tensor.

# 3. Dynamics

Dynamical analyses based on the streamline formalism have assumed that only one mode is present at any given location, with the exception of Longaretti (1989). We make the same single mode assumption in the present work as a first approximation.

The B ring outer edge is maintained by the $m = 2$ resonance with Mimas, and the edge location may be assumed to have achieved equilibrium on the time-scales on interest here. One can therefore assume that the semimajor axes $a$ of all streamlines are constant on these time scales. The confinement process also constrains the shape of the surface density close to the edge. Here, instead, as the confinement process is not included in the dynamical analysis, the shape of the surface density profile is assumed; in practice, only constant density profiles have been explored (except briefly in section 6.2.3), motivated by the fact that the B ring edge is sharp.

Furthermore, in a first approximation, the effect of the ring stress tensor is neglected, as it is usually small compared to the ring self-gravity. In other words, self-gravity is the only collective effect in this limit, and the apsidal shift $m\Delta$ becomes constant throughout the edge mode and in time, as can
be self-consistently checked with the help of the dynamical equations. This constant apsidal shift can be chosen to be zero by a change of origin of angles. Consequently, in this approximation, $Z = \epsilon$, $\gamma = 0$ and $q = ad\epsilon/da$.

In a second approximation, a small stress tensor can be introduced, which produces a small shift in $m\Delta$ across the extent of the edge mode. The observed shift — $\approx m \times 1^{\circ}$ to $m \times 3^{\circ}$ for the forced $m = 2$ mode — provides a direct measure of the ratio of the stress tensor to the ring self-gravity time scales (and therefore of their relative magnitude), as discussed in section 4 [see Eq. (44)]. To leading order in this ratio, the equation specifying the eccentricity profile of the free and forced mode is unchanged; the only difference between the first and second approximation is the introduction of an equation for the apsidal shift profile $\Delta(a)$.

These two levels of approximation are formalized in this section. For completeness, the basis of the perturbation approach (epicyclic perturbation equations) is recalled first.

### 3.1. Perturbation equations

The streamline formalism borrows from celestial mechanics an osculating motion approach, transposed to epicyclic motion instead of the more familiar elliptic setting (Borderies and Longaretti, 1987; Longaretti and Borderies, 1991; Borderies-Rappaport and Longaretti, 1994; Longaretti, 1992). In other words, perturbed motions are represented by the same relations for the ring fluid particle position and velocity, except that the epicyclic elements must depend on time. The shift from elliptic to epicyclic theory is necessary to remove a (non-negligible) inconsistency of absolute order $J_2$ in the original formulation of the streamline formalism, and makes the associated osculating elements nearly identical to the observed ones (e.g., in the determination of ring orbits), as was erroneously implicitly assumed in the original formulation.

The relevant perturbation equations read

$$ \frac{d\epsilon}{dt} = \frac{1}{na_e} [R \sin M_e + 2S \cos M_e], \qquad (22) $$

$$ \frac{d\varpi_e}{dt} = \dot{\varpi}_p + \frac{1}{na_e\epsilon} [-R \cos M_e + 2S \sin M_e], \qquad (23) $$

to leading order in $\epsilon$ and $J_2$; $R$ and $S$ are the radial and azimuthal components of the perturbing acceleration, and $\dot{\varpi}_p = \Omega - \kappa$ represents the effect of the
planet oblateness on the precession of the apses. Note that $d\varphi_e/dt - \Omega = d(\bar{\omega}_e + M_e)/dt - \Omega$ is negligible compared to $d\bar{\omega}_e/dt$ (higher order in $\epsilon$) and is therefore neglected.

These equations are nearly identical to their elliptic counterparts, thanks to an appropriate definition of epicyclic elements. This in turn ascertains that the results obtained from the analysis of streamline dynamics in the original BGT publications is correct in spite of the inconsistencies pointed out above, once the osculating elements are interpreted as epicyclic elements. For more details see Longaretti (1992) and Longaretti (2018) or the papers mentioned above.

In the streamline formalism, these equations are used only to determine the secular evolution of the epicyclic elements. The short period terms are averaged over, and are anyway of negligible amplitude. The perturbation equation for $Z$ used below follows by combining these relations with Eqs. (14) and (17), and results from these short-period averages. For simplicity, the short-period averaged quantities (eccentricity, apsidal shift, etc) are denoted in the same way as the instantaneous ones; this is a common convention in the streamline formalism literature.

## 3.2. Stationary self-gravitational dynamics

In this limit, stationary edge modes can be described in the frame rotating at the pattern speed and the perturbation equations reduce to (see section 3 of Longaretti 2018)

$$ \frac{dZ}{dt} = \left(\frac{dZ}{dt}\right)_{sg} + \left(\frac{dZ}{dt}\right)_{pl} + \left(\frac{dZ}{dt}\right)_{sat} = 0, \qquad (24) $$

where

$$ \left(\frac{dZ}{dt}\right)_{sg} = -i \int da' \frac{2na^3\sigma_0'}{M_p} H(q_{aa'}^2) \frac{\epsilon - \epsilon'}{(a - a')^2}, \qquad (25) $$

$$ \left(\frac{dZ}{dt}\right)_{pl} = -i[m(\Omega - \Omega_p) - \kappa]\epsilon, \qquad (26) $$

$$ \left(\frac{dZ}{dt}\right)_{sat} = -in \frac{a\Psi_{m_s k_s}}{2GM_p} \delta_{m,2} \qquad (27) $$

which represent respectively the effect of the ring self-gravity, the planet-induced relative drift of the streamlines and the satellite resonant forcing (see ). In these relations, $m\Delta = 0$ has been used, so that $Z = \epsilon$.
In the last relation, $\delta_{m,2}$ is a Kronecker-like symbol, used as a reminder that this applies only to the forced $m = 2$ mode: $\delta_{m,2} = 1$ if $m = 2$ and $= 0$ otherwise. The other quantities are defined as follows: $M_p$ is the mass of the planet (Saturn), $q_{aa'} = a(\epsilon - \epsilon')/(a - a')$, $H(q^2) = [1 - (1 - q^2)^{1/2}]/[q^2(1 - q^2)^{1/2}]$, and $n = (GM_p/a^3)^{1/2}$ is the elliptic mean motion. $H(q^2)$ diverges as $q \to 1$; this prevents streamline crossing, which would occur at $q = 1$. Finally, the satellite forcing term $\Psi_{m_s k_s}$ is explicitly given in Appendix B.

The problem is further simplified for numerical purposes by recasting the self-gravity integral into a sum over $N$ streamlines, so that Eq. (24) reduces to$^2$:

$$ \sum_{j \neq i} \frac{\bar{n} m_j}{\pi M_p} a^2 H(q_{ij}^2) \frac{\epsilon_i - \epsilon_j}{(a_i - a_j)^2} - (\kappa_i - m\Omega_i)\epsilon_i = m\Omega_p\epsilon_i, \quad (28) $$

$$ \sum_{j \neq i} \frac{\bar{n} m_j}{\pi M_p} a^2 H(q_{ij}^2) \frac{\epsilon_i - \epsilon_j}{(a_i - a_j)^2} - [\kappa_i - m(\Omega_i - \Omega_p)]\epsilon_i = -\bar{n} \frac{a \Psi_{mk}}{2GM_p}. \quad (29) $$

The first equation applies to free edge modes, the second to the forced $m = 2$ one; $m_i$ is the mass of streamline $i$. By assumption, all $N$ streamlines are of equal width $\delta a$ and the surface density $\sigma_0$ is assumed constant in the edge region so that $m_i = 2\pi\sigma_0 a \delta a$ for all streamlines. For simplicity, the semimajor axis $a$ is assumed constant except in the difference $a_i - a_j$ (the error introduced by this simplification is negligible), and $\bar{n} = (GM_p/a^3)^{1/2}$ is the associated elliptic mean motion.

These equations can be recast as$^3$:

$$ \sum_{j=1}^{N} L_{ij} \epsilon_j = \frac{m \Omega_p}{\bar{n}} \epsilon_i, \quad (30) $$

$$ \sum_{j=1}^{N} \left( L_{ij} - \frac{m \Omega_p}{\bar{n}} I_{ij} \right) \epsilon_j = - \frac{a \Psi_{mk}}{2GM_p} \equiv F, \quad (31) $$

<footnote>$^2$These relations make use of an implicit origin of time and angles for the $m = 2$ forced mode, so that the forcing satellite potential has zero phase at $t = 0$ and $m\Delta = 0$. Such choices are always possible when dealing with a single mode. Using an arbitrary origin of angles requires to add a factor $\exp im\Delta_s$ to the right-hand side of Eq. (27); in the absence of dissipation, $\Delta_s = \Delta$ and this common phase disappears from the equation.</footnote>

<footnote>$^3$I am indebted to Peter Goldreich (private communication) for pointing out the usefulness of this reorganization.</footnote>
where $L_{ij}$ is implicitly defined by reorganizing the terms on the left-hand side of Eq (28) (after division by $\bar{n}$) and $I$ stands for the identity matrix.

In the linear limit $q \to 0$, $H(q^2) = 1/2$. Eq. (30) then reduces to an eigenvalue problem for the vector made up by collecting the eccentricities $\epsilon_i$ of all streamlines, and with eigenvalue $m\Omega_p/\bar{n}$. In the trapped wave picture of edge modes, this equation describes the structure of a standing density wave. One expects therefore that the eigenvalue problem will exhibit $N$ solutions differing by the number of radial nodes (from 0 to $N-1$). This expectation is borne out in the numerical solutions that are presented later on.

As the number of streamlines can be arbitrarily chosen and arbitrarily large, one may ask what limits in practice the number of radial nodes. Several such limits can be expected on physical grounds. First, modes of short wavelengths will be strongly damped, and eventually, not detectable at smaller and smaller wavelengths. Second, no wavelength smaller than the scale of the granularity of the rings (a few meters or tens of meters) can exist. Finally, on scales smaller than the scale height (itself of the order of the granularity scale), the physical approximations leading to the two-dimensional problem solved here break down.

### 3.3. Stress tensor correction to the dynamics

Introducing a stress-tensor correction to the dynamics modifies Eq. (24) into

$$ \frac{dZ}{dt} = \left(\frac{dZ}{dt}\right)_{sg} + \left(\frac{dZ}{dt}\right)_{pl} + \left(\frac{dZ}{dt}\right)_{sat} + \left(\frac{dZ}{dt}\right)_{vis} = 0, \quad (32) $$

where the various contributions now read

$$ \left(\frac{dZ}{dt}\right)_{sg} = -i \int da' \frac{2na^3\sigma_0'}{M_p} H(q_{aa'}^2) \frac{Z-Z'}{(a-a')^2}, \quad (33) $$

$$ \left(\frac{dZ}{dt}\right)_{vis} = \frac{i}{an\sigma_0} \frac{\partial}{\partial a} \left( (t_2 + it_1) \frac{a}{q} \frac{\partial Z}{\partial a} \right), \quad (34) $$

$$ \left(\frac{dZ}{dt}\right)_{pl} = -i[m(\Omega - \Omega_p) - \kappa]Z, \quad (35) $$

and where $t_1$ and $t_2$ are viscous-like and pressure-like (respectively) effective components of the pressure tensor. These will be specified below. The satellite term is unchanged. Eq. (34) is more compact but strictly equivalent to the expression given in Eq. (10.58) of Longaretti (2018).
Only $t_2$ enters the corrected Eqs. (28) and (29), which corresponds to the real part<fnref n="4" /> of $Z$. To leading order in stress-tensor correction, this contribution can be neglected, leaving these two relations unchanged. The imaginary part of $Z$ then obeys the following equation, in the stationary limit:

$$ \sum_{j=1}^{N} (L_{ij} - \frac{m\Omega_p}{\bar{n}} I_{ij}) \epsilon_j \sin m\Delta_j = \frac{2\pi}{\bar{n}^2 m_i} \Delta^{\pm} \left( \frac{a}{q} \frac{d\epsilon}{da} t_1 \right). \quad (36) $$

where $\Delta^{\pm}X$ stands for the variation of $X$ from the current streamline inner boundary to its outer boundary. The other quantities are the same as in subsection 3.2. This equation applies to both free and forced edge modes. In the absence of dissipation, $t_1 = 0$ and $m\Delta_j = 0$ (consistently with the choice $m\Delta = 0$ far in the mode evanescent region). Thus, in the absence of dissipation, the mode apsidal shifts are all strictly aligned. Conversely, a small difference of apsidal shifts between the inner mode region and the edge requires the presence of viscous dissipation, in the stationary dynamics limit. As this apsidal shift difference is small, $\cos m\Delta_j = 1$ to a very high level of approximation; this completes the self-consistency check of the approximation used in Eqs. (28) and (29).

To sum up, from a practical point of view, Eq. (36) is solved by first solving the stress-free dynamics of the previous subsection and using the resulting eccentricity profile in this equation to determine the apsidal shift profile. The error introduced by this two-step process is negligible because the stress tensor ($t_1$ and $t_2$) is small compared to the ring self-gravity.

### 3.4. Stress tensor model

Only a simple approximation of the stress-tensor is needed for the purposes of this work. Namely

$$ t_1(q) = T_1^d(q) \nu_{eff} \sigma_0 \bar{n}, \quad (37) $$

$$ T_1^d(q) = q \frac{q_1 - q}{(1 - q^2)^{3/2}}, \quad (38) $$

where $\nu_{eff} \simeq c^2/\bar{n}$ ($c$ is the local velocity dispersion<fnref n="5" />) and $q_1$ is a parameter. This model is motivated both by the simple hydrodynamic approximations

<footnote n="4">⁴This ignores the negligible apsidal shift difference across the mode due to the ring stress tensor. For simplicity, the apsidal shift at $a \to 0$ has been set to zero.</footnote>

<footnote n="5">⁵This expression is appropriate for dense enough rings, which is the case of the B ring.</footnote>
described in Appendix A and the empirical family of models of Borderies et al. (1986). The simple hydrodynamic models of Appendix A help to understand the physical nature of $t_1$ and $t_2$, which embody respectively the azimuthally averaged viscous damping and pressure in this hydrodynamic limit. The dimensional scaling of the above model is correct for all contexts where the ring optical depth $\tau \gtrsim 1$ as is the case at the B ring edge (with possible reformulations of the velocity dispersion in terms of other ring parameters); the dimensionless $q$ dependence has been somewhat simplified to encompass both viscously stable and overstable contexts, while satisfying the physical constraints identified in Borderies et al. (1986), and discussed now.

Stable rings have $t_1 < 0$ for all $q$ (i.e., $q_1 \le 0$) while viscously overstable ones have $t_1 > 0$ at small $q$ and $t_1 < 0$ at large $q$, coloredri.e., $q > q_1$ (Borderies et al., 1985, 1986; Longaretti, 2018). In this second context, $q_1$ varies from $\lesssim 0.5$ (Mosqueira, 1996) to 0.7 in all known models (including the models of Appendix A) but the relevant point for the present purpose is the sign of $t_1$, as $T_1$ is of order unity for the moderate values of $q$ appropriate for the B ring edge modes. When $t_1 > 0$ edge modes usually undergo libration motions, making them non-stationary, but this does not change the results that will be deduced from Eq. (36), as these apply to the equilibrium point around which librations may or may not take place (Longaretti and Rappaport, 1995). When rings are stable ($t_1 < 0$ for all $q$), either $q_1 = -1$ or $q_1 = 0$ is assumed, in line with the stable, moderately dissipative model<fnref n="6" /> of Borderies et al. (1986). Note finally that both $t_1$ and $t_2$ (the latter being neglected here) $\propto q$ (or higher powers of $q$) when $q \to 0$ and diverge when $q \to 1$ (in order to prevent streamline crossing, which cannot occur in a fluid medium with a large enough collision frequency).

For future use, it is also useful to provide an expression for the viscous flux of angular momentum $L_H^{vis} = 2\pi a^2 a_{r\theta}$:

$$a_{r\theta}(q) = A(q)\nu_d\sigma_0 n, \qquad (39)$$

$$A(q) = \frac{q_a - q}{(1 - q^2)^{3/2}}, \qquad (40)$$

where $q_a$ specifies the level of nonlinearity required to achieve angular momentum luminosity reversal. The physical motivation of this model is the

<footnote n="6">⁶This makes the dependence on $q$ as $q \to 0$ more realistic than for the stable hydrodynamic model of Appendix A where the leading order terms in the Taylor expansion of $q$ cancel out somewhat fortuitously.</footnote>
same as for Eqs. (37) and (38); typically, $q_a \simeq 0.8$ in the dense granular flow limit, and this seems to apply in more realistic, less dense flows as well (Mosqueira 1996, Fig. 18). The hydrodynamic models of Appendix A display even larger (and probably unrealistic) values of $q_a$. The fact that collisions are dissipative implies that $q_1 < q_a$ (Borderies et al. 1985 for the granular flow model of dense rings; Longaretti 1992 for a more general setting).

The standard angular momentum flux expression $L_H = 3\pi a^2 \nu_0 \sigma_0 n$ is recovered in the $q=0$ limit by defining $\nu_0 = 2q_a\nu_d/3$ where $\nu_0$ is the ring effective viscosity in unperturbed flows. If one further requires consistency between this hybrid model for $a_{r\theta}$ and the models of Appendix A in the limit $q \to 0$, one must have $2q_a\nu_d/3 = \nu_{eff}$. This is assumed here for simplicity.

The model used above implicitly assumes that transport in the perturbed region close to the edge is dominated by the ring collisional stress tensor and not by self-gravity wakes. This seems at odds with the by-now commonly admitted fact that self-gravity wakes most probably dominate viscous angular momentum transport in the B ring unperturbed regions (away from the outer edge), as they are detected there (Colwell et al., 2007; Hedman et al., 2007; Nicholson and Hedman, 2010). In fact, self-gravity wakes are most probably quenched in the edge perturbed region (see section 4.2.3), which justifies the implicit assumption made here.

## 4. Semi-quantitative dynamics in the two-streamline approximation

The two-streamline approximation considerably simplifies the analysis of the dynamics. Such an approximation is analytically tractable and provides useful insights into the behavior of the more precise N-streamline numerical solutions that will be presented next. The two-streamline approximation was introduced for eccentric narrow rings by Borderies et al. (1983). It has been formulated for a single edge mode in Longaretti (2018) in a simplified manner, and the present section elaborates on this earlier analysis.

The two streamlines are characterized by their semimajor axes $a_1, a_2$, eccentricities $\epsilon_1, \epsilon_2$ and apsidal shifts, $m\Delta_1$ and $m\Delta_2$; the last two are collected in the complex eccentricities $Z_1$ and $Z_2$; $\delta a = a_2 - a_1$ is the characteristic width of the perturbed region; by construction, this is half the full extent between the inner and outer boundary of the perturbed region, with similar definitions for $\delta\epsilon$ and $\delta m\Delta$. The two streamlines are assumed to have equal
masses $m_1 = m_2 = m_r/2$, i.e., a constant surface density profile is assumed in the perturbed region, for simplicity.

Defining $\epsilon_B$ as the eccentricity at the edge, the total width of the perturbed region $\Delta a$ is constrained by two requirements. First, in order of magnitude, $q \simeq a\epsilon_B/\Delta a \lesssim 1$, but this requirement is rather crude. Second, one may assume $\Delta a \sim 2 \times (a_B - a_r)$ where $a_r$ is the resonance location of each mode and $a_B$ the B ring outer edge location. One expects the second requirement to be more relevant and precise than the first, so that each streamline has a width $\delta a = \Delta a/2 \sim (a_B - a_r)$.

The most important difference between a narrow ring and an edge mode bears on the relation between the mean eccentricity difference $\delta\epsilon_0$ across the structure and the mean eccentricity $\epsilon_0$. In the major narrow rings, $\delta\epsilon_0 \ll \epsilon_0$ while in an edge mode, $\delta\epsilon_0 \approx \epsilon_0$, because the eccentricity vanishes far enough from the edge. For simplicity $\epsilon$ is assumed to vary linearly from the inner boundary $a_i$ of the first streamline (where it is assumed to vanish) to the outer boundary of the second streamline $a_B$. The first streamline is therefore assumed to have $\epsilon_1 = \epsilon_B/4$ while the second one has $\epsilon_2 = 3\epsilon_B/4$. In this approximation,

$$ \delta\epsilon_0 \equiv \epsilon_2 - \epsilon_1, \qquad (41) $$

$$ \epsilon_0 \equiv \frac{\epsilon_1 + \epsilon_2}{2} = \delta\epsilon_0, \qquad (42) $$

The last equality is essentially a consequence of the inner boundary condition, discussed in more detail in Appendix C.2.

In the stationary limit, the two streamline approximation allows us to obtain two important ring diagnostics: the ring surface density and the effective viscosity in the mode region. The first one follows from the role of the ring self-gravity which prevents streamlines from losing phase coherence under the action of the planet forcing potential. The second one follows from relating the mode apsidal phase shift across its radial extent to the ring effective viscosity, as mechanical energy dissipation prevents all streamline apsidal shifts from being exactly aligned.

This analysis encompasses the whole dynamics when the ring viscosity effectively damps any possible motion around the stationary one. This is the case when the ring is viscously stable. In presence of viscous overstabilities, the analysis of this section applies to the stationary state around which librations take place. Consequently this analysis provides relevant semi-quantitative constraints, whether the ring is viscously stable or not.
The ring diagnostics just mentioned will be derived in a more precise way in sections 5 and 6; the semi-quantitative analysis exposed here is intended to provide insights into the underlying physics of these numerical solutions.

### 4.1. Free and forced modes stationary dynamics

<table><thead><tr><th>m</th><th>δa (km)</th><th>σ<sub>min</sub> (g/cm<sup>2</sup>)</th><th>σ<sub>max</sub> (g/cm<sup>2</sup>)</th></tr></thead><tbody><tr><td>2</td><td>33</td><td>120</td><td>240</td></tr><tr><td>3</td><td>25</td><td>135</td><td>270</td></tr><tr><td>4</td><td>19</td><td>120</td><td>240</td></tr><tr><td>5</td><td>23</td><td>250</td><td>500</td></tr></tbody></table>

Table 1: Estimates of the free nodeless edge mode surface density from Eq. (47), based on the data of Table 4 of French et al. (2023). In these estimates, the nodeless mode is identified with the mode with resonance closest to the edge, ignoring the forced $m = 2$ mode. The minimum estimate takes collective edge mode dynamical effects into account along the lines discussed in the text, while the maximum corresponds to a linear limit. The minimum value is more likely to be correct (see text).

The free and forced mode stationary dynamics can be analyzed along the lines of section 4.1 of Longaretti (2018). The equation for $d(Z_2 - Z_1)/dt = 0$ is the most important for the present analysis. It is obtained from Eqs. (28) and (36) and leads to, by separating the real and imaginary parts<fnref n="7" />:

$$ \frac{\delta a}{a} \approx \frac{2\Omega_{sg}}{3(m-1)\bar{n}}, \qquad (43) $$

$$ \delta(m\Delta) \approx -\frac{\lambda_1}{\Omega_{sg}}, \qquad (44) $$

where

$$ \begin{aligned} \Omega_{sg} &= \frac{\bar{n} m_r}{\pi M_p} \left(\frac{a}{\delta a}\right)^2 H(q^2) \\ &\approx 4H(q^2)\bar{n} \frac{\sigma_0 a^2}{M_p} \frac{a}{\delta a}, \end{aligned} \qquad (45) $$

$$ \lambda_1 = \frac{8\pi at_1}{q\bar{n}m_r\delta a} \approx \operatorname{sgn}(t_1) \frac{4\nu_{eff}}{\delta a^2}, \qquad (46) $$

<footnote n="7">⁷The quantity $\delta(m\Delta)$ has no relation to the quantity $\delta m\Delta$ of Appendix C.1.1.</footnote>
and $m_r$ is the mass in the radial extent of the perturbed region associated with the edge mode. Also, the relation $\delta\epsilon \simeq \epsilon$ has been used. These equations are simplified versions of Eqs. (10.95) and (10.96) of Longaretti (2018). Note that by construction, the satellite forcing disappears from these equations for forced modes, which is why they apply to both free and forced modes.

The first of these equations, Eq. (43), provides a relation between the ring surface density and the size of the resonant cavity $\delta a$. Eq. (44) on the other hand allows us to relate the apsidal shift difference across the edge mode radial extent to the effective viscosity of the ring in the perturbed region. The sign of the apsidal shift change across the edge mode region $\delta(m\Delta)$ is related to the sign of $\lambda_1$.

More precisely, Eq. (43) can be recast to deduce the local surface density in the edge mode region with the help of Eq. (45):

$$ \sigma_0 \sim \frac{3(m-1) M_p}{8H(q^2)} \left(\frac{\delta a}{a}\right)^2, \qquad (47) $$

This relation is not very precise due to the fact that $\delta a$ is not precisely constrained, as a consequence of the fact that the eccentricity $\epsilon$ decreases slowly with decreasing $a$ inside the mode resonance location. The scaling with $m$ and $\delta a$ is robust, however.

This estimate applies to the free nodeless $m \neq 1$ edge modes$^8$; these are identified with the mode whose resonance is closest to the edge for $m = 3, 4, 5$ and the second closest to the edge for $m = 2$ (the closest to the edge being the mode forced by Mimas), in Table 4 of French et al. (2023). For the forced mode, this relation is even less precise, as both $\delta a$ and $q$ are smaller than for the free $m = 2$ mode. It is therefore not used for the forced mode; instead, a more precise estimate is derived in the next subsection for the forced mode, from the observed amplitude of this mode at the edge.

For the time being, and from $\delta a \sim (a_B - a_r)$, Eq. (47) leads to the results collected in Table 1 for the surface density estimated from nodeless edge modes. Two different assumptions for $H(q^2)$ are adopted in this estimate. First, the linear limit ($q \to 0$) is assumed, i.e., $H(q^2) = 1/2$; this is in line with the fact that quite a number of the nodeless edge modes seem to have relatively modest amplitude at the edge, suggesting that their maximum $q$ is rather small, an expectation borne out in the more precise numerical $N$

<footnote>$^8$A two-streamline approximation can only capture nodeless edge modes.</footnote>
streamlines estimated made later on in this paper. Second, one may argue on the contrary that the relevant value for the mean $q$ in this relation is a collective one and furthermore, cannot be very different, on average in time, from the value required by viscous angular momentum flux reversal $q_a \simeq 0.8$, which gives $H(q_a^2) \simeq 1$. The two estimates of the surface density resulting from these two choices are noted $\sigma_{max}$ and $\sigma_{min}$ respectively in Table 1. Considering that the requirement $q = q_a$ must eventually be met for the edge to be actually confined by Mimas, the minimum value is more likely to be correct than the maximum one; a look at the last column of Table 2 shows that this $q = q_a$ requirement should be easily met when taking into account the presence of multiple edge modes at the B ring edge. Indeed, in first approximation, this changes $q$ into an effective value $q^2 \simeq \sum_{\alpha} q_{\alpha}^2$, where $\alpha$ is an index for each individual mode<fnref n="9" />.

The values obtained for $m = 2, 3, 4$ are consistent with each other. Although the approximation used here is in principle uncertain within a factor of at least a few, it turns out that these estimates are surprisingly precise when compared to numerical results. However, the surface density obtained for $m = 5$ seems off. This is related to the breaking of the trend of decreasing $\delta a$ with increasing $m$, which suggests that the nodeless mode may not have been detected for this $m$ value.

One can derive a similar result for $m = 1$ by using the correct form of $m(\Omega - \Omega_p) - \kappa$ for this mode. This leads to

$$ \sigma_0 \sim \frac{21}{8} J_2 \frac{M_p}{a^2} \left(\frac{R_p}{a}\right)^2 \left(\frac{\delta a}{a}\right)^2, \qquad (48) $$

and yields $\sigma_0 \sim 75$ g/cm² for the $m = 1$ nodeless free mode, with $\delta a \sim 150$ km from Table 4 of French et al. (2023). No collective correction for $H(q^2)$ has been applied here, for two reasons. First, the $m = 1$ mode has a very large cavity, while the collective edge mode effect is relevant much closer to the edge; second, and precisely because of this extensive cavity, the mean value of $q$ for this mode is expected to be rather small, even taking into account that several $m = 1$ modes are detected at the B ring edge. Both arguments suggest that the linear limit $H(q^2) = 1/2$ should be appropriate for this mode.

<footnote n="9">⁹A partial but sufficient justification of this statement can be found in Longaretti 1989. A more extended analysis will be given in a future paper.</footnote>
Let us conclude this section by examining the implication of Eq. (44). Combined with Eqs. (45) and (46), this yields an estimate of the effective viscosity in the perturbed region from the observed apsidal shift of the forced mode:

$$ \nu_{eff} \sim \frac{n\delta a^2 \sigma_0 a^2 a}{2 M_p \delta a} \delta(m\Delta) \\ \sim 8 \times 10^4 \left( \frac{\Delta_B}{2.9^\circ} \right) \left( \frac{\sigma_0}{200 \text{ g/cm}^2} \right) \left( \frac{\delta a}{15 \text{ km}} \right) \text{cm}^2/\text{s}, \quad (49) $$

where $m\Delta_B$ is the observed apsidal shift of the forced mode with respect to Mimas mean longitude. The forced mode $m = 2$ and cavity size $\delta a$ have been used in the second relation. The implication of this result will be discussed in the next subsection.

## 4.2. Forced mode amplitude and edge torque

For forced modes, one can relate the mode amplitude at the edge $\epsilon_B$ to the magnitude of the forcing on the one hand (a relation involving the ring surface density, due to the dependence of forced density wave wavelengths on the surface density), and the apsidal shift at the edge $m\Delta_B$ — and correlative effective viscosity — to the satellite torque on the other hand. Further relating the viscosity in the satellite perturbed region to the unperturbed viscosity allows for self-consistency checks on the assumptions involved and of the magnitude of the observed apsidal shift, with respect to the expectation that viscous transport in the unperturbed B ring region is dominated by self-gravity wakes.

### 4.2.1. Forced amplitude and related effective viscosity

The only mode of interest in this subsection is the edge mode forced by the Mimas 2:1 resonance. The relevant new relations one can obtain for this mode follow from Eqs. (29) [or (31)] and (36) by focusing on the outermost streamline in the two streamline approximation. Eq. (29) leads to

$$ \begin{aligned} & \frac{m_r}{2\pi M_p} \left(\frac{a}{\delta a}\right)^2 H(q^2) \frac{\epsilon_B}{2} + \frac{\omega(a_2) 3\epsilon_B}{n} \frac{4}{2} \\ &= \left(\frac{a_r}{a_s}\right) \left(\frac{M_s}{M_p}\right) \left(\frac{H_{mk}}{2}\right) = \left(\frac{M_s}{M_p}\right) \left(\frac{m-1}{m}\right)^{5/3} H_{mk} \\ & \approx \frac{3M_s}{4M_p}, \end{aligned} \quad (50) $$
where $\omega(a_2)/n = [m(\Omega(a_2) - \Omega_p) - \kappa(a_2)]/n \simeq -3(m-1)\delta a/4a$ and $H_{mk}$ is defined by Eq. (B.7). The last approximate relation is coincidental for $m=2$.

Recalling that $m_r = 4\pi a\delta a\sigma_0$ and $a\epsilon_B \sim 2\delta a$, one finds

$$ \sigma_0 \sim \frac{3M_s/4 + 9M_p(\delta a/a)^2/8}{2a_r^2 H(q^2)}. \qquad (51) $$

This leads to $\sigma_{min} \simeq 130 \text{ g/cm}^2$ and $\sigma_{max} \simeq 260 \text{ g/cm}^2$ for $H(q^2) = 1/2$ and $H(q^2) = H(q_a^2) \simeq 1$, respectively. This is broadly consistent with the results collected in Table 1, considering the uncertainty involved.

Turning to Eq. (36) yields

$$ \begin{aligned} & \left[ \frac{m_r}{2\pi M_p} \left(\frac{a}{\delta a}\right)^2 H(q^2) \frac{\epsilon_B}{2} + \frac{\omega(a_2) 3\epsilon_B}{n} \frac{1}{4} \right] \sin m\Delta_B \\ & = -\operatorname{sgn}(t_1) \frac{4\pi}{nm_r} \nu_{eff}\sigma_0. \end{aligned} \qquad (52) $$

Taking into account Eq. (50), this equation reduces to

$$ \begin{aligned} \nu_{eff} &\simeq -\operatorname{sgn}(t_1) \frac{M_s}{M_p} m n a \delta a \sin m \Delta_B \\ &\simeq -2 \operatorname{sgn}(t_1) \frac{M_s}{M_p} m n a \delta a \delta(m \Delta), \end{aligned} \qquad (53) $$

where the large $m$ limit has been taken $[(m-1)/m \simeq 1, H_{mk} \simeq m]$. This limit overestimates the actual result by a factor of order $\sim 2$ for the Mimas $m=2$ forced edge mode.

The joint constraints of the exchanges of angular momentum and energy between the satellite and the ring in the forced mode region implies that the effective viscosity in the forced edge mode region, $\nu_{eff}$, and in the inner unperturbed region, $\nu_0$, are related by (Borderies et al., 1982; Longaretti, 2018)

$$ \nu_{eff} \simeq \frac{\nu_0}{m} \left(\frac{M_p}{M_s}\right)^{1/2} \simeq 2 \times 10^3 \nu_0. \qquad (54) $$

To the same level of approximation, $\delta a \sim a(M_s/M_p)^{1/2}$ (Borderies et al., 1982). Injecting these last two relations in Eq. (53) finally yields

$$ \nu_0 \simeq -\operatorname{sgn}(t_1) \left(\frac{M_s}{M_p}\right)^2 n m^2 a^2 \sin m \Delta_B, \qquad (55) $$
Note that this equation implies that $m\Delta_B$ and $t_1$ have opposite signs as $\nu_0 > 0$.

### 4.2.2. Satellite torque and effective viscosity

Yet another relation between the effective unperturbed viscosity and the apsidal shift at the edge has been derived by Borderies et al. (1982) from the requirement of balance of the unperturbed viscous angular momentum flux and the satellite torque [proportional to $\delta(m\Delta)$ for a forced edge mode]. Let us recall this derivation, due to its relevance for the present discussion.

The satellite torque reads (see, e.g., Longaretti 1992, 2018 and references therein)

$$ 
\begin{aligned}
T_s &= \int da \pi m a \sigma_0 \Psi_{mk} \epsilon \sin m\Delta \\
&\approx -(mna)^2 \sigma_0 \frac{M_s}{M_p} (a\epsilon_B) \Delta a \sin m\Delta_B \\
&\approx -m^2 n^2 a^4 \sigma_0 \left(\frac{M_s}{M_p}\right)^2 \sin m\Delta_B.
\end{aligned} 
\quad (56) $$

In the last relation, $a\epsilon_B \approx \Delta a \approx a(M_s/M_p)^{1/2}$ has been used. This makes the torque expression identical to the approximation $T_{BGT}$ used by Borderies et al. (1982).

Balancing this torque with the unperturbed angular momentum flux $L_h^{vis} = 3\pi\sigma_0 na^2\nu_0$ yields:

$$ \nu_0 \approx \frac{m^2}{3\pi} \left(\frac{M_s}{M_p}\right)^2 na^2 \sin m\Delta_B. \quad (57) $$

Comparing Eqs. (55) and (57) shows that there is about an order of magnitude of discrepancy between the two estimates, due to the crudeness of the approximations involved.

The fact that $\nu_{eff} \gg \nu_0$ has a consequence that was first pointed out by Borderies et al. (1984). From Eq. (39), the viscous angular momentum flux in the perturbed region would largely exceed the satellite torque, unless $q - q_a \ll 1$ or $\sigma_0$ is much smaller than in the unperturbed region. Because the outer edge of the B ring is rather sharp, $\sigma_0$ remains large very close to the edge and $q \approx q_a$ necessarily applies in the perturbed region of the forced mode. This was used above to take into account collective effects in the estimate of the ring surface density in the narrow edge mode region.
The satellite torque is necessarily negative for the satellite to confine the edge as the angular momentum flux due to the stress tensor torque in the unperturbed region (a positive quantity) must be balanced by the satellite torque, i.e., angular momentum must flow from the ring to the satellite. This in turn implies that $m\Delta_B > 0$ with the conventions on the origin of time and angles adopted earlier, i.e., $\Delta_B - \varphi_s > 0$ for an arbitrary choice of origins.

Consequently, Eq. (55) and Eq. (57) imply that $t_1 < 0$ in the perturbed region, and libration motions (oscillations in eccentricity and apsidal shift) are damped (see Longaretti and Rappaport 1995 and section 4.1 of Longaretti 2018). Possible ways around this conclusion will be discussed in a forthcoming paper on multiple modes time-dependent dynamics.

### 4.2.3. Viscosity constraints and wake transport

Observations imply that $\Delta_B \sim 1^\circ$ (French et al., 2023), so that the two estimates above for the unperturbed viscosity give $\nu_0 \sim 2 - 20 \text{ cm}^2/\text{s}$. However, values as high as $\Delta_B \sim 2.9^\circ$ have been reported in the literature (Spitale and Porco, 2010) although a later analysis of the B ring edge modes concludes that this might be an upper limit (Nicholson et al., 2014). From Eqs. (55) and (57), this gives $\nu_0 \sim 12 - 120 (\Delta_B/2.9^\circ) \text{ cm}^2/\text{s}$. Note finally that the rather crude estimate of Eq. (49) combined with Eq. (54) gives $\nu_0 \sim 50 (\Delta_B/2.9^\circ) \text{ cm}^2/\text{s}$ for the forced mode, in the relevant range.

The lower bound on $\nu_0$ just derived is safely above the ring minimal viscosity $\nu_m \sim nd^2$ where $d$ is the ring particles diameter ($d \sim 1 \text{ m}$ so that $\nu_m \sim 1 \text{ cm/s}^2$). We now show that the estimate on $\nu_0$ is also consistent with the now prevalent idea that transport in the unperturbed region is dominated by self-gravity wakes. Indeed,

$$ \nu_{sg} \approx C \frac{G^2 \sigma_0^2}{n^3}, \qquad (58) $$

where $C$ is a dimensionless quantity with $C \simeq 20 - 40$ in the vicinity of the B ring edge, for pure ice particle density (0.9 g/cm$^3$) (Daisaka et al., 2001; Yasui et al., 2012). Reducing the particle density by a factor of 2 reduces $C$ by a factor of 3, finally bringing the relevant range down to $C \simeq 7 - 14$; this lower particle density is favored by the analysis of the ring brightness asymmetry (French et al., 2007). It is adopted in this paper as well as the related, lower range of $C$.

From Eq. (48), setting $\sigma_0 = 75 \text{ g/cm}^2$ leads to $\nu_{sg} \simeq 20 - 40 \text{ cm}^2/\text{s}$ from Eq. (58) with $C \simeq 7 - 13$. Thus, the constraints on $\nu_0$ obtained in the
previous subsection are compatible with the existence of wake transport. The discussion of this point in section 6.2, based on more precise numerical results for the various estimates provided here, essentially confirms this inference.

Let us conclude this section by pointing out that self-gravity wakes are most probably quenched in the forced mode perturbed region (the last 15 to 30 km inside the B ring edge). Indeed, wakes are present when the Toomre parameter is low enough$^{10}$, $Q = c\kappa/\pi G\sigma_0 \lesssim 2$ (Toomre, 1964). On the one hand, the required enhanced dissipation in the edge region implied by Eq. (54) most probably translates into a velocity dispersion driving $Q$ sub-stantially above the instability threshold; on the other hand, if close-packing of the ring particles is sufficiently strong, the enhanced dissipation implied by Eq. (54) may result from an increased ring thickness instead (and hence a larger number of collisions in the edge region), possibly maintaining Toomre's criterion. However, the surface density derived in this region (here and in section 6.2) implies that a close-packed many particle thick ring configuration is unlikely in the forced mode perturbed region; in this case, the quenching of self-gravity wakes is a likely outcome.

### 4.2.4. Summary

The number of different estimates of the ring viscosity derived above and their relations to one another is somewhat labyrinthine. For this reason, a short summary of the preceding analysis is presented here.

* First, an estimate of the forced mode effective viscosity in the near-edge perturbed region $\nu_{eff}$ is obtained [Eq. (53)], from the dynamical equations. Its magnitude is observationally constrained by the observed apsidal shift of the forced $m = 2$ mode.

* Next, two different estimates of the viscosity in the unperturbed region $\nu_0$ are derived. The first one [Eq. (55)] comes from the estimate of $\nu_{eff}$ and from the consistency of the budgets of energy and angular momentum exchanges with the satellite, which implies a large increase

<footnote>$^{10}$This same criterion implies that wakes are possibly unstable at surprisingly low surface densities. Indeed, $\nu_{sg} > \nu_m$ implies that $\sigma_0 \gtrsim 15$ g/cm$^2$ (for $C$ in the middle of the admissible range) while assuming the minimum velocity dispersion $c_m \simeq nd$ in Toomre criterion leads to $\sigma_0 \gtrsim 10$ g/cm$^2$ at the B ring edge. Both constraints are roughly compatible (as they should), and indeed suggest that wakes exist at low surface densities in the B ring edge region.</footnote>
in dissipation in the perturbed region with respect to the unperturbed one [Eq. (54)] (Borderies et al., 1982). The second one [Eq. (57)] comes from the required balance of the satellite torque with the unperturbed viscous angular momentum flux for the satellite to be able to confine the ring. Both estimates are also constrained from the observed apsidal shift and should in principle be identical if they were obtained in an exact manner. In practice, due to differences in approximations between the two estimates, there is a factor of about 10 of difference, leading to $v_0 \sim 12 - 120 (\Delta_B/2.9^\circ) \text{ cm}^2/\text{s}$.

• Finally, these two estimates are compared with the viscosity derived under the assumption that gravitational wakes dominate angular momentum transport in the unperturbed region [Eq. (58)]. This last viscosity estimate is constrained by the surface density derived above from the $m=1$ mode. This leads to $v_{sg} \approx 20 - 40 \text{ cm}^2/\text{s}$.

This analysis shows that the ring transport properties derived from the data through the analysis of the forced mode dynamics is consistent with the idea that self-gravity wakes dominate angular momentum transport in the unperturbed near-edge region. A similar but more precise discussion based on numerical results will be provided in section 6.2.3.

**5. Numerical solution of single mode stationary dynamics: free modes**

The numerical results presented in this work bear only on stationary modes. The time-dependent dynamics is much more involved from a numerical point of view, and will be dealt with elsewhere. Also, only single modes are considered.

As a consequence, the surface density estimates used to best approximate the observed mode properties (cavity size and amplitude at the edge) ignores the collective $q$ effect discussed in relation to Table 1, although this approximation probably leads to overestimate the surface density by a factor $\sim 2$, as discussed in section 4.1.

In such numerical solutions, streamlines are equally spaced between the inner boundary $a_{in}$ and the ring edge $a_B$. In practice, $N = 200$ streamlines are used for almost all simulations as a compromise between the speed of
computations and the precision of the results (relative precision of a few percent at worse for the eccentricities $\epsilon$ of the edge modes of interest). Unless otherwise specified (in particular when exploring the first twelve eigenmodes for $m = 1$), $a_B - a_{in} = 2\lambda_{dw}$ is used in all numerical solutions ($\lambda_{dw}$ is the first wavelength of a freely propagating density wave; see Longaretti 1992, Eq. 7.43). The boundary conditions are discussed in Appendix C. For the search of the twelve first $m = 1$ free modes performed below, the domain has been extended to $4\lambda_{dw}$ and the number of streamlines to $N = 400$, to the same level of precision. This domain is sufficient to capture a large enough evanescent zone for all modes.

A few runs have been made for $N = 400$ and 800 to quantify the precision of the numerical results for the first three modes for $m \neq 1$. The least precise quantity is the value of $q$ at the edge, which does not quite converge to zero for $N = 200$ (as it should from the discussion of Appendix C.1.1 and Appendix C.1.2), but this affects only marginally the precision of the resonance location ($a_B - a_r$) and the derived eccentricities and pattern speeds (a few percent at most of relative errors in both cases).

## 5.1. *Numerical procedure*

Eq. (30) is solved for the eigenvectors (eccentricities of the *N* streamlines) and eigenvalue ($m\Omega_p$) in the linear limit [$H(q^2) = 1/2$] for a given *m* and mass profile by standard dedicated linear algebra routines. In practice a constant mass profile is used and only the associated constant ring surface density must be specified. The resulting linear solution is then used as a guess for an iterative search of the nonlinear solution; the iterations proceed until the solution is converged with respect to a relative change in the eccentricity profile of the order of $10^{-6}$.

Both linear and nonlinear solutions are first normalized by fixing the maximum $q$ obtained along the profile in the linear solution. The amplitudes $a_B\epsilon_B$ at the edge quoted in tables below are obtained for a maximum nonlinearity parameter $q$ that is adjusted in order to match the observed free mode amplitudes $a_B\epsilon_B$ for $m = 2$ to $m = 5$. Next, the ring surface density is found by successive approximations until the resonance locations of the nodeless modes match the observations. Finally, the parameter $\Gamma$ of Eq. (62) is found from the common surface density that best matches the resonance locations of all these nodeless modes. A separate search of the same type is performed
for the $m = 1$ edge mode<fnref n="11" />.

These best matches of surface density were not determined in a very precise way (not better than ~ 20%), considering that the assumptions of constant density profile and mode independence are not very precise.

The overall arbitrariness of the normalization of the eccentricities is obvious in the linear regime. It is somewhat less obvious in the nonlinear one. In fact, numerical solutions show that the nonlinearity, embedded in $H(q^2)$, is mild enough as to produce only a mild distortion of the linear solution, for all normalizations. In particular, the nonlinearity does not specify a unique solution for a given number of radial nodes, but a continuum, as in the linear regime. One can choose, e.g., the amplitude at the edge $a_B\epsilon_B$ (within a bound ensuring that $q < 1$ everywhere), or the maximum $q$ of the mode as the free parameter of each family of solutions. In practice this maximum $q$ is set in a rather direct way for the linear solutions, but the reported values pertain to the nonlinear ones obtained from these linear solutions; the difference between the two is negligible except for the free $m = 2$ nodeless mode.

The actual observed amplitudes of edge modes are ultimately determined by processes not discussed here. If the rings are viscously stable, free edge modes will eventually be completely damped on relatively short time-scales (of the order of a few hundred years). Conversely, the fact that free modes are observed strongly suggests that the rings are viscously overstable, at least part of the time. In this case, the amplitude of a single mode is specified by the requirement that the viscous coefficient $t_1$ vanishes on average over the edge mode region. The consequences of this requirement for edge mode amplitudes are however not straightforward due to the simultaneous presence of quite a few different edge modes at the B ring edge.

In Eq. (28), the eccentricities and the pattern speed are not independent quantities. Multiplying this equation by $m_i$ and summing over $i$ produces a cancellation of the self-gravity terms and leads to

$$m\Omega_p = \frac{\sum_i m_i \epsilon_i (m\Omega_i - \kappa_i)}{\sum_i m_i \epsilon_i} \qquad (59)$$

This relation is satisfied by necessity by the linear eigenvalues. In the nonlinear regime, this relation is enforced in the iterative process establishing the

<footnote n="11">¹¹As this mode extends much farther inside the B ring, due to the fact that its pattern speed is $O(J_2\bar{n})$, the characteristic surface density may be different.</footnote>
solution from the linear one.

## 5.2. Numerical results

### 5.2.1. Analysis in the independent mode limit

For illustrative purposes, the first free eigenmodes are shown for $m = 3$, $\sigma_0 = 200$ g/cm$^2$ and a maximum $q = 0.2$ in the linear regime on Fig. 1. The figure shows the nonlinear and linear solutions are almost identical for such a small maximum value of $q$. Note that the sign of the amplitude $a\epsilon$ is arbitrary in these solutions (see also section 2.2 on this point).

For the same purpose, Fig. 2 shows the influence of the number $N$ of streamlines on the solution for the same modes under the same choice of physical parameters. The nonlinear solutions are nearly identical. The only (barely) noticeable difference is that at the edge, $q$ for $N = 400$ is about a factor of two smaller than for a $N = 100$, in line with the discussion of Appendix C.1.2.

This exploration of basic mode structure is concluded with Fig. 3, which displays the $m = 1$ and $m = 2$ mode amplitudes and nonlinearity parameters in the linear and nonlinear regimes. The $m = 1$ mode is substantially more extended, as already pointed out. The $m = 2$ mode exhibits a much more pronounced difference between the linear and nonlinear solutions, due to the larger maximum $q$ chosen to fit its observed edge amplitude, while the linear and nonlinear solutions of the $m = 1$ mode are indistinguishable due to the low $q$ needed to match the edge amplitude.

A rather large number of numerical solutions (a few tens) has been explored to find the best matching solutions for free nodeless modes, when compared to the results of the data analysis of French et al. (2023). The match is searched for two quantities, the mode cavity size $\Delta a = a_B - a_r$ and the mode amplitude at the edge, $a_B\epsilon_B$; for each $m$, the nodeless mode is assumed to be the free mode with smallest cavity size$^{12}$. This assumption is natural in the limit of independent, stationary edge modes, but is less obvious when time-dependent dynamical effects are taken into account. This is briefly discussed in section 7.

In the linear limit, only the cavity size depends on the surface density and the problem boils down to finding a single surface density making the nu-

<footnote>$^{12}$Note that for $m = 2$, the smallest cavity pertains to the forced mode, but this mode is excluded by definition from the present analysis. It is analyzed in the next section.</footnote>
<figure><img src="image_2.png" /><figcaption>Figure 1: The $m = 3$ edge mode amplitude profiles $a\epsilon(a)$ and $q(a)$, for the first three eigenmodes, $n_r = 0$ (blue), $n_r = 1$ (red) and $n_r = 2$ (black), where $n_r$ is the number of radial nodes of the mode (number of half wavelengths of the trapped wave). Nonlinear solution: plain line, linear solution: dashed line; $\lambda_{dw} = 94$ km and $\sigma_0 = 200$ g/cm$^2$. The origin of radii is the outer edge, and the resonance location of each mode is materialized by a vertical line in the first graph.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>Figure 2: The $m = 3$ edge mode amplitude profiles $a\epsilon(a)$ and $q(a)$, for the first three eigenmodes, $n_r = 0$ (blue), $n_r = 1$ (red) and $n_r = 2$ (black). Number of streamlines: 400 (plain line), and 100 (dashed line); $\lambda_{dw} = 94$ km and $\sigma_0 = 200$ g/cm$^2$.</figcaption></figure>
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>Figure 3: The $m = 1$ (left) and $m = 2$ (right) edge modes amplitude profiles $a\epsilon(a)$ (top) and $q(a)$ (bottom), for the first three eigenmodes, $n_r = 0$ (blue), $n_r = 1$ (red) and $n_r = 2$ (black). Nonlinear solutions: plain lines, linear solutions: dashed line; $\lambda_{dw} = 133$ km and $\sigma_0 = 200$ g/cm$^2$ ($m = 2$); $\lambda_{dw} = 594$ km and $\sigma_0 = 60$ g/cm$^2$ ($m = 1$).</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 4: The dependence on surface density (top) and maximum *q* (bottom) of the amplitude (left) and cavity size (right) of the nodeless *m* = 3 edge mode. For varying *q*<sub>max</sub>, σ<sub>0</sub> = 200 g/cm<sup>2</sup>; for varying σ<sub>0</sub>, *q*<sub>max</sub> = 0.2. The various black dashed-lines represent physically motivated fits to the numerical results (see text).</figcaption></figure>
merical solution cavity sizes match at best the observed ones for all nodeless modes at once. Indeed, in the linear limit, the eccentricity at the edge can be adjusted independently by scaling the eccentricity profile. As pointed out semi-quantitatively in section 4.1, a mode cavity size increases with increasing surface density $\sigma_0$ and this has been used to find the best fitting common surface density for the cavity sizes of the nodeless $m = 2$ to 5 modes in the linear regime.

This is no longer true in the nonlinear regime, where $q$ also influences the cavity size. Although most modes are found to be only weakly nonlinear, this effect cannot be completely ignored, in particular for the $m = 2$ free mode.

Qualitatively, a change of $\sigma_0$ has substantial effect on both the cavity size $\delta a$ and $a_B\epsilon_B$, while a change of $q_{max}$ has substantial effect on $a_B\epsilon_B$ but only a weak effect on $\delta a$. These trends are displayed in a more quantitative manner on Fig. 4. The dependencies displayed on this figure can be understood from the two following scaling relations:

$$ \delta a \propto [\sigma_0 H(\tilde{q}^2)]^{1/2}, \qquad (60) $$

$$ a_B \epsilon_B \propto \tilde{q} \delta a \propto \tilde{q} [\sigma_0 H(\tilde{q}^2)]^{1/2}, \qquad (61) $$

where $\tilde{q} \lesssim q_{max}$ is a characteristic value of $q$ in the mode cavity. The first scaling comes from Eq. (47) and the second one from the definition of $q$ which can be approximated by $\tilde{q} \simeq a_B\epsilon_B/2\delta a$. Note that holding $q_{max}$ or $\tilde{q}$ constant is equivalent.

Let us focus first on the dependence of the cavity size and mode amplitude at the edge on $\sigma_0$ at constant $q_{max}$, i.e., on the two top sub-graphs of Fig. (4). The two preceding scaling relations imply that $\delta a = \delta a^*(\sigma_0/\sigma_0^*)^{1/2}$ and $a_B\epsilon_B = a_B\epsilon_B^*(\sigma_0/\sigma_0^*)^{1/2}$ where $(\sigma_0^*, \delta a^*)$ and $(\sigma_0^*, a_B\epsilon_B^*)$ stand for a reference point in each of these two graphs. The blacked dashed lines show these two relations, where $\sigma_0^* = 50 \text{ g/cm}^2$ has been chosen.

Turning now to the two bottom sub-graphs (dependence on $q_{max}$ at constant background surface density), the same argument leads to $\delta a = \delta a^*[H(\tilde{q}^2)/H(\tilde{q}^{*2})]^{1/2}$ and $a_B\epsilon_B = a_B\epsilon_B^*(\tilde{q}/\tilde{q}^*)[H(\tilde{q}^2)/H(\tilde{q}^{*2})]^{1/2}$. The black dashed lines show again these two relations with $\sigma_0^* = 50 \text{ g/cm}^2$. In these scalings $\tilde{q} = q_{max}/\sqrt{2}$ has been used (i.e. $\tilde{q}^2 = q_{max}^2/2$; this value was found by trial and error). The fit of the mode amplitude is not strictly linear due to the weak dependence of $H$ on $\tilde{q}$ that is visible on the fit of the cavity
size. Note in this respect that $q_{max} < 1$ implies that $\tilde{q} < 0.71$. For such a value, the $(H/H^*)^{1/2}$ scaling factor is still quite modest. Because of this, the independence of the cavity size on $q$ in the linear limit translates into a weak dependence on the nonlinear regime.

This exercise shows that although Eq. (47) is only semi-quantitatively correct, the scalings implied by this relation are rather precise. The least precise is the amplitude one as a function of $q_{max}$ (lower left graph), and this is clearly due to the use of the rather crude approximation of $\tilde{q}$ from the definition of $q$ instead of some more precise but unknown average on the whole $q$ profile; this also explains why the deviation increases with increasing $q_{max}$ on the lower left graph.

Considering the limitation of the present analysis, mostly due to the neglect of the modes mutual interactions through the ring self-gravity, but also to possible variations of the surface density profile in the edge mode region, the best matching numerical solutions of Table 2 have not been determined to better than ~5 - 10% in relative precision in $\sigma_0$ and $q_{max}$. These best matching solutions have been found by a trial and error process taking into account the dependencies just described. This process converges rather fast at the level of precision just mentioned.

The results are collected in Table 2. The best fitting surface density for these modes is $\sigma_0 \simeq 200$ g/cm$^2$ (in the independent mode dynamics limit discussed here). The best match for the $m=1$ mode is obtained for a lower surface density, $\sigma_0 \simeq 60$ g/cm$^2$. This follows because this mode has a much wider cavity, in line with the fact that the first wavelength of the underlying trapped density wave is much longer (see footnote 11). Also, the surface density at the edge may be enhanced compared to deeper into the ring due to, e.g., accretion of material from the Cassini Division$^{13}$. For this reason, a specific surface density has been searched for to match the $m=1$ nodeless mode cavity size.

The nodeless $m=5$ mode is the outlier in this table. The mismatch between the expected and measured cavity size is substantially larger than for all other modes, both in relative and absolute value. Also this mode breaks the expected trend of decreasing cavity size with increasing $m$, an oddity already pointed out in section 4.1. A possible explanation is that the

<footnote>$^{13}$As, for example, drift velocities due to ballistic transport increase with decreasing surface density.</footnote>
<table><caption>Table 2: Comparison of the nodeless modes resonance location ($\Delta a = a_B - a_r$) and amplitude at the edge ($a_B\epsilon_B$) determined from the Cassini data analysis (French et al., 2023) and the best matching numerical solutions, for a constant surface density profile $\sigma_0 = 200 \text{ g/cm}^2$ ($m = 2$ to 5) and $\sigma_0 = 60 \text{ g/cm}^2$ ($m = 1$). The numerical maximum value of $q = q_{max}$ corresponding to the best fitting value of the edge amplitude in the nonlinear regime is also given. Distances are rounded to the nearest kilometer, $q_{max}$ to the first decimal place.</caption><thead><tr><th>m</th><th>Δa (km)<br>(data)</th><th>Δa (km)<br>(theory)</th><th>a<sub>B</sub>ε<sub>B</sub> (km)<br>(data)</th><th>a<sub>B</sub>ε<sub>B</sub> (km)<br>(theory)</th><th>q<sub>max</sub></th></tr></thead><tbody><tr><td>2</td><td>33</td><td>35</td><td>38</td><td>35</td><td>0.6</td></tr><tr><td>3</td><td>26</td><td>23</td><td>7</td><td>7</td><td>0.2</td></tr><tr><td>4</td><td>19</td><td>19</td><td>3</td><td>3</td><td>0.1</td></tr><tr><td>5</td><td>23</td><td>16</td><td>5</td><td>5</td><td>0.2</td></tr><tr><td>1</td><td>151</td><td>144</td><td>24</td><td>23</td><td>0.1</td></tr></tbody></table>

$m = 5$ nodeless mode had too small an amplitude to be detected at the time of the Cassini mission.

It is convenient for data analysis to define a quantity $\Gamma$ such that Eq. (47) can be put in a more precise form (for $m > 1$):

$$ \sigma_0 = \frac{3(m-1)(a_B - a_r)^2 M_p}{8\pi^2\Gamma a_r^4} \frac{1}{2\langle H(q^2) \rangle}. \quad (62) $$

This definition is adopted for consistency with French et al. (2023). In this expression, and in line with Eq. (47), $\langle H(q^2) \rangle$ stands for some characteristic (average) value of $H$ over the mode profile. In the present context, and considering the modest maximum $q$ value obtained for these modes (see Table 2) and the rather weak dependency on $q$ shown on the lower right graph of Fig. 4, the $H(q^2) = 1/2$ linear limit may be applied (even for $m = 2$) for the purpose of an overall fit of $\sigma_0$ to a sufficient level of approximation, in which case $2\langle H(q^2) \rangle \simeq 1$ and$^{14}$ $\Gamma \simeq 1/16$ must be chosen in order to best reproduce the results obtained from the numerical solutions.

The quantity $\Gamma$ has a geometric interpretation that is discussed in Appendix C.1.1. Namely, this quantity embodies the constraint that the total

<footnote>$^{14}$The fact that this estimate of $\Gamma$ is the approximately the inverse of a power of 2 is coincidental.</footnote>
phase change along a complete outgoing/ingoing wave cycle must be a mul-
tiple of $2\pi$. This can also be seen as a form of resonant condition for the
ingoing and outgoing wave to add up constructively and not destructively,
in order for the mode to exist.

It is of direct interest to the present investigation to look into the reso-
nance locations of free modes of higher number of radial nodes $n_r$. This has
been done for the best matching solution for $\sigma_0$, and conventionally keeping
the same $q_{max}$ for all eigenmodes. This last choice does not affect the con-
clusion below, due to the limited impact of this choice on the resonant cavity
sizes of the modes.

These computations have been made under the hypothesis that the vari-
ous mode-like features found by French et al. (2023) for each *m* correspond to
actual edge modes of higher radial wavenumbers. To this effect, the observed
mode-like features are reordered in increasing cavity size order, and this is
assumed to correspond to consecutive values of the number of radial nodes
*n*<sub>r</sub> (the treatment of the *m* = 1 mode is again different, as discussed right
below). The results of these numerical computations are collected in Table
3. The nodeless solutions already collected in Table 2 are reproduced here,
for comparison purposes.

The results are unequivocal. For *m* = 2 to 5, it is impossible that the
various mode-like features detected besides the nodeless one correspond to
higher numbers of radial nodes, as the observationally determined cavity sizes
are always too narrow by a wide margin.

Conversely, the *m* = 1 modes may be compatible with the idea that
they correspond to different radial wavenumbers, if one assumes that not
all possible modes have been detected. To test this idea, the best matching
number of nodes for each detected mode-like feature has been looked for, up
to *n*<sub>r</sub> = 12. The results of this search are given in Table 3. It is legitimate
to ask whether the matches found are coincidental or not. In this respect
one should note that the consecutive distance between resonance radii of
successive number of radial nodes is ~ 100 km for *n*<sub>r</sub> = 5 and ~ 60 km for
*n*<sub>r</sub> = 10. The match in both instances would seem rather coincidental if the
dynamics producing these resonance locations were unrelated to the number
of radial nodes. Also, *m* = 1 modes are much easier to excite than all other
*m* > 1 modes, because *m* = 1 corresponds to the natural epicyclic motion
of test particles around an oblate planet, i.e., only weak collective effects are
needed to maintain the mode coherence, compared to the other modes (this
is also the reason behind the difference in cavity sizes between *m* = 1 and
<table><caption>Table 3: Comparison of the modes resonance location (as measured by the mode cavity size $δa = a_B - a_r$) for various radial node numbers $n_r$, and associated mode amplitude at the edge ($a_Bϵ_B$), determined from the *Cassini* data analysis (French et al., 2023) and from numerical solutions, for a constant surface density profile $σ_0 = 200$ g/cm$^2$ ($m = 2$ to 5) and $σ_0 = 60$ g/cm$^2$ ($m = 1$). Distances are rounded to the nearest kilometer.</caption><thead><tr><th>m</th><th>n<sub>r</sub></th><th>δa (km)<br>(data)</th><th>δa (km)<br>(theory)</th></tr></thead><tbody><tr><td>2</td><td>0</td><td>33</td><td>35</td></tr><tr><td rowspan="3">3</td><td>0</td><td>26</td><td>23</td></tr><tr><td>1</td><td>28</td><td>62</td></tr><tr><td>2</td><td>35</td><td>88</td></tr><tr><td rowspan="3">4</td><td>0</td><td>19</td><td>19</td></tr><tr><td>1</td><td>26</td><td>51</td></tr><tr><td>2</td><td>32</td><td>72</td></tr><tr><td rowspan="2">5</td><td>0</td><td>16</td><td>23</td></tr><tr><td>1</td><td>35</td><td>44</td></tr><tr><td rowspan="3">1</td><td>0</td><td>151</td><td>147</td></tr><tr><td>5</td><td>879</td><td>871</td></tr><tr><td>10</td><td>1208</td><td>1232</td></tr></tbody></table>

other modes). This is true for all such modes, independently of their radial number of nodes. This may explain why $m = 1$ may be the only azimuthal wavenumber for which free modes with different numbers of radial nodes may be found. However, if these modes are indeed associated with the different $n_r$, the absence of all the other $n_r$ values (i.e., 9 modes up to $n_r = 12$) is rather intriguing, but possibly not so unlikely as these modes may have too small amplitudes to be detected. On the intermediate time-scales of interest here, at least two processes may lead to varying mode amplitudes: overstable librations and nonlinear mode interactions. This may possibly make some of these modes undetectable at some times, but it still seems intriguing that so many of them would coincidentally have low amplitudes precisely at the time of the *Cassini* mission.

In any case, an alternative interpretation of the existence of this multiplicity of mode-like features must be found for $m = 2$ to 5. This is discussed
in section 7.3.

<figure><img src="image_9.png" /><figcaption>Figure 5: The $m = 3$ nodeless edge mode amplitude profiles $a\epsilon(a)$ (top) and distance between the resonances of the first two modes (bottom) as a function of the pseudo-collective nonlinearity parameter $q_c$; $\sigma_0 = 200$ g/cm$^2$; $q_{max}$ is slightly adjusted as $q_c$ increases in order to maintain the mode amplitude $a_B\epsilon_B$ at its observed value. The black dashed-lines represent physically motivated fits to the numerical results (see text).</figcaption></figure>
### 5.2.2. Semi-quantitative discussion of collective mode effects

As pointed out in section 4.1, the presence of a collection of modes at the B ring edge has some impact on the dynamics of individual modes. Indeed, one is also motivated to assume $H(q^2) = H(q_a^2) \simeq 1$ in Eq. (62) for the edge to be confined by Mimas; this requirement should hold over the cavity size of the forced $m = 2$ mode (see Longaretti 2018 for a discussion of this point), which is comparable for the other low $m$ nodeless modes cavity sizes. This results in halving the surface density needed to obtain matches between theoretical and observational results in Table 2. Indeed, solving Eq. (28) under the assumption $H(q^2) = H(q_a^2) \simeq 1$ leads to the same formal solution as for $H(q^2) = 1/2$, provided that the surface density is halved; for the same reason, $\Gamma \simeq 1/16$ is again the correct value to match the various modes amplitudes and Eq. (62) reduces to

$$ \sigma_0 = \frac{3(m-1)(a_B - a_r)^2 M_p}{8\pi^2 \Gamma a_r^4}. \qquad (63) $$

Investigating the collective effects on mode spacing presented in Table 3 is however less straightforward. A very simple attempt has been made in this direction by exploiting in a different form the idea that $q$ must be of the order of $q_a$. In practice, some constant $q_c^2$ contribution to $q^2$ has been added in the self-gravity term, modifying the self-gravity nonlinear function $H(q^2)$ into $H(q_c^2 + q^2)$ with various values of $q_c$. This rather crude expedient is suggested by the short discussion of the collective effect of multiple mode in section 4.1.

It turns out that the spacing of the various successive $n_r$ modes resonance location for a given $m$ increases with $q_c$, as shown on the bottom graph of Fig. 5. This behavior can be explained in the following way. First note that the cavity size of the nodeless mode $\delta a$ (top graph of the same figure) is a fraction of the first wavelength of the underlying standing density wave, and the same remark applies to the distance $\Delta r_s$ between the resonance location of the first two modes (which should be comparable to the second half of the first wavelength). This argument suggests that $\Delta r_s \propto \delta a$. Next, the scaling arguments used in interpreting Fig. 4 applies *mutatis mutandis*, i.e., Eq. (60) applies with $\tilde{q}^2 = q_c^2 + q_{max}^2/2$. This argument implies that $\delta a \propto \Delta r_s \propto [H(\tilde{q}^2)]^{1/2}$. Thus $\delta a = \delta a^*[H(\tilde{q}^2)/H(\tilde{q}^{*2})]^{1/2}$ and $\Delta r_s = \Delta r_s^*[H(\tilde{q}^2)/H(\tilde{q}^{*2})]^{1/2}$ with the same meaning as before for the starred quantities (now chosen for $q_c = 0$). The resulting behavior is displayed by the dashed lines on Fig. 5.
This makes the possibility that the modes associated with a given azimuthal wavenumber $m$ correspond to modes of higher number radial of nodes even less realistic, but such a conclusion needs to be confirmed by true multimode dynamical analyses.

## 6. Numerical solution: forced mode

The simulation domain is extended to $4\lambda_{dw}$ compared to the $2\lambda_{dw}$ choice used for most free mode numerical solutions, for unessential technical reasons. Only the last $3\lambda_{dw}$ are shown in the figures, to focus on the informative part. The number of streamlines is now fixed at $N = 400$, to take into account the increase in size of the computational domain without losing precision in the computations.

### 6.1. Numerical procedure

For the Mimas 2:1 resonance, the problem is solved in two steps:

1. The linear solution ($H = 1/2$) of Eq. (31) for the vector of eccentricity values along the mode is first determined through matrix inversion:

$$ \epsilon = (L - m\Omega_p I/\bar{n})^{-1} F \cdot \mathbf{1} \qquad (64) $$

where $\mathbf{1}$ is the unit vector and $I$ the identity matrix, $\Omega_p$ the pattern speed forced by the satellite resonance and $\bar{n}$ some characteristic value of the elliptic mean motion, e.g., $n(a_r)$ or $n(a_B)$ (the difference between the two is negligible). The (assumed constant) surface density is adjusted iteratively until the mode amplitude at the edge corresponds to the observations within better than 1 km.

2. The nonlinear solution is then found iteratively with the linear solution as a guess, until the solution is converged with respect to a relative change in the eccentricity profile of the order of $10^{-6}$ (see section 5.1 for this part of the procedure).

### 6.2. Numerical results

The quantities of interest in this section are often normalized by the equivalent quantity in the linear theory of density waves. Also the $m = 2$ forced mode is compared to the equivalent linear density wave that would propagate in the absence of the edge. The relevant expressions in the streamline formalism framework can be found in section 7.2 of Longaretti (1992), in particular Eqs. (7.43), (7.50) and (7.53).
<figure><img src="image_10.png" /><figcaption>Figure 6: The linear edge mode profile $\epsilon(a)$ (blue) normalized by the linear density wave asymptotic amplitude $\epsilon_0 = 2.5 \times 10^{-4}$ (left axis); $a_B\epsilon(a)$ in km (right axis); $\lambda_{dw} = 139.5$ km and $\sigma_0 = 220$ g/cm$^2$. The origin of radii is the resonance location, indicated by a green vertical dotted line, as well as the ring edge. Also shown are the inviscid density wave solution profiles $\epsilon \cos m\Delta^{dw}$ (black), $\epsilon \sin m\Delta^{dw}$ (red) and $\epsilon$ (dotted black). The crosses represent the test particle solution. The vertical black dotted line specifies the limit of validity of the test particle approximation [i.e., a distance $a_r(M_s/M_p)^{1/2}$ inside the resonance location].</figcaption></figure>

### 6.2.1. Eccentricity profile of the B ring forced edge mode

It is informative to contrast first an edge mode and a density wave. To this effect, Eq. (31) is solved with a constant surface density profile ($\sigma_0 = 220$ g/cm$^2$ assumed; see below the origin of this value). Fig. 6 compares the behavior of the linear version of Eq. (31) (i.e., where $H(q^2) = 1/2$ has been set) with the behavior of $\epsilon \cos m\Delta^{dw}$, $\epsilon \sin m\Delta^{dw}$ and $\epsilon$ for the underlying forced outgoing density wave in the linear regime. The test particle solution is also shown on this figure. All amplitudes are normalized by the density wave asymptotic amplitude $\epsilon_0 = 2.5 \times 10^{-4}$, and the radial coordinate is
<figure><img src="image_11.png" /><figcaption>Figure 7: Top: The linear (red) and nonlinear (blue) edge mode profiles $\epsilon(a)$; bottom: associated $q$ profiles; $\epsilon_0 = 2.5 \times 10^{-4}$, $\lambda_{dw} = 139.5$ km and $\sigma_0 = 220$ g/cm$^2$ (see text). The origin of radii is the resonance location. The crosses represent the test particle solution. The vertical black dotted line specifies the limit of validity of the test particle approximation.</figcaption></figure>

$(a - a_r)/\lambda_{dw}$ where $\lambda_{dw} = 139.5$ km is the forced linear density wave first wavelength. The linear theory of density waves underestimates the amplitude at the edge by a factor $\simeq 2$. Both the linear density wave and the test particle $\epsilon$ profiles provide reasonably good approximations to the linear edge mode
profile in the evanescent region. This is related to the fact that collective
effects become negligible in the evanescent region. Fig. 7 compares the forced
linear and nonlinear $\epsilon$ and $q$ profiles for the same setting. The nonlinear
correction is rather small as the nonlinearity parameter $q$ never reaches very
large values.

The value of $\sigma_0$ has been obtained from the insight gained through Eq. (50). Neglecting the $\omega(a_2)$ term in this equation (which is only a fraction of the other term) yields a relation of the form $m_r\epsilon_B = \text{cst}$ with $m_r = 4\pi\sigma_0 a \delta a$. This approximate relation indicates that $m_r = m_r(\epsilon_B)$ with a monotonically decreasing behavior: when the perturbed region mass $m_r$ is too small, the ring edge amplitude $\epsilon_B$ is too large with respect to the observed amplitude, and conversely. Reverting to the exact equation, this constraint is reexpressed in an iterative form $m_r^{(i)} = m_r^{(i-1)} \times (\epsilon_e^{(i)}/\epsilon_B)$ where $m_r^{(i-1)}$ is the mass of the edge perturbed region at iteration $(i-1)$ and $\epsilon_e^{(i)}$ the computed edge mode amplitude at iteration $(i)$ from this assumed mass. The mass of the perturbed region is adjusted through this relation until the computed value of $a\epsilon_e^{(i)}$ matches the observed one (a tolerance of 1 km has been adopted as stopping criterion). The iteration is started with Eq. (51) providing an initial guess for the surface density. Note that this process will converge as long as $m_r$ and $\epsilon_B$ vary in opposite directions and the functional dependence $m_r = m_r(\epsilon_B)$ is sufficiently well-behaved [i.e., $m_r(\epsilon_B)$ sufficiently flat] which is the case here. The value $\sigma_0 = 220 \text{ g/cm}^2$ has been obtained for the nonlinear solution, but as the difference between the linear and nonlinear solutions is small, it is nearly valid also for the linear one.

### 6.2.2. Apsidal shift profile of the B ring forced edge mode

Let us now turn to the solution of Eq. (36). The resulting profile of $\epsilon \sin m\Delta$ is shown in Fig. 8.

In this computations, the magnitude of $\sin m\Delta$ at the edge is adjusted to match $\Delta_B \sim 2.9^\circ$ (the largest of the values obtained for this quantity from the various analyses of the *Cassini* mission data). An iterative procedure has been implemented to this effect. This procedure is the same as in section 6.2.1, *mutatis mutandis*: the effective viscosity $\nu_{eff}$ is proportional to the apsidal shift at the edge $\Delta$; it is modified until convergence to the observed apsidal shift is obtained within a relative tolerance of $10^{-2}$. As Eq. (36) is linear in $\epsilon \sin m\Delta$ [$q_{ij}$ and $q$ are fixed by the solution of Eq. (31)], this requires only one iteration. For the same reason, the figure is easily rescaled for a different choice of $\Delta_B$. The maximum change in the profile of $q$ due to $m\Delta$
<figure><img src="image_12.png" /><figcaption>Figure 8: The edge mode $\epsilon \sin m\Delta$ profiles for the stress tensor model of section 3.3 and a constant surface density profile (blue) and a linearly increasing one (red). The profiles are normalized by the linear density wave asymptotic amplitude ($\epsilon_0 = 2.5 \times 10^{-4}$); $\lambda_{dw} = 139.5$ km and $\sigma_0 = 220$ g/cm$^2$. The origin of radii is the resonance location.</figcaption></figure>

does not exceed 1%, self-consistently justifying the ordering between the ring self-gravity and stress tensor assumed in section 3.

The effective viscosities obtained from this procedure are given by

$$ \nu_{eff} = 3 \times 10^4 \left( \frac{\Delta_B}{2.9^\circ} \right) \text{cm}^2/\text{s}, \qquad (65) $$

$$ \nu_0 \sim 15 \left( \frac{\Delta_B}{2.9^\circ} \right) \text{cm}^2/\text{s}, \qquad (66) $$

where Eq. (54) has been used in the second expression, relevant for the viscous angular momentum transport in the unperturbed region.
### 6.2.3. Torque and wake transport

We are now in position to perform a numerical integration of the torque density and scale the result against the Borderies et al. (1982) approximation $T_{BGT}$ of Eq. (56). For this scaling the torque is evaluated for the constant surface density $\sigma_0 = 220 \text{ g/cm}^2$ (see section 6.2.1). For the B ring edge parameters, this leads to a reference torque value of $T_{BGT} = -1.7 \times 10^{20} (\Delta_B/2.9^\circ) \text{g.cm}^2/\text{s}^2$.

The satellite torque due to the forced mode is computed from the numerical solution with the help of its formal exact integral expression:

$$ 
\begin{aligned}
T_s &= \int da \, \pi m a \sigma_0 \Psi_{mk} \epsilon \sin m \Delta \\
&= 1.4 T_{BGT} \left( \frac{\Delta_B}{2.9^\circ} \right) = -2.4 \times 10^{20} \left( \frac{\Delta_B}{2.9^\circ} \right) \text{g.cm}^2/\text{s}^2.
\end{aligned} 
\quad (67) $$

The BGT approximation of the torque appears to be surprisingly precise, considering the number of factors of order unity or more ignored in its derivation.

Balancing the torque $T_s$ with the unperturbed viscous angular momentum flux $L_h^{vis} = 3\pi\sigma_0 na^2\nu_0$ gives another, more precise determination of $\nu_0$

$$ \nu_0 = \frac{|T_s|}{3\pi\sigma_0 na^2} \approx 15 \left( \frac{\Delta_B}{2.9^\circ} \right) \text{cm}^2/\text{s}, \quad (68) $$

where the estimate $\sigma_0 \approx 60 \text{ g/cm}^2$, obtained in section 5.2 from the analysis of the free $m = 1$ nodeless mode, has been used as characterizing the unperturbed region. This is in remarkable agreement with the cruder estimate of Eq. (66).

On the other hand, for this value of $\sigma_0$, Eq. (58) with $C \approx 7 - 13$ yields $\nu_{sg} \approx 13 - 27 \text{ cm}^2/\text{s}$, which is consistent with $\nu_0$, but not comfortably so as $\Delta_B$ is scaled by what is likely a maximum value. The overlap may be improved in several different ways. First, $\Delta_B$ may be affected by librations and this may in turn provide biased measurements for the long term time-average of this quantity that is of interest here. Second, the estimate of the torque may be too low. It turns out that increasing linearly the ring surface density by a factor of 3 in the region of extent $\lambda_{dw}$ inside the edge increases the torque by a factor of 2, i.e., covering the range of uncertainty of $\nu_{sg}$. Such alternative surface density profiles have been somewhat explored numerically for the discussion of the present paragraph, but not in an extensive way. It
remains to be seen if the related enhancement in surface density is compat-
ible with the actual optical depth profile close to the edge (which is as yet
unknown, due to the complex changes in compression with time due to the
multiplicity of modes observed there). Conversely, changes in the stress ten-
sor model (e.g., changing $q_1$ from -1 to 0, or adopting a granular flow model
for the stress tensor instead) do not change the torque by more than a few
percents.

This short discussion suggests that the results of this analysis are compat-
ible with the dominance of wake transport in the unperturbed B ring region,
a conclusion whose validity should however await a detailed analysis of the
edge time-dependent dynamics to be more reliably assessed.

7. Summary and discussion

In this work, the structure and properties of edge modes at the B ring
edge have been quantified following the theoretical formulation of the trapped
density wave picture of edge mode briefly outlined in Longaretti (2018),
following an unpublished lead by Peter Goldreich. This analysis relies on
the streamline formalism (Borderies et al., 1983, 1985, 1986), recast in the
fluid context in epicyclic variables (Longaretti and Borderies, 1986, 1991;
Borderies-Rappaport and Longaretti, 1994; Longaretti, 1992, 2018).

One expects that the rather large collection of free edge modes identified
in this work among the mode-like features found by Nicholson et al. (2014)
and French et al. (2023) have some mutual and non-resonant influence on
their dynamics. However, as a first approximation, it is both informative
and useful to explore this dynamics while ignoring these interactions. The
kinematics and dynamics of the trapped wave picture has been more thor-
oughly formulated here in this independent mode limit, as well as the related
formal description of edge modes as solutions of an eigenvalue problem in the
linear limit; in this view, a collection of free modes can be associated with a
given azimuthal wavenumber *m*, each corresponding to a different number of
radial nodes *n*<sub>*r*</sub>. The issue of the correct boundary conditions to be applied
to such modes in numerical solutions has also been addressed in Appendix
C.

The resulting stationary dynamics of the B ring edge modes has first been
explored through a semi-quantitative analysis relying on a two-streamline ap-
proximation which provides insight into the physical origin of the connection
between mode structural properties (resonance location, edge amplitude, etc)
and important ring physical properties, most notably the magnitude of the surface density in the edge mode region and the level of energy dissipation through the ring effective viscosity. These structural properties and ring diagnostics have been quantified next more precisely through numerical solutions of the mode structure equations.

### 7.1. Stationary free single modes

Free mode-like features have been detected for a variety of azimuthal wavenumbers, namely $m = 1$ (three features), $m = 2$ (one), $m = 3$ (three), $m = 4$ (three), and $m = 5$ (two). Whether all these features are actual modes in the meaning specified by the stationary mode equations, or some other type of time-dependent feature associated with a given $m$ mode, is not a priori obvious. In this respect, an important issue is to assess whether the various modes detected for each azimuthal wavenumber can be associated with edge modes of different radial wavenumbers. For $m \ge 2$, the present investigation concludes that only the features with associated resonance location closest to the edge which can be positively interpreted as a mode for the various azimuthal wavenumbers (to the exception of $m = 5$); in all cases, this is a nodeless mode. This conclusion follows from the fact that the observed resonance location spacing of mode-like features is systematically much smaller than the theoretical prediction for modes of successive number of radial nodes, for $m \ne 1$. For $m = 1$, on the contrary, all mode-like features can be interpreted as actual free modes and many such modes are missing in the observations.

The surface density in the edge region can be estimated from free mode resonance locations. This is applied to the nodeless modes identified for each $m$. In the independent mode dynamics limit, this gives $\sigma_0 \simeq 200$ g/cm$^2$ for the nodeless $m = 2$ to $m = 4$ modes; this does not apply to $m = 5$, as the nodeless solution seems to be absent from the data. Collective mode effects are expected to bring this down to $\sim 100$ g/cm$^2$ for $m = 2$ to $m = 4$ (see section 4.1). The $m = 1$ nodeless mode, which extends substantially deeper in the B ring than $m \ne 1$ modes, gives a more modest surface density, $\sim 60$ g/cm$^2$, which indicates that the edge itself (i.e., in the last 20 or 30 km to the B ring edge) is denser than the region inside it.

Besides the intrinsic uncertainty on the magnitude of the collective effect correction just mentioned, the next major source of uncertainty comes from the assumption of a constant surface density profile. Preliminary tests with
non-constant profiles suggest that the associated uncertainty is ~30% on the determination of the surface density.

Another source of uncertainty comes from the identification of the nodeless mode for each *m* with the feature whose associated resonance is closest to the ring edge. This is necessarily the case in the stationary dynamics limit of independent modes. It will be shown elsewhere that viscous overstabilities do not change this conclusion. However, mode non-resonant interactions through the ring self-gravity may lead to the production of mode-like features that may kinematically resemble an additional mode closer to the edge, if one focuses only on the edge streamline, as in data analyses. This possibility as well as its implications for the present conclusions will be discussed elsewhere as well.

## 7.2. Stationary forced single mode

The *m* = 2 mode forced by Mimas has an externally determined resonance location. However, for this mode, and contrarily to free modes, the amplitude is forced as well. This amplitude can again be related to the ring surface density through the action of the ring self-gravity and planet differential rotation, leading to $\sigma_0 \sim 220 \text{ g/cm}^2$ for this mode in the independent mode dynamics limit. Again, a reduction by a factor of $\sim 2$ is expected under the collective effect of edge modes, bringing this estimate in the $\sim 100 \text{ g/cm}^2$ vicinity; these estimates are consistent with the results of the analysis of the free modes summarized above.

In the dissipationless limit, the forced $m = 2$ mode periapses should be aligned with Mimas. However, a small shift of the order of 1 to 3° has been detected (Nicholson et al., 2014; French et al., 2023; Spitale and Porco, 2010), and is related to the level of viscous dissipation in the ring edge perturbed (last 20 to 30 km) and in the unperturbed region (further inside the ring). The sign of this shift is consistent with the ring outer edge being in a normal, viscous damping regime, in agreement with both the requirement of angular momentum transfer from the ring to Mimas due to the torque exerted by this satellite through the forced mode and the requirement of viscous angular momentum flux being nearly cancelled close to the edge for this confinement to be possible at all (Borderies et al., 1984, 1989; Longaretti, 2018). This last point implies in turn that potential viscous overstabilities should be quenched (Borderies et al., 1985; Longaretti, 1992), at least in the outermost part of the perturbed edge region and in the stationary dynamics limit. A way around
this conclusion in a time-dependent dynamics context is proposed in the final discussion, section 7.3.

The magnitude of the associated effective viscosity close to the edge is quite large and uncertain within an order of magnitude due to various modeling uncertainties (in particular the shape of the surface density profile and the magnitude of the observed shift between Mimas and the forced edge mode), $\nu \sim 1 - 6 \times 10^4$ cm$^2$/s. This is however consistent with, on the one hand, the large increase in dissipation in the edge perturbed region required by the balance of energy and angular momentum transfers from the ring to the satellite (Borderies et al., 1982) and with a related magnitude of the effective viscosity in the unperturbed region of $\nu \sim 5 - 30$ cm$^2$/s. This last upper limit in particular is roughly consistent with the fact that angular momentum transport is dominated by self-gravity wakes in unperturbed regions; indeed, for a surface density of the order of 60 g/cm$^2$ deduced from the analysis of the $m=1$ mode (this should be characteristic of the unperturbed region, considering the rather small level of nonlinearity of this mode and deeper penetration of this mode in the ring compared to $m \neq 1$ modes) the wake transport viscosity is estimated at $\nu_{sg} \simeq 13 - 27$ cm$^2$/s. This analysis has assumed that self-gravity wakes are quenched in the outer, perturbed region of the B ring edge. This possibility is supported by the very large excess dissipation just mentioned and the related expected increase in the ring velocity dispersion, which should prevent Toomre's criterion of instability of self-gravity wakes from being satisfied (Toomre, 1964).

The $q$ profile presented above for the forced mode (see Fig. 7) never reaches values as large as the ones quoted in Hahn et al. (2009) (see in particular their figures 7b and 9b). The explanation for this behavior is related to the behavior of the self-gravity integral$^{15}$, which tends, as a rule, to compensate changes of $\sigma_0$ by changes of $q$. The model of these authors produces a drop in the surface density close to the edge, with a nearly constant surface density inside the last 50 km or so, whereas observational constraints exclude the existence of a surface density drop until the last kilometer or last few kilometers to the edge, as assumed here. Consequently, the maximum $q$ in Hahn et al. (2009) is necessarily larger than the one obtained here.

Quite clearly, the modest $q$ values ($0.5$ at most) obtained here seem at

<footnote>$^{15}$There is no formal proof of this behavior, but it is easily noticed in numerical solutions and derives from the fact that both quantities appear in the combination $\sigma_0 H(q_{ij}^2)$.</footnote>
odd with the requirement of viscous angular momentum flux reversal ($q \simeq q_a \simeq 0.8$ for a variety of stress tensor models). In spite of their surface density profile, Hahn et al. (2009) do not reach this requirement either for their stress tensor model. The failure to meet this requirement is at the root of the failure of their analysis of the B ring edge confinement by Mimas for realistic stress tensor models, and also raises a similar issue for the original confinement process of Borderies et al. (1982).

The only option at this point for this confinement to occur under realistic physical conditions is through the collective edge mode effects already mentioned. Note in particular from Table 2 that $(\sum_{\alpha} q_{\alpha, max}^2)^{1/2} \sim 0.85$ where $\alpha$ is an index identifying the various nodeless modes. Such an expression of the effective $q$ of a collection of modes will be justified in detail elsewhere, but can also be deduced from the analysis of Longaretti (1989). It provides a time average estimate of a number of modes’ collective $q$, and is close to the expected value $q_a \simeq 0.8$ considering the approximate nature of the argument. Note that although only the nodeless mode contributions are taken into account in this sum, the other mode-like features have modest amplitudes and should have negligible contributions to this sum; if these come out of nonlinear non-resonant couplings, they should average out on intermediate time-scales anyway. Although the independent mode dynamics analysis performed in this work is a first approximation, these remarks are encouraging; a more complete explanation of the B ring edge confinement dynamics along these lines will be explored in more detail elsewhere.

### 7.3. Perspectives

The most important caveat in this investigation relates to mode mutual interactions, and possibly as well, the existence of overstable librations of edge modes. Both processes are mostly mediated by the ring self-gravity. This issue is addressed in a forthcoming paper, and only a few important points will be mentioned here, which should be considered as indicative but provisional until this analysis is published. They are included here because they help forming a self-consistent picture of the B ring edge dynamics.

A first, rather formal question, concerns the possibility of representing each edge mode in the usual, single mode approximation adopted here, in spite of the presence of other modes. As will be shown elsewhere, it turns out that the dynamics of each mode is formally very similar to a single mode dynamics, except for a formal redefinition of effective quantities, such as $H(q^2)$ for the self-gravity, and $t_1$, $t_2$ and $a_{r\theta}$ for the stress tensor. The formal change
is rather mild — substitution of a general, long-period dependent $Q$ parameter to the mode $q$ in all these quantities, whose time average expression has been used in the forced mode discussion right above. The kinematic justification of this result can be found in Longaretti (1989), as well as preliminary elements of the dynamics.

The most important consequence of this substitution of the generalized nonlinearity parameter $Q$ to the usual $q$ is that the amplitudes of all modes (eccentricities in particular) undergo slow quasi-periodic variations, on time scales of $\sim 10$ yrs. This in itself can explain why the nodeless modes for some $m$ may not be detected at some time (in particular for the B ring edge, $m = 5$). Again, some partial justification of this result can be found in Longaretti (1989).

Another important consequence is that the edge confinement requirement $Q = q_a$ cannot be satisfied at all times, but only in some time-average sense. If at the same time the damping coefficient $t_1$ is negative on average, all free modes and librations should eventually be damped, which suggests that, on the contrary, some process is at work preventing this. On the one hand, the simple model of Eqs. (37) to (40) does indeed imply that if $a_{r\theta}$ vanishes on average (therefore ensuring that the confinement of the B ring edge is satisfied), $t_1 < 0$ on average. On the other hand, this last constraint cannot apply throughout the entire edge region, as $q$ is bound to vanish far enough from the edge. In the inner region, $t_1 > 0$ on average is expected if $q_1 > 0$ ($q_1$ being the characteristic value of $q$ for which $t_1$ changes sign). Finally, a viscous overstable inner edge region will produce some libration in the outer region due to coupling through the ring self-gravity, with an amplitude defined by a balance between self-gravity forcing and the outer edge region unavoidable damping as $t_1 < 0$ applies there. The whole issue therefore revolves around the expected amplitude of this forced libration, and correlatively, on the relative size of the inner overstable region with respect to the stable outer one.

In any case, these time-dependent phenomena may offer a natural explanation for the prolonged existence of edge modes. They also provide a potential explanation as to why some modes are detected and not others: if partial viscous overstabilities do indeed take place, they may all be overstable, but due to the time-dependence of their amplitudes (itself a combined effect of viscous overstabilities and mode long-period interactions), some may not be detectable all the time (see Longaretti and Rappaport 1995 for some justification of this statement and of the qualitative discussion of this para-
graph). This long period time dependence also offers a dynamical cause for the other mode-like features associated with the various azimuthal wavenumbers detected at the B ring edge. In this perspective, one can show that the features associated with the nodeless mode should be identified with the time-independent part of overstable modes, which does indeed conform to the stationary dynamics explored here, except for the long-period variations due to mode interactions, which transpire mostly through varying eccentricities $\epsilon$ and associated edge amplitudes $a_B\epsilon_B$. This justifies or at least motivates the relevance of the approximation studied in the present work, in the more general context of multiple edge mode coupled dynamics. Also, this rather complex dynamics is most likely the source of the yet unexplained residuals obtained in the mode fitting procedure of Nicholson et al. (2014) and French et al. (2023). Improving on this procedure may be quite challenging though.

Finally, the long-period mode interaction dynamics briefly described here makes the edge confinement issue a time-dependent problem — and a considerably more complex one to analyze than previously suspected.

This paper is devoted to the B ring edge modes and mode-like features. Similar, and much more numerous features are detected at the A ring outer edge (Nicholson et al., ????). The dynamics there is substantially more complex due to the fact that the forcing by Janus is intermittent as a consequence of the Janus-Epimetheus shared horseshoe orbit. Consequently, and much more than for the B ring edge, a time-dependent dynamical analysis will be required to provide insight into the peculiar features observed at the A ring outer edge.

## Appendix A. Simple analytic models of the ring stress tensor

This Appendix expands upon the short discussion of the stress tensor in section 3.3. To this effect a simple hydrodynamic-like approximation is discussed, with appropriate simplified choices of the pressure and viscosity model. Namely,

$$P_{ik} = p\delta_{ik} - 2\sigma\nu u_{ik}, \qquad (A.1)$$

where $u_{ik}$ is the strain tensor defined from the fluid velocity $\mathbf{u}$

$$u_{ik} = \frac{1}{2} \left( \frac{\partial u_i}{\partial x_k} + \frac{\partial u_k}{\partial x_i} \right), \qquad (A.2)$$
and where<fnref n="16" /> $p = \sigma c^2$ and $\nu = \alpha(c^2/n)(\sigma/\sigma_0)^\beta$ ($c$ is the ring velocity dispersion, but can also be interpreted as the ring sound speed). In this model, $\alpha$ is a constant of order unity, and $\beta = 0$ for a simple viscously stable model, while $\beta \gtrsim 1$ allows for viscous overstabilities to develop; $\beta \simeq 1$ from N-body simulations (Salo et al., 2001). Such a model is often used in the analysis of axisymmetric viscous overstabilities (e.g., Schmidt et al. 2001; Latter and Ogilvie 2009), which are now believed to form most of the small-scale structure of the A and B rings. This model has the same dependence on dimensional quantities for $\beta \ge 1$ as the model used in section 3.4 and a similar dependence on $q$. This is explicitly shown below for $\beta = 1$; in particular the two most important physical properties of this $q$ dependence are preserved ($q_1 < q_a$ and divergence when $q \to 1$).

In a cylindrical coordinate system centered on Saturn, one has

$$u_{rr} = \frac{\partial u_r}{\partial r}, \qquad (A.3)$$

$$u_{r\theta} = \frac{1}{2} \left( \frac{1}{r} \frac{\partial u_r}{\partial \theta} + \frac{\partial u_\theta}{\partial r} - \frac{u_\theta}{r} \right), \qquad (A.4)$$

$$u_{\theta\theta} = \left( \frac{1}{r} \frac{\partial u_\theta}{\partial \theta} + u_r \right). \qquad (A.5)$$

Only $P_{rr}$ and $P_{r\theta}$ are necessary in the dynamical equations of section 3.3. Combining Eqs. (A.3) and (A.4) with Eqs. (10) and (11), one finds

$$P_{rr} = \frac{\sigma_0 c^2}{J} - 2\nu_{eff}\sigma_0 n \frac{q \sin M'}{J^{\beta+2}}, \qquad (A.6)$$

$$P_{r\theta} = -2\nu_{eff}\sigma_0 n \left[ \frac{1}{4J^{\beta+2}} - \frac{1}{J^{\beta+1}} \right], \qquad (A.7)$$

to leading order in $\epsilon$, $J_2$ and $\delta a/a$. In this expression, $M' = m(\varphi - \Omega_p t) + m\Delta + \gamma$ and $\nu_{eff} = \alpha c^2/n$ have been defined.

More precisely, we only need to evaluate $t_1$, $t_2$ and $a_{r\theta}$, defined by:

$$t_1 = s_{rr} + 2c_{r\theta}, \qquad (A.8)$$

$$t_2 = 2s_{r\theta} - c_{rr}, \qquad (A.9)$$

$$c_{\alpha\beta} + i s_{\alpha\beta} = \langle \exp(iM') P_{\alpha\beta}(M') \rangle, \qquad (A.10)$$

$$a_{r\theta} = \langle P_{r\theta}(M') \rangle, \qquad (A.11)$$

<footnote n="16">¹⁶A vertical integration has been performed under an assumption of vertical hydrostatic equilibrium to obtain these relations.</footnote>
where the bracket notation stands for the azimuthal average over $M' = m(\varphi - \Omega_p t) + m\Delta + \gamma$.

For $\beta = 0$ (hydrodynamic model with kinematic viscosity independent of the ring surface density) and $\beta = 1$ (viscous overstability model), $t_1$, $t_2$ and $a_{r\theta}$ have analytic expressions.

For $\beta = 0$:

$$t_2 = -\sigma_0 c^2 q H(q^2), \qquad (A.12)$$

$$t_1 = -\sigma_0 \nu_{eff} n \left[ \frac{3q^2 - 2 + 2(1 - q^2)^{3/2}}{2q(1 - q^2)^{3/2}} \right]$$

$$\simeq -\sigma_0 \nu_{eff} n \left[ \frac{3q^3}{4(1 - q^2)^{3/2}} \right], \qquad (A.13)$$

$$a_{r\theta} = \sigma_0 \nu_{eff} n \left[ \frac{3 - 4q^2}{2(1 - q^2)^{3/2}} \right]. \qquad (A.14)$$

The second expression for $t_1$ is accurate to better than 25% in the whole range $0 \le q \le 1$. For this model, both $t_1$ and $t_2$ are negative. Viscous angular momentum flux reversal takes place when $a_{r\theta} < 0$ i.e., $q > q_a = \sqrt{3}/2 \simeq 0.87$ (Borderies et al., 1982).

For $\beta = 1$:

$$t_2 = -\sigma_0 c^2 q H(q^2), \qquad (A.15)$$

$$t_1 = -\sigma_0 \nu_{eff} n \left[ \frac{q(2q^2 - 1)}{2(1 - q^2)^{5/2}} \right], \qquad (A.16)$$

$$a_{r\theta} = \sigma_0 \nu_{eff} n \left[ \frac{6 - 7q^2}{4(1 - q^2)^{5/2}} \right]. \qquad (A.17)$$

The expression for $t_2$ does not depend on $\beta$ and is unchanged. Angular momentum flux reversal now occurs for a substantially higher value of $q_a \simeq (6/7)^{1/2} \simeq 0.93$. It is unclear that this model is realistic in this respect, as such high values of $q_a$ are not found in microphysically motivated models or in numerical simulations. The sign of the viscous term $t_1$ now depends on $q$. For $q < q_1 = 1/2^{1/2} \simeq 0.71$, $t_1 > 0$, leading to the excitation of the complex eccentricity $Z$ instead of its damping. For $q > q_1$, the usual damping behavior is recovered.
## Appendix B. Satellite resonant forcing

Consider a satellite with mass $M_s$ whose orbit is characterized by its epicyclic semimajor axis $a_s$ and epicyclic eccentricity $\epsilon_s$. The relevant component of the satellite potential giving rise to a Lindblad resonance reads (Goldreich and Tremaine, 1980, 1982; Shu, 1984)

$$ \phi_s(r, \theta, t) = \Phi_{m_s, k_s}(r/a_s) \cos m_s(\theta - \Omega_{p,s}t). \quad (B.1) $$

This relation follows from a double Fourier expansion of the satellite exact potential: one in azimuth (sum over $m_s$ terms) and one with respect to small deviations from purely circular motion (sum over $k_s$ terms). Consequently the satellite potential Fourier component $\Phi_{m_s, k_s} \propto \epsilon_s^{|k_s|}$ and the strongest resonances correspond to $k_s = 0$. For such resonances (the only ones of interest in this work), $\Omega_{p,s} = \Omega_s$ (the satellite epicyclic angular velocity).

The satellite resonant forcing is characterized by (setting $r = a$ in $\Phi_{m_s, k_s}$)

$$ \Psi_{m_s, k_s}(a) \equiv a \frac{d\Phi_{m_s, k_s}}{da} + 2m_s \Phi_{m_s, k_s}. \quad (B.2) $$

By convention, for an inner Lindblad resonance (ILR), $m > 0$. One has, for $m > 1$ (Goldreich and Tremaine, 1980; Shu, 1984) (the only case of interest here)

$$ \Phi_{m_s, k=0}(r/a_s) = -\frac{GM_s}{a_s} b_{1/2}^{m_s}(r/a_s). \quad (B.3) $$

In this relation, $b_{1/2}^j$ is the usual Laplace coefficient. Its derivative is usually computed with the help of the following identity

$$ \frac{db_{\gamma}^{m_s}(\beta)}{d\beta} = \gamma [-2\beta b_{\gamma+1}^{m_s}(\beta) + b_{\gamma+1}^{m_s+1}(\beta) + b_{\gamma+1}^{m_s-1}(\beta)]. \quad (B.4) $$

$\beta = r/a_s < 1$ and $m_s > 1$ at an inner Lindblad resonance (ILR); $\beta > 1$ and $m_s < 0$ at an outer Lindblad resonance (OLR).

Let us introduce an auxiliary dimensionless quantity $H_{m_s, k_s}$ such that

$$ \Psi_{m_s, k_s} = -\frac{GM_s}{a_s} H_{m_s k_s}(\beta), \quad (B.5) $$

where $\beta = a_r/a_s$ for an ILR$^{17}$. For a leading order ILR ($k_s = 0$), one has

<footnote>$^{17}$The expression for an OLR is obtained by switching $a_r$ and $a_s$, with $m_s < 0$ in Eq. (B.2).</footnote>
$$H_{m_s,0} = \left( 2m_s + \beta \frac{d}{d\beta} \right) b_{1/2}^m(\beta). \qquad (B.6)$$

In the large $m$ limit, this reduces to

$$H_{m_s,0} = \frac{2m_s}{\pi} \left[ 2K_0(2/3) + K_1(2/3) \right] \approx \frac{5m_s}{3}, \qquad (B.7)$$

where $K_0$ and $K_1$ are modified Bessel functions (Goldreich and Tremaine, 1980). The last approximate expression is surprisingly accurate even for low $m_s$. For example, in the case of the Mimas 2:1 resonance, $5m_s/3 = 3.3$ and $H_{m_s,k_s} = 2.38$. This approximation is therefore quite sufficient for quick estimates, for all $m > 1$ values.

## Appendix C. Numerical solution: boundary conditions

Detailed mode structures are obtained numerically only for stationary single mode dynamics in this work. Detailed time-dependent phenomena (librations, multiple mode interactions and evolution) require a much more sophisticated numerical approach and will be explored elsewhere.

The numerical solution of Eqs. (30) and (31) requires us to specify boundary conditions. The two boundaries of an edge mode are of a very different nature. The outer boundary is sharp and free, i.e., no stress is applied on the ring from the outside, which is essentially empty. The inner one is by design more fuzzy, as there is no sharp radius $a_i$ where $q = 0$ for $a < a_i$ and $q > 0$ for $a > a_i$. In fact, $q \to 0$ as $a \to 0$ in a rather smooth way, although in practice $q$ becomes negligible a few wavelengths of the trapped wave inside the resonance location.

Therefore the outer boundary discussion is a physical one, whereas the inner boundary condition is a numerical one. Also, these two discussions are somewhat different for the ring self-gravity and the ring stress tensor. These considerations underlie the structure of the present section.

Only single modes are discussed in this section. However, the results established here apply to a collection of edge modes as well, *mutatis mutandis*.

### Appendix C.1. Outer edge

#### Appendix C.1.1. Ring stress tensor

The stress tensor components are constant at the boundary as the fluid is bounded by a vacuum (i.e., stress-free) and the stress must be continuous for
the acceleration per unit volume to remain finite. A necessary consequence of this vanishing is that $\Delta P_{ik} = P_{ik}(\mathbf{r} + \boldsymbol{\xi}) - P_{ik}(\mathbf{r}) = 0$ at the boundary, where $\boldsymbol{\xi}$ is the fluid particle displacement from its unperturbed location (the stress tensor must remain constant during the displacement of a fluid particle at the free surface). This fluid displacement occurs at constant $a$, as fluid particles move on streamlines of constant $a$, at least on the intermediate time-scales of interest here, so that only $\varphi$ varies in such displacements.

Only the $P_{rr}$ and $P_{r\theta}$ components of the pressure tensor matter<fnref n="18" />. From Eqs. (A.6) and (A.7), one finds that a necessary condition for $\Delta P_{rr} = 0$ and $\Delta P_{r\theta} = 0$ is that $q \to 0$ as $a \to a_B$ ($c^2$ has no reason to vanish for fluid particles at the boundary, and $\sigma_0$, however small, cannot be strictly zero for any fluid particle moving on the boundary to be part of the fluid and not the vacuum). This conclusion, although based on the specific example of the stress tensor model of Appendix A, is easily generalized to all stress tensor models used in ring studies.

Note that this requirement on $q$ is also sufficient for the right-hand side of Eq. (36) to remain finite. Indeed, computing the mass average of $(dZ/dt)_{vis}$ over the last streamline before the edge (of width $\delta a$) leads to $\langle dZ/dt \rangle_{vis} = i(t_2 + it_1) \exp(i(m\Delta + \gamma)/(an\sigma_0\Delta a)) [t_{1,2}(a_B^+) = 0$ has been used]. Because $t_{1,2} \propto q\sigma_0$, $q \to 0$ as $a \to a_B$ is a necessary for this to remain finite as $\Delta a \to 0$. It is also sufficient if $q$ converges to 0 at least as fast as $\Delta a$.

It is of some interest to explore the consequences of this constraint. In the asymptotic linear WKB analysis of forced density waves ($|ka| \gg 1$), radial derivatives are dominated by the derivative of the phase $m\Delta_0^{dw}$ defined in Eq. (20) so that $q = ad\epsilon/da = -ka\epsilon \sin m\Delta_0^{dw}$. Therefore, requiring $q = 0$ at the edge leads to $m\Delta_0^{dw} + \delta m\Delta = n_r\pi$ ($n_r$ is an integer) or equivalently

$$ \int_0^{x_B} kdx = n_r\pi - \delta m\Delta, \qquad (C.1) $$

where the fact that $n_r$ stands for the number of radial nodes of the edge mode has been anticipated in the notation. Note that the overall change in phase in the “phase loop” (one complete wave cycle from, e.g., the resonance to the edge and back to the resonance) is<fnref n="19" /> $m\Delta_+^{dw} - m\Delta_-^{dw} = 2 \int_0^{x_B} kdx + 2\delta m\Delta = 2n_r\pi$ (the notations have been defined in section 2.2).

<footnote n="18">¹⁸See, e.g., Eqs. (10.48) and (10.49) of Longaretti 2018, taking into account that $\partial/\partial\theta$ terms vanish upon the phase averaging performed to obtain secular contributions.</footnote>

<footnote n="19">¹⁹A difference instead of a sum is involved due to the sign convention of $m\Delta_\pm$.</footnote>
The linear dispersion relation of density waves reads (for $m > 1$):

$$|k| = \frac{3}{2\pi}(m-1)\left(\frac{M_p}{\sigma_0 a_r^4}\right)(a - a_r). \qquad (C.2)$$

Combining Eq. (C.2), (C.1) and (62) for nodeless modes ($n_r = 0$) leads to $\Gamma = |\delta m \Delta| / 2\pi^2$. The quantity $\delta m \Delta$ is not specified by the linear theory of free density waves, and must be quantified numerically. It will be shown in section 5.2 that $\Gamma \simeq 1/16$ for independent modes, but is most likely closer to $1/8$ when mode interactions are taken into account.

Expressions for $m = 1$ involve the gravitational harmonic coefficient $J_2$ (to leading order) and are not given, but they are easily derived from the relevant expansion of $m(\Omega - \Omega_p) - \kappa$ close to the resonance.

This outer boundary condition of a standing wave on the stress tensor is very similar, e.g., to the resonant amplification loop requirement of the SLING amplifier mechanism of Shu et al. (1990). The term “resonant cavity” refers to this feature, and not simply to the fact that the cavity is bounded by a resonant radius. This point is further commented upon in Appendix C.1.3

### Appendix C.1.2. Ring self-gravity: free modes

The outer boundary condition for the ring self-gravity is intriguing. The B ring outer edge is sharp, i.e., the surface density $\sigma_0$ drops to zero on very short scales at the edge (most probably of the order of tens or hundred of meters, i.e., the ring scale height). This outer edge is a free boundary, and the force density due to self-gravity must also remain finite there. In fact, Poisson’s equation ensures that this requirement is satisfied (the force is continuous) at the ring edge and upper and lower surface before the vertical integration and azimuthal averaging leading to Eq. (25) have been performed.

On the other hand, this requirement raises a conundrum as it implies that the self-gravity integral in this equation remains finite as $a \to a_B$ where $a_B$ is the ring edge. Indeed, this integral is defined only in principal value, and the principal value is no longer defined right at the edge. This is clearly related to the fact that the approximations used in the formulation of this integral cannot be valid all the way to the edge in the continuous limit. More precisely, this integral expression fails on distances smaller than the ring thickness; similarly, the streamline width (in a discrete approximation) cannot be smaller than the ring thickness. However, requiring that this
integral is nevertheless relevant or at least usable requires to have either
$\sigma_0 \to 0$ as $a \to a_B$ or $(\epsilon - \epsilon')/(a - a')^2 = q_{aa'}/(a - a') \to 0$ as $a, a' \to a_B$.

Numerical solutions avoid this conundrum if the discretization length $\delta a$ is larger than the zone where the self-gravity integral approximation fails, a feature that is always satisfied in this work. Interestingly, numerical solutions do naturally enforce $q \to 0$ in the $\delta a \to 0$ limit at the edge for the constant surface density profiles examined in this work (in fact, faster than $\delta a$), in line with the argument just made, so that the stress tensor outer boundary condition is naturally enforced by the self-gravitational dynamics as approximated through the self-gravity integral (but the $q$ profile close to the edge may not be the same).

This discussion shows that outer boundary condition ($q \to 0$) is not a
true boundary condition for the self-gravitational dynamics, but rather an
emergent dynamical property at the edge, although it is a true constraint for
the stress tensor.

Appendix C.1.3. Ring self-gravity: forced modes

This discussion involves some of the dynamics of the edge confinement
enforced by Mimas. Although the details of this dynamics are not part of the
numerical solutions explored here, this aspect of the discussion is required in
order to set up the simplified problem examined in this work in a way that
will be minimally changed when the full dynamics is solved for.

Ignoring the details of the edge confinement physics

Forced modes also exist in a cavity between the resonance and the outer
edge, but the underlying physics is different from free modes. On the one
hand, the edge position adjusts so that the viscous angular momentum flux
and the satellite torque balance each other; on the other hand, the resonance
location is produced by the forcing potential. The wider the cavity, the larger
the satellite torque (up to the first wavelength or so) so there is some edge
location satisfying the equilibrium constraint (provided that the full density
wave torque is larger than the viscous torque, but if this were not the case,
the edge would not be located so close to the satellite resonance). This is
what defines the size of the cavity for forced modes.

For free modes, the size of the cavity follows by adjusting the resonance
location so that the “boundary condition” ($q = 0$ at the edge) is satisfied.
This is necessary for the self-gravity integral to be finite at the edge at finite
surface density, as discussed right above. Thus for free modes, the outer
boundary is given (in fact, fixed by the forced mode dynamics as just pointed
out) and the cavity size adjusts, consistently with the boundary condition,
which expresses that the mode is a stationary trapped mode.

There is no physical reason why the resonance cavity sizes of forced modes
should be the same as for the free nodeless mode of same *m*. This can only
happen by coincidence, and, in fact, this is not the case at the B ring edge.

However, the self-gravity integral of the forced mode *must* also be finite
at the edge. How can this be possible with the forced mode cavity being dif-
ferent from the free mode one? In fact, the loop resonant condition Eq. (C.1)
is relevant for both free and forced modes but for different values of δ*m*Δ,
as shown from the numerical results presented below. In forced modes, and
for this reason, the wave phase *m*Δ$_{±}$ will significantly deviate in the outer-
most edge region from the linear forced density wave asymptotic form where
δ*m*Δ = π/4 and ε = cst, as well as from the actual forced density wave
behavior at the edge location. This is explicitly shown in section 6.2. This
deviation is ultimately what makes the ring self-gravity integral finite.

To conclude this part of the discussion, it may be worth pointing out
that in the constant surface density profile model, and when modifying the
forced edge mode satellite torque so that one changes the location of the ring
outer edge, the forcing becomes truly resonant in the usual sense (very large
amplitude of the forced mode) when the edge location nearly matches the free
mode requirement Eq. (C.1) (with the independent free mode requirement
δmΔ ≃ π/8), i.e., when the forced edge mode cavity nearly matches the free
edge mode one<fnref n="20" />. This provides possibly a more direct reason why the free
edge mode cavity is called resonant. In fact, and in the absence of dissipation,
it is impossible to converge numerically on a stationary forced mode structure
when this matching is exact. If such a situation did occur in reality (with
or without dissipation), the configuration would not be stationary and the
satellite/edge distance would quickly adjust in order to break this resonance
condition.

*Accounting for the edge confinement requirements*

For actual ring outer edges, the large enhancement in ring dissipation in

<footnote n="20">²⁰This matching can be realized, e.g., by changing the satellite mass. Also, the ability to change the cavity size implies that, for forced modes, $\delta m \Delta$ is not a fixed quantity; this in itself is sufficient to conclude that the linear density wave asymptotic result, $\delta m \Delta = \pi/4$, has no reason to apply in practice.</footnote>
the forced mode region (Borderies et al. 1982) implies that the viscous angular momentum flux must be nearly canceled for the local balance between the viscous angular momentum flux and the satellite angular momentum extraction to be satisfied. As discussed in section 4.2, this may occur in only one of two ways: $\sigma_0 \to 0$, or $q \to q_a$ as $a \to a_B$ ($q_a$ is the level of non-linearity ensuring the reversal of the viscous flux of angular momentum introduced in section 3.3). Requiring the edge to be sharp makes $q \simeq q_a$ the only viable option in the perturbed edge region. But this apparently conflicts with the fact that $q \to 0$ is also required at the edge, for the self-gravity integral and the stress tensor contribution to remain finite.

The fact that the edge cannot be infinitely sharp is the most simple (and actually only) way out of this conundrum. The surface density must eventually drop to zero close to the edge. In the present analysis, the edge is assumed sharp at scales smaller than a streamline width, whatever the number of streamlines (i.e., the edge is assumed infinitely sharp for all practical purposes). But this is possible only because the physics of the edge confinement is not accounted for. In any case, this is not expected to affect the results in any significant way, as the region of breaking of this assumption is quite narrow.

In fact, the reason why an edge is sharp is somewhat coincidental. Borderies et al. (1989) looked at this for the Encke gap edge and they concluded this was a property of the stress tensor which had no compelling reason to be correct, it just happened to be satisfied for all stress tensor models that had been devised to this date. So it is not exceptional, but not necessary either. If this property did not obtain, the cavity size would still be defined by the torque balance constraint for forced modes, but the edge would not be sharp.

## Appendix C.2. Inner boundary

The inner boundary condition is also suggested by the trapped wave picture of edge modes. Far inside the resonance (in the evanescent region), $\epsilon \propto 1/(a_r - a)$ for forced density waves (see, e.g., Shu 1984). This is also true of pure test particle forced solutions, consistently with the fact that collective effects play little role in the evanescent region, so that this feature is expected to hold for forced edge modes. In practice, numerical solutions show that this asymptotic behavior is already quite precise a fraction of the first wavelength inside the evanescent region for both forced and free modes. It turns out that imposing $\epsilon = 0$ at the inner boundary [i.e., ignoring the role
of streamlines inside the first one in the sum of Eq. (28)] leads to errors in
the percent range if the first streamline lies at least $2\lambda_{dw}$ inside the edge (at
least for the first three eigenmodes) ($\lambda_{dw}$ is the first wavelength of a freely
propagating density wave). This has been checked numerically by doubling
this extent. Therefore, this boundary condition is applied to all numeri-
cal solutions presented in the present work. This boundary condition has
been modified to include enough wavelengths when the first 12 eigenmodes
of $m = 1$ have been computed in section 5.2.

Acknowledgements

I am indebted to Peter Goldreich for many insightful and illuminating
discussions over the years on all aspects of rings dynamics. More specific
inputs of his on the present work have been acknowledged where appropriate
in the text.

This work has also benefited from Phil Nicholson’s careful and relevant
feedback.

I thank the two anonymous referees for their detailed and useful comments
that helped improve the presentation of the work discussed in this paper.

References

Borderies, N., Goldreich, P., Tremaine, S., 1982. Sharp edges of planetary rings. Nature 299, 209-211.

Borderies, N., Goldreich, P., Tremaine, S., 1983. The dynamics of elliptical rings. Astron. J. 88, 1560-1568.

Borderies, N., Goldreich, P., Tremaine, S., 1984. Unsolved problems in planetary ring dynamics. In: Greenberg, R., Brahic, A. (Eds.), Planetary Rings. University of Arizona Press, pp. 713-734.

Borderies, N., Goldreich, P., Tremaine, S., 1985. A granular flow model for
dense planetary rings. Icarus 63, 406-420.

Borderies, N., Goldreich, P., Tremaine, S., 1986. Nonlinear density waves in
planetary rings. Icarus 68, 522-533.

Borderies, N., Goldreich, P., Tremaine, S., 1989. The formation of sharp
edges in planetary rings by nearby satellites. Icarus 80, 344-360.
Borderies, N., Longaretti, P. Y., 1987. Description and behavior of streamlines in planetary rings. Icarus 72, 593-603.

Borderies-Rappaport, N., Longaretti, P.-Y., 1994. Test particle motion around an oblate planet. Icarus 107, 129-141.

Colwell, J. E., Esposito, L. W., Sremčević, M., Stewart, G. R., McClintock, W. E., 2007. Self-gravity wakes and radial structure of Saturn's B ring. Icarus 190, 127-144.

Daisaka, H., Tanaka, H., Ida, S., 2001. Viscosity in a Dense Planetary Ring with Self-Gravitating Particles. Icarus 154, 296-312.

French, R. G., Nicholson, P. D., McGhee-French, C. A., Longaretti, P.-Y., Hedman, M. M., Colwell, J. E., Marouf, E. A., Rappaport, N., Flury, S., Fong, J., Maguire, R., Steranka, G., 2023. Normal modes and their librations at the outer edge of Saturn's B Ring, as observed in Cassini stellar and radio occultation data. Icarus 405, 1156-1178.

French, R. G., Salo, H., McGhee, C. A., Dones, L., 2007. Hst observations of azimuthal asymmetry in saturn's rings. Icarus 189 (2), 493-522.

Goldreich, P., Tremaine, S., 1979a. Precession of the epsilon ring of Uranus. Astron. J. 84, 1638-1641.

Goldreich, P., Tremaine, S., Nov. 1979b. The excitation of density waves at the Lindblad and corotation resonances by an external potential. Astrophys. J. 233, 857-871.

Goldreich, P., Tremaine, S., 1980. Disk-satellite interactions. Astrophys. J. 241, 425-441.

Goldreich, P., Tremaine, S., 1982. The dynamics of planetary rings. Ann. Rev. Astron. Astrophys. 20, 249-283.

Hahn, J. M., Spitale, J. N., Porco, C. C., 2009. Dynamics of the Sharp Edges of Broad Planetary Rings. Astrophys. J. 699, 686-710.

Hedman, M. M., Nicholson, P. D., Salo, H., Wallis, B. D., Buratti, B. J., Baines, K. H., Brown, R. H., Clark, R. N., 2007. Self-Gravity Wake Structures in Saturn's A Ring Revealed by Cassini VIMS. Astron. J. 133, 2624-2629.
Latter, H. N., Ogilvie, G. I., 2009. The viscous overstability, nonlinear wave-
trains, and finescale structure in dense planetary rings. Icarus 202, 565–583.

Longaretti, P.-Y., 1989. Uranian ring dynamics — An analysis of multimode motions. Icarus 82, 281–287.

Longaretti, P.-Y., 1992. Planetary ring dynamics: from Boltzmann's equation to celestial mechanics. In: Benest, D., Froeschle, C. (Eds.), Interrelations Between Physics and Dynamics for Minor Bodies in the Solar System. Editions Frontières. Available on arXiv: astro-ph/1606.00759.

Longaretti, P.-Y., 2018. Theory of Narrow Rings and Sharp Edges. In: Tiscareno, M., Murray, C. (Eds.), Planetary Ring Systems. Cambridge University Press. Available on arXiv: astro-ph/1702.02079.

Longaretti, P.-Y., Borderies, N., 1986. Nonlinear study of the Mimas 5:3 density wave. Icarus 67, 211–223.

Longaretti, P.-Y., Borderies, N., 1991. Streamline formalism and ring orbit determination. Icarus 94, 165–170.

Longaretti, P.-Y., Rappaport, N., 1995. Viscous overstabilities in dense narrow planetary rings. Icarus 116, 376–396.

Mosqueira, I., 1996. Local Simulations of Perturbed Dense Planetary Rings. Icarus 122.

Nicholson, P. D., French, R. G., Hedman, M. M., Marouf, E. A., Colwell, J. E., 2014. Noncircular features in Saturn's rings I: The edge of the B ring. Icarus 227, 152–175.

Nicholson, P. D., French, R. G., McGhee-French, C. A., Longaretti, P.-Y., Hedmann, M. M., El Moutamid, M., Colwell, J., Marouf, E. A., Flury, S., Fong, J., Maguire, R., Steranka, G., ???? The Seven-lobed Shape of Saturn's A Ring Edge. Icarus 390, 115287.

Nicholson, P. D., Hedman, M. M., 2010. Self-gravity wake parameters in Saturn's A and B rings. Icarus 206, 410–423.

Papaloizou, J. C. B., Lin, D. N. C., 1988. On the pulsational overstability in narrowly confined viscous rings. Astrophys. J. 331, 838–860.
Salo, H., Schmidt, J., Spahn, F., 2001. Viscous Overstability in Saturn's B Ring. I. Direct Simulations and Measurement of Transport Coefficients. Icarus 153, 295-315.

Schmidt, J., Salo, H., Spahn, F., Petzschmann, O., 2001. Viscous Overstability in Saturn's B-Ring. II. Hydrodynamic Theory and Comparison to Simulations. Icarus 153, 316-331.

Shu, F. H., 1984. Waves in planetary rings. In: Greenberg, R., Brahic, A. (Eds.), Planetary Rings. University of Arizona Press, pp. 513-561.

Shu, F. H., Tremaine, S., Adams, F. C., Ruden, S. P., 1990. SLING Amplification and Eccentric Gravitational Instabilities in Gaseous Disks. Astrophys. J. 358, 495-514.

Shu, F. H., Yuan, C., Lissauer, J. J., 1985. Nonlinear spiral density waves - An inviscid theory. Astrophys. J. 291, 356-376.

Spitale, J. N., Porco, C. C., 2010. Detection of Free Unstable Modes and Massive Bodies in Saturn's Outer B Ring. Astron. J. 140, 1747-1757.

Toomre, A., 1964. On the gravitational stability of a disk of stars. Astrophys. J. 139, 1217-1238.

Yasui, Y., Ohtsuki, K., Daisaka, H., 2012. Viscosity in Planetary Rings with Spinning Self-gravitating Particles. Astron. J. 143, 110.