ABSTRACT

Using the linearized Boltzmann equation, we investigate how grooves carved in the phase space of a half-mass Mestel disc can trigger the vigorous growth of two-armed spiral eigenmodes. Such grooves result from the self-induced dynamics of a disc subject to finite-N shot noise, as swing-amplified noise patterns push stars towards lower angular momentum orbits at their inner Lindblad radius. Supplementing the linear theory with analytical arguments, we show that the dominant spiral mode is a cavity mode with reflections off the forbidden region around corotation and off the deepest groove. Other subdominant modes are identified as groove modes. We provide evidence that the depletion of near-circular orbits, and not the addition of radial orbits, is the crucial physical ingredient that causes these new eigenmodes. Thus, it is possible for an isolated, linearly stable stellar disc to spontaneously become linearly unstable via the self-induced formation of phase-space grooves through finite-N dynamics. These results may help explain the growth and maintenance of spiral patterns in real disc galaxies.

1 INTRODUCTION

Despite the many tantalizing hints gleaned from detailed observations (Meidt, Rand & Merrifield 2009; Salo et al. 2010; Beckman et al. 2018), analytical calculations (Goldreich & Lynden-Bell 1965; Julian & Toomre 1966; Lynden-Bell & Kalnajs 1972; Mark 1976; Omurkanov & Polyachenko 2014; De Rijcke & Voulis 2016), and numerical simulations (Sellwood 2011; D’Onghia, Vogelsberger & Hernquist 2013; Sellwood & Carlberg 2014; Saha & Elmegreen 2016; Semczuk, Łokas & del Pino 2017), the cause(s) and the life expectancy of spiral structures in disc galaxies are still uncertain. Even restricting the general problem to that of the gravitational dynamics of an axially symmetric, razor-thin stellar disc still leaves a wealth of dynamical processes to be explored and understood. Given the observed variety of the spiral galaxy zoo, with specimens ranging from the beautifully symmetric grand-design spirals to the patchy and chaotic flocculent disc galaxies, a unique explanation may even seem elusive. This is unfortunate, given the apparent importance of spiral structures for the secular evolution of the galaxies that host them (Zhang 1998; Sellwood & Binney 2002; Roškar et al. 2008, 2012; Sellwood 2014; Fouvry et al. 2015; Vaghmare et al. 2015; Daniel & Wyse 2018) and for regulating the rate and location of their star formation (Seigar & James 2002; Aramyan et al. 2016; Hart et al. 2017).

Early on in the history of this topic, it was understood that stellar discs can respond vigorously to disturbances whose wavelength is much smaller than the host disc. Such localized perturbations can be shown to be strongly amplified as they shear from leading to trailing (Goldreich & Lynden-Bell 1965; Julian & Toomre 1966; Toomre 1981). This ‘swing amplification’ mechanism appears to be a plausible explanation for the ragged appearance of the flocculent disc galaxies. It cannot, at first sight, explain the open, long-wavelength patterns observed in grand-design systems. In the original quasi-steady-state theory for grand-design spirals (Lin & Shu 1964), the gravitational maintenance of an imposed spiral pattern was expounded. However, this theory lacked a generating mechanism for the spiral patterns.

Barring those cases where a spiral pattern can be linked to an external pertuber’s tidal forces, an internal cause must be searched. It was soon realized that such patterns could originate from growing eigenmodes in linearly unstable stellar discs (Toomre 1964; Hunter 1965; Kalnajs 1977; Vauterin & Dejonghe 1996). If a galaxy is linearly unstable, even the smallest perturbation is sufficient to trigger its eigenmodes (Toomre 1964; Hunter 1965; Kalnajs 1977; Vauterin & Dejonghe 1996; Pichon & Cannon 1997) and begin its evolution towards higher entropy states (Lynden-Bell & Kalnajs 1972). This leads to secular radial migration of stars and disc heating (Zhang 1999; Sellwood & Binney 2002; Roškar et al. 2008, 2012).

In this paper, we show that such linearly unstable states need not be exceptional and that they are, in fact, a natural outcome of the self-induced dynamics of a finite-N stellar disc (with N ≫ 1 the number of stars). This point was already argued by Sellwood (2012) and Fouvry et al. (2015), who investigated the evolution of a linearly stable half-mass Mestel disc (Mestel 1963; Toomre 1981) using N-body simulations and the integration of the inhomogeneous Balescu–Lenard equation (Heyvaerts 2010; Chavanis 2012), respectively. The latter equation, which is valid at the order 1/N, accounts for the self-driven orbital secular diffusion of a self-gravitating system induced by the intrinsic shot noise due to its discreteness. This diffusion results in the formation of a resonant groove in phase space along a specific resonant direction. In other words, the finite-N dynamics of an initially linearly stable stellar disc can lead to its spontaneous destabilization and the subsequent growth of spiral-shaped eigenmodes. The origin of these modes lies in the gravitational amplification of waves trapped between the groove and corotation.

This paper is organized as follows. In Section 2, we introduce the half-mass Mestel disc model, both without and with phase-space grooves. This is followed by our linear mode analysis of this model in Section 3. We discuss the significance of our results in Section 4 and we conclude with Section 5. The linear stability tool we employ is shortly described in Appendix  A.

2 THE HALF-MASS MESTEL DISC

2.1 The linearly stable half-mass Mestel disc

The Mestel disc (Mestel 1963) has a surface mass density given by
(1)
which self-consistently generates a gravitational field with a binding potential
(2)
Here, v0 is the constant circular velocity of this stellar disc model and r0 is a scale length. Clearly, |$v_0^2 = 2\pi G \Sigma _0 r_0$| for consistency (this relation defines the density scale factor Σ0). A self-consistent distribution function, or DF, for this model exists, of the form
(3)
with E binding energy, Jϕ angular momentum, and q a real number that links the radial velocity dispersion σ to the circular velocity v0 via the relation |$\sigma ^2 = v_0^2/(1+q)$| (Toomre 1977). We adopt the numerical values q = 11.4 and r0 = 20, similarly to Sellwood (2012).
In order to obtain a stellar disc with a finite total mass and with a non-diverging central density, we multiply this DF with two cut-out functions, Hinner(Jϕ) and Houter(Jϕ), of the form
(4)
(5)
Here, we choose n = 4, m = 5, router = 11.5, and we adopt units such that v0 = G = rinner = 1. As in Sellwood (2012), we use a Plummer softening length ε = 1/8 in equation (A5) for the inter-particle interaction potential.
The DF is further multiplied with an active fraction ξ = 1/2, such that the DF describing the dynamics of the stellar component of the half-mass Mestel disc has the form
(6)
The rest of the matter, making up the deficit between the stellar density generated by the DF (6) and the Mestel disc density (1), is supposed to reside in a rigid, unresponsive halo and bulge. As reported by Toomre (1981), Evans & Read (1998), Sellwood (2012), and Fouvry et al. (2015), this stellar disc is linearly stable, a fact that we confirmed with our own stability analysis.

2.2 The grooved half-mass Mestel disc

While this half-mass Mestel disc is linearly stable, this does not mean that it cannot evolve secularly through finite-N effects. In order to maximally elucidate the underlying dynamics, Sellwood (2012) and Fouvry et al. (2015) consider only gravitational forces up to the m = 2 harmonic in their studies of the evolution of this disc model. Therefore, we will likewise only consider m = 2 spiral patterns in the remainder.

As shown by Sellwood (2012) and Fouvry et al. (2015), swing-amplified finite-N stellar density fluctuations in this Mestel disc transport angular momentum away from their inner Lindblad resonances (ILRs) and scatter stars away from (near-)circular orbits around certain angular momenta. If a spiral wave is time-independent in a corotating coordinate frame then each star evolves along a track of constant Jacobi integral, EJ, given by
(7)
with Ωp the wave’s pattern speed. The interaction between the stars and transient swing-amplified patterns induces stellar migration through action space, changing the stars’ binding energy and angular momentum but preserving their Jacobi integral. At a pattern’s ILR, this pushes stars away from circular orbits towards orbits with lower angular momentum Jϕ and higher values of the radial action Jr. Hence, the DF develops what we here refer to as grooves (loci around constant EJ-values where stars have diffused away from high-Jϕ orbits, leaving behind a depleted DF) and ridges (loci around constant EJ-values where stars have diffused towards high-Jr orbits, causing an enhanced DF).

In Fig. 1, we show the DF of the half-mass Mestel disc, grooved on the long term by the non-linear evolution of swing-amplified particle-shot noise, as obtained in simulation 50M of Sellwood (2012)1 [see also fig. 7 in Sellwood (2012) and fig. 4 in Fouvry et al. (2015)]. Three grooves can be identified where circular orbits are depopulated around angular momentum values Jϕ ≈ 1.2, 1.55, and 1.90. In order to rapidly convert between binding energy E, angular momentum Jϕ, and the radial action Jr, we used the formulae from Williams, Evans & Bowden (2014).

The phase-space distribution of the particles in simulation 50M of Sellwood (2012) at time 1400, presented as a function of angular momentum, denoted by Jϕ, and the radial action, denoted by Jr. Circular orbits are characterized by Jr = 0. The contours of constant stellar phase-space density are evenly spaced between 2.5 and 95 per cent of the maximum value of the DF (estimated as the mean of the 10 highest valued pixels). Three grooves can be identified where circular orbits are depopulated around the angular momentum values Jϕ ≈ 1.2, 1.55, and 1.90.
Figure 1.

The phase-space distribution of the particles in simulation 50M of Sellwood (2012) at time 1400, presented as a function of angular momentum, denoted by Jϕ, and the radial action, denoted by Jr. Circular orbits are characterized by Jr = 0. The contours of constant stellar phase-space density are evenly spaced between 2.5 and 95 per cent of the maximum value of the DF (estimated as the mean of the 10 highest valued pixels). Three grooves can be identified where circular orbits are depopulated around the angular momentum values Jϕ ≈ 1.2, 1.55, and 1.90.

Concurrent with the appearance of these grooves in action space, the stellar disc develops a set of vigorously growing spiral patterns whose amplitudes exponentiate with time (see fig. 2 in Sellwood 2012). The most vigorous among these modes has a rotation frequency ω = mΩp ≈ 0.55 (see fig. 4 in Sellwood 2012). It cannot be an eigenmode of the ungrooved half-mass Mestel disc because it is known to be linearly stable. Since the mode already appears in the linear regime, it is unlikely to be caused by non-linear mode coupling (Masset & Tagger 1997). Randomizing the azimuthal particle positions does not prevent the pattern’s growth, quite on the contrary, so any non-axisymmetric features that grew during the disc’s first phase of evolution cannot have caused this growing mode (see fig. 5 in Sellwood 2012). The fact that it grows exponentially with time suggests that it is a true eigenmode particular to the grooved half-mass Mestel disc.

The DF of the fiducial grooved half-mass Mestel disc. The contours are the same as in Fig. 1. The dotted line traces a constant Jacobi integral along the main groove.
Figure 2.

The DF of the fiducial grooved half-mass Mestel disc. The contours are the same as in Fig. 1. The dotted line traces a constant Jacobi integral along the main groove.

Using linear stability analysis, we now show that the vigorously exponentiating patterns observed by Sellwood (2012) are indeed eigenmodes of the grooved half-mass Mestel disc. For this, we use pystab, a python/C++ computer code to analyse the stability of a responsive, self-gravitating razor-thin stellar disc embedded in the gravitational field of a rigid axisymmetric or spherically symmetric central bulge and dark-matter halo. The details of the mathematical formalism behind this code and of its implementation can be found in Vauterin & Dejonghe (1996), Dury et al. (2008), and De Rijcke & Voulis (2016), so we will not repeat these here in detail. We provide a brief overview of the formalism and of the employed values of some numerical parameters of the code in Appendix  A.

3 MODE-ANALYSIS OF THE GROOVED HALF-MASS MESTEL DISC

3.1 The fiducial grooved half-mass Mestel disc

We first mimicked the grooves visible in the DF of simulation 50M (see Fig. 1) of Sellwood (2012) by multiplying the DF F0(E, Jϕ) with ‘groove functions’ of the form
(8)
with x varying along a line of constant Jacobi integral, normalized such that x = 0 at the circular orbit and x = 1 at the radial orbit. Each groove function differs from 1 only inside a narrow region of width 2w centred on a given Jacobi integral. If fgroove(x) > 1, the phase-space density is increased; if fgroove(x) < 1, the phase-space density is decreased.

Based on Fig. 1, we identify three angular momenta where circular orbits are significantly depopulated: at Jϕ ≈ 1.2 (first groove), 1.55 (second groove), and 1.9 (third groove), corresponding to the ILR of waves with pattern speeds Ωp = 0.244, 0.189, and 0.154, respectively. We will henceforth refer to the half-mass Mestel disc with these three grooves as the fiducial grooved half-mass Mestel disc.

At the first groove, at Jϕ ≈ 1.2, the value of the DF decreases by 57 per cent while, at the second and third grooves, at Jϕ ≈ 1.55 and 1.9, the value of the DF decreases by 5 per cent. The first groove leads to a marked population increase of the lower angular momentum/higher radial action orbits, which is mimicked here by setting α = 11 in equation (8). The second and, especially, the third grooves show much less pronounced DF increases towards low angular momentum orbits, so we choose α = 6 for these. For each groove, the parameters w, a, b, c, and d are chosen such that the grooves visible in Fig. 1 are adequately reproduced, cf. Figs 2 and 3, and Table 1. We took care to conserve the total mass of the stellar disc to within a few tenths of a per cent.

Top panel: a cut through the DF of the fiducial grooved half-mass Mestel disc along the Jr = 0 line, where the circular orbits reside, as a function of angular momentum Jϕ (full line). For a comparison, the dotted line traces the ungrooved DF. The grey histograms show the distribution of the particles in simulation 50Mc of Sellwood (2012) over circular orbits (based on simulation data kindly provided to us by Prof. J. Sellwood). Bottom panel: a cut through the DF of the fiducial grooved half-mass Mestel disc along the main ridge, plotted as a function of the radial action Jr (full line). The dotted line traces the ungrooved DF along the main ridge. In particular, we note the depletion of stars on quasi-circular orbits, and the addition of stars on more radial orbits.
Figure 3.

Top panel: a cut through the DF of the fiducial grooved half-mass Mestel disc along the Jr = 0 line, where the circular orbits reside, as a function of angular momentum Jϕ (full line). For a comparison, the dotted line traces the ungrooved DF. The grey histograms show the distribution of the particles in simulation 50Mc of Sellwood (2012) over circular orbits (based on simulation data kindly provided to us by Prof. J. Sellwood). Bottom panel: a cut through the DF of the fiducial grooved half-mass Mestel disc along the main ridge, plotted as a function of the radial action Jr (full line). The dotted line traces the ungrooved DF along the main ridge. In particular, we note the depletion of stars on quasi-circular orbits, and the addition of stars on more radial orbits.

Table 1.

The parameter values employed in the ‘groove function’ fgroove that defines the profile of each groove (cf. equation 8).

αwabcd
First groove110.290.474−0.044−5.219−44.796
Second groove60.100.5620.388−3.371−15.619
Third groove60.120.0810.869−0.486−3.242
αwabcd
First groove110.290.474−0.044−5.219−44.796
Second groove60.100.5620.388−3.371−15.619
Third groove60.120.0810.869−0.486−3.242
Table 1.

The parameter values employed in the ‘groove function’ fgroove that defines the profile of each groove (cf. equation 8).

αwabcd
First groove110.290.474−0.044−5.219−44.796
Second groove60.100.5620.388−3.371−15.619
Third groove60.120.0810.869−0.486−3.242
αwabcd
First groove110.290.474−0.044−5.219−44.796
Second groove60.100.5620.388−3.371−15.619
Third groove60.120.0810.869−0.486−3.242

3.2 Linear stability analysis of the fiducial grooved Mestel disc

As explained in Appendix  A, each eigenmode corresponds to a complex frequency ω for which the matrix |${\mathcal {C}}(\omega)$|⁠, defined in equation (A7), has a unity eigenvalue. In Fig. 4, we plot the quantity min |λ(ω) − 1| in the complex frequency plane of the fiducial grooved half-mass Mestel disc. This quantity is zero at the loci of the eigenmodes.

The quantity min |λ(ω) − 1| in the complex frequency plane of the fiducial grooved half-mass Mestel disc. An eigenmode lives wherever this quantity is zero (dark regions).
Figure 4.

The quantity min |λ(ω) − 1| in the complex frequency plane of the fiducial grooved half-mass Mestel disc. An eigenmode lives wherever this quantity is zero (dark regions).

The complex frequencies ω of the dominant eigenmodes retrieved with pystab for the fiducial grooved half-mass Mestel disc are listed in Table 2, along with the radii where the main resonances (ILR, CR, OLR) occur, with CR the corotation resonance radius and OLR the outer Lindblad resonance radius. Resonances that (approximately) overlap with the position of a groove are printed in boldface. For the Mestel disc, the resonance radii rres, ℓ of an m-armed mode with pattern speed mΩp = ωreal, with ωreal the real part of the mode frequency, are given by
(9)
Here, ℓ = 0 gives the CR while ℓ = ±1 gives the Lindblad resonances.
Table 2.

The eigenmodes of the half-mass Mestel disc with three grooves at Jϕ = 1.2, 1.55, and 1.9. The complex mode frequency is denoted by ω. The radii of the main resonances (ILR, CR, OLR) are indicated and resonances that (approximately) overlap with the position of a groove are printed in boldface. The physical nature of each mode (global or groove mode) is given in the last column.

ωILRCROLRType
0.597 + 0.013i0.983.355.72global mode
1.662 + 0.004i0.351.202.01groove mode
0.465 + 0.001i1.264.307.34global mode
1.101 + 0.000i0.531.823.10groove mode
ωILRCROLRType
0.597 + 0.013i0.983.355.72global mode
1.662 + 0.004i0.351.202.01groove mode
0.465 + 0.001i1.264.307.34global mode
1.101 + 0.000i0.531.823.10groove mode
Table 2.

The eigenmodes of the half-mass Mestel disc with three grooves at Jϕ = 1.2, 1.55, and 1.9. The complex mode frequency is denoted by ω. The radii of the main resonances (ILR, CR, OLR) are indicated and resonances that (approximately) overlap with the position of a groove are printed in boldface. The physical nature of each mode (global or groove mode) is given in the last column.

ωILRCROLRType
0.597 + 0.013i0.983.355.72global mode
1.662 + 0.004i0.351.202.01groove mode
0.465 + 0.001i1.264.307.34global mode
1.101 + 0.000i0.531.823.10groove mode
ωILRCROLRType
0.597 + 0.013i0.983.355.72global mode
1.662 + 0.004i0.351.202.01groove mode
0.465 + 0.001i1.264.307.34global mode
1.101 + 0.000i0.531.823.10groove mode

3.3 The ω = 0.597 + 0.013i mode: a cavity mode

The most rapidly growing global spiral mode has a pattern speed mΩp = 0.597, which places its ILR at a radius of 0.98, very near the inner edge of the first groove, while its OLR occurs at a radius of 5.72. The other global spiral mode, with a pattern speed mΩp = 0.465, has a negligible growth rate (this is actually the fastest growing member of a cluster of eigenmodes with pattern speeds between 0.45 and 0.5, cf. Fig. 4). The surface density of the dominant global spiral mode is plotted in Fig. 5 (top panel) where we compare it with the dominant spiral pattern found in the N-body simulations reported in Sellwood (2012) (bottom panel). The frequency of the emerging pattern is estimated2 at ω ≈ 0.577 + 0.018i in simulation 50M and at ω ≈ 0.599 + 0.016i in simulation 50Mc. We find ω = 0.597 + 0.013i for the dominant eigenmode of the fiducial grooved Mestel disc. This good agreement between the pattern speeds leads us to confirm that we are, in fact, seeing the same mode in the linear stability analysis and in the N-body simulations.

Top panel: the positive part of the surface density perturbation of the m = 2 global mode with frequency 0.597 + 0.013i of the fiducial grooved half-mass Mestel disc. The contour values increase linearly between 1 and 99 per cent of the maximum surface density value. The full line indicates the position of the first, main groove and dashed lines mark the mode’s ILR, CR, and OLR. Bottom panel: the surface density perturbation of this mode taken from N-body simulation 50Mc by Sellwood (2012). Dotted lines mark the Lindblad resonances while the full line marks the CR radius.
Figure 5.

Top panel: the positive part of the surface density perturbation of the m = 2 global mode with frequency 0.597 + 0.013i of the fiducial grooved half-mass Mestel disc. The contour values increase linearly between 1 and 99 per cent of the maximum surface density value. The full line indicates the position of the first, main groove and dashed lines mark the mode’s ILR, CR, and OLR. Bottom panel: the surface density perturbation of this mode taken from N-body simulation 50Mc by Sellwood (2012). Dotted lines mark the Lindblad resonances while the full line marks the CR radius.

3.3.1 Mode generating mechanism: linear stability analysis

As noted by Sellwood & Carlberg (2014), a phase-space groove induces a spatially localized change of the velocity with which spiral waves can be radially propagated through the disc. A groove locally changes the ‘impedance’ of the stellar disc. Using the suggestive analogy of waves being partially reflected at the junction between two strings with different mass densities (and hence different wave velocities), these authors argue that stellar density waves impacting radially inward on a phase-space groove will, likewise, be partially reflected outward again. This creates the possibility of setting up a resonant cavity within which inwards moving waves are (partially) reflected back out again by the groove and outwards moving waves are reflected back inwards at corotation, thus creating a standing wave pattern.

In the WASER feedback cycle (Toomre 1969; Mark 1976), an ingoing short-wavelength trailing wave is reflected from the unresponsive, dynamically hot disc interior [if the Toomre Q parameter (Toomre 1964) rises sufficiently rapidly towards small radii] into an outgoing long-wavelength trailing wave. At CR, this wave overreflects into an ingoing and an outgoing short-wavelength trailing wave. Hence, the WASER only involves trailing waves. Moreover, it leads to rather modest wave amplifications and then only for low Q values. Much more impressive amplifications can be obtained by the swing amplifier feedback cycle (Toomre 1977, 1981; Bertin et al. 1989). In this theory, an ingoing trailing wave is reflected from the galaxy centre as an outgoing leading wave. Around CR, this wave is overreflected into an outgoing and an ingoing trailing wave at which point amplification factors of 1–2 orders of magnitude can be achieved. Significant amplification occurs for discs with Q ≲ 2. Interference between the leading and trailing waves within the resonant cavity is expected to produce a growing eigenmode with a ‘lumpy’ density distribution (Athanassoula 1984; Binney & Tremaine 2008).

If the inner Q-profile of the half-mass Mestel disc were steep enough to reflect ingoing waves outwards again, it would already have done so in the ungrooved model and that model is linearly stable. Therefore, it is unlikely that the WASER feedback cycle is responsible for the growth of the global modes in the grooved half-mass Mestel disc. On the other hand, the Toomre Q parameter is smaller than 2.0 in a large part of disc, between radii in the range ≈1.4–9.0, and hovers around a value of 1.5 between radii in the range ≈2.0–6.0. From the famous ‘dust to ashes’ figure of Toomre (1981), it is clear that swing amplification in the Q = 1.5, X = kcritr/m = 2 half-mass Mestel disc without inner or outer cut-outs can boost wave amplitudes by a factor of ∼30 (which significantly hastens the disc’s secular dynamics). Unfortunately, those amplified waves are subsequently absorbed at the ILR. This is, however, not an issue for the grooved half-mass Mestel disc. In the case of waves with a pattern speed mΩp = 0.597, corresponding to the dominant global mode of the fiducial grooved half-mass Mestel disc, the ILR (at a radius of 0.98) is shielded by the reflective first groove (around a radius of 1.2). Moreover, this dominant global mode clearly has a ‘lumpy’ density distribution (cf. Fig. 5), indicative of interference between leading and trailing waves.

3.3.2 Mode generating mechanism: WKB arguments

If an eigenmode is indeed a standing wave built from wave packets travelling inside a resonance cavity between radii rinf and rsup then it must obey the quantum condition
(10)
where k is the radial wavenumber of the waves. The mode’s growth rate ωim = ℑ{ω} then follows from the relation
(11)
with vg(r) the radial group velocity of the travelling wave patterns (Mark 1977; Bertin 2014). Hence, the narrower the cavity and the larger the group velocity, the stronger the growth. These integrals are to be performed over a full circle between their lower and upper bounds.
In the WKBJ limit, the dispersion relation for travelling waves is given by
(12)
Here,
(13)
with Ω and κ the angular and epicycle frequency of near-circular orbits, respectively (Kalnajs 1965; Lin & Shu 1966). The reduction factor on the right-hand side of this equation is the product of two factors: an exponential factor caused by the Plummer softening with softening length ε (Romeo 1994), for which we adopt the value ε = 1/8, and the well-known reduction factor
(14)
with In a modified Bessel function of the first kind, due to the stellar velocity dispersion. This latter expression is only correct for a Schwarzschild DF so the following results must be regarded as approximate. The group velocity vg of a wave packet centred on wavenumber k is given by
(15)
From the dispersion relation stated above, we can derive the relation
(16)
for the group velocity. The partial derivatives of the reduction factor |$\mathcal {F}$| that figure in this equation can be evaluated analytically as
(17)
where all Bessel functions have χ as argument. We terminate the infinite sums in the evaluation of this reduction factor and its derivatives when a relative accuracy of 10−6 has been obtained.

We now wish to compute the frequency of a growing standing wave in the half-mass Mestel disc with a groove at a radius rgroove, which is supposed to reflect incoming wave packets back out again, using this analytical theory. For a given real frequency ωreal, we first compute the inner radius of the forbidden region around CR, where the dispersion relation (12) admits no real solutions for k. This we equate to the outer radius of the resonance cavity, rsup. For the inner radius, rinf, we use the maximum of the groove radius and the wave’s ILR radius. In practice, this always turned out to be the groove radius.

For any radius r inside the resonance cavity, we then find the two positive solutions kshortreal, r) and klongreal, r) of the dispersion relation (12). The quantum condition (10) can be rewritten as
(18)
Numerical root finding finally allows us to solve this equation for the real frequency ωreal in case the disc supports a cavity mode. Equation (16) yields the group velocities, vg, short(r) and vg, long(r), of the short and long branch waves, respectively. These are then inserted into equation (11), which can be rewritten as
(19)
This, finally, provides us with an estimate for the growth rate ωim of the growing mode.

For the fiducial value of the groove radius, rgroove = 1.2, this yields the analytical frequency estimate ω = 0.510 + 0.017i. This is tantalizingly close to the value ω = 0.597 + 0.013i that we derived from the full linear theory. Together with all the arguments given above, this can be taken as evidence for the idea that the dominant global mode of the grooved Mestel disc is caused by swing-amplified travelling wave packets inside a resonance cavity between the groove – rendered reflective by the depopulation of near-circular orbits – and the mode’s CR.

3.4 The ω = 1.662 + 0.004i mode: a groove mode

Glancing back at Table 2, the second most unstable mode of the grooved Mestel disc, depicted in Fig. 6, has properties that suggest that it is a groove mode caused by the first groove,3 depopulating circular orbits around angular momentum Jϕ = 1.2. There actually appears to exist a cluster of eigenmodes with pattern speeds very close to this value, cf. Fig. 4, of which this is the most rapidly growing member.

The positive part of the surface density perturbation of the m = 2 mode of the fiducial grooved half-mass Mestel disc, with frequency ω = 1.662 + 0.004i. The contour values increase linearly between 1 and 99 per cent of the maximum surface density value. The dashed lines mark the modes’ ILR, CR, and OLR and the full line indicates the position of the first, main groove.
Figure 6.

The positive part of the surface density perturbation of the m = 2 mode of the fiducial grooved half-mass Mestel disc, with frequency ω = 1.662 + 0.004i. The contour values increase linearly between 1 and 99 per cent of the maximum surface density value. The dashed lines mark the modes’ ILR, CR, and OLR and the full line indicates the position of the first, main groove.

Sellwood & Kahn (1991) describe growing instabilities with their CR at a sharp groove, i.e. a local depression in the stellar density distribution, that are an instance of the negative-mass instability (Lovelace & Hohlfeld 1978). These authors show that such groove modes originate from the coupling of two waves, one on each side of the groove, with their corotation resonance inside the groove. The wavy displacements of disc material at the groove’s edges couple across the groove through gravity. Thus, local overdensities can exchange angular momentum in such a way that material on the groove’s inner edge gains angular momentum and is pulled outwards while material on the groove’s outer edge loses angular momentum and is pulled inwards. This counteracts the disc’s shear and facilitates mass clumping. Thus, the amplitude of these mass displacements grows exponentially with time. The rest of the pattern can then be explained as the disc’s response to these growing overdensities via the mechanism described by Julian & Toomre (1966).

The CR of the ω = 1.662 + 0.004i mode indeed sits squarely inside the first groove, cf. Fig. 7. It is also clear from this figure that the mode, as expected for a groove mode, consists of two density peaks sitting on either side of the groove with a gravitational potential peak between them. This suggests this is a groove mode caused by the first groove. There also appears to exist a neutral eigenmode with a pattern speed mΩp = 1.101 that can, likewise, be interpreted as a groove mode caused by the third groove. Its vanishing growth rate explains why it did not appear in the simulations reported by Sellwood (2012).

Close-up view of the density peak of the ω = 1.662 + 0.004i groove mode from Fig. 6. Full-line contours trace the surface density perturbation (in 10 equidistant steps between 1 and 99 per cent of its maximum value) while the dotted contours trace the mode’s gravitational binding potential (in 10 equidistant steps between 1 and 99 per cent of its maximum value). Dashed circles mark the loci of the mode’s ILR, CR, and OLR. The two full-line circles demarcate the region where the groove reaches 50 per cent of its maximum depth, which is a measure for its width.
Figure 7.

Close-up view of the density peak of the ω = 1.662 + 0.004i groove mode from Fig. 6. Full-line contours trace the surface density perturbation (in 10 equidistant steps between 1 and 99 per cent of its maximum value) while the dotted contours trace the mode’s gravitational binding potential (in 10 equidistant steps between 1 and 99 per cent of its maximum value). Dashed circles mark the loci of the mode’s ILR, CR, and OLR. The two full-line circles demarcate the region where the groove reaches 50 per cent of its maximum depth, which is a measure for its width.

4 DISCUSSION

4.1 The number of grooves

The presence of the minor second and third grooves has little bearing on the existence of the modes although it does slightly impact their frequencies. For instance, if we remove these two minor grooves, keeping only the deep and wide first groove, the main global mode is retrieved with a frequency ω = 0.603 + 0.016i, with its ILR sitting at a radius of 0.97, still just inside the inner edge of the groove. The groove mode now has a frequency ω = 1.695 + 0.003i. Its CR, now at a radius 1.18, again falls inside the groove. The remaining groove is obviously instrumental to the existence of these modes. Without it, the grooveless half-mass Mestel disc is, as already mentioned, linearly stable.

4.2 The groove profile

A groove consists both of a depletion of circular orbits and an overpopulation of more radial orbits. In order to test which of these two features is actually causing the modes in the grooved Mestel disc, we first computed the eigenmode spectrum of a model with the fiducial groove function modified according to the prescription
(20)
so that no stars are ever removed but only added. This removes the grooves at the loci of near-circular orbits but leaves the ridges at more radial orbits in place. This particular model turned out to have a single (barely) growing global eigenmode at ω = 0.611 + 0.001i. The high-frequency mode of the fiducial model is completely absent.
Next, we computed the eigenmode spectrum of a model with the fiducial groove function modified according to the prescription
(21)
so that no stars are ever added but only removed. This removes the ridges at more radial orbits but leaves the grooves at the loci of near-circular orbits in place. This model has two growing eigenmodes, at ω = 0.593 + 0.013i and ω = 1.663 + 0.003i. These are virtually the same as the eigenmodes of the fiducial grooved Mestel disc. Likewise, we analysed a model in which the depth of the three grooves of the fiducial model is kept constant along lines of constant Jacobi integral. In other words, the first groove depresses the DF by 57 per cent and the second and third grooves by 5 per cent everywhere along the groove in question. The stellar mass of this model is a bit smaller than in the fiducial model but this has no apparent consequences. Indeed, this model’s mode spectrum is almost identical to that of the fiducial model.

Putting these results together, we are led to the conclusion that it is the depopulation of the near-circular orbits and not the overpopulation of the radial orbits that is crucial to the existence of the growing eigenmodes of the grooved Mestel disc. This corroborates the results from Sellwood & Kahn (1991), who use an approximate analytical computation to demonstrate that grooves are destabilizing while ridges are not.

4.3 The groove depth

We varied the depths of the grooves, keeping their respective depths in proportion. This leads to the eigenmode spectra listed in Table 3. Clearly, the deeper the grooves, the more rapidly the mode with a pattern speed mΩp ≈ 0.6 grows. This confirms the finding by Sellwood (2012) and Sellwood & Carlberg (2014) that new instabilities only start to grow when the phase-space grooves are sufficiently deep.

Table 3.

Reprise of Table 2, with the eigenmodes of the half-mass Mestel disc with three grooves at Jϕ = 1.2, 1.55, and 1.9 with various depths. The depth of the first groove, expressed relative to its fiducial depth, is indicated in the first column. The depths of the other two grooves are varied proportionally. The complex mode frequency is denoted by ω. The radii of the main resonances (ILR, CR, OLR) are indicated and resonances that (approximately) overlap with the position of a groove are printed in boldface. The physical nature of each mode (global or groove mode) is given in the last column.

DepthωILRCROLR Type
155 per cent0.582 + 0.019i1.003.445.86global mode
1.709 + 0.016i0.341.172.00groove mode
1.078 + 0.010i0.541.863.17groove mode
0.471 + 0.001i1.244.257.25global mode
100 per cent0.597 + 0.013i0.983.355.72global mode
1.662 + 0.004i0.351.182.01groove mode
0.465 + 0.001i1.264.317.35global mode
80 per cent0.604 + 0.008i0.973.315.65global mode
60 per cent0.611 + 0.002i0.963.275.59global mode
DepthωILRCROLR Type
155 per cent0.582 + 0.019i1.003.445.86global mode
1.709 + 0.016i0.341.172.00groove mode
1.078 + 0.010i0.541.863.17groove mode
0.471 + 0.001i1.244.257.25global mode
100 per cent0.597 + 0.013i0.983.355.72global mode
1.662 + 0.004i0.351.182.01groove mode
0.465 + 0.001i1.264.317.35global mode
80 per cent0.604 + 0.008i0.973.315.65global mode
60 per cent0.611 + 0.002i0.963.275.59global mode
Table 3.

Reprise of Table 2, with the eigenmodes of the half-mass Mestel disc with three grooves at Jϕ = 1.2, 1.55, and 1.9 with various depths. The depth of the first groove, expressed relative to its fiducial depth, is indicated in the first column. The depths of the other two grooves are varied proportionally. The complex mode frequency is denoted by ω. The radii of the main resonances (ILR, CR, OLR) are indicated and resonances that (approximately) overlap with the position of a groove are printed in boldface. The physical nature of each mode (global or groove mode) is given in the last column.

DepthωILRCROLR Type
155 per cent0.582 + 0.019i1.003.445.86global mode
1.709 + 0.016i0.341.172.00groove mode
1.078 + 0.010i0.541.863.17groove mode
0.471 + 0.001i1.244.257.25global mode
100 per cent0.597 + 0.013i0.983.355.72global mode
1.662 + 0.004i0.351.182.01groove mode
0.465 + 0.001i1.264.317.35global mode
80 per cent0.604 + 0.008i0.973.315.65global mode
60 per cent0.611 + 0.002i0.963.275.59global mode
DepthωILRCROLR Type
155 per cent0.582 + 0.019i1.003.445.86global mode
1.709 + 0.016i0.341.172.00groove mode
1.078 + 0.010i0.541.863.17groove mode
0.471 + 0.001i1.244.257.25global mode
100 per cent0.597 + 0.013i0.983.355.72global mode
1.662 + 0.004i0.351.182.01groove mode
0.465 + 0.001i1.264.317.35global mode
80 per cent0.604 + 0.008i0.973.315.65global mode
60 per cent0.611 + 0.002i0.963.275.59global mode

Moreover, new eigenmodes appear as the grooves are deepened. For instance, the mode with a frequency ω = 1.078 + 0.010i, which appears in the spectrum of the model with a first groove that is 1.55 times as deep as in the fiducial model, is a groove mode with its CR at the position of the third groove, at angular momentum Jϕ = 1.9. This corroborates the results from Sellwood & Kahn (1991), who find that a groove of a given width needs to exceed a critical depth before it can create a groove mode. Using approximate analytical calculations, these authors predict the critical depth to be proportional to the groove’s width, which is in rough agreement with the required groove depths we observe here.

4.4 The position of the grooves

An investigation of the role of the positions of the grooves is not possible with the techniques of Sellwood (2012) and Fouvry et al. (2015), where the grooves grow due to stellar dynamics and their position is determined by the physics of the problem. Here, on the contrary, we have full control over the positions of the grooves. We simply change the angular momentum of the circular orbits depleted by the first groove between 0.8 and 3.5 while maintaining a constant pattern speed ratio between the three grooves and keeping all other groove properties constant. This causes the grooves to jointly sweep across the phase-space region with the highest stellar density.

The result is plotted in Fig. 8. Here, the black bullets indicate the frequency of the dominant global eigenmode in the complex frequency plane for different positions of the grooves. Each data point is labelled with the angular momentum Jϕ of the circular orbits depleted by the first groove as it sweeps through phase space. Growing eigenmodes are triggered only if the grooves fall in a ‘responsive’ region of phase space (see also De Rijcke & Voulis (2016)). In this case, this means a groove carving into the most densely populated region of phase space: that of the circular orbits with radii between 1.5 and 2.5. A groove falling outside this region has a much less marked effect, as can be expected from removing stars from sparsely populated phase-space surroundings.

Position of the dominant m = 2 eigenmode in the complex frequency plane of the half-mass Mestel disc as a function of groove position, computed using the full linear theory (black bullets) and using the approximate calculation from paragraph 3.3.2 (grey stars). The angular momentum of the first groove takes on values between Jϕ = 0.8 and 3.5, as indicated next to each data point. The ratio of the angular momentum of the second and third grooves to that of the first groove is kept constant, causing the grooves to jointly sweep through phase space. The fiducial grooved Mestel disc corresponds to Jϕ = 1.2.
Figure 8.

Position of the dominant m = 2 eigenmode in the complex frequency plane of the half-mass Mestel disc as a function of groove position, computed using the full linear theory (black bullets) and using the approximate calculation from paragraph 3.3.2 (grey stars). The angular momentum of the first groove takes on values between Jϕ = 0.8 and 3.5, as indicated next to each data point. The ratio of the angular momentum of the second and third grooves to that of the first groove is kept constant, causing the grooves to jointly sweep through phase space. The fiducial grooved Mestel disc corresponds to Jϕ = 1.2.

The grey stars in Fig. 8 show the analytical estimate for the frequency of the dominant global eigenmode as a function of groove position, using the analytical formalism detailed in paragraph 3.3.2. Each star is labelled with the angular momentum Jϕ of the circular orbits depleted by the first groove. This estimate, derived from an analytical calculation based on the assumption that we are dealing with a cavity mode between the groove radius and the forbidden region around corotation, is clearly in rough agreement with the full linear theory result. Given the approximations underlying this calculation, this level of agreement is actually surprisingly good and can be taken as strong evidence for our interpretation of the dominant global eigenmode as a cavitymode.

4.5 Gravitational softening

Gravitational softening with a Plummer kernel with a softening length ε = 1/8, as employed by Sellwood (2012), has a very strong impact on the eigenmode spectrum of the grooved Mestel disc. In Table 4, we list the dominant eigenmodes of the unsoftened grooved Mestel disc model. Without gravitational softening, the groove modes grow much faster than in the softened model. For instance, the groove mode caused by the first groove is now by far the fastest growing mode. The global modes around mΩp ≈ 0.45 and mΩp ≈ 0.6 appear to be relatively unaffected by softening but they are now surpassed in growth rate by the groove modes. Moreover, new eigenmodes appear in the spectrum of the unsoftened model: the very rapidly growing mode at ω = 0.6740 + 0.040i and the slowly growing mode at ω = 0.5420 + 0.003i.

Table 4.

The eigenmodes of the unsoftened half-mass Mestel disc with three grooves at Jϕ = 1.2, 1.55, and 1.9. The complex mode frequency is denoted by ω. The radii of the main resonances (ILR, CR, OLR) are indicated and resonances that (approximately) overlap with the position of a groove are printed in boldface. The physical nature of each mode (global or groove mode) is given in the last column.

ωILRCROLRType
1.7260 + 0.068i0.341.161.98groove mode
0.6740 + 0.040i0.872.975.06global mode
0.5470 + 0.011i1.073.656.24global mode
1.0650 + 0.004i0.551.883.21groove mode
0.5420 + 0.003i1.083.696.29global mode
0.4200 + 0.001i1.394.768.13global mode
ωILRCROLRType
1.7260 + 0.068i0.341.161.98groove mode
0.6740 + 0.040i0.872.975.06global mode
0.5470 + 0.011i1.073.656.24global mode
1.0650 + 0.004i0.551.883.21groove mode
0.5420 + 0.003i1.083.696.29global mode
0.4200 + 0.001i1.394.768.13global mode
Table 4.

The eigenmodes of the unsoftened half-mass Mestel disc with three grooves at Jϕ = 1.2, 1.55, and 1.9. The complex mode frequency is denoted by ω. The radii of the main resonances (ILR, CR, OLR) are indicated and resonances that (approximately) overlap with the position of a groove are printed in boldface. The physical nature of each mode (global or groove mode) is given in the last column.

ωILRCROLRType
1.7260 + 0.068i0.341.161.98groove mode
0.6740 + 0.040i0.872.975.06global mode
0.5470 + 0.011i1.073.656.24global mode
1.0650 + 0.004i0.551.883.21groove mode
0.5420 + 0.003i1.083.696.29global mode
0.4200 + 0.001i1.394.768.13global mode
ωILRCROLRType
1.7260 + 0.068i0.341.161.98groove mode
0.6740 + 0.040i0.872.975.06global mode
0.5470 + 0.011i1.073.656.24global mode
1.0650 + 0.004i0.551.883.21groove mode
0.5420 + 0.003i1.083.696.29global mode
0.4200 + 0.001i1.394.768.13global mode

5 CONCLUSIONS

We have studied the linear stability of a half-mass Mestel disc, with a distribution function modified with phase-space grooves like the ones found in the works of Sellwood (2012), Fouvry et al. (2015), and Fouvry (2017). The ungrooved model is linearly stable so that any modes found in the grooved models must somehow be related to the presence of the grooves. Since we compare our linear stability computations with N-body simulations, we take into account the effects of gravitational softening (De Rijcke, Fouvry & Dehnen 2018).

In the fiducial grooved model, we mimic the grooves visible in the DF of the N-body simulation 50M of Sellwood (2012) and show that this model has m = 2 eigenmodes with properties very close to the two-armed spiral patterns visible in its simulated analogue. We argue that, at least in this case, the dominant new eigenmode is a cavity mode, living inside a resonance cavity between the groove – rendered reflective by the removal of stars from near-circular orbits – and the forbidden region around corotation (Mark 1977; Bertin 2014). Other, more slowly growing new eigenmodes can be interpreted as groove modes (Sellwood & Kahn 1991).

This confirms the suggestion put forward in Sellwood (2012) and Sellwood & Carlberg (2014) that grooves can be a potent source of new, rapidly growing eigenmodes. These grooves are produced by a series of uncorrelated swing-amplified waves, sourced by the finite-N nature of any self-gravitating stellar system. In other words, it is possible for an isolated, linearly stable stellar disc, initially without recourse to any growing eigenmodes that can transport its angular momentum outwards as desired by the second law of thermodynamics (Lynden-Bell & Kalnajs 1972), to spontaneously become linearly unstable via the formation of grooves in phase space through finite-N dynamics. The simple model of the half-mass Mestel disc with inner and outer cut-outs serves as a template for more realistic galaxy models, equipped with more plausible bulge and halo components. The idea is that the dynamical processes leading to the vigorously growing spiral modes in this Mestel disc also explain the spirals in N-body simulations of more complex galaxy models, something which needs to be investigated further.

Ultimately, one hopes that understanding these dynamical processes in simulated galaxies furthers our understanding of the growth and maintenance of spiral structure in real galaxies. Real spiral galaxies evidently contain much more stars than the ∼106–108 particles that are routinely used in N-body simulations and their stellar bodies are, therefore, much smoother than those of their simulated counterparts. As shown by Sellwood (2012) and Fouvry et al. (2015), the time at which the disc’s dynamics, driven by finite-N shot noise, has carved sufficiently deep phase-space grooves and the new eigenmodes start to manifest themselves grows with the number N of particles, so one can wonder whether the mechanism discussed here actually applies to real galaxies. However, whereas the stellar bodies of real galaxies are perhaps much smoother than those of simulated galaxies, they possess other potent sources of gravitational ‘noise’, such as giant molecular clouds or globular clusters moving in or through the stellar disc (D’Onghia et al. 2013), orbiting satellite dwarf galaxies (Martínez-Delgado et al. 2008; Laporte et al. 2018), and dark-matter subhaloes (Dubinski et al. 2008; Chequers, Widrow & Darling 2018; Hu & Sijacki 2018), which could play the same role as the particle noise present in simulations. Moreover, stellar discs are dynamically cold systems so that through self-gravity any perturbation is dressed by a polarization cloud, which further hastens the disc’s evolution.

ACKNOWLEDGEMENTS

We dedicate this paper to the memory of Donald Lynden-Bell (1935–2018). The study of spiral patterns in disc galaxies is but one of many fields of astronomy that he enriched with inspiring and crucial contributions.

We wish to thank Prof. Jerry Sellwood for kindly providing us with additional information about the N-body simulations presented in Sellwood (2012). We also thank the referee for a fast and helpful report. SDR acknowledges financial support from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 721463 to the SUNDIAL ITN network. JBF acknowledges support from Program number HST-HF2-51374 which was provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555. This research is carried out in part within the context of Spin(e) (ANR-13-BS05-0005, http://cosmicorigin.org).

Footnotes

1

Based on simulation data kindly provided to us by Prof. J. Sellwood.

2

Private communication with J. Sellwood.

3

Private communication with J. Sellwood: power spectra of the simulations described in Sellwood (2012) show the existence of a slowly growing pattern with a pattern speed mΩp ≈ 1.65 when extended to higher pattern frequencies than was reported in fig. 4 of that work.

REFERENCES

Aramyan
L. S.
et al. ,
2016
,
MNRAS
,
459
,
3130

Athanassoula
E.
,
1984
,
Phys. Rep.
,
114
,
319

Beckman
J. E.
,
Font
J.
,
Borlaff
A.
,
García-Lorenzo
B.
,
2018
,
ApJ
,
854
,
182

Bertin
G.
,
2014
,
Dynamics of Galaxies
.
Cambridge Univ. Press
,
Cambridge

Bertin
G.
,
Lin
C. C.
,
Lowe
S. A.
,
Thurstans
R. P.
,
1989
,
ApJ
,
338
,
104

Binney
J.
,
Tremaine
S.
,
2008
,
Galactic Dynamics
, 2nd edn.
Princeton University Press
,
Princeton, New Jersey

Chavanis
P.-H.
,
2012
,
Phys. A: Stat. Mech. Appl.
,
391
,
3680

Chequers
M. H.
,
Widrow
L. M.
,
Darling
K.
,
2018
,
MNRAS
,
480
,
4244

D’Onghia
E.
,
Vogelsberger
M.
,
Hernquist
L.
,
2013
,
ApJ
,
766
,
34

Daniel
K. J.
,
Wyse
R. F. G.
,
2018
,
MNRAS
,
476
,
1561

De Rijcke
S.
,
Voulis
I.
,
2016
,
MNRAS
,
456
,
2024

De Rijcke
S.
,
Fouvry
J. B.
,
Dehnen
W.
,
2018
,
MNRAS
,
preprint (arxiv.1812.07001)

Dubinski
J.
,
Gauthier
J.-R.
,
Widrow
L.
,
Nickerson
S.
,
2008
, in
Funes
J. G.
,
Corsini
E. M.
, eds,
ASP Conf. Ser. Vol. 396, Formation and Evolution of Galaxy Disks (Spiral and Bar Instabilities Provoked by Dark Matter Satellites)
.
Astron. Soc. Pac
,
San Francisco
, p.
321

Dury
V.
,
de Rijcke
S.
,
Debattista
V. P.
,
Dejonghe
H.
,
2008
,
MNRAS
,
387
,
2

Evans
N. W.
,
Read
J. C. A.
,
1998
,
MNRAS
,
300
,
106

Fouvry
J.-B.
,
2017
,
Secular Evolution of Self-gravitating Systems Over Cosmic Age
.
Springer International Publishing AG
,
Switzerland

Fouvry
J. B.
,
Pichon
C.
,
Magorrian
J.
,
Chavanis
P. H.
,
2015
,
A&A
,
584
,
A129

Goldreich
P.
,
Lynden-Bell
D.
,
1965
,
MNRAS
,
130
,
125

Hart
R. E.
,
Bamford
S. P.
,
Casteels
K. R. V.
,
Kruk
S. J.
,
Lintott
C. J.
,
Masters
K. L.
,
2017
,
MNRAS
,
468
,
1850

Heyvaerts
J.
,
2010
,
MNRAS
,
407
,
355

Hu
S.
,
Sijacki
D.
,
2018
,
MNRAS
,
478
,
1576

Hunter
C.
,
1965
,
MNRAS
,
129
,
321

Julian
W. H.
,
Toomre
A.
,
1966
,
ApJ
,
146
,
810

Kalnajs
A. J.
,
1965
,
PhD thesis
,
Harvard Univ

Kalnajs
A. J.
,
1977
,
ApJ
,
212
,
637

Laporte
C. F. P.
,
Johnston
K. V.
,
Gómez
F. A.
,
Garavito-Camargo
N.
,
Besla
G.
,
2018
,
MNRAS
,
481
,
286

Lin
C. C.
,
Shu
F. H.
,
1964
,
ApJ
,
140
,
646

Lin
C. C.
,
Shu
F. H.
,
1966
,
Proc. Natl. Acad. Sci. USA.
,
55
,
229

Lovelace
R. V. E.
,
Hohlfeld
R. G.
,
1978
,
ApJ
,
221
,
51

Lynden-Bell
D.
,
Kalnajs
A. J.
,
1972
,
MNRAS
,
157
,
1

Mark
J. W.-K.
,
1977
,
ApJ
,
212
,
645

Mark
J. W. K.
,
1976
,
ApJ
,
205
,
363

Martínez-Delgado
D.
,
Peñarrubia
J.
,
Gabany
R. J.
,
Trujillo
I.
,
Majewski
S. R.
,
Pohlen
M.
,
2008
,
ApJ
,
689
,
184

Masset
F.
,
Tagger
M.
,
1997
,
A&A
,
322
,
442

Meidt
S. E.
,
Rand
R. J.
,
Merrifield
M. R.
,
2009
,
ApJ
,
702
,
277

Mestel
L.
,
1963
,
MNRAS
,
126
,
553

Omurkanov
T. Z.
,
Polyachenko
E. V.
,
2014
,
Astron. Lett.
,
40
,
724

Pichon
C.
,
Cannon
R. C.
,
1997
,
MNRAS
,
291
,
616

Romeo
A. B.
,
1994
,
A&A
,
286
,
799

Roškar
R.
,
Debattista
V. P.
,
Quinn
T. R.
,
Stinson
G. S.
,
Wadsley
J.
,
2008
,
ApJ
,
684
,
L79

Roškar
R.
,
Debattista
V. P.
,
Quinn
T. R.
,
Wadsley
J.
,
2012
,
MNRAS
,
426
,
2089

Saha
K.
,
Elmegreen
B.
,
2016
,
ApJ
,
826
,
L21

Salo
H.
,
Laurikainen
E.
,
Buta
R.
,
Knapen
J. H.
,
2010
,
ApJ
,
715
,
L56

Seigar
M. S.
,
James
P. A.
,
2002
,
MNRAS
,
337
,
1113

Sellwood
J. A.
,
2011
,
MNRAS
,
410
,
1637

Sellwood
J. A.
,
2012
,
ApJ
,
751
,
44

Sellwood
J. A.
,
2014
,
Rev. Mod. Phys.
,
86
,
1

Sellwood
J. A.
,
Binney
J. J.
,
2002
,
MNRAS
,
336
,
785

Sellwood
J. A.
,
Carlberg
R. G.
,
2014
,
ApJ
,
785
,
137

Sellwood
J. A.
,
Kahn
F. D.
,
1991
,
MNRAS
,
250
,
278

Semczuk
M.
,
Łokas
E. L.
,
del Pino
A.
,
2017
,
ApJ
,
834
,
7

Toomre
A.
,
1964
,
ApJ
,
139
,
1217

Toomre
A.
,
1969
,
ApJ
,
158
,
899

Toomre
A.
,
1977
,
ARA&A
,
15
,
437

Toomre
A.
,
1981
, in
Fall
S. M.
,
Lynden-Bell
D.
, eds,
Structure and Evolution of Normal Galaxies: What Amplifies the Spirals
,
Cambridge University Press
,
Cambridge
, p.
111

Vaghmare
K.
,
Barway
S.
,
Mathur
S.
,
Kembhavi
A. K.
,
2015
,
MNRAS
,
450
,
873

Vauterin
P.
,
Dejonghe
H.
,
1996
,
A&A
,
313
,
465

Williams
A. A.
,
Evans
N. W.
,
Bowden
A. D.
,
2014
,
MNRAS
,
442
,
1405

Zhang
X.
,
1998
,
ApJ
,
499
,
93

Zhang
X.
,
1999
,
ApJ
,
518
,
613

APPENDIX A: GRAVITATIONALLY SOFTENED LINEAR STABILITY ANALYSIS

In short, an axially symmetric disc galaxy model is characterized by a distribution function F0(E, Jϕ), with E the specific binding energy and Jϕ the specific angular momentum of a stellar orbit, and a (positive) gravitational binding potential V0(r). pystab retrieves those complex frequencies ω for which a spiral-shaped perturbation of the form
(A1)
constitutes an eigenmode with an infinitesimal amplitude. Here, (r, ϕ) are polar coordinates in the stellar disc, Vpert(r) is a complex function quantifying the mode’s amplitude and phase, m is its multiplicity, Ωp = ℜ{ω}/m its pattern speed, and ℑ{ω} its growth rate. A general perturbing potential can always be expanded in such spirals and, for infinitesimal amplitudes, these can be studied independently from each other.
In essence, pystab solves the first-order collisionless Boltzmann equation to find the response distribution function fresp(r, ϕ, vr, vϕ, t) produced by a given perturbation Vpert(r, ϕ, t). This response distribution function generates a response density given by
(A2)
which, in turn, corresponds to a response gravitational potential given by
(A3)
with ψ the inter-particle interaction potential. For Newtonian gravity, one evidently has
(A4)
while Plummer softened interactions are quantified by the choice
(A5)
with ε the softening length.

Since we want to validate our approach by comparing particular results with the numerical simulations presented in Sellwood (2012), we mimic the strategies employed in that work, in particular the softening method. There, the initial condition is generated by sampling stellar particles from the distribution function F0(E, J) evaluated using the mean Newtonian gravitational potential V0(r), independent of the gravitational softening that is employed to evolve the particles through time. Moreover, the axially symmetric force field of the base state is evaluated correctly, i.e. without softening, while the non-axisymmetric force field of the growing waves is softened. Therefore, in our search for unstable models, we only implement Plummer-type softening in the response potential Vresp(r, θ, t), but not in the axially symmetric base state potential V0(r). Using this strategy, equation (A3) is the only place where the softened gravitational interaction enters the computation of the modes. This approach to including gravitational softening in linear stability analysis is detailed in De Rijcke et al. (2018).

Eigenmodes are identified by the fact that
(A6)
and pystab employs a matrix method (Kalnajs 1977; Vauterin & Dejonghe 1996) to find them. To do so, the perturbing potential Vpert is expanded in a basis of potentials, V. The response to each basis potential, denoted by |$V_{\ell ,\rm resp}$|⁠, can likewise be expanded in this basis as
(A7)
If the perturbation is an eigenmode, then the ω-dependent |${\mathcal {C}}$| matrix can be shown to possess a unity eigenvalue (Vauterin & Dejonghe 1996). This feature is exploited by pystab to identify the eigenmodes. Since it is assumed that the mode amplitude is zero at time t = −∞, we only consider growing modes, i.e. with frequencies with a positive imaginary part.

The formalism contains a number of technical parameters, such as the number of orbits on which phase space is sampled [here we use norbit(norbit + 1)/2 orbits with norbit = 600 in the allowed triangle of turning points – or pericentre/apocentre – space], the number nFourier of Fourier components in which the periodic part of the perturbing potential is expanded (here we use nFourier = 80), the number of potential-density pairs (PDPs) that is used for the expansion of the radial part of the perturbing potential and density (we typically use 40 PDPs), and the shape and extent of the PDP density basis functions. Here, we use PDP densities Σ with compact radial support. They cover the relevant part of the stellar disc and are evenly spaced on a logarithmic scale so the resolution is highest in the inner regions of the disc. Their radial widths are automatically chosen such that consecutive basis functions are unresolved and can be used to represent any smooth radial function. The corresponding PDP potentials are computed numerically using equation (A3).

Author notes

Hubble fellow.

This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://dbpia.nl.go.kr/journals/pages/open_access/funder_policies/chorus/standard_publication_model)