Saturation of spiral instabilities in disk galaxies

Spiral density waves can arise in galactic disks as linear instabilities of the underlying stellar distribution function. Such an instability grows exponentially in amplitude at some fixed growth rate $\beta$ before saturating nonlinearly. However, the mechanisms behind saturation, and the resulting saturated spiral amplitude, have received little attention. Here we argue that one important saturation mechanism is the nonlinear trapping of stars near the spiral's corotation resonance. Under this mechanism, we show analytically that an $m$-armed spiral instability will saturate when the libration frequency of resonantly trapped orbits reaches $\omega_\mathrm{lib} \sim \mathrm{a\,\, few}\times m^{1/2} \beta$. For a galaxy with a flat rotation curve this implies a maximum relative spiral surface density $\vert \delta\Sigma/\Sigma_0\vert \sim \mathrm{a\,\,few} \times (\beta/\Omega_\mathrm{p})^2 \cot \alpha$, where $\Omega_\mathrm{p}$ is the spiral pattern speed and $\alpha$ is its pitch angle. This result is in reasonable agreement with recent $N$-body simulations, and suggests that spirals driven by internally-generated instabilities reach relative amplitudes of at most a few tens of percent; higher amplitude spirals, like in M51 and NGC 1300, are likely caused by very strong bars and/or tidal perturbations.


INTRODUCTION
After decades of work there is still no consensus on the origin(s) and character of spiral structure in galaxies.Opinions differ over the mechanisms that provoke spiral responses, the lifetimes of individual spiral patterns, the importance of gas flows and star formation, and so on -for an array of perspectives, see e.g.reviews by Athanassoula (1984); Dobbs & Baba (2014); Sellwood & Masters (2022).However, one area in which there are some robust results -and which we hope bears some relation to real galaxies -is the study of isolated, razor-thin, gas-free stellar disks.In this highly simplified context, both N -body simulations and linear perturbation theory have demonstrated clearly that spiral structure can arise as a linear instability (i.e. an exponentially growing, rigidly rotating Landau mode) of the underlying axisymmetric stellar distribution function (DF) (Sellwood & Kahn 1991;Jalali & Hunter 2005).For example, a 'grooved' disk, in which there is a deficit of stars near a particular galactocentric radius (or more precisely, a localized deficit in the DF of orbital angular momenta), can be linearly unstable to spiral modes which exponentiate out of the finite-N noise (Sellwood & Carlberg 2014;De Rijcke et al. 2019).
The basic premise of the theory promoted by Sellwood & Carlberg over recent decades -and initiated by Sellwood & Lin (1989) -is that once a spiral mode has exponentiated to a large enough amplitude, it will resonantly scatter stars at its associated Lindblad resonances.This resonant scattering both drains the amplitude of the mode and changes the underlying axisymmetric DF (e.g. by carving a new groove, see ⋆ E-mail: chamilton@ias.eduSellwood & Carlberg 2019).The newly-modified axisymmetric DF is itself unstable to a new set of spiral modes, which scatter stars at their own Lindblad resonances, and so on in a 'recurrent cycle of groove modes' (see also Dekker 1976;Sridhar 2019).
On the other hand, a key observable associated with spirals is their amplitude.Observational studies of spiral amplitudes have historically suffered from small sample sizes, with at most ∼ 100 galaxies per study and often much less -see e.g.Rix & Zaritsky (1995); Grosbøl et al. (2004); Block et al. (2004); Elmegreen et al. (2011).But this situation has changed dramatically in the past decade or so due to large-scale spectroscopic surveys like SDSS (York et al. 2000;Strauss et al. 2002), combined with classification efforts from citizen science (Lintott et al. 2008) and machine learning (Domínguez Sánchez et al. 2018), which facilitate statistical studies many thousands of galaxies (Masters et al. 2019;Savchenko et al. 2020;Porter-Temple et al. 2022;Yu et al. 2022).Recent analyses have found that in galaxies without close companions, observable spiral amplitudes range from a few percent to several tens of percent as a fraction of the axisymmetric background; these amplitudes are most strongly (inversely) correlated with the central concentration of the galaxy, and correlate only weakly with other parameters like bar strength, pitch angle, star formation rate, and environmental characteristics (Yu & Ho 2020;Smith et al. 2022).Arms with relative overdensities ≳ 100% seem to exist only in galaxies that are undergoing strong tidal interactions like in M51, or in very strongly barred galaxies like NGC 1300 (Elmegreen et al. 1989;Rix & Rieke 1993;Kendall et al. 2011Kendall et al. , 2015)).
With regard to this observable, the groove-cycle paradigm is incomplete, since it does not answer the question: what sets the maximum amplitude of the recurrent spirals?The answer cannot simply be that the scattering of stars at Lindblad resonances drains energy and angular momentum from the spiral modes -in fact, Sellwood & Carlberg (2022) showed that the mode amplitudes saturate before any significant changes to the DF near Lindblad resonances has occured 1 .Thus, some other process must be halting the exponential growth.One natural candidate is the nonlinear trapping of stars at the spiral's corotation resonance, which transports stars back and forth in angular momentum across the resonance, while leaving their radial actions unchanged.Sellwood & Carlberg (2022) showed that the saturation amplitude of the surface density perturbation scaled as |δΣ| ∝ β 2 , where β is the linear instability's growth rate.This scaling is consistent with equating the growth timescale of the linear mode (∼ β −1 ) with the libration period of nonlinearly trapped orbits (∝ |δΣ| −1/2 ).The possibility that trapping is somehow responsible for saturation has also been mentioned in passing by several other authors (e.g.Kalnajs 1977;Donner & Thomasson 1994), but the precise mechanism behind the saturation, and a quantitative estimate of the resulting amplitude, has been lacking.
In this paper we argue quite generally that saturation of spiral instabilities may well involve nonlinear trapping at corotation, which crucially flattens the angular momentum DF in the vicinity of the resonance, erasing the sharp features that allowed for instability.We derive an analytical formula for the saturated amplitude of a spiral mode whose growth rate is very small -i.e.β/Ωp ≪ 1, where Ωp is the pattern speed -and recover the scaling δΣ ∝ β 2 .While Sellwood & Carlberg (2022)'s spirals do not satisfy the condition β/Ωp ≪ 1 particularly well (in a sense we will define precisely below), the physics behind their saturation is essentially the same as that captured by our calculation, so their spirals obey the same scaling albeit with a modified prefactor.
Inspiration for our analytic calculation was primarily drawn from the classic plasma papers by O'Neil (1965) and Dewar (1973).In the galactic context, the only comparable calculation we know of is by Morozov & Shukhman (1980) (later reproduced as §2.2 of the book by Fridman & Poliachenko 1984), who studied the saturation of amplified spiral perturbations in a stable galactic disk.
This paper is structured as follows.In §2 we explain the physics behind spiral mode amplification and saturation, and thereby calculate the resulting saturation amplitude and estimate the corresponding relative surface density perturbation (some mathematical details are relegated to the Appendix).In §3 we compare our results with those of Sellwood & Carlberg (2022) and with the classic problem of bar-halo friction, and discuss briefly the astrophysical implications of our work.We summarize in §4.
1 More precisely, these authors investigated saturation of spiral instabilties in N -body simulations of their aforementioned grooved razor-thin disks, mostly restricting the sectoral harmonics of the non-axisymmetric gravitational field to azimuthal waveumbers m = 2 or m = 3. See also Sellwood & Lin (1989); Donner & Thomasson (1994).

CALCULATING THE SATURATION AMPLITUDE
To establish the basic idea behind our calculation, we first present Figure 1.The green line in panel (a), which is adapted from Figure 6 of Sellwood & Carlberg (2022), shows the logarithmic amplitude of a spiral density wave from their N -body simulation of an initially grooved, razor thin disk.The gray line illustrates schematically the simplification we will use here: we pretend that the evolution consists of a linear phase in which the spiral amplitude grows exponentially, followed immediately by a nonlinear phase in which the saturated spiral has constant amplitude.In panel (b) we sketch the distribution function (DF) of angular momenta in the disk, F (L).
The initial unstable distribution is shown with a black line.Saturation occurs because nonlinear trapping tends to flatten the DF in the vicinity of corotation (red line), erasing the feature that drove the instability.By calculating the angular momentum content of the spiral in the linear ( §2.2) and nonlinear ( §2.3) regimes, and equating the two, we will be able to read off an approximate expression for the saturation amplitude ( §2.4).
It is worth noting that unstable modes in many different physical systems exhibit very similar behavior to that shown in Figure 1.This is true not only for other unstable stellar disks (Donner & Thomasson 1994) but also for stellar clusters prone to radial orbit instabilities (Palmer 1994;Petersen et al. 2023), for various unstable plasmas (Fried et al. 1971;Berk et al. 1995;Lesur 2020;Bierwage et al. 2021), and in grooved/ridged Hamiltonian Mean Field (HMF) models (Campa & Chavanis 2017;Modak & Hamilton, in prep).Thus, our calculation can, with minor adjustments, be used to predict saturation amplitudes in a whole class of unstable kinetic systems.

Model
We consider a two-dimensional, razor-thin stellar disk, which we embed in a rigid dark matter halo so as to render the disk only weakly unstable at t = 0. Let the gravitational potential of the whole system be Φ, which we decompose into a dominant axisymmetric background part Φ0, and a non-axisymmetric perturbation δΦ which is generated by the non-axisymmetric part of the stellar distribution.
In the disk plane (φ, R) we can describe stellar orbits using angle-action coordinates (Binney & Tremaine 2008): where JR is the radial action and Jφ = R2 dφ/dt ≡ L is the angular momentum.A star orbiting in this disk has Hamiltonian where H0 = v 2 /2 + Φ0(x) and (x, v) are the star's position and velocity.The orbital frequencies are Note that stars on nearly-circular orbits -which is approximately true of most stars in a rotationally-supported, axisymmetric disk -have JR ≈ 0, θφ ≈ φ and Ω ≈ Ω(L).In particular, when making various estimates in this paper we will consider the specific case of a Mestel disk, which is also the model used by Sellwood & Carlberg (2022).In this model, stars on circular orbits have a flat rotation curve: It follows that for circular orbits, L = RV0 and so We also define the distribution function (DF) of stars in the disk, f , such that d 2 v f = Σ(x, t) is the stellar surface density.We further decompose this into an axisymmetric part and a non-axymmetric part: where and the corresponding density perturbation is δΣ(x, t) = d 2 v δf .

Linear regime
Now, at t = 0 we suppose the axisymmetric part of the DF is F (J, 0) = f0(J), and that this DF admits an unstable, m-armed, global spiral Landau mode with growth rate β, which rotates rigidly in the φ direction with pattern speed Ωp = ω/m (so the complex mode frequency is ω + iβ).We will not attempt to calculate β or ω self-consistently in this paper, but this can be done using linear response theory (e.g.Evans & Read 1998;Jalali & Hunter 2005;De Rijcke & Voulis 2016;De Rijcke et al. 2019;Petersen et al. 2023).This spiral mode will exponentiate out of the initial non-axisymmetric fluctuation 2 δf (t = 0), which can either be imposed by hand or can be due to some initial (e.g.finite-N ) noise.After an early transient phase the exponentially growing mode will dominate the potential, so for t ≳ β −1 we can write for some shape function mg(R).Expanding this as a Fourier series in angles, δΦ Now, the rate of angular momentum transfer to a population of stars by a potential perturbation δΦ is Assuming δΦ is weak and F = f0 is time-independent, we can solve the linearized collisionless Boltzmann equation ∂δfn/∂t + in • Ωδfn − i(n • ∂f0/∂J)δΦn = 0 for δfn in terms of δΦn.Plugging the result into (11), using the form (10) for δΦn, and again ignoring terms that are exponentially subdominant for t ≳ β −1 , we arrive at (see Dekker 1976 for details): where ω = mΩp and and we used the shorthand nR = n.We recognise Mm(z) as closely related to the response matrix (e.g.Kalnajs 1977) which encodes the self-consistent response of the system to m-armed perturbations of frequency z (with Im z > 0).In particular, if (8) is a true Landau mode of the system then Im [Mm(ω+iβ)] = 0. From ( 12), this guarantees total angular momentum conservation dL/dt = 0 in the linear approximation.
To gain further insight let us expand the right hand side of (12) for small β/ω ≪ 1.This gives where and M ′ m (ω) ≡ [dMm/dz]z=ω.The terms Lres and Lres correspond to angular momentum transferred to the resonant and non-resonant stars respectively.For instance, writing Im Mm(ω) = Im lim β→0 + Mm(ω + iβ) and using Plemelj's formula it is easily to show that This expression involves only those stars whose orbital frequencies are exactly resonant with the spiral mode.On the contrary, Lnon−res does not involve the resonant stars (its J integral involves a principal value at each resonance location).
Since the total angular momentum is conserved, we see from ( 14) that to first order in β/ω the angular momentum gained/lost by the resonant stars is lost/gained by the non-resonant stars, Lnon−res ≈ − Lres.Physically, the nonresonant stars, which make up the bulk of the system, oscillate in such a way as to sustain the spiral mode (8), while the resonant stars, which make up a small subset of the system, feed angular momentum to the mode and cause it to grow (Nelson & Tremaine 1999;Vandervoort 2003).We can therefore identify L mode (t) ≡ t 0 dt ′ Lnon−res(t ′ ) ≈ − t 0 dt ′ Lres(t ′ ) as the 'angular momentum of the spiral mode'.4Integrating (17) forward in time we get for t ≫ β −1 : It is found empirically that the dominant resonance responsible for angular momentum transfer amplifying spiral modes (at least for groove-induced modes, see Sellwood & Kahn 1991;Palmer 1994 and the discussion in §3.1) is the corotation resonance Ωφ = Ωp.Therefore we will keep only the contributions from n = 0 in (18).Performing the resulting integral over L we find where LCR(JR) is the angular momentum value corresponding to the corotation resonance at a given JR: and is the current magnitude of the m-th Fourier component of the potential fluctuation (10).
The expressions ( 18)-( 19) are of zeroth order in the small parameter β/ω.More generally, if we kept higher order terms in ( 14) we would find that the finite value of β broadens the δ-function in (18) so that resonant stars occupy a width in orbital frequency space of |∆Ω| ∼ β.Converting this to angular momentum space for circular orbits in the Mestel disk gives a corotation resonance width according to linear theory of where Ωφ(LCR, 0) = Ωp.This will be important when interpreting the results of Sellwood & Carlberg (2022) in §3.1.

Nonlinear regime
The above results rely on linear theory, i.e. they assume that stellar orbits are only slightly -and sinusoidally -nudged off their unperturbed trajectories by the spiral potential δΦ.This implies that any changes to the angle-averaged DF F are of second order in the spiral amplitude and so can be ignored to the accuracy we require here.However, sufficiently close to resonance and on sufficiently long timescales, stars will in fact be trapped onto qualitatively different orbits, such that when viewed in the rotating frame of the spiral they librate around a fixed angle.Mathematically, for the corotation resonance this motion can be described using a pendulum equation in the variables (θφ − Ωpt, L) -see the Appendix for details.
In particular, the period of the nonlinear librations is t lib = 2π/ω lib , where ω lib is given in equation (A8).The maximum extent in L of the 'island' of trapped orbits at corotation is given by the island half-width L h -see equation (A9).Suppose for a moment that the spiral amplitude is held constant.Then in and around the librating island the full DF f (θ, J, t) gets progressively more sheared (phase-mixed) along the pendulum Hamiltonian contours.As a result, once t ≫ t lib the DF is a function only of the pendulum Hamiltonian itself (so is uniform along any given Hamiltonian contour), and is independent of oscillation phase (for trapped orbits, this means independent of libration angle) -see e.g.§3.3.3 of Hamilton et al. 2023.Figure 1b shows the corresponding effect that this phase mixing has on the angleaveraged DF of resonant stars F (L, t).For t → ∞ there is a significant change to F in the libration region |L−LCR| ≲ L h compared to t = 0, but little change outside of this region.Most importantly, F has flattened in the vicinity of the resonance.
The above assumes the spiral amplitude is fixed, whereas in truth it is growing, as is the width of the trapped island: L h ∝ |Φm| 1/2 ∝ exp(βt/2).However, the phase mixing occurs on the timescale ∼ ω −1 lib , and as we will see, around the time that the spiral saturates this is actually ≲ β −1 , so the spiral amplitude may be considered quasi-static 5 .Then, the flattening of the axisymmetric part of the DF in the vicinity of the resonance quenches the gradients that were responsible for the mode amplification (e.g. by filling in the groove, see Figure 1b), and so the perturbations cease to grow and the system reaches a stationary state. 6Just as in §2.2, we can calculate the angular momentum excess contained in the nonresonant stars in this saturated state (i.e. the angular momentum of the spiral mode) by setting it equal to minus the change in the angular momentum of the resonant stars compared to t = 0.In the Appendix we show that this is approximately . (23)

Saturation amplitude
Following Figure 1a, we now equate the angular momentum in the spiral immediately before saturation ( 19) with that after saturation ( 23).This requirement leads to an equation of the form dJR Q(Φm, JR) = 0.However, neither linear or nonlinear perturbations at corotation are able to change a star's radial action, so we must have Q = 0 at every JR.This allows us to read off an expression for the final saturation amplitude of the spiral potential at corotation: This result can be made more physically enlightening if we use (A7) and (A8) to express it as In other words, saturation occurs when the nonlinear libration frequency ω lib is a few times larger than the linear growth rate β.Using equations (A8)-(A9) we can estimate the final island width to be which is the same as the linear width ∆L (equation ( 22)).Thus, the librating island of nonlinearly trapped orbits simply grows until it engulfs the original resonance width implied by time-dependent linear theory.Spiral instabilities with higher growth rates thus saturate at larger amplitudes.Note that we did not have to consider explicitly any 'groove' or other feature in the DF in order to arrive at equation (24).The only requirement was that some feature in the DF causes it to be unstable, and that this instability can be quenched by flattening the DF in the vicinity of the corotation resonance (Figure 1).Note also that for one-armed spirals, equation ( 25) is identical to equation ( 12) of Dewar (1973), who studied the saturation of linear momentum weak (logarithmic) divergence, and applies to such a small fraction of stars that we do not expect it to affect our simple estimates.
6 More precisely, one can show that the resulting phase-mixed DF f (which contains both axisymmetric and non-axisymmetric parts) is an approximate soliton; that is, a self-consistent, nonlinear, rigidly-rotating solution to the combined Vlasov-Poisson system of equations, much like a stellar bar (Vandervoort 2003).transfer between electrons and a single longitudinal plasma wave.
Let us evaluate equation ( 24) for the Mestel disk.Using (5) we find We can further relate this to a relative surface density perturbation at corotation, |δΣm(RCR)/Σ0|.Using equation ( 27) we can easily estimate this to be of order ∼ |Φm|/V 2 0 ∼ 4m −1 (β/Ωp) 2 .To make this estimate more precise, we must relate the potential of the spiral to the corresponding surface density perturbations.This can be done in the WKB approximation provided the radial wavelength λ of the spiral is not too large; in this case |Φm| ≈ Gλ| δΣm| (Binney & Tremaine 2008, equation (6.30)).We now divide this by the surface density of the Mestel disk Σ0 = ξV 2 0 /(2πGR), where ξ is the disk mass fraction (so ξ = 1/2 for a half-mass disk).Evaluating at the corotation radius, we find where α ≡ arctan[λm/(2πRCR)] is the pitch angle at corotation.Note that when written in this form, (28) does not depend explicitly on m.

Comparison with Sellwood & Carlberg (2022)
Sellwood & Carlberg (2022) considered spiral saturation in N -body simulations of initially axisymmetric half-mass Mestel disks (Binney & Tremaine 2008).They generated the instabilities by carving a groove in the DF of these disks centered on a particular location in angular momentum space at t = 0, namely L * ≈ 6.5R0V0 (where R0 is their length unit) -see Figure 1b for a schematic illustration.By varying the width of the groove or the level of random motion in the disk, they were able to generate spirals with very similar pattern speeds Ωp but different growth rates β.Their Figure 11 shows the saturation amplitude of the resulting spiral instabilities as a function of β/Ωp.To compare our results with those of Sellwood & Carlberg (2022), in Figure 2 we plot |δΣm/Σ0| as a function of the dimensionless growth rate β/Ωp, following equation ( 28), for different pitch angles α.With red markers we show the 6 data points for the m = 2 spirals extracted by Sellwood & Carlberg (2022) in their Figure 11.We see that the theory is an excellent fit to the data for α ≈ 75 • .However, inspection of Sellwood & Carlberg (2022)'s figures suggests their m = 2 spirals have pitch angles closer to α ∼ 35 • , in which case the prediction ( 28) is consistently too large by a factor of a few, although the scaling |δΣm/Σ0| ∝ (β/Ωp) 2 remains accurate.
Why might our the theory be systematically overestimating Sellwood & Carlberg (2022)'s spiral amplitudes?The main reason is that (28) was derived in the limit of very small β/Ωp, and this is not a good approximation for Sellwood & Carlberg (2022)'s spirals.One way to see this clearly is to note that in this limit, all nonzero contributions to the angular momentum transfer between resonant and non-resonant stars in the linear regime occur at exact resonant locations in action space (see equation ( 18)).This cannot be true for Sellwood & Carlberg (2022)'s spirals, since the groove (centered on L * = 6.5R0V0 and with a width of wL ∼ 0.2R0V0) does not coincide with the corotation resonance (LCR ≈ 7.1R0V0 for m = 2)7 or any other major resonance of the half-mass Mestel disk (see the top panel of Sellwood & Carlberg 2022's Figure 4).In other words, if we were to calculate angular momentum transfer using equation (18) using their grooved DF, the grooved region itself would contribute almost nothing, which is clearly absurd.Instead, a more accurate calculation of the mode's angular momentum content in the linear regime should involve an integral over each broadened resonance, with width (equation ( 22)): This broadened function, when centered on the corotation resonance, engulfs the groove region for all six of Sellwood & Carlberg (2022)'s β/Ωp values, so the mode can feed off the sharp gradients in the DF created by the groove, as illustrated in Figure 1b.Despite this subtlety, the physics of mode saturation is essentially the same as in our simplified calculation: nonlinear trapping at corotation and subsequent phase mixing flattens the DF, filling in the grooved region and hence quenching the instability.This explains why Sellwood & Carlberg (2022) still find the scaling |δΣm/Σ0| ∝ β 2 predicted by ( 24).The appropriate prefactor in their case is smaller, likely because to quench the instability it is not necessary for the DF to be completely phase mixed in the resonant region; it is sufficient merely to fill in a significant portion of the groove.

Comparison with bar-halo friction
It is tempting to draw an analogy between our calculation and the classic problem of bar-halo friction in galaxies (Tremaine & Weinberg 1984;Chiba & Schönrich 2022;Chiba 2023;Hamilton et al. 2023).These problems are indeed closely related: in both cases, angular momentum is transferred to/from a global, rigidly rotating 'mode' (spiral or bar), and the equations governing this transfer in the linear regime ( §2.2) look very similar.Moreover, in both cases the angular momentum transfer can be quenched due to the trapping of orbits at (e.g.) corotation.However, the precise physics involved is subtly different.To illustrate this, let us consider the simple case of a constant-amplitude bar whose pattern speed Ωp is decreasing monotonically.The decrease in Ωp causes e.g. the corotation resonance location to sweep to larger radii (larger angular momenta LCR).In particular, if the resonance location moves by a resonance width ∼ L h on a timescale short compared to the libration time of trapped orbits t lib (the so-called 'fast regime', see Tremaine & Weinberg 1984), then the bar will always be encountering fresh, untrapped material with which to resonate.Mathematically, in this regime nonlinear trapping can be ignored and linear theory will continue to be valid (see Chiba 2023 for a general treatment).On the other hand, in the case of spiral Landau modes, the phase space location of the corotation resonance is fixed, and only its width grows with time until eventually it is comparable to the width implied by the finite β in linear theory (compare equations ( 22)-( 26)).Thus in the absence of diffusive processes, secular evolution, etc. (see §3.3), amplifying spirals never encounter 'fresh material' at all.They simply grow until the trapping region has exhausted the initial population of resonant orbits (within |∆Ω| ∼ β −1 of the resonance) off of which the spiral was feeding.

Astrophysical implications
Saturation amplitudes of spirals are of key astrophysical importance because they may allow one to distinguish between different theories of spiral structure (D'Onghia et al. 2013;Sellwood & Carlberg 2021).More generally, spiral strength is a key parameter in determining the secular evolution of galactic disks, their associated gas flows, star formation rates, etc. (Roberts 1969;Binney & Lacey 1988;Sellwood & Binney 2002;Slyz et al. 2003;Kim & Kim 2014).
Taken at face value, the result (28) suggests that spiral structure driven by instabilities will only ever reach a relative amplitude of several tens of percent compared to the axisymmetric component of a galaxy.This is because (i) most spirals have cot α ≲ a few (Lingard et al. 2021), and (ii) it is difficult to imagine any sustainable cycle of instabilities whose individual mode strengths amplify faster than ∼ Ωp, so typically (β/Ωp) 2 ≪ 1.This result is consistent with observations of isolated galaxies; overdensities of 100% or greater are mostly confined to galaxies undergoing strong tidal interactions (Elmegreen et al. 1989;Rix & Rieke 1993;Kendall et al. 2011Kendall et al. , 2015)).
There may, however, exist a slight tension between theory and observation given that equation ( 28) predicts that more tightly-wrapped spirals will have larger amplitudes, whereas Díaz-García et al. ( 2019) and Yu & Ho (2020) found (very weak) correlations in the opposite sense in S 4 G and SDSS data respectively.On the other hand, β, Ωp and α are emphatically not independent variables, and since self-gravity is weaker for tighter spirals, a decrease in α will correspond to a decrease in β/Ωp.
It is enlightening to estimate the timescale tsat over which saturation of realistic spirals might be achieved.To do this, first recall that saturation occurs when ω lib (tsat) ∼ β.Second, a basic scaling of equation (A8) tells us that at corotation, ω lib (t) ∼ ϵ(t) 1/2 Ωp, where ϵ(t) ≡ |Φm(t)/Φ0| is the dimensionless strength of the perturbation.Third, in the linear regime we have ϵ(t) = ϵ(0) exp(βt).Putting these three ingredients together, we get a saturation time In particular, if the initial fluctuation level is set by Poisson noise alone then ϵ(0) ∼ N −1/2 , and we have (very roughly): where TCR = 2π/Ωp is the orbital period at corotation.Thus, if it has to exponentiate out of Poisson noise, a rather vigorously growing spiral with β/Ωp ≳ 0.1 will saturate after a few Gyr.More realistic galactic noise tends to be of higher amplitude than is implied by Poisson statistics (Sellwood 2012;Fouvry et al. 2015), though this will only change the timescale tsat by a logarithmic factor.A caveat to our work is that we have considered a highly idealized system in which we isolated the interaction of a single spiral mode with an ambient stellar population.In reality, there can be multiple modes at play simultaneously.Strong coupling between different modes can cause spirals to saturate (Sygnet et al. 1988;Laughlin et al. 1997), although Sellwood & Carlberg (2022) found that mode-coupling was not a significant effect in their simulations.Still, superposing many additional weak modes on top of our primary one will have some diffusive effect on stellar orbits, as will as other noise sources like molecular clouds.Moreover, spiral waves can be damped by hydrodynamic interactions with the interstellar gas (Kalnajs 1972;Roberts & Shu 1972;Dekker 1976).One can get an idea of the impact these effects will have on our primary spiral mode by characterising them all, very crudely, by a simple 'dissipation timescale' ν −1 .Then in the linear regime, the spiral growth rate β will be replaced by β − ν, meaning that if ν ≳ β the mode in question will never grow out of the noise.In the more realistic case that ν ≪ β, the linear physics will proceed essentially unchanged from the case ν = 0. but as the mode amplitude approaches the nonlinear regime, one effect of the additional forces will be to de-trap stellar orbits (Johnston 1971;Auerbach 1977;Hamilton et al. 2023).This will counteract the flattening of the DF near corotation and hence allow the spiral amplitude to continue to grow.Thus, we speculate that more noisy/gas rich galaxies may allow for somewhat higher amplitude spirals.

SUMMARY
A complete theory of spiral structure should explain not only the origin, but also the amplitude, of spiral waves.It is well known that, at least in isolated, razor-thin, gas-free stellar disks, spiral structure can arise as a linear instability of the underlying DF.In this paper we have calculated the saturation amplitude of such an instability under the assumption that transferral of angular momentum between stars and spiral modes primarily occurs at the corotation resonance.We find that (i) saturation occurs when the nonlinear resonant libration frequency of the stars is comparable to the initial linear growth rate of the instability β, and (ii) for galaxies with flat rotation curves, the resulting relative saturation amplitude in surface density is ∼ (β/Ωp) 2 cot α, where Ωp is the spiral pattern speed and α is its pitch angle.

Figure 1 .
Figure 1.(a) Growth and saturation of a spiral mode in an Nbody simulation of a razor thin stellar disk (adapted from Sellwood & Carlberg 2022, green line).The gray line illustrates the simplification used in this paper: we separate the spiral's life into a linear phase of perfectly exponential growth, and a nonlinearly saturated phase of constant amplitude.(b) Illustration of the change to the axisymmetric DF F (L, t) in the vicinity of a corotation resonance.The initial instability is driven by some sharp feature in the DF (here a groove, centered on L = L * with width w L , see Sellwood & Carlberg 2022).In the saturated state the DF is flattened around the resonance L = L CR .

Figure 2 .
Figure 2. Saturation amplitude of spiral instabilities at corotation as a function of dimensionless growth rate.The theoretical result (28) is shown for ξ = 1/2 and different values of pitch angle α (solid lines).The results of Sellwood & Carlberg (2022)'s simulations of m = 2 spirals are shown with red circles (c.f.their Figure 11).