Spiral instabilities: Mechanism for recurrence

TL;DR

This paper proposes that recurrent spiral patterns in stellar discs are caused by self-excited instabilities, specifically groove-type instabilities triggered by scattering at Lindblad resonances, explaining persistent spiral activity.

Contribution

It demonstrates that scattering at Lindblad resonances seeds new instabilities, providing a comprehensive mechanism for recurrent spiral patterns in stellar discs.

Findings

Spiral patterns result from superposition of transient modes.

Resonance scattering at Lindblad resonances triggers new instabilities.

The mechanism explains persistent spiral activity in galaxies.

Abstract

We argue that self-excited instabilities are the cause of spiral patterns in simulations of unperturbed stellar discs. In previous papers, we have found that spiral patterns were caused by a few concurrent waves, which we claimed were modes. The superposition of a few steadily rotating waves inevitably causes the appearance of the disc to change continuously, and creates the kind of shearing spiral patterns that have been widely reported. Although we have found that individual modes last for relatively few rotations, spiral activity persists because fresh instabilities appear, which we suspected were excited by the changes to the disc caused by previous disturbances. Here we confirm our suspicion by demonstrating that scattering at either of the Lindblad resonances seeds a new groove-type instability. With this logical gap closed, our understanding of the behaviour in the simulations is almost complete. We believe that our robust mechanism is a major cause of spiral patterns in the old stellar discs of galaxies, including the Milky Way where we have previously reported evidence for resonance scattering in the recently released Gaia data.

Tables

Table 1. Table 1: Numerical parameters

Grid size	106 $\times$ 128
Active sectoral harmonic	2 or 3
$R_0$	8 grid units
Softening length	$R_0/8$
Number of particles	$5\times {10}^7$
Basic time-step	$R_0/\left(40V_0\right)$
Time step zones	5
Guard zones	4

Table 2. Table 2: Summary of the modes fitted to data from simulations described in § 3 , including the fitted complex frequencies, ω 𝜔 \omega with β 𝛽 \beta being the growth rate of the mode, and the radii of the principal resonances for circular orbits given by ( m + l ​ 2 ) / ℜ ⁡ ( ω ) 𝑚 𝑙 2 𝜔 (m+l\surd 2)/\Re(\omega) , although we leave uninteresting entries blank for clarity.

Line #	Model	$m$	Mode #	time range	$\omega =m{\Omega }_p+i\beta$	$R_{\mathrm{ILR}}$	$R_{\mathrm{CR}}$	$R_{\mathrm{OLR}}$
1	G	2	1	$100\le t\le 200$	$0.656\pm 0.005+\left(0.017\pm 0.003\right)i$	$\sim 0.89$	$\sim 3.05$	$\sim 5.20$
2	G	2	2	$200\le t\le 500$	$0.393\pm 0.003+\left(0.022\pm 0.002\right)i$		$\sim 5.09$
3	GR	2	1	$400\le t\le 700$	$0.398\pm 0.001+\left(0.019\pm 0.002\right)i$		$\sim 5.03$
4	GR	2	2	$400\le t\le 700$	$0.672\pm 0.001+\left(0.013\pm 0.004\right)i$		$\sim 2.98$
5	GSOLR	2	1	$400\le t\le 700$	$0.400\pm 0.008+\left(0.019\pm 0.005\right)i$		$\sim 5.00$
6	GSCR	2	1	$400\le t\le 700$	$0.662\pm 0.030+\left(0.010\pm 0.008\right)i$		$\sim 3.02$
7	E	3	1	$300\le t\le 600$	$0.233\pm 0.001+\left(0.005\pm 0.001\right)i$	$\sim 6.8$	$\sim 12.8$	$\sim 18.9$
8	E	3	2	$500\le t\le 780$	$0.417\pm 0.001+\left(0.020\pm 0.006\right)i$		$\sim 7.2$
9	ER	3	1	$600\le t\le 800$	$0.423\pm 0.001+\left(0.014\pm 0.001\right)i$		$\sim 7.1$
10	ER2	2	1	$600\le t\le 900$	$0.271\pm 0.002+\left(0.013\pm 0.002\right)i$		$\sim 7.4$
11	ER2	2	2	$600\le t\le 900$	$0.136\pm 0.006+\left(0.022\pm 0.004\right)i$		$\sim 14.7$
12	GS3OLR	3	1	$400\le t\le 680$	$0.614\pm 0.009+\left(0.013\pm 0.001\right)i$		$\sim 4.9$

Equations

$$\delta \Sigma (R,\varphi ,t)=A_m(R)e^{i(m\varphi -\omega t)},$$

$$m\left[{\Omega }_p-{\Omega }_c(R)\right]=l\kappa (R),$$

$$m\approx \frac{R_{\mathrm{CR}}{\kappa }^2}{2\pi XG\Sigma },$$

$$f(E,L_{z})=\cases{xFL_{z}^{q}e^{-E/\sigma_{R}^{2}}&$L_{z}>0$\cr 0&otherwise,\cr}$$

$$F=\frac{1}{G{R}_0(R_0{V}_0{)}^q}\frac{(q/2+0.5{)}^{q/2+1}}{{\pi }^{3/2}(q/2-0.5)!}.$$

$$T(L_z)={\left[1+{\left(\frac{R_0{V}_0}{L_z}\right)}^{\nu }\right]}^{-1}{\left[1+{\left(\frac{L_z}{R_1{V}_0}\right)}^{\mu }\right]}^{-1},$$

Full text

††pubyear: 2019††pagerange: Spiral instabilities: Mechanism for recurrence–References

1 Introduction

Theorists have long dreamed that the spiral patterns gracing most disc galaxies, and simulations thereof, reflect the normal modes of the disc. That is, they are mildly non-linear manifestations of self-excited, linear instabilities of the stellar disc that are uniformly rotating and exponentially growing density waves. Long ago, this grand program was set back by two early discoveries from normal mode analyses of apparently reasonable models of featureless discs: models where the rotation curve rose gently from the centre were dominated by bar-forming instabilities (Hohl, 1971; Kalnajs, 1978), while models having a dense (bulge-like) centre had no modes whatsoever (Toomre, 1981).

Here we begin by reviewing why we think spiral modes are in fact the most promising mechansim, despite this early setback, and add to our case by demonstrating how spiral modes recur.

1.1 Modes in galactic discs

A normal mode of any system is a self-sustaining, sinusoidal disturbance of fixed frequency and constant shape, save for a possible uniform rotation; the frequency would be complex if the mode were to grow or decay. In the case of galaxy discs, the perturbed surface density of a mode is the real part of

[TABLE]

where $m$ is the angular periodicity, $\omega=m\Omega_{p}+i\beta$, $\Omega_{p}$ is the angular rate of rotation, usually called the pattern speed, and $\beta$ is the growth rate. The complex function $A_{m}(R)$, which is independent of time, describes the radial variation of amplitude and phase of the mode.

Stability analysis of a system supposes small amplitude perturbations about the equilibrium state and is linearized when any terms that involve products of small amplitude terms are discarded – see Kalnajs (1971) for a careful formulation. The equilibrium is linearly unstable if any of the resulting normal modes have a positive growth rate, since its amplitude will exponentiate out of the noise until the neglected 2nd and higher order terms become no longer negligible.111The word instability is sometimes used to imply a purely imaginary frequency, and a mode with a complex frequency is then described as an overstability. Here we adopt the more usual convention of simply describing all frequencies with $\beta>0$ as instabilities.

Normal modes can be standing wave oscillations of the system that exist between two reflecting barriers, as in organ pipes and guitar strings, which are generally described as cavity modes in galaxy discs. The prime example in galaxies is the bar-forming mode, for which reflections take place at the centre and at corotation (Toomre, 1981; Binney & Tremaine, 2008). An instability of this type is possible only if the disturbance has no inner Lindblad resonance (hereafter ILR), since linear theory (Mark, 1974) predicts that any disturbance that encounters an ILR will be absorbed, and therefore damped. An ILR must be present for any reasonable pattern speed when the centre is dense, and therefore no small-amplitude cavity mode is possible in a featureless disc of this kind. But $\Omega_{c}-\kappa/2$ has a maximum value in mass models having quasi-harmonic cores, and bar-forming instabilities avoid resonance damping by having pattern speeds that exceed this maximum. Here $\Omega_{c}(R)$ is the angular frequency of circular motion at radius $R$ in the disc mid-plane, and $\kappa(R)$ is the usual frequency of small-amplitude radial oscillations about a circular orbit (Binney & Tremaine, 2008). The dominant mode of several bar-unstable models has been identified in simulations, with excellent quantitative agreement of both the frequency and mode shape (Sellwood & Athanassoula, 1986; Earn & Sellwood, 1995). Other types of cavity mode have been proposed by Mark (1977) and by Sellwood & Carlberg (2014) and are described below.

Galactic discs can also support another class of mode. The best known examples are edge modes (Toomre, 1981; Papaloizou & Lin, 1989) and groove modes (Sellwood & Lin, 1989; Sellwood & Kahn, 1991). They are not standing waves, however, and maintain a fixed shape and frequency by other means. In the case of the edge mode, a small non-axisymmetric distortion of the disc where the density decreases steeply, moves high density material out to places where the equilibrium density was lower, and conversely at other azimuthal phases. On their own, such co-orbiting distortions would be neutrally stable and therefore of no interest. But as Julian & Toomre (1966) taught, a cool surrounding disc responds vigorously to a co-orbiting mass excess, creating a trailing wake that extends radially far into the shear flow on either side of the perturbing mass; this behaviour is a consequence of swing amplification (Goldreich & Lynden-Bell, 1965; Julian & Toomre, 1966; Toomre, 1981; Binney, 2019). The distorted edge therefore excites a strong supporting response from the interior disc that is not balanced by the exterior response because the equilibrium density drops rapidly with radius at the edge. The attraction of the interior wake on the original density excess increases its angular momentum, causing it to grow exponentially as it rotates. The necessary conditions for instability were summarized by Toomre (1989).

A groove in a disc is effectively two closely spaced edges, which however give rise to a single mode because the distortions on each edge are gravitationally coupled. Note that it is a steep gradient in the angular momentum density that matters for both edge and groove modes; epicyclic blurring can mask the steepness in the actual surface density profile. As for edge modes, it is the supporting response of the surrounding disc that causes the groove mode to have a substantial radial extent and to grow rapidly (Sellwood & Kahn, 1991), and these authors were able to obtain good quantitative agreement between their analytic predictions and simulations.

Note that swing amplication, famously illustrated in the dust-to-ashes figure of Toomre (1981), is not a mode both because the shape changes with time and the amplitude variation is not a simple exponential. Also the wake response to an imposed co-orbiting mass clump (Julian & Toomre, 1966) is not a mode because, to first order, it would disperse if the clump were removed, and it therefore is not self-sustaining. Both are simply responses of the disc to externally imposed disturbances. However, they are both very helpful concepts when trying to understand the mechanisms of self-sustaining modes.

1.2 Noise, resonances, and heating

The collisionless fluid of stars invoked in these analytical treatments is an idealization that assumes the stars to be infinitely finely divided so that phase space is smooth. The numbers of stars in galaxy discs is large enough that this assumption holds quite well (see Sellwood, 2014a, for caveats), but galaxies contain mass clumps such as star clusters and giant molecular clouds, and the number of particles employed in a simulation is generally several orders of magnitude fewer than the number of stars. Thus shot noise in both real galaxies and in simulations gives rise to significant inhomogeneities in the disc.

The spectrum of shot noise in a shearing distribution of randomly distributed gravitating masses inevitably contains leading wave components that will be strongly amplified as the shear carries them from leading to trailing. This behaviour has two important consequences. The first consequence of swing-amplified shot noise is that each heavy particle develops a wake (Julian & Toomre, 1966; Binney, 2019). The two point correlation function of the particles becomes greater along the direction of the wake and lower in other directions, and the particle distribution is said to be polarized. Since the distribution of particles is no longer perfectly random, the amplitude of all components of the noise spectrum is enhanced, causing subsequent noise-induced fluctuations to be stronger, although linear theory predicts this cycle should asymptote in a few epicycle periods to a mean steady excess over the level expected from uncorrelated noise (Julian & Toomre, 1966; Toomre & Kalnajs, 1991).

Second, the collective amplified response of any one component of the noise orbiting at the angular frequency $\Omega_{p}$, creates a coherent trailing wave in the disc that propagates away from corotation (Toomre, 1969, 1981) until it reaches a Lindblad resonance where it is absorbed (Mark, 1974). For near circular orbits the Lindblad resonances occur where the Doppler shifted frequency at which the stars encounter the wave is equal to the epicyclic frequency of their small radial excursions, i.e., the radii at which

[TABLE]

with $l=\pm 1$ at Lindblad resonances. The negative sign is for the ILR, where stars overtake the wave, and the positive is for the OLR where the wave overtakes the stars, at the local epicycle frequency in both cases. Wave-particle interactions at the resonance cause localized irreversible changes to the energy and angular momenta of stars. Jacobi’s invariant (Binney & Tremaine, 2008) implies that changes are related as $\Delta E=\Omega_{p}\Delta L_{z}$. On average and to second order, particles lose $L_{z}$ at the ILR and gain at the OLR (Lynden-Bell & Kalnajs, 1972; Carlberg & Sellwood, 1985), and this outward transfer of $L_{z}$ allows the wave to extract free energy from the galactic potential enabling the scattered particles to acquire additional random energy at both resonances. The resulting depopulation of stars on near circular orbits over the narrow region of each resonance creates a “scratch” in the disc that may alter its stability properties.

It is important to realize that linear theory neglects this second order effect by assumption, i.e. it does not allow for changes to the equilibrium state. In fact, Sellwood (2012) found that the amplitudes of successive episodes of uncorrelated swing amplified noise rose slowly, but continuously, because the consequent scratches to the disc from each episode caused partial reflections of subsequent disturbances that allowed further amplification (Sellwood & Carlberg, 2014; Fouvry & Pichon, 2015). This process continued until the the partial reflections became strong enough that the disc was able to support an unstable mode (Sellwood, 2012; De Rijcke et al., 2019), and coherent growth to large amplitude began. We discuss this behaviour further in §1.4.

Here we have used the word “scratch” to describe quite mild changes to the distribution function (hereafter DF) that can cause partial reflections of a wave propagating radially within the disc. But scattering at a Lindblad resonance could also carve a similar feature that seeds a groove mode instead, and we will show below that this appears to be the more usual behaviour.

Note that the Lindblad resonances are closer to corotation, where $\Omega_{p}=\Omega_{c}(R_{\rm CR})$ ($l=0$ in eq. 2), when $m$ is large than for waves of lower $m$, implying that more free energy can be extracted from the potential, causing more rapid heating when the disc supports larger-scale waves. The value of $m$ that is amplified most strongly (Julian & Toomre, 1966; Toomre, 1981) is

[TABLE]

with $1<X<2$ in a flat rotation curve. Thus the preferred $m$ varies inversely with the disc surface density $\Sigma$ (Sellwood & Carlberg, 1984; Athanassoula et al., 1987). Whatever the origin of disturbances, a sub-maximal disc will prefer higher $m$, and therefore heat more slowly than would heavier discs.

Since spiral activity heats collisionless particles, or stars, an uncooled disc must become less able to support collective disturbances over time, as is well known. Sellwood & Carlberg (1984) established that dynamical cooling by gas and star formation is needed to counter secular heating and to maintain spiral activity, a result that has been confirmed in many more recent simulations (e.g. Roškar et al., 2008; Agertz et al., 2011; Aumer et al., 2016). It nicely accounts for the observation (e.g. Oort, 1962) that almost all spiral patterns are seen in galaxies that contain gas and are forming stars.

1.3 The spiral challenge

As already noted, a galaxy model having a dense centre and no sharp features, such as an edge or a groove, should not support any normal modes at all. Since spirals are ubiquitous in disc galaxies containing a modest gas fraction, and also develop spontaneously in simulations of isolated discs, some mechanism is needed to excite them. Currently, there are at least four proposed mechanisms:

1. i

Following Mark (1977), Bertin & Lin (1996) suggest that spirals result from a cavity type mode in a low-mass disc that is dynamically cool over most of the disc, but which also possess an inner “$Q$ barrier” to shield the ILR. These authors imagine that most galaxies support a single, long-lived, mildly-unstable mode that persists for many tens of galactic rotations and becomes “quasi-steady” due to dissipative shocks in the gas, but allow that superposition of a second mode may be needed in some cases. As given by eq. (3), strong swing-amplification for $X\sim 2$ sets a preference for multi-arm instabilities over 2-arm modes in low mass discs. By considering only bi-symmetric disturbances in low-mass discs, Bertin et al. (1989) exploited the mild amplification when $X>3$ in order to obtain slowly-growing spiral modes in their global stability analysis of many galaxy models.

Simulations by Sellwood (2011) of one of the cases presented by Bertin et al. (1989) confirmed that a single, slowly-growing mode was present when disturbance forces were restricted to $m=2$. Not surprisingly, however, he also found much more vigorous instabilities appeared when higher sectoral harmonics contributed to disturbance forces, and the contrived basic state of the disc that was designed to support the $m=2$ mode was rapidly changed by these more vigorous instabilities. This evidence alone should have ruled out the theory, although Shu (2016) ignored it as he continued to advocate for their picture. 2. ii

Toomre (1990) and Toomre & Kalnajs (1991) abandoned the idea of spirals as normal modes, and advocated instead that a collection of massive clumps in the disc, each of which becomes dressed with its own wake, would create a “kaleidoscope” of shearing spiral patterns. Their simulations of this process were all confined to a single shearing patch with a modest number of particles. D’Onghia et al. (2013) conducted global simulations of a low-mass disc, embedded in a rigid halo, employing $10^{8}$ star particles to which they added a sprinkling of heavy particles that induced evolving multi-arm spiral patterns in the stars. In a separate experiment they also tried a single perturber of mass $10^{7}\;$M⊙, which they removed again after it had completed one orbit; the simulation continued to manifest spiral activity in response to the non-axisymmetric density distribution created by the original imposed mass, which they described as non-linear behaviour.

These authors suggest that “ragged” spiral activity in galaxies results from such responses to co-orbiting giant molecular clouds, massive star clusters, etc., and to the lingering disc responses should any disperse. Although we do not doubt their numerical results, we remain unconvinced of their importance for spiral activity in galaxies. First, the $10^{7}\;$M⊙ particle that D’Onghia et al. (2013) employed produced only a modest, and not very extensive wake in their low-mass disc. Spirals in real galaxies generally have greater amplitude and radial extent than this, suggesting that yet more massive clumps would be needed. Second, a collection of randomly placed heavy particles does not seem likely to produce a net response that is predominantly 2- or 3-armed, as are observed in the overwhelming majority of galaxies (Davis et al., 2012; Hart et al., 2016; Yu et al., 2018). Third, clumps massive and numerous enough that their associated wakes produce large-amplitude and radially-extensive spiral patterns would scatter disc stars, and heat the disc rapidly so that the responses would fade quickly unless the disc were cooled aggressively, and the necessary cooling (Toomre, 1990) seems rather extreme. A fourth, and overriding, reason is that the idea is unnecessary, because discs readily support unstable spiral modes, as we discuss next. 3. iii

Sellwood & Carlberg (2014) demonstrated that spiral activity in simulations resulted from superposition of a number of coherent, uniformly rotating waves. They found that each wave grew and decayed, but was detectable over a period of some ten rotations at its corotation radius, and they presented substantial evidence that the waves were in fact modes. The vigorous modes in their new picture differ substantially from those invoked by Bertin & Lin (1996) because they work best where swing-amplication is strongest, do not last for nearly as long, and fresh instabilities develop to maintain spiral activity.

Note that the spiral appearance changes still more rapidly than do the modes in their simulations. This is because the superposition of several modes, each having a different pattern speed and perhaps also angular periodicity, as well as time varying amplitude, causes the pattern of visible spiral arms to change radically in less than one orbit.222An animation showing the time evolution of the net density when two notional patterns are superposed is at

http://www.physics.rutgers.edu/$\sim$sellwood/spirals.html

Sellwood & Carlberg (2014) also postulated a new cavity mode that relied upon partial reflection of travelling waves off an impedance variation created by a previous perturbation within the disc. They suggested that ILR scattering by a past wave would have “scratched” an originally smooth distribution function (see §1.2) to create a deficiency of stars on near circular orbits at some radius. A subsequent ingoing trailing wave encountering the scratch before reaching its own ILR, would find that the scratch in the otherwise smooth disc represents an abrupt change of impedance, and is therefore partially reflected into a outgoing leading wave. A second reflection of the wave by swing amplification at corotation creates a cavity that will support standing waves having frequencies allowed by the usual phase closure condition. They dubbed this new type of cavity mode a “mirror mode.” As for the bar mode that also includes swing amplification from leading to trailing, a reasonable reflected fraction at the partial mirror will cause a mirror mode to grow rapidly, and the instability must run its course within a modest number of pattern rotations. Nevertheless, spiral activity can be maintained because other modes develop in rapid succession. We present evidence in §1.4 below that the first real instability in the tests of the Mestel discs reported in Sellwood (2012) was almost certainly of this type. Some mirror modes may also have been present in the simulations of Sellwood & Carlberg (2014), but we here (§3) argue for a more probable recurrence cycle. 4. iv

Grand et al. (2012a, b), Baba et al. (2013), Roca-Fàbrega et al. (2013), and others present simulations of mostly sub-maximal discs and report that spiral patterns are shearing structures, in which the density maxima wind as would material arms, or nearly so. Grand et al. (2012a, b), Kawata et al. (2014), and other authors have described spiral arm streaming motions that are consistent with expectations set out by Kalnajs (1973); however, such a flow pattern is required for any self-consistent spiral, irrespective of its nature or origin. Baba et al. (2013) and Michikoshi & Kokubo (2018) find that spirals in their simulations behave as predicted by swing amplification theory (Toomre, 1981). Grand et al. (2013) reported a correlation between the mean pitch angle of spiral arms in their simulations and the shear rate in the disc, which they also suggested was consistent with swing amplification theory. Baba (2015) reported similar shearing and evolving spirals in the outer discs of barred simulations. Kumamoto & Noguchi (2016) excited regular spiral patterns that were particle wakes in the disc forced by a few heavy particles equally spaced around rings, finding additional swing amplification as the wakes lined up. The earlier part of this body of work was summarized in the review by Dobbs & Baba (2014), and all these authors propose that shearing patterns are the fundamental character of spiral arms. We discuss these findings in §1.4.

The Gaia satellite has revealed the local phase space structure of the Milky Way in unprecedented detail (Gaia collaboration: Katz et al., 2018). Sellwood et al. (2019) used these data to try to discriminate among the different theories just described. They calculated the changes to a smooth distribution function that would be caused by a single episode of each of the last three spiral models and compared the predictions with the Gaia data. They concluded that the features in action space seemed more consistent with the transient spiral mode model (Sellwood & Carlberg, 2014) than with either of the other two. Note that Monari et al. (2018) found evidence in the same data for scattering by the Milky Way bar, but some features remained that they could not attribute to the bar. However, Hunt et al. (2018) concluded that other features in the same data were consistent with shearing spiral arms.

Sellwood et al. (2019) also argued that the long-lived spiral mode model advocated by Bertin & Lin (1996) would not cause any pronounced changes, because the ILR, where the largest changes generally occur, is shielded by the $Q$ barrier, no net change is expected at corotation, and the OLR would probably lie too far outside the solar circle to affect the local distribution. Thus if their delicate mechanism for spiral generation were indeed to operate in the Milky Way, it would contribute little to the observd extensive sub-structure in phase space and would implausibly have to survive in the mild disequilibrium state of the disc revealed in the Gaia data.

1.4 Modes or shearing patterns?

The shearing spiral behaviour that is apparent in almost all simulations of cool, isolated discs possesses many aspects of swing-amplification and wakes, as the above cited papers have reported. We have also reported such behaviour: e.g., Fig. 3 of Sellwood & Carlberg (1984) indicated the time evolution of a single 3-arm spiral disturbance that apparently sheared and amplified to an open trailing pattern before decaying as it continued to wind and, in simulations of much improved numerical quality, Sellwood & Carlberg (2014, their Fig. 3) illustrated constantly changing patterns, which is qualitatively similar to the behaviour that all authors report.

However, we have also found (e.g. Sellwood, 1989, 2011; Sellwood & Carlberg, 2014, and later in this paper) that the constantly changing appearance of spiral patterns results from the superposition of a modest number of longer lived waves (see also Quillen et al., 2011). Even though a mode (eq. 1) has a fixed pattern speed at all radii and a constant shape function, $A_{m}(R)$, shearing patterns are the inevitable consequence of superposed modes. All that is required to produce the appearance of a shearing transient spiral (see footnote 2) from two or more superposed patterns of fixed shape is that the modes closer to the centre have higher pattern speed, which is always true.

The proposition that the shearing patterns are the fundamental behaviour is not a satisfactory explanation for the origin of the spirals, because it is unable to answer a key question: how can the spiral amplitudes in simulations be largely independent of the number of particles? No recent paper that argues for this interpretation addresses this issue, yet there is substantial evidence to support it: Sellwood (2011), Sellwood (2012), and Sellwood & Carlberg (2014) all reported similar final amplitudes in experiments in which the number of particles ranged over orders of magnitude. As $N$ is increased, the amplitude of shot noise fluctuations must decrease as $N^{-1/2}$, and some kind of growth mechanism is required to produce final amplitudes that are independent of $N$. Swing-amplification and/or wakes do not lead to indefinite growth, and therefore cannot account for final amplitudes that are independent of $N$, at least in linear theory. We would agree with an argument that non-linear scattering (§1.2) is important, but that plays into our case for unstable modes (Sellwood & Carlberg, 2014, and this paper).

It is likely that most models adopted in simulations are unstable, because they allow some feedback through the centre and/or possess outer edges sharp enough to excite normal modes. Mildly unstable modes will be seeded at low amplitude when $N$ is large and may not cause visible changes for several rotations. Both Sellwood (2011) and Sellwood & Carlberg (2014) present cases in which the first visible features took longer and longer to appear as the number of particles was increased.333Pronouncements of stability after a short period of evolution in large-$N$ models (e.g. D’Onghia et al., 2013) are unlikely to hold in longer integrations. Once the first instability has appeared, we find that activity takes off, as the particle distribution becomes more and more structured by the previous evolution; the new results presented in §3 below show in more detail how a recurrent cycle of modes can occur.

However, it is not neccessary to start from a disc that possesses one or more global instabilities. Toomre (1981) claimed that the half-mass Mestel disc is linearly stable, and this case was studied in simulations by Sellwood (2012). Fig. 1 illustrates the evolution of the $N=5\times 10^{7}$ simulation from his paper; details of the model and the numerical method are given in §2 below. The top panel reproduces the cyan curve from Fig. 2 of Sellwood (2012), which reported at each instant the greatest value of the ratio $\Sigma_{2}(R)/\Sigma_{0}(R)$ within the radius range $1.2<R<12$, where $\Sigma_{2}$ is the amplitude of the bi-symmetric disturbance density. After an initial surge by a factor of a few as each particle created its own wake (§1.2), these $m=2$ features manifested slow secular growth. During this phase, the amplitudes of successive swing-amplified episodes increased slowly because each created mild scratches in the previously smooth disc that enabled partial reflections. The partial mirrors so created at first reflected small fractions of the incident waves, but as the disturbance amplitudes rose with each episode, the new scratches became more reflective. No coherent modes were detectable during this secular growth phase, and the power spectrum, such as shown in the middle panel, was characterized by multiple uncorrelated frequencies arising from whichever noise features were strongest at the time; the corotation radius of each event determines the frequency of the wave that propagates to its ILR where it is absorbed. These findings were consistent with Toomre’s linear theory prediction of global stability, since the secular growth was caused by the non-linear scattering terms that his analysis neglected. The largest amplitude disturbances reached an overdensity of $\sim 2$% by $t\sim 1000$, at which point the reflections became strong enough to cause true unstable modes in this originally stable disc. The amplitude in the top panel of Fig. 1 began to rise more rapidly and the power spectrum, in the bottom panel, that has a different amplitude scale shows much stronger and more coherent waves, which are modes.

Fouvry & Pichon (2015) and Fouvry et al. (2015) successfully developed a quasi-linear theoretical description of this behaviour. De Rijcke et al. (2019) conducted a global mode analysis of this model at $t=1400$, and estimated the frequency of the dominant instability to be $\omega=0.597+0.013i$ in excellent agreement with Sellwood’s own estimate (his paper has factor 10 typo in the growth rate). Although both analyses were of the model denoted 50c, in which the particles coordinates had been scrambled at $t=1400$, the faster wave in the bottom panel of Fig. 1 confirms a mode of similar pattern speed when the evolution is uninterrupted. Fig. 2 shows the logarithmic spiral decomposition of the mode Sellwood (2012) fitted to the data from his model 50c. The leading/trailing bias reflects the trailing appearance of the mode, but the leading side has significant amplitude indicating that it was a mirror, or cavity, mode that operated by an inner reflection off the partial mirror to a leading wave, as Sellwood & Carlberg (2014) described.

Sridhar (2019) calculates the changes to the DF at the principal resonances caused by a transient disturbance, and attempts to show how they lead to further instabilities. However, the correspondence between the modes he calculates and the groove and mirror modes that we have reported is not immediately apparent from his analysis.

Thus we consider that every simulation of a cool, unperturbed disc must inevitably possess, or develop, unstable modes. The unbounded growth of instabilities, until second order terms become important, provides the only viable mechanism to give rise to large-amplitude spirals no matter how large a number of particles is employed. With a recurrence mechanism, such as that described here (§3), the simulations must support multiple unstable modes and their superposition naturally leads to the shearing transient activity that is almost universally reported. Note that externally perturbed models, such as those presented by D’Onghia et al. (2013) and Kumamoto & Noguchi (2016), may not manifest modes if the initial responses develop before any modes could grow to significant amplitude.

1.5 Objective of this paper

The complicated mechanism just described would be required to account for the origin of spirals only in models that do not possess any linear instabilities. The purpose of this paper is to show how a recurring cycle of instabilities can arise in models that do possess a single mild instability at the outset. We argue that such an instability cycle is the origin of spiral activity in all simulations, and hopefully in real galaxies also.

Sellwood & Carlberg (2014) also reported that coherent waves appear and decay at increasing radii and lower frequency over time. It is possible that the later instabilities are independent modes caused by uncorrelated noise fluctuations scratching the DF at larger radii. However, Sellwood & Carlberg (2014) preferred the idea that the decay of one mode created conditions to seed a new instability, although they did not provide a detailed mechanism for recurrence. Here we describe the recurrence mechanism in detail.

In fact, Sellwood & Lin (1989) presented a recurrence mechanism for self-gravitating spiral modes in low mass particle discs orbiting around a central mass. The mechanism they proposed was that scattering at the OLR created a groove in phase space that excited a new instability. We thought it unlikely that the same recurrence mechanism would work in heavy discs, for reasons that we give in §4, but we here demonstrate that this expectation was wrong.

2 Technique

Since the dynamical behaviour of fully self-consistent simulations can be very complicated, we find it fruitful to run simplified simulations that can capture the phenomena we wish to study without them being obscured by unrelated activity. Once understood, it will naturally be important to show that the behaviour persists under more general conditions.

Accordingly, we adopt the razor-thin Mestel disc used in the studies by Zang (1976) and Toomre (1981), which is characterized by a constant circular speed $V_{0}$ at all radii. Note that since $\Omega_{c}=V_{0}/R$ and $\kappa=\surd 2V_{0}/R$, $\Omega_{c}-\kappa/m$ increases indefinitely towards the centre except for $m=1$, when it is negative at all radii. Thus an ILR will intervene to damp all disturbances having $\Omega_{p}>0$ when $m\geq 2$, prohibiting possible cavity modes in a smooth disc, except for $m=1$.

The axisymmetric surface density $\Sigma_{0}(R)=V_{0}^{2}/(2\pi GR)$ would self-consistently yield the appropriate central attraction for centrifugal balance. If the surface density is reduced to fraction $x$ ($0<x\leq 1$), with the removed mass added to a rigid halo to maintain centrifugal balance, then eq. (3) implies that the most vigorously amplified disturbances ($X\simeq 1.5$) will have $m=1.33/x$. Since Zang (1976) had found that the full-mass disc was prone to $m=1$ cavity modes, Toomre (1981) preferred a model with $x=0.5$ in order to avoid lop-sided modes. Note that although swing-amplification will be near its most vigorous for $m=2$ disturbances for this choice of $x$, bisymmetric cavity modes in a smooth disc are disallowed because an ILR would block the feedback loop.

Toomre employed the DF (Toomre, 1977; Binney & Tremaine, 2008)

[TABLE]

where $q=V_{0}^{2}/\sigma_{R}^{2}-1$ and the normalization constant is

[TABLE]

[Sellwood (2012) erroneously omitted the factor $R_{0}^{-(q+1)}$.] Choosing $q=11.44$ yields a Gaussian distribution of velocities such that the $x=0.5$ disc has $Q=1.5$. Toomre further multiplied the DF $f$ by the double taper function

[TABLE]

to create a central cut out centered at $R_{0}$ and an outer taper centered at $R_{1}$, while maintaining the centripetal acceleration $-V_{0}^{2}/R$ everywhere. Setting the taper indices $\nu=4$ and $\mu=5$ yielded an idealized, smooth disc model that Toomre claimed possessed no small amplitude unstable modes. We choose $R_{1}=11.5R_{0}$, and limit the radial extent of the disc by an energy cut-off that eliminates particles having sufficient energy to pass $R=20R_{0}$. Here we adopt units such that $V_{0}=R_{0}=G=1$.

The evolution of one of the simulations of this model reported by Sellwood (2012) was illustrated in Fig. 1 above. The simulations in the present paper use the same code and disc model, have the same number of particles (see Table 1), but are integrated to $t=1000$ only and so are unaffected by the later rapid growth phase.

The particles in our simulations are constrained to move in a plane over a 2D polar mesh at which the accelerations are calculated, and interpolated to the position of each particle. In most cases, disturbance forces are restricted to a single sectoral harmonic, $m=2$ or $m=3$. Since the gravitational field is a convolution of the mass density with a Green function that is most efficiently computed by Fourier transforms, it is easy to restrict the sectoral hamonics that contribute to the field when using a polar grid. A full description of our numerical procedures is given in the on-line manual (Sellwood, 2014b) and the code itself is available for download.

Table 1 gives the values of the numerical parameters for most simulations reported here, and we note when they are varied in a few cases. Sellwood (2012) reported, and we have reconfirmed in this study, that all our results are insensitive to reasonable changes to grid resoution, time step and zones, and number of particles. Changes to the softening length do affect the frequencies of the modes, but not the qualititative behaviour.

To set up each model we draw particles from the tapered DF (eqs. 4 and 6) as described in the appendix of Debattista & Sellwood (2000). However, in some experiments we modified the distribution of selected particles in order to introduce a single additional feature into the DF, as set out in §3.1.

We also restarted some simulations after “scrambling” a copy of the particle distribution. By this we mean that we changed $(R,\phi,v_{R},v_{\phi})\rightarrow(R,\phi^{\prime},v_{R},v_{\phi})$ for every particle, with $\phi^{\prime}$ being chosen at random from a distribution that is uniform in 0 to $2\pi$. Clearly scrambling resets the amplitudes of all non-axisymmetric disturbances back to the shot noise level of the initial disc, while preserving the radius and both velocity components (in polar coordinates), so that any features in the action distribution of the particles that had been introduced during prior evolution would be preserved. This is therefore not equivalent to a fresh start with a different random seed.

Note that the radial action, $J_{R}\equiv\oint\dot{R}dR/(2\pi)$ (Binney & Tremaine, 2008), has dimensions of angular momentum, is zero for a circular orbit, and increases with the eccentricity of the orbit. In an axisymmetric potential, the angular momentum, $L_{z}$, is the other action for orbits confined to a plane.

3 Results

In this section we study the effect of introducing two distinct features into the otherwise unperturbed disc, both designed to excite an initial mode of a specific type, namely: a groove mode (Sellwood & Lin, 1989; Sellwood & Kahn, 1991) and an outer edge mode (Toomre, 1981; Papaloizou & Lin, 1989). In both cases, we follow its evolution and study the further modes that are excited. We show for each case that the second mode is a true instability of the modified disc by erasing all non-axisymmetric structure at the time the second mode began to grow, as just described, and demonstrating that the scrambled particle distribution possessed a very similar instability to that which developed in the continued original run. Furthermore we identify the scattering feature in action space that was responsible for the second instability, as reported below in §3.1.3.

3.1 Initial groove mode

We create a deficiency of low-$J_{R}$ particles over a narrow inclined range in $(L_{z},J_{R})$ space by giving some particles a larger $J_{R}$, as shown in Fig. 3. The slope of this feature, $-1/2$, approximately traces the locus of the ILR of a bisymmetric disturbance having a pattern speed $\Omega_{n}=0.1$, although the disturbance was, in fact, purely hypothetical. The ILR of a disturbance having this $\Omega_{n}$ would lie at $L_{z}=2.93$ for circular orbits, for which $J_{R}=0$.444Note that all the major resonances with any pattern speed and any angular periodicity in the self-similar Mestel disc lie on lines of slope $\simeq-1/2$ in the space of these actions, for $J_{R}\ll L_{z}$.

We describe this feature as a groove, since it is a deficiency at low $J_{R}$ over a narrow range in $L_{z}$. Since orbits in this self-similar disc have typical epicyclic radii of $\sigma_{R}/\kappa\simeq 0.2R_{g}$ about the guiding centre radius, this narrow groove in action space causes no noticeable change to the surface density profile.

Our, perhaps somewhat clumsy, procedure for shifting particles was motivated by our observations of resonance scattering in previous simulations. We calculated the frequency distance from the resonance $\delta\omega=2\Omega_{\phi}-\Omega_{R}-2\Omega_{n}$ for each particle, where $\Omega_{\phi}$ and $\Omega_{R}$ are respectively the angular frequencies of the guiding centre and radial oscillation for orbits of arbitrary eccentricity; note that these frequencies tend to the familiar quantities $\Omega_{\phi}\rightarrow\Omega_{c}$ and $\Omega_{R}\rightarrow\kappa$ as the eccentricity $\rightarrow 0$. We then select particles that have $g^{2}=1-(\delta\omega/w_{\omega}\Omega_{n})^{2}>0$ as candidates to be scattered. The relative frequency width, $w_{\omega}=0.06$, of this feature is deliberately rather narrow. Of these, we further select only those particles whose energy of random motion $E_{r}=E-E_{c}(L_{z})<gE_{r,{\rm lim}}$, where the groove extends to the limiting random energy $E_{r,{\rm lim}}=0.025$ at $g=1$, i.e. along the resonance locus. The probability that these selected particles are given larger random energy decreases from $d$ to zero as $E_{r}$ increases from zero to this maximum, where the fractional depth of the groove $d=0.25$. The new $E_{r}=E_{r,{\rm lim}}(1+3p)$ where $p$ is the solution of $\exp(-p^{2})=s$, with $s$ being randomly drawn from a distribution uniform in 0–1, and finally $\Delta L_{z}=\Delta E/\Omega_{n}$.

A simulation that employed this modified DF, model G, supported a groove instability at first. The red curve in the top panel of Fig. 4 presents the time evolution of the ratio $\langle\Sigma_{2}(R)/\Sigma_{0}(R)\rangle$ averaged over the radius range $2<R<8$. Since this is a different measure from that presented in Fig. 1, the dotted curve indicates this measure for the same model without the groove. The green and blue lines are explained below. Power spectra over two periods of growth are shown in the lower two panels. The middle panel is for the same time period as the middle panel in Fig. 1, although the grey scales differ, and the contrast in the appearance is quite striking.

3.1.1 First mode

The disturbance in model G grew slowly until $t\sim 400$, and the power spectrum over this time interval in the middle panel of Fig. 4 is dominated by a coherent disturbance of angular frequency $2\Omega_{p}\approx 0.65$ that is mostly localized between the ILR and corotation. A couple of other frequencies that are probably due to swing-amplified noise are also faintly visible, and these uncorellated mild disturbances are responsible for the small, but rapid, amplitude fluctuations in the red line in the top panel.

We employed the mode-fitting procedure devised by Sellwood & Athanassoula (1986) to the data from model G over the time interval $100\leq t\leq 400$ and give its estimated eigenfrequency in line 1 of Table 2. The shape of the fitted mode, with its principal resonances for circular orbits ($J_{R}=0$) marked by the circles, is shown in the upper panel of Fig 5, and the decomposition of the mode into logarithmic spirals is in the lower panel.

The mechanism for the groove instability was fully explained by Sellwood & Kahn (1991) and outlined above in §1.1. It is not a cavity-type mode, and the frequency is determined by the gradients in angular momentum density for orbits of small radial action. The density changes caused by non-axisymmetric disturbances at the groove “edges” excite a supporting response from the surrounding disc, making it a large-scale spiral mode. In principle, the circular angular frequency of the groove centre should be the pattern speed of the mode, but geometric factors in low-$m$ modes and finite growth rate cause the corotation radius to lie slightly farther out. The shape of the mode is determined by the shear rate and surface density profile of the disc, so that in the present self-similar disc, groove modes should have closely similar shapes and pitch angles and their spatial scale will vary with the corotation radius.

Evidence that this mode is a groove instability (Sellwood & Lin, 1989) is two-fold. First, corotation (line 1 of Table 2) lies close to, but just outside (Sellwood & Kahn, 1991), the groove centre at $L_{z}=2.93$. Second, the leading component in the mode transform (bottom panel of Fig. 5) is much weaker than that in Fig. 2 because groove modes (§1.1) do not operate by a feedback cycle of combined leading and trailing waves.

Our choice of a very narrow frequency width for the groove ensured a low growth rate for the mode (Sellwood & Kahn, 1991). We found in other simulations, not described in detail here, that wider grooves excited a mode that grew more vigorously to higher amplitude. We preferred to study an instability that caused mild non-linear changes in order to make it easier to follow further activity that was a consequence of the initial mode.

3.1.2 Second mode

After the slow growth to $t=400$, the disturbance amplitude in model G, shown by the red line in the top panel of Fig. 4, rose more rapidly over the time range $400<t<700$. The power spectrum over this time range in the bottom panel revealed a new coherent wave having a frequency $\sim 0.45$. As the mode had clearly saturated by the end of this period, we fitted a mode over the time interval $200<t<500$, obtaining good fits with eigenfrequencies in the range given in line 2 of Table 2.

In order to verify that this second mode was a true linear instability of the modified disc, we made a copy of the particle distribution at $t=400$ and evolved the scrambled (see §2) particle distribution in a separate simulation, designated GR. The amplitude evolution in the simulation starting from this scrambled disc is shown by the green line in the top panel of Fig. 4. The very first value was measured from the particles before they were scrambled, after which the amplitude rose out of the shot noise steadily to $t=700$ when the run was stopped.

We found that the evolution of the scrambled disc, model GR, was again dominated by an exponentially growing disturbance of fixed pattern speed. The eigenfrequency fitted to the data from this run, given in line 3 of Table 2, is within the uncertainties the same as that in line 2 – the second mode found in model G. A second, milder instability was also present that had the fitted frequency given in line 4. At first we thought it could be the original groove mode, although the pattern speed is distinctly higher, but we identify its progeny at the end of §3.1.3.

3.1.3 Cause of the second mode

The second more vigorous mode, which was not present in the early evolution of the same model G (middle panel of Fig. 4), was caused by changes to the DF produced by the first mode. Fig. 7 reveals the changes to the density of particles in action space between times 0 and 400. The colour scale in this figure shows where the density of particles at $t=400$ is both greater and less than that at $t=0$, and reveals several features. The dashed line marks the locus of the corotation resonance, and the dotted lines the Lindblad resonances, for the fitted frequency of the first mode in this simulation. These lines indicate that the most prominent feature is the result of ILR scattering and the faint pair of features at larger $L_{z}$ are associated with the OLR. A paired deficiency and excess have appeared near the location of the original groove at $L_{z}\simeq 2.93$, reflecting the changes that must have been caused by the instability. As such changes occur near corotation, marked by the dashed line, they should not cause any significant change to $J_{R}$ (Sellwood & Binney, 2002), which is largely true. Since the groove instability is so mild, Fig. 7 also shows the effects of scattering by weak noise features that generally had lower frequencies than the excited mode.

Note the differences in the nature of the changes caused by scattering at the ILR from those at the OLR. In action space, the scattering vector at a resonance has slope $\Delta J_{R}/\Delta Lz=l/m$ (Sellwood & Binney, 2002), where $m$ is the angular multiplicity of the pattern and $l=0,\;\mp 1$ at respectively corotation and the inner and outer Lindblad resonances. Formally the slope $l/m$ is exact only as $J_{R}\rightarrow 0$ but we find it is an excellent approximation over the range of $J_{R}$ of interest here. As noted before (Sellwood, 2012; Sellwood & Carlberg, 2014; Sellwood et al., 2019), the scattering vector slope, $-1/2$, is almost perfectly aligned with the ILR locus for $m=2$ waves, implying that particles stay on resonance as they gain $J_{R}$, allowing large changes to build up. Thus ILR scattering moves particles along the resonance locus, creating a deficiency of particles at low $J_{R}$, and an excess for larger $J_{R}$, just as we created artificially at the start of the first simulation (Fig. 3). On the other hand, scattering vectors at the OLR have positive slope $+1/m$, while the resonance locus again has slope $-1/2$. Thus we see particles are scattered across the resonance from a sloping deficiency to an almost parallel line where there is a slight excess. It is perhaps interesting that the lines are approximately parallel, indicating that the action changes are roughly independent of the initial $J_{R}$. More importantly, $J_{R}$ is an intrinsically positive quantity, and therefore Lindblad resonance scattering always creates a deficiency at $J_{R}\ga 0$.

We suspected that one of the features in Fig. 7 was responsible for the new instability in model G. In order to determine which, we extracted all the particles in one of the features and spliced them into a pristine undisturbed particle distribution that we then evolved. In detail, we began the splicing procedure by creating a fresh set of particles from the doubly tapered Toomre-Zang DF described in §2 that lacked any additional features; recall that Toomre (1981) had predicted, and Sellwood (2012) had confirmed, that this model has no coherent, small amplitude instabilties. We then extracted all the particles from the evolved distribution at $t=400$ that lay within a trapezium in action space of width $\Delta L_{z}=1$ that was defined by $|L_{z}+0.5J_{R}-L_{z,0}|<0.5$, for some value of $L_{z,0}$. We spliced the extracted particles into the undisturbed set by substituting the values of $(x,y,v_{x},v_{y})$ of these selected particles from the evolved set for the values in the undisturbed set of particles that lay within the same trapezium in action space.555Since the central attraction, $-V_{0}^{2}/R$, is an unchanging function throughout these experiments, we do not have to worry about possible changes to the axisymmetric potential. This process required attention to two points in order to obtain a near seamless splice. First, almost all the particles in the evolved model had changed their $L_{z}$ values during the evolution to $t=400$, making it impossible to just replace particle $n$ from the pristine set with particle $n$ from the evolved set, even though they were both originally drawn from the same DF by the same algorithm. We therefore had to compute $J^{\prime}=L_{z}+0.5J_{R}$ for every particle in both sets and rank order them to identify those in the desired range. Second, the large numbers of particles in the desired range of $J^{\prime}$, typically some 4.5 million, was not precisely the same in the two sets, so that we had to skip the occasional particle, usually one in several hundred, in the larger set to ensure equal numbers that spanned the entire range. Fig. 7 shows the difference between the spliced and undisturbed sets for $L_{z,0}=5$. It is clear that the OLR scattering feature from Fig. 7 has been satisfactorily isolated from all other features in action space at that time, without introducing sharp edges to the splice.

The time evolution of the disturbance density in a simulation started from the spliced distribution, model GSOLR, is indicated by the blue line in the top panel of Fig. 4. The mild non-axisymmetry of the spliced particles at $t=400$, was erased by scrambling before the evolution commenced, and the disturbance density grew out of the noise, closely following the evolution from the simply scrambled disc, model GR, traced by the green line.666Other modes are excited after $t\sim 700$, and interference between disturbances rotating at different rates can cause large temporary decreases in the measured amplitude. The eigenfrequency fitted to the data from this spliced run is given in line 5 of Table 2 and is, again within the uncertainties, the same as those in lines 2 and 3, the second mode in model G and the mode in the scrambled disc, model GR. This evidence conclusively shows that the OLR scattering feature of the first mode in model G was responsible for creating the second instability.

Note that the corotation resonance for near circular orbits of the modes given in lines 2, 3, and 5 of Table 2, all lie near the OLR feature created by the first mode (Fig. 7) strongly suggesting that the second mode was also a groove instability. It seems likely that the deficit of particles near $J_{R}\ga 0$ was the cause.

We were quite surprised to find it was the OLR scattering feature that was responsible for the second instability for the reasons set out in the discussion below (§4). Yet our experiment with the isolated weak scattering feature of the OLR leaves no doubt that it did indeed provoke a new instability that was more vigorous than the original mode.

The near coincidence of the OLR for the first mode (line 1 of Table 2) with CR for the second (lines 2, 3 and 5 of Table 2) suggests non-linear mode coupling. Resonance scattering to create the new groove is a non-linear effect, but we have strong evidence against the hypothesis of coupling at resonances between large amplitude waves that was first proposed by Tagger et al. (1987). A technical reason is that their theory requires the interaction of three waves, with the third wave having rotational symmetry that is the sum or difference of the two other waves, and our simulation in which disturbance forces were restricted to $m=2$ could not support disturbances having $m=0$ or 4. However, the more compelling argument is that scrambling the particles erased all non-axisymmetric features in the density distribution, which need to have significant amplitude for the mode coupling mechanism to work. Instead, our simulations provide clear evidence that the second mode is a true, linear instability of the modified disc at $t=400$.

We tried other experiments for which we chose other splice centres, but do not illustrate the results. Choosing $L_{z,0}=3$, model GSCR supported a mild groove mode with frequency given in line 6 of Table 2; the deficiency of the original DF had been shifted to lower $L_{z}$ by the evolution of the first groove mode, causing the slightly higher pattern speed than that reported in line 1. It is satisfying to note that this instability approximately matches that of the second mode in the scrambled model GR, reported in line 4 of the table. Simulations of particle distributions that resulted from splicing with $L_{z,0}=1$ or from $L_{z,0}=2$ appeared not to support any coherent instability. The strong ILR scattering feature from the first mode lies within the inner taper, and is therefore unable to excite a significant instability because the tapered surounding disc would provide only a weak supporting response of (Sellwood & Kahn, 1991).

3.2 Initial edge mode

All we needed to do to excite an edge mode was to increase the value of $\mu$ for the outer index taper (eq. 6). Setting $\mu=20$ created a vigorous instability that saturated at large amplitude, strongly distorting a large fraction of the disc. After some experimentation, we found that setting $\mu=12$, model E, yielded an outer edge mode with a moderate growth rate that enabled us to understand the subsequent evolution. Note that for the sequence of runs reported in this section, in which the initial disturbances were in the outer disc where the grid becomes coarser, we used a 4 times finer grid ($404\times 512$). We also restricted disturbance forces to $m=3$, rather than $m=2$, in order that the modes have smaller radial extent, but the other numerical parameters were unchanged.

The evolution of the non-axisymmetric density of model E is shown in Fig. 8. A 3-arm spiral grows slowly in the outer disc until $t\sim 600$, after which additional patterns appear closer to the disc centre. The amplitude evolution of the $m=3$ disturbance density over the radius range $3<R<12$ and power spectra over two periods of growth are shown in Fig. 9.

3.2.1 First mode

The evolution of model E to $t=600$, is dominated (middle panel of Fig. 9) by a single, exponentially growing disturbance of frequency $3\Omega_{p}\ga 0.2$, that extends from its ILR to OLR, with corotation perhaps near $R\sim 13$. The mode fitting software yielded an eigenfrequency given in line 7 of Table 2 and the fitted mode and its logarithmic spiral transform are presented in Fig. 10.

It has the hallmarks of an edge mode. As expected (Toomre, 1989), $R_{\rm CR}$ (line 7 of Table 2) lies outside the radius of the steepest gradient, which is at $L_{z}=11.5$, the centre of the outer taper. Again, the logarithmic spiral transform of the mode, lower panel of Fig. 10, has no significant leading component as expected for an edge mode, since the mechanism does not require feedback through leading waves (see §1.1).

Fig. 11 displays the changes to the density of particles in action space between $t=0$ and $t=600$ in model E, and the dashed and dotted lines mark the loci of the principal resonances for the measured pattern speed of the original edge mode. It reveals scattering features at corotation and the ILR, but there are no significant changes at the OLR because there is very little disc left at those radii, and angular momentum transport by this edge mode is from the ILR to CR only. Since the power spectra (Fig. 9) and the fitted mode shape (Fig. 10) were drawn using the relative overdensity at each radius, they exaggerate the apparent amplitude of the mode in the very outer part of the disc. The very low frequency features in the bottom panel of Fig. 9 have a small absolute amplitude, and probably arise from shot noise due to the low particle density near the disc edge.

As noted above, scattering vectors in action space have slopes $=l/m$. Thus the near perfect alignment of the scattering vector with the locus of the ILR occurs only for $m=2$ disturbances, and therefore the excess and deficiency in Fig. 11 at the ILR are slightly offset from the resonance locus for this $m=3$ pattern, and the groove carved by the mode at $L_{z}=6.9$ is fractionally exterior to its ILR.

3.2.2 An inwardly propagating cascade of groove modes

The power spectrum of the later part of this simulation, in the bottom panel of Fig. 9, indicates several additional coherent waves, in addition to the continuing presence of the first mode, which had not even finished growing by $t=600$.

We have fitted two modes to the data from the simulation over the period $500<t<780$, recovering both the edge mode and the second mode, which has an estimated frequency given in line 8 of Table 2. Corotation for this second mode, also given in line 8, is somewhat outside the groove centre carved by the original edge mode, as is usual groove mode (Sellwood & Kahn, 1991).

The green line in Fig. 9 shows the amplitude evolution we obtain in a new simulation, model ER, started from the scrambled particles of model E at $t=600$. Once again we find the perturbation amplitude grows rapidly with the estimated frequency given in line 9 of Table 2, which is in reasonable agreement with the frequency of the second mode given in line 8. Note the second mode outgrows the original edge mode, which is also detectable in the scrambled simulation. The true uncertainties in the measured frequencies and resonance radii in this case are probably greater then the formal values given in the table because we were attempting to fit two vigorously growing modes over a short time interval.

The numerical evidence just presented implies that the second mode in model E is a groove mode seeded by ILR scattering by the first mode (Fig. 11). Once again, the near coincidence of the ILR of mode 1 (line 7 of Table 2) is CR of its daughter (line 8) is not evidence for non-linear mode coupling, for the same reasons as given above. We obtained the daughter mode (line 9) even when we eliminated all pre-existing waves by scrambling the disc particles, and we have strong evidence that the new mode is a linear instability caused by scattering at the ILR of the first mode.

The power spectrum in the bottom panel of Fig. 9 strongly suggests a cascade of groove instabilities to ever higher frequency, with corotation of each subsequent instability lying at the ILR of the previous. Since the initial outer edge mode had a low growth rate, the inner disc had to wait for a succession of destabilising grooves to be carved before each mode could develop. This delay is all the more remarkable because the dynamical clock runs faster at smaller radii. The long absence of coherent disturbances in the inner disc could therefore happen only if the single instability the disc possesses is at its outer edge.

3.3 Changes of angular symmetry

Disturbance forces in the two sets of experiments described above were each confined to a single sectoral harmonic: $m=2$ for the initial groove mode (§3.1) and $m=3$ for the initial edge mode (§3.2). Furthermore, the simulations with scrambled or spliced particle sets were also restricted to the same symmetry as in the parent simulation. Here we show that the features in action space created by the first mode can also give rise to instabilities having other rotational symmetries.

3.3.1 From 3 to 2

The initial $m=3$ edge mode of model E presented in §3.2 scattered particles to produce the features shown in Fig. 11 by $t=600$. We have already demonstrated that the ILR scattering feature causes a new $m=3$ groove instability that we identified both in the continued run, model E, and in a separate simulation, model ER, started from a scrambled copy of the particles at $t=600$. Here we report another simulation, model ER2, with the same scrambled copy of the particles at $t=600$, but with disturbance forces restricted to $m=2$ instead.

The disturbance amplitude exponentiates rapidly, as shown in the top panel of Fig.12, and the power spectrum in the bottom panel reveals two simultaneous instabilities. Our mode fitting software finds two modes, a groove instability given in line 10 of Table 2 and an edge mode having a frequency given in line 11.

As noted above (§3.2.2), the $m=3$ edge mode in model E was still growing at $t=600$, and it is therefore hardly surprising to find that the edge of this disc remains steep enough to excite a new edge mode this time. The new $m=2$ mode, line 11 of Table 2, has a higher growth rate than the $m=3$ edge mode because the swing amplifier is at full strength in this case, $X=2$, whereas $X=4/3$ for $m=3$. Because the instability is more vigorous, corotation is at the larger radius, as also indicated in the Table.

The $m=2$ groove mode in model ER2 also has a higher growth rate, line 10, than at $m=3$ instability from the same file of particles (line 9 of Table 2), consistent with more vigorous swing amplification, and corotation is still farther from the groove created by ILR scattering by the original $m=3$ edge mode at $R\simeq 6.9$, as expected for a more vigorous and larger-scale groove instability (Sellwood & Kahn, 1991).

3.3.2 From 2 to 3

In this next case, model GS3OLR, we start from the spliced set of particles shown in Fig. 7, but this time restrict perturbation forces to the $m=3$ sectoral harmonic. Our estimated frequency of the dominant instability is given in line 12 of Table 2, and corotation is again almost exactly at the location of the groove in action space, as for the modes in lines 2, 3 and 5. Notice that the growth rate is a little lower than for the corresponding $m=2$ case (line 5 of Table 2), consistent with expectations since swing amplification is less vigorous for $m=3$ than for $m=2$.

These last two experiments have demonstrated that a groove created by Lindblad resonance scattering from a pattern of one sectoral harmonic can excite an instability of a different $m$.

3.4 A very slightly more general case

Perturbation forces in all the preceding simulations were restricted to a single sectoral harmonic, which we now relax. We start a simulation, model GU, from the particle distribution used in the groove mode case, shown in Fig. 3, which we know from model G will excite at least a mild $m=2$ instability. However, perturbation forces in this new experiment now include $0\leq m\leq 8$ terms, with the exception that $m=1$ is excluded to avoid possible unbalanced forces from the fixed rigid halo component.

Multiple patterns develop over time, as shown in Fig. 13, with perhaps the first feature to appear being an $m=3$ mode. Fig. 14 presents power spectra from model GU over the same time interval for both $m=2$ and $m=3$, revealing a number of coherent waves with pattern speeds independent of radius. A number of strong scattering features in action space are visible at $t=300$ in Fig. 15 (top). The bottom panel shows that changes become so intricate in the later evolution that we are no longer able to isolate the effects of individual modes, but there seems little doubt that the disc is being repeatedly destabilized by scattering at resonances of successive disturbances, each of which removes particles from near-circular orbits over narrow ranges of $L_{z}$.

4 Discussion

The simulations of Sellwood & Carlberg (2014) manifested an apparently outwardly propagating cascade of spiral modes to ever larger radii over time. Although they presented evidence of resonance scattering and argued that the changes to the DF excited subsequent modes, they did not identify the mechanism. Here we have presented compelling evidence that the recurrence mechanism identified long ago by Sellwood & Lin (1989) for $m=4$ spirals a low-mass disc with a Keplerian rotation curve can, in fact, operate in massive discs having flat rotation curves. We doubted that this would be viable for three principal reasons that, in hindsight, no longer seem particularly cogent.

1. 1

Probably the strongest reason was that the radial extent of spiral modes, which are expected to extend from the ILR to the OLR (Lin & Shu, 1964), is much greater for $m=2$ waves in a disc having a flat rotation curve. Both the lower $m$ value and the shape of the rotation curve cause the ratio $R_{\rm OLR}/R_{\rm ILR}$, which was 1.67 for $m=4$ waves in a Kepler potential, to rise to $\sim 5.8$ for $m=2$ patterns in a flat rotation curve. If the patterns were required not to overlap, it would be hard to fit many into any disc of reasonable extent. But in fact, for a cascade of outwardly propagating $m=2$ modes, successive disturbances have radii in the ratio $R_{\rm OLR}/R_{\rm CR}\simeq 1.7$. Even for an inwardly propagating $m=2$ cascade initiated by an edge mode, say, the ratio is $R_{\rm CR}/R_{\rm ILR}\simeq 3.4$. These smaller ratios, which also allow extensive overlap of successive modes, would be even further reduced for larger values of $m$. 2. 2

It is well known that $m=2$ spiral disturbances have much greater amplitude inside corotation than outside – examples include the famous dust-to-ashes figure from Toomre (1981), as well as the power spectra and modes presented here, e.g. the middle panel of Fig. 4. However, the fitted mode in that case (Fig. 5) does have a weak spiral that extends past the OLR circle. This general disparity suggested that scattering at the OLR would not be as important as at the ILR. Since self-excited modes in a disc must conserve angular momentum (there is no externally applied torque), the gains and losses must balance and the angular momentum emitted at the ILR must be absorbed at other resonances. With the exception of edge modes, where corotation takes up almost all the angular momentum lost at the ILR (Fig. 11), the OLR must absorb a significant amount. 3. 3

The third reason we expected that OLR scattering would be less consequential than at the ILR has already been stated. The close alignment of the ILR locus and the scattering vector causes particles to stay on resonance as they are scattered allowing large changes to build up. This does not happen at the OLR, where scattered particles are quickly moved off resonance and the angular momentum gain must therefore be shared among many more particles. Again this is true and consistent with our results presented in Fig. 7, which led us to doubt that the feeble OLR scattering seen there could excite the second mode. However, it seems that the depopulation of the low-$J_{R}$ region, as many particles are shifted away from near-circular orbits, is the dominant source of excitation of the next mode.

These arguments make it clear that we were wrong to doubt the possibility of a cascade of groove instabilities, running either radially outward or inward, could occur in a massive disc having a flat rotation curve. The evidence presented in §3 demonstrates that cascades of groove modes, with each causing changes to excite the next, arise naturally.

5 Conclusions

In this paper, we have presented a new part of our picture that spiral patterns result from true instabilities in galaxy discs. We have demonstrated that scattering at either Lindblad resonance by any one wave carves a new groove in the disc that excites a further mode. We presented examples of both inwardly and outwardly propagating cascades of instabilities in simulations restricted to a single sectoral harmonic. We also showed that a groove created by one wave of a certain angular symmetry can excite a mode of different angular symmetry. Thus, in more general simulations that include disturbance forces from multiple sectoral harmonics, such as are inevitable when tree codes are employed, it quickly becomes very difficult to follow the causal chain from one wave to the next.

The mirror mode mechanism, proposed by Sellwood & Carlberg (2014), seemed to be responsible for the first true mode in models 50 of Sellwood (2012), as evidenced by Fig. 2 shown here. However, we suspect it is important only in models that lack strong instabilities. Once a recurrent cycle of groove modes becomes established, mild variations of impedance to cause partial reflections of travelling waves must become physically less important as scattering features in phase space increase in both the strength and number – see e.g. the bottom panel of our Fig. 15 and Fig. 6 of Sellwood & Carlberg (2014), where 3D motion was allowed. The recurrent cycle of groove modes reported here is then far more likely to be the main driver of evolution.

The additional evidence we have presented here strengthens our argument that spirals in simulations are manifestations of global instabilities. The modes in our simulations are each detectable as coherent waves of fixed frequency for a number of rotations, but the visual pattern of spirals changes on a much shorter time scale. The superposition of two or more modes inevitably leads to the apparently shearing spiral features that are widely reported. But only unbounded growth of true instabilities in the linear regime can account for the final amplitudes being independent of the number of particles, as we have repeatedly shown in our previous work (Sellwood, 2011, 2012; Sellwood & Carlberg, 2014). None of the many papers arguing for the fundamental nature of shearing features has so far presented any mechanism that could deliver unbounded growth. Swing amplification and wakes, which are invoked in these papers, undoubtedly enhance collective features, but only by a fixed factor in linear theory (Julian & Toomre, 1966; Toomre, 1981). D’Onghia et al. (2013) describe non-linear consequences of responses to the mass clumps they introduced into their disc models, which do appear to cause more persisting non-axisymmetric structures, but it is unclear as yet whether they could account for indefinite growth. The non-linear scattering at resonances, that is a crucial part of our picture, readily excites new instabilities, leading to a robust cycle of recurring modes that accounts for the kind of behaviour that is almost universally observed in simulations of disc galaxies.

Naturally, resonant scattering causes irreversible changes, and the consequent secular increase in random motion would, if unchecked, degrade the ability of the disc to support continuing collective instabilities, and spiral activity must fade. But Sellwood & Carlberg (1984) were the first to show that allowing for dissipation and star formation in a reasonable gas fraction could allow spiral activity to continue “indefinitely,” an idea that has been supported in much subsequent work.

In their analysis of the local Gaia data, Sellwood et al. (2019) argued that an action-space representation of local phase space revealed evidence for Lindblad resonance scattering, that was not all attributable to the bar (Monari et al., 2018). This evidence suggests that the Milky Way has supported global spiral modes of the type we describe in our simulations, and we are reasonably confident that spirals in the old stellar discs of many other galaxies also result from such instabilities.

Acknowledgements

We thank an anonymous referee for thoughtful comments on the paper, and also Elena D’Onghia and especially Daisuke Kawata for helpful conversations during the Gaia19 program at KITP, where this paper was advanced significantly. KITP is supported in part by NSF grant PHY-1748958. JAS also acknowledges the hospitality of Steward Observatory.