The shearing sheet and swing amplification revisited

TL;DR

This paper revisits the classic shearing sheet and swing amplification analysis, providing clearer derivations and new insights into stellar disc dynamics, including wavepacket behavior and limitations of existing dispersion relations.

Contribution

It offers a more accessible derivation of key results, explores the group velocity and wavepacket structure, and highlights limitations of the Lin-Shu-Kalnajs dispersion relation in disc dynamics.

Findings

Two wavepackets emerge inside corotation, not one.

Disturbances persist longer in velocity space than in real space.

Limitations of the Lin-Shu-Kalnajs dispersion relation are demonstrated.

Abstract

The principal results of the classic analysis of the shearing sheet and swing amplification by Julian & Toomre (1966) are re-derived in a more accessible way and then used to gain a better quantitative understanding of the dynamics of stellar discs. The axisymmetric limit of the shearing sheet is derived and used to re-derive Kalnajs' 1965 dispersion relation and Toomre's 1964 stability criterion for axisymmetric disturbances. Using the shearing sheet to revisit Toomre's important 1969 paper on the group velocity implied by Lin-Shu-Kalnajs dispersion relation, we discover that two rather than one wavepackets emerges inside corotation: one each side of the inner Lindblad resonance. Although LSK dispersion relation provides useful interpretations of both wavepackets, the shearing sheet highlights the limitations of the LSK approach to disc dynamics. Disturbances by no means avoid an annulus around corotation, as the LSK dispersion relation implies. While disturbances of the shearing sheet have a limited life in real space, they live on much longer in velocity space, which Gaia allows us to probe extensively. C++ code is provided to facilitate applications of winding spiral waves.

Equations

$${\cal L}={\textstyle{1\over 2}}\biggl{[}\dot{x}^{2}+(R+x)^{2}\biggl{(}\Omega+{\dot{y}\over R}\biggr{)}^{2}\biggr{]}-\Phi(R+x)$$

$$p_x$$

$$p_y$$

$$\simeq R\Omega +2\Omega x+\dot{y}.$$

$$H$$

$${\Delta }_y\equiv p_y-R\Omega =2\Omega x+\dot{y}.$$

$$\frac{\partial \Phi }{\partial x}$$

$$\frac{{\partial }^2\Phi }{\partial x^2}$$

$$=\Omega \left(\Omega -4A\right),$$

$$\Phi (R+x)\simeq \Phi (R)+R{\Omega }^{2}x+\frac{1}{2}\Omega (\Omega -4A)x^2,$$

$$H$$

$$+R{\Omega }^{2}x+\frac{1}{2}\Omega (\Omega -A)x^2$$

$$\simeq \frac{1}{2}\left[p_x^2+{\Delta }_y^2-R^2{\Omega }^2\right]+\Phi (R)-x\Omega {\Delta }_y+\frac{1}{2}{\kappa }^2{x}^2$$

$$\overline{x}\equiv 2\Omega {\Delta }_y/{\kappa }^2.$$

$$A$$

$${\kappa }^2$$

$$\dot{y}$$

$$=-2Ax(circular\; orbit).$$

$$v_{\varphi }=2(A-\Omega )x+\frac{{\kappa }^2\overline{x}}{2\Omega }=2B(x-\overline{x}),$$

$$x$$

$$v_{\varphi }$$

$$\dot{y}$$

$$\Rightarrow y(t^{\prime })$$

$$=y_0+{\Delta }_y\frac{A}{B}{t}^{\prime }-2\frac{\Omega }{\kappa }X\sin {\theta }_r,$$

$$H$$

$$=\frac{1}{2}\left[p_x^2+{\Delta }_y^2\left(1-\frac{4{\Omega }^2}{{\kappa }^2}\right)-R^2{\Omega }^2\right]+\Phi (R)+\frac{{\kappa }^2}{8B^2}{v}_{\varphi }^2$$

$$=\frac{1}{2}\left[p_x^2+{\Delta }_y^2\frac{A}{A-\Omega }-R^2{\Omega }^2\right]+\Phi (R)+\frac{\Omega }{-2B}{v}_{\varphi }^2$$

$$=H_x(p_x,v_{\varphi })+H_y({\Delta }_y),$$

$$H_x\equiv \frac{1}{2}\left(p_x^2+\frac{\Omega }{-B}{v}_{\varphi }^2\right)=\frac{1}{2}{\kappa }^2{X}^2.$$

$$H_x=\frac{1}{2}\left[p_x^2+{\kappa }^2(x-\overline{x}{)}^2\right].$$

$$f_0=F{\mathrm{e}}^{-H_x/{\sigma }^2},$$

$${\Sigma }_0=\int \mathrm{d}{p}_x\mathrm{d}{v}_{\varphi }{f}_0=\frac{4\pi {\sigma }^{2}BF}{\kappa }.$$

$$k_{\mathrm{crit}}\equiv \frac{2\pi }{{\lambda }_{\mathrm{crit}}}=\frac{{\kappa }^2}{2\pi G{\Sigma }_0}.$$

$$Q\equiv \frac{\kappa \sigma }{3.36G{\Sigma }_0},$$

$$\omega =k_y\dot{y}=-2Ax{k}_y.$$

$$x_{\mathrm{L}}=\pm \frac{\kappa }{2A{k}_y}.$$

$${\Sigma }_1(\mathbf{x},t)={\tilde{\Sigma }}_1(t){\mathrm{e}}^{\mathrm{i}\mathbf{k}\cdot \mathbf{x}},$$

Full text

1 Introduction

Spiral structure has fascinated astronomers since its discovery two centuries ago. Although our quantitative grasp of this phenomenon remains inadequate, we can now claim a good qualitative understanding of the phenomenon, and we know that it plays a key role in the dynamical and chemical evolution of galaxies like ours (e.g. Aumer et al., 2016b).

The key ingredients of spiral structure are: the Lin-Shu-Kalnajs (LSK) dispersion relation111This relation, first given by Lin & Shu (1966), generalises the relation obtained for axisymmetric disturbances by Kalnajs (1965). for running density waves; the swing amplifier (Goldreich & Lynden-Bell, 1965; Julian & Toomre, 1966; Toomre, 1981); and resonant absorption at Lindblad resonances (Lynden-Bell & Kalnajs, 1972; Sellwood & Carlberg, 2014). The picture is as follows. Noise from any source will contain a packet of leading spiral waves. The LSK dispersion relation is such that the packet will travel away from the nearest Lindblad resonance towards corotation, unwinding as it goes (Toomre, 1969). Eventually the tight-winding approximation on which the LSK relation is based fails, and following Julian & Toomre (1966, hereafter JT66) the packet’s evolution must be understood in the context of the shearing sheet. As the packet swings from leading to trailing, it is amplified, and two stronger packets of trailing waves propagate away from corotation towards the Lindblad resonances. Should any part of these waves be reflected back towards corotation as leading waves, repeated swing-amplification can permit insignificant noise to grow into a manifest trailing spiral pattern and ultimately a strong bar.

One way in which a trailing wave moving away from corotation can convert into a leading wave that returns to corotation, is refraction around the galactic centre (Toomre, 1981). However, this will not occur if a Lindblad resonance is encountered before the centre is reached because the wave will be resonantly absorbed there (Lynden-Bell & Kalnajs, 1972) rather than proceeding to the centre. Moreover, when the circular-speed curve is rather flat, any wave will have a Lindblad resonance between corotation and the centre.

Sellwood & Carlberg (2014), however, pointed out when a wave is resonantly absorbed at its Lindblad resonance, a barrier is formed that is liable to reflect some part of any subsequent wave as it approaches the barrier en route to its Lindblad resonance. This is because absorption of the wave permanently modifies the distribution function (DF) in a narrow zone around the resonance, and it is generically the case that waves are partially reflected when the impedance of the medium in which they are travelling changes in a distance smaller than their wavelength. Hence, if initially a disc’s DF is a smooth function of the actions, the first wavepacket reaches its Lindblad resonance, where it is resonantly absorbed by Landau damping. As a result of this absorption, the DF develops large gradients in the vicinity of the Lindblad resonance and when a subsequent packet of swing-amplified waves attempts to pass through this region of modified DF, it is partially reflected back towards corotation, where it is again swing amplified. Consequently this second wave packet modifies the DF more strongly at its Lindblad resonances than did the first packet, and the scope for reflection back to corotation for repeated amplification grows as the disc ages (Sellwood & Carlberg, 2014). Eventually the effectiveness of swing amplification is such that the spiral structure becomes an order-unity phenomenon and refashions the centre of the disc into a bar. The bar may later buckle into a bulge/bar (Combes & Sanders, 1981; Raha et al., 1991).

The swing amplifier and resonant absorption are the stand-out pieces of physics in this beautiful mechanism by which galaxies like ours evolve. The first aim of this paper is to present an accessible account of the swing amplifier, which Binney & Tremaine (2008) shied away from attempting. They did so because the derivation in JT66 is subtle and hard to follow. It relies on the theory of characteristics, which will be unfamiliar to most students, and employs a bewilderingly large number of symbols and changes of variables. This paper re-derives the key results of JT66 in a way that will be more accessible to the average member of the community. With this pedagogic aim in view, more intermediate steps in the algebra are given than are strictly necessary.

The second aim of the paper is to illustrate the value of the key result in JT66 by plotting a variety of figures that help us to understand how stellar discs work. Many of these figures appeared already in either JT66, Toomre (1969) or Toomre (1981), but with limited explanation of the computational details.

Section 2 presents the shearing sheet from a Hamiltonian perspective. Section 3 sets up the key equation by linearising the collisionless Boltzmann equation (CBE). Section 3.1 computes the signature of a winding spiral in velocity space. Section 3.2 derives the equation that governs axisymmetric disturbances. This simpler equation admits modes and their dispersion relation is derived and show to be identical with the axisymmetric limit of the LSK dispersion relation. Section 4 uses the key equation to investigate the response of the disc to an impulsive jolt. Section 5 uses it to determine the increase in density around an orbiting mass that arises because the disc is highly ‘polarisable’, and Section 6 uses it to study the dynamics of wavepackets. Section 7 discusses the insights we gain from the formalism and its limitations. Section 8 sums up. Appendix A evaluates a required Gaussian integral, while Appendix B describes a C++ class that computes the dynamics of winding waves.

2 Shearing sheet

We consider waves characterised by a large value of the azimuthal quantum number $m$ and study a small near-Cartesian patch that moves around the disc at the frequency $\Omega$ of a circular orbit at the mean radius $R$ of the patch. We define a radial coordinate $x$ by taking a general radius to be $R+x$, and we define a complementary coordinate by $y\equiv R\phi$, where $\phi$ is the angle between the lines from the disc’s centre to the patch’s centre and the point $(x,y)$ (Fig. 1). In this rotating frame, the Lagrangian is

[TABLE]

because in an inertial frame the angular velocity of the point at $R+x$ is $\Omega+\dot{\phi}=\Omega+\dot{y}/(R+x)\simeq\Omega+\dot{y}/R$. From ${\cal L}$ we can read off the momenta:

[TABLE]

It follows that the Hamiltonian is

[TABLE]

Since $H$ doesn’t contain $t$ or $y$, we have two constants of motion: $H$ and $p_{y}$ or

[TABLE]

We will use $\Delta_{y}$ rather than $p_{y}$ because, unlike $p_{y}$, it is first-order in the small quantities $x,y$.

The derivatives of $\Phi$ at the origin are

[TABLE]

where $A$ is Oort’s first constant (Table 1). Hence

[TABLE]

so

[TABLE]

where $\kappa^{2}\equiv 4\Omega(\Omega-A)$. From the way $x$ and $p_{x}$ appear in $H$ if follows that $x$ oscillates harmonically at the epicycle frequency $\kappa$ about

[TABLE]

2.1 Orbits

Circular orbits are ones on which $x=\overline{x}$, so from equation (6)

[TABLE]

This relation describes how the sheet shears in the absence of random velocities.

Let $v_{\phi}\equiv 2Ax+\dot{y}$ be the azimuthal speed relative to the local circular orbit. Then from equation (14) and the definition (6) of $\Delta_{y}$ we have

[TABLE]

where $B$ is Oort’s second constant (Table 1). Hence $v_{\phi}$ tells us how far a star is from its guiding centre. Given that $x$ oscillates harmonically, the coordinates of any star can be written

[TABLE]

where $\theta_{r}=\kappa t+\hbox{constant}$. Similarly,

[TABLE]

where $y_{0}$ is a constant of integration.

When we eliminate $x$ from $H$ in favour of $v_{\phi}$ we get

[TABLE]

where

[TABLE]

Since $B<0$ and $p_{x}$ is the radial component of velocity, in velocity space lines of constant $H_{x}$ form ellipses. From equation (11) we obtain an alternative expression for $H_{x}$:

[TABLE]

2.2 Characteristic scales

By Jeans’ theorem a suitable equilibrium DF is

[TABLE]

where $F$ and $\sigma$ are free constants. By equation (28) this DF generates a biaxial Maxwellian velocity distribution with dispersion $\sigma$ in $p_{x}$ and $\sqrt{-B/\Omega}\sigma$ in $v_{\phi}$. The surface density it generates is

[TABLE]

It is convenient to specify the constants $F$ and $\sigma$ through the values of two functions of them. Toomre (1964) showed that even a stellar disc in which all stars are on perfectly circular orbits (so $\sigma=0$) is stable to axisymmetric perturbations with wavenumbers smaller than the critical wavenumber

[TABLE]

We use $k_{\rm crit}$ as proxy for the surface density $\Sigma_{0}$ or, equivalently, the normalisation of the DF, $F$. Our proxy for $\sigma$ is

[TABLE]

which is the ratio of the actual velocity dispersion $\sigma$ to the minimum value that ensures that the disc is stable to axisymmetric disturbances of any wavenumber (Toomre, 1964).

If the disc is disturbed by a potential $\Phi_{1}\propto\cos(k_{y}y)$, a star that is on a circular orbit at $x$ experiences the disturbance at frequency

[TABLE]

At the locations $\pm x_{\rm L}$ of the Lindblad resonances, the absolute value of this frequency coincides with the star’s radial frequency $\kappa$. Clearly the Lindblad resonances of this disturbance lie at

[TABLE]

If we adopt plausible values $\Sigma_{0}=49\,{\rm M}_{\odot}\,\mathrm{pc}^{-2}$ and $\kappa=37\,\mathrm{km\,s}^{-1}\,\mathrm{kpc}^{-1}$ for the solar neighbourhood (Binney & Tremaine, 2008), then $\lambda_{\rm crit}\simeq 6\,\mathrm{kpc}$ and $k_{\rm crit}\sim 1\,\mathrm{kpc}^{-1}$.

2.3 Shearing pattern

We assume that our coordinate patch is centred on the corotation radius of a spiral wave as sketched on the right of Fig. 2. Hence the pattern shears with particles on circular orbits. That is, if the perturbed surface density has the form

[TABLE]

the phase ${\bf k}\cdot{\bf x}$ is constant at the location of every particle that moves on a circular orbit. For this to be so, ${\bf k}$ has to be a function of time. We determine that function by equating phases at a particle’s locations at two different times:

[TABLE]

where we have used equation (16). Henceforth we take the origin of time to be the instant

[TABLE]

at which $k_{x}=0$. At times $t<0$, $k_{x}<0$ and the spiral is leading, while at $t>0$ it is trailing. In general

[TABLE]

and

[TABLE]

goes through a minimum as $k_{x}$ passes through zero.

So long as $|k|\ll R$, the gravitational potential generated by a sinusoidal density variation in a razor-thin disc is (e.g. Binney & Tremaine, 2008)

[TABLE]

The strength of this potential peaks as $k_{x}$ passes through zero.

Using our solution (19) and (21) for a general orbit we now have that

[TABLE]

If we define the phase

[TABLE]

this simplifies to

[TABLE]

It follows that if $t$ and $t^{\prime}$ are two instants along a given orbit, then

[TABLE]

3 Linearizing the CBE

With $f=f_{0}(H_{0})+f_{1}$ and $H=H_{0}+\Phi_{1}$, the CBE

[TABLE]

where $[,]$ is a Poisson bracket, linearises to

[TABLE]

On integration we have

[TABLE]

where the integral is along unperturbed orbits and $f_{1}$ has been assumed to vanish for $t^{\prime}<t_{\rm i}$.

The unperturbed DF (30) is a function of not only $p_{x},x$ but through $\overline{x}$ or $v_{\phi}$ also of the momentum $\Delta_{y}$. Using equation (29) for $H_{x}$ with $\Phi_{1}({\bf x})\propto\mathrm{e}^{{\rm i}{\bf k}\cdot{\bf x}}$, we have

[TABLE]

Since we have assumed that $f_{1}$ vanishes at $t_{\rm i}$, it will vanish at all subsequent time unless something contributes to $\Phi_{1}$ in addition to the surface density associated with $f_{1}$. Hence, we take

[TABLE]

where the second equality uses equation (48), $\widetilde{\Sigma}_{\mathrm{e}}(t^{\prime})\mathrm{e}^{{\rm i}{\bf k}\cdot{\bf x}|_{t^{\prime}}}$ is a fictitious surface density that provides a gravitational driving force, and $\widetilde{\Sigma}_{1}$ is defined by

[TABLE]

In view of equations (51) to (53), with equations (19) and (41) we have

[TABLE]

Next we have to compute $\widetilde{\Sigma}_{1}$ by integrating $f_{1}$ over ${\bf v}$, but before integrating we eliminate $X,\theta_{r}$ in favour velocity components $p_{x},v_{\phi}$ at $t^{\prime}=0$. Using $\theta_{r}=\theta_{0}+\kappa t^{\prime}$ we write

[TABLE]

where

[TABLE]

Similarly

[TABLE]

The Jacobian of the (time-dependent) transformation in velocity is simply

[TABLE]

With ${\bf U}$ replacing $(p_{x},v_{\phi})$ the phase factor of the integrand in equation (56) becomes

[TABLE]

where $C(t)\equiv\cos\kappa t$, $S(t)\equiv\sin\kappa t$ and $C^{\prime}$ and $S^{\prime}$ are defined the same way but with $t^{\prime}$ replacing $t$.

Using equation (31) to eliminate $F$ from equation (56), we see that the coefficient of $\mathrm{e}^{{\rm i}{\bf k}\cdot{\bf x}}$ in the expression for the perturbed surface density can be written

[TABLE]

where

[TABLE]

is a dimensionless kernel that involves the Gaussian integral

[TABLE]

The integral is of the form

[TABLE]

where

[TABLE]

Hence

[TABLE]

In terms of the critical wavenumber (32), $I$ simplifies to

[TABLE]

Now we replace ${\bf b}$ by the dimensionless vector

[TABLE]

and note that

[TABLE]

where we have introduced $Q$ (eqn. 33). With equations (77) and (79), the kernel (68) can be written

[TABLE]

The vectors $\widehat{\bf b}$ and ${\bf c}$ are functions of the dimensionless numbers $At$, $At^{\prime}$, $\kappa t$, $\kappa t^{\prime}$ and $\Omega/\kappa$, while $\widehat{\bf b}$ has an additional (linear) dependence on $k_{y}/k_{\rm crit}$. The argument $(Q\widehat{b})^{2}$ of the exponential is proportional to $\sigma^{2}$ and has no dependence on $\Sigma_{0}$.

In the following we shall refer to equation (66) as the JT equation and $K$ as the JT kernel. Although JT66 normalised their kernel differently, equation (66) is equivalent to their equation.

3.1 Impact on velocity space

Once we have computed the evolution of $\widetilde{\Sigma}_{1}$ from the JT equation, we can use equation (56) to compute the evolution of $f_{1}$. Since we can directly observe the phase space of our Galaxy, $f_{1}({\bf x},{\bf v},t)$ contains observationally significant information that was erased when we integrated over ${\bf U}$ to compute $\widetilde{\Sigma}_{1}$.

Equation (56) is explicit about the ${\bf x}$ dependence of $f_{1}$, but its ${\bf v}$ dependence is buried in the variables $X$ and $\theta_{r}$. Above we expressed $X$ and $\theta_{r}$ as functions of the velocity at the special moment when $k_{x}=0$. Now we want to express them as functions of the velocity at the current time $t$. We do this by now writing $\theta_{r}=\theta_{0}+\kappa(t^{\prime}-t)$, where $\theta_{0}$ is a star’s current radial phase. With this new definition of $\theta_{0}$,

[TABLE]

where ${\bf U}$ is related to the current velocity components by equation (60). Equation (64) becomes

[TABLE]

where $C^{\prime}=\cos[\kappa(t^{\prime}-t)]$ and $S^{\prime}=\sin[\kappa(t^{\prime}-t)]$. Completing the transition from $X,\theta_{r}$ to ${\bf U}$, we find

[TABLE]

3.2 Axisymmetric limit

We now recover the form of the JT equation that describes the response of a disc to axisymmetric disturbances. That is, we take the limit $k_{y}\to 0$ of the JT equation. Since $\widehat{\bf b}\propto k_{y}$ and both the prefactor and the exponent in equation (80) vanish with $k_{y}$, superficial analysis yields the conclusion that in this limit $K\to 0$, making the JT equation useless. But we have to bear in mind that $k_{x}=2Atk_{y}$ and we require $k_{x}\neq 0$, so we must let $|t|\to\infty$ as we let $k_{y}\to 0$ such that the product $tk_{y}$ is constant. By allowing $t$ to diverge we are recognising that the instant when $k_{x}=0$ and the wave crests run radially lies in the remote past (or future if $k_{x}<0$).

A glance at equations (72) shows that in this limit $k_{y}\to 0$ but $k_{y}t=\hbox{constant}$, ${\bf b}$ and ${\bf c}$ simplify such that

[TABLE]

where $k_{x}^{\prime}\equiv k_{x}(t^{\prime})$, etc. Now as $k_{y}\to 0$ ${\rm d}k_{x}/{\rm d}t\to 0$, so over any finite interval $t-t^{\prime}$ there will be negligible change in $k_{x}$. Moreover, while we are obliged to let $|t,t^{\prime}|\to\infty$, it suffices to consider forces that acted a finite time in the past. That is, we need to evaluate $K(t,t^{\prime})$ for finite $|t-t^{\prime}|$, so we can neglect the difference between $k_{x}$ and $k_{x}^{\prime}$. Simplifying the right sides of equations (87) thus and using the trigonometrical identities $CC^{\prime}+SS^{\prime}=\cos[\kappa(t-t^{\prime})]$ and $SC^{\prime}-CS^{\prime}=\sin[\kappa(t-t^{\prime})]$, we find that in the axisymmetric limit $k_{y}\to 0$, the JT kernel reduces to (JT66)

[TABLE]

where

[TABLE]

The JT equation now reads

[TABLE]

We replace $t^{\prime}$ by $t^{\prime\prime}\equiv t-t^{\prime}$. As $t^{\prime}$ goes from $-\infty$ to $t$, $t^{\prime\prime}$ goes from $\infty$ to 0. Hence the equation becomes

[TABLE]

The right side has the form of a Laplace convolution of ${\cal K}$ with the sum of surface densities, so the Laplace transforms $\overline{\Sigma}_{1}(p)$, etc are related by

[TABLE]

Values of $p$ for which $\kappa\overline{{\cal K}}=1$ are of particular interest: for these frequencies the JT equation has a non-zero solution in the absence of a stimulating density $\Sigma_{\mathrm{e}}$. That is, these are the frequencies of the sheet’s axisymmetric normal modes.

Using the identity (Abramowitz & Stegun, 1965, eqn. 9.6.34)

[TABLE]

it is easy to show that the Laplace transform of ${\cal K}(\tau)$ is

[TABLE]

When $n=-|n|$ the second term in the big bracket above takes the value that the first term takes when $n=|n|$, so we can drop the second term and double the remaining sum:

[TABLE]

We relate this to the standard LSK dispersion relation by defining the dimensionless frequency $s\equiv{\rm i}p/\kappa$. In terms of $s$ the dispersion relation $\kappa\overline{{\cal K}}(p)=1$ can be written

[TABLE]

When the LSK dispersion relation is written in analogous form (Binney & Tremaine, 2008, eqn. K.25), and divided by $1-s^{2}$, one obtains

[TABLE]

The recurrence relation

[TABLE]

allows us put the LSK relation into the form

[TABLE]

Gathering together the coefficients of $I_{n}$ makes the sum precisely the sum that appears in equation (96), so the same manoeuvre to simplify the coefficient of $I_{n}$ establishes that the LSK dispersion relation can be written in the form (99) that was obtained from the JT kernel.

Passing to the axisymmetric limit removes all meaning to the concepts ‘pattern speed’ and ‘corotation’. Before we take this step, $x=0$ is singled out as the location at which stars move at a fixed phase in the perturbation. After going axisymmetric, every value of $x$ is equivalent. Hence the axisymmetric modes of sheets are standing sinusoidal waves of unlimited extent in $x$.

When the right-hand side of equation (99) is plotted as a function of $k_{x}/k_{\rm crit}$ at a fixed values of $\sigma$ and $s^{2}$, one obtains a curve that has a minimum at $k/k_{\rm crit}\sim 2$. The smaller the value of $\sigma$ or the larger the value of $s^{2}$, the lower the curve reaches. If the curve crosses the line $y=0$, the system has a mode of frequency $s$ with the associated value of $k_{x}/k_{\rm crit}$. The system is stable if the curve crosses $y=0$ only for $s^{2}>0$, and Toomre’s stability criterion $Q>1$ emerges by finding the value of $\sigma$ at which the curve just touches $y=0$ for infinitesimal positive $s^{2}$. The sheet’s modes oscillate indefinitely if stable ($s^{2}>0$) or they grow exponentially when unstable: there are no overstable modes.

4 Impulsive excitation of the disc

For our first application of the JT equation we follow JT66 by exciting the disc impulsively. We take

[TABLE]

where $\Sigma_{\rm s}$ sets the magnitude of the impulse and the denominator $\kappa$ ensures dimensional soundness. With this choice, the JT equation yields

[TABLE]

Equation (104) is straightforwardly solved by establishing a grid of about a thousand equally spaced times in the interval $(t_{\rm i},5\pi/\kappa)$, approximating the integral over $t^{\prime}$ by the trapezium rule on this grid, and then for each possible end time $t$ solving the resulting algebraic equation for $\widetilde{\Sigma}_{1}(t)$.

Fig. 3 shows solutions obtained in this way for four values of $t_{\rm i}$ in a disc with a perfectly flat rotation curve, $Q=1$ and $k_{y}=k_{\rm crit}$. Since the applied stimulus perturbs velocities but not positions, time is required for an overdensity to develop that can be amplified as $k_{x}$ passes through zero at $t=0$. Hence earlier excitation produces a bigger response up to times later than $t_{\rm i}\simeq-3\pi/2\kappa$. However, moving $t_{\rm i}$ further back decreases the peak in $\widetilde{\Sigma}(t)$ because the disturbance then has time to phase mix away prior to $t=0$. As a reflection of this tradeoff between phase mixing and the need for time for fluctuations in density to emerge, in Fig. 3 the red dotted curve for $\kappa t_{\rm i}=-2\pi$ almost exactly coincides with the black, long-dashed curve for $\kappa t_{\rm i}=-\pi$.

The black curves in Fig. 4 show, for a fixed value of $k_{y}$, the evolution of $\widetilde{\Sigma}_{1}$ for $t_{\rm i}=-3\pi/2\kappa$ and three values of $Q$, namely $1.1,\,1.2$ and $1.3$. As expected, as self-gravity becomes less important, the amplitude of the response declines steeply. Nonetheless, self-gravity greatly enhances the disc’s response even for the largest plotted value of $Q$. This fact is established by the red curve in Fig. 4, which plots $K(t,t_{\rm i})$ for $Q=1.3$, which is what the response would be in a disc of test particles.

The full black and dashed red curves in Fig. 5 show, respectively, the responses with and without self-gravity for the case $Q=1.5$ and $k=k_{\rm crit}/2$. In this warmer sheet and at this longer wavelength, swing amplification is less powerful and the complex responses can be studied on a common scale. Notwithstanding the comparative weakness of self-gravity in this case, it dramatically modifies the structure of the response, not just its amplitude. In the absence of self-gravity, the response consists of a quickly achieved maximum, followed by series of blips, each of which comprises a minimum quickly followed by a maximum. The first blip is centred on $\kappa t=5\pi/4$ and the second on and $10\pi/3$. With self-gravity, the response comprises a rapid rise to a plateau at a value slightly larger than the peak response achieved without self-gravity, followed by a steep rise to a high peak at $\kappa t=5\pi/4$ that commences just after $k_{x}$ passes through zero. This peak is then rapidly followed by an almost double-bottomed minimum. The self-gravitating response begins to rise out this minimum at essentially the same moment that the response without self-gravity starts to rise from its second minimum. The self-gravitating response has a second maximum at $\kappa t\simeq 4\pi$, when without self-gravity the response is essentially zero.

The key feature of Fig. 5 is the strongly anharmonic nature of the response in the absence of gravity: Fourier analysis of the dashed red curve would clearly show significant power at many frequencies. Given that the JT kernel $K(t,t^{\prime})$ is obtained in the epicycle approximation, in which all orbits are perfectly harmonic, and is driven by a perfectly sinusoidal gravitational field, the anharmonic nature of the response is remarkable. It arises through the integration of the sinusoidal gravitational field along unperturbed orbits that carry stars in and out and around in azimuth at a rate that speeds up or slows down as $x$ decreases or increases. The self-gravitating response is less anharmonic, but the complex structure of its minimum strongly suggests an interference pattern.

Fig. 6 plots for seven values of $Q$ the amplification factor

[TABLE]

versus the lengthscale $\lambda/\lambda_{\rm crit}$. Here the outer max operator involves a search over the time $t_{\rm i}$ at which the disc is jolted, while the inner operators range over the times at which the magnitudes of the responses are measured. The value of $t_{\rm i}$ that produces the maximum response depends strongly on $\lambda/\lambda_{\rm crit}$. If the latter is significantly less than unity, the disturbance phase mixes rapidly, so the largest response arises when the disc is jolted shortly before $k_{x}$ passes through zero, and $\kappa t_{\rm i}<\pi/2$. If $\lambda/\lambda_{\rm crit}\ga 1$, phase mixing is slower and the largest response is obtained by jolting the disc early on, so there is time for a significant overdensity to emerge before $k_{x}$ passes through zero. This is especially true if $Q-1$ is small. Since maximising the response involves arranging for an oscillating system to have a favourable phase at a particular instant, the optimum value of $t_{\rm i}$ is not a continuous function of $Q$ and $\lambda/\lambda_{\rm crit}$. This fact gives rise to kinks and bumps in the curves of Fig. 6 that might be thought indicatons of numerical error but are not.

Fig. 6 closely resembles the central panel of Fig. 7 in Toomre (1981) although the quantity it plots is quite different. Toomre (private communication) plotted the ratio of the rms responses obtained when the disc was stimulated by leading and trailing waves with wavevectors that differed only in the sign of $k_{x}$. The wavevector of the leading wave would later become identical to the initial wavevector of the trailing wave, but not before it had been swing amplified. Thus Toomre’s ratio is a very clean measure of the effectiveness of the amplifier. By plotting the ratio of the peak responses with and without self-gravity Fig. 6 is a measure of the effectiveness of self-gravity. The similarity of the two figures indicates that the principal importance of self-gravity is in driving the swing amplifier.

Fig. 6 shows that at every value of $Q$, the amplification vanishes with the lengthscale as a consequence of fast phase mixing, then rises to a peak $\lambda/\lambda_{\rm crit}\simeq 2$ before falling to a plateau that extends from $\lambda/\lambda_{\rm crit}\simeq 3$ to the longest lengthscales. Discontinuities in the optimum value of $t_{\rm i}$ sometimes generate bumps in the plateau.

In N-body simulations of discs with sustained star formation, $Q$ settles to a value around $1.4$ (e.g. Aumer et al., 2016a). Fig. 6 indicates that in such a disc perturbations are swing-amplified by a factor that hovers around 14 for $1.3<\lambda/\lambda_{\rm crit}<2.1$.

Fig. 7 shows the evolution of the perturbation to velocity space in a sheet with $Q=1.2$ when a wave that is stimulated by the surface density (103) at $t=-3\pi/2\kappa$ is swing amplified. Successive panels show velocity space at intervals $\pi/\kappa$ apart, starting with $t=0$, when $k_{x}=0$. The dotted curve in Fig. 4 shows the evolution of this wave’s density in real space. This peaks before $t=3\pi/2\kappa$ and has long vanished by $t=4\pi/\kappa$. The central panel in the bottom row of Fig. 7 shows that at this late time the velocity-space signature of the wave is still growing. That is, the wave persists for much longer in velocity space than it does in real space.

In velocity space the wave manifests as a complex pattern of maxima and minima that is constantly shrinking in scale velocity while increasing in amplitude at a declining rate. The pattern varies with the location ${\bf x}$ of the plotted velocity space, becoming less symmetric as one moves away from corotation. In general it looks like a network of cells.

5 Response to a cloud

We now compute the response of a stellar disc to the insertion at time $t_{\rm i}$ of a mass such as a molecular cloud that moves on a circular orbit. The mass provides an endless succession of the $\delta$-function stimuli we considered in the last section.

Previously we considered stimuli that had a well defined circumferential wavenumber $k_{y}$. Now we are considering a succession of broad-band stimuli since when we Fourier transform the surface density

[TABLE]

we find that every wavenumber has non-vanishing amplitude:

[TABLE]

Notice that $\widehat{\Sigma}$ has dimensions of mass, unlike the quantity $\widetilde{\Sigma}_{1}(t)$ defined above as the coefficient of $\mathrm{e}^{{\rm i}{\bf k}\cdot{\bf x}}$ in the expression for $\Sigma_{1}({\bf x},t)$. This difference reflects the fact that to recover $\Sigma({\bf x})$ from $\widehat{\Sigma}({\bf k})\mathrm{e}^{{\rm i}{\bf k}\cdot{\bf x}}$ we have to integrate over ${\bf k}$. Fortunately, for given ${\bf k}$, $\widehat{\Sigma}_{1}$ will evolve in time in just the same way that $\widetilde{\Sigma}_{1}$ does, so in the JT equation (66) we can replace $\widetilde{\Sigma}_{1}$ by $\widehat{\Sigma}_{1}$. Then we have

[TABLE]

Fig. 8 shows solutions to equation (108) with $\Delta=0$ for $k_{y}=k_{\rm crit}/2$ and several values of the time $t_{\rm i}$ at which the mass $M$ is added. The corresponding values of $\kappa t_{\rm i}/\pi$ are given in the legend at top left. The quantity plotted is the amplitude of a swinging wave as a function of time, with the origin of time taken as the instant at which $k_{x}=0$, when the crests run radially. The curve with the smallest amplitude is that for the latest time of mass insertion $t_{\rm i}=0$. The other curves show that as the moment of insertion of the mass is pushed back, the amplitude of the wave grows until it reaches peak amplitude for $\kappa t_{\rm i}/\pi\sim 2.5$, and then settles to a steady value for even earlier insertion.

The mass $M$ will disturb the density in our coordinate patch with a superposition of waves like that shown by the red curve in Fig. 8 that differ in the times at which $k_{x}=0$. If the moment of mass insertion lies far in the past, every wave will be described by a curve that closely resembles the red curve in Fig. 8 but shifted to the right or left according as its crests are vertical later or earlier than the wave for which the red curve was computed. For this wave, $k_{x}=2Ak_{y}t$ (eqn. 40). For a wave that has $k_{x}=0$ at time $t_{x}$, the relation between $t$ and $k_{x}$ is different: $k_{x}=2Ak_{y}(t-t_{x})$. Whereas in Fig. 8 there is one value of $k_{x}$ at each time, at any time after insertion of the mass, the patch will contain waves with infinitely many values of $k_{x}$ on account of the presence of waves that swing through $k_{x}=0$ at every possible time.

Whereas in the section on impulsive excitation only one value of $k_{y}$ was in play, the mass excites waves with every value of $k_{y}$. The discussion of the last paragraph shows that from Fig. 8 we can infer the value of $\widehat{\Sigma}_{1}({\bf k})$ along a line of constant $k_{y}$ in the $(k_{x},k_{y})$ plane, and to determine $\widehat{\Sigma}_{1}({\bf k})$ over the the half plane $k_{y}>0$ we just need to solve equation (108) for positive values of $k_{y}$. Since $\widehat{\Sigma}_{1}({\bf k})$ is the Fourier transform of a real function, it satisfies

[TABLE]

and this relation gives the values we need in the other half plane.

Fig. 9 shows contours of constant overdensity

[TABLE]

for a cloud with characteristic size $\Delta=0.05\lambda_{\rm crit}$ in a disc with a flat rotation curve and $Q=1.4$. The contour values are in units of $M/(0.1\lambda_{\rm crit})^{2}$, so if the mass were spread at uniform density it would contribute one unit within a square $0.1\lambda_{\rm crit}$ on a side. Since the innermost contour is at overdensity 2 in these units and has an area that is roughly twice that of such a square, the excess stellar mass within this contour alone amounts to four times the cloud’s mass. From this fact it is clear that a stellar disc like that of our Galaxy is very polarisable: the effective mass of an object is several times its actual mass on account of its tendency to cause passing stars to linger in its vicinity.

6 Wave packets

We now excite the disc with the potential of the mass distribution

[TABLE]

The Fourier transform of the spatial part of this distribution is

[TABLE]

When this transform is used in the JT equation (66) with $\widetilde{\Sigma}$ replaced by $\widehat{\Sigma}$ as in the last section, both $\widehat{\Sigma}_{\mathrm{e}}$ and $\widehat{\Sigma}_{1}$ become explicit functions of ${\bf k}$. Hitherto the initial value of $k_{x}$ has been encoded via equation (40) in the variable $t_{\rm i}$. Now we treat $t_{\rm i}$ as a function of $k_{x}(0)$:

[TABLE]

As a consequence of this dependence of $t_{\rm i}$ on ${\bf k}$, the upper limit of the time integral in the JT equation loses its status as the current time. Instead we write

[TABLE]

where $\tau$ is the current time. When evaluating $\widehat{\Sigma}_{\mathrm{e}}({\bf k},\tau^{\prime})$ we use the current value of $k_{x}$:

[TABLE]

When we inverse Fourier transform $\widehat{\Sigma}_{1}({\bf k},\tau)$ to recover $\Sigma_{1}({\bf x},\tau)$ we likewise set $k_{x}$ to its value at time $\tau$.

In the previous applications of the JT equation, $\widehat{\Sigma}_{\mathrm{e}}$ has been real, with the consequence that $\widehat{\Sigma}_{1}$ has remained real. Now $\widehat{\Sigma}_{\mathrm{e}}$ is complex, so $\widehat{\Sigma}_{1}={\cal R}+{\rm i}{\cal I}$ also becomes complex. Using equation (109) we reason that

[TABLE]

where in the last line ${\cal R}+{\rm i}{\cal I}$ stands for the coefficient of $\pi\delta(k_{y}-k_{\rm i})$ in the transform $\widehat{\Sigma}_{1}({\bf k})$ (cf. eqn. 112).

Fig. 10 shows the result of exciting the disc with the external density (111) when $Q=1.2$, $\Delta=x_{\rm L}/\surd 2$ and $k_{\rm i}/k_{\rm crit}=\kappa/(4\pi A)=0.225$, where $x_{\rm L}$ is the distance of the wave’s Lindblad resonance from the corotation resonance (eqn. 35). Although the exciting density is centred on $x=-1.5x_{\rm L}$ (marked by a red vertical line), the excitation is concentrated between the Lindblad resonances because the external density corotates with particles at $x=0$, so stars that lie near the centre of the exciting density pass it rather rapidly, and are less strongly perturbed than more distant, but slower passing stars.

The envelope of excitation first grows in amplitude and then bifurcates into packets that propagate away from corotation inwards and outwards. The crests of the waves are essentially stationary, but they wax and wane in such a way that the packet moves quite coherently. That is, the phase velocity is much smaller than the group velocity. The smallness of the phase velocity is a direct consequence of the excitation being dominated by particles that nearly corotate with the exciting density, which is stationary in our reference frame. As the packet moves and fades, $k_{x}$ gradually increases. Toomre (1969) showed that the trajectories of the packets’ centres can be accurately predicted from the group velocity implied by the LSK dispersion relation.

If Fig. 10 a subsidiary packet develops inside the radius $x_{\rm i}$ (marked in red) on which the exciting potential is centred. The left edge of this packet lies near the value $x=-2x_{\rm L}$ at which stars on circular orbits perceive the wave at a frequency $\omega=2\kappa$. Our explanation of this feature relies on the structure of the axisymmetric limit of the JT kernel, which is the subject of the next section.

Smaller values of $k_{y}/k_{\rm crit}$ than applies in Fig. 10 yield similar wavepackets but ones that decay more slowly. The example shown in Fig. 11 has $k_{y}/k_{\rm crit}=0.15$. The largest packet moves briskly to the inner Lindblad resonance but seems to get stuck with about a quarter of its length over the resonance. In this stationary position, it slowly decays. Fig. 13 shows how this disturbance looks in the $xy$ plane at six of the times plotted in Fig. 11. The growth in the wave’s amplitude up to $t\sim 3P_{0}=6\pi/\kappa$ is evident, as is the steady increase in $k_{x}$. One can also see the tendency of both ingoing and outgoing packets to stick when they reach a Lindblad resonance.

Fig. 12 shows the effect of increasing $k_{y}$ to $k_{\rm crit}/2$. At this shorter azimuthal wavelength, the disc is significantly more responsive (Fig. 6) and the disturbance stimulated by the same mass is $\sim$twice as big. The stimulated spiral decays significantly faster: the response has almost extinguished by $t=5P_{0}$ rather than clearly persisting to $t=7P_{0}$ in Fig. 10 or to $t=10P_{0}$ in Fig. 11. The Lindblad resonances of this wave lie much closer to corotation than do the Lindblad resonances of the previous waves with lower $k_{y}$. In consequence the generated wavepacket fills the region between the Lindblad resonances, and, perhaps because it touches both Lindblad resonances from the outset, it shows little tendency to move.

6.1 Application of the LSK dispersion relation

In Section 3.2 we showed that the axisymmetric limit of the JT kernel gives rise to a dispersion relation (99) that is identical to the axisymmetric limit of the LSK dispersion relation. We now use this relation to elucidate the wavepackets evident in Figs. 10 and 11.

Fig. 14 plots the right-hand side of equation (99) as a function of $k_{x}$ for $Q=1.2$ and several values of $s$. The full curves correspond to $s=0.4,0.6,0.8$ and $0.9$, while the dashed curves correspond to $s=1.1,1.7,1.8$ and $1.9$. If the curve for some $s$ crosses the dashed line $y=0$, then there are modes with that frequency at the two values of $k_{x}$ at which the curve crosses $y=0$. On account of the denominators $(n+1)^{2}-s^{2}$ in equation (99), the curves change discontinuously as $s$ increases past an integer. Thus while the curve for $s=0.99$ (not plotted) drops far below $y=0$ and thus yields modes at both small and large $k_{x}$, the curve for $s=1.01$ (also not plotted) stays far above $y=0$ for all $k_{x}$. Thus we have two bands of $s$ for which modes exist: $1>s>0.41$ and $2>s>1.78$. Since the dispersion relation is a function of $s^{2}$, there are equivalent bands with $s\to-s$.

We now argue that a disc will respond to tightly wound spirals very much as it does to axisymmetric disturbance with the same value of $k_{x}$. On this understanding, we now consider

[TABLE]

to be a function of $x$, namely the frequency at which an $m$-armed spiral is perceived by a star on a near-circular orbit. The two bands in $s<0$ for which modes exist thus correspond to two bands in $x$ within which self-sustaining oscillations are possible. In Figs. 10 and 11, the band $-1<s<-0.41$ occupies territory to the right of the middle vertical line, which marks the inner Lindblad resonance, while the band $-2<s<-1.78$ occupies territory to the right of the left vertical line, which marks the harmonic of the Lindblad resonance. The band $-1<s<-0.41$ should extend 60% of the distance from the middle vertical line to $x=0$, while the band $-2<s<-1.78$ should extend 56% of the distance from the left black vertical line to the red line. In Figs. 10 and 11 the wavepackets somewhat exceed these bounds, but qualitatively the agreement is good.

Hence we can interpret the wavepackets seen in Figs. 10 and 11 as follows. The disturbance created by the external driving potential $\Phi_{\mathrm{e}}$ is largest and lives longest at locations where a self-sustaining mode is possible at the frequency at which stars on nearly circular orbits perceive $\Phi_{\mathrm{e}}$. At $x<0$ these locations are bounded on the left by $x=-2x_{\rm L}$ and $x=-x_{\rm L}$, while at $x>1$ they are bounded on the right by $x=x_{\rm L}$. Since $\Phi_{\mathrm{e}}$ has finite width in frequency space, modes with finite ranges in $k_{x}$ are excited. Interference between these modes first causes the region of largest net amplitude to migrate towards the bounding lines $x=\pm x_{\rm L}$, $x=\pm 2x_{\rm L}$, etc, and then causes the net disturbance to fade as different modes drift more and more out of phase and cancel ever more completely.

7 Discussion

A priori it is not obvious that in a warm, self-gravitating disc a perturbation with the initial form $\Sigma_{1}({\bf x})\propto\mathrm{e}^{{\rm i}{\bf k}\cdot{\bf x}}$ will evolve such that $\Sigma({\bf x})$ remains proportional to $\mathrm{e}^{{\rm i}{\bf k}\cdot{\bf x}}$ with ${\bf k}$ evolving so the phase is constant at particles that are on circular orbits. In our treatment this is established a posteriori by showing that we can construct a solution to the linearised CBE under the ansatz that ${\bf k}$ evolves by simple shear. JT66 established this fact more cleanly by Fourier transforming the CBE and then identifying its characteristics.

The existence of the JT66 solutions to the linearised perturbation problem clarifies the status of LSK running waves. The non-existence in a disc with $Q>1$ of modes with values of $s^{2}$ smaller than a threshold value is interpreted as an absence of waves in a region around corotation. As we have seen, the JT equation establishes that axisymmetric modes exist at all radii. The LSK dispersion relation has no solutions in a region around corotation precisely because any non-axisymmetric wave must wind up: the wavevector must be an explicit function of time, a possibility that one excludes a priori in the derivation of a dispersion relation. In light of this remark, the existence of solutions to the WKB dispersion relation away from corotation becomes puzzling. The puzzle is resolved by the dynamics of wavepackets, which Toomre (1969) showed move towards the Lindblad resonances. Consequently, a packet’s central ${\bf k}$ value becomes a function of time through the dependence of the LSK dispersion relation on radius. In this way the LSK theory manages after a fashion to encompass the growth in $k_{x}$ that is a categorical imperative in a shearing system.

Nevertheless, LSK waves do not play in discs a role analogous to Maxwell’s electromagnetic waves in electrodynamics. Electromagnetic waves are typically generated by shaking a charge, just as a solution of the JT equation can be generated by gravitationally jolting the sheet. Once generated, the waves of Section 4 are free oscillations of the system, exactly as a light wave is as it moves from emitting to absorbing atom. The standing wave in a laser is generally imagined to comprise counter-propagating travelling electromagnetic waves that are inter-converted the the laser’s end mirrors. If spiral structure can be understood in terms of normal modes, it will be by imagining the disc to be buzzing with JT66 winding waves, not LSK waves. The JT66 waves are simply as close as one can get to normal modes in a shearing sheet. True modes don’t exist because the general JT kernel is not time-translation invariant.

In the special case of axisymmetry, waves don’t wind up, so $k_{x}$ is time-independent and the JT equation admits normal modes. Their dispersion relation is identical with that of LSK waves with $m=0$.

We have explored the dynamics of wavepackets predicted by the JT equation by using the method Toomre (1969) devised to excite wavepackets comprising a coherent group of waves that share a common value of $k_{y}$. When the stimulating density has a small value of $k_{y}$ ($\la 0.5k_{\rm crit}$), packets form either side of corotation and move inwards and outwards to the Lindblad resonances just as the LSK dispersion relation predicts. If the stimulating density has a larger value of $k_{y}$, the sheet responds more energetically and in a less localised way. The stimulated structure decays quite rapidly without significant motion away from corotation. LSK theory provides no insight into this behaviour.

A feature of wave patterns with reasonable small values of $k_{y}/k_{\rm crit}$ that Toomre (1969) overlooked, is a wavepacket that fills the region at $-2<x/x_{\rm L}<-1$, which lies inside the Lindblad resonance at $x=-x_{\rm L}$. This newly identified wavepacket can be comprehended on the basis of the LSK dispersion relation to the same extent as the wavepacket in $-1<x/x_{\rm L}<0$ that Toomre (1969) first displayed. The new wavepacket is a nonetheless a curious phenomenon physically, because its surface density oscillates faster than can any of the stars from which it is built. Moreover, it is unlikely that any real galactic disc displays an excitation of this type because in a real disc the inner Lindblad resonance is much further removed from corotation than is the outer Lindblad resonance, and the analogue of $x=-2x_{\rm L}$ would lie at an exceedingly small radius. The shearing sheet, unlike a real disc, has resonances that are symmetrically arranged around corotation because the epicycle frequencies $\kappa$ of its particles are independent of $x$, whereas those of stars are roughly proportional to $1/R$.

Any wavepacket ultimately ceases to be visible in real space. We have seen, however, that in velocity space it lives on. However, because first-order perturbation theory does not encompass resonant trapping, it provides an inadequate account of the manner in which waves decay. N-body simulations (Sellwood, 2012) and calculations based on matrix mechanics (Fouvry et al., 2015) show that a wave is absorbed by stars that resonate with the wave. Such stars are concentrated around the radii of the Lindblad resonances. The concentration of the wave’s energy on altering the orbits of a small fraction of the disc’s stars plays a key role in the long-term drift of the disc to bar formation (Sellwood & Carlberg, 2014), so it is a pity that this physics is missed by the current theory. To capture this phenomenon one needs recognise that the perturbation modifies a star’s frequencies in parallel with its orbit. When this fact is omitted, each star’s phase relative to the wave increments steadily, so the star spends a pre-defined time at phases that cause it to absorb energy from the wave, followed by an equal time during which its phase leads to energy being surrendered. After a star becomes resonantly trapped, its phase relative to the wave only librates, and normally is such that energy is absorbed from the wave rather than being emitted.

Equation (108) for the perturbation caused by an orbiting mass implies non-linear dependence on $\Sigma_{0}$ but linear dependence on $M$: we can create a new solution by doubling both $M$ and $\widehat{\Sigma}_{1}$. It follows that any object that moves on a circular orbit will have an effective mass that is significantly greater than its actual mass; this property is not restricted to massive objects such as GMCs. It applies also to stars so long as they move on circular orbits. While few stars do move on perfectly circular orbits, radial oscillations with amplitude significantly smaller than the extent of the wake shown in Fig. 9 should not prevent the formation of a wake. Indeed, Toomre & Kalnajs (1991) showed that in N-body simulations, regions of enhanced density like that shown in Fig. 9 can be detected around individual stars if you stack images of the sheet such that a different star always lies exactly at the centre of the image. Fouvry et al. (2015) showed that the enhancement of the masses of individual stars through wake formation accelerates the relaxation of discs by increasing the level of Poisson noise. Much earlier, Julian (1967) pointed out that GMCs would stochastically accelerate disc stars faster than one would naively expect because their effective masses are larger than their physical masses by virtue of the wakes they raise in the stellar disc.

Giant molecular clouds (GMCs) are thought to be destroyed by outflows from the massive stars that they bring into the world. Hence it is interesting to ask about the persistence of the wake that a GMC generates in the stellar disc. After all, a GMC gathers around it a stellar entourage significantly more massive than itself. Will not this entourage seed an overdensity after the GMC has been dispersed? Somewhat counter-intuitively, the current theory predicts that the wake will quickly disperse rather than survive the GMC’s destruction. The argument is that, as JT66 demonstrated, each spiral wave evolves in isolation, and the wake comprises waves that are already trailing. For the wake to persist, a source is required of leading waves that can be swing amplified. So long as the GMC lives, its gravitational field furnishes that source, but the wake cannot self seed.

The dynamical independence of waves with different ${\bf k}$ values must derive from the orthogonality of $\mathrm{e}^{{\rm i}{\bf k}\cdot{\bf x}}$ and $\mathrm{e}^{{\rm i}{\bf k}^{\prime}\cdot{\bf x}}$ for ${\bf k}\neq{\bf k}^{\prime}$. From a physical perspective this orthogonality is questionable because it hinges on an infinite domain of integration over ${\bf x}$, which is entirely unphysical. Consequently, the implication of the JT equation that GMC’s after images are short-lived should be accepted cautiously.

The central principle of the linear theory used here is integration along unperturbed trajectories. Such an integration will yield useful results so long as the perturbation has not changed the orbit qualitatively. In the case of a mass $M$, some orbits will be qualitatively changed by being trapped by the corotation resonance. For the results computed to be useful, $M$ must be small enough that only a small fraction of stars are trapped by the corotation resonance (e.g. Binney, 2018).

8 Conclusions

Spiral structure is a complex phenomenon and decades have been required to achieve a reasonable understanding of it. That understanding has been assembled by patching together insights from several different approaches. Matrix mechanics (Kalnajs, 1977; Toomre, 1981; Fouvry et al., 2015) and N-body simulations (e.g. Sellwood & Carlberg, 2014; Fouvry et al., 2015) provide the most trustworthy information, but taken alone they yield insufficient insight into the physics of the phenomenon. Consequently, our understanding of the physics of spiral structure is heavily dependent on linear theory (Toomre, 1981). This comes in two variants and requires the use of approximations that are not strictly justifiable.

The LSK theory of running waves is quite well known but it requires the tight-winding approximation, which inevitably fails as leading waves approach corotation. Moreover, taken on its own, the LSK theory is profoundly misleading in that it draws an exclusion zone around corotation and suggests that waves that approach this zone bounce off it. Matrix mechanics and N-body simulations clearly show that far from being a quiet zone, the corotation region lies at the heart of spiral-structure.

The shearing sheet that JT66 introduced to stellar dynamics following the seminal paper of Goldreich & Lynden-Bell (1965) on gas discs, provides the only tractable model of this beating heart of the system, so deserves to be widely understood. It provides a tractable model of a spiral feature that includes the vital time dependence of its wavevector ${\bf k}$. We have re-derived and slightly extended the key results of JT66, and illustrated their value by both reproducing important figures from JT66, Toomre (1969) and Toomre (1981) and by plotting some things that cannot be found in those papers. Our aim has been twofold: (i) to derive the JT equation in a way that is as self-contained and elementary as possible, and then (ii) to explain more fully than JT66 and Toomre (1969) did how this equation can be used to compute a variety of phenomena. These include, the wake that an orbiting body assembles around it, the movement of wavepackets between corotation and Lindblad resonances, and the criterion for a disc’s stability to axisymmetric disturbances. All these applications are quite subtle and merit explicit explanation.

Axisymmetric disturbances cannot wind up, so the wavevector ${\bf k}$ is not a function of time. Consequently, the axisymmetric limit of the JT equation admits modes, and their dispersion relation is just the appropriate LSK relation.

Application of the theory to wavepackets clarifies the standing of the LSK dispersion relation. For sufficiently small azimuthal wavenumbers ($k/k_{\rm crit}\la 0.5$), wavepackets move towards Lindblad resonances while becoming more tightly wound, much as the LSK relation predicts. But for larger wavenumbers, the LSK relation is of no use. The fundamental problem with LSK theory is that in a shearing system, waves will wind up, so ${\bf k}$ should be an explicit function of time, yet such a functional dependence is explicitly excluded in the derivation of a dispersion relation.

Everything here flows from the Volterra integral equation that JT66 first derived, which admits inexpensive numerical solution. Only in the axisymmetric limit $k_{y}\to 0$ is the JT kernel $K(t,t^{\prime})$ translationally invariant, i.e., a function ${\cal K}(t-t^{\prime})$. From this it follows that in the non-axisymmetric case a modal analysis such as that pursued by Lin & Shu (1966) is impossible. The ingredients required to derive the JT equation are (i) the linearised Boltzmann equation, (ii) the approximate solution to Poisson’s equation for short wavelength waves, and (iii) formulae for general unperturbed orbits that are simple enough to permit analytic integration over velocities to obtain the perturbed density from the perturbed DF. Within the epicycle approximation, all of these ingredients appear to be as available for a true disc as the shearing sheet. If it proves impossible to generalise the JT kernel to a disc, the realism of the kernel might still be significantly enhanced by making the epicycle frequency $\kappa$ a decreasing function of $x$ and thus moving corotation towards the outer Lindblad resonance from the midpoint of the gap between the Lindblad resonances.

A striking result obtained here is the clarity with which a swing-amplified wave can be seen in velocity space long after it has disappeared from real space. The Gaia mission (Gaia Collaboration et al., 2018) makes it straightforward to probe in great detail the structure of velocity space at locations up to $\sim 2\,\mathrm{kpc}$ from the Sun. It has been known since the Hipparcos mission that the local $UV$ plane is rich in structure (Dehnen, 1998), but universally accepted explanations of this structure are still lacking. The shearing sheet will surely play a prominent role in efforts to understand velocity-space structures in the new data.

Acknowledgements

I thank Jerry Sellwood and Scott Tremaine for comments on an early draft and Alar Toomre for being a most helpful, meticulous and understanding referee. This work has been supported by the UK Science and Technology Facilities Council under grant number ST/N000919/1.

Appendix A Gaussian integral

We have

[TABLE]

Hence

[TABLE]

Appendix B Example code

The online version of this paper includes a file jt_class.h. Once this header file has been included in a C++ program, the statements

JulianT S(ky,Q); double A=S.Afactor(1000);

will place in A, for a wave with $k_{y}={\tt ky*}k_{\rm crit}$ in a disc with stability parameter Q, the amplification factor analogous to those plotted in Fig. 6. By default the disc has a flat rotation curve, but this can be changed by the statement S.set(Omega,kappa) with the desired frequencies. The current frequencies are extracted by the statement S.fill(Omega,kappa,A)

The statements

int np=1000;

double Sigma[np], theta[np];

S.Afactor(Sigma,theta,np);

will place in the arrays the quantities $\widetilde{\Sigma}_{1}(t)$ and $\kappa t$ shown by the black curves in Fig. 4. The statement

S.Afactor(Sig,Sigma,theta,np);

will additionally place in Sig[np] data for the red curve in Fig. 4. All the above amplification factors are for the case that the disc is jolted at time $t_{\rm i}$ such that $\kappa t_{\rm i}=-3\pi/2$. The statement

Afactor(Sig0,Sigma,theta,theta0,np);

will return the amplification factor, etc, for a disc jolted at $\kappa t_{\rm i}={\tt theta0*}\pi$.

The statement S.K(theta,thetap) will return the value of the JT kernel at the given values of $\kappa t$, while S.Kt(t,tp) will yield the value of the kernel at the given times.

The statement S.Igrand(theta,thetap,Ux,Uy,dpsi) returns the value of the fraction on the second line of equation (85) at the given values of ${\tt theta}=\kappa t$, ${\tt thetap}=\kappa t^{\prime}$ and velocity components, and places in double dpsi the value of $\psi(t^{\prime})-\psi(t)$. Hence the statements

double Skr=sin(kxx+kyy),Ckr=cos(kxx+kyy);

double th=theta[850],dth=theta[1]-theta[0];

double sum=-S.Igrand(th,0,Ux,Uy,dpsi);

sum =(Skrcos(dpsi)+Ckr*sin(dpsi));

for(int it=0;it<=nt;it++){

 double thp=it*dth;

 double I=S.Igrand(th,thp,Ux,Uy,dpsi);

 sum-=Sigma[it]I(Skrcos(dpsi)+Ckrsin(dpsi))*dth;

}

double f1=sumexp(-.5(UxUx+UyUy));

will place in f1 the perturbation to the density at (Ux,Uy) in the velocity space of (x,y) at the time corresponding to theta[850].