SELF-GRAVITY, RESONANCES AND ORBITAL DIFFUSION IN STELLAR DISKS

Since we are imagining that stars are conserved, the equation that governs the secular evolution of the DF takes the form

$$\begin{eqnarray}&&{\displaystyle \frac{\partial f}{\partial t}}=-{\displaystyle \frac{\partial }{\partial \boldsymbol{J}}}\cdot \boldsymbol{F},\end{eqnarray} \tag{ 1 }$$

where $\boldsymbol{F}$ is the diffusive flux of stars in action space. If $\dot{P}(\boldsymbol{J},\bf{\Delta })$ is the rate of increase with time of the probability that a star scatters from $\boldsymbol{J}$ to $\boldsymbol{J}+\bf{\Delta }$, then $\boldsymbol{F}$ is given by (Binney & Lacey 1988)

$$\begin{eqnarray}&&{F}_{i}=f\bar{{\rm{\Delta }}_{i}}-{\displaystyle \frac{1}{2}}{\displaystyle \frac{\partial f\bar{{\rm{\Delta }}_{{ij}}^{2}}}{\partial {J}_{j}}},\end{eqnarray} \tag{ 2 }$$

where the first- and second-order diffusion coefficients are

$$\begin{eqnarray}\begin{array}{ccc}\bar{{\rm{\Delta }}_{i}}(\boldsymbol{J}) & = & {\displaystyle \int }{d}^{3}\bf{\Delta }\;{\rm{\Delta }}_{i}\dot{P}(\boldsymbol{J},\bf{\Delta })\\ \bar{{\rm{\Delta }}_{{ij}}^{2}}(\boldsymbol{J}) & = & {\displaystyle \int }{d}^{3}\bf{\Delta }\;{\rm{\Delta }}_{i}{\rm{\Delta }}_{j}\dot{P}(\boldsymbol{J},\bf{\Delta }).\end{array}\end{eqnarray} \tag{ 3 }$$

Binney & Lacey (1988) showed that in the relevant circumstances the first- and second-order diffusion coefficients are related by

$$\begin{eqnarray}&&\bar{{\rm{\Delta }}_{i}}={\displaystyle \frac{1}{2}}{\displaystyle \frac{\partial \bar{{\rm{\Delta }}_{{ij}}^{2}}}{\partial {J}_{j}}},\end{eqnarray} \tag{ 4 }$$

so the diffusive flux can be written entirely in terms of the second-order coefficients:

$$\begin{eqnarray}&&{F}_{i}=-{\displaystyle \frac{1}{2}}\bar{{\rm{\Delta }}_{{ij}}^{2}}{\displaystyle \frac{\partial f}{\partial {J}_{j}}}.\end{eqnarray} \tag{ 5 }$$

By expanding $\psi (\boldsymbol{x},t)$, the fluctuating part of the gravitational potential, in angle-action variables,

$$\begin{eqnarray}&&\psi (\boldsymbol{x},t)=\psi (\boldsymbol{\theta },\boldsymbol{J},t)={\displaystyle \sum }_{\boldsymbol{m}}{\psi }_{\boldsymbol{m}}(\boldsymbol{J},t){e}^{i\boldsymbol{m}\cdot \boldsymbol{\theta }},\end{eqnarray} \tag{ 6 }$$

Binney & Lacey (1988) showed that the second-order diffusion coefficients are related to the fluctuations in the potential by

$$\begin{eqnarray}&&\bar{{\rm{\Delta }}_{{ij}}^{2}}(\boldsymbol{J})={\displaystyle \sum }_{\boldsymbol{m}}{m}_{i}{m}_{j}{\tilde{c}}_{\boldsymbol{m}}(\boldsymbol{J},\boldsymbol{m}\cdot \bf{\Omega }(\boldsymbol{J})),\end{eqnarray} \tag{ 7 }$$

where an overline indicates an ensemble average and ${\tilde{c}}_{\boldsymbol{m}}$ is the Fourier transform with respect to time of the autocorrelation of the $\boldsymbol{m}$ Fourier component of the potential

$$\begin{eqnarray}&&{\tilde{c}}_{\boldsymbol{m}}(\boldsymbol{J},\nu )={{\displaystyle \int }}_{-\infty }^{\infty }d\tau \;{e}^{i\nu \tau }\;\bar{{\psi }_{\boldsymbol{m}}(\boldsymbol{J},t){\psi }_{\boldsymbol{m}}^{*}(\boldsymbol{J},t-\tau )}.\end{eqnarray} \tag{ 8 }$$

By the Wiener–Khinchin theorem, ${\tilde{c}}_{\boldsymbol{m}}(\boldsymbol{J},\nu )$ is the power spectrum of the stationary random variable ${\psi }_{\boldsymbol{m}}(\boldsymbol{J},t)$. Equation (7) tells us that diffusion is driven by resonances because it implies that the rate at which stars diffuse from action $\boldsymbol{J}$ is proportional to the power that the fluctuating field has at any of the orbit's characteristic frequencies $\boldsymbol{m}\cdot \bf{\Omega }(\boldsymbol{J})$. Hence, if the fluctuations are confined to a narrow frequency range, perhaps because they are associated with spiral arms, stars that respond strongly to them will be located at only a few points in action space.

We will find it expedient to expand the fluctuating potential $\psi (\boldsymbol{x},t)$ in a set of basis functions, the members of which are enumerated by an index q:

$$\begin{eqnarray}&&\psi (\boldsymbol{x},t)={\displaystyle \sum }_{q}{b}_{q}(t){\psi }^{(q)}(\boldsymbol{x}),\end{eqnarray} \tag{ 9 }$$

where ${b}_{q}(t)$ is a random variable. Following Kalnajs (1971), we require our basis potentials to be orthonormal to the densities ${\rho }^{(p)}(\boldsymbol{x})$ that generate them, so we have

$$\begin{eqnarray}&&{\rm{\nabla }}^{2}{\psi }^{(p)}=4\pi G{\rho }^{(q)}\ \mathrm{and}\int {d}^{3}\boldsymbol{x}\;{\rho }^{(p)}(\boldsymbol{x}){[{\psi }^{(q)}(\boldsymbol{x})]}^{*}=-{\delta }_{{pq}}.\end{eqnarray} \tag{ 10 }$$

Now we have

$$\begin{eqnarray}\begin{array}{ccc}{\psi }_{\boldsymbol{m}}(\boldsymbol{J},t) & = & {\displaystyle \frac{1}{{(2\pi )}^{3}}}{\displaystyle \int }{d}^{3}\boldsymbol{\theta }\;{e}^{-i\boldsymbol{m}\cdot \boldsymbol{\theta }}\psi (\boldsymbol{\theta },\boldsymbol{J},t)\\ & = & {\displaystyle \sum }_{p}{b}_{p}(t){\psi }_{\boldsymbol{m}}^{(p)}(\boldsymbol{J}),\end{array}\end{eqnarray} \tag{ 11 }$$

where

$$\begin{eqnarray}&&{\psi }_{\boldsymbol{m}}^{(p)}(\boldsymbol{J})\equiv {\displaystyle \frac{1}{{(2\pi )}^{3}}}\int {d}^{3}\boldsymbol{\theta }\;{e}^{-i\boldsymbol{m}\cdot \boldsymbol{\theta }}{\psi }^{(p)}[\boldsymbol{x}(\boldsymbol{\theta },\boldsymbol{J})].\end{eqnarray} \tag{ 12 }$$

Hence the required power spectrum is

$$\begin{eqnarray}&&{\tilde{c}}_{\boldsymbol{m}}(\boldsymbol{J},\nu )={\displaystyle \sum }_{{pq}}{B}_{{pq}}(\nu ){\psi }_{\boldsymbol{m}}^{(p)}(\boldsymbol{J}){{\psi }_{\boldsymbol{m}}^{(q)}}^{*}(\boldsymbol{J}),\end{eqnarray} \tag{ 13 }$$

where

$$\begin{eqnarray}&&{B}_{{pq}}(\nu )\equiv {{\displaystyle \int }}_{-\infty }^{\infty }d\tau \;{e}^{i\nu \tau }\;\bar{{b}_{p}(t){b}_{q}^{*}(t-\tau )}\end{eqnarray} \tag{ 14 }$$

is the Fourier transform of the cross-correlation of the amplitudes of the p and q basis functions.

Below we shall require an expression for ${B}_{{pq}}(\nu )$ in terms of the Fourier transform

$$\begin{eqnarray}&&{\tilde{b}}_{p}(\nu )=\int {dt}\;{e}^{i\nu t}{b}_{p}(t)\end{eqnarray} \tag{ 15 }$$

of ${b}_{p}(t)$. If $\psi (t)$ is a stationary random process, then it is straightforward to show that

$$\begin{eqnarray}&&\bar{{\tilde{b}}_{p}(\nu ){\tilde{b}}_{q}^{*}({\nu }^{\prime })}=2\pi \delta (\nu -{\nu }^{\prime }){B}_{{pq}}(\nu ).\end{eqnarray} \tag{ 16 }$$

As we indicated in the Introduction, a stellar disk is exposed to several sources of fluctuations. The issue that we now have to confront is that the fluctuation ψ in the potential that stars experience, which is what appears in the above formulae, differs from the original stimulation, ${\psi }^{e}$, because the disk has non-negligible mass, so through Poisson's equation it makes a contribution ${\psi }^{s}$ to the actual gravitational potential ψ (Weinberg 2001a). We shall refer to ${\psi }^{e}$ as the "bare" stimulus and to

$$\begin{eqnarray}&&\psi (t)={\psi }^{e}(t)+{\psi }^{s}(t)\end{eqnarray} \tag{ 17 }$$

as the "dressed" stimulus. We now seek a relationship between the dressed and bare stimuli.

Let $\psi \prime$ be the change in the potential that the disk would generate if its particles moved in the sum of the unperturbed potential and the stimulating potential ${\psi }^{e}$. Then ${\psi }^{\prime }(\boldsymbol{x})$ is linearly related to ${\psi }^{e}(\boldsymbol{x})$ so for each time-lapse τ there is a linear response operator $M(\tau )$ that connects these functions

$$\begin{eqnarray}&&{\psi }^{\prime }(t)={{\displaystyle \int }}_{-\infty }^{t}{dt}\prime \;M(t-t\prime ){\psi }^{e}(t\prime ).\end{eqnarray} \tag{ 18 }$$

Since the mass of the disk actually contributes to the potential in which its particles move, changes in the disk's potential at a early time t' contribute alongside ${\psi }^{e}(t\prime )$ to the disturbance of disk particles at later times, so the fluctuating component of the disk's potential ${\psi }^{s}$ satisfies

$$\begin{eqnarray}&&{\psi }^{s}(t)={{\displaystyle \int }}_{-\infty }^{t}{dt}\prime \;M(t-t\prime )[{\psi }^{s}(t\prime )+{\psi }^{e}(t\prime )].\end{eqnarray} \tag{ 19 }$$

Inserting this expression into Equation (17), we obtain

$$\begin{eqnarray}&&\psi (t)={\psi }^{e}(t)+{{\displaystyle \int }}_{-\infty }^{t}{dt}\prime \;M(t-t\prime )\psi (t\prime ).\end{eqnarray} \tag{ 20 }$$

In this equation the potentials are functions of $\boldsymbol{x}$ as well as t and $M(t-t\prime )$ is an operator that maps one function of space onto another. The basis functions introduced above reduce this operator to a matrix, so when we write

$$\begin{eqnarray}\begin{array}{ccc}{\psi }^{e}(\boldsymbol{x},t) & = & {\displaystyle \sum }_{p}{a}_{p}(t){\psi }^{(p)}(\boldsymbol{x})\\ \psi (\boldsymbol{x},t) & = & {\displaystyle \sum }_{p}{b}_{p}(t){\psi }^{(p)}(\boldsymbol{x}),\end{array}\end{eqnarray} \tag{ 21 }$$

Equation (20) can be written

$$\begin{eqnarray}\begin{array}{ccc}{\displaystyle \sum }_{p}{b}_{p}(t){\psi }^{(p)}(\boldsymbol{x}) & = & {\displaystyle \sum }_{q}\left[{a}_{q}(t){\psi }^{(q)}(\boldsymbol{x})\right.\\ & & +\left.{{\displaystyle \int }}_{-\infty }^{t}{dt}\prime M(t-t\prime ){b}_{q}(t\prime ){\psi }^{(q)}(\boldsymbol{x})\right].\end{array}\end{eqnarray} \tag{ 22 }$$

We multiply both sides of this equation by $-\int {d}^{3}\boldsymbol{x}\;{[{\rho }^{(r)}(\boldsymbol{x})]}^{*}$ and with Equation (10) obtain

$$\begin{eqnarray}&&{b}_{r}(t)={a}_{r}(t)+{\displaystyle \sum }_{q}{{\displaystyle \int }}_{-\infty }^{t}{dt}\prime {M}_{{rq}}(t-t\prime ){b}_{q}(t\prime ),\end{eqnarray} \tag{ 23 }$$

where

$$\begin{eqnarray}&&{M}_{{rq}}(t-t\prime )=-\int {d}^{3}\boldsymbol{x}\;{[{\rho }^{(r)}(\boldsymbol{x})]}^{*}\;M(t-t\prime )\;{\psi }^{(q)}(\boldsymbol{x}).\end{eqnarray} \tag{ 24 }$$

The temporal convolution can be reduced to a multiplication by taking a Fourier transform: multiplication of Equation (23) by $\int {dt}\;{e}^{i\nu t}$ yields

$$\begin{eqnarray}&&{\tilde{b}}_{r}(\nu )={\tilde{a}}_{r}(\nu )+{\displaystyle \sum }_{q}{\tilde{M}}_{{rq}}(\nu ){\tilde{b}}_{q}(\nu ).\end{eqnarray} \tag{ 25 }$$

Hence

$$\begin{eqnarray}&&\tilde{\boldsymbol{b}}(\nu )={[\boldsymbol{I}-\tilde{\boldsymbol{M}}(\nu )]}^{-1}\tilde{\boldsymbol{a}}(\nu ),\end{eqnarray} \tag{ 26 }$$

where boldface implies vectors and matrices indexed with p.

Equations (16) and (26) enable us to relate B_pq to the basis coefficients of the stimulating field

$$\begin{eqnarray}&&\begin{array}{c}\boldsymbol{B}(\nu )={\displaystyle \frac{1}{2\pi }}{\displaystyle \int }d\nu \prime \;{[\boldsymbol{I}-\tilde{\boldsymbol{M}}(\nu )]}^{-1}\bar{\tilde{\boldsymbol{a}}(\nu )\otimes {\tilde{\boldsymbol{a}}}^{*}({\nu }^{\prime })}\\ \times {[\boldsymbol{I}-{\tilde{\boldsymbol{M}}}^{\dagger }({\nu }^{\prime })]}^{-1}.\end{array}\end{eqnarray} \tag{ 27 }$$

We will show below that analogously to Equation (16)

$$\begin{eqnarray}&&\bar{{\tilde{a}}_{p}(\nu ){\tilde{a}}_{q}^{*}({\nu }^{\prime })}=2\pi \delta (\nu -{\nu }^{\prime }){A}_{{pq}}(\nu ).\end{eqnarray} \tag{ 28 }$$

Hence Equation (13) can be written

$$\begin{eqnarray}\begin{array}{ccc}{\tilde{c}}_{\boldsymbol{m}}(\boldsymbol{J},\nu ) & = & {\displaystyle \sum }_{{pq}}{\psi }_{\boldsymbol{m}}^{(p)}(\boldsymbol{J}){{\psi }_{\boldsymbol{m}}^{(q)}}^{*}(\boldsymbol{J})\\ & & \times {\left\{{[\boldsymbol{I}-\tilde{\boldsymbol{M}}(\nu )]}^{-1}\boldsymbol{A}(\nu ){[\boldsymbol{I}-{\tilde{\boldsymbol{M}}}^{\dagger }(\nu )]}^{-1}\right\}}_{{pq}}.\end{array}\end{eqnarray} \tag{ 29 }$$

Our derivation of the dressed secular diffusion coefficients sketched previously is based on the master equation (Binney & Tremaine 2008) which led to the first- and second-order diffusion coefficients from Equation (3). One can also recover these diffusion coefficients via a timescale decoupling of the collisionless Boltzmann equation (Weinberg 2001a; Pichon & Aubert 2006; Chavanis 2012; Fouvry et al. 2015). Various sources of external perturbations can then be considered to induce secular evolution (Weinberg 2001b; Aubert & Pichon 2007).

Since we are working in two dimensions, the basis functions ${\psi }^{(p)}(\boldsymbol{x})$ become functions ${\psi }^{(p)}(R,\phi )$ of plane polar coordinates that are orthonormal to the generating surface densities ${\rm{\Sigma }}^{(p)}(R,\phi )$. Our problem is simplified if we can choose the basis ${\psi }^{(p)}$ such that the response operator $\tilde{M}$ is diagonal and we now show that this is possible.

It is well known that the potential generated via Poisson's equation by a tightly wound spiral wave is itself a spiral wave with the same wavevector $\boldsymbol{k}=({k}_{r},{k}_{\phi })$ (e.g., Binney & Tremaine 2008, Section 6.2.2). Moreover, the computation of the dynamical response of particles within our disk to a tightly wound perturbing potential that oscillates at temporal frequency ν is covered by standard texts (see, for example, Binney & Tremaine 2008; Section 6.2.2(d) or Binney 2013; Section 4.2 for two different approaches). The result is that a spiral potential ${\psi }^{(p)}(\boldsymbol{x}){e}^{i\nu t}$ creates a spiral perturbation in the surface density of test particles, which through Poisson's equation creates a response potential ${\psi }^{\prime }(\boldsymbol{x}){e}^{i\nu t}$ that differs from the original stimulating potential only in magnitude. In fact

$$\begin{eqnarray}&&{\psi }^{\prime }(R,\phi )={\lambda }_{\boldsymbol{k}}{\psi }^{(p)}(R,\phi ),\end{eqnarray} \tag{ 30 }$$

where

$$\begin{eqnarray}&&{\lambda }_{\boldsymbol{k}}={\displaystyle \frac{2\pi G\rm{\Sigma }| {k}_{r}| }{{\kappa }^{2}(1-{s}^{2})}}\mathcal{F}(s,\chi ).\end{eqnarray} \tag{ 31 }$$

Here Σ is the disk's surface density,

$$\begin{eqnarray}&&s\equiv {\displaystyle \frac{\nu -{k}_{\phi }\;{\rm{\Omega }}_{\phi }}{\kappa }}\end{eqnarray} \tag{ 32 }$$

is the ratio of the frequency at which a star experiences the perturbation to the epicycle frequency, and $\mathcal{F}$ is the reduction factor (Kalnajs 1965; Lin & Shu 1966)

$$\begin{eqnarray}&&\mathcal{F}(s,\chi )\equiv 2\;(1-{s}^{2}){\displaystyle \frac{{e}^{-\chi }}{\chi }}{\displaystyle \sum }_{{m}_{r}=1}^{+\infty }{\displaystyle \frac{{\mathcal{I}}_{{m}_{r}}(\chi )}{1-{s}^{2}{/m}_{r}^{2}}},\end{eqnarray} \tag{ 33 }$$

where ${\mathcal{I}}_{{m}_{r}}$ is a modified Bessel function and the dimensionless quantity

$$\begin{eqnarray}&&\chi \equiv {\displaystyle \frac{{\sigma }_{r}^{2}\;{k}_{r}^{2}}{{\kappa }^{2}}}\end{eqnarray} \tag{ 34 }$$

is a measure of how warm the disk is. In cases of interest the reduction factor is a number slightly smaller than unity and of little interest.

The proportionality (30) suggests that in a basis formed of tightly wound spiral waves the Fourier transformed response operator $\tilde{M}(\nu )$ is diagonal with ${\lambda }_{\boldsymbol{k}}$ the diagonal element associated with the given spiral wave. Hence the natural procedure might seem to be the adoption of the complete set of logarithmic spirals (e.g., Binney & Tremaine 2008, Section 2.6.3) as the basis ${\psi }^{(p)}$. Unfortunately, $\tilde{M}(\nu )$ is not, in fact, diagonal in the basis formed by logarithmic spirals for the following reason. The demonstration that a spiral perturbation generates a spiral response scaled by ${\lambda }_{\boldsymbol{k}}$ is a local result: the disk is analyzed in just an annulus, and in the spirit of WKB analysis, the wave considered is a packet of finite length. Since the frequencies κ and ${\rm{\Omega }}_{\phi }$ that appear in Equation (32) are functions of radius, ${\lambda }_{\boldsymbol{k}}$ is also a function of radius whereas a diagonal element of $\tilde{M}(\nu )$ should be a constant. Hence only a short packet of spiral waves provides a good approximation to an eigenfunction of $\tilde{M}(\nu )$. Such a packet is a non-trivial superposition of logarithmic spirals, so $\tilde{M}(\nu )$ cannot be diagonal in the basis provided by logarithmic spirals. Physically, the dynamics of the disk is inherently local on account of the existence of resonant radii, so basis functions such as logarithmic spirals that extend from the disk's center to infinity cannot make $\tilde{M}(\nu )$ diagonal. Given that we want $\tilde{M}(\nu )$ to become diagonal, we must work with basis functions ${\psi }^{(p)}$ that are local.

Fouvry et al. (2015) show how to construct a biothorgonal basis of localized spirals. They divide the range $({R}_{\mathrm{min}},{R}_{\mathrm{max}})$ of relevant radii into intervals of width σ centered on R₀, and then for any given wavevector $\boldsymbol{k}=({k}_{r},{k}_{\phi })$ create a basis function for each interval. Specifically, their basis potentials are

$$\begin{eqnarray}&&{\psi }^{(\boldsymbol{k},{R}_{0})}(R,\phi )=\sqrt{{\displaystyle \frac{G}{| {k}_{r}| {R}_{0}}}}{\displaystyle \frac{{e}^{i({k}_{r}R+{k}_{\phi }\phi )}}{{(\pi {\sigma }^{2})}^{1/4}}}\mathrm{exp}\left[-{\displaystyle \frac{{(R-{R}_{0})}^{2}}{2{\sigma }^{2}}}\right],\end{eqnarray} \tag{ 35 }$$

with ${k}_{r}{R}_{0}\gg 1$. The corresponding surface densities are

$$\begin{eqnarray}&&{\rm{\Sigma }}^{(\boldsymbol{k},{R}_{0})}=-{\displaystyle \frac{| {k}_{r}| }{2\pi G}}{\psi }^{(\boldsymbol{k},{R}_{0})}.\end{eqnarray} \tag{ 36 }$$

They show that two basis functions ${\psi }^{({\boldsymbol{k}}^{1},{R}_{0}^{1})}$ and ${\psi }^{({\boldsymbol{k}}^{2},{R}_{0}^{2})}$ will be biorthogonal only when $\rm{\Delta }{R}_{0}\equiv {R}_{0}^{1}-{R}_{0}^{2}$ and $\rm{\Delta }{k}_{r}={k}_{r}^{1}-{k}_{r}^{2}$ satisfy

$$\begin{eqnarray}\left\{\begin{array}{ccc}\rm{\Delta }{R}_{0} & \gg \sigma , & \rm{or}\;\rm{\Delta }{R}_{0}=0,\\ \rm{\Delta }{k}_{r} & \gg {\displaystyle \frac{1}{\sigma }}, & \rm{or}\;\rm{\Delta }{k}_{r}=0.\end{array}\right.\end{eqnarray} \tag{ 37 }$$

That is, the centers of neighboring bands have to be separated by more than the width of a band, and within any band, adjacent wavenumbers have to differ by enough to give a significant phase difference across the band.

Now that the basis potentials have been chosen, one can use the mapping $(\boldsymbol{x},\boldsymbol{v})\mapsto (\boldsymbol{\theta },\boldsymbol{J})$ provided by the epicycle approximation (e.g., Binney 2013, Equation (82)) to compute their Fourier transforms with respect to the angle variables:

$$\begin{eqnarray}\begin{array}{ccc}{\psi }_{\boldsymbol{m}}^{(\boldsymbol{k},{R}_{0})}(\boldsymbol{J}) & = & \;{\delta }_{{m}_{\phi }}^{{k}_{\phi }}\;{e}^{{{im}}_{r}{\theta }_{R}^{0}}\;\sqrt{{\displaystyle \frac{G}{| {k}_{r}| {R}_{0}}}}{\displaystyle \frac{1}{{(\pi {\sigma }^{2})}^{1/4}}}\;{e}^{{{ik}}_{r}{R}_{g}}\\ & & \times {\mathcal{J}}_{{m}_{r}}\left(\sqrt{\frac{2{J}_{r}}{\kappa }}{k}_{r}\right)\;\mathrm{exp}\left[-{\displaystyle \frac{{({R}_{g}-{R}_{0})}^{2}}{2{\sigma }^{2}}}\right],\end{array}\end{eqnarray} \tag{ 38 }$$

where ${R}_{g}({J}_{\phi })$ is the radius of the circular orbit with angular momentum ${J}_{\phi }$ and ${\mathcal{J}}_{{m}_{r}}$ is a Bessel function of the first kind. On account of the tight-winding condition ${k}_{r}{R}_{g}\gg 1$, the phase shift ${\theta }_{R}^{0}$ is given by

$$\begin{eqnarray}&&{\theta }_{R}^{0}\simeq -\pi /2.\end{eqnarray} \tag{ 39 }$$

In this basis, the response matrix $\tilde{M}$ is diagonal, having diagonal elements

$$\begin{eqnarray}&&{\tilde{M}}_{\left({\boldsymbol{k}}^{1},{R}_{0}^{1}\right)\left({\boldsymbol{k}}^{2},{R}_{0}^{2}\right)}={\delta }_{{k}_{r}^{1}}^{{k}_{r}^{2}}\;{\delta }_{{k}_{\phi }^{1}}^{{k}_{\phi }^{2}}\;{\delta }_{{R}_{0}^{1}}^{{R}_{0}^{2}}\;{\lambda }_{\boldsymbol{k}},\end{eqnarray} \tag{ 40 }$$

where ${\lambda }_{\boldsymbol{k}}$ is defined by Equation (31). This expression for $\tilde{\boldsymbol{M}}$ allows us to rewrite Equation (29) in the form

$$\begin{eqnarray}&&{\tilde{c}}_{\boldsymbol{m}}(\boldsymbol{J},\nu )={\displaystyle \sum }_{{pq}}{\psi }_{\boldsymbol{m}}^{(p)}(\boldsymbol{J}){{\psi }_{\boldsymbol{m}}^{(q)}}^{*}(\boldsymbol{J}){\displaystyle \frac{{A}_{{pq}}(\nu )}{(1-{\lambda }_{{\boldsymbol{k}}^{p}})(1-{\lambda }_{{\boldsymbol{k}}^{q}})}}.\end{eqnarray} \tag{ 41 }$$

In the chosen basis the expansion coefficients of the stimulating field are

$$\begin{eqnarray}\begin{array}{ccc}{a}_{p}(\nu ) & = & -{\displaystyle \int }{d}^{2}\boldsymbol{x}{\left[{\rm{\Sigma }}^{(p)}(\boldsymbol{x})\right]}^{*}{\tilde{\psi }}^{e}(\boldsymbol{x},\nu )\\ & = & \sqrt{{\displaystyle \frac{| {k}_{r}^{p}| }{{{GR}}_{0}^{p}}}}{\displaystyle \frac{1}{{(\pi {\sigma }^{2})}^{1/4}}}{\displaystyle \int }{dR}\;R\\ & & \times \mathrm{exp}\left[-{\displaystyle \frac{{(R-{R}_{0}^{p})}^{2}}{2{\sigma }^{2}}}\right]{e}^{-{{iRk}}_{r}^{p}}{\tilde{\psi }}_{{k}_{\phi }^{p}}^{e}(R,\nu ),\end{array}\end{eqnarray} \tag{ 42 }$$

where

$$\begin{eqnarray}&&{\tilde{\psi }}_{{k}_{\phi }^{p}}^{e}(R,\nu )\equiv {\displaystyle \frac{1}{2\pi }}\int d\phi \;{e}^{-{{ik}}_{\phi }^{p}\phi }{\tilde{\psi }}^{e}(R,\phi ,\nu )\end{eqnarray} \tag{ 43 }$$

is the Fourier transform in azimuthal angle and time of the stimulating potential. We further define the local radial Fourier transform of ${\psi }^{e}$ within the segment centered on ${R}_{0}^{p}$ by (Gabor 1946)

$$\begin{eqnarray}\begin{array}{ccc}{\tilde{\psi }}_{{\boldsymbol{k}}^{p}}^{e}({R}_{0}^{p},\nu ) & \equiv & {\displaystyle \frac{1}{2\pi }}{\displaystyle \int }{dR}\;\mathrm{exp}\left[-{\displaystyle \frac{{(R-{R}_{0}^{p})}^{2}}{2{\sigma }^{2}}}\right]\\ & & \times {e}^{-i(R-{R}_{0}^{p}){k}_{r}^{p}}{\tilde{\psi }}_{{k}_{\phi }^{p}}^{e}(R,\nu ).\end{array}\end{eqnarray} \tag{ 44 }$$

This definition is motivated by the consequence that thus defined the local radial Fourier transform of a uniform potential ${\psi }^{e}=1$ is independent of ${R}_{0}^{p}$. If in Equation (42) we approximate the leading factor R in the integrand by R₀, we then have

$$\begin{eqnarray}&&{\tilde{a}}_{p}(\nu )=\sqrt{{\displaystyle \frac{| {k}_{r}^{p}| {R}_{0}^{p}}{G}}}{\displaystyle \frac{2\pi }{{(\pi {\sigma }^{2})}^{1/4}}}{e}^{-{{iR}}_{0}^{p}{k}_{r}^{p}}{\tilde{\psi }}_{{\boldsymbol{k}}^{p}}^{e}({R}_{0}^{p},\nu ).\end{eqnarray} \tag{ 45 }$$

We require the ensemble average $\bar{{\tilde{a}}_{p}(\nu ){\tilde{a}}_{q}^{*}({\nu }^{\prime })}$ (Equations (28) and (29)), which is related to the ensemble average $\bar{{\psi }^{e}(\boldsymbol{x},t){\psi }^{e}(\boldsymbol{x}\prime ,t\prime )}$. We assume that stimulating fluctuations are quasi-stationary in the sense that

$$\begin{eqnarray}&&\bar{{\psi }_{{k}_{\phi }}^{e}(R,t){\psi }_{{k}_{\phi }}^{e*}(R\prime ,t\prime )}={C}_{{k}_{\phi }}(t-t\prime ,R-R\prime ,(R+R\prime )/2),\end{eqnarray} \tag{ 46 }$$

with the dependence on $R+R\prime$ being weak. With this assumption that the process ${\psi }_{{k}_{\phi }}^{e}$ is stationary in time and "locally stationary" in space, it follows that

$$\begin{eqnarray}\begin{array}{ccc}\bar{{\tilde{\psi }}_{{\boldsymbol{k}}^{p}}^{e}(R,\nu ){\tilde{\psi }}_{{\boldsymbol{k}}^{q}}^{e*}(R\prime ,{\nu }^{\prime })} & = & 2\pi \delta (\nu -{\nu }^{\prime })\delta ({k}_{r}^{p}-{k}_{r}^{q})\\ & & \times C({\boldsymbol{k}}^{p},\nu ,(R+R\prime )/2),\end{array}\end{eqnarray} \tag{ 47 }$$

where $C(\boldsymbol{k},\nu ,R)$ is the spatio-temporal power spectrum of the stimulating noise in the neighborhood of R. Now we can write

$$\begin{eqnarray}\begin{array}{ccc}\bar{{\tilde{a}}_{p}(\nu ){\tilde{a}}_{q}^{*}({\nu }^{\prime })} & = & {\displaystyle \frac{| {k}_{r}^{p}| {R}_{0}^{p}}{G}}{\displaystyle \frac{{(2\pi )}^{3}}{{(\pi {\sigma }^{2})}^{1/2}}}{e}^{-{{ik}}_{r}^{p}({R}_{0}^{p}-{R}_{0}^{q})}\delta (\nu -{\nu }^{\prime })\\ & & \times \delta ({k}_{r}^{p}-{k}_{r}^{q})C({\boldsymbol{k}}^{p},\nu ,({R}_{0}^{p}+{R}_{0}^{q})/2),\end{array}\end{eqnarray} \tag{ 48 }$$

so by Equation (28)

$$\begin{eqnarray}\begin{array}{ccc}{A}_{{pq}}(\nu ) & = & {\displaystyle \frac{| {k}_{r}^{p}| {R}_{0}^{p}}{G}}{\displaystyle \frac{{(2\pi )}^{2}}{{(\pi {\sigma }^{2})}^{1/2}}}{e}^{-{{ik}}_{r}^{p}({R}_{0}^{p}-{R}_{0}^{q})}\\ & & \times \delta ({k}_{r}^{p}-{k}_{r}^{q})C({\boldsymbol{k}}^{p},\nu ,({R}_{0}^{p}+{R}_{0}^{q})/2).\end{array}\end{eqnarray} \tag{ 49 }$$

To obtain the diffusion coefficients from Equation (7) we substitute into Equation (41) for ${\tilde{c}}_{\boldsymbol{m}}$ our expressions (38) for ${\psi }_{\boldsymbol{m}}^{(p)}(\boldsymbol{J})$ and the above expression for $\boldsymbol{A}$. We have

$$\begin{eqnarray}\begin{array}{ccc}{\tilde{c}}_{\boldsymbol{m}}(\boldsymbol{J},\nu ) & = & {\displaystyle \frac{{(2\pi )}^{2}}{\pi {\sigma }^{2}}}{\displaystyle \sum }_{{pq}}{\delta }_{{m}_{\phi }}^{{k}_{\phi }^{p}}{\delta }_{{m}_{\phi }}^{{k}_{\phi }^{q}}{\mathcal{J}}_{{m}_{r}}^{2}\left(\sqrt{{\displaystyle \frac{2{J}_{r}}{\kappa }}}{k}_{r}^{p}\right){e}^{-{{ik}}_{r}^{p}({R}_{0}^{p}-{R}_{0}^{q})}\\ & & \times \mathrm{exp}\left[-{\displaystyle \frac{{({R}_{g}-{R}_{0}^{p})}^{2}+{({R}_{g}-{R}_{0}^{q})}^{2}}{2{\sigma }^{2}}}\right]\\ & & \times \delta ({k}_{r}^{p}-{k}_{r}^{q}){\displaystyle \frac{C({\boldsymbol{k}}^{p},\nu ,({R}_{0}^{p}+{R}_{0}^{q})/2)}{{(1-{\lambda }_{{\boldsymbol{k}}^{p}})}^{2}}}.\end{array}\end{eqnarray} \tag{ 50 }$$

The sums over p and q expand into sums over ${\boldsymbol{k}}^{p}$, ${\boldsymbol{k}}^{q}$, ${R}_{0}^{p}$ and ${R}_{0}^{q}$. On account of the factors ${\delta }_{{m}_{\phi }}^{{k}_{\phi }^{p}}$ and ${\delta }_{{m}_{\phi }}^{{k}_{\phi }^{q}}$ the sums over the ϕ components of the ${\boldsymbol{k}}^{i}$ are trivial. The sums over the ${k}_{r}^{i}$ and ${R}_{0}^{i}$ we approximate with integrals by the substitutions

$$\begin{eqnarray}\begin{array}{ccc}{\displaystyle \sum }_{{k}_{r}}f({k}_{r}) & \to & {\displaystyle \frac{1}{\rm{\Delta }{k}_{r}}}{\displaystyle \int }{{dk}}_{r}\;f({k}_{r}),\\ {\displaystyle \sum }_{{R}_{0}}g({R}_{0}) & \to & {\displaystyle \frac{1}{\rm{\Delta }{R}_{0}}}{\displaystyle \int }{{dR}}_{0}\;g({R}_{0}),\end{array}\end{eqnarray} \tag{ 51 }$$

where $\rm{\Delta }{k}_{r}$ is the difference between successive values of k_r in the sum, and similarly for $\rm{\Delta }{R}_{0}$. Then the Dirac delta function in Equation (50) allows us to integrate over ${k}_{r}^{q}$. We assume that σ is small enough that we can approximate each Gaussian exponential by $\sqrt{2\pi }\sigma \delta ({R}_{g}-{R}_{0}^{i})$ so we can trivially integrate over the ${R}_{0}^{i}$. Finally, the presence in Equation (48) of a rapidly oscillating complex exponential ${e}^{{{ik}}_{r}{R}_{0}}$ imposes that the intervals $\rm{\Delta }{k}_{r}$ and $\rm{\Delta }{R}_{0}$ must satisfy the critical-sampling condition $\rm{\Delta }{k}_{r}^{i}\rm{\Delta }{R}_{0}^{i}=2\pi$ (Daubechies 1990). With this condition, Equation (50) simplifies to

$$\begin{eqnarray}&&{\tilde{c}}_{\boldsymbol{m}}(\boldsymbol{J},\nu )=2\int {{dk}}_{r}\;{\displaystyle \frac{{\mathcal{J}}_{{m}_{r}}^{2}\left(\sqrt{2{J}_{r}/\kappa }\;{k}_{r}\right)}{{(1-{\lambda }_{\boldsymbol{k}})}^{2}}}\;C(\boldsymbol{k},\nu ,{R}_{g}),\end{eqnarray} \tag{ 52 }$$

where by hypothesis the dependence of C on ${R}_{g}$ is weak, and we have ${k}_{\phi }={m}_{\phi }$. This is our principal result. It enables us to compute the diffusion tensor at any point in the action space of a thin, self-gravitating disk given the power spectrum of the noise that excites spiral structure in the disk.

To proceed further we need to assume some form for the power spectrum $C(\boldsymbol{k},\nu ,{R}_{g})$ that appears in Equation (52). An inevitable source of noise is shot noise caused by the the finite number of stars in the disk, and massive, compact gas clouds are a source of spectrally similar noise, so let us investigate the impact that shot noise has. In this case the power spectrum is independent of ν and k_r, and varies with radius like $\surd \rm{\Sigma }(R)$. Then $C\propto \rm{\Sigma }$, so to within a normalization that depends on particle number, we have

$$\begin{eqnarray}&&{\tilde{c}}_{\boldsymbol{m}}(\boldsymbol{J},\nu )={\rm{\Sigma }}_{t}({J}_{\phi })\int {{dk}}_{r}\;{\displaystyle \frac{{\mathcal{J}}_{{m}_{r}}^{2}\left(\sqrt{2{J}_{r}/\kappa }\;{k}_{r}\right)}{{(1-{\lambda }_{\boldsymbol{k}})}^{2}}}.\end{eqnarray} \tag{ 63 }$$

The eigenvalues ${\lambda }_{\boldsymbol{k}}$ have to be evaluated at $\nu =\boldsymbol{m}\cdot \bf{\Omega }$, and then $s={m}_{r}$ by Equation (32). In order to handle the singularity of the Equation (31) when $s=\pm 1$, one adds a small imaginary part to the frequency of evaluation, so that $s={m}_{r}+i\eta$. As long as η in modulus is small compared to imaginary part of the least damped mode of the disk, adding this complex part makes a negligible contribution on the expression of $\mathrm{Re}(\lambda )$.

In Equation (63) the integral over k_r should formally be over the full range 0 to $\infty$. However, small values of k_r are unphysical and violate our assumption of tightly wound spirals. Values of k_r that are larger than $\sim 2\pi$ divided by the thickness of a galactic disk are also unphysical, and in the case ${m}_{r}=0$ of the CR the integral diverges at ${J}_{r}=0$ since then the Bessel function remains non-zero to arbitrarily high k_r and ${(1-{\lambda }_{\boldsymbol{k}})}^{-2}$ is always greater than unity. Hence we must determine appropriate upper and lower limits to the integration on k_r.

At any point in action space the biggest contribution to the diffusion tensor will come from waves that yield the largest value of ${\lambda }_{\boldsymbol{k}}$. Hence we now examine the structure of the function ${k}_{r}\mapsto {\lambda }_{\boldsymbol{k}}$ for given $\boldsymbol{m}$ and $\boldsymbol{J}$. Figure 4 shows that ${\lambda }_{\boldsymbol{k}}$ has a well-defined peak at a value ${k}_{\mathrm{max}}$ of k_r that decreases as ${J}_{\phi }$ increases. In fact ${k}_{\mathrm{max}}\propto 1/{J}_{\phi }$ because the radius of a near-circular orbit is $R\propto {J}_{\phi }$ and in a scale-free model we expect ${k}_{\mathrm{max}}\propto 1/R$. For the same reason we expect the width ${\rm{\Delta }}_{k}$ of the peak in ${\lambda }_{\boldsymbol{k}}$ to be proportional to $1/{J}_{\phi }$. In light of these observations we adopt as lower and upper bounds on the k_r integral the wavenumbers ${k}_{\mathrm{inf}}$ and ${k}_{\mathrm{sup}}$ defined by

$$\begin{eqnarray}&&{\lambda }_{{k}_{\mathrm{inf}}}={\lambda }_{{k}_{\mathrm{sup}}}={\displaystyle \frac{1}{2}}{\lambda }_{{k}_{\mathrm{max}}}.\end{eqnarray} \tag{ 64 }$$

These limiting values of k_r are marked on Figure 4.

Zoom In Zoom Out Reset image size

At each point in action space there are contributions to the diffusion coefficients from several values of $\boldsymbol{m}$. The contribution from ${m}_{r}=-1$ is driven by waves that have their inner Lindblad resonance at that point, while that from ${m}_{r}=+1$ is driven by waves that have their outer Lindblad resonance there, and the contribution from ${m}_{r}=0$ is driven by waves locally in corotation. Since ${\lambda }_{\boldsymbol{k}}$ depends on ${s}^{2}={m}_{r}^{2}$, the value of ${\lambda }_{\boldsymbol{k}}$ is the same for ${m}_{r}=\pm 1$. At a given point in action space different values of the frequency ν are associated with ${m}_{r}=\pm 1$, but since we are considering shot noise, the fluctuations have the same power at all frequencies. Hence the values ${m}_{r}=\pm 1$ contribute equally to the diffusion coefficients. The lower curve in Figure 5 shows the extent to which this contribution is amplified by the disk's self-gravity and we see that the amplification is modest. The upper curve in Figure 5 shows the amplification by self-gravity in the case ${m}_{r}=0$ of corotation: it is much larger.

Zoom In Zoom Out Reset image size

In a cool disk, $| \partial f/\partial {J}_{r}| \gg | \partial f/\partial {J}_{\phi }|$ so by Equation (5) a reasonable approximation to the diffusive flux, when ${m}_{r}\ne 0$, is

$$\begin{eqnarray}&&{F}_{i}=-{\displaystyle \frac{1}{2}}{\displaystyle \sum }_{\boldsymbol{m}}{m}_{i}{\tilde{c}}_{\boldsymbol{m}}{m}_{r}{\displaystyle \frac{\partial {f}_{0}}{\partial {J}_{r}}}.\end{eqnarray} \tag{ 65 }$$

It follows that waves at inner Lindblad resonance, which couple through $({m}_{r},{m}_{\phi })=(-1,2)$ drive a diffusive flux

$$\begin{eqnarray}&&\boldsymbol{F}=-(-1,2){\displaystyle \frac{1}{2}}{\tilde{c}}_{(-\mathrm{1,2})}\left|{\displaystyle \frac{\partial {f}_{0}}{\partial {J}_{r}}}\right|,\end{eqnarray} \tag{ 66 }$$

while waves at outer Lindblad resonance drive the flux

$$\begin{eqnarray}&&\boldsymbol{F}=(1,2){\displaystyle \frac{1}{2}}{\tilde{c}}_{(\mathrm{1,2})}\left|{\displaystyle \frac{\partial {f}_{0}}{\partial {J}_{r}}}\right|.\end{eqnarray} \tag{ 67 }$$

These formulae show that waves with a Lindblad resonance at $\boldsymbol{J}$ drive a flux toward increasing J_r: these waves are heating the disk by increasing the eccentricities of orbits. Wave that have their ILR at $\boldsymbol{J}$ drive stars to lower angular momentum, while those with their OLR at $\boldsymbol{J}$ drive stars to higher angular momentum. We shall see below that at low ${J}_{\phi }$ waves with a local ILR dominate, so on average angular momenta decrease, while at large ${J}_{\phi }$ the waves with a local OLR dominate and angular momenta increase on average. Thus spiral structure conveys angular momentum outwards, thus liberating energy that is used to heat the disk.

Waves at corotation couple through $\boldsymbol{m}=(0,2)$ so (a) they can only drive diffusion parallel to the ${J}_{\phi }$ axis, and (b) that diffusion is proportional to the gradient's small component $\partial f/\partial {J}_{\phi }$. In practice this diffusion is important only in so far as the disk has a metallicity gradient, when the DF of stars of any particular metallicity can have a significant derivative with respect to ${J}_{\phi }$ and the large value of ${\tilde{c}}_{\boldsymbol{m}}$ at corotation can have a big impact on the metallicity distribution—radial migration at corotation is important for chemical evolution (Sellwood & Binney 2002; Schönrich & Binney 2009).

In this section we examine the extent to which our analytic results can explain the simulations of tapered Mestel disks described by S12. We do not expect perfect agreement between our results and the numerical experiments because we have employed several approximations. In particular we have assumed that the driving fluctuations are white noise that has an amplitude squared that is proportional to the disk's surface density. We have assumed, moreover, that these fluctuations drive tightly wound waves. S12 restricted disturbing forces to ${m}_{\phi }=2$, so we impose this same restriction on $\boldsymbol{m}$.

The dark contours in Figure 6 show the magnitude of the diffusive flux that is generated by the two Lindblad resonances and corotation when one adopts the same numerical pre-factors as S12. The gray arrow shows the direction of the diffusive flux at the location of the peak flux; the direction is similar at neighboring points. The thin contours in Figure 6 show the value of the DF in S12 after diffusion has taken place. Originally the DF peaked along the ${J}_{\phi }$ axis at ${J}_{\phi }\simeq 1$, becoming small to the left of that point on account of the inner taper. Above the location of the peak DF, there is a prominent horn of enhanced stellar density that extends roughly in the direction of the diffusive flux. Thus our analytic results are in good qualitative agreement with the numerical experiments of S12.

Zoom In Zoom Out Reset image size

In Figure 6 waves at ILR drive diffusion parallel to vectors that are inclined by 153° to the ${J}_{\phi }$ axis, whereas the net flux makes an angle of 111° with the axis. That these angles are similar attests to the dominant role of waves that have ILR near where the DF peaks. This dominance arises because $\partial {f}_{0}/\partial {J}_{r}$ is always negative and in this region $\partial {f}_{0}/\partial {J}_{\phi }\gt 0$, so $| \boldsymbol{m}\cdot \partial {f}_{0}/\partial \boldsymbol{J}|$ is larger for the waves at ILR, which have ${m}_{r}=-1$, than for the waves at OLR, while ${\tilde{c}}_{\boldsymbol{m}}$ is the same for both values of m_r.

Given that we have assumed that the driving fluctuations are white noise, it is perhaps surprising that the magnitude of the diffusive flux is as sharply peaked in action space as the dark contours in Figure 6 show it to be. This localization surely reflects the sharp tapering of the DF at small radii. Just outside this taper the disk becomes significantly self-gravitating and most effective at amplifying the stimulating white noise. The amplified waves tend to have their ILRs further in, where the taper is cutting into the disk.

We have seen that self-gravity amplifies stimuli most strongly at corotation (Figure 5). In fact, the dominance of corotation increases without limit as Q decreases because the value of ${\lambda }_{\mathrm{max}}$ associated with corotation approaches unity, and the associated amplification diverges, while the value of ${\lambda }_{\mathrm{max}}$ associated with the LRs remains significantly below unity. Figure 7 illustrates this phenomenon by comparing the velocity of diffusion in action space,

$$\begin{eqnarray}&&\boldsymbol{v}={\displaystyle \frac{\boldsymbol{F}}{f}},\end{eqnarray} \tag{ 68 }$$

for two values of ξ. For the value $\xi =0.5$ assumed by S12 (left panel) the velocity vectors are significantly non-zero only within the narrow strip at ${J}_{\phi }\sim 1$ associated with the ILR and point up and to the left. For $\xi =0.73$ the velocity is also quite large at ${J}_{\phi }\gt 2$ near the ${J}_{\phi }$ axis and is there directed up and to the right.

Zoom In Zoom Out Reset image size