© The Authors 2025

1. Introduction

Spiral density perturbations are ubiquitous in late-type galaxies (by definition), and they are an important ingredient for the formation and evolution of these galaxies because the resulting non-axisymmetric gravitational forces alter the otherwise simple dynamics of the interstellar medium (ISM) and of the stars. The ISM is often strongly affected and forced to high densities, which triggers star formation, while the perturbation of the stellar orbits leads to radial migration and radial heating (Sellwood & Binney 2002), which correspond to diffusion in the orbital actions J_ϕ and J_R, respectively.

Studies of these processes often require some model for the gravitational potential of the spiral perturbations, and in many situations, it is only necessary to model the dynamics in two dimensions without considering the vertical motions. Unfortunately, the only known gravitational potential for a flat (razor-thin) spiral density perturbation is that developed by Kalnajs’ (1971), which has an amplitude that declines with radius like R^−3/2 and is unsuitable as a model for galactic spirals. Instead, the tight-winding approximation (often also called ‘WKB approximation’¹) allows reasonably accurate estimates for the gravitational potential of a spiral pattern with a pitch angle α ≲ 20°.

Dynamical modellers typically employ simple models for the gravitational potential (e.g. Roberts et al. 1979; Contopoulos & Grosbøl 1986; Eilers et al. 2020; Hamilton et al. 2024) that are often based on this approximation. However, the actual density corresponding to these models has never been computed. Similarly, the actual potential generated by reasonably realistic spiral density perturbations is largely unknown.

The purpose of this study is to fill these gaps. To this end, I implement in Appendix B an efficient numerical tool for finding the gravitational potential for any flat density and, vice versa, the density from the potential known in the equatorial plane only. This tool is then applied to assess the accuracy of the tight-winding approximations, including a more accurate novel version, for logarithmic spirals with amplitude that declines with radius either as a power law or exponentially (like the surface density of disc galaxies) with and without radial truncation. I also obtain the actual density of models for spiral potentials that have been used in previous studies, and find that the density sometimes behaves in unexpected ways. For example, the pitch angle can deviate significantly from the target.

Finally, I consider the gravitational torques of a spiral perturbation and derive an approximate expression for their local mean at each radius. For trailing spirals, this is always negative at small and positive at large radii, implying an outward angular-momentum transport (Zhang 1996), in addition to the advective angular momentum transport by diffusion driven by the standard deviation of the torques.

I begin the presentation in Section 2 by introducing my notations and giving a brief overview over known results, including the tight-winding approximation. In Section 3, I extend Kalnajs’ (1971) potential-density pairs to logarithmic spirals whose density amplitude is an arbitrary power law in radius. I also develop a novel more accurate tight-winding approximation. Section 4 considers the non-local gravitational torque induced by a spiral perturbation. In Section 5, I consider various models for spiral perturbations and compare their numerically computed properties to the predictions of the tight-winding approximations. Finally, Sections 6 and 7 discuss my findings and present my conclusion.

2. Spiral density perturbations

This section introduces basic concepts and notations and summarises the important tight-winding approximation, but presents no novel results.

2.1. Poisson’s equation for razor-thin systems

This study exclusively considers the situation of razor-thin mass distributions with some surface density Σ(R), where R ≡ (x, y), and the spatial density ρ(r) = Σ(R) δ(z), where r ≡ (x, y, z) = (R, z). The gravitational potential

$$\begin{array}{c}\hfill \mathrm{\Phi }(\mathit{r})=-G{\displaystyle \int \int \int \frac{\rho ({\mathit{r}}^{\prime }){\mathrm{d}}^3{\mathit{r}}^{\prime }}{|\mathit{r}-{\mathit{r}}^{\prime }|}}\end{array}$$(1)

satisfies the Poisson equation ΔΦ(r) = 4πG ρ(r), which implies

$$\begin{array}{c}\hfill \mathrm{\Sigma }(\mathit{R})=\frac{1}{2\pi G}{\displaystyle \underset{z\to 0^{+}}{lim}}\frac{\partial \mathrm{\Phi }}{\partial z}.\end{array}$$(2)

However, often only the potential in the equatorial plane Ψ(R)≡Φ(R, z = 0) is required or even known, when Eq. (1) becomes

$$\begin{array}{c}\hfill \mathrm{\Psi }(\mathit{R})=-G{\displaystyle \int \int \frac{\mathrm{\Sigma }({\mathit{R}}^{\prime }){\mathrm{d}}^2{\mathit{R}}^{\prime }}{|\mathit{R}-{\mathit{R}}^{\prime }|}}.\end{array}$$(3)

The inverse operation, obtaining Σ(R) from Ψ(R), is not as straightforward as for three-dimensional models because Eq. (2) cannot be used. Nonetheless, a useful relation can be inferred from the fact that a razor-thin plane density wave Σ ∝ e^ik ⋅ R has the potential Φ(r) = Ψ(R)e^−k|z|, where k = |k| and

$$\begin{array}{c}\hfill \mathrm{\Psi }(\mathit{R})=-2\pi G\mathrm{\Sigma }(\mathit{R})/k.\end{array}$$(4)

It follows that in Fourier space, Eq. (3) becomes $\widehat{\mathrm{\Psi }}=-2\pi G\widehat{\mathrm{\Sigma }}/k$ or $2\pi G\widehat{\mathrm{\Sigma }}=-k\widehat{\mathrm{\Psi }}=\widehat{\mathrm{\Delta }\mathrm{\Psi }}/k$, where a hat denotes the Fourier transform. Thus, −ΔΨ relates to 2πGΣ in exactly the same way as 2πGΣ relates to Ψ, and therefore,

$$\begin{array}{c}\hfill \mathrm{\Sigma }(\mathit{R})=\frac{1}{4{\pi }^{2}G}{\displaystyle \int \int \frac{\mathrm{\Delta }\mathrm{\Psi }({\mathit{R}}^{\prime }){\mathrm{d}}^2{\mathit{R}}^{\prime }}{|\mathit{R}-{\mathit{R}}^{\prime }|}}\end{array}$$(5)

(H. Dejonghe in Binney & Tremaine 2008, problem 2.22; see Appendix A for a proof without Fourier transform). Thus, obtaining Σ(R) from Ψ(R) requires a similar effort as the inverse relation and can be achieved with the same numerical or analytic tools. Appendix B details the numerical tool I used in this study.

2.2. Spiral patterns

The density and potential of the disc can be uniquely decomposed into azimuthal Fourier components of the general form

$$\begin{array}{cc}\hfill \mathrm{\Sigma }(\mathit{R})& =S_m(R){\mathrm{e}}^{\mathrm{i}m[\varphi -{\psi }_{\mathrm{\Sigma }}(R)]},\hfill \end{array}$$(6a)

$$\begin{array}{cc}\hfill \mathrm{\Psi }(\mathit{R})& =-P_m(R){\mathrm{e}}^{\mathrm{i}m[\varphi -{\psi }_{\mathrm{\Psi }}(R)]},\hfill \end{array}$$(6b)

with amplitudes P_​m, S_​m ≥ 0 and phases ψ_Ψ, ψ_Σ, which are real-valued functions of the radius R = |R|. As usual, the physical potential and density are given by the real parts of Eqs. (6), and we can assume m ≥ 0 without loss of generality. For a spiral pattern, the phase is a monotonously increasing (or decreasing) function of radius. The pitch angle α > 0 between the wave crest at ϕ = ψ(R) and concentric circles relates to the phase via

$$\begin{array}{c}\hfill \lambda \equiv \frac{\mathrm{d}\psi }{\mathrm{d}lnR}=\frac{\sigma }{tan\alpha },\end{array}$$(7)

where σ = ±1 gives the sense in which the spiral twists, such that α > 0. Galactic spirals are trailing, implying $\sigma =-\mathit{sign}({\overline{v}}_{\varphi })$, where ${\overline{v}}_{\varphi }$ is the mean azimuthal (rotational) velocity. The pitch angles of galactic spirals tend to vary with radius (Kennicutt 1981; Savchenko & Reshetnikov 2013; Chugunov et al. 2025), but a widespread model is that of a constant pitch angle, which is obtained for

$$\begin{array}{c}\hfill \psi (R)=\lambda ln(R/R_0)\end{array}$$(8)

and referred to as a ‘logarithmic spiral’. In this case, the phase factors in Eqs. (6) become e^{im[ϕ − λln(R/R₀)]}, which is a plane wave in both ϕ and lnR.

Except for models with power-law amplitudes (see Section 3), no analytic solutions to Eq. (3) are known for spirals of the general form (6a). Instead, the following approximation is often employed.

2.3. The tight-winding (or WKB) approximation

Each annulus of a disc with, say, an m = 2 perturbation generates a quadrupolar gravitational field. For a bar-like perturbation, these quadrupoles have the same orientation, and their contributions to the field at some position accumulate. For a spiral, on the other hand, the orientations vary with radius, and most of the contributions from distant annuli to the local field cancel (destructive interference), such that it is dominated by the local pattern. To estimate the field from the local pattern, the curvature of the spiral can be ignored and it can be approximated as a plane wave. In analogy to the relation (4), this gives (Toomre 1964; Lin & Shu 1964)

$$\begin{array}{cc}\hfill \mathrm{\Psi }(\mathit{R})& \approx -\frac{2\pi G}{k(R)}\mathrm{\Sigma }(\mathit{R}),\hfill \end{array}$$(9)

where k(R) is a locally appropriate horizontal wavenumber. Following Lin & Shu (1964), k(R) is traditionally taken to be just the radial wavenumber

$$\begin{array}{c}\hfill k(R)=|k_R|=\frac{m|\lambda |}{R}=\frac{m}{Rtan\alpha }.\end{array}$$(10a)

However, the analogy to plane waves suggests the absolute value of the horizontal wavevector, which also includes the azimuthal wavenumber k_ϕ = m/R, i.e.

$$\begin{array}{c}\hfill k(R)=\sqrt{k_R^2+k_{\varphi }^2}=\frac{m\sqrt{{\lambda }^2+1}}{R}=\frac{m}{Rsin\alpha }.\end{array}$$(10b)

Thus, the potential is weaker for more arms (larger m) and more tightly wound spirals (smaller α). Apparently, Eq. (10b) has not been used explicitly, but is implicit in shearing-sheet approaches for modelling spiral structure (e.g. Julian & Toomre 1966, Eq. (11)). Both Eqs. (10) are first-order accurate in the sense that the fractional error made in Eq. (9) is 𝒪(sin α) and the total error 𝒪(sin²​α), but Eq. (10b) is significantly better than Eq. (10a), as demonstrated in Figs. 2 and 3.

In general, the phases ψ_Ψ and ψ_Σ of the potential and density of a spiral perturbation differ, but the first-order approximations (10) give ψ_Ψ = ψ_Σ and hence do not predict the phase offset δψ ≡ ψ_Ψ − ψ_Σ. Shu’s (1970) second-order approximation

$$\begin{array}{c}\hfill \mathrm{\Sigma }(\mathit{R})=-\frac{\mathrm{i}\sigma }{2\pi G\sqrt{R}}\frac{\partial \sqrt{R}\mathrm{\Psi }}{\partial R}\times [1+O({sin}^2\alpha )],\end{array}$$(11)

predicts δψ, but does not provide Ψ(R) for a given Σ(R).

3. Scale-invariant spirals

A simple situation is that of a logarithmic spiral whose density amplitude follows a power law in radius,

$$\begin{array}{cc}\hfill \mathrm{\Sigma }(R,\varphi )& =S_{m,0}{(R/R_0)}^{-\gamma }{\mathrm{e}}^{\mathrm{i}m[\varphi -\lambda ln(R/R_0)]},\hfill \end{array}$$(12)

such that the total R dependence, Σ ∝ R^{−γ − imλ}, is still a power law and hence scale-invariant. The gravitational potential of these models is (see Appendix D for a derivation)

$$\begin{array}{cc}\hfill \mathrm{\Phi }(\mathit{r})& =-2\pi G{f}_m(|z|/r)r\mathrm{\Sigma }(r,\varphi ),\hfill \end{array}$$(13)

with spherical radius r = |r| and

$$\begin{array}{cc}\hfill f_m(x)& =-{(\mathrm{d}{P}_{-1/2-\mathrm{i}\zeta }^{|m|}/\mathrm{d}x)}_{x=0}^{-1}{P}_{-1/2-\mathrm{i}\zeta }^{|m|}(x),\hfill \end{array}$$(14)

where

$$\begin{array}{cc}\hfill \zeta & \equiv m\lambda -\mathrm{i}(\gamma -\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$2$}\right.)\hfill \end{array}$$(15)

and P_ν^m(x) denotes the associated Legendre function of the first kind. The function f_m describes the run of Φ(r) with latitude and at x → 1 declines like (1 − x)^m/2.

This study focusses on the potential in the z = 0 plane,

$$\begin{array}{c}\hfill \mathrm{\Psi }(\mathit{R})=-2\pi G{K}_m(\zeta )R\mathrm{\Sigma }(\mathit{R}),\end{array}$$(16)

where

$$\begin{array}{cc}\hfill K_m(\zeta )& \equiv f_m(0)=\frac{1}{2}\frac{\mathrm{\Gamma }(\frac{1}{4}+\frac{m}{2}+\mathrm{i}\frac{\zeta }{2})\mathrm{\Gamma }(\frac{1}{4}+\frac{m}{2}-\mathrm{i}\frac{\zeta }{2})}{\mathrm{\Gamma }(\frac{3}{4}+\frac{m}{2}+\mathrm{i}\frac{\zeta }{2})\mathrm{\Gamma }(\frac{3}{4}+\frac{m}{2}-\mathrm{i}\frac{\zeta }{2})},\hfill \end{array}$$(17)

which satisfies K_m(−ζ) = K_m(ζ), K_m(ζ^*) = K_m^*(ζ) (where a star denotes complex conjugation), and

$$\begin{array}{cc}\hfill K_m(\zeta )K_{m+1}(\zeta )& ={[{\zeta }^2+(m+\frac{1}{2}{)}^2]}^{-1},\hfill \end{array}$$(18)

$$\begin{array}{cc}\hfill K_m(\zeta )K_m(\zeta +\mathrm{i})& ={[(\zeta +\frac{\mathrm{i}}{2}{)}^2+m^2]}^{-1}.\hfill \end{array}$$(19)

Interpolating these relations gives the asymptotic behaviour

$$\begin{array}{cc}\hfill K_m(\zeta )& \approx {[{\zeta }^2+m^2-\frac{1}{4}]}^{-1/2},\hfill \end{array}$$(20)

which is an excellent approximation even for pitch angles as large as 40°, as Fig. 1 demonstrates.

Fig. 1. Assessing the accuracy of the approximation (20) for m = 2 (for m > 2 the accuracy is better). The real (imaginary) parts are even (odd) functions of $\gamma -\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$ (hence the relations for $\gamma <\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$ are hidden in the top panel). I also show (dashed) the relation given by Kalnajs (1971).

For these scale-invariant spirals, the potential Ψ and density Σ are proportional to the same phase factor e^{im[ϕ − λln(R/R₀)]}, such that they have the same constant pitch angle α. For $\gamma =\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$, ζ is real-valued, and hence, so are f_m(x) and K_m(ζ), such that Ψ and Σ even have the same phase ψ(R). In this case, S_​m ∝ R^−3/2 and P_​m ∝ R^−1/2, and we obtain the potential-density pairs of Kalnajs (1971), who for real-valued ζ gave Ψ(R) with Eqs. (17,18) and an inferior version of Eq. (20).

For $\gamma \ne \raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$, however, K_m(ζ) is complex-valued and causes a constant phase offset δψ = −arg(K_m(ζ))/m of Ψ with respect to Σ (see also Fig. 2), such that for a trailing spiral, the potential lags the density for $\gamma <\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$ and leads it for $\gamma >\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$ (Zhang 1996).

Fig. 2. Assessing the tight-winding (WKB) approximation for scale-invariant spirals. The relative error of the density amplitude and the error in the phase offset δψ ≡ ψΨ − ψΣ are plotted vs. pitch angle α for m = 2 and various values of the exponent γ (for these scale-invariant models, the errors are the same at each radius). Since the first-order approximations give δψapprox = 0, their error reflects the actual phase offset of the models.

As shown in the next section, this constant offset would imply a non-vanishing total torque in contradiction to Noether’s theorem. However, when the total torque is expressed as an integral over pairs of annuli torquing each other, it is found to be zero. This contradiction arises because of the unphysical nature of these models at R → 0 (see Appendix E for details).

3.1. Comparison with the tight-winding approximation

The scale-invariant spirals have no inner or outer edge, such that there is a local spiral field everywhere. This makes these models the ideal situation for the tight-winding approximation. Fig. 2 plots the errors in amplitude and phase made by the first- (10) and second-order approximation (11), respectively. The traditional tight-winding approximations (10a) and (11) underestimate the density or, conversely, overestimate the potential and forces, but the first-order approximation (10b) only makes very small errors in amplitude. The order of the approximations is only obvious in the phase error (a peculiarity of these power-law models), where the first-order methods predict δψ = 0, while the phase error of the second-order approximation is quadratic in the pitch angle, as expected.

For $\gamma =\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$ (red; Kalnajs’ 1971 spirals), ψ_Ψ = ψ_Σ and the approximation using Eq. (10b) is identical to Eq. (20) and hence extremely accurate. Away from the z = 0 plane, however, the tight-winding approximation deteriorates, as it predicts ψ_Ψ(R, z) = ψ_Σ(R), while actually ${\psi }_{\mathrm{\Psi }}(R,z)={\psi }_{\mathrm{\Sigma }}(\sqrt{R^2+z^2})$, i.e. the phase is in fact constant on spheres, not on cylinders. A more detailed assessment of the potential and its approximation away from the plane is beyond the scope of this study.

3.2. An improved tight-winding approximation

In the traditional tight-winding approximation, a spiral is locally approximated as plane wave, i.e. assumed to have constant amplitude and pitch angle α ∝ R. The amplitude might instead be assumed to decline locally like a power law and the pitch angle be constant. In other words, the spiral is approximated locally as a scale-invariant spiral. When its potential is in turn approximated via Eq. (20), it can be expressed in the standard form (9) using the complex-valued wavenumber

$$\begin{array}{c}\hfill k(R)=\frac{1}{R}\sqrt{{\zeta }^2(R)+m^2-\frac{1}{4}},\end{array}$$(21)

where the principal value of the square root is taken and

$$\begin{array}{cc}\hfill \zeta (R)& =\mathrm{i}[\frac{3}{2}+\frac{R}{\mathrm{\Sigma }}\frac{\partial \mathrm{\Sigma }}{\partial R}]=m\lambda (R)-\mathrm{i}(\gamma (R)-\frac{3}{2}),\hfill \end{array}$$(22)

with λ(R) as defined in Eq. (7) and γ(R)≡ − dlnS_​m/dlnR.

It is also possible to set $\zeta (R)=\text{i}[\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.+(R/\mathrm{\Psi })(\partial \mathrm{\Psi }/\partial R)]$, which allows the approximation of Σ(R) given Ψ(R), and, when also approximating $\sqrt{{\zeta }^2+m^2-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$4$}\right.}\approx \sigma \zeta$, gives the second-order approximation (11). A somewhat more simply computable and manipulable form than Eq. (21) is provided by the Taylor expansion,

$$\begin{array}{cc}\hfill \frac{1}{k(R)}\approx \frac{R}{\sqrt{{\zeta }^2+m^2}}& =\frac{R}{m\sqrt{{\lambda }^2+1}}[1+\frac{\mathrm{i}\lambda ϵ}{\sqrt{{\lambda }^2+1}}+\mathcal{O}{ϵ}^2],\hfill \end{array}$$(23)

where $ϵ\equiv (\gamma -\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$2$}\right.)/m\sqrt{{\lambda }^2+1}=sin\alpha (\gamma -\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$2$}\right.)/m$. When this is truncated at zeroth order in ϵ (or for $\gamma =\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$ when ϵ = 0), this collapses to the first-order approximation (10b). These relations suggest that the novel approximation (21) is at least second-order accurate. In terms of the pitch angle α (instead of λ), the approximation (9) resulting from the expansion (23) is conveniently expressed as

$$\begin{array}{cc}\hfill \mathrm{\Psi }(\mathit{R})& \approx -2\pi G\frac{sin\alpha }{m}[1+\mathrm{i}\sigma \frac{sin2\alpha }{2m}(\gamma -\frac{3}{2})]R\mathrm{\Sigma }(\mathit{R}).\hfill \end{array}$$(24)

Section 5 assesses this together with the first-order approximations for various spiral models.

4. The gravitational torque of spiral perturbations

A spiral perturbation exerts a gravitational torque per unit mass of −∂Φ/∂ϕ. This induces changes in the stellar angular momenta, which cause radial migration (Sellwood & Binney 2002) and are particularly pronounced for stars that co-rotate with the spiral pattern (and for stars on the Lindblad resonances) because they experience coherent torquing, while the torques largely average away for other stars. The distribution of angular-momentum changes δJ_ϕ induced by these torques per unit time is characterised by its first two moments ⟨δJ_ϕ⟩ and ⟨(δJ_ϕ)²⟩. The second moment determines the local diffusion of J_ϕ, which causes advective angular-momentum transport.

The first moment, ⟨δJ_ϕ⟩, describes non-local angular-momentum transport and is caused by the average torque at fixed J_ϕ. While the azimuthal average of −∂Φ/∂ϕ vanishes, that of the torque density τ ≡ −Σ(∂Ψ/∂ϕ) does not in general vanish because the distribution of stars is not axially symmetric owing to the spiral pattern itself. For potential and density perturbations of the form (6), the azimuthal average of τ is (Zhang & Buta 2007)

$$\begin{array}{c}\hfill {\langle \tau \rangle }_{\varphi }(R)=\frac{1}{2}m{S}_m(R)P_m(R)sin(m\delta \psi (R)),\delta \psi \equiv {\psi }_{\mathrm{\Psi }}-{\psi }_{\mathrm{\Sigma }}\end{array}$$(25)

if m_Ψ = m_Σ = m and ⟨τ⟩_ϕ = 0 otherwise. Thus, if the potential and density are out of phase at radius R, a spiral perturbation exerts a net torque over the annulus. For a trailing spiral, if σ δψ < 0 the potential leads the density, and stars gain angular momentum on average. Conversely, if σ δψ > 0, the potential trails the density, and stars lose angular momentum on average.

According to Noether’s theorem, the total angular momentum is conserved. The total torque over the whole disc, T  =  2π∫₀^∞⟨τ⟩_ϕR dR, therefore vanishes (as shown explicitly in Appendix E), and these gains and losses at different radii balance, but constitute a non-local transport of angular momentum from regions where σδψ > 0 to those where σδψ < 0.

The azimuthally averaged torque density ⟨τ⟩_ϕ can be computed from the numerically computed potential for any spiral density perturbation, but a simple way to approximate it would be insightful.

4.1. Torque in the tight-winding approximation

For a density component of the form (6a), the approximation (9) gives for the azimuthally averaged torque density

$$\begin{array}{cc}\hfill {\langle \tau \rangle }_{\varphi }& \approx -m\pi G{S}_m^2\mathfrak{I}\{k^{-1}\}.\hfill \end{array}$$(26)

The first-order approximations (10) have real-valued k(R) and hence give ⟨τ⟩_ϕ = 0. This is not surprising: As the product of the first-order quantities P_​m and sin δψ, the torque is a second-order quantity and cannot be predicted by first-order approximations. For the novel approximation (24), we find

$$\begin{array}{cc}\hfill {\langle \tau \rangle }_{\varphi }& \approx -\frac{\sigma \pi }{2m}GR{S}_m^{2}sin\alpha sin2\alpha (\gamma -\frac{3}{2}).\hfill \end{array}$$(27)

4.2. Torque from the potential

Instead of estimating ⟨τ⟩_ϕ from the target density and its second-order approximate potential as above, I now compute it for a first-order approximate potential and its (unknown) actual density using the exact expression (Binney & Tremaine 2008, Eq. 6.18)

$$\begin{array}{c}\hfill T(R)=-\frac{R}{4\pi G}{\displaystyle {\int }_0^{2\pi }}\mathrm{d}\varphi {\displaystyle {\int }_{-\infty }^{\infty }}\mathrm{d}z\frac{\partial \mathrm{\Phi }}{\partial R}{\frac{\partial \mathrm{\Phi }}{\partial \varphi }|}_R\end{array}$$(28)

for the cumulative disc torque T(R)≡2π∫₀^RdR′R′⟨τ⟩_ϕ(R′). In order to apply this, we must evaluate the integral over z. For the first-order tight-winding approximation based on plane waves, the natural approximation is Φ(r) = Ψ(R)e^−k(R)|z|, when the integral over z gives a factor 1/k(R), and we obtain after inserting Eqs. (6a) and (9), taking the derivatives, integrating over ϕ (and recalling that only the real part of Ψ has meaning)

$$\begin{array}{cc}\hfill T(R)& \approx \frac{\sigma \pi G{m}^2}{tan\alpha }\frac{S_m^2(R)}{k^3(R)}.\hfill \end{array}$$(29)

The azimuthally averaged torque density can now be estimated as ⟨τ⟩_ϕ = (dT/dR)/(2πR), which for the approximation (10b) again gives Eq. (27). As T(R) in Eq. (29) vanishes for R → ∞ (provided S_​m decays faster than R^−3/2 in this limit), this also implies that the estimate (27) obtains zero total torque exactly.

5. Application: The potential of spirals

I now consider various models for spiral perturbations and study the associated gravitational potentials, but also assess the ability of the tight-winding approximations to predict them.

5.1. Spirals with an exponential amplitude profile

The unperturbed (m = 0) surface density of disc galaxies generally follows an exponential decline,

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_0(R)={\mathrm{\Sigma }}_{\mathrm{c}}{\mathrm{e}}^{-R/R_{\mathrm{d}}}.\end{array}$$(30)

A simple and natural model for a spiral is one with a constant relative amplitude A_m ≡ |Σ_m|/Σ₀ (Section 5.2.3 considers A_m to depend on radius). Logarithmic spirals with this amplitude have the density

$$\begin{array}{c}\hfill \mathrm{\Sigma }(\mathit{R})={\mathrm{\Sigma }}_{\mathrm{sp}}{\mathrm{e}}^{-R/R_{\mathrm{d}}+\mathrm{i}m[\varphi -\lambda ln(R/R_0)]},\end{array}$$(31)

where Σ_sp = A_mΣ_c. I now consider various tight-winding approximations for the gravitational potential as well as that computed numerically.

5.1.1. Assessing the tight-winding approximations

A common model is

$$\begin{array}{c}\hfill \mathrm{\Psi }(\mathit{R})=-A_{\mathrm{\Psi }}R{\mathrm{e}}^{-R/R_{\mathrm{d}}+\mathrm{i}m[\varphi -\lambda ln(R/R_0)]},\end{array}$$(32)

with some (often unspecified) constant A_Ψ, which corresponds to the first-order approximation (9) with Eq. (10a) if

$$\begin{array}{c}\hfill A_{\mathrm{\Psi }}=2\pi G{\mathrm{\Sigma }}_{\mathrm{sp}}{m}^{-1}tan\alpha \end{array}$$(33)

(Contopoulos & Grosbøl 1986, and several subsequent studies) or Eq. (10b) if A_Ψ = 2πGΣ_spm⁻¹sin α. For constant and real-valued A_Ψ, the potential has the same phase as the target density (31), such that α_Ψ = α. My novel approximation can be expressed in this form by setting A_Ψ to a complex-valued function of radius. For five values of α, I compute numerically (via Eq. 5 and the gravity solver of Appendix B) the surface densities Σ(R) corresponding to these potential approximations and compare in Fig. 3 their properties to those of the target (31).

Fig. 3. Assessing the ability of the first-order (left) and our novel (right) tight-winding approximations to provide a gravitational potential for the target density (31) of a logarithmic spiral with pitch angle α and exponentially declining amplitude. From the density that is actually generated by the approximate potential, I plot the radial profiles of the relative amplitude error (top) and the errors in phase (middle) and pitch angle (bottom), which for the first-order methods are also the respective offsets between the potential and density because to first order, ψΨ = ψtarget.

The first-order method using the non-standard k(R) in Eq. (10b) approximates the amplitude significantly better than the standard form (10a). However, the two first-order methods give the same phase for the density as for the potential, and hence, an incorrect phase offset and pitch angle for the density, while the novel approximation gives an excellent match for R ≲ 3R_d or α ≲ 20 and still a reasonably accurate density phase at higher values. The full version of my novel approximation (Eq. (21)) is somewhat better than its reduced form (Eqs. (23) or (24)), but the difference is only significant in the regime where neither is very good any more.

The failure of all the tight-winding approximations at larger radii is not surprising because owing to the exponentially declining amplitude, the local field (on which these approximations are based) becomes ever less relevant compared to the far field from the inner parts (ignored by these approximations).

5.1.2. The exact gravitational potential

Finally, I also numerically computed the actual potential generated by the density (31) and plot in Fig. 4 its amplitude relative to that of the traditional first-order tight-winding approximation, as well as the phase offset and the deviation of the pitch angle between the density and potential. While these relations are apparently not easy to predict by tight-winding approximations for large radii and/or pitch angles, they are sufficiently smooth for interpolation, allowing a cost-efficient numerical implementation.

Fig. 4. Properties of a spiral with density (Eq. (31)) and its (numerically computed) potential. Radial profiles of the amplitude ratio to the first-order tight-winding approximation (Eqs. (9) and (10b)), the offsets of phase and pitch angle between the potential and density, and the torque per annulus (< 0 if angular momentum is lost for stars in a galaxy with a trailing spiral).

Again, the pitch of the potential is larger than that of the density and increases with radius at large α_Σ and R even more so than predicted by our novel second-order approximation.

In the bottom panel Fig. 4, I plot the contribution −σdT/dR of each annulus to the total torque. Owing to the increasing phase offset δψ between the potential and density changing sign at $R=\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$2$}\right.R_{\mathrm{d}}$, the torque is negative (for a trailing spiral) at $R<\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$2$}\right.R_{\mathrm{d}}$ and positive at $R>\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$2$}\right.R_{\mathrm{d}}$ (such that total angular momentum remains conserved).

The approximation (27) only slightly underestimates the local torque (dashed in the bottom panel of Fig. 4) and makes a much better impression than the approximations for the potential on which it is based. Moreover, the differences between Eq. (27) and its parent (Eq. (26) with k from Eq. (21)) are very small (not shown).

5.2. Radially tapered spirals

All the spirals considered so far extend radially over the whole disc, but real spirals may be limited to a radial range or, equivalently, an azimuthal range. Theoretical analysis of orbits perturbed by a spiral pattern suggest that the response of the system is in phase with the spiral only in a region around the co-rotation resonance (e.g. Contopoulos & Grosbøl 1986), and N-body simulations also show spiral patterns to be confined to a region around co-rotation and limited by the Lindblad resonances (Sellwood & Carlberg 2019). Moreover, in barred galaxies, spiral structures are mostly absent from the bar region.

5.2.1. The phase factor of Roberts et al.

With the goal to model a spiral that begins just outside a bar, Roberts et al. (1979) applied the phase

$$\begin{array}{c}\hfill {\psi }_{\mathrm{R}79}(R)\equiv \frac{\sigma }{Ntan{\alpha }_{\infty }}ln[1+{(R/R_{\mathrm{sp}})}^N]\end{array}$$(34)

to the potential

$$\begin{array}{c}\hfill \mathrm{\Psi }(\mathit{R})=-v_0^2\frac{a^{m+1}{R}^m}{{(a^2+R^2)}^{m+1/2}}{\mathrm{e}}^{\mathrm{i}m(\varphi -\psi )},\end{array}$$(35a)

with m = 2. At R ≪ R_sp, this phase function evaluates to zero and the potential perturbation stays aligned (α = 90°) as for a bar, but at R > R_sp, the potential perturbation becomes spiral-like, reaching a pitch angle of α → α_∞ at R ≫ R_sp with the sharpness of the transition determined by the parameter N. In case of a constant phase ψ, (35a) is the potential of a razor-thin perturbation with density²

$$\begin{array}{c}\hfill \mathrm{\Sigma }(\mathit{R})=\frac{v_0^2}{2\pi G}\frac{(2m+1)a^{m+2}{R}^m}{{(a^2+R^2)}^{m+3/2}}{\mathrm{e}}^{\mathrm{i}m(\varphi -\psi )},\end{array}$$(35b)

but for the phase (34) no formula for the associated density exists (and Roberts et al. 1979 did not attempt to estimate it). I computed it numerically (for the same parameter values as used by Roberts et al.) and plot in Fig. 5 its amplitude, phase, and pitch angle. The latter switches indeed from α = 90° at small R to α = α_∞ at large R as intended, albeit perhaps not as cleanly as for the potential: The bar region (R < R_sp = 1.8s indicated by grey lines in Fig. 5) is significantly pitched, and most of the torque exchange occurs within this region.

Fig. 5. Properties of the (numerically computed) density that generates the potential (35a) with the phase ψ = ψR79(R) (Eq. 34) for Rsp = 1.8a (grey vertical line) and N = 5 (as used by Roberts et al. 1979) and various values for α∞ (Roberts et al. used α∞ = 20°). For α∞ = 90°, ψ = 0, and the density is given by Eq. (35b).

The density declines like Σ ∝ Ψ_m/R ∝ R⁻⁴, as expected from the tight-winding approximation (9), and less steeply than for the bar-like case (α_∞ = 90°), i.e. Σ ∝ R⁻⁵ (Eq. 35b). The unperturbed (m = 0) component of the model used by Roberts et al. was also of type (35), such that the relative amplitude of the spiral declined like R⁻⁴/R⁻³ = 1/R.

At least two recent studies applied the phase (34) with N = 100 to the potential (32), with the intention of having a spiral pattern only at R > R_sp, but unaware that at R < R_sp this generates a bar-like perturbation (in addition to any other explicitly modelled bar) with significant σψ′_Σ < 0 (twisting the wrong way) near R = R_sp. I abstain from plotting these nonsensical models.

5.2.2. The spirals of Hamilton et al.

Hamilton et al. (2024) used a model with the potential being a logarithmic spiral,

$$\begin{array}{cc}\hfill \mathrm{\Psi }(\mathit{R})& =\frac{\eta }{m}{v}_0^2{B}_{\mathrm{H}24}(R,\beta ){\mathrm{e}}^{\mathrm{i}m[\varphi -\lambda ln(R/R_0)]}\hfill \end{array}$$(36a)

with

$$\begin{array}{cc}\hfill B_{\mathrm{H}24}(R,\beta )& \equiv exp(-\frac{{(R-R_{\mathrm{cr}})}^2}{2R_{\beta }^2}),\hfill \end{array}$$(36b)

where $R_{\beta }=\beta \sqrt{2}{R}_{\mathrm{cr}}/m$ with R_cr the co-rotation radius of the spiral pattern. Hamilton et al. modelled the m = 0 component of the total potential as Ψ = v₀²lnR, such that η is the (maximum of the) ratio of the tangential force to the unperturbed force. Moreover, for this model the epicycle and circular frequencies are related as $\kappa =\sqrt{2}\mathrm{\Omega }$, such that the Lindblad resonances are located at $R=(1\pm \sqrt{2}/m)R_{\mathrm{cr}}$, and for β = 1 the standard deviation of the potential envelope B_H24 equals the distance to the Lindblad resonances. Hamilton et al. considered $\beta =\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$, β = 1, and β = ∞. The latter case is identical to the scale-invariant spirals (Eq. (16)) with γ = 1, while for finite β, I computed Σ(R) numerically.

A common measure of the spiral strength is the ratio A_m = |Σ_m/Σ₀| of the density amplitudes of the perturbation and the unperturbed (m = 0) disc. I obtained (a lower limit for) A_m by assuming a maximum disc for the m = 0 component, i.e. Eq. (30) with Σ_c = 0.411 v₀²/GR_d, such that the maximum circular speed generated by the disc alone (without the dark-matter halo) equals v₀.

In Fig. 6, the resulting A_m = 2 is plotted (in the second panel from top) for the parameter choices used by Hamilton et al. (2024) and R_cr = 1.5R_d, while the lower panels show the phase offset and density pitch angle as in previous figures. For β = ∞, which is their favourite model, A_m = 2 has an inverted profile with a minimum at R = R_d. Towards smaller and larger radii, it increases and eventually exceeds unity, implying negative density. The finite-β models do not suffer from this issue (except at R ≪ R_d), but their spiral strength A_m = 2 is not symmetric with respect to R = R_cr (unlike the potential amplitude) and peaks at R > R_cr.

Fig. 6. Assessing the model (36) used by Hamilton et al. (2024) for the same parameters as used by these authors. Plotted as function of radius from top to bottom, I show the relative perturbation strength of the forces, the relative amplitude Am = 2 of the density perturbation (assuming the unperturbed disc is maximal with scale radius related to the co-rotation radius of the spiral as Rcr = 1.5Rd), the phase offset between the potential and density, and the density pitch angle.

For the model $\beta =\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$ with α = 30°, the density pitch angle decreases at R > R_cr, implying that the spiral becomes tighter than anticipated by Hamilton et al. (2024). The implications of all these properties for the results of the study by Hamilton et al. (2024) are discussed in Section 6.5.

5.2.3. Spirals with simple density models

All the models considered so far (with the exception of Fig. 4) were specified via the spiral potential, which is required for dynamical modelling. Unfortunately, as we have seen, simple potentials tend to correspond to spiral density patterns that deviate from the target, in particular, if the spiral pattern is restricted in radius. Instead, I now consider simple models in the density and numerically compute the corresponding potential.

I still assumed that the unperturbed (m = 0) density is of the form (30) and that the relative amplitude of the spiral perturbation is A_m(R) = |Σ_m|/Σ₀. The spiral perturbation then has the density

$$\begin{array}{c}\hfill \mathrm{\Sigma }(\mathit{R})={\mathrm{\Sigma }}_{\mathrm{c}}{A}_m(R){\mathrm{e}}^{-R/R_{\mathrm{d}}+\mathrm{i}m[\varphi -\lambda ln(R/R_0)]}.\end{array}$$(37a)

For constant A_m, this is identical to the target density (31) from the previous section, whose potential was presented in Fig. 4. Instead, I now consider simple non-constant amplitude functions of the form

$$\begin{array}{c}\hfill A_m(R)=A_{max}{B}_{\mathrm{H}24}(R,\beta )\end{array}$$(37b)

and plot in Fig. 7 for $\beta =\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$ and 1 (β = ∞ gives the model shown in Fig. 4) the same quantities as for the models of Hamilton et al. in Fig. 6, plus the torque per annulus.

Fig. 7. Similar to Fig. 6, but for the density models of Eq. (37a) with constant pitch angle α and relative perturbation density-amplitude Am(R) = AmaxBH24(R, β), i.e., the same profile as the relative perturbation force-amplitude for the model of Hamilton et al. (2024, shown in Fig. 6) for $\beta =\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$ (left) and β = 1 (right). Differently from that figure, the pitch angle of the potential is plotted. The estimates for the potential from the second-order tight-winding approximation (24) are also shown, as is the torque per annulus (bottom) and its approximation (27).

A comparison of these two figures shows many similarities: (i) The relative force perturbation peaks at smaller radii than the relative density perturbation, although this effect is weaker for our more realistic density models, where it is hardly present for $\beta =\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$. (ii) At radii where A_m is significant (not much smaller than its maximum), the phase offset δψ ≡ ψ_Ψ − ψ_Σ decreases slowly with radius (already seen in Fig. 4), and the potential pattern becomes more open with an outward-increasing pitch angle (although this effect is much weaker for the density models). Finally, (iii) at large radii where A_m is insignificant, the potential decays like R^{−1 − m} with α_Ψ = 90°. This is simply the expected behaviour for the outer field of the inner spiral (e.g. Binney & Tremaine 2008, Eq. (2.93)) without any significant local perturbation, and it is therefore also impossible to predict via the tight-winding approximations (dashed in Fig. 7). In the transition regions, at the edges of the spiral, the potential perturbation even twists backwards, as indicated by the pitch angle of the potential exceeding 90°.

In the bottom panels of Fig. 7, I also plot the net torque per annulus. The general pattern, already seen in Fig. 4 for a spiral without radial truncation, of outward angular-momentum transport by a trailing spiral is also present for these radially tapered spirals. The approximation (27) (dashed) is not perfect, but gives a reasonable estimate for this torque, in particular, for spirals with a larger radial extent and/or smaller pitch angles. Interestingly, the total torque at −σ⟨τ⟩ < 0, i.e. the total amount of angular momentum that is transported non-locally per unit time, is comparable for the models of $\beta =\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$ and β = 1 (and the same α): For the shorter spiral pattern, the torque is limited to a smaller region, which is compensated for by a higher torque density.

6. Discussion

6.1. How good is the tight-winding (WKB) approximation?

Many of the early results on disc stability rely on the first-order tight-winding (or an equivalent) approximation (e.g. Toomre 1964; Goldreich & Lynden-Bell 1965; Julian & Toomre 1966; Toomre 1981), and various later studies have employed models for the potential of spirals based on this approximation (e.g. Contopoulos & Grosbøl 1986; Grosbøl 1993; Grosbøl & Carraro 2018; Eilers et al. 2020). Little is known about its accuracy, however. Binney & Tremaine (2008) claimed that “in many situations [it] works fairly well even for |kR| as small as unity” (i.e. pitch angles as large as 60° for m = 2), but failed to give evidence of this claim. When prompted, Scott Tremaine (private communication) pointed to the fact that it works very well for Kalnajs’ (1971) spiral models, as I have shown in Section 3.1. These models are exceptional, however, in that the density and potential have the same phase, which the (first-order) tight-winding approximation always predicts, but which does not hold in general, nor for realistic spirals.

I assessed the tight-winding approximation for the important case of a spiral perturbation with an amplitude that declines exponentially, like the unperturbed surface density itself, such that the relative strength of the spiral is constant with radius. For this case, the amplitude of the perturbing potential is fairly well approximated, in particular, if the non-standard form (10b) for k is used. However, the first-order approximation does not give a phase offset δψ between the potential and density, and hence, it fails to predict the radially increasing pitch of the potential. If the spiral perturbation has a finite radial extent, the tight-winding approximation only works within the spiral region and fails at its edges and outside.

6.2. An improved tight-winding approximation

Kalnajs (1971) provided the exact potential of logarithmic spirals whose density amplitude declines like S_​m ∝ R^−3/2. I generalised these models to S_​m ∝ R^−γ, where the exponent γ is a free parameter. These appear to be the only flat spiral models whose potential is known in closed form. For $\gamma \ne \raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$, the phase of the potential in the equatorial plane, Ψ(R), differs from that of the density by a constant offset δψ = ψ_Ψ − ψ_Σ (such that the pitch angle is still the same), which is approximately proportional to $\gamma -\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$.

By locally approximating an arbitrary spiral pattern by such a scale-invariant model, I obtained the second-order tight-winding approximation (24)

$$\begin{array}{cc}\hfill \mathrm{\Psi }(\mathit{R})& \approx -2\pi G\frac{sin\alpha }{m}[1+\mathrm{i}\sigma \frac{sin2\alpha }{2m}(\gamma -\frac{3}{2})]R\mathrm{\Sigma }(\mathit{R}),\hfill \end{array}$$

where α = α(R) is the local pitch angle of the density and γ = γ(R) =  − dlnS_​m/dlnR the power-law slope of its amplitude, while σ = sign(dψ_Σ/dR) gives the sense of the spiral twist (for a trailing spiral $\sigma =-\mathit{sign}({\overline{v}}_{\varphi })$). This novel tight-winding approximation is much better than the traditional first-order approach. In particular, it predicts the phase offset δψ and pitch of the potential with reasonable accuracy.

6.3. The potential of typical spiral perturbations

As the surface density of spiral galaxies typically declines nearly exponentially with radius, so does the amplitude of its spiral perturbation (except at its inner and outer edges). In particular, the amplitude declines faster than a power law, and as a consequence, the phase offset δψ is no longer constant. We may estimate it from the second-order approximation (24) as

$$\begin{array}{cc}\hfill \delta \psi (R)& \approx -\frac{\sigma }{m}arctan\left[\frac{sin2\alpha }{2m}(\gamma -\frac{3}{2})\right]\simeq -\frac{\sigma }{m^2}\alpha (\gamma -\frac{3}{2}).\hfill \end{array}$$(38)

Typically, the density amplitude is shallower than S_​m ∝ R^−3/2 at small radii and steeper at larger radii, such that the potential trails the density at small R and leads it at large R (for a trailing spiral). For a purely exponential profile with scale radius R_d, the transition occurs at $R=\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$2$}\right.R_{\mathrm{d}}$.

I demonstrate this in Fig. 8, which plots the crest lines (the positions of the density maxima and potential minima at each radius) for logarithmic spirals with an exponentially declining amplitude and various pitch angles (for the density). The solid red curves show the exact potential, and the dashed curves were obtained by the approximation (24), which out to R = 5R_d predicts the actual potential very well. Beyond this radius, this local approximation fails because the potential depends ever more on the inner rather than the local spiral.

Fig. 8. Crest lines of the density Σ (blue) and potential Ψ (red, partly concealed by those of the density) for logarithmic m = 2 spirals with a pitch angle α as indicated (top) and an exponentially declining radial density-amplitude (i.e. constant relative amplitude Am = |Σm|/Σ0), either radially unlimited (top) or with sharp inner and outer edges (bottom). I also plot (dashed) crest lines of the potential obtained via the second-order approximation (24) (to first order, the crest lines of Σ and Ψ are identical).

The bottom panels of Fig. 8 plot the same, but for a spiral that has a non-zero (and exponentially declining) amplitude only at R_d < R < 4R_d. In most of this range, the potential closely follows the density, as predicted by the tight-winding approximations, but towards the edges and outside of the spiral range, the potential aligns itself to a constant phase and decays like Ψ ∝ R^m at R → 0 and Ψ ∝ R^{−m − 1} at R → ∞, as expected for the outer potential of an m-fold perturbation.

Galactic spirals tend to have an outward-decreasing pitch angle (Savchenko & Reshetnikov 2013), which implies via Eq. (38) that the phase offset does not increase as much as for α =  const.

6.4. Net torque and non-local angular momentum transport

The potential of a spiral perturbation exerts a torque −∂Φ/∂ϕ per unit mass. For most stars, this torque averages out, but for stars in orbital resonance, the torque accumulates and causes significant angular-momentum changes δJ_ϕ. These changes cause radial migration (Sellwood & Binney 2002) by (i) advective angular-momentum transport caused by local diffusion of J_ϕ, as characterised by ⟨(δJ_ϕ)²⟩, and (ii) non-local angular-momentum transport described directly by ⟨δJ_ϕ⟩.

The sign of δJ_ϕ depends on the relative orbital phase of the star to the spiral perturbation, which is not completely random, as the spiral pattern is made up of the very same stars. As a result, there remains a net torque at each radius, which is second order in both the spiral amplitude and pitch angle, and the sign of which depends on whether the potential of the perturbation trails (σδψ > 0) or leads (σδψ < 0) the density. For trailing spirals, $\sigma =-\mathit{sign}({\overline{v}}_{\varphi })$ and the net torque reduces stellar angular momenta at small radii and increases them at large radii, constituting a non-local outward transport of angular momentum by gravitational torques (Zhang 1996).

Thus, every trailing spiral perturbation transports angular momentum outwards, not only via local diffusion, but also by a net gravitational torque at each radius. Remarkably, this conclusion can be derived analytically (in Sections 4.1 and 4.2) without considering orbital dynamics. It implies that spirals alter the state of the disc irreversibly.

6.5. Implications of poor spiral models

The results of models using inadequate gravitational potentials for the spiral pattern are obviously compromised. Unfortunately, these shortcomings are not as uncommon as might be hoped. I give two recent examples.

With the intention to model Milky Way spirals, Eilers et al. (2020) used a gravitational potential that was meant to be the first-order tight-winding approximation (32), but deviated by a factor Rtanα/mh_z with h_z = 1 kpc. This, however, gives the tight-winding approximation for a density amplitude that differs from the intended exponential model by the same factor, i.e. the relative amplitude A_m = |Σ_m|/Σ₀ of the spiral perturbation increases linearly with radius. This is certainly unphysical and implies that the associated dynamical model of the non-axisymmetric stellar kinematics is inadequate Eilers et al. (2020) also neglected the gravity from the Milky Way bar, which renders their model even less adequate).

As already shown in Section 5.2.2, the models employed by Hamilton et al. (2024) also used unrealistic profiles of the density amplitude of the spiral perturbation, in particular, their model with β = ∞ in Eqs. (36), which Hamilton et al. favoured because it would resemble the model of Eilers et al. (2020). Their models $\beta =\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$ and β = 1 have more reasonable density profiles and induce much less radial heating, which agrees well with the situation that is observed for the Milky Way. The claim of Hamilton et al. that the Milky Way stellar disc is dynamically colder than expected is therefore unfounded.

6.6. The vertical dimension

In this study, only flat (razor-thin) spiral perturbations have been considered, and of these, only their gravitational potential in the equatorial plane. This is often sufficient for modelling horizontal disc dynamics, at least to first order.

While the discs and their spiral perturbation of real galaxies are thin in the sense that the scale height h_z is much smaller than the scale length R_d, h_z can become significant compared to the radial wavelength of a spiral perturbation. In the first-order tight-winding approximation, the potential decays as e^−k|z| away from z = 0, where k = m/Rsin α, and convolving this with the vertical profile ∝e^−|z|/h_z reduces the central potential by the factor (1 + h_zk)⁻¹. This reduction is larger at smaller R and smaller α (more tightly wound spirals), and implies that razor-thin models overestimate the gravity of spiral perturbations by ∼ (1 + h_zk).

To avoid this and model the vertical extent of galactic discs and their spiral perturbations more accurately, full three-dimensional models for these perturbations and their potentials are required. These are the subject of an ongoing study.

7. Conclusions

I studied the gravitational potential of flat (razor-thin) spiral perturbations by means of an accurate numerical solution and compared them to the widely used tight-winding (or WKB) approximation. I also extended Kalnajs’ (1971) scale-invariant spirals to general power-law amplitudes |Σ_m|∝R^−γ (Kalnajs’s spirals have $\gamma =\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$). I studied various simple models for spiral potentials or densities, most of which were logarithmic spirals (constant pitch angle α). My main results are as follows.

• For typical m = 2 spirals with exponentially declining amplitude (and hence, with constant relative amplitude), the traditional (first-order) tight-winding approximation gives ≲10% errors for pitch angles α ≲ 20°, and deteriorates beyond that.

• Models for spiral potentials used in the literature are typically based on the first-order tight-winding approximation and inherit all its weaknesses, implying that the actually implied density deviates from the intended target, sometimes significantly so.

• A second-order tight-winding approximation (Eq. 24), based on locally approximating the spiral perturbation by a scale-invariant spiral (rather than a plane wave), is significantly better and can predict the phase offset δψ of the potential from the density (unlike the first-order approximations, which always give δψ = 0).

• Towards the inner or outer edge of a spiral perturbations, all tight-winding approximations fail. Beyond the edges, the gravitational potential of the spiral resembles that of a weak bar. Hence, to model these radially limited spiral perturbations, approximate techniques cannot be used, but numerical treatment is required.

• Generally, the gravitational potential for trailing spirals trails the density at small radii and leads it at large radii, with the transition at the radius at which the density amplitude declines like R^−3/2. These phase offsets result in gravitational torques that transport angular momentum outwards and change the state of the disc irreversibly. Eq. (27) gives a simple but reasonably accurate approximation for the average torque density at each radius.

• I developed an efficient numerical method for computing the gravitational potential of spiral perturbations for arbitrary multiplicity m, phase functions ψ(R), and amplitudes S_​m(R). Implementations in C++ and python are made available.

Acknowledgments

I thank Howard Cohl for discussing the Laplace coefficients and for sharing Kalnajs’ letter to him, Scott Tremaine for discussing the computation of the Laplace coefficients and the tight-winding approximation, Agris Kalnajs for opening my eyes to gravitational torques by pointing out the related problem of the scale-invariant spirals with $\gamma \ne \raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$, and the anonymous reviewer for a thorough reading and a prompt, encouraging, and helpful report.

References

Appendix A: Proof of equation (5)

I begin with the observation that the potential of a razor-thin disc can always be expressed as Φ(r) = F(R, |z|), where F(r) is a smooth function satisfying ΔF = 0 for z ≥ 0 (one can construct F by setting F = Φ for z ≥ 0 and by analytic continuation for z < 0). Then Ψ(R) = F(R, 0) and Poisson’s equation at z = 0 implies

$$\begin{array}{cc}\hfill 2\pi G\mathrm{\Sigma }(\mathit{R})& ={\displaystyle \underset{z\to 0^{+}}{lim}}\frac{\partial \mathrm{\Phi }}{\partial z}=\frac{\partial F}{\partial z}{|}_{z=0},\text{and}\hfill \end{array}$$(A.1)

$$\begin{array}{cc}\hfill \mathrm{\Delta }\mathrm{\Psi }& =-\frac{\partial F}{\partial z}{|}_{z=0}.\hfill \end{array}$$(A.2)

Because $\stackrel{\sim }{F}(\mathit{r})\equiv -\partial F/\partial z$ also satisfies $\mathrm{\Delta }\stackrel{\sim }{F}=0$, another razor-thin model can be constructed, whose potential in the plane

$$\begin{array}{c}\hfill \stackrel{\sim }{\mathrm{\Psi }}(\mathit{R})=\stackrel{\sim }{F}(\mathit{R},0)=-\frac{\partial F}{\partial z}{|}_{z=0}=-2\pi G\mathrm{\Sigma }(\mathit{R}).\end{array}$$(A.3)

For the surface density of the new model, Eq. (A.1) applied to $\stackrel{\sim }{F}$ gives

$$\begin{array}{cc}\hfill 2\pi G\stackrel{\sim }{\mathrm{\Sigma }}(\mathit{R})& =\frac{\partial \stackrel{\sim }{F}}{\partial z}{|}_{z=0}=-\frac{\partial F}{\partial z}{|}_{z=0}=\mathrm{\Delta }\mathrm{\Psi },\hfill \end{array}$$(A.4)

where the last equality follows from Eq. (A.2). Finally, when the Poisson integral (3) relating $\stackrel{\sim }{\mathrm{\Psi }}$ and $\stackrel{\sim }{\mathrm{\Sigma }}$ is expressed in terms of Σ and ΔΨ using Eqs. (A.3, A.4), Eq. (5) follows. ▪

Appendix B: Evaluating the Poisson integral numerically

Assume that the density has been decomposed into azimuthal Fourier components

$$\begin{array}{cc}\hfill \mathrm{\Sigma }(\mathit{R})& ={\displaystyle \sum _{m=-\infty }^{\infty }}{\mathrm{\Sigma }}_m(R){\mathrm{e}}^{\mathrm{i}m\varphi }\text{with}{\mathrm{\Sigma }}_m(R)=\frac{1}{2\pi }{\displaystyle {\int }_0^{2\pi }}\mathrm{d}\varphi \mathrm{\Sigma }(\mathit{R}){\mathrm{e}}^{-\mathrm{i}m\varphi },\hfill \end{array}$$(B.1)

where the Σ_m(R) are generally complex-valued and often given in polar form Σ_m(R) = S_​m(R) e^−imψ(R) with real-valued amplitude S_​m and phase ψ. When this is inserted into the Poisson integral (3), the azimuthal Fourier component of the potential is found to be

$$\begin{array}{ccc}\hfill {\mathrm{\Psi }}_m(R)& =-G{\displaystyle {\int }_0^{\infty }}\mathrm{d}{R}^{\prime }{R}^{\prime }{\mathrm{\Sigma }}_m(R^{\prime }){\displaystyle {\int }_0^{2\pi }}\frac{\mathrm{d}{\varphi }^{\prime }{\mathrm{e}}^{\mathrm{i}m({\varphi }^{\prime }-\varphi )}}{\sqrt{R^2-2R{R}^{\prime }cos|\varphi -{\varphi }^{\prime }|+R^{\prime 2}}}\hfill & \hfill =-\pi G{\displaystyle {\int }_0^{\infty }}\mathrm{d}{R}^{\prime }{\mathrm{\Sigma }}_m(R^{\prime })b_{1/2}^m(R/R^{\prime }),\end{array}$$(B.2)

where

$$\begin{array}{c}\hfill b_s^m(\alpha )=\frac{1}{\pi }{\displaystyle {\int }_0^{2\pi }}\frac{{\mathrm{e}}^{\mathrm{i}m\phi }\mathrm{d}\phi }{{(1-2\alpha cos\phi +{\alpha }^2)}^s}=\frac{1}{2\pi }{\displaystyle {\int }_0^{\pi }}\frac{cosm\phi \mathrm{d}\phi }{{(1-2\alpha cos\phi +{\alpha }^2)}^s}\end{array}$$(B.3)

are the ‘Laplace coefficients’. Their computation is detailed in Appendix C, while Fig. B.1 plots them for $s=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$. At α = 1, b_1/2^m has a logarithmic singularity for all m. Other important properties are

$$\begin{array}{cccccccc}\hfill b_s^{-m}(\alpha )& =b_s^m(\alpha ),\hfill & \hfill b_s^m({\alpha }^{-1})& ={\alpha }^{2s}{b}_s^m(\alpha ),\hfill & b_s^m(\alpha )& \sim {\alpha }^m\text{as}\alpha \to 0,\text{and}& b_s^m(\alpha )& \sim {\alpha }^{-2s-m}\text{as}\alpha \to \infty .\end{array}$$(B.4)

In order to compute the derivative dΨ_m/dR needed for the force, the derivative of b_1/2^m in Eq. (B.2) could be taken via the relation (Hagihara 1972, eq. 7.4.13)

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{b}_s^m}{\mathrm{d}\alpha }& =\frac{[m+(m+2s){\alpha }^2]b_s^m(\alpha )-2[m+1-s]b_s^{m+1}(\alpha )}{\alpha (1-{\alpha }^2)}.\hfill \end{array}$$(B.5)

Alternatively, because Eq. (B.2) is a convolution, the derivative can be shifted onto Σ_m:

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{\mathrm{\Psi }}_m}{\mathrm{d}R}& =-\frac{\pi G}{R}{\displaystyle {\int }_0^{\infty }}\mathrm{d}{R}^{\prime }\frac{\mathrm{d}(R^{\prime }{\mathrm{\Sigma }}_m)}{\mathrm{d}{R}^{\prime }}{b}_{1/2}^m(R/R^{\prime }).\hfill \end{array}$$(B.6)

Fig. B.1. The Laplace coefficients b1/2m(α) (solid) and their asymptotes at α ≪ 1 and α ≫ 1 (dotted).

My aim is to employ Gauß-Legendre quadrature to evaluate the integrals in (B.2) and (B.6). This requires to transform to integrals with finite boundaries and finite integrand. To avoid the infinite integration interval and shift the singularity to x = 0, I use the substitution

$$\begin{array}{c}\hfill x=\frac{R^{\prime }-R}{R^{\prime }+R+a},R^{\prime }=\frac{R+x(R+a)}{1-x}\end{array}$$(B.7)

with some scale length a ≥ 0 (for small R it is often better to set a > 0), such that the integral (B.2) becomes

$$\begin{array}{cc}\hfill {\mathrm{\Psi }}_m(R)& =\pi G{\displaystyle {\int }_{-R/(R+a)}^1}\mathrm{d}x\frac{2R+a}{{(1-x)}^2}{\mathrm{\Sigma }}_m(R^{\prime }(x))b_{1/2}^m(R/R^{\prime }(x)).\hfill \end{array}$$(B.8)

The logarithmic singularity at x = 0 can now be dealt with by the second substitution x = u³ or x = u⁵, such that the integrand is smooth at u = 0. Alternatively, the singularity can be treated by a technique used by Huré (2005) in a very similar context. Suppose we already know a potential-density pair ${\overline{\mathrm{\Psi }}}_m(R)$, ${\overline{\mathrm{\Sigma }}}_m(R)$ that is an exact solution. Then

$$\begin{array}{c}\hfill {\mathrm{\Psi }}_m(R)={\mathrm{\Sigma }}_m(R)\frac{{\overline{\mathrm{\Psi }}}_m(R)}{{\overline{\mathrm{\Sigma }}}_m(R)}-\pi G{\displaystyle {\int }_0^{\infty }}\mathrm{d}{R}^{\prime }[{\mathrm{\Sigma }}_m(R^{\prime })-{\mathrm{\Sigma }}_m(R)\frac{{\overline{\mathrm{\Sigma }}}_m(R^{\prime })}{{\overline{\mathrm{\Sigma }}}_m(R)}]b_{1/2}^m(R/R^{\prime }).\end{array}$$(B.9a)

Thus, the integral only computes the difference between Ψ_m and $({\mathrm{\Sigma }}_m/{\overline{\mathrm{\Sigma }}}_m)$ times ${\overline{\mathrm{\Psi }}}_m$. The point of this construct is that the integrand now vanishes at R′=R and the previous method without the second substitution can be used. A suitable choice for ${\overline{\mathrm{\Psi }}}_m$ and ${\overline{\mathrm{\Sigma }}}_m$ is the model (35) with parameter a set such that ${\overline{\mathrm{\Sigma }}}_m(R\prime )$ is maximal at R′=R, giving

$$\begin{array}{cccc}\hfill \frac{{\overline{\mathrm{\Psi }}}_m(R)}{{\overline{\mathrm{\Sigma }}}_m(R)}& =-\frac{2\pi GR}{\sqrt{m(m+3)}}\frac{2m+3}{2m+1},\hfill & \hfill \frac{{\overline{\mathrm{\Sigma }}}_m(R^{\prime })}{{\overline{\mathrm{\Sigma }}}_m(R)}& =\frac{{(2m+3)}^{m+3/2}{(R^{\prime }/R)}^m}{{(m+3+m{[R^{\prime }/R]}^2)}^{m+3/2}}.\hfill \end{array}$$(B.9b)

Appendix C: Computing the Laplace coefficients

The Laplace coefficients b_s^m(α), defined in Eq. (B.3), have been studied by Laplace, Lagrange, Legendre, Euler, Cauchy, and Gauß amongst many others and are covered in every textbook on celestial mechanics (e.g. Hagihara 1972; Tremaine 2023). Because of their properties (B.4), their computation can be limited to m ≥ 0 and 0 ≤ α ≤ 1, and I assume that these conditions are satisfied. The Laplace coefficients are related to the associated Legendre functions of the second kind, Q_ν^μ, with half-integer degree ν and integer order μ (also known as toroidal functions) via

$$\begin{array}{cc}\hfill b_s^m(\alpha )& =\frac{2}{\mathrm{\Gamma }(s)\sqrt{\pi \alpha }{({\alpha }^2-1)}^{s-1/2}}{Q}_{m-1/2}^{s-1/2}(\chi )\text{for}\alpha <1\text{with}\chi \equiv \frac{1}{2}(\alpha +{\alpha }^{-1}).\hfill \end{array}$$(C.1)

A representation due to Lagrange (Hagihara 1972, eq. 7.1.2) is

$$\begin{array}{cc}\hfill b_s^m(\alpha )& =\frac{2\mathrm{\Gamma }(m+s)}{\mathrm{\Gamma }(s)\mathrm{\Gamma }(m+1)}{\alpha }^m{}_2{F}_1(s,m+s;m+1;{\alpha }^2)=\frac{2{\alpha }^m}{\mathrm{\Gamma }{(s)}^2}{\displaystyle \sum _{n=0}^{\infty }}\frac{\mathrm{\Gamma }(n+s)\mathrm{\Gamma }(m+n+s)}{\mathrm{\Gamma }(n+1)\mathrm{\Gamma }(m+n+1)}{\alpha }^{2n},\hfill \end{array}$$(C.2a)

where ₂F₁ is the hypergeometric function. One possibility to compute the Laplace coefficients is via the hypergeometric series (the second expression in Eq. C.2a), which gives the algorithm

$$\begin{array}{cccccc}\hfill b_s^m(\alpha )& ={\displaystyle \sum _{n=0}^{\infty }}{t}_{s,n}^m,\hfill & \hfill t_{s,0}^m& =2\frac{{(s)}_m}{m!}{\alpha }^m,\hfill & \frac{t_{s,n+1}^m}{t_{s,n}^m}& =\frac{(n+s)(m+n+s)}{(n+1)(m+n+1)}{\alpha }^2,\end{array}$$(C.2b)

where (a)_n = a(a + 1)…(a + n − 1) = Γ(a + n)/Γ(a) is the (rising) Pochhammer symbol. This algorithm converges for α < 1, but only very slowly for α → 1, where b_s^m(α) is singular. For $s=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$, the singularity at α = 1 is logarithmic, i.e. b_1/2^m(α)∼ln(1 − α), and a suitable alternative expression based on Eq. (15.3.10) of Abramowitz & Stegun (1972) is

$$\begin{array}{cc}\hfill b_{1/2}^m(\alpha )& =\frac{2{\alpha }^m}{{\pi }^{3/2}\mathrm{\Gamma }(m+\frac{1}{2})}{\displaystyle \sum _{n=0}^{\infty }}\frac{\mathrm{\Gamma }(n+\frac{1}{2})\mathrm{\Gamma }(m+n+\frac{1}{2})}{n{!}^2}[2\psi (n+1)-\psi (n+\frac{1}{2})-\psi (m+n+\frac{1}{2})-lnϵ]{ϵ}^n,\hfill \end{array}$$(C.3a)

where ϵ = 1 − α² and ψ(z)≡Γ′(z)/Γ(z) is the digamma function. The resulting algorithm is

$$\begin{array}{cccccc}\hfill b_{1/2}^m(\alpha )& ={\displaystyle \sum _{n=0}^{\infty }}{p}_n^m{q}_n^m,\hfill & \hfill p_0^m& =\frac{2{\alpha }^m}{\pi },\hfill & \frac{p_{n+1}^m}{p_n^m}& =\frac{(n+\frac{1}{2})(m+n+\frac{1}{2})}{n^2}{ϵ}^2,\\ \hfill & q_0^m\hfill & \hfill =4ln2-lnϵ-{\displaystyle \sum _{k=0}^{m-1}}\frac{1}{k+\frac{1}{2}},& q_{n+1}^m\hfill & =q_n^m+\frac{2}{n+1}-\frac{1}{n+\frac{1}{2}}-\frac{1}{m+n+\frac{1}{2}}.\end{array}$$(C.3b)

Fig. C.1. Relative error (top) and costs (bottom) of the computation of the Laplace coefficient bsm via four different methods described in the text for $s=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$ and m = 2 (left) or m = 4 (right).

Potentially more efficient algorithms are based on the recursion relations

$$\begin{array}{cc}\hfill (m+1-s)b_s^{m+1}(\alpha )& =m[\alpha +{\alpha }^{-1}]b_s^m(\alpha )-(m-1+s)b_s^{m-1}(\alpha )\hfill \end{array}$$(C.4a)

(Euler in Hagihara 1972, eq 7.7) to increment or decrement m and

$$\begin{array}{cc}\hfill s{(1\pm \alpha )}^2[b_{s+1}^m(\alpha )\mp b_{s+1}^{m+1}(\alpha )]& =(m+s)b_s^m(\alpha )\pm (m-s+1)b_s^{m+1}(\alpha )\hfill \end{array}$$(C.4b)

(Hagihara 1972, eq. 7.11) to increment or decrement s. From these and some starting values for b_s^m and b_s^m + 1, the values of the Laplace coefficients for any s and m can, in theory, be computed. Possible starting values are b_1/2⁰(α) = 4K(α)/π and b_1/2¹(α) = 4[K(α)−E(α)]/πα, where K and E are complete elliptical integrals of the first and second kind, respectively. These can be efficiently computed simultaneously from their relation to Gauß’ arithmetic-geometric mean (NIST 2024, §https://dlmf.nist.gov/19.8i), which gives the algorithm

$$\begin{array}{cccccc}\hfill a_0,g_0& =1,\sqrt{1-{\alpha }^2},\hfill & \hfill a_{n+1},g_{n+1}& =\frac{1}{2}(a_n+g_n),\sqrt{a_n{g}_n},\hfill & b_{1/2}^0(\alpha )& =2/a_{\infty },\end{array}$$(C.5a)

$$\begin{array}{cccccc}\hfill c_0& =\frac{1}{2}{\alpha }^2,\hfill & \hfill c_{n+1}& =c_n+2^{n-3}\sqrt{a_n-g_n},\hfill & b_{1/2}^1(\alpha )& =b_{1/2}^0(\alpha )c_{\infty }/\alpha ,\end{array}$$(C.5b)

which converges quadratically to a_∞ = g_∞.

Unfortunately, the recursion (C.4a) is unstable in the upwards direction³ and hence cannot be used to compute the Laplace coefficients for α^m ≪ 1 with sufficient accuracy when starting from m = 0, 1. In that regime, one may either use the series expression (algorithm C.2b) or apply the recursion (C.4a) downwards (Miller’s algorithm). To this end, one starts from some value m′≫m with initial values b_s^m′+1 = α and b_s^m′ = 1, corresponding to the asymptotic behaviour b_s^m + 1/b_s^m ∼ α at small α. Next, the recursion is applied to obtain corresponding values for b_s^m′−1,  b_s^m′−2, …, b_s⁰. At small m, these are only wrong by the same overall normalisation (for the same reason that the recursion is unstable in the upwards direction), which for $s=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$ can be obtained from b_1/2⁰(x) via (C.5a) and in general from the relation Tremaine (2023, problem 4.6)

$$\begin{array}{cc}\hfill \frac{1}{2}{b}_s^0(\alpha )+{\displaystyle \sum _{m=1}^{\infty }}{b}_s^m(\alpha )& ={(1-\alpha )}^{-2s}.\hfill \end{array}$$(C.6)

The asymptotic b_s^m(α)∼α^m for small α implies that reaching a precision ϵ requires α^m′−m = ϵ or m′=m + lnϵ/lnα = m − 16/log₁₀(α) for the usual 64-bit floating-point arithmetic. These recursive methods are particularly useful for computing several Laplace coefficients with different s or m, but even when only a single coefficient is required, as in our application, they are competitive.

This gives four different methods to compute the Laplace coefficients. For $s=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$, the only value actually required by the methods of Appendix B, I assessed them by comparison to a calculation based on the first equality in (C.2a) computed using implementations for the hypergeometric and Gamma functions from the Gnu Scientific Library (gsl, Gough 2011). At α → 1 and m ≫ 1, the gsl routines become inaccurate or even fail altogether (exit with an error), in which case I took the mean of the two least deviating values from my four methods as estimate for the correct value. Fig. C.1 plots the resulting relative errors (top) and the required CPU time (bottom) for b_1/2² (left) and b_1/2⁴ (right). It appears that for small α the direct summation of the hypergeometric series (blue) is accurate and most efficient, while at larger α the upward recursion (green) is the best method.

Appendix D: The potential of the scale-invariant spirals

Following the ansatz of Appendix A, the potential can be expressed as Φ(r) = F(R, |z|) with a smooth function F(r) that satisfies

$$\begin{array}{cc}\hfill \mathrm{\Delta }F& =0\text{at}z\ge 0,F\to 0\text{as}z\to \infty ,\text{and}\frac{\partial F}{\partial z}{|}_{z=0}=2\pi G\mathrm{\Sigma }(\mathit{R}).\hfill \end{array}$$(D.1)

Since Σ(R) = S_​m, 0(R/R₀)^{−γ − imλ}e^imϕ ∝ R^{−3/2 − iζ}e^imϕ is scale-invariant, a natural ansatz for F(r) is a scale-invariant form, too:

$$\begin{array}{cc}\hfill F(\mathit{r})& =-2\pi G{R}_0{S}_{m,0}{(r/R_0)}^{-1/2-\mathrm{i}\zeta }{f}_m(cos\theta ){\mathrm{e}}^{\mathrm{i}m\varphi },\hfill \end{array}$$(D.2)

with spherical coordinates (r, θ, ϕ). The function f_m must be determined by the conditions (D.1), in particular (with x = cos θ)

$$\begin{array}{c}\hfill \mathrm{\Delta }F\propto {(r/R_0)}^{-5/2-\mathrm{i}\zeta }{\mathrm{e}}^{\mathrm{i}m\varphi }\{[(-\frac{1}{2}-\mathrm{i}\zeta )(\frac{1}{2}-\mathrm{i}\zeta )-\frac{m^2}{1-x^2}]f_m(x)+\frac{\mathrm{d}}{\mathrm{d}x}(1-x^2)\frac{\mathrm{d}{f}_m}{\mathrm{d}x}\}=0.\end{array}$$(D.3)

For this to hold and f → 0 as x → 1 (so that F → 0 as z → ∞), $f_m(x)=C{P}_{-1/2-\text{i}\zeta }^{|m|}(x)$, where P_ν^μ is the associated Legendre function of the first kind. The constant C is determined by the last of the conditions (D.1), giving $C^{-1}=-\text{d}{P}_{-1/2-\text{i}\zeta }^{|m|}/\text{d}x{|}_{x=0}$. Thus,