The shearing sheet and swing amplification revisited
James Binney



=2pt

    ∂ı iR_ b
        
        
        =cmmib10
        =cmmib7
        ==#1
        	#1#1#1#11112131415161718192021222324252627282930313233ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz.⊤#1#2.#1#1#2.#2. v⊤
        ⊤.ḍv_ c J R I Kk_ BMNRASAnnRA&AApJApJSApJLNatAJA&AA&ASPASPJCAPPASJPhys. Rev. EdfΛCDMLOSIℒS Gyr Myr kpcmixrescutconst km s^-1 km mas yr^-1 mas GeV cm^-3  M_⊙v_∥ pcdeR_cR_LR_d#1#2#1#2HItmr_ sz_ stol_Jrms^ T^ Tθ^ TJ^ TΦ^ T^ T^ Tr^ Tϑ^ Tp^ T tm sgn#1Fig. <ref> sos^2^-1_ pω_ pRΩ
         LSwing amplification revisited]
        The shearing sheet and swing amplification revisitedJames Binney]
          James Binney^1E-mail: binney@thphys.ox.ac.uk
        ^1Rudolf Peierls Centre for Theoretical Physics, Clarendon Laboratory,
          Parks Road, Oxford, OX1 3PU, UK
        [
            [
            December 30, 2023
        =====================
        
        
        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.
        
          Galaxy:
          kinematics and dynamics – galaxies: kinematics and dynamics – methods:
          numerical
        
        
        § 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 <cit.>.
        
        The key ingredients of spiral structure are: the Lin-Shu-Kalnajs (LSK)
        dispersion relation[This relation, first given by <cit.>,
        generalises the relation obtained for axisymmetric disturbances by
        <cit.>.] for running density waves; the swing amplifier
        <cit.>; and resonant absorption at Lindblad
        resonances <cit.>. 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 <cit.>.
        Eventually the tight-winding approximation on which the LSK relation is based
        fails, and following <cit.> 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 <cit.>. However, this will not occur if a Lindblad
        resonance is encountered before the centre is reached because the wave will
        be resonantly absorbed there <cit.> 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.  
        
        
        <cit.>, 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
        <cit.>. 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 <cit.>.
        
        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 <cit.> 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,
        <cit.> or <cit.>, but with limited explanation of the
        computational details. 
        
        Section <ref> presents the shearing sheet from a Hamiltonian
        perspective. Section <ref> sets up the key equation by linearising
        the collisionless Boltzmann equation (CBE). Section <ref> computes
        the signature of a winding spiral in velocity space.
        Section <ref> 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 <ref> uses the key equation to
        investigate the response of the disc to an impulsive jolt.
        Section <ref> uses it to determine the increase in density around an
        orbiting mass that arises because the disc is highly `polarisable', and
        Section <ref> uses it to study the dynamics of wavepackets.
        Section <ref> discusses the insights we gain from the formalism
        and its limitations. Section <ref> sums up. Appendix <ref>
        evaluates a required Gaussian integral, while Appendix <ref> describes
        a C++ class that computes the dynamics of winding waves.
        
        
        
        
        § 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 Ω 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
        +x, and we define a complementary coordinate by y≡ϕ,
        where ϕ is the angle between the lines from the disc's centre to the
        patch's centre and the point (x,y) (Fig. <ref>). In this
        rotating frame, the Lagrangian is
        
            =12[ẋ^2+(+x)^2(+ẏ)^2]-Φ(+x)
        
        because in an inertial frame the angular velocity of the point at  R+x is
        Ω+ϕ̇=Ω+ẏ/(+x)≃Ω+ẏ/. From  we can read off 
        the momenta:
        
            p_x   =ẋ
            p_y   =(+x)^2(+ẏ)1   ≃+2 x+ẏ.
        
        It follows that the Hamiltonian is
        
            H   =12( p_x^2+p_y^2(1+x/)^2)- p_y+Φ.
        
         Since H doesn't contain t or y, we have two constants of motion: H and p_y or 
        
            Δ_y≡ p_y-=2
            x+ẏ.
        
         We will use Δ_y rather than p_y because, unlike p_y, it is
        first-order in the small quantities x,y.
        
        The derivatives of Φ at the origin are
        
            Φ x   =RΩ^2^2Φ x^2   =^2+2   =(-4
            A),
        
         where A is Oort's first constant (Table <ref>). Hence
        
            Φ(+x)≃Φ()+^2 x+12(-4A)x^2,
        
        so
        
            H   ≃12[p_x^2+p_y^2(1-2x+3x^2^2)]-
            p_y+Φ()   3cm+^2x+12(-A)x^2   ≃12[p_x^2+Δ_y^2-^2^2]+Φ()-xΔ_y+12κ^2x^2
        
        where κ^2≡4(-A).
        From the way x and p_x appear in H if follows that x oscillates
        harmonically  at the epicycle frequency κ about
        
            x≡2Δ_y/κ^2.
        
        
         §.§ Orbits
        
        
         Circular orbits are ones on which x=x, so from equation
         (<ref>)
        
            ẏ   =Δ_y-2 x=(κ^22-2)x   =-2Ax
             (circular orbit).
        
        This relation describes how the sheet shears in the absence of random
        velocities.
        
         Let v_ϕ≡2Ax+ ẏ be the azimuthal speed relative to the
        local circular orbit. Then from equation (<ref>) and the definition
        (<ref>) of Δ_y we have
        
            v_ϕ=2(A-)x+κ^2x2=2B(x-x),
        
         where B is Oort's second constant (Table <ref>).
        Hence v_ϕ tells us how far a star is from its guiding centre.
        Given that x oscillates harmonically, the coordinates of any star can be
        written
        
            x   =x+Xcosθ_r  ⇒   p_x=-κ Xsinθ_r
            
            v_ϕ   =2B(x-x)=2BXcosθ_r.
        
        where θ_r=κ t+constant. Similarly,
        
            ẏ   =Δ_y -2 x=Δ_y -2(x+Xcosθ_r)⇒ 
            y(t')   =y_0+Δ_y(1-4Ω^2κ^2)t'-2κXsinθ_r   =y_0+Δ_yA Bt'-2κXsinθ_r,
        
        where y_0 is a constant of integration.
        
         When we eliminate x from H in favour of v_ϕ we get
        
            H   ≃12[p_x^2+Δ_y^2-^2^2]+Φ()
            +12κ^2[(x-x)^2-x^2]   =12[p_x^2+Δ_y^2(1-4^2κ^2)-^2^2]+Φ()
            +κ^2 8B^2v_ϕ^2   =12[p_x^2+Δ_y^2A A--^2^2]+Φ()
            + -2Bv_ϕ^2   =H_x(p_x,v_ϕ)+H_y(Δ_y),
        
        where
        
            H_x≡12(p_x^2+Ω-Bv_ϕ^2)=12κ^2X^2.
        
        Since B<0 and p_x is the radial component of velocity, in velocity space
        lines of constant H_x form ellipses.  From equation (<ref>) we obtain
        an alternative expression for H_x:
        
            H_x=12[p_x^2+κ^2(x-x)^2].
        
        
         §.§ Characteristic scales
        
        
        By Jeans' theorem a suitable equilibrium DF is
        
            f_0=F^-H_x/σ^2,
        
        where F and σ are free constants. By equation (<ref>) this
        DF generates a biaxial Maxwellian velocity distribution with dispersion
        σ in p_x and √(-B/Ω)σ in v_ϕ.  The
        surface density it generates  is
        
            Σ_0=∫p̣_xṿ_ϕ f_0=4πσ^2BFκ.
        
        
        It  is convenient to specify the constants F and σ through the
        values of two functions of them. <cit.> showed that even a stellar disc in which all
        stars are on perfectly circular orbits (so σ=0) is stable to
        axisymmetric perturbations with wavenumbers smaller than  the critical wavenumber 
        
            k_ crit≡2πλ_ crit=κ^22π
            GΣ_0.
        
        We use k_ crit as proxy for the surface density Σ_0 or,
        equivalently, the
        normalisation of the DF, F. Our proxy for σ is
        
            Q≡κσ3.36GΣ_0,
        
        which is the ratio of the actual velocity dispersion σ to the minimum value
        that ensures that the disc is stable to axisymmetric disturbances of any
        wavenumber <cit.>.
        
        If the disc is disturbed by a potential Φ_1∝cos(k_yy), a star
        that is on a circular orbit at x experiences the disturbance at frequency
        
            ω=k_yẏ=-2Axk_y.
        
        At the locations ± x_ L of the Lindblad resonances, the absolute
        value of this frequency coincides with the star's radial frequency κ.
        Clearly the Lindblad resonances of this disturbance lie at
        
            x_ L=±κ2Ak_y.
        
        
        If we adopt plausible values Σ_0=49^-2 and
        κ=37^-1 for the solar neighbourhood <cit.>, then
        λ_ crit≃6 and k_ crit∼1^-1.
        
        
        
        
         §.§ 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. <ref>. Hence the
        pattern shears with particles on circular orbits. That is, if the perturbed
        surface density has the form
        
            Σ_1(,t)=Σ_1(t)^ı·,
         
         the phase · is constant at the location of every particle that
        moves on a circular orbit. For this to be so,  has to be a function of
        time. We determine that function by equating phases at a particle's locations at two
        different times:
        
            k_x(t_ı)x+k_yy_0=k_x(t)x+k_yy(t)⇒     
            k_x(t)=k_x(t_ı)+2k_yA(t-t_ı),
        
         where we have used equation (<ref>).
        Henceforth we take the origin of time to be the instant 
        
            t_c≡-k_x(t_ı) 2k_yA,
         
        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
         
            k_x(t)=2Ak_yt,
        
         and 
        
            |k|=k_y√(1+4A^2t^2)
        
         goes through a minimum as k_x passes
        through zero.
        
        So long as |k|≪ R, the gravitational potential generated by a
        sinusoidal density variation in a razor-thin disc is <cit.>
            Φ_1(,t)=P(t)^ı· =-2π GΣ(,t)|k|.
        
        The strength of this potential peaks as k_x passes through zero.
        
        
        Using our solution (<ref>) and (<ref>) for a general orbit we now have that
        
            k_xx+k_yy   =2k_yAt(x+Xcosθ_r)       +k_y(y_0+Δ_yA
            Bt-2κXsinθ_r)   =k_y[y_0+2X(Atcosθ_r-κsinθ_r)
            ].
        
        If we define the phase 
        
            ψ(t)≡2k_yX(Atcosθ_r-(/κ)sinθ_r).
        
        this simplifies to
        
            ·=k_yy_0+ψ(t).
        
        It follows that if t and t' are two instants along a given
        orbit, then
        
            ·|_t'=·|_t+ψ(t')-ψ(t).
        
        
        § LINEARIZING THE CBE
        
        
        With f=f_0(H_0)+f_1 and H=H_0+Φ_1, the CBE
        
            f t+[f,H]=0,
        
         where [,] is a Poisson bracket, linearises to
        
            f̣_1ṭ≡ f_1 t+[f_1,H_0]=[Φ_1,f_0].
        
         On integration we have
        
            f_1=∫_t_ı^tṭ' Φ_1· f_0,
        
         where the integral is along unperturbed orbits and f_1 has been assumed to
        vanish for t'<t_ı.
        
        The unperturbed  DF (<ref>) is a function of not only p_x,x but
        through x or v_ϕ also
        of the momentum Δ_y. Using equation (<ref>) for H_x
        with
        Φ_1()∝^ı·, we have 
        
            Φ_1· f_0=
            -ı
            FΦ_1^-H_x/σ^2σ^2[k_xp_x-k_y2(x-x)].
         
        
        Since we have assumed that f_1 vanishes at t_ı, it
        will vanish at all subsequent time unless something contributes to Φ_1
        in addition to the surface density associated with f_1. Hence,
        we take
        
            Φ_1(t')=-2π G|k|[Σ_(t')+Σ_1(t')]^ı·|_t'   =-2π G|k|[Σ_(t')+Σ_1(t')]^ı[·|_t+ψ(t')-ψ(t)],
        
         where the second equality uses equation (<ref>),
        Σ_(t')^ı·|_t' is a fictitious surface density that provides a
        gravitational driving
        force, and
        Σ_1 is defined by
        
            Σ_1(,t)=∫^̣2 f_1=Σ_1(t)^ı·.
        
        
        In view of equations (<ref>) to (<ref>), with
        equations (<ref>) and (<ref>)  we have
        
            f_1(t)   =2π GFı^-H_x/σ^2σ^2^ı·∫_t_ı^tṭ' (Σ_+Σ_1)   ×^ı[ψ(t')-ψ(t)]-2At'κ Xsinθ_r-2
            Xcosθ_r√(1+4A^2t^'2).
        
        
        Next we have to compute Σ_1 by integrating  f_1 over ,
        but before integrating we eliminate
        X,θ_r in favour velocity components p_x,v_ϕ at t'=0. Using θ_r=θ_0+κ t' we write
        
            κ Xsinθ_r   =κ X(sinθ_0cosκ
            t'+cosθ_0sinκ t')   =-U_xcosκ t'+U_ysinκ t',
        
        where
        
            U_x   ≡-κ Xsinθ_0=p_x
            U_y   ≡κ
            Xcosθ_0=(κ/2B)v_ϕ.
        
        Similarly
        
            κ Xcosθ_r=U_ycosκ t'+U_xsinκ t'.
        
        The Jacobian of the (time-dependent) transformation in velocity is simply
        
            (p_x,v_ϕ)(U_x,U_y)=2Bκ.
        
        With  replacing (p_x,v_ϕ) the phase factor of the integrand in equation
        (<ref>) becomes
        
            ψ(t')-ψ(t)   =2(k_y/κ)[U_x{A(t'S'-tS)+Ωκ(C'-C)}    +U_y{A(t'C'-tC)-κ(S'-S)}],
        
        where C(t)≡cosκ t, S(t)≡sinκ t and C' and S'
        are defined the same way but with t' replacing t.
        
        Using equation (<ref>) to eliminate F from equation (<ref>), we
        see that the
        coefficient of ^ı· in the expression for the
        perturbed surface density  can be written
        
            Σ_1(t)   =∫p̣_xṿ_ϕ ^-ı· f_1   =∫_t_ı^tκ ṭ' K(t,t')[Σ_(t')+Σ_1(t')],
        
        where
        
            K(t,t')≡I(t,t')√(1+4A^2t^'2)
        
        is a dimensionless kernel that involves the Gaussian integral
        
            I(t   ,t')≡-2ıGΣ_0κσ^4∫^̣2 ^-U^2/2σ^2^ı[ψ(t')-ψ(t)]   ×[At'(-U_xC'+U_yS')+(Ω/κ)(U_yC'+U_xS')].
        
        The integral is of the form
        
            ∫^̣2 · ^-a^2U^2+2ı·=ıπ·
            a^4 ^-b^2/a^2,
        
        where
        
            a^2   =12σ^2
            b_x   =k_yκ[A(t'S'-tS)+(/κ)(C'-C)]
            b_y   =k_yκ[A(t'C'-tC)-(/κ)(S'-S)]
            c_x   =-At' C'+κ S'    
            c_y=At' S'+κ C'.
        
        Hence
        
            I(t,t')=8πGΣ_0κ · ^-b^2/a^2.
        
        In terms of the critical wavenumber (<ref>),
        I simplifies to
        
            I(t,t')=4κ k_ crit · ^-b^2/a^2.
        
        Now we replace  by the dimensionless vector
        
            ≡κ k_ crit
        
        and note that
        
            b  a =√2σ k_ critκb=3.36√22πQb,
        
        where we have introduced Q (eqn. <ref>).
        With equations (<ref>) and (<ref>), 
        the kernel (<ref>) can be written
        
            K(t,t')=4· exp(-0.572Q^2b^2)√(1+4A^2t^'2).
        
        The vectors  and  are functions of the dimensionless numbers At, At',
        κ t, κ t'
        and Ω/κ, while  has an additional (linear) dependence on
        k_y/k_ crit. The argument (Qb)^2 of the exponential is
        proportional to σ^2 and has no dependence on Σ_0.
        
        In the following we shall refer to equation (<ref>) as the JT
        equation and K as the JT kernel. Although JT66 normalised their kernel
        differently, equation (<ref>) is equivalent to their equation.
        
        
        
        
         §.§ Impact on velocity space
        
        
        Once we have computed the evolution of Σ_1 from the JT
        equation, we can use equation (<ref>) to compute  the evolution
        of f_1. Since we can directly observe the phase space of our Galaxy,
        f_1(,,t) contains observationally significant information that was erased when
        we integrated over  to compute Σ_1.
        
        Equation (<ref>) is explicit about the  dependence of
        f_1, but its  dependence is buried in the variables X and
        θ_r. Above we expressed X and θ_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
        θ_r=θ_0+κ(t'-t), where θ_0 is a star's current
        radial phase. With this new definition of θ_0, 
        
            κ Xsinθ_r   =-U_xcos[κ(t'-t)]+U_ysin[(κ(t'-t)]κ Xcosθ_r   = U_ycos[κ(t'-t)]+U_xsin[κ(t'-t)],
        
        where  is related to the current velocity components by equation
        (<ref>). Equation (<ref>) becomes
        
            ψ(t')-ψ(t)   =2(k_y/κ)[U_x{A(t'-t)S'+Ωκ(C'-1)}    +U_y{A(t'-t)C'-κS'}]
        
        where C'=cos[κ(t'-t)] and S'=sin[κ(t'-t)]. Completing the
        transition from X,θ_r to , we find
        
            f_1(t)   =-ıκ^-H_x/σ^23.36QB^ı·∫_t_ı^tκṭ' (Σ_+Σ_1)^ı[ψ(t')-ψ(t)]   ×A(t'-t)(-U_xC'+U_yS')+(Ω/κ)(U_yC'+U_xS')
            σ^3√(1+4A^2t^'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→0
        of the JT equation. Since ∝ k_y and both the prefactor and the
        exponent in equation (<ref>) vanish with k_y, superficial
        analysis yields the conclusion that in this limit K→0, making the JT
        equation useless. But we have to bear in mind that k_x=2Atk_y and we
        require k_x0, so we must let |t|→∞ as we let k_y→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 (<ref>) shows that in this limit k_y→0
        but k_yt=constant,  and  simplify such that
        
            →1 2k_ crit(k_x'S'-k_xS,k_x'C'-k_xC)4√(1+A^2t'^2)   →2(-C',S'),
        
        where k_x'≡ k_x(t'), etc. Now as k_y→0ḳ_x/ṭ→0, so over
        any finite interval t-t' there will be negligible change in k_x.
        Moreover, while we are obliged to let |t,t'|→∞, it suffices to
        consider forces that acted a finite time in the past. That is, we need to
        evaluate K(t,t') for finite |t-t'|, so we can neglect the difference
        between k_x and k_x'. Simplifying the right sides of equations
        (<ref>) thus and using the trigonometrical identities
        CC'+SS'=cos[κ(t-t')] and SC'-CS'=sin[κ(t-t')], we find that
        in the axisymmetric limit k_y→0, the JT kernel reduces to (JT66)
        
            (t-t')   ≡lim_k_y→0K(t,t')   =k_x k_ critsin[κ(t-t')]^(cos[κ(t-t')]-1)χ,
        
        where
        
            χ≡Q^2k_x^23.497 k_
            crit^2=(k_xσκ)^2.
        
        The JT equation now reads
        
            Σ_1(t)=∫_-∞^tκ ṭ' (t-t')[Σ_(t')+Σ_1(t')].
        
        We replace t' by t”≡ t-t'. As t' goes from -∞ to t,
        t” goes from ∞ to 0. Hence the equation becomes
        
            Σ_1(t)=∫_0^∞κ ṭ” (t”)[Σ_(t-t”)+Σ_1(t-t”)].
        
        The right side has the form of a Laplace convolution of  with the sum of
        surface densities, so the Laplace transforms Σ_1(p), etc
        are related by
        
            Σ_1(p)=κ(p)Σ_(p)
            1-κ(p).
        
        Values of p for which κ=1 are of particular interest: for
        these frequencies the JT equation has a non-zero solution in the absence of
        a stimulating density Σ_. That is, these are the frequencies of the
        sheet's axisymmetric normal modes. 
        
        
        Using the identity <cit.>
            ^zcosθ=∑_n=-∞^∞ I_n(z)cos(nθ),
        
        it is easy to show that the Laplace transform of (τ) is
        
            (p)=   k_x^-χ 2k_ crit∑_n=-∞^∞ I_n(χ)((n+1)κ
            p^2+(n+1)^2κ^2           -(n-1)κ p^2+(n-1)^2κ^2).
        
        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:
        
            (p)=k_x^-χ k_ crit∑_n=-∞^∞ I_n(χ)n+1κ p^2+(n+1)^2κ^2.
        
        We relate this to the standard LSK dispersion relation by defining the
        dimensionless frequency s≡ı p/κ. In terms of s the dispersion
        relation κ(p)=1 can be written
        
            0=1-k_x k_ crit^-χ∑_n=-∞^∞ I_n(χ)(n+1) (n+1)^2-s^2.
        
        When the LSK dispersion relation is written in analogous form
        <cit.>, and  divided by 1-s^2, one obtains
        
            0=1-|k| k_ crit^-χ∑_n=-∞^∞I_n(χ)χn^2
            n^2-s^2.
        
         The recurrence relation
        
            I_n-1(χ)-I_n+1(χ)=2nI_n(χ)χ
         
        allows us put the LSK relation into the form
        
            0=1-|k| 2k_ crit^-χ∑_n=-∞^∞ [I_n-1(χ)-I_n+1(χ)] n
            n^2-s^2.
        
        Gathering together the coefficients of I_n makes the sum precisely the sum
        that appears in equation (<ref>), so the same manoeuvre to simplify
        the coefficient of I_n establishes that the LSK dispersion relation can be
        written in the form (<ref>) 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 (<ref>) is plotted as a
        function of k_x/k_ crit at a fixed values of σ and s^2, one
        obtains a curve that has a minimum at k/k_ crit∼2. The smaller the
        value of σ 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_ 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 σ 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.
        
        
        
        
        § IMPULSIVE EXCITATION OF THE DISC
        
        
        For our first application of the JT equation we follow JT66 by exciting the
        disc impulsively. We take 
        
            Σ_(t)=Σ_ sκδ(t-t_ı),
        
         where Σ_ s sets the magnitude of the impulse and the
        denominator κ ensures dimensional soundness.
        With this choice, the JT equation yields
        
            Σ_1(t)=K(t,t_ı)Σ_ s+
            ∫_t_ı^tκṭ' K(t,t')Σ_1(t').
        
        Equation (<ref>) is straightforwardly solved by establishing a grid
        of about a thousand equally spaced times in the interval (t_ı,5π/κ),
        approximating the integral over t' by the trapezium rule on this grid, and
        then for each possible end time t solving the resulting algebraic equation
        for Σ_1(t).
        
        
        
        Fig. <ref> shows solutions obtained in this way for four values
        of t_ı in a disc with a perfectly flat rotation curve, Q=1 and
        k_y=k_ 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_ı≃-3π/2κ.  However, moving t_ı further back decreases the
        peak in Σ(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. <ref> the red dotted curve for κ t_ı=-2π almost
        exactly coincides with the black, long-dashed curve for κ t_ı=-π.
        
        The black curves in Fig. <ref> show, for a fixed value of
        k_y, the evolution of
        Σ_1 for t_ı=-3π/2κ 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. <ref>,
        which plots K(t,t_ı) 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. <ref> show,
        respectively, the responses with and without self-gravity for the case
        Q=1.5 and k=k_ 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 κ t=5π/4 and the second on and 10π/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 κ t=5π/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 κ t≃4π,
        when without self-gravity the response is essentially zero.
        
        The key feature of Fig. <ref> 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') 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. <ref> plots for seven values of Q the amplification
        factor 
        
            A_ max≡max_t_ı{max_t[Σ_1(t)/Σ_ s]max_t[K(t,t_ı)]}
        
         versus the lengthscale λ/λ_ crit. Here the outer max
        operator involves a search over the time t_ı 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_ı that produces the maximum
        response depends strongly on λ/λ_ 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 κ t_ı<π/2. If λ/λ_ crit1,
        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_ı
        is not a continuous function of Q and λ/λ_ crit.
        This fact gives rise to kinks and bumps in the curves of
        Fig. <ref> that might be thought indicatons of numerical error
        but are not.
        
        Fig. <ref> closely resembles the central panel of Fig. 7 in
        <cit.> 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. <ref> 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. <ref> 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
        λ/λ_ crit≃2 before falling to a plateau that extends
        from λ/λ_ crit≃3 to the longest lengthscales.
        Discontinuities in the optimum value of t_ı sometimes generate bumps in
        the plateau. 
        
        In N-body simulations of discs with sustained star formation, Q
        settles to a value around 1.4<cit.>.
        Fig. <ref> indicates that in such a disc perturbations are
        swing-amplified by a factor that hovers around 14 for
        1.3<λ/λ_ crit<2.1.
        
        Fig. <ref> 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
        (<ref>) at t=-3π/2κ is
        swing amplified. Successive panels show velocity space at intervals
        π/κ apart, starting with t=0, when k_x=0. The dotted curve in
        Fig. <ref> shows the evolution of this wave's density in real space.
        This peaks before t=3π/2κ and has long vanished by t=4π/κ. The
        central panel in the bottom row of Fig. <ref> 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  of
        the plotted velocity space, becoming less symmetric as one moves away from
        corotation. In general it looks like a network of cells.
        
        
        
        § RESPONSE TO A CLOUD
        
        
        We now compute the response of a stellar disc to the insertion at time t_ı
        of a mass such as a molecular cloud that moves on a circular orbit.  The mass
        provides an endless succession of the δ-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
        
            Σ_ e()=M ^-||^2/2Δ^22πΔ^2
        
        we find that every wavenumber has non-vanishing amplitude:
        
            Σ_ e()≡∫^̣2 ^-ı· Σ_
            e() =M ^-||^2Δ^2/2.
        
        Notice that Σ has dimensions of mass, unlike the quantity
        Σ_1(t) defined above as the coefficient of
        ^ı· in the expression for Σ_1(,t). This difference
        reflects the fact that to recover Σ() from
        Σ()^ı· we have to integrate over .
        Fortunately, for given , Σ_1 will evolve in time in just
        the same way that Σ_1 does, so in the JT equation
        (<ref>) we can replace Σ_1 by Σ_1.
        Then we have
        
            Σ_1(,t)
            =∫_t_ı^tṭ' K(t,t')[M^-||^2Δ^2/2+Σ_1(,t')].
        
        
        Fig. <ref> shows solutions to equation (<ref>) with
        Δ=0 for
        k_y=k_ crit/2 and several values of the time t_ı at which the
        mass M is added. The corresponding values of κ t_ı/π 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_ı=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 κ t_ı/π∼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. <ref>
        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. <ref> 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_yt
        (eqn. <ref>). 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. <ref> 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. <ref> we can infer the value of Σ_1() along a
        line of constant k_y in the (k_x,k_y) plane, and to determine
        Σ_1() over the the half plane k_y>0 we just need to solve
        equation (<ref>) for  positive values of k_y. Since
        Σ_1() is the Fourier transform of a real function, it
        satisfies
        
            Σ_1(-)=Σ_1^*(),
        
        and this relation gives the values we need in the other half plane.
        
        Fig. <ref> shows contours of constant overdensity
        
            Σ_1()=∫^̣2(2π)^2 ^ı· Σ_1()
        
        for a cloud with characteristic size Δ=0.05λ_ crit in a
        disc with a flat rotation curve and Q=1.4. The contour values are in units
        of M/(0.1λ_ crit)^2, so if the mass were spread at uniform
        density it would contribute one unit within a square 0.1λ_ 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.
        
        
        
        
        § WAVE PACKETS
        
        
        We now excite the disc with the potential of the mass
        distribution
        
            Σ_(,t)=M√(2π)Δexp((x-x_ı)^22Δ^2)
            cos(k_ı y)^-t^2/t_0^2.
        
        The Fourier transform of the spatial part of this distribution is
        
            Σ_()=πM^-k_x^2Δ^2/2-ı
            k_xx_ı[δ(k_y+k_ı)+δ(k_y-k_ı)].
        
        When this transform is used in the JT equation (<ref>) with
        Σ replaced by Σ as in the last section, both
        Σ_ and Σ_1 become explicit functions of
        . Hitherto  the initial value of k_x has been 
        encoded via equation (<ref>) in the variable t_ı.   Now we
        treat t_ı as a function of k_x(0):
        
            t_ı=k_x(0)2Ak_y.
        
        As a consequence of this dependence of t_ı on , the upper limit of
        the time integral in the JT equation loses its status as the current time. Instead we write
        
            t=t_ı+τ,
        
        where τ is the current time. When evaluating
        Σ_(,τ') we use the current value of k_x:
        
            k_x(τ')=k_x(0)+2Aτ' k_y.
        
        When we inverse Fourier transform Σ_1(,τ) to recover
        Σ_1(,τ) we likewise set k_x to its value at time τ.
        
        
        
        In the previous applications of the JT equation, Σ_ has been
        real, with the consequence that Σ_1 has remained real. Now
        Σ_ is complex, so Σ_1=+ı also becomes complex.
        Using equation (<ref>) we reason that
        
            Σ_1(x,t)   =∫_0^∞ḳ_y2π∫_-∞^∞ḳ_x2π(^ı·Σ_1(,t)
            +^-ı·Σ_1(-,t)
            )   =2∫_0^∞ḳ_y2π∫_-∞^∞ḳ_x2π[cos(·)-sin(·))]   =∫_-∞^∞ḳ_x2π{[cos(k_xx)-sin(k_xx)]cos(k_ı y)       -[sin(k_xx)+cos(k_xx)]sin(k_ı y)},
        
        where in the last line +ı stands for the coefficient of
        πδ(k_y-k_ı) in the transform Σ_1()
        (cf. eqn. <ref>).
        
        Fig. <ref> shows the result of exciting the disc with the external
        density (<ref>) when Q=1.2, Δ=x_ L/√2 and
        k_ı/k_ crit=κ/(4π A)=0.225, where x_ L is the distance
        of the wave's Lindblad resonance from the corotation resonance
        (eqn. <ref>). Although the exciting density is centred on
        x=-1.5x_ 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.
        <cit.> showed that the trajectories of the packets' centres can be
        accurately predicted from the group velocity implied by the LSK dispersion
        relation.
        
        If Fig. <ref> a subsidiary packet develops  inside the radius
        x_ı (marked in red) on which the exciting potential is centred. The left
        edge of this packet lies near the value x=-2x_ L at which stars on
        circular orbits
        perceive the wave at a frequency ω=2κ. 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_ crit than applies in Fig. <ref>
        yield similar wavepackets but ones that decay more slowly. The example shown
        in Fig. <ref> has k_y/k_ 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. <ref> shows how this
        disturbance looks in the xy plane at six of the times plotted in
        Fig. <ref>. The growth in the wave's amplitude up to
        t∼3P_0=6π/κ 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. <ref> shows the effect of increasing k_y to k_
        crit/2. At this shorter azimuthal wavelength, the disc is significantly
        more responsive (Fig. <ref>) and the disturbance stimulated by
        the same mass is ∼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. <ref> or to t=10P_0
        in Fig. <ref>. 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.
        
         
        
        
         §.§ Application of the LSK dispersion relation
        
        
        In Section <ref> we showed that the axisymmetric limit of the JT
        kernel gives rise to a dispersion relation (<ref>) that is identical to the
        axisymmetric limit of the LSK dispersion relation.
        We now use this  relation to elucidate the wavepackets
        evident in Figs. <ref> and <ref>.
        
        
        
        Fig. <ref> plots the right-hand side of equation (<ref>)
        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 (<ref>), 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→-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 
        
            s=mΩ_ p-Ωκ
        
        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. <ref> and
        <ref>, 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. <ref> and
        <ref> the wavepackets  somewhat exceed these bounds, but qualitatively the
        agreement is good. 
        
        Hence we can interpret the wavepackets seen in
        Figs. <ref> and <ref> as follows. The disturbance
        created by the external driving potential Φ_ 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 Φ_. At x<0 these
        locations are bounded on the left by x=-2x_ L and x=-x_ L,
        while at x>1 they are bounded on the right by x=x_ L. Since
        Φ_ 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=± x_ L, x=±2x_ 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.
        
        
        
        
        
        
        
        
        § DISCUSSION
        
        
        A priori it is not obvious that in a warm, self-gravitating disc a
        perturbation with the initial form Σ_1()∝^ı·
        will evolve such that Σ() remains proportional to
        ^ı· with  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  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 <cit.> showed move towards the Lindblad resonances.
        Consequently, a packet's central  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 <ref> 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 <cit.> 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 (0.5k_ 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_ crit
        that <cit.> overlooked, is a wavepacket that fills the region at
        -2<x/x_ L<-1, which lies inside the Lindblad resonance at x=-x_
        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_
        L<0 that <cit.> 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_ 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 κ 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
        <cit.> and calculations based on matrix mechanics
        <cit.> 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 <cit.>, 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 (<ref>) for the perturbation caused by an orbiting mass
        implies non-linear dependence on Σ_0 but linear dependence on M: we
        can create a new solution by doubling both M and Σ_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. <ref> should not
        prevent the formation of a wake. Indeed, <cit.> showed that in
        N-body simulations, regions of enhanced density like that shown in
        Fig. <ref> 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. <cit.> 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,
        <cit.> 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  values must derive
        from the orthogonality of ^ı· and ^ı'· for
        '. From a physical perspective this orthogonality is questionable
        because it hinges on an infinite domain of integration over , 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 <cit.>.
        
        
        
        § 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 <cit.> and N-body
        simulations <cit.> 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 <cit.>.
        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 <cit.> 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 .  We have re-derived and
        slightly extended the key results of JT66, and illustrated their value by both
        reproducing important figures from JT66, <cit.>  and
        <cit.> 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
        <cit.> 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  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_ crit0.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
         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→0 is the JT kernel K(t,t') translationally
        invariant, i.e., a function (t-t'). From this it follows that in the
        non-axisymmetric case a modal analysis such as that pursued by
        <cit.> 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 κ 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 <cit.> makes it straightforward to
        probe in great detail the structure of velocity space at locations up to
        ∼2 from the Sun. It has been known since the Hipparcos mission that
        the local UV plane is rich in structure <cit.>, 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.
        
        
        mn2e
        
        § GAUSSIAN INTEGRAL
        
        We have
        
            ∫ụ (u+c)^-a^2u^2+2ı bu   =^-b^2/a^2∫ụ (u+c)^-(au-ı
            b/a)^2   =√π^-b^2/a^2 a^2(ca+ı b/a).
        
        Hence
        
            ∫ụ_y    ^-a^2u_y^2+2ı b_yu_y∫ụ_x ^-a^2u_x^2+2ı
            b_xu_x(c_xu_x+c_yu_y)   =
            √π^-b_x^2/a^2
            a^2∫ụ_y ^-a^2u_y^2+2ı b_yu_y(c_yu_ya+ı c_xb_x/a)   =ıπ
            a^4^-(b_x^2+b_y^2)/a^2( c_xb_x+ c_yb_y).
        
        
        § 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
        
        =10pt 
        JulianT S(ky,Q); double A=S.Afactor(1000);
         will place in  A, for a wave with k_y= ky*k_ crit
        in a disc with stability parameter  Q, the amplification factor analogous to
        those plotted in Fig. <ref>. 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
        
        =10pt 
        int np=1000;
        double Sigma[np], theta[np];
        S.Afactor(Sigma,theta,np);
         will place in the arrays the quantities Σ_1(t)
        and κ t shown by the black curves in Fig. <ref>.
        The statement
        
        =10pt 
        S.Afactor(Sig,Sigma,theta,np);
         will  additionally place in  Sig[np] data
        for the red curve in Fig. <ref>. All the above amplification
        factors are for the case that the disc is jolted at time t_ı such that
        κ t_ı=-3π/2. The statement
        
        =10pt 
        Afactor(Sig0,Sigma,theta,theta0,np);
         will return the amplification factor, etc, for a disc jolted at
        κ t_ı= theta0*π.
        
        The statement  S.K(theta,thetap) will return the value of the JT kernel
        at the given values of κ 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 (<ref>) at the given values of
        theta=κ t, thetap=κ t' and velocity components, and places
        in  double dpsi the value of ψ(t')-ψ(t). Hence the statements
         
        =10pt 
        double Skr=sin(kx*x+ky*y),Ckr=cos(kx*x+ky*y);
        double th=theta[850],dth=theta[1]-theta[0];
        				double sum=-S.Igrand(th,0,Ux,Uy,dpsi);
        				sum *=(Skr*cos(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*(Skr*cos(dpsi)+Ckr*sin(dpsi))*dth;
        				}
        				double f1=sum*exp(-.5*(Ux*Ux+Uy*Uy));
         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].
        
         =0pt
        #include <stdio.h>
        #include <math.h>
        #define MAX(A,B) ((A)>(B)?(A):(B))
        
        class JulianT{  	private:
            		double ky,Q;
            		double Omega,kappa,A,B,sigma,sig2,Pi;
            		double K(double theta,double thetap){      			double S=sin(theta),C=cos(theta),
                			Sp=sin(thetap),Cp=cos(thetap);
              			double bx=ky/(kappa*kappa)
                			*(A*(thetap*Sp-theta*S)+Omega*(Cp-C));
              			double by=ky/(kappa*kappa)
                			*(A*(thetap*Cp-theta*C)-Omega*(Sp-S));
              			double cx=(-A*thetap*Cp+Omega*Sp)/kappa;
              			double cy=( A*thetap*Sp+Omega*Cp)/kappa;
              			return -Pi*8*sigma*exp(-(bx*bx+by*by)*2*sig2)
                			*(cx*bx+cy*by)/sqrt(1+pow(2*A/kappa*thetap,2));
            }  	public:
            		JulianT(double kyi,double Qi){      			Q=Qi; Omega=1; kappa=sqrt(2)*Omega;
              			A=Omega*(1-.25*pow(kappa/Omega,2)); 
              			B=A-Omega;
              			sigma=1; sig2=sigma*sigma; Pi=acos(-1);
              			double kcrit=3.36*kappa*Q/(2*Pi*sigma);
              			ky=kyi*kcrit; 
            }    		double Afactor(){      			int np=500;
              			double Sigma[np],theta[np];
              			Sigma[0]=0; theta[0]=-1.5*Pi;
              			double dtheta=6.5*Pi/(double)(np-1);
              			double fac=1/(3.36*Q),Sigmax=0,Kmax=0;
              			for(int i=1;i<np;i++){        				theta[i]=theta[i-1]+dtheta;
                				double sum=0,Sig0=-fac*K(theta[i],theta[0]);
                				for(int j=1;j<i;j++)
                  					sum+=K(theta[i],theta[j])*Sigma[j];
                				Sigma[i]=(Sig0-fac*dtheta*sum)
                  			/(1-.5*dtheta*fac*K(theta[i],theta[i]));
                				Kmax=MAX(Kmax,fabs(Sig0));
                				Sigmax=MAX(Sigmax,fabs(Sigma[i]));
              }      			return Sigmax/Kmax;
            }};
        
        
        When we consider  non-axisymmetric LSK waves (ones with m0), we associate each value of
        s=m(Ω_ p-Ω)/κ with a given radius. Corotation (x=0)
        is associated with s=0, and the Lindblad resonances are associated with
        s=±1.
