A new approach to the problem of modes in mestel disks [original]

INAUGURAL-DISSERTATION

zur

Erlangung der Doktorw urde

der

Naturwissenshaftlih-Mathematishen Gesamtfakultat

der

Rupreht-Karls-Universitat Heidelb erg

vorgelegt von

Dipl.-Phys. Markus Demleitner

aus Erlangen

Tag der m undlihen Pr ufung: 17. Mai 2000

Uber einen neuen Zugang

zum Problem von

Moden in Mestelscheiben

Gutahter: Prof. Dr. Burkhard Fuhs

Prof. Dr. Josef Fried

Dissertation

Submitted to the

Combined Faulties for the Sienes and Mathematis

of Rup erto Carola University, Heidelb erg, Germany

for the Degree of

Do tor of Sienes

A New Approach to

the Problem of Modes

in Mestel Disks

presented by

Markus Demleitner, b orn in Erlangen

Heidelb erg, May 17, 2000

Gutahter: Prof. Dr. Burkhard Fuhs

Prof. Dr. Josef Fried



Ub er einen neuen Zugang zum Problem von Mo den in Mestelsheib en

Diese Arbeit untersucht, wann Scheiben mit global flacher Rotationskurve (Mestelscheiben) sta-

tion¨are, harmonische St¨orungen (Normalmoden) zulassen und wie diese beschaffen sind. W¨ahrend

die existierenden Analysen dieses Problems (Zang 1976, Read 1997) die Orbits in diesen Scheiben

seminumerisch bestimmen, werden hier von vorneherein gen¨aherte Orbits verwendet. Dies er-

laubt eine Bestimmung des Kerns der das Problem beschreibenden Integralgleichung in fast

geschlossener Form.

Zus¨atzlich zur selbstkonsistenten Scheibe wird der Fall einer Scheibe, in der ein Teil der

Matrie immobilisiert wird (cut-out disk, gekappte Scheibe), untersucht. Die gekappte Scheibe

hat eine finite aktive Masse und ist nicht mehr selbst¨ahnlich, was zu einem qualitativ anderen

Verhalten f¨uhrt. Im selbstkonsistenten Fall ist der Kern relativ ¨uberschaubar, wenn auch immer

noch zu kompliziert, um auf analytische L¨osungen hoffen zu k¨onnen. Die Terme f¨ur die gekappte

Scheibe hingegen sind ausgesprochen unhandlich.

Im Wesentlichen reproduziert die N¨aherung die Ergebnisse der vorangehenden Analysen be-

merkenswert gut. Im Fall der selbstkonsistenten Scheibe erm¨oglicht der geschlossene Ausdruck des

Kerns ein verbessertes Verst¨andnis des Eigenschaften der neutralen (nichtrotierenden, nichtwach-

senden) Moden. Zus¨atzlich bietet er einen Zugang zu der nach den Arbeiten von Zang und Read

noch offenen Frage der nichtaxialsymmetrischen wachsenden Moden. Solche Moden k¨onnen nicht

¨uber der Stabilit¨atsgrenze der neutralen Moden existieren, und unterhalb dieser Grenze k¨onnen

sie keine beschr¨ankte Mellintransformierte besitzen. Die Existenz solcher Moden konnte nicht

gezeigt werden.

A New Approah to the Mo des in Mestel Disks

In this work I examine the modes admitted by the Mestel disk, a disk with a globally flat rotation

curve. In contrast to previous analyses of this problem by Zang (1976) and Read (1997), I

approximate the orbits to obtain almost closed expressions for the kernel of the integral equation

governing the behaviour of the modes.

I investigate the modes admitted by both the self-consistent and a cut-out Mestel disk, the

difference being that in the latter case a part of the matter in the disk is immobilized. This breaks

the self-similarity and produces a pronouncedly different picture. While the expressions for the

kernel in the self-consistent disk are quite managable (though still beyond the reach of analytic

techniques), the kernels for cut-out disks tends to rather complicated indeed.

In general, my approximation reproduces the results of the previous works remarkably well.

Due to the sheer size of the terms, examining the solution behaviour in the approximation does not

save computing time compared to Zang’s method at least for the cut-out disks. The more handy

expressions in the self-consistent disk, on the other hand, allow an intuitive understanding of

most of the properties of neutral (nonrotating, nongrowing) modes there. Also, non-axisymmetric

modes of finite growth rate and pattern speed in the self-consistent disk become almost treatable.

I can prove that there are no such modes above the velocity dispersions at which the neutral

modes appear, and that any modes that exist below these thresholds cannot possess a bounded

Mellin transform. Unfortunately, I still cannot prove the existence of such modes.

Contents

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2

Closing in on Modes . . . . . . . . . . . . . . . . . . . . . . . . . 2

Modes in the Mestel disks according to Read . . . . . . . . . . . . . . 5

This work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

The Matrix Equation . . . . . . . . . . . . . . . . . . . . . . . . . . 8

The Physics of the Problem . . . . . . . . . . . . . . . . . . . . . . 8

Derivation of the Matrix Equation . . . . . . . . . . . . . . . . . . . 14

Simplifying the Kernel . . . . . . . . . . . . . . . . . . . . . . . . . 23

Integration Over the Azimuthal Angle . . . . . . . . . . . . . . . . . 23

Integration Over the Epicyclic Angle . . . . . . . . . . . . . . . . . . 23

Integration Over the Second Integral . . . . . . . . . . . . . . . . . . 24

Integration Over the Auxiliary Angle . . . . . . . . . . . . . . . . . . 25

Integration Over the Angular Momentum . . . . . . . . . . . . . . . . 27

The Kernel of the Integral Equation . . . . . . . . . . . . . . . . . . 33

The Cut-Out-Disk . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

The Cut-Out Distribution Function . . . . . . . . . . . . . . . . . . . 39

The Boltzmann Equation for the Cut-Out Disk . . . . . . . . . . . . . 40

The Matrix Equation for the Cut-Out Disk . . . . . . . . . . . . . . . 40

Integration Over the Angular Momentum in the Cut-Out Disk . . . . . . . 41

Numerics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45

The Classic Modes . . . . . . . . . . . . . . . . . . . . . . . . . . . 50

Axisymmetric Stability . . . . . . . . . . . . . . . . . . . . . . . . 50

One-Armed Modes . . . . . . . . . . . . . . . . . . . . . . . . . . 53

Two-Armed Modes . . . . . . . . . . . . . . . . . . . . . . . . . . 58

Three-Armed Modes . . . . . . . . . . . . . . . . . . . . . . . . . 63

Four-Armed Modes . . . . . . . . . . . . . . . . . . . . . . . . . . 67

More Modes? . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71

The Singular Disk revisited . . . . . . . . . . . . . . . . . . . . . . 71

Summary and Conclusions . . . . . . . . . . . . . . . . . . . . . . . . 78

Appendix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82

References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86

2Introduction

1. Intro dution

1.1. Closing in on Modes

The riddle why some of the nebular patches seen on the sky appear to have a clearly discern-

able spiral structure puzzles astronomers ever since this spirality was first noticed in by Lord

Rosse (1852).

The first serious attempt at a theory of spiral structure was undertaken by Bertil Lindblad

(1926, 1940). Although his theories were badly flawed, his basic assumption still is the

predominant explanation of spiral structure: A spiral arm is a density wave in which the spiral

is not material as, e.g., the spiral formed by milk being poured into a cup of tea. Instead, the

stars and in particular the interstellar medium (ISM) are slowed down while passing through

the density wave but finally manage to leave it again, so that the constituent stars and clouds

of a spiral arm constantly change. For most density wave theories, a convenient model is a

traffic jam that slowly moves but is made up of different cars at different times since vehicles

slow down at the rear of the congestion and leave it at its head.

It should be noted that there are theories of spiral structure without a density wave, the most

prominent of which probably is stochastic self-propagating star formation (e.g., Jungwiert

and Palouˇs 1994). However, kinematical studies and investigation of the older populations

within spiral galaxies (supposed to be dominating the infrared K-band) show that even in

flocculent galaxies without much global structure some dynamical processes usually are at

work (Thornley and Mundy 1997, Kuno et al 1997), so that while SSPSF might provide a

framework to understand what is going on in irregular galaxies, bona fide spirals are most

likely linked to density waves.

Despite the empirical evidence in favour density waves as the driving agent for spiral structure

and the relatively long time astrophysicists have been developing their theory, dynamical

theories about them are still a muddy ground in which progress is slow and hard. Lindblad’s

original ideas were quickly superseded by the density wave theory connected with the names

of Lin and Shu (Lin and Shu 1964, for a review see Rohlfs 1977 or Bertin and Shu 1996).

Lin, Shu, and their followers simply inserted a density wave into a disk and tried to find

conditions under which this spiral would rigidly rotate, thereby deriving a dispersion relation

linking angular frequency and wave length of a spiral wave.

Although the Lin-Shu theory appeared to explain quite a few observations (Burton 1971,

Rots 1975), its rather phenomenological approach kept the astrophysicists looking for deeper

insight, even more so since gradually more and more shortcomings of the theory surfaced.

Probably the most severe of them was that due to the non-constant dispersion relation a

pure Lin-Shu mode moves inwards or outwards (Toomre 1969). Also, it was found out that

Introduction 3

such modes are absorbed at the Lindblad resonances, i.e., the radii within a disk at which

the angular frequency of the mode is in resonance with the epicyclic frequency (Lynden-Bell

and Kalnajs 1972).

These problems indicated that the global structure of a galaxy, largely disregarded in the

first inceptions of the Lin-Shu theory, plays a significant role in the formation of spiral arms.

During the seventies and the eighties, several theories were worked out that strived to enrich

the Lin-Shu theory with such global features (Bertin et al 1989). A common trait of those

“modal” theories is a feedback loop, in which incoming waves are reflected at the inner border

of the disk or some other feature near the centre. This might be facilitated by the increasing

velocity dispersion near the bulge. By the dispersion relation of Lin-Shu density waves, they

will bend back outwards quite analogous to the bending of short radio waves in the ionosphere

of the earth. This inner reflecting zone is commonly referred to as the Q-barrier, where the

letter Qrefers to Toomre’s stability parameter Q(Toomre 1964), the variation of which plays

the central role in the reflection process.

To close the feedback loop, an outer reflector is required. The most popular candidate

suggested for this outer reflector is the corotation resonance, i.e., the radius at which the

pattern speed of the mode is equal to the angular frequency of a star on a circular orbit. At

corotation, one can reflect (and even amplify) density waves back inwards by a mechanism

dubbed overreflection by the proponents of this model (Mark 1976), or by swing amplification

(Julian and Toomre 1966).

While modal theories do indeed give a framework for the interpretation of observations, they

frequently lack analytic tractability, particularly when the spiral patterns are rather open

and the WKBJ approximation they usually employ breaks down. In the introduction to their

1996 book, Bertin and Shu describe their approach as “semiempirical, that is, our dynamical

studies will always attempt to address issues directly raised by empirical facts [...] We shall

avoid too much speculation on dynamical mechanisms and physical scenarios that do not

have sufficient empirical basis.”

It was felt by many, however, that daring some “speculation” might nevertheless pay off. Pi-

oneering works attempting stability analyses “from first principles” include the local analyses

by Toomre (1964), Goldreich and Lynden-Bell (1965) and Julian and Toomre (1966), and

the global analyses in a series of papers by Kalnajs (1971, 1976a, 1976b, 1977) and in Zang’s

Ph.D. thesis (1976).

The hope of these works was that by simply considering the stability properties of the set

of equations given by the Poisson equation and either the Jeans equations or the Boltzmann

equation one would find that disks develop spiral structure as a part of their nature, as it were,

at least when certain conditions are met, most notable the absence of too much unordered

motion, or, in gas-dynamical terms, pressure support.

The first of the works listed above, Toomre (1964), focused on the axisymmetric stability

of disks in a local analysis, i.e., just looking at a patch of a larger disk, thus developing a

variant of Jeans’ (1919) analysis for a rotating disk, or, more precisely, a shearing sheet. It

was this work that introduced the concept of the stability parameter Qas the ratio between

the velocity dispersion of the disk σand a critical velocity dispersion σcrit below which the

disk becomes susceptible to growing axisymmetric (ringlike) disturbances.

4Introduction

Goldreich and Lynden-Bell (1965) also restricted themselves to a local analysis but examined

non-axisymmetric perturbations. They described their patch by the Jeans equations (“gas-

dynamically”) and found that violent albeit transient instabilities occur in this shearing

medium even when it is stable to axisymmetric perturbations, whereas quasistationary modes

in the sense of Lin and Shu are absent. Julian and Toomre (1966) chose to investigate

the same problem using the Boltzmann equation (“stellar-dynamically”) and found a rather

similar behaviour.

Later, Toomre (1981) coined the term “swing amplification” for this sort of instability and

speculated that at least flocculent galaxies might be shaped by a recurrence of swing am-

plification events. As a matter of fact, Fuchs (1991) showed that by nonlinear coupling

instabilities might trigger new instabilities themselves, thus producing some sort of feedback

loop again. Models along these lines got additional credibility form simulations like those

of Sellwood and Carlberg (1984), who also found the “swirling hotch-potch” of perpetual

growth and decay of spiral fragments.

However, there are galaxies with a beautiful spiral structure that appears to be “global” and

thus cannot be expected to turn up in the local analysis. In particular, disk models based

on plastering the disk with copies of a patch do not possess any borders like the Q-barrier

that appeared so central in the modal picture. Although it was found that some of the most

beautiful spirals can be explained without regarding their internal dynamics as simply being

the result of interactions (Toomre and Toomre, 1972), it still appears that the local analyses

miss important ingredients to a real disk. As Toomre (1977) wrote, it appears likely that

“even in M51 the main spiral structure is in essence just a badly-dented version of what

would have been there anyway.”

At this point the global analyses mentioned above come into play. After the rather abstract

work of Kalnajs mentioned above, Zang (1976) was the first to succeed in a complete stability

analysis of a galactic disk using the Boltzmann equation. Zang chose the Mestel disk (Mestel

1963) as basic state, a disk the rotation curve of which is flat from the centre out to infinity.

While this certainly is a grossly unrealistic model because of its infinite mass and a density

singularity in the centre, it is simple enough to allow detailed calculations that require only

moderate numerical effort. What is more, the part of real galaxies in which spiral structure

is prominent is frequently characterised by a flat rotation curve, so that Zang could indeed

hope that at least some of his results should be readily generalisable to more realistic galaxy

models. One of the remarkable results of Zang was that he could not find “interesting”—in

a sense that will be clarified below—modes in the Mestel disk unless he immobilised a part

of the disk matter in such a way that he essentially introduced a Q-barrier in the disk.

The axisymmetric stability of the self-consistent disk was later reconsidered by Lemos,

Kalnajs and Lynden-Bell (1991) for a gaseous disk. They basically reproduced Zang’s find-

ings, except that the gaseous disk may also admit breathing modes (i.e., modes with no radial

node) if the adiabatic exponent is larger than 3/2, which is impossible for a stellar disk.

One drawback of Zang’s work was that it was never published in generally accessible liter-

ature and that the presentation of the material is not always as clear as one might wish.

Furthermore, it made heavy use of the fact that the Mestel disk is self-similar, i.e., that the

orbits near the centre of the disk are just scaled copies of those further out. Those were

the reasons to embark on this research: To find an easily apprehensible, conceptually and

Introduction 5

computational simple and possibly readily generalisable way to reproduce Zang’s results and

to follow the behaviour of the Mestel disk even further than Zang without having to take

recourse to numerical integration.

While I was working on this project, the accessibility of Zang’s results was greatly enhanced

by Read’s (1997) Ph.D. Thesis and the papers by Evans and Read (1998a, 1998b) based on

it. Read generalised Zang’s method to power law disks of almost arbitrary exponent. The

following chapter is intended to supply the reader with enough background from those works

to see the parallels and differences between their and my approach.

1.2. Modes in the Mestel disks according to Read

In her PhD thesis, Read (1997—henceforth referred to as Read) investigates the stability of

power law disks (her results are also published in Evans and Read (1998a, 1998b)), of which

the Mestel disk studied here and earlier by Zang (1976 – henceforth referred to as Zang) is

a special case. Since Read’s work is published and much more readable than Zang’s, I will

usually refer to her work instead of Zang’s, although most of the results on Mestel disks were

already present in the earlier work. Since these previous results provide the framework in

which to judge mine, it seems appropriate to shortly summarise them here.

Quite as I will do in the present work, Read starts out with computing the orbits of mass

points in her disks. Read extended her investigation to the entire family of physically feasible

power law disks with a surface density Σ(r) = Σ0(r/r0)β−1, of which the Mestel disk is the

β= 0 special case. In contrast to the present work, she chooses to keep exact expressions.

This forces her to introduce an auxiliary integral not reducible to well-known functions.

Using this integral, all the orbit parameters can be expressed in rather short terms. From

this analysis, Read picks the maximum radial velocity of the orbit, U, and the home radius

RHas the integrals of motion in which to perform the calculation.

The distribution function used by Read in the Mestel disk case is the classical one also

employed here, and she investigates “cut-outs” by letting only a part of the matter participate

in the perturbations. Two cases are worked out: (a) the “inner cut out”, where active matter

is removed only in the centre, and (b) the “double cut-out”, where active matter is removed

both in the centre and in the outer parts of the disk. When used in this work, the term

“cut-out disk” is usually supposed to mean doubly cut out disks with inner and outer cut-out

indices N=M= 2 (see chapter 4.1 for the meaning of these parameters and the shape of

the cut-out).

Read goes on to derive the equation that governs the evolution of the power law disks using

the projection method suggested by Kalnajs (1971). Starting out at the Poisson and the

linearised Boltzmann equations, she finds expressions for the response of a disk to an imposed

perturbation. By projection of both an imposed perturbation and the response to it onto the

basis of logarithmic spirals, she arrives at one integral equation for each azimuthal harmonic

m, where different harmonics are decoupled due to the linearity of the problem.

The kernel of this integral equation, called the transfer function Sm(α, α′) by Read (Eq. 3.80),

decomposes into an integral over the angular momentum Lz, into which the dependence on

the size of the orbit is condensed, and a sum over Fourier components of a given orbit with

6Introduction

respect to the logarithmic spirals that control to what extent a given orbit is susceptible to

the perturbation. This decomposition is possible because of the self-similarity of the disk and

is particularly convenient because the integral over Lzcan be carried out analytically for the

distribution functions chosen by Read.

To compute the kernel, one still has to perform an integration over the eccentric velocity U

and determine the Fourier coefficients. Read does this numerically and finds, as the most

striking feature, that the kernel exhibits a strong trailing bias, i.e., Sm(α, α′) is much larger

for α > α′than on the other side of the α=α′diagonal. This means that power flows almost

exclusively from leading to trailing spiral components in these disks.

The solutions of those integral equations can be found by transforming the equation A(α) =

RSm(α, α′)A(α′)dα′to λA(α) = RSm(α, α′)A(α′)dα′, thus writing the problem as an eigen-

value problem (note that αand α′will be used with their roles exchanged in the rest of this

work, and Read’s Smis Kor ¯

Khere; I apologise for the confusion). In that formulation,

one has found solutions to the original equations when λ, the“mathematical eigenvalue” in

Zang’s nomenclature, is unity. By discretizing Sm, one can use standard numerical software

to obtain those eigenvalues.

At this point Read can examine the stability properties of her disks, starting with the axisym-

metric perturbations. She finds that the Mestel disk becomes unstable at σu= 0.378069 vc

with the rotation speed of the disk vc, in close agreement with the local analysis (Toomre

1964). For the doubly cut-out disk with N=M= 2, axisymmetric instabilities set in at

σu= 0.377 vc. As expected, in no case breathing modes are found.

Since the most beautiful galaxies feature prominent two-armed patterns, the first non-

axisymmetric modes one will look for are the bisymmetric ones. For the self-consistent disk

(i.e., the one with no cut-out applied), Read can only find neutral bisymmetric modes with

zero pattern speed. Due to the lack of a time scale in the Mestel disk, the existence of one

growing mode would lead to the expectation that Mestel disks should be unstable to modes

of every pattern speed and growth rate at a time when their temperature was low enough.

However, Read cannot conclusively rule out such a continuum of modes since the integral

equation describing this problem is very hard to treat numerically. In the Mestel disk, neutral

bisymmetric modes set in at σu= 0.17 vc.

By introducing the cut-out, a power-law disk attains a length-, and therfore a time scale.

Thus isolated modes are possible and are indeed found. But like Zang before her, Read finds

that disks in which the density of the active matter declines only slowly are surprisingly stable

against bisymmetric modes. In the doubly cut-out Mestel disk, they set in at σu= 0.205 vc

which is far less than the velocity dispersion needed to stabilise the disk against axisymmetric

perturbations. This changes only when the inner cut-out is made significantly steeper, thus

allowing bisymmetric modes in axisymmetrically stable disks at the expense of probably

unrealistically sharp inner edges.

The dependence of the minimum velocity dispersion on the steepness of the cut-out is consis-

tent with qualitative reasoning based on the Lin-Shu theory. The cut-out acts as a Q-barrier

shielding waves travelling towards the galactic centre from the inner Lindblad resonance (ILR)

that would otherwise absorb them and additionally reflects the waves back, thus opening the

possibility for a feedback loop. Indeed, for the modes found by Read, the ILR always lies

safely within the cut-out (i.e., in zones with little active matter).

Introduction 7

Zang (1976) had already found that Mestel disks strongly favour m= 1-modes over their

m= 2 counterparts. Indeed, m= 1-modes are already growing when axisymmetric modes

are just marginally stable. Unfortunately, however, the marginal stability is quite problematic

in Read’s m= 1 analysis, since the mathematical eigenvalues appear to follow the real axis in

the plane of complex eigenvalues when the growth rate vanishes and the pattern speed tends

to zero. Thus, growing one-armed modes seem to exist at any temperature. This suggestion

is also backed by results for the self-consistent disk.

A concern for m= 1-modes is that they might shift the barycentre of the disk. However, Read

finds that at least for the Mestel disk, this is not the case, since in cut-out disks the density

profile of the mode tends to cancel out the drags to the barycentre. Even the logarithmic

spiral in the self-consistent case has no effect on the barycentre since the perturbation falls

off faster than the unperturbed surface density.

Read also investigates three- and four-armed modes. As a rule, these are even move difficult

to maintain than the two-armed modes. However, the stability behaviour of the m= 4

mode in the cut-out Mestel disk has mathematical eigenvalues at very low pattern speeds

and vanishing growth rate along the positive real axis. For m= 3, this does not occur.

Analogously, the self-consistent Mestel disk has neutral modes for m= 4 even at high velocity

dispersions. Based on analogous behaviour for varying disk index β, Read speculates that

this phenomenon is linked to the emergence of closed orbits with m= 3 symmetry which are

not present in the m= 4 case and which might somehow damp the modes.

1.3. This work

This work set out as an attempt to reproduce Zang’s (1976) results in a more straightforward

way. Even if they have become far more accessible in the meantime due to Evan’s and Read’s

efforts, it is still worthwhile to try and get analytic approximations for the kernel of the

integral equation describing the modes in Mestel disks, not the least in the hope of trying to

unravel sume of the vexating properties of the Mestel disk.

After briefly laying out the basic physics of the problem, I will derive analytic approximations

for the orbits in the Mestel disk and further simplify them in order to find convenient action-

angle coordinates in which I will express the equations. The derivation of the kernel then

requires a few quadratures, the most involved of which is the one over the angular momentum.

Still, all of these quadratures can be carried out analytically.

Having established the kernels for the self-consistent and the cut-out disks, I will show that the

results obtained using this simplified theory are in satisfying accord with the previous results

Zang and Read obtained,where I will go from axisymmetric to four-armed perturbations.

Finally, I tackle the problem Read could not conclusively solve, the issue of growing modes

in the self-consistent disk.

8The Matrix Equation

2. The Matrix Equation

In this section I will state the basic assumptions of the model, derive an approximated Hamil-

tonian and combine it with the Boltzmann equation to arrive at an integro-differential equa-

tion. Instead of writing this equation down, I will immediately convert it into what Kalnajs

(1971) calls a matrix equation by projection onto the space spanned by the logarithmic spirals.

2.1. The Physics of the Problem

2.1.1. The Mestel Disk

The Mestel disk is named after Leon Mestel, who in a work published in 1963 investigated

the gravitational collapse of a spherical cloud of gas, back then being mainly concerned with

star formation. He found that under appropriate assumptions the collapse will lead to a disk

with a surface density profile

Σ(r) = Σ0r0/r, (1)

where r0is a normalising length, at which, obviously, Σ(r0) = Σ0, and ris the radial coordi-

nate with the centre of the disk at r= 0. Using the Poisson equation for a razor-thin disk,

this can immediately be integrated to yield the potential

Φ(r) = 2πGΣ0r0ln r(2)

with the gravitational constant G.

As can be seen from (1), r0is not really a free parameter of the disk, since disks with different

r0can also be viewed as having identical r0but different Σ0. I exploit this liberty to choose

r0to one unit length.

This type of disk is interesting for the investigation of spiral structure not only because of

its analytic simplicity, but also because it has a flat rotation curve, i.e. the circular velocity

vcof a mass point a Mestel disk is independent from its distance to the disk’s centre. At

least nearly flat rotation curves are fairly universal for the off-centre regions of giant spiral

galaxies (e.g., Rubin et al. 1985 or Persic and Salucci 1991 and references therein), with indeed

flatter rotation curves apparently associated with a more dominant global mode (Biviano et

al. 1991).

In this work, I choose the unit velocity such that this circular velocity

vc= 2πGΣ0= 1.(3)

The Matrix Equation 9

In this way, the equations of motion in the Mestel disk in polar coordinates (r,ϑ) are

¨r=r˙

ϑ2−1/r (4a)

Lz=r2˙

ϑ= const,(4b)

where Lzdenotes the angular momentum per unit mass and the dot signifies the temporal

derivative.

2.1.2. Derivation of the Hamiltonian in Ation-Angle Co ordinates

In my programme of solving the Boltzmann and Poisson equations of the perturbed Mestel

disk, I will start with the Boltzmann equation. This equation takes its simplest form when it is

formulated in the action-angle coordinates of the problem. Therefore, it is favourable to start

out with deriving a Hamiltonian of the perturbed Mestel disk in action-angle coordinates.

The most straightforward way to find the transformation equations between the (non-

canonical sets of coordinates (r, ϑ, ˙r, ˙

ϑ) and the yet-to-be-determined action-angle coordinates

(w1, w2, Lz, J) is to solve the equations of motion. Unfortunately, an exact closed solution of

(4) does not exist in elementary functions.

The path most frequently taken to get approximate solutions, the epicyclic approximation,

yields a Hamiltonian that does not contain any terms originating from non-circular motions

and thus is useless for the investigation of perturbations. I therefore employ a different

approximation based on the method of harmonic balance. The method of harmonic balance

belongs to the large family of asymptotic methods, where nonlinear differential equations are

solved by expanding the nonlinearity into a series and dropping higher terms. The choice of

the basis functions for this expansion is crucial for the quality of the solutions.

In the case of primarily oscillatory motion the most straightforward basis functions are sines

and cosines, and thus the coefficients of the series are just the Fourier coefficients of the

nonlinearity. Keeping only the constant and the two lowest order terms (those proportional

to sin(x) and cos(x)), one is left with a linear differential equation that can readily be solved.

More on asymptotic expansions can be found in Bogoliubov and Mitropolski (1961), Kuypers

(1990) gives a very nice and concise treatment of the method of harmonic balance itself.

Inserting (4b) into (4a) yields

¨r=L2

z/r3−1/r. (5)

In order to perform the development of the nonlinearity in (5), I substitute

r(t) = R01 + ζcos(κt)(6)

with an epicyclic frequency κand a parameter ζthat gives the ellipticity of an orbit such

that ζ= 0 indicates a circular orbit. R0is the orbit’s reference radius. The parameter ζis

the small parameter of this asymptotic expansion, so ζshould not exceed 0.1, say, for the

majority of the constituents of the model galaxy. I will return to the issue how much of a

restriction this is in a moment.

10 The Matrix Equation

The first Fourier terms of the right hand side of (5) with respect to tare

a0=L2

z1 + ζ2/2−R2

01−ζ22

01−ζ25/2

a1=−2R2

01−ζ22p1−ζ2−1+ 3L2

zζ2

ζR3

01−ζ25/2

b1= 0,

(7)

where aiand bidenote the cosine and sine coefficients, respectively. Dropping all higher

Fourier terms Eq. (5) now reads

¨r=a0+a1

ζr

R0−1.(8)

I parametrise the orbits by R0and ζ. Inserting (6) and comparing the terms independent of

tone immediately obtains

Lz=√2R0(1 −ζ2)

pζ2+ 2 .(9)

By comparing the terms with cos(κt) the epicyclic frequency κis found to be

κ=√2

ζR0s1−2p1−ζ2

ζ2+ 2 .(10)

To determine the circular frequency of the guiding centre Ω in this approximation one inserts

(6) and (9) into Eq. (4b) and Fourier analyses with respect to tagain. This yields

Ω = Lz

0(1 −ζ2)3/2.(11)

0.4

0.5

0.6

0.7

0.8

0.9

1.1

1.2

0 1 2 3 4 5 6 7

0.95

1.05

1.1

1.15

1.2

orbital phase

Fig. 1: A comparison of nu-

merically integrated and approx-

imated orbits. The upper pair is

for ζ= 0.15 and has to be read

with the left scale, the lower one

for ζ= 0.05 and corresponds to

the right scale. In both cases,

R0= 1.

The Matrix Equation 11

This completely determines the orbits. For what kind of orbit parameters will those approxi-

mated orbits be useful? The temporal derivative of (6) shows that the maximal radial velocity

attained on a given orbit (the eccentric velocity in Zang’s nomenclature) is U=R0ζκ. In

the epicyclic approximation of the Mestel disk, κ=√2Ω, so that U=√2ζvcfor small ζ.

As can be seen in Fig. 1, the method of harmonic balance gives at least useful estimates of

both κand the radial extent of the orbits up to ζ≈0.15. However, the shape of the real orbits

is quite different from a sine already when ζis that high. If our Galaxy were a Mestel disk,

in the solar neighbourhood (vc= 200 ···250 km s−1) that would correspond to an eccentric

velocity of roughly U= 50 km s−1, which again means that its orbit should be confined to

a torus roughly 2 kpc wide. It is probably justified to assume that most of the stars that

significantly participate in a spiral instability will more or less meet this criterion. Of course,

this does not decide the question whether the quite severe distortion of the shape of the orbits

will greatly degrade the results, and it is hard to do this a priori. Still, the fact that the times

spent within and out of the reference radius are quite comparable for the approximated and

the true orbits may restore some confidence shaken by the optical appearance of the ζ= 0.15

orbit.

However, even this approximation is not enough to make the problem of the transformation

to action-angle coordinates tractable, since unfortunately determining the Hamiltonian in

action coordinates from Eqs. (10) and (11) results in rather messy expressions. In particular,

solving (10) for ζgets quite involved. However, for small ζone may linearise (10) and (11)

κ=√2/R0Ω = 1/R0(12)

(this is identical to the result of classic epicyclic theory) and finds on developing the Hamil-

tonian to the second order and dropping higher Fourier coefficients

H0(J, Lz) = 1

2+ ln(Lz) + √2J

,(13)

where J=κζ2R2

0/2. This Jsimply is the energy in the epicyclic motion divided by κ, which,

in an epicyclic approximation, is the second integral and thus the second action besides Lz.

Clearly, this Hamiltonian is the simplest one that could possibly be related to the Mestel

disk. The justification of the rather brutal simplifications I was forced to adopt can (and

will) only be given a posteriori by showing that the disks described by the Hamiltonian (13)

do indeed behave quite like the “true” Mestel disks studied by Zang and Read.

The action variables Jand Lzare conjugate to the angle variables

w1=κt ∂H0

∂J = ˙w1=κ

w2= Ωt∂H0

∂Lz

= ˙w2= Ω

(14)

Although the Hamiltonian (13) is of order ζ2due to the presence of a term with J/Lz, the

canonical equations (14) are not correct to ζ2; in fact, while ∂H0/∂Lz= 1/Lz−√2J/L2

one has Ω = 1/Lz−(3/2)J/L2

zto second order in ζ, a 5% difference in the second-order

12 The Matrix Equation

coefficient (the corresponding relations for w1and Jare correct to second order). However,

the term with Jis important for the derivation of the Boltzmann equation, and terms of

order J/Lz=O(ζ2) will be disregarded after its derivation.

Now, if I had derived the exact action-angle coordinates, the transformation from Carte-

sian coordinates and their canonically conjugate momenta into the space of action-angle

coordinates (w1, w2, J, Lz) would be a canonic transformation. Since (6) and (12) are ap-

proximations, it is almost certainly not. However, if the ”non-canonicity” caused by the

approximations is small enough, i.e., of order O(ζ2), the error made by treating this trans-

formation as canonic is negligible. I will now check that this indeed the case by computing

the Jacobian of the transformation from the set of action-angle coordinates (w1, w2, Lz, J) to

the Cartesian coordinates (x, y) and the momenta per unit mass ( ˙x, ˙y).

The transformation from the action-angle coordinates to the set of polar coordinates and

their temporal derivatives is

r=Lz+q√2JLzcos(w1)

˙r=−√2s√2J

sin(w1)

ϑ=w2−√2s√2J

sin(w1)

ϑ=−1

+ 2p√2JLz

cos(w1).

Inserting this into the transformation form polar to Cartesian coordinates yields

x=Lz(1 + ζcos(w1)) cos w2−√2ζsin(w1)

y= (Lz+ζcos(w1)) sin w2−√2ζsin(w1)

˙x=−√2ζsin(w1) cos w2−√2ζsin(w1)

−(1 + ζcos(w1)) 1−2ζcos(w1)sin w2−√2ζsin(w1)(16)

˙y=−√2ζsin(w1) sin w2−√2ζsin(w1)

−(1 + ζcos(w1)) 1−2ζcos(w1)cos w2−√2ζsin(w1),

where I have kept ζ=q√2J/Lzto save space.

The Jacobian ∂(x, y, ˙x, ˙y)/∂(w1, w2, J, Lz) of this transformation turns out to be

J= 1 + 6s√2J

cos(w1) + 9√2J

cos2(w1) + 4s√2J

cos3(w1).(17)

The Matrix Equation 13

Developing this into a series in ζ, one finds

J= 1 + 6ζcos(w1) + 9ζ2cos2(w1) + O(ζ3).(18)

Thus, only after averaging over one epicyclic period is the Jacobian really of the desired

1+O(ζ2) order. This should not be a major issue in limiting the accuracy of my calculations,

though, in particular because the kernel I will calculate is a quantity averaged over both

periodic motions.

2.1.3. The Distribution Funtion

Binney and Tremaine (1987) give a distribution function of the Mestel disk (without counter-

rotating stars) as

f0(E, Lz) = FBTLq

zeE/σ2, Lz>0

0, Lz<0(19)

with the relative energy E= Ψ −v2/2, the relative potential Ψ, and

FBT =Σ0˜r0

2q/2√πq−1

2!σq+2 , q =σ−2−1.(20)

Noting that Eis just the negative of the Hamiltonian, I can express the distribution function

in Jand Lzas

f0(Lz, J) = F

exp −√2J

σ2Lz!,(21)

where I have collected the terms not dependent on the actions in

F=e−1/2σ2Σ0˜r0

2q/2√πq−1

2!σq+2 .(22)

The equilibrium distribution function fonly depends on actions, as it should by the Jeans

theorem. Again, this is not much of a test for the approximated transformation equations,

but it provides some reassurance.

The parameter qusually is the range 10 ···100 in disks responsive to spiral perturbations,

corresponding to σ≈(0.1···0.3)vcor (20 ···70)km s−1for a disk like our own Galaxy’s in

the regions that show spiral structure. Thus, an application of Stirling’s formula in Fis well

justified and significantly reduces later numerical hassles. The result is

F=√2√1−2σ21−1/σ2

4eπ2σ2G(23)

where (3) has been used to eliminate Σ0and I have set r0= ˜r0.

14 The Matrix Equation

2.1.4. The Potential Basis

Following Kalnajs (1971), I use logarithmic spirals

ϕα,m(r, ϑ) = eiαln(r)+imϑ−iωt

√r(24)

as the basis into which the potential is developed. The factor r−1/2has been added to ensure

definiteness of the Lz-integration in the matrix equation and will silently be replaced by a

L−1/2

zwhen linearising in the actions. This is not exact for a linearisation in ζand has in

that sense to be justified a posteriori. However, the potential basis is just that, a set of basis

functions, and thus slightly changing the shape of the basis functions while retaining their

central features (orthogonality, span) should not be critical anyway and would only influence

the expansion coefficients but not the expansion itself.

The sign of the term with ωhas been chosen so that spirals with ω,m, and αall positive are

trailing. While the main impetus of this work is the study of marginal (i.e., neither growing

nor decaying) modes, it should be noted that growing modes correspond to an ωwith positive

imaginary part with this distribution of signs. This will be of particular importance in the

application of Landau’s rule further below.

The quantity α, appearing here as the tangent of the inclination of a spiral arm, can also be

viewed as something like a radial wave number. Indeed, the larger |α|, the tighter wrapped

the spirals are. However, since αis multiplied with ln(r) instead of rin (24), it is not true

that the spiral pattern is periodic with period 2π/α in r. Instead, 2π/α is the period of ln r,

and thus αshould be called logarithmic wave number. I will nevertheless talk about wave

numbers in this text.

It should be noted that with this definition the pattern speed of the wave is Ωp=ω/m, since

at a given point in the galaxy a potential trough has to pass mtimes for the spiral pattern

to complete an entire rotation around the centre.

The expression (24) has to be reduced to order O(ζ2). Some care is required in doing this since

I want to keep terms like sin(ζsin(w1)), in accord with my programme of keeping oscillatory

terms of order exp(iζ). This is done by linearising just the exponent to order ζ, not the entire

exponential.

After the linearisation the potential basis has the form

ϕα,m(w1, w2, J, Lz) = Liα−1/2

zexp 

imw2+ i s√2J

Lzαcos(w1)−√2msin(w1)−iωt

.

(25)

In this work, I assume all possible potential perturbations have an expansion like

∞

m=0 Z∞

−∞

Φα,mϕα,m(w1, w2, J, Lz)dα. (26)

This is certainly true for all continuous (and therefore bounded for finite r) perturbations

that vanish quickly enough as r→ ∞ to be acceptable in a linear perturbation analysis.

The Matrix Equation 15

2.2. Derivation of the Matrix Equation

2.2.1. The Linearised Boltzmann Equation

The response of the Mestel disk to a perturbation of the type (25) can be computed using

the Boltzmann equation,

dt =∂f

∂t + [f, H] = 0,(27)

where the symbol [ .,.] denotes Poisson brackets, His the total Hamiltonian of the system,

and fits total distribution function. For brevity, I do not indicate the dependence on the

angles and actions in the distribution function and the potential (and their perturbations)

except where it seems necessary for reasons of clarity.

Since I am only interested in small deviations from the equilibrium, I may significantly sim-

plify this equation by linearising it. To this end, I write H=H0+ Φ1and f=f0+f1,

where H0and f0are given by (13) and (21), respectively. I immediately specialise on the

case of a perturbation potential Φ1of the form (25). I furthermore demand that Φ1≪H0

and f1≪f0. Using the linearity of the Poisson brackets, (27) then becomes

∂f1

∂t + [f0, H0] + [f1, H0] + [f0,Φ1] + [f1,Φ1] = 0 (28)

The second summand on the left side of this equation is zero by definition of the equilibrium,

whereas the last summand is neglected since it is quadratic in the small quantities f1and

Φ1. Hence the linearised Boltzmann equation reads

∂f1

∂t + [f1, H] = −[f0,Φ1].(29)

After inserting the canonical equations and (12), the left hand side of (29) becomes

l.h.s.=∂f1

∂t +∂f1

∂w1

˙w1+∂f1

∂w2

˙w2

=∂f1

∂t +√2

∂f1

∂w1

∂f1

∂w2

(30)

The right hand side of (29) is somewhat more involved. Evaluating at first only the derivatives

of f0one gets

r.h.s. = f0

(σLz)2√2Lz

∂Φ1

∂w1

+σ2Lz−√2J∂Φ1

∂w2.(31)

In this expression the term with Jis of order ζ2and is therefore discarded.

Before inserting the potential basis (25), I have to decide on the sign of ϕα,m in the Hamilton

function. Since I want to save me the trouble of maintaining a sign in the Poisson equation

(where it enters with the opposite sign), I use Φ1=−Φα,mϕα,m here, where Φα,m is a complex

16 The Matrix Equation

number and corresponds to the coefficients of the series (26). Carrying out the remaining

derivatives, one arrives at

r.h.s.=−i Φα,m

FLiα−3

σ2q2√2Jαsin(w1) + √2mcos(w1)−mσ2pLz

exp 

imw2+ is√2J

Lzαcos(w1)−√2msin(w1)−√2J

σ2Lz−iωt

.

(32)

2.2.2. Computation of the Resp onse

The solution of the inhomogeneous linear partial differential equation (30)=(32) gives the

response f1(J, Lz, w1, w2, t) of the disk to an elementary forcing (25) in the limit of vanishing

perturbation amplitude. To find a solution, I first convert the Boltzmann equation to an

ordinary differential equation which, due to the linearity of the equation, is easily solved.

The structure of (32) invites a separation of variables by

f1(w1, w2, t) = g1(w1)eimw2−iωt,(33)

which results in the inhomogeneous ordinary linear differential equation

g′

1(w1) + i(ωLz+m)

√2g1(w1) = C αsin(w1) + √2mcos(w1)−mσ2sLz

2√2J!

expiq√2J/Lzαcos(w1)−√2msin(w1)−√2J

σ2Lz.

(34)

Here, the prime denotes the derivative with respect to w1and

C=−i Φα,mq2√2JFLiα−3

σ2.(35)

This inhomogeneous equation is solved in the usual way by finding the solution of the ho-

mogenous equation and then applying the method of variation of constants. The solution of

the homogenous equation is

˜g1(w1) = ei(ωLz−m)w1/√2.(36)

By variation of constants one immediately has

g1(w1) = eiηw1D+Zw1

S(w′

1)dw′

1,(37)

where I have abbreviated η= (ωLz−m)/√2 and Dis an integration constant. The lower

integration bound is arbitrary, but the choice of 0 is quite convenient for the calculations

below. The integrand S(w′

1) is given by

S(w′

1) =C αsin(w1) + √2mcos(w1)−mσ2sLz

2√2J!

exp−iηw′

1+ iq√2J/Lzαcos(w′

1)−√2msin(w′

1)−√2J

σ2Lz

(38)

The Matrix Equation 17

and has the property that S(w′

1+ 2π) = e2πiηS(w′

1). Therefore,

Z2π+w1

S(w′

1)dw′

1=Z2π

S(w′

1)dw′

1+ e2πiηZw1

S(w′

1)dw′

1.(39)

I now have to find the particular solution appropriate to my problem, i.e., fix D. Because (34)

is a first-order differential equation, I only need one condition that in this case is generated

from a periodicity requirement: fhas to be uniquely defined and thus 2π-periodic in w1.

Hence fmust satisfy f(w1) = f(w1+ 2π), which in turn implies g(w1) = g(w1+ 2π).

Inserting (37) and applying (39) one arrives at

D=1

e−2πiη−1Z2π

S(w′

1)dw′

1.(40)

Note that, in contrast to the usual theory of linear differential equations with initial or

boundary conditions, this requirement does not necessarily cut down the solution space to

a single point. For example, requiring f(w1) = f(w1+π) will lead to different functions

that also satisfy f(w1) = f(w1+ 2π). It might be a worthwhile exercise to investigate such

solutions.

The integral in (37) can also be somewhat simplified—at least as far as computability is

concerned—by using periodicity arguments. Here one uses

Zw1

S(w′

1)dw′

1=Z2π

S(w′

1)dw′

1+ e2πiηZw1

S(w′

1)dw′

1−Z2π

S(w′

1+w1)dw′

implying Rw1

0S(w′

1)dw′

1=R2π

0S(w′

1)dw′

1−R2π

0S(w′

1+w1)dw′

1/1−e2πiηand finds on

inserting (40)

g1(w1) = eiηw1

e−2πiη−1Z2π

S(w′

1+w1)dw′

1(42)

This integral can be further simplified by noting that the trigonometric functions of w1

in the coefficient appear in the same linear combination as those in the exponent derived

w.r.t. w1, which is hardly surprising, since they came into the coefficient by a derivation of

the exponent. This invites partial integration in order to get rid of the trigonometric terms

outside the exponential. Setting

u=−i

Cs√2J

exp i (w1+w′

1)mσ2

√2+η

v= exp 

imσ2

√2(w1+w′

1) + is√2J

Lzαcos(w1+w′

1)−√2msin(w1+w′

1)−√2J

σ2Lz



(43)

one may write S=v′/u (the prime now denotes derivation w.r.t. w′

1).

18 The Matrix Equation

Assembling the result one is left with

g1(w1) = CsLz

√2Jexp −√2J

σ2Lz! i expis√2J

Lzαcos(w1)−√2msin(w1)

−1

e−2πiη−1√2η+mσ2

Z2π

expis√2J

Lzαcos(w1+w′

1)−√2msin(w1+w′

1)−iηw′

1dw′

1!.

(44)

Substituting this equation into (33) gives the desired response of the disk.

2.2.3. The Poisson Equation

After I have computed the response to an elementary forcing, I proceed to ask when this

response is self-consistent, i.e., for which parameters the density perturbation resulting from

the change in the distribution function causes just the potential perturbation that I started

with. To answer this question, the Poisson equation linking the potential Φ and the density

△Φ = −4πGρ

with Newton’s gravitational constant G must be solved simultaneously to the Boltzmann

equation. Since the equation is linear, it is immediately valid for the perturbation quantities

just like for the total potential and density.

Intuitively, one expects that the elementary forcing (25) will as a rule have a rather different

shape after having been mangled through the Boltzmann and Poisson equations, i.e., the

response to (25) in potential space contains many contributions of the form (25) and has the

form (26).

Indeed, the Laplace operator in the Poisson equation does not commute with the operator

in the Boltzmann equation, and thus they will not possess common eigenfunctions. What

is more, the Laplacian expressed in action-angle coordinates is quite cumbersome. Thus,

the integro-differential equation resulting from a simple insertion of the Boltzmann into the

Poisson equation would be very hard to handle. To overcome this difficulty, I follow Kalnajs

(1971) in avoiding this integro-differential equation by converting the equation system into a

matrix equation. This can be done by projecting the potential perturbation and the density

response onto a set of suitable basis functions.

The basis functions in potential space have been chosen above to the be the logarithmic

spirals (24). To find the corresponding basis in density space, I employ Toomre’s device

(Clutton-Brock 1972). One has

△e−k|z|=k2e−k|z|−2kδ(z),(46)

and, by definition of the Bessel function,

△Jm(kr)eimϑ=−k2Jm(kr)eimϑ.(47)

The Matrix Equation 19

Since the Laplace operator in cylinder coordinates,

△=1

∂

∂r r∂

∂r+1

∂2

∂ϑ2+∂2

∂z2,(48)

decomposes into independent operators in (r, ϑ) and z, using (46) and (47) one can easily

verify that the functions Φkand ρkdefined by

Φk(r, ϑ, z) = Jm(kr)eimϑe−k|z|ρk(r, ϑ, z) = k

2πGJm(kr)eimϑδ(z) (49)

form a potential-density pair. Multiplying by a suitable function g(k)—it is completely

unrelated to the function g1(w1) introduced in the last chapter—and integrating over kone

sees that

ρ(r, ϑ, z) = eimϑ

2πG Z∞

Jm(kr)g(k)kδ(z)dk

Φ(r, ϑ, z) = eimϑ Z∞

Jm(kr)g(k)e−k|z|dk

(50)

form a potential-density pair as well.

Multiplying the second line of (50) with rJm(k′r), integrating over rand applying the Fourier-

Bessel-Theorem Z∞

Jm(k′r)Jm(kr)r dr =δ(k′−k)

k(51)

yields

g(k) = kek|z|−imϑ Z∞

Φ(r, ϑ, z)Jm(kr)r dr (52)

(k′has been renamed khere).

The resulting g(k) can now be inserted into the first line of (50). After an integration over z

one has

µ(r, ϑ) = 1

2πG Zk2Jm(kr)ZΦ(r′, ϑ, 0)Jm(kr′)r′dr′dk, (53)

where µdenotes a surface density. Inserting (24), for the moment omitting the dependence

on t, one has integrals of the type R∞

0xµJν(ax). The integrals can be carried out using

Gradshteyn and Ryzhik (1979), 6.561.14,

Z∞

xµJν(ax)dx = 2µa−µ−1Γ(1 + ν+µ)/2

Γ(1 + ν−µ)/2,(54)

where Gradshteyn and Ryzhik require that −(Re ν+ 1) <Re µ < 1

2and a > 0, which

given the asymptotics of the Bessel function is of course necessary for the definiteness of the

integral. In the present context, however, it is permissable to carry this formula to the limits

Re µ→1/2 and a→0.

20 The Matrix Equation

For the r′-integration, one sets

x=r′ν=m µ =1

2+ i α a =k,

and for the integration over k

x=k ν =m µ =1

2−iα a =r.

This gives the basis functions in surface density,

µα,m(r, ϑ) = Km(α)eimϑr−3/2+iα,(55)

where

Km(α) = Γ3

4+1

2m+1

2iα

πGΓ1

4+1

2m+1

2iα

2(56)

is related to the Kalnajs gravity factor as used by Read, K(α, m), through Km(α) =

2πGK(α, m)−1. In writing down Km(α) in terms of moduli, I have already exploited

the fact that Γ(¯z) = Γ(z).

The two sets of functions µα,m and ϕα,m both form a basis for two function spaces. If

both a potential perturbation Φ(r, ϑ) and the density response computed by projection of the

distribution function resulting from the Boltzmann equation into the (r, ϑ)-space are elements

of those spaces, they can be expanded into two series,

Φ(r, ϑ) = X

mZΦα,mϕα,m(r, ϑ)dα

µ(r, ϑ) = X

mZMα,mµα,m(r, ϑ)dα.

(57)

By definition of ϕα,m and µα,m as an orthonormal potential density pair, the Poisson equation

for the coefficients of the series (57) takes the simple form

Φα,m =Mα,m.(58)

For later reference, it should be noted that the scalar product of two corresponding basis

functions is

Z2π

0Z∞

ϕ∗

α,mµα,mr dr dϑ = 4π2Km(α)δm,m′δ(α−α′).(59)

The Matrix Equation 21

2.2.4. The Matrix Equation

I am now ready to write down the matrix equation. To this end (42) can be written as

fresponse

1=LΦimposed (60)

with a linear operator L. Assuming for the moment that this equation is formulated in

Cartesian coordinates, the perturbation of the distribution function f1can be transformed

into a density perturbation µ1by integration over the momenta per unit mass ˙xand ˙y.

Inserting the expansions for µresponse and Φimposed, for the moment assumed to be expressed

in Cartesian coordinates as well, one finds

mZdα µα,m(x, y)Mresponse

α,m =ZdpxdpyX

mZdα Φimposed

α,m Lϕα,m(x, y).(61)

In this chapter only, I write the differentials directly behind their integral signs to improve

the readability of the expressions.

A projection of this equation onto ϕα′,m′yields

4π2Km′(α′)Mresponse

α′,m′=Zdx dy ZdpxdpyX

mZdα Φimposed

α,m (Lϕα,m)ϕ∗

α′,m′.(62)

One advantage of Kalnajs’ method is that the projection can be carried out in any set of

coordinates. Since canonical transformations have a unity Jacobian, I may directly replace

the impulse-space coordinates with the action-angle coordinates derived above.

Inserting the Poisson equation (58) yields

4π2Km′(α′)Φresponse

α′,m′=Z∞

−∞

dα Z∞

dLzZ∞

dJ Z2π

dw1

∞

m=0 Z2π

dw2Φimposed

α,m (Lϕα,m)ϕ∗

α′,m′

.(63)

I have silently swapped the order of integrations into something convenient for me. The

integration limits on Lzwould be −∞ and ∞in the general case. Here, however, I demanded

that the distribution function was identically zero for Lz<0, hence the smaller integration

interval.

I am looking for self-consistent solutions, and thus the stability problem is solved when

Φimposed = Φresponse. Before looking for such solutions, I first carry out the integrals over

the action and angle coordinates. The structure of the integrand in (63) invites a separation

of this task into the integration of a part independent of w′

1and one containing the integral

over the auxiliary angle w′

1. After solving (63) for Φα′,m′, these integrands take the forms

I1=Φα,m

FLi(α−α′)−2

4π2σ2Km′(α′)

exp 

is√2J

Lz(α−α′) cos(w1)−√2(m−m′) sin(w1)+ i(m−m′)w2−√2J

σ2Lz

,

(64)

22 The Matrix Equation

and

I2=iΦα,m

√2FLi(α−α′)−2

zmσ2+√2η

8π2σ2Km′(α′)(e−2πiη−1) exp i(m−m′)w2−iηw′

1−√2J

σ2Lz

+ is√2J

Lzαcos(w′

1+w1)−α′cos(w1)−√2msin(w′

1+w1) + √2m′sin(w1)!,

(65)

respectively. Here, I have written the complete right-hand sides, not indicating for the four

(or, in the case of I2, five) quadratures and the one summation required to compute the

matrix elements of the kernel. I will use this convention throughout the following and call

the right-hand sides emerging after the integration over the variable xas I(x)

1,2, not writing

down the remaining integral and summation signs.

Again adopting a sloppy stance towards exchanging integration orders, the complete integral

equation now is

Φα′,m′=Z∞

−∞

dα ∞

m=0 Z∞

dLzZ2π

dw′

1Z∞

dJ Z2π

dw1Z2π

dw2(I1+I2).(66)

Simplifying the Kernel 23

3. Simplifying the Kernel

In this section I will carry out the quadratures in (66). This will result in an integral-free

expression for the kernel in elementary and special functions. One series introduced in this

process will not be representable in well-known functions. Concluding this section, I will

discuss some general features of the kernel.

3.1. Integration Over the Azimuthal Angle

The first integration, R2π

0dw2, is straightforward. In both I1and I2, the only term dependent

on w2is

Z2π

ei(m−m′)w2dw2= 2πδm,m′.(67)

The Kronecker delta decouples modes of different m, which was to be expected in a linear

analysis. I will, however, postpone the summation over mand thus its elimination until all

integrations have been carried out, since keeping mdoes not overly complicate the compu-

tations and might give some hints what should happen in the mildly nonlinear regime when

modes of different mstart to couple.

3.2. Integration Over the Epicyclic Angle

The next step is the integration over w1. The integrals to be computed are

Z2π

expiq√2J/Lz(α−α′) cos(w1)−√2 (m−m′) sin(w1)dw1(68)

for I1and

Z2π

exp iq√2J/Lzαcos(w′

1)−√2msin(w′

1)−α′cos(w1)

−αsin(w′

1) + √2mcos(w′

1)−√2m′sin(w1)!dw1

(68)

for I2.

24 Simplifying the Kernel

This type of integral can be carried out by noting that the integrands can be rewritten

eAsin(w1)+Bcos(w1)= exp pA2+B2A

√A2+B2sin(w1) + B

√A2+B2cos(w1)!

= exppA2+B2cos(v) sin(w1) + sin(v) cos(w1)

= exppA2+B2sin(w1+v),

(70)

where Aand Bare arbitrary coefficients and vhas to be chosen such that sin(v) =

B/√A2+B2and cos(v) = A/√A2+B2. Since

Z2π

esin(w1+v)dw1=Z2π

esin(w1)dw1+Z2π+v

2π

esin(w1)dw1

=Z2π

esin(w1)dw1+Zv

esin(w1)dw1

=Z2π

esin(w1)dw1,

(71)

one may drop the vfrom the argument of the sine in (70). By expansion into Bessel functions

and using the fact that R2π

0sin(sin(w1)) dw1vanishes for symmetry reasons one finds

Z2π

eiAsin(w1)dw1= 2πJ0(A).(72)

Actually, this result is readily found in any integral table; however, since I am going to apply

the very same technique for more complex terms later, a broader presentation in this simple

case seemed justified.

An application of (70), (71), and (72) to (68) and (69) yields

I(w1)

1= Φα,m

Fδm,m′Li(α−α′)−2

σ2Km′(α′)J0

s√2J

Lzp(α−α′)2+ 2(m−m′)2

exp −√2J

σ2Lz!

(73)

and

I(w1)

2=iΦα,m

Fδm,m′Li(α−α′)−2

zωLz−m(1 −σ2)

√2σ2Km′(α′)

eiw′

1(ωLz−m)/√2

(e−√2πi(ωLz−m)−1)

exp −√2J

σ2Lz!J0s√2J

Lzα2+α′2+ 2(m2+m′2) + 2√2 sin(w′

1)(mα′−m′α)−

2 cos(w′

1)(αα′+ 2mm′)1/2.

(74)

Simplifying the Kernel 25

3.3. Integration Over the Second Integral

I next perform the integration over J. The terms in question are

Z∞

J0q√2J/Lzp(α−α′)2+ 2(m−m′)2exp −√2J

σ2Lz!dJ (75)

for I1and

Z∞

J0q√2J/Lzα2+α′2+ 2(m2+m′2) + 2√2 sin(w′

1)(mα′−m′α)

−2 cos(w1)(αα′+ 2mm′)1/2exp −√2J

σ2Lz!dJ.

(76)

for I2.

Again, the two integrals are structurally equivalent and can be computed using

Z∞

xν+1e−αx2Jν(βx)dx =βνe−β2/(4α)

(2α)ν+1 (77)

(Gradshteyn and Ryzhik 1979, 6.631.4) by letting x=√Jand ν= 0. The results are

I(J)

1= Φα,m

Fδm,m′Li(α−α′)−1

√2Km′(α′)exp −σ2

4(α−α′)2+ 2(m−m′)2(78)

and

I(J)

2=iΦα,m

Fδm,m′Li(α−α′)−1

zωLz−m(1 −σ2)

2Km(α) e(σ2/4(α2+α′2)+2(m2+m′2)

e−iw′

1(ωLz−m)/√2

e−√2πi(ωLz−m)−1

exp−σ2

2√2 sin(w′

1)(mα′−m′α)−cos(w′

1)(αα′+ 2mm′).

(79)

26 Simplifying the Kernel

3.4. Integration Over the Auxiliary Angle

The integration over w′

1is only necessary in I2. The term in question is

Z2π

expσ2

2αα′+2mm′cos(w′

1)+ σ2

√2mα′−m′αsin(w′

1)+i w′

1(ωLz−m)/√2dw′

1.(80)

This integrand resembles the one treated in the integration over w1, where the case is more

complicated here because of the presence of a linear term in the exponent that unfortunately

breaks the periodicity of the integrand that could be exploited in (71). Applying the method

outlined above, the integral attains the structure

Z2π

cos(τw′

1)eνcos(v−w′

1)+ i sin(τw′

1)eνcos(v−w′

1)dw′

1.(81)

Here,

τ=−(ωLz−m)/√2

v= arctan√2(αm′−α′m), αα′+ 2 mm′

ν=−σ2

2q(α2+ 2m2)(α′2+ 2m′2).

(82)

One has to use the two-argument arctan in the definition of vsince the parameters can put

vin any quadrant. To get rid of the arctan, I set

Λ = αα′+ 2mm′+√2 i (αm′−α′m)

q(α2+ 2m2)(α′2+ 2m′2),

(83)

after which I can write v=−ln(Λ). The expansion of expνcos(v−w′

1)into Bessel functions

is (Abramowitz and Stegun 1968, 9.6.34)

eνcos(v−w′

1)= I0(ν) + 2 ∞

k=1

Ik(ν) cosk(v−w′

1).(84)

After expanding the cosine, one can again exchange the order of integration and summation,

and using the usual symmetry arguments one finds

Z2π

cos(τw′

1)eνcos(v−w′

1)+ i sin(τw′

1)eνcos(v−w′

1)dw′

=−i I0(ν)e2πiτ−1

τ+ 2 ∞

k=1

Ik(ν)e2πiτ−1

τ2−k2ksin(v) + i τcos(v).

(85)

Simplifying the Kernel 27

Using this, the integral over w′

1computes to

I(w′

2=−Φα,m

Fδm,m′Li(α−α′)−1

zωLz−m(1 −σ2)

√2Km′(α′) e(σ2/4)(α2+α′2)+2(m2+m′2)

I0(σ2

2√(α2+2m2)(α′2+2m′2))1

ωLz−m+

∞

k=1

Ik(σ2

2√(α2+2m2)(α′2+2m′2))Λk

ωLz−m+√2k+Λ−k

ωLz−m−√k!.

(86)

3.5. Integration Over the Angular Momentum

In addition to splitting the integration into I(J)

1and I(w′

2, it is advantageous to carry out

the Lz-quadrature of I(w′

2separate for the part without and the part with the sum.

3.5.1. Integration of

I(J)

For the integration over Lzin I(J)

1one only has to compute R∞

0Li(α−α′)−1

zdLz, which after

the substitution u= ln(Lz) is

Z∞

−∞

eI(α−α′)udu = 2πδ(α−α′).(87)

Thus,

I(Lz)

1=√2πΦ(α, m)Fδm,m′δ(α−α′)

Km′(α′)exp −σ2

4(α−α′)2+ 2(m−m′)2(88)

3.5.2. Integration of

, Part Without Sum

For the term not containing the sum in I(w′

2, one has to compute

Z∞

Liγ−1

ωLz−m(1 −σ2)

ωLz−mdLz,(89)

where for brevity γ=α−α′. After a substitution u= ln(ωLz/m) this becomes

m

ωiγZ∞

−∞

eiγu

eu−1eu−(1 −σ2)du (90)

28 Simplifying the Kernel

2π

3π

−RR

Fig. 2: Integration path for the

Lz-Integration over I(J)

1. The

crosses mark the poles relevant

for the integration.

when ω > 0. Although this expression is not overly complex, it pays off to first analyse the

simpler problem of computing

Z∞

−∞

e(iγ+λ)u

eu−1du (91)

for λ∈ {0,1}, since (90) nicely decomposes into integrals of this type and later integrals,

though more complex, can be evaluated along the same lines. I will perform this integration

using the residue theorem.

The first thing to note about (91) is that the integral does not exist in a conventional sense,

since the integrand does not vanish as u→ ±∞. However, (91) does exist as a distribution

limit, i.e., it is meaningful to compute an integral over, say, γof a suitable function multiplied

with (91). The distribution limit is well known from integrals like R∞

−∞ exp(iωt)dt frequently

occurring in Fourier theory. For the moment, I will pretend that the integrand vanishes faster

than 1/u for u→ ∞ and return to the question how to handle the tails later.

The second complication in the computation of (91) is that the integrand diverges at u= 0.

The most straightforward way to handle this would be to take a Cauchy principal value,

which, however, leads to nonphysical results. This was noted in a plasma-physical context by

Landau (1946), who, in the same publication, also gave the correct prescription for handling

this sort of situation. Landau’s idea was that by requiring that in the very distant past the

system was in equilibrium, one may conclude that any mode must have started out as a

growing mode. If it has evolved continuously, any integration contours must also have been

deformed continuously during the evolution of the system that has lead to the quasi-steady

state characterised by a purely real ω. Consequently, a real ωhas to be seen as the limit of

ω+ iǫ(for my sign convention in (24)) as ǫ→0 and ǫ > 0. This prescription is generally

known as Landau’s rule. A welcome side-effect of applying it is that the expressions to be

derived shortly will also be valid for growing modes (complex ωwith positive imaginary part).

Giving ωa minute positive imaginary part results in a shift of the poles of the integrand in

(101) into the lower half plane, Lz=m(ω−iǫ)/(ω2+ǫ2). Equivalently, the substitution

leading to (90) would put the integration path into the upper half plane if ωhas a positive

imaginary part. Thus, the correct application of Landau’s rule here is to integrate around

Simplifying the Kernel 29

the pole on the real axis into the upper half-plane. The resulting integration path for γ > 0

is shown in Fig. 2.

A third problem complicates the use of the residue theorem in the present case. On letting R

in Fig. 2 to infinity, the path should be continuously deformable. This is not the case here,

where a pole is in the way every 2πi. However, the poles of the integrand get very sharp

indeed, so that ignoring just “a few” points the contour sweeps over during the deformation,

one can continuously deform the contour and still gets a rapid decay of the modulus of the

integrand on the upper line of the path in Fig. 2. A bit more precisely put, one can show that

for any ǫ > 0 there is an R0such that the integral R∞+iR

−∞+iRis smaller than ǫfor all R > R0

and |R−2πij|> ǫ.

Another way to view this is to notice that the residues (93) rapidly decrease with growing

j. Thus, even when one fixes the upper horizontal line in the integration path Fig. 2 and

only lets the vertical parts move towards ±∞, the error one makes is quite small when the

upper horizontal line is put far enough into the complex plane avoiding the poles, since both

(1) the contribution from the contour integral along the upper horizontal line, and (2) the

contribution from the residues above the line will be very small. Numerical checks show

that this somewhat sloppy treatment delivers correct results. In these numerical checks, the

periodic tails have to be taken care of. This can be done by damping the integrand with a

Gaussian of width aand then letting agrow until the results converge. Usually, however,

it works better to numerically compute a few (30, say) definite integrals with boundaries

equally distributed over a period of the tail, where the boundaries should be far enough from

zero to ensure they are in the region where the integrand is purely periodic for all practical

purposes. The mean value of those definite integrals then is an approximation for the integral

without the contribution form the periodic tail.

The denominator in (91) has poles of order 1 at

u(1)

j= 2πij. (92)

The residues at these points are easily computed using res(f, z) = limu→zf(u)(u−z). An

application of de l’Hospial’s rule yields

res e(iγ+λ)u

eu−1, u(1)

j= e−2πij(iγ+λ−1).(93)

When γ > 0, the residues from (93) have to be summed up for j= 1 ···∞ (starting at 1

ensures that Landau’s rule is followed). Because the terms with λin the series evaluate to 1,

this leads to a single geometric series,

∞

j=1

e−2πγj =−1 + ∞

j=0

e−2πγj =1

1−e−2πγ −1 = 1

e2πγ −1(94)

and one finds

Z∞

−∞

e(iγ+λ)u

eu−1du =2πi

e2πγ −1−border terms.(95)

30 Simplifying the Kernel

The same result is found for γ < 0, where this time the summation limits are j= 0 ···−∞:

−∞

j=0

e−2πγj =∞

j=0

e2πγj =1

1−e2πγ .(96)

The “border terms” in (95) arise from the harmonic tails of the integrand in (91), and they

enter negatively since they contribute to the contour integral off the real axis. The harmonic

tails of the integrand are

e(iγ+λ)u

eu−1∼−eiγu;λ= 0 and u→ −∞

eiγu;λ= 1 and u→ ∞.(97)

In the other two cases the integrand falls off exponentially.

Thus, the vertical parts of the integration contour in Fig. 2 contribute

lim

x→−∞Zx

x+i∞−eiγu du =−lim

x→∞

ieiγx

γat −∞, and

lim

x→∞Zx+i∞

ieiγu du = lim

x→∞

ieiγx

γat +∞,

(98)

where I have substituted u=x+ iy.

As briefly mentioned above, the limits in (98) have to be taken in the sense of a distribution

with the following integration over αin mind. In the distribution limit,

lim

x→−∞

eiγx

γ=−iπδ(γ)

lim

x→∞

eiγx

γ= iπδ(γ).

(99)

Expressions of this type have to be read as a shorthand for

lim

x→∞Z∞

−∞

f(γ)eiγx

γdγ =−iπf(0),(100)

say, as is customary in distribution theory. For a concise and not too technical treatment of

distribution theory, see Wladimirow (1972).

Inserting this result and (95) into (90) finally yields

Z∞

Liγ−1

ωLz−m(1 −σ2)

ωLz−mdLz=πm

ωiγ2iσ2

e2πγ −1+δ(γ)(2 −σ2).(101)

Simplifying the Kernel 31

3.5.3. Integration of

, Part With Sum

In the part of I(w′

2with the sum, one has to compute

Z∞

Liγ−1

ωLz−m(1 −σ2)

√2Λ2k

ωLz−m+√2k+1

ωLz−m−√2kdLz(102)

Again substituting u= ln(ωLz/m), this takes the form

m

ωiγZ∞

−∞

eiγu eu+σ2−1

√22ξ2−(1 −eu)2Λ2k(1 + √2ξ−eu) + 1 −√2ξ−eudu, (103)

where ξ=k/m. As above, I will first tackle the technically easier but conceptually equivalent

problem of evaluating

Z∞

−∞

e(iγ+λ+1)u

22ξ2−(1 −eu)2du (104)

for λ∈ {−1,0,1}.

The integral (104) shares all the problems that were to be solved in the previous paragraph:

it converges only in the sense of a distribution limit, it has a pole on the integration interval,

and the integration contour is not continuously deformable to infinity. Luckily, the case is

rather analogous, so the basic arguments and techniques remain the same.

I start out with an analysis of the poles, i.e., the zeroes of the denominator of the integrand

of (104). These lie at

u(2)

j= ln(1 + √2ξ) + 2πijand u(3)

j= ln(1 −√2ξ) + 2πij. (105)

Of course, when ξ > 1/√2, the poles at u(3)

jare no longer situated on the real axis but

are at ˜u(3)

j= ln(1 −√2ξ) + (2j+ 1)πi. I will handle this in a moment. The case ξ= 1/√2

cannot occur since ξby definition is a rational number. Note, however, that the √2 originates

from the ratio of epicyclic and azimuthal frequency of the orbit. Had I kept the more exact

expressions (9), (10), and (11), there could be values of ζfor which that ratio is rational

which in turn could displace the poles at u(3)

jto −∞ for some kand m. In Read’s formalism,

this situation does occur; at least in her case, it did not require significant changes to the

procedures to solve the integral equation. Quite possibly simply ignoring the problem of poles

on the integration contour would be sufficient in a higher-order analysis with my formalism

as well.

The residues at the poles u(2)

jand u(3)

jare again of order 1 and can be computed following

the prescription given for u(1)

j. The results are

res e(iγ+λ+1)u

22ξ2−(eu−1)2;u(2)

j!=−eln(1+√2ξ)(iγ+λ)e2πij(iγ+λ)

4√2ξ

res e(iγ+λ+1)u

22ξ2−(eu−1)2;u(3)

j!= eln(1−√2ξ)(iγ+λ)e2πij(iγ+λ)

4√2ξ

(106)

32 Simplifying the Kernel

2π

3π

−RR

ln 1 + √2ζ

ln 1−√2ζ

Fig. 3: Integration path for the

Lz-Integration over I(w′

2. The

crosses mark the poles relevant

for the integration. The dia-

monds indicate where the poles

lie in the case ξ > 1/√2. While

it does not really matter, one

could let the integration contour

follow the real axis at ln(1−√2ξ)

in this case.

Both of those residues have to be summed up for j= 1 ···∞ when γ > 0. This is straight-

forward and yields

Z∞

−∞

e(iγ+λ+1)u

22ξ2−(eu−1)2du = iπeln(1−√2ξ)(iγ+λ)−eln(1+√2ξ)(iγ+λ)

2√2ξ(e2πγ −1) + border terms.(107)

Now, in the case ξ > 1/√2 the poles u(3)

jare replaced by poles at ˜u(3)

j. Here, for γ > 0

all residues at poles in the (open) upper half plane have to be added up. To carry over the

results from above, one can set the branch cut of the logarithm in the definition of u(3)

jin

(105) such that they ”shift down” by −iπwhen going from ξ < 1/√2 to ξ > 1/√2. To do

this, one sets ln(−1) = −iπ, which, when ln(1) = 0 is to be preserved, suggests a branch cut

at the positive imaginary axis (π

2≥ln(eiϕ)>−3πi

2). With this choice, (107) remains valid

for all rational and positive ξ.

The border terms this time arise from the λ=−1 contribution at u=−∞ and the λ= 1

contribution at u=∞. The λ= 0 contribution is well behaved at both u=±∞. An

evaluation of the contributions following the path lined out in the treatment of I1

border terms = 









πδ(γ)

4ξ2−2λ=−1

πδ(γ)λ= +1

(108)

Thus the contribution from the part with the series of I2is

Z∞

Liγ−1

z2m+√2σ2(ωLz−√2m)Λ2k

ωLz−√2m+k+1

ωLz−√2m−k=

−i√2πω−iγ

e2πγ −1(√2k−mσ2)Λ2ke(iγ−1) ln(m−√2k)

−(√2k+mσ2)e(iγ−1) ln(m+√2k)

+πδ(γ)

√2m

ωiγ 2√2k−m(1 + σ2)

√2k−mΛ2k+2√2k+m(1 + σ2)

√2k+m!.

(109)

Simplifying the Kernel 33

Since the I(Lz)

2is a rather lengthy expression and furthermore naturally decomposes into a

part with and without the Dirac delta, I postpone assembling it until I write down the kernel

in chapter 3.6.

3.5.4. The ase

ω= 0

The reasoning in the previous three paragraphs breaks down when ω= 0. Indeed, as Read

(1997) has already noted, the ω= 0 case cannot be obtained by taking the limit of the

ω6= 0 expressions, mainly, of course, because this limit does not exist—expiγln(ω)has an

essential singularity at ω= 0 when γ6∈ Q, so according to Picard’s theorem the image of an

arbitrarily small neighbourhood of zero under the kernel as an analytic function of ωis at

least dense in the complex plane.

To correctly evaluate the kernel when ω= 0, one therefore has to go gack to an expression

before the ω−iγemerges and take the limit there; in the present problem the most convient

choice for these expressions are I(J)

1and I(w′

2. While I(J)

1does not depend on ω,I(w′

2does.

Setting ω= 0 there leaves Li(α−β)−1

zas the integrand, and thus one has

I2(Lz,ω=0) =−√2πΦα,mδ(α−α′)Fδm,m′m(1 −σ2)

Km′(α′) e(σ2/4)(α2+α′2)+2(m2+m′2)

I0(σ2

2√(α2+2m2)(α′2+2m′2))

m+∞

k=1

Ik(σ2

2√(α2+2m2)(α′2+2m′2))Λk

m−√2k+Λ−k

m+√k!.

(110)

3.6. The Kernel of the Integral Equation

After having carried out the quadratures over the actions and angles in (64) and (65), I can

now assemble the integral equation in the form I will use to examine it.

The first observation to make is that there are contributions proportional to δ(α−α′) in I(Lz)

as well as in I(Lz)

2coming from both the part without the series (101) and the part with the

series (109). One can immediately integrate these contributions over α, and thus the integral

equation resulting will have the form

Φα′,m′=F(α′)Φα′,m′+Z∞

−∞ K(α, α′)Φα,m′dα, (111)

where the summand with Foriginates from the terms proportional to δ(α−α′).

If Kwas a compact operator, this would be a homogenous Fredholm equation of the second

kind. According to the Fredholm Alternative, this would imply that solutions exist uncondi-

tionally, unless one happens to hit an eigenfunction of K. In that case, additional solvability

conditions would have to be satisfied. However, the Fredholm Alternative is only valid for

integral equations with compact integral operators. Now, because of the infinite integration

region the operator in (111) certainly is not compact. What is worse, it will shortly be seen

34 Simplifying the Kernel

that Khas a Cauchy-type singularity on the α=α′diagonal, i.e. K(α, α′)∼1

α−α′as α→α′

Under those circumstances, the solution theory of integral equations is significantly different.

I will return to this issue in chapter 6.

Carrying out the sum over mand the integral over αin I(J)

1, one arrives at

F1=√2πF

Km′(α′).(112)

A similar procedure for the parts of (101) and (109) after a substitution into I(w′

2yields

F2=πF(2 −σ2)

√2Km′(α′)I0σ2

2(α′2+ 2m′2)e−σ2(α′2+2m′2)/2

F3=√2πF

Km′(α′)e−σ2(α′2+2m′2)/2∞

k=1 m′2(σ2−1)

2k2−m′2+ 1Ikσ2

2(α′2+ 2m′2).

(113)

In F3, one can evaluate a part of the series (that converges absolutely because of the Ik) using

(84) written as 2 P∞

k=1 Ik(z) = ez−I0(z). This leads to

F3=√2πF

Km′(α′) −1 + e−σ2(α′2+2m′2)/2

I0σ2

2(α′2+ 2m′2)+∞

k=1 m′2(σ2−1)

2k2−m′2Ikσ2

2(α′2+ 2m′2)!.

(114)

In the parts of (101) and (109) without the Dirac delta, one can sum over m′and has

K1=−m′

ωi(α−α′)i√2πσ2F

Km′(α′)(e2π(α−α′)−1)I0(σ2

2√(α2+2m′2)(α′2+2m′2)) e−σ2

4(α2+α′2+4m′2)

(115)

from (101) and

K2=i√2πω−i(α−α′)F

Km′(α′)(e2π(α−α′)−1)e−σ2(α2+α′2+4m′2)/4∞

k=1

Ik(σ2

2√(α2+2m′2)(α′2+2m′2))

(√2k−m′σ2)Λke(i(α−α′)−1) ln(m′−√2k)−(√2k+m′σ2)Λ−ke(i(α−α′)−1) ln(m′+√2k)

(116)

from (109).

Collecting (112), (113), and (114) and substituting (23) and (56), one is left with

F=√1−2σ21−1/σ2

2σ2eσ2(α′2+2m′2)/2+1 Γ1

4+1

2m′+1

2iα′

Γ3

4+1

2m′+1

2iα′

4−σ2

2I0σ2

2(α′2+ 2m′2)+∞

k=1 m′2(σ2−1)

2k2−m′2Ikσ2

2(α′2+ 2m′2).

(117)

Simplifying the Kernel 35

-30

-20

-10

-30

-20

-10

α’

-2

-1.5

-1

-0.5

0.5

(a)

-2

-1.5

-1

-0.5

0.5

-60 -40 -20 0 20 40 60

K(α−d, α +d)

(b)

-2

-1.5

-1

-0.5

0.5

1.5

-100 -80 -60 -40 -20 0 20

K(α, −α)

()

Fig. 4: The kernel of the self-consistent disk. (a) shows the real part of the kernel K(α, α′)

in (111) of an ω= 0.5,σ= 0.2, and m= 2-disk; (b) is the real part of the same kernel

along the main diagonal (black line) and along parallels K(α−d, α′+d)with d= 0.5(upper

dotted line) and d= 1 (lower dotted line); (c) is the same kernel along the secondary diagonal

α=−α′, this time with the imaginary part drawn as a dashed line. This is the “spur” visible

in (a).

36 Simplifying the Kernel

The kernel itself computes to

K=−iω−i(α−α′)√1−2σ21−1/σ2

2σ2eσ2(α2+α′2+4m′2)/4+1e2π(α−α′)−1Γ1

4+1

2m+1

2iα′

Γ3

4+1

2m+1

2iα′

m′i(α−α′)σ2I0(σ2

2√(α2+2m′2)(α′2+2m′2))−∞

k=1

Ik(σ2

2√(α2+2m′2)(α′2+2m′2))

(√2k−m′σ2)Λke(i(α−α′)−1) ln(m′−√2k)−(√2k+m′σ2)Λ−ke(i(α−α′)−1) ln(m′+√2k)!,

(118)

where, from (83),

Λ = αα′+√2m′(√2m′+ i α−α′)

p(2m′2+α2)(2m′2+α′2).(119)

Interestingly, the only term in (118) dependent on ω, i.e., ωi(α−α′), can be factored out of

the kernel, and Fis entirely independent of ω. Under these circumstances, the value of ωis

completely irrelevant for the solvability of the integral equation. This means that as long as

this kernel is valid, one will have either no solution at all or a two-dimensional continuum of

solutions (corresponding to the real and imaginary parts of ω). Of course, the case ω= 0

handled below is an exception here, since the kernel (118) is not valid for ω= 0. Actually, by

considering the physics, this result is not too surprising and is required by the scale-freeness

of the disk. Without a time scale to ‘cling to’ a disk has no way to know what growth rate

or pattern speed it should prefer. However, see chapter 6.1 for a more thorough discussion.

The kernel Klinks the input spectrum of a perturbation with respect to the potential basis,

Φm′(α), to the output spectrum, Φm′(α′), which is why Zang called it the “transfer function”.

This interpretation is useful when one wants to attain some qualitative understanding of the

various features of K. In the case of the self-consistent disk, this interpretation is not literally

correct, since the coefficient of Φm′(α′) depends on α′. Still, to gain some qualitative insight a

consideration of the kernel is valuable, and in the cut-out case the kernel does indeed directly

link the input and the output spectrum of the perturbation. Since the real parts of the kernels

are not qualitatively different for the cut-out and self-consistent disks (at some distance from

the diagonal, this statement is true for the imaginary part as well), much of what is said here

helps in understanding the cut-out disks.

By just looking at a picture of Klike the one in Fig. 4a, one immediately notices the strong

asymmetry in the kernel: Where α > α′, the kernel practically vanishes, except very close to

the diagonal, whereas on the other side of the diagonal, quite a few features are visible. The

interpretation is evident when one recalls that for ω > 0 a spiral is more trailing when αis

larger. Now, the kernel for α > α′gives the power flow from more trailing to more leading

spirals—which is minute. For this reason, this asymmetry is known as the trailing bias.

Going from α > α′towards the diagonal, one next encounters an extended contribution

streching out along the α=α′line. In the real part of the kernel shown in Fig. 4a, this

diagonal is rather inconspicuous; still, the large extent of the diagonal makes it one of the

dominant features of the kernel. Kernels for axisymmetric (m= 0) perturbations basically

Simplifying the Kernel 37

consist only of the diagonal. The imaginary part of the self-consistent kernel, on the other

hand, has a singularity on the diagonal which severely influences the entire character of the

eigenvalue problem (111). The diagonal roughly gives the feedback of a given logarithmic

spiral onto itself.

-30

-20

-10

-30

-20

-10

30 α’

-0.5

0.5

1.5

(a) -30

-20

-10

-30

-20

-10

α’

-1

(b)

Fig. 6: The kernel of the self-consistent disk for (a) m= 2,σ= 0.12 and (b) m= 1,

σ= 0.2. Again, ω= 0.5.

Right behind the diagonal on its trailing side, one or a few major peaks dominate the visual

appearance. A comparison between Fig. 4a and Fig. 6a suggests that these peaks widen

when the velocity dispersion decreases, and comparing Fig. 4a and Fig. 6b, one suspects that

lower mcorrespond to more concentrated peaks. A closer investigation reveals that this is a

definite trend.

In Fig. 6b one sees that while the m= 1-kernel is of course still symmetric on the diagonal, it

quickly becomes almost antisymmetric with respect to the secondary diagonal, which sets it

apart from the approximate symmetry of the m= 2 kernels (and these for the other azimuthal

harmonics). This, and a quite pronounced peak on the diagonal, is special to the one-armed

case and does not occur for other azimuthal harmonics.

The peaks are the main means of transport of power between spirals of different wave number

α. The first peak on the secondary diagonal transports power from leading to trailing waves

around α= 0, which suggests a relation to the well-known swing amplifier, in that the latter

is most effective when the wave pattern is almost radial.

Finally, a “spur” continues down the secondary diagonal, with small ripples on top of an os-

cillation that is only slowly damped. This spur, reaching farther out for smaller σ, transports

power from strongly leading to strongly trailing spiral components and might be related to

either reflection at the centre or possibly waves travelling through the centre along the lines

of the scheme hinted at in Binney and Tremaine (1987), Fig. 6-20. However, compared to this

scheme, an explanation of the spur would require the signs to be reversed, and of course it is

not obvious how these schemes should work in the presence of an inner Lindblad resonance,

when the centre becomes unreachable for density waves. Still, the fact that the coupling

mediated by the spur prefers waves of comparable wave numbers but different sign strongly

suggests some mechanism along these lines.

38 Simplifying the Kernel

Also note that in the kernel, there is little sign of a feedback loop that would transport power

from trailing to leading waves. There is, however, a nonzero contribution on the trailing-to-

leading side of the diagonal close to it. Taking a look forward, it should be noted that the

kernels of the cut-out disk do indeed have a small peak around the secondary diagonal on the

leading side of the main diagonal. It is tempting to connect this feature with the reflection

at the Q-barrier, even more so since it is the only qualitative difference between the cut-out

and the self-consistent kernels except for the singularity of the latter on the main diagonal.

3.6.1. The ase

ω= 0

As pointed out above, the case ω= 0 requires a special treatment. Fortunately, the expres-

sions become much more compact in this case.

Adding together (86) and (88), one has

K=√1−2σ21−1/σ2

√2σ2Γ1

4+1

2m+1

2iα′

Γ3

4+1

2m+1

2iα′

2 1−(1 + σ2)e−σ2(α′2+2m′2)/2

I0σ2

2(α′2+ 2m′2)+∞

k=1

2m′2

m′2−2k2Ikσ2

2(α′2+ 2m′2)!.

(120)

Remarkably, the Dirac and Kronecker deltas in the kernel have collapsed the integral operator

Φ(α)→RK(α, α′)Φ(α)dα to a simple multiplicative operator, and the integral equation

becomes an ordinary nonlinear equation,

Φm′(α′) = K(α′)Φm′(α′),or K(α′)−1 = 0.(121)

The solutions of (121) will always be marginal modes, since ω= 0 of course implies a vanishing

growth rate. To find these modes, one has to solve (121) as a nonlinear equation in the two

variables α′and σ. Since the kernel is quickly evaluated and I will not mass-produce such

results, bisection works well as a solution finder.

The Cut-Out-Disk 39

4. The Cut-Out-Disk

4.1. The Cut-Out Distribution Function

Two unphysical properties of the Mestel disk are the singularity of the density in the centre

of the disk and its infinite total mass. Those features are responsible for the “border terms”

in the Lzintegration above, and they cause the singularity on the diagonal. What is more,

the perfectly self-similar Mestel disk has no length or time scales, which leads to the weird

phenomenon that either a continuum of modes should be unstable or the disk is entirely

stable, with the possible exception of neutral modes with ω= 0 that do not need a time

scale.

A catch-all for those problems is the introduction of a cut-out. The idea, already promoted

by Zang (1976) is that a part of the matter in the disk is immobilised, i.e., it is only visible

in the rotation curve but not in the distribution function used in the Boltzmann equation.

Physically, this could mean that a stellar component with high velocity dispersion is present—

like bulges near the centre of a spiral galaxy—, or that a gravitationally unresponsive dark

matter halo strongly contributes to the rotation curve. Both possibilities appear to occur in

nature, and so immobilising a part of the mass is not an entirely unphysical procedure.

Mathematically, the immobilisation is effected by multiplying the distribution function with

another function H(Lz) with the properties that 0 ≤H(Lz)≤1 and preferably H(Lz) = 0

for Lz→0,∞strongly enough to ensure that the contributions proportional to δ(α−α′)

in I(Lz)

1,2are suppressed. Additionally, one drops the requirement of having a self-consistent

model, i.e., the potential remains the same as before. Thus, while some of the matter remains

dynamically important, it no longer takes part in a perturbation.

In the choice of the cut-out I do not strive to mimic any realistic models of haloes, disks,

or bulges. My primary concern is analytic simplicity of the resulting expressions. A cut-out

that leads to relatively compact expressions would be

H(Lz) = 2LzL0

0+L2

.(122)

with a positive L0. While I would recommend this functions for further analytic exploration

of a cut-out Mestel disk, my primary motive here is to check my approximation by comparing

my results to those of Read (1997), and therefore I will follow her choice of cut-out functions.

To damp down both harmonic tails in (90) and (103), I need to immobilise matter both at

the centre and in the outer regions. Therefore, I will use Read’s doubly cut-out function

H(Lz) = LN

zLM

(LN

z+LN

0)(LM

z+LM

c).(123)

40 The Cut-Out-Disk

0.2

0.4

0.6

0.8

1.2

1.4

0 2 4 6 8 10

f0(Lz, J)

J= 0, self-consistent

J= 0 cut-out

J= 0.01, self-consistent

J= 0.01 cut-out

Fig. 7: Effect of the cut out in a

σ= 0.2disk for two values of J.

Note that the distribution func-

tion of the self-consistent disk re-

ally diverges as J= 0 and Lz→

Since Mand Nwill occur as the orders of poles and the computation of high-order residues

leads to really messy expressions, I will only investigate the M=N= 2 case here. I do not

expect the dependence of the behaviour of the disk on Mand Nto be very different from

what Read found.

4.2. The Boltzmann Equation for the Cut-Out Disk

The modified distribution function does only introduce changes into the right hand side of

the linearised Boltzmann equation (29). The equivalent of (31) now reads

r.h.s. = f0 σ2Lz−√2J

σ2L2

+ 2 L4

z−L2

0L2

Lz(L2

z+L2

0)(L2

z+L2

c)∂Φ1

∂w1

+√2

σ2Lz

∂Φ1

∂w2!.(124)

Again, the term with Jis of order ζ2and is therefore discarded. Inserting the linearised

potential basis (25) now yields

r.h.s.=−i√2 Φα,mFLiα−1

σ2(L2

z+L2

0)(L2

z+L2

c)eimw2+ip√2J

Lz(αcos(w1)−√2msin(w1))−√2J

σ2Lz−iωt

q√2Jαsin(w1) + √2mcos(w1)−σ2m√Lz

√21 + 2L4

z−2L2

0L2

(L2

z+L2

0)(L2

z+L2

c)!

(125)

The Cut-Out-Disk 41

4.3. The Matrix Equation for the Cut-Out Disk

Just as the cut-out influences Read’s equations only via the angular momentum function, the

changes from (32) to (125) only involve Lz. Thus, all derivations up to the integration over Lz

may be carried out for the cut-out disk following the path outlined for the self-consistent disk

above. The two contributions to the kernel of the integral equation describing asymptotic

modes in the cut-out disk are, with no quadratures yet carried out,

I1=Φα,m

FLi(α−β)−2

4π2σ2Km′(α′)

zL2

(L2

z+L2

0)(L2

z+L2

exp 

is√2J

Lz(α−α′) cos(w1)−√2(m−m′) sin(w1)



exp i(m−m′)w2−√2J

σ2Lz!

(126)

and

I2=iΦα,m

FLi(α−β)−2

zmσ2+√2η

8π2σ2Km′(α′)(e2πiη−1)

zL2

(L2

z+L2

0)(L2

z+L2

1 + 2mσ2(L4

z−L2

0L2

(mσ2+√2η)(L2

z+L2

0)(L2

z+L2

c)exp i(m−m′)w2+ iηw′

1−√2J

σ2Lz

+is√2J

Lzαcos(w′

1+w1)−α′cos(w1)−√2msin(w′

1+w1) + √2m′sin(w1)!.

(127)

The conventions stated for the notation of (64) and (65) apply here, too, with the difference

that I denote the integrands for the cut-out disks with ¯

I1,2.

To highlight the terms that are added in going from the equivalent expressions for the self-

consistent disk, (64) and (65), the contributions already present in the self-consistent disk

are grayed out.

The black factors in ¯

I1,2just pull through all the quadratures I performed in the self-consistent

disk up until the Lzintegration. There, matters get significantly more complicated than

before.

42 The Cut-Out-Disk

4.4. Integration Over the Angular Momentum in the

Cut-Out Disk

4.4.1. Integration of

I(J)

The integrand of ¯

I(J)

1is the black term in (126) times the Li(α−α′)−1

zfrom the self-consistent

disk, which on transferring an L2

cto the rest of the term, leaves

Z∞

Li(α−α′)+1

(L2

z+L2

0)(L2

z+L2

c)dLz(128)

as the integral to evaluate. On substituting u= ln(Lz), one gets

Z∞

−∞

eiuγ+2

e2u+L2

02e2u+L2

c2du, (129)

where again I abbreviate γ=α−α′.

3π

5π

7π

−Rln(L0) ln(Lc)R

Fig. 8: Integration path for

the first integration of ¯

I(J)

1. The

crosses mark the poles relevant

for the integration.

The integrand has poles of order 2 at

¯u(1)

j= ln(L0) + iπ(j+ 1/2)

¯u(2)

j= ln(Lc) + iπ(j+ 1/2)

Assuming γ > 0 again, one can use the integration path sketched in Fig. 8. Of the three

problems noted with the integration of (91), only one is present here: Going towards +i∞, the

poles of the integrand go all the way up, so that a continuous deformation of the integration

path as we go to infinity is not possible. However, the remarks given above about why the

residue theorem may still be applied hold for the present integration as well.

The Cut-Out-Disk 43

Otherwise, the integrand is much tamer than the one I considered in the self-consistent case,

since it falls off exponentially as u→ ±∞. As if to make good for this, the residues are

substantially harder to compute. For poles of order two one has

res(f(z), u) = lim

u→z

duf(u)(z−u)2;

The determination of the limit requires a triple application of de l’Hospital’s rule. The results

are

res eu(iγ+2)

e2u+L2

02e2u+L2

c2; ¯u(1)

j!=−Liγ

e−(πγ/2)(2j+1)

2(L2

0−L2

c)(130)

for ¯u(1)

jand

res eu(iγ+2)

e2u+L2

02e2u+L2

c2; ¯u(2)

j!=Liγ

e−(πγ/2)(2j+1)

2(L2

0−L2

c)(131)

for ¯u(2)

The two residues (130) and (131) must be summed up for j= 0 ···∞, which results in an

alternating geometric series of the form Pj(−1)jqj. Such a series is easily determined to be

∞

j=0

(−1)jqj=∞

j=0

q2j−∞

j=0

q2j+1 =1−q

1−q2=1

1 + q.(132)

Using this to sum up (130) and (131), applying the residue theorem and inserting the resulting

expression into (128) delivers

I(Lz)

1=−IπΦ(α, m′)√2FL2

4Km′(α′)(L2

0−L2

Li(α−α′)

0−Li(α−α′)

sinh(π/2)(α−α′)e−σ2

4(α−α′)2.(133)

Here, I have already performed the summation over m(cf. Eq. 67). This is even more

advisable in ¯

Iw′

2, since the expressions are complicated enough even without keeping mand

m′. Note that (133) is regular at α=α′despite the sinhπ(α−α′)/2in the denominator,

because the Li(α−α′)

0−Li(α−α′)

cin the numerator also has a zero of order 1 at α=α′.

4.4.2. Integration of

I(w′

, Part Without Sum

As I did in the self-consistent case, I split up the integration of I(w′

2into two parts. At

first I consider the terms without the series. After again substituting u= ln(Lz) and writing

γ=α−α′, one has to compute

Z∞

−∞

e(iγ+2)u

(m−euω)(L2

0+ e2u)2(L2

c+ e2u)2mL2

0L2

c(1 + σ2)−ωL2

0L2

ceu

−m(σ2−1)(L2

0+L2

c)e2u−ω(L2

0+L2

c)e3u+m(1 −3σ2)e4u−ωe5udu.

(134)

44 The Cut-Out-Disk

2π

3π

−Rln m

ω

ln(L0) ln(Lc)R

Fig. 9: Integration path for the

part without the sum in the in-

tegration the first part of ¯

I(w′

with γ > 0. The crosses mark

the poles relevant for the integra-

tion.

This problem is essentially solved when one knows the integral

Z∞

−∞

eiγueλu

(m−euω)(L2

0+ e2u)2(L2

c+ e2u)2du (135)

for λ∈ {2,3,4,5,6,7}. This integrand has poles at ¯u(1)

jand ¯u(2)

j, and, in addition to that, at

¯u(3)

j= 2πij+ ln(m/ω),(136)

where the poles at ¯u(1,2)

jare of order two, the poles at ¯u(3)

jare of order one. The situation in

this case is sketched in Fig. 9.

At this point, the computations really get messy for Read’s cut-out, and it seems wise to take

recourse to computer algebra systems. In the appendix, I give the results from this and the

following quadratures as exported by Maple. Apart from the sheer size of the expressions,

the calculations can be carried out in close parallel to those in the case without the sum and

those in the singular disk.

The result of having still rather clumsy expressions does not seem to warrant the effort

of further cleaning up these formulae. As mentioned above, if goal is a cut-out disk with

relatively compact expressions, I recommend using a cut-out of the type (122).

4.4.3. Integration of

I(w′

, Part With Sum

After the usual substitutions, the terms containing the Lzintegration of the part of ¯

I(w′

with the series can be computed evaluating

Z∞

e(iγ+2)u2ksin(kln Λ) −i√2 cos(kln Λ)m+ωeu

(ωeu−m)2−2k2e2u+L2

02e2u+L2

c2mL2

0L2

c(1 + σ2)−ωL2

0L2

ceu

−m(σ2−1)(L2

0+L2

c)e2u−ω(L2

0+L2

c)e3u+m(1 −3σ2)e4u−ωe5u(137)

The Cut-Out-Disk 45

Now the reduced problem is finding

Z∞

e(iγ+λ)u

(ωeu−m)2−k2L2

0+ e2u2L2

c+ e2u2du (138)

for λ∈ {2,3,4,5,6,7,8}. The zeroes the denominator of this integrand are at ¯u(1)

j, ¯u(2)

j, and

¯u(4)

j= ln m−√2k

ω!+ 2πij

¯u(5)

j= ln m+√2k

ω!+ 2πij

.(139)

The pole at u(4)

jplays the role of the pole at u(3)

jin the Iw′

2integration in the self-consistent

disk, i.e., when √2m < k they are shifted vertically, and again a proper choice of the loga-

rithm’s branch cut allows the use of identical expressions regardless of the ratio between m

and k.

2π

3π

ln(L0) ln(Lc)

−Rln √2m+k

ω

ln √2m−k

ω

Fig. 10: Integration path for

the part with the sum in the in-

tegration of ¯

I(w′

2with γ > 0and

k < √2m. The crosses mark the

poles relevant for the integration.

The situation found here is depicted in Fig. 10. Basically, the integral can be computed as

the sum of four series over residues at poles of orders one (u(3)

j,u(4)

j) and two (u(1)

j,u(2)

j). But

while the expressions for the residues at u(3)

jand u(4)

jare quite like those in the self-consistent

disk with the denominator adapted, those for the residues at u(1)

jand u(2)

jare really messy,

and I once more refer the reader to Appendix A.

Suffice it to say that the kernel of the integral equation for the cut-out disk is quite longish

to write down but much easier to handle than the self-consistent disk’s kernel: it is bounded

at the diagonal α=α′and it falls off rather quickly with increasing distance to the diagonal.

So, the cut-out disk it is much more accessible to a numerical study of its behaviour.

46 The Cut-Out-Disk

4.5. Numerics

After the final integration over Lzthe integral equation for modes in the cut-out disk has the

form

Φα′,m′=Z∞

−∞

K(α, α′)Φα,m′dα with

K(α, α′) = ¯

I(Lz)

1+¯

I(Lz,1)

2+¯

I(Lz,2)

2(140)

where ¯

I(Lz)

1is given by (133) and ¯

I(Lz,1,2)

2can be found in the appendix. Formally, this is

a homogenous Fredholm integral equation of the second kind. However, due to the infinite

integration interval, the operator induced by ¯

Kis not compact, and therefore the well-known

results from the theory of integral equations (e.g., Pipkin 1991), and in particular the Fred-

holm Alternative, are not immediately applicable. Thanks to the exponential decay of the

kernel, however, the operator is “compact for all practical purposes”, meaning that as long

the functions the operator is applied to do not diverge too quickly as α→ ∞, the Fredholm

Alternative still holds. Indeed, functions that do not vanish as α→ ∞ can probably not be

considered solutions to the given problem since they do not satisfy the condition that the

perturbation of the surface density is much smaller than the surface density of the basic state

unless they were finely crafted to do so.

This “almost-compactness” of ¯

Kwas already noted by Zang, who simply cut off the inte-

gration region when increasing it did not change the solutions noticeably any more. I will

adopt this scheme for purposes of the numerical investigation of the solutions. As for closed

solutions, the integral equation (140) is far too complex to undertake an attempt to find

them, even if there was any hope that a kernel that complex might have any.

The kernel of ¯

Kis quite easy to compute, since the series usually converges rather quickly and

high-ksummands are small. However, the diagonal elements of the kernel are evaluated with

poor accuracy, even diverging in the imaginary as α→α′at the limited accuracy provided

by the C library or the numerical coprocessor. The problems surface when |α−α′|<

∼0.05

for IEEE double precision floats and parameters in the interesting range (say, 0.05 < σ < 0.5

and 0.001 < ω < 10).

The reason for this behaviour is of course the presence of the sinh−1term in (133) and the

corresponding terms for ¯

I2. Although it is balanced by zeroes of order 1 in the numerators,

rounding errors invalidate the results from the computation. To overcome this difficulty, one

could rewrite the kernel isolating the zeroes in both numerator and denominator and compute

them separately. With this approach the rounding errors are more confined, and indeed the

GNU C library on iA86 and PPC architectures handles terms like x/ sinh(x) quite well down

to very small x. Another approach would be to use a suitable expansion of the kernel in the

vicinity of the diagonal.

Both these approaches are unattractive because of the sheer size of the expressions. Therefore,

I resorted to polynomial interpolation to compute the kernel for α≈α′. When a kernel value

K(α, α′) is requested with |α−α′|<0.05, I compute the four points ¯

K(α′−0.18 + 0.12 j, α′),

j= 0 ···3 and determine ¯

K(α, α′) from the third-order polynomial fitting through these four

points. For the values of ω,σ,L0, and Lcused in this work, the characteristic periods in the

kernel are generally larger than 0.1, and thus this approach works very well.

The Cut-Out-Disk 47

Following Read and Zang, the integral equation itself is solved by re-writing it as an eigenvalue

equation,

λΦm′(α′) = Z¯

K(α, α′)Φm′(α)dα. (141)

The complex numbers λfor which a Φα,m exists so that this equation can be satisfied are

called mathematical (or artificial) eigenvalues. Although those eigenvalues have no physical

significance in themselves, they will allow to judge on the solvability of the integral equation

and have the advantage of being quite easily computed (see chapter 5.2). To find the Φα,m′

that really are modes, one has to look for parameters for which ¯

Khas an eigenvalue 1.

4.5.1. Quadrature

To compute the eigenvalues of ¯

K, I employ the method of Fredholm discretisation (also known

as Nyst¨om method): Given a quadrature formula Rf(x)dx =PN

i=1 wif(xi) with equidistant

grid points xione can approximate the integral equation (141) with the system of linear

equations

λΦ(αj, m′) =

i=1

wi¯

K(αi, αj)Φ(αi, m′).(141)

Since the kernel oscillates rather quickly where it is large, a large number of grid points

would be necessary to ensure useful accuracy of a simple quadrature formula without special

adaption to the problem at hand. Thus, using an advanced quadrature formula is highly

desirable.

In general applications of numerical integration, the most popular advanced quadrature for-

mulae belong to the family of Gauss quadratures. These give quadratures exact for a weight

function times a polynomial, but in general require an unevenly spaced grid. When solving

integral equations, this feature makes them quite useless, since one has to derive a system

of linear equations, and thus the discretisation has to be the same in αiand αj. There

are, however, adaptions of the Gauss quadrature on a grid (the most common of which is

the Gauss-Legendre quadrature) and related methods to compute quadrature formulae with

weight functions on an evenly spaced grid (e.g., Press et al 1992, p. 798).

Unfortunately, choosing a weight function for the present kernel is highly nontrivial, because

the three periodic functions present in the kernel (with periods 2π/ω, 2π/L0, and 2π/Lc) are

overlaid in a rather complicated way, and quasiperiodic terms further worsen the situation. I

have tried implementing Press’ prescription mentioned above with the dominant oscillation

(ωi(α−α′)) as a weight function, which did not yield a very satisfying result. This did not really

come as a surprise, since factoring out ωi(α−α′)does not smoothen the kernel too significantly.

Trying more complicated weight functions did not improve the situation.

A relatively proven tool to handle the discretisation of the integral equation at hand is the

quadrature already used by Read (1997) and Zang (1976), which is based on an approximation

of the kernel by Lagrangian polynomials. However, I found this scheme to be effective only

when the kernel already is rather smooth. and it did not work too well for me. It is somewhat

hard to understand why it worked that much better for Read, who only factored out the

Kalnajs factor.

48 The Cut-Out-Disk

I finally settled for a quadrature based on local spline interpolation. In this scheme, one first

determines the kernel on a grid of grid spacing ∆α, yielding N2numbers ¯

Ki,j. One then

one puts a cubic spline s(α) through all the points ( ¯

Kij)i=1···N, i.e. one row of the resulting

matrix. This spline usually is a good representation of the true kernel unless the kernel is

grossly undersampled. Hence an integral over this spline is a reasonable approximation to

the integral RαN

α1

K(α, αj)dα. This integral could in principle be computed analytically, but

for convenience I evaluate in by the trapezoidal rule. Even with 30 abscissae per ∆αinterval,

the computing time is still dominated by the evaluation of the kernel.

In terms of the matrix ¯

Kij, this scheme invites setting ¯

Kij =Rαi+h/2

αi−h/2s(α)dα, where at α1

and αNone only integrates from α1and to αN, respectively. Strict numerics would tell one

that this tapering with the first and last step does not mak much sense—in practice, it does.

After this procedure, one transposes ¯

Kij and repeats it. Since the eigenvalues of a matrix

and its transpose are identical, one does not need to undo the transposition.

4.5.2. Auray

The accuracy of the eigenvalues obtained by this discretisation depends on a number of

parameters. Given that the purpose of my numerical calculations mainly is checking whether

or not the approximations that lead the kernel in (140) deliver results comparable to those

of Read, an accuracy of a few percent in the eigenvalues will suffice. What parameters are

needed to attain this accuracy?

Firstly, the accuracy with which the kernel can be computed. Apart from the problems near

the diagonal mentioned above, this is not a very serious issue. To check this, I have evaluated

the expressions for the kernel in Maple with various settings of Maple’s digits parameter and

found that the expressions are not particularly prone to round-off errors. The values resulting

from the C++ programme that delivered the results were in every case correct to at least

five digits, in most cases correct to four digits less than machine accuracy.

Secondly, the numerical condition of matrix with respect to Schur’s algorithm that was used

to compute the eigenvalues. If this problem was ill-conditioned, one would expect slightly

different discretisation sizes to lead to significantly different results. This is not the case;

for example, computing the largest mathematical eigenvalue for a m= 2 perturbation with

ω= 0.5 in a σ= 0.2 disk, one finds that with N= 250, the eigenvalue is 0.8932459−1.253385i,

whereas for N= 251, one finds 0.8932604 −1.253254i.

Finally, the parameters of the discretisation, Nand L, influence the accuracy of the eigenval-

ues. For the ω= 0.5, σ= 0.2 disk, one finds that the largest mathematical eigenvalues for an

m= 2 perturbation are 0.892403 −1.26264i when N= 150 and L= 30, whereas for N= 300

and L= 60 it is 0.892419 −1.26247i. So, for this moderately high σ,L= 30 would easily

suffice for precision in the range of one per cent. As σdecreases, however, the Lrequired

to attain a given precision increases, and thus I use L= 40 in the calculations below unless

otherwise stated.

The accuracy of the solution depends on Nthrough the grid point spacing L/(N−1). Having

N= 200 at L= 30, the largest mathematical eigenvalue is 0.89342−1.251445, with N= 300,

it is 0.89413 −1.243457, and with N= 400, it is 0.89437 −1.240663. It turns out that this

The Cut-Out-Disk 49

result is fairly generic, in that neither σnor ωsignificantly influence the dependence of the

eigenvalues on L/(N−1). Thus, with L= 40, N= 250 is a good choice.

50 The Classic Modes

5. The Classi Mo des

The term “classic modes” here refers to the modes already well investigated by Read and,

before her, Zang. The main purpose of this chapter is to show that in my approximation

all the major features of the exact calculations show up, and that the quantitative results

obtained from the two methods usually agree to within a few per cent.

5.1. Axisymmetric Stability

This chapter investigates the most fundamental stability property a disk can have, the sta-

bility against ring-like perturbations. In the potential basis (24), ring-like perturbations

corresponds to setting m= 0. After this, 2π/ ln(α) is the distance between two rings in the

perturbation.

Axisymmetric stability is a fundamental stability property not only because it was the first

to be closely investigated (Toomre 1964), but also because it is the equivalent of the Jeans

instability in rotating disks. A long tradition of theoretical evidence and numerical exper-

iments indicates that a disk susceptible to axisymmetric instabilities quickly fragmentates

(e.g., Hockney and Hohl 1969).

To have a quick assessment of the stability of a disk, one frequently gives its Toomre-Q, defined

as the ratio of its velocity dispersion to the velocity dispersion required for axisymmetric

stability, Q=σ/σmin. The denominator in this definition is usually computed from the local

analysis Toomre gave in 1964. As Zang showed, the local and global analyses are in excellent

agreement for the Mestel disk. Whereas in real galactic disks Qusually varies dramatically

with the position within the disks, the self-consistent Mestel disk has one global Q.

5.1.1. Axisymmetri Stability of the Self-Consistent Disk

As Read remarks, m= 0 perturbations are not orientable and can therefore not rotate. Thus

one may set ω= 0 (or use a purely imaginary ωif one wants to investigate growing modes).

However, in letting m′= 0 some care is required as once more, and one has to set m=m′= 0

in (78) and (86) to get the limit right—just setting m′= 0 in (110) leads to erroneous results.

The two integrations over Lzare as straightforward as in the general ω= 0 case and yield

K(m=0) =√1−2σ21−1/σ2

2eσ2|Γ(1/4 + iα′/2)|2

|Γ(3/4 + iα′/2)|21−e−σ2α′2/2I0(σ2α′2/2).(143)

To find modes, one has to look for (σ, α′) pairs that solve K(σ, α′) = 1. The special functions

in (143) prevent an analytic solution of this problem, and so one has to take recourse to

The Classic Modes 51

0.2

0.4

0.6

0.8

1.2

1.4

1.6

1.8

0 2 4 6 8 10

K(α′)

α′

σ= 0.2

σ= 0.3

σ= 0.38597

σ= 0.5

Fig. 11: Plots of the kernel

(143) governing m= 0 modes in

the self-consistent Mestel disk for

four values of σ.

numerics. The result is shown in Fig. 11. For cool disks, K(α′) crosses unity twice, which

means that there are two marginal modes in those disks. This will be of minor interest, since

between them there are growing modes which will dominate those marginal modes. The

axisymmetric modes set in at σ= 0.3860 with αcrit = 3.51. Below this threshold, the disks

are stable against axisymmetric perturbations.

This result is in remarkable agreement with the findings of Read. She found that σ= 0.3780

was required to stabilise the disk, and the critical wavelength is αcrit = 3.46. The figures

agree to within 2% in both σand αcrit, which in a way is surprising given that those disks

are already quite hot, and thus a large number of stars will be on orbits for which the

approximations in chapter 2.1 are very poor. On the other hand, this stability limit is fairly

generic and is already given correctly in a local analysis without assuming a specific disk

model. Therefore, while it is reassuring that it is reproduced, the agreement is not really an

indicator for the quality of the approximation.

I will use Q=σ/0.3860 in the rest of the work for consistency, although it does not make

much of a difference.

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

-40 -30 -20 -10 0 10 20 30-0.16

-0.14

-0.12

-0.1

-0.08

-0.06

-0.04

-0.02

0.02

K(α, α)

-0.06

-0.04

-0.02

0.02

0.04

0.06

-10 -5 0 5 10

K(α, −α)

Fig. 12: The kernel for axisym-

metric perturbations of the cut-

out disk (σ= 0.2). In the upper

figure, the real (solid line) and

the imaginary part (dotted, the

right scale applies) along the di-

agonal is shown. In this case,

the imaginary part is plotted a

little bit off the diagonal since it

vanishes exactly on the diagonal.

In the lower plot, the real (solid

line) and the imaginary part of

the kernel are shown along the

secondary axis.

52 The Classic Modes

5.1.2. Axisymmetri stability of the Cut-Out Disk

For determining the axisymmetric stability of the cut-out disks, the kernel derived for the

general case is again unsuitable since letting m= 0 in the finished expression for the kernel

delivers incorrect results. Instead, one has to back up to perform the Lzintegration with

ω= 0 and m= 0. Again, this setting greatly simplifies the calculations that finally yields

K(m′=0)(α, α′) =iLi(α−α′)

c−Li(α−α′)

0L2

c√1−2σ21+1/σ2

4σ2(L2

0−L2

c)(e−π(α−α′)−1) Γ(1 + 2iα′)/4

Γ(3 + 2iα′)/4

I0σ2αα′/2−eσ2αα′/2(144)

As Zang (1976) first noted (and one might suspect from the visual appearance of the cuts

through the kernel in Fig. 12), the kernel of the integral equation governing axisymmetric

perturbations of the Mestel disk can be brought to Hermitian form. The transformation that

does this is identical to the one used by Zang. The basic observation is that (144) would

be Hermitian but for the terms with the Gamma function. Looking back at (140), one sees

that by writing the integral equation for a new function ˜

Φα,m =pKm(α)Φα,m, one attains

Hermitian symmetry in the resulting kernel pK′

m(α)/K′

m(α′)¯

K(m′=0)(α, α′). Now, if the

kernel is equivalent to a Hermitian kernel, its eigenvalues will be real, and, as it turns out

when using the machinery developed in chapter 4.5 on this kernel, this is indeed the case.

0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5

λ1(σ)

Lc= 100

Lc= 10

Fig. 13: The largest mathemat-

ical eigenvalue of the m= 0 cut-

out disk as a function of σfor

two values of the outer cut-out

Lc.

After that, finding the limit for marginal stability is just a matter of determining for which

σthe largest mathematical eigenvalue is unity. For the doubly cut-out disk with Lc= 100

(I deviate from my standard Lc= 10 following Read in this) one finds σmin = 0.349 or

Q= 0.904, where Read has σmin = 0.327 or Q= 0.866. These figures agree to within about

5%, which is somewhat worse than the agreement in the self-similar disk. After the discussion

at the end of the last subsection this does not come as a surprise. In the case of the cut-out

disk the global structure of the disk plays a major role, and it is clear that at σ≈0.35 this

global structure is already quite distorted by the approximations. Still, the qualitative result

is that the cut-out slightly destabilises the disk, and also quantitatively the agreement is quite

satisfactory.

The Classic Modes 53

On reducing Lc, the disks become more stable. For example, at Lc= 10, σmin = 0.274. This

is easily explained by the reduced surface density of the disk.

5.2. One-Armed Modes

Both the works of Read and Zang found that the m= 1 modes are the most unstable in

Mestel disks. This came as a surprise at least to Zang, since the classic picture of a spiral

galaxy invariably involved a two-armed pattern. As a matter of fact, a marginally stable

(Q= 1) cut-out disk admits only one-armed growing modes (at least with my choice of the

index of the inner cut-out). In this chapter, I will show that the present formalism reproduces

this result. The more exotic features of the m= 1 perturbation in the Mestel disk are also

recovered.

5.2.1. Neutral Mo des in the Self-Consistent Disk

The m > 0 cases are more complicated than the axisymmetric case was, since now one has to

decide on the solvability of the full integral equation (111) with the kernel (118). Making this

decision is quite involved indeed and requires techniques different from those applied when

ω= 0 or with a cut-out. Furthermore, neither Read nor Zang could solve this problem, but

both were rather convinced from indirect arguments that growing modes should not exist in

the self-consistent disk.

For these reasons, I postpone pondering this issue to chapter 6 and for now stick to the “classic

modes” known to Read. In the self-consistent disk, those are the neutral modes, i.e. single

logarithmic spirals with ω= 0. Since a growth rate or a pattern speed would introduce

time or length scales which the Mestel disk does not provide, those neutral modes are the

only type of isolated modes that are possible in the self-consistent disk without breaking its

self-similarity.

0 5 10 15 20 25 30 35 40 45 50

K(α′)

α′

Fig. 14: The kernel for neutral m= 1 modes. I show curves for (top to bottom) σ= 0.2,

σ= 0.3,σ= 0.38597, and σ= 0.5.

54 The Classic Modes

To find modes, one has to solve (121) for a given σ. Again, this can only be done numerically

(or graphically). At least qualitatively, however, one can already say that as long as my

approximations are valid, the self-consistent disk will admit neutral m= 1 modes at any

temperature. Evaluating the kernel at α′= 0, one sees that at least for σ < 1/√2, the kernel

is larger than 1, whereas Ktends to zero for α′→ ∞. This latter asymptotics can be seen by

recalling that Iν(z)∼ez/√2πz for large z, which, combined with the factor exp(−σ2α′2/2)

and a 1/α′form the Gamma functions gives the kernel an asymptotics of α′−3/2. Thus, for

large α′, the kernel vanishes, and because it is a continuous function of α′, it must cross the

K= 1-line at some α′. Indeed, as Fig. 14 indicates, this seems to be the case. As σincreases,

the wave number of the mode decreases.

Of course, it does not seem likely that a disk admits modes at any temperature form an

astrophysical point of view. While one could maintain that at some point the linearised

analysis will certainly no longer be applicable and higher contributions would suppress the

modes, the behaviour found here is probably related to the fact that m= 1 perturbations can

shift the centre of mass of the disk and that density cusps in the centres of galaxies tend to be

unstable. This phenomenon is well known from numerical studies and was noticed as early as

1963 by von Hoerner as a nuisance and later re-discovered by Miller and Smith (1992) who

investigated it on its own right. Recently, Taga (1999) reported his simulations lead him to

believe that a black hole embedded in a dense stellar cluster would become unstable to to

displacements when its mass was below about 10% of the mass of the cluster; for the Mestel

disk, the central singularity is not associated with a compact object at all, and thus this 10%-

rule would certainly make the cusp eligible for the displacement instability. In the formalism

used here, the centre of mass is artificially kept fixed. Given these results, this is a serious

forcing, and it appears entirely possible that the instability of the centre cusp—that is not

dynamically present in the cut-out disk—resurfaces as the non-quenchable m= 1 instability.

On the other hand, all this might as well be interpreatble along the lines of what I discuss

in the context of the very similar behaviour is found for m= 4 in chapter 5.5. Possibly this

issue might be decided by the appearence of the true marginal mode (as the limit of growing

modes) arising from the equations laid out in chapter 6.1.

On a side note, the kernel (120) is symmetric in α′. This means that in addition to “trailing”

modes with positive α′(the disk stars enter the spiral arm from the convex side), there is

also a “leading” one (the stars enter the arm from the concave side). This could have been

anticipated since both modes are steady state solutions of a time-reversible physical system—

both the Poisson and the Boltzmann equation are ignorant of the time arrow—, and thus the

anti-spiral theorem (Lynden-Bell and Ostriker 1967) applies. The situation for the modes in

the cut-out disk is different, since by application of Landau’s rule in their construction they

are limits of growing modes and thus not steady-state solutions even when their growth rate

vanishes.

5.2.2. Mo des in the Cut-Out Disk

In contrast to the self-consistent disk, the cut-out disk has a characteristic length scale (L0,

say), and ωis present in many terms in the kernel. Thus, discrete and growing modes should

be possible. However, finding them is more difficult than in the axisymmetric case. The

The Classic Modes 55

0.2

0.4

0.6

0.8

0 0.5 1 1.5 2 2.5

Re(λ)

Im(λ)

(a) (b)

-0.2

0.2

0.4

0.6

0.8

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6

Re(λ)

Im(λ)

Im(ω) = 0.000

Im(ω) = 0.010

Im(ω) = 0.100

Im(ω) = 0.200

Im(ω) = 0.400

()

-0.6

-0.4

-0.2

0.2

0.4

0.6

0 0.5 1 1.5 2

Re(λ)

Im(λ)

σ= 0.150

σ= 0.200

σ= 0.400

σ= 0.600

(d)

Fig. 15: Eigenvalues of the cut-out m= 1 kernel. (a) shows spectra for kernels with ω=

0.001 ···0.6in 20 steps, where each symbol corresponds to one value of ω(for most of them,

actually two values, the symbols repeat after the filled circle); (b) is a plot of the largest math-

ematical eigenvalues of a Q= 1 disk versus the Re(ω)and growth rate s= Im(ω)for Read’s

calculations for the inner cut-out disk, taken from Read (1997) with kind permission of the au-

thor; (c) is a plot of the eigenvalues with the largest real part form the approximated kernel ver-

sus Re(ω) = (0.0001,0.001,0.005,0.01,0.05,0.1,0.15,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.2)

and the indicated values of Im(ω); (d) is the same plot for a fixed growth rate Im(ω) = 0.1

and the indicated values of the velocity dispersion σ.

56 The Classic Modes

symmetry of the kernel that resulted in purely real eigenvalues is broken for nonzero m, and

therefore the eigenvalues will in general be complex.

Pondering the kernel depicted in Fig. 17, one recognises the features discussed in chapter 3.6:

A contribution reaching out to large αon the diagonal, quite pronounced in this case, and

a spur along the secondary diagonal. The trailing bias discussed above is of course present

here, too. As the m= 1 kernel of the self-consistent disk (Fig. 6b), the kernel, though

still symmetric on the diagonal, quickly attains some distinctly antisymmetrical (w.r.t. the

secondary diagonal) component, but has a roughly symmetrical appearance farther away from

the diagonal.

Figure 15a shows the eigenvalue spectra of the m= 1 kernels for a few values of real ω(i.e.,

the plot only shows kernels describing marginal modes). The eye-catcher in this diagram is

the arc that spans from about 0.3 + 0.7i to 1.9 + 0i far above the bulk of the eigenvalues.

Those “rogue” eigenvalues do not show up for any azimuthal harmonic except m= 1.

Speculating for the moment, the physical significance of this feature can be expected to lie

in one of the two features special to m= 1 modes, namely either the fact that the one-armed

modes do not possess an inner Lindblad resonance, or that they can move the centre of

gravity in violation of my assumptions. It seems hard to decide this matter handwavingly.

Fortunately, this issue is not relevant for the basic stability properties of the disks because

the eigenvalue with the largest real part will always be in the bulk part of the spectrum.

However, in the plots of the largest eigenvalues I use the ones with the largest real part when

m= 1, whereas for the other azimuthal harmonics I use the ones with the largest moduli.

-30

-20

-10

-30

-20

-10

α’

-0.2

0.2

0.4

0.6

Fig. 17: The real part of the

kernel for m= 1 perturbations

of the cut-out disk. This kernel

is for a velocity dispersion σ=

0.3. The imaginary part is quite

similar.

This bulk part of the spectrum is quite simple: for each kernel, the eigenvalues are grouped

along a dented curve from the origin to some point in the upper half plane, where this point

lies the farther out the in the upper half plane the larger ωis. For very small ω, the spectrum

is almost aligned along the real axis. Apparently, no eigenvalues are located below the real

axis as long as Im(ω) = 0, but the curves approach the real axis as ωdecreases.

Most eigenvalues come in pairs, which, given the limited accuracy of the calculation, may well

indicate that each pair actually corresponds to exactly one eigenvalue with a two-dimensional

eigenspace. I will return to this issue when discussing the m= 2 modes.

The Classic Modes 57

Those eigenvalues as such are, however, only of limited interest, since all the physics is in

the eigenvalue λ= 1 + 0i; only the eigenvector(s) belonging to this eigenvalue correspond

to modes in the disk described by the kernel under scrutiny. The obvious question now is

whether these exist. Before plunging into the numerics, one should first approach the problem

from a point of view of the theory of complex functions. Reformulating (141), one has

λ(ω) = 1

Φm′α′(ω)Z¯

Kα, α′(ω)Φm′(α)dα. (145)

Solutions Φm′(α) that I am interested in will be depend at least continuously on ω, and

I require that they do not vanish except on countably many isolated points (with these

assumptions, Φm′(α) must even be an analytic function of ωbecause of its definition via

an integral). Now (145) can be read as the definition of a function λ(ω), and since ¯

Kis

an analytic function of ωin the sliced plane (taken as a complex number with a mode’s

growth rate in the imaginary part) and this property carries through the integration, λ(ω)

is analytic as well. This means that one can exploit the rich theory of conformal mappings

when exploring the mapping between ωand the mathematical eigenvalues.

Fig. 15c shows a part of this mapping by plotting the largest mathematical eigenvalues for

various complex ωin the (Re λ, Im λ) plane. Points of equal Im(ω) are connected by spline-

interpolated lines. The decisive point in this diagram is λ= 1+0i, since when a line crosses the

real axis there, the disk admits a mode with the corresponding growth rate. However, since

λ(ω) is (as good as) analytic, it depends on ωcontinuously, i.e., when the line corresponding

to a growth rate Im(ω1) crosses the real axis left of 1 + 0i and the another of growth rate

Im(ω2) right of it, the midpoint theorem states that there has to be a growth rate Im(ω0) in

between that for some Re(ω) has an eigenvalue 1 + 0i.

A similar reasoning works for σwith the difference that of course the kernel is not an analytic

function of ω+ iσ. Still, using that λ(σ) is continuous, one may again use the midpoint

theorem to state that when the eigenvalue line corresponding to various Re(ω) at a given

growth rate Im(ω) crosses the real axis right of unity and if for large velocity dispersions all

eigenvalues are left of the Re(λ) = 1 line, then there will be at least one σ0which allows a

unity eigenvalue for at least one ωwith the given growth rate as its imaginary part.

So, while the only point in the diagram relevant for the physics in the Figure 15b is 1 + 0i,

the mappings λ(ω) and λ(σ) do contain valuable information. The first thing to note is that

at the velocity dispersion at which this plot was made (σ= 0.3860, corresponding to Q= 1),

there should indeed be a growth rate for which the kernel has a unity eigenvalue, and thus,

it does admit a growing mode.

This issue is much more difficult for the marginal (Im(ω) = 0) mode. From the graph, it

is hard to tell whether or not it eventually reaches (or crosses) the real axis. In the above

discussion of the spectrum of the kernel governing marginal modes (Fig. 15a), it was already

pointed out that a very accurate calculation would probably show that all eigenvalues have

positive imaginary parts. But even when using 800 discretisation points on a grid of 120×120

in αand α′, the largest eigenvalue of a Q= 1 kernel is at 1.7485 −4.053 ×10−9—certainly

attempting a decision on the positions of the true eigenvalues with respect to the real axis

are very daring given the uncertainties involved with computing eigenvalues of matrices of

that size.

58 The Classic Modes

A related question is whether the slow march of the ω=ǫ+ 0i eigenvalue along the real axis

as ǫ→0 apparent in Fig. 15c continues or stops. It appears that it does stop at ω= 1.47;

at least going from ω= 10−6to ω= 10−8does not change the leading five figures of the

eigenvalue 1.4753 −7.4499 ×10−7when evaluating the kernel on a N= 200 grid.

Comparing Read’s results in Fig. 15b (for an inner cut-out disks, so one has to expect her

disks to be somewhat more stable) and my results in Fig. 15c, one immediately notices that

although the curves roughly match, there are some differences in details. Read’s curves

approach the real axis somewhat faster, and they reach out to far larger values of Re(λ) as

ω→0. While the first disagreement is not very worrisome, particularly because for nonzero

growth rate there is as good an agreement in the points where the λ(ω) curves cross the real

axis as can be expected given that Read’s disk only has an inner cut-out, the second one

is serious since it will lead to different stability predictions in the limit of vanishing growth

rate; while it appears uncertain if cut-out Mestel disks could ever be stabilised against m= 1

perturbations from Fig. 15b and indeed Read does not even attempt to find a σthat would

do this, Fig. 15c nourishes hope that the disks might eventually become stable as one raises

the velocity dispersion.

In fact, this is a point in which the doubly cut-out disk behaves significantly different from

the inner-cut out disk. Increasing the outer cut-out angular momentum Lcto 100 raises the

largest mathematical eigenvalue to 1.889 for ω= 0.0001 and to 1.890 for ω= 0.00001 when

Q= 1. For Lc= 1000 the respective numbers are 1.926 and 1.928. Thus, the less important

the outer cut-out becomes, the larger are the largest mathematical eigenvalues. Obviously,

the infinite mass in the remote parts of the disk destabilises the inner parts with respect to

m= 1 perturbations, and the behaviour of the doubly cut-out disk should be more applicable

to real galaxies that do not have an infinite mass.

For an Lc= 10 doubly cut-out disk, one may attempt to determine a stability limit for the

m= 1 perturbation. Evaluating the m= 1 kernel for ω= 10−8on an N= 200 grid one

finds that for σ= 0.57 the second-largest mathematical eigenvalue is unity within 0.001, and

raising σfurther, it drops below unity. The value of this result might be somewhat debatable

since σ= 0.57 is really a high velocity dispersion, at which not only my assumption that most

stars are on almost circular orbits breaks down spectacularly but also the idea of a razor-thin

disk becomes completely meaningless unless the approximate equipartition between radial

and vertical degrees of freedom observed in spiral galaxies would be severely violated. Still,

one can give a limit at which hypothetical disks would be stable against m= 1 perturbations.

The susceptibility of the cut-out disks to m= 1 modes might have been expected from a

modal picture, since m= 1 modes do not possess an inner Lindblad resonance that could

absorb density waves, which again means that as soon as some feedback loop is present, the

instabilities will grow much easier than when an ILR is present that has to be hidden by

some kind of Q-barrier. On a related note, Sellwood (1985) found in a series of simulations

of various models of the Galaxy that without a halo, strong m= 1 perturbations develop

on relatively short time scales. An outer cut-out effectively simulates such a halo and thus

works to stabilise the disk.

I postpone the discussion of the solutions of the integral equation to the next chapter, in

which I treat the aesthetically more satisfying two-armed spirals. Suffice it to say that the

solutions do not look qualitatively different from the ones Read found.

The Classic Modes 59

5.3. Two-Armed Modes

The two-armed modes are of particular interest to investigations of spiral structure, mainly

because it were the beautiful patterns of galaxies like M51 who inspired the first researchers

to tackle the problem of the origin of this structure, and consequently, most of the early works

on the subject were heavily biased towards two-armed patterns. Indeed, in the case of M51

and a few more of the classic galaxies, the m= 2 contribution dominates the decomposition

of an image of the galaxy into angular harmonics in both the V- (supposed to trace gas and

young stars) and K-band (which is thought to trace the bulk of the mass). In the case of

M51, this was demonstrated quite convincingly by Rix and Rieke (1993).

However, the galaxies with the most impressive bisymmetric spiral structure almost invariably

are part of a galaxy pair or group, and already the work of Toomre and Toomre (1972)

had clearly shown that the interaction itself can form magnificent two-armed patterns even

without self-gravity and thus instabilities in the disks themselves. Of course, the Toomre’s

own disclaimer cited in the introduction about the spiral structure that would have been

there anyway cautions one not to overinterpret these findings.

Taking a less biased sample than the galaxies humans perceive as the most beautiful ones, the

situation is far more complex (e.g., Naim et. al. 1995) with all kinds of azimuthal harmonics

present in galaxies, including strong m= 1 components. So, maybe Zang’s (1976) discomfort

at finding that Mestel disks are quite stable against m= 2 perturbations was not so justified

after all, even when Mestel disks are supposed to be more than a mere thought experiment.

-2.5

-2

-1.5

-1

-0.5

0.5

1.5

0 5 10 15 20 25 30 35 40 45 50

K(α′)

α′

σ= 0.1

σ= 0.155

σ= 0.386

Fig. 18: The kernel for neutral m= 2 perturbations in the self-consistent disk.

5.3.1. Neutral Mo des in the Self-Consistent Disk

Comparing Figure 18 with Figure 14, one first notices that the neutral m= 2 modes show a

behaviour completely different from the neutral m= 1 modes. Their kernels start out below

zero, decline a bit and then rise. After going through a maximum, they then settle down to

zero for large α′. The rise of the Q= 1 curve (σ= 0.386) is not nearly large enough to let it

cross unity. Thus axisymmetrically stable disks do not admit neutral m= 2 modes.

60 The Classic Modes

The neutral m= 2 modes in the self-consistent disk set in at σ= 0.155 with α′= 25.18. Read

found that for her β=−0.1 disks, the modes set in at σ= 0.156 and for her β= 0.1 disks,

at σ= 0.149 (Read 1999); unfortunately, she had no data available for the Mestel disk itself.

Still, she expects that an interpolation would lead to a useful result (in her own words: “In

general, nothing special happened at β= 0. We wrote code specially for β= 0, using Zang’s

formulae, and the results always agreed with those from the general code with β=±0.001 or

something similarly small”). It thus seems that here, again, my approximations deliver an

excellent agreement with the more exact semi-analytic approach.

Beyond the limit of σ= 0.155, there are two branches of neutral modes since the kernel

crosses the K(α′) line twice. Of course the remarks about the anti-spirality of the neutral

modes apply here as well, so in total there are four branches of neutral m= 2 modes in the

self-consistent Mestel disk.

5.3.2. Mo des in the Cut-Out Disk

The kernel of the m= 2 modes in the cut-out disk (Fig. 21) has all the features familiar by

now, set apart from the m= 1 kernel by its closer adherence to symmetry with respect to the

secondary diagonal. Also, the diagonal appears less important in relation to the off-diagonal

spur. Turning to the eigenvalue spectra in Fig. 19a now, one notices that in contrast to the

m= 1 kernels, there usually are two branches of eigenvalues near the origin. Going towards

eigenvalues of larger moduli, the branches merge to form a single arc of pairs of eigenvalues.

In the superposition of the spectra shown in Fig. 19a, one sees that with decreasing ω, the

spectrum of a kernel is basically rotated clockwise and somewhat elongated. This leads to

a series of arcs of eigenvalues. The most relevant of these arcs is the one with the largest

moduli.

Those lines of eigenvalues of largest modulus are shown for three different values of the growth

rate Im(ω) in Fig. 19c. This figure can be be compared with Read’s results in Fig. 19b, keeping

in mind that these results are for a inner cut-out disk. While generally the approximated

doubly cut-out disks are a little more more stable and the descent towards the real axis is a

little steeper, the results match quite well (the match could be improved increasing Lc, thus

lessening the impact of the outer cut-out).

Furthermore, with the outer cut-out an overtaking of eigenvalues occurs, i.e., the largest

mathematical eigenvalue jumps from one family to a different one. This can be seen in the

plot of the largest eigenvalues as a kink and is more clearly visible in the spectra of Fig. 19a.

For high ω, the largest eigenvalue belongs to the upper arc, whereas for lower ωit belongs

the rightmost arc of pairs of eigenvalues.

The one most interesting number in the stability analysis is the velocity dispersion at which

marginal modes set in; in the M=N= 2 cut-out disk, this happens at σ= 0.214 (Q= 0.554)

with ω= 1.096, corresponding to a pattern speed of Ωp= 0.548. Read’s values (this time for

the doubly cut-out disk) are σ= 0.205 (Q= 0.542) and ω= 1.164, which translates into a

discrepancy of 4% in both σand ω; normalised onto the limit for axisymmetric stability, the

difference in Qis 2%. Given that the velocity dispersion in those disks is still quite high, this

is a surprisingly good agreement.

The qualitative behaviour of the lines of largest eigenvalues as a function of σand ωfor

two-armed modes can be seen in Fig. 19d that can be compared to Read’s Fig. 7.22 (where

The Classic Modes 61

-0.6

-0.4

-0.2

0.2

0.4

0.6

-0.2 0 0.2 0.4 0.6 0.8 1 1.2

Re(λ)

Im(λ)

(a) (b)

-0.15

-0.1

-0.05

0.05

0.1

0.15

0.2

0.25

0.3

-0.1 0 0.1 0.2 0.3 0.4 0.5 0.6

Re(λ)

Im(λ)

Im(ω) = 0.0

Im(ω) = 0.2

Im(ω) = 0.5

()

-0.4

-0.3

-0.2

-0.1

0.1

0.2

0.3

0.4

0 0.2 0.4 0.6 0.8 1

Re(λ)

Im(λ)

σ= 0.232

σ= 0.309

σ= 0.386

(d)

Fig. 19: Eigenvalues of the cut-out m= 2 kernel. (a) shows the spectra for kernels from

ω= 0.1···1.2in 20 steps, the upmost points belonging to the ω= 1.2kernel; (b) is a plot

of the largest mathematical eigenvalues of a Q= 1 disk versus Re(ω)and growth rate sfor

Read’s exact calculations, taken from Read (1997) with kind permission of the author; (c) is

a plot of the largest mathematical eigenvalue versus the indicated values of the growth rate

Im(ω)and Re(ω) = 0.8···2.2in 5 steps; (d) is the same plot varying the velocity dispersion

σinstead.

62 The Classic Modes

-30

-20

-10

-30

-20

-10

α’

-1.4

-1.2

-1

-0.8

-0.6

-0.4

-0.2

0.2

0.4

ΚFig. 21: The real part of the

kernel for m= 2 perturbations

of the cut-out disk. To facili-

tate comparisons, as in the cor-

responding plot for m= 1 this

kernel is for a velocity dispersion

of σ= 0.3.

the lines I computed are for Q= 0.6, Q= 0.8, and Q= 1), still bearing in mind that Read’s

disk has no outer cut-out. Again the approximation reproduces the previous results quite

satisfactorily, with the doubly cut-out disk being a little more stable than its inner cut-out

counterpart.

-0.2

-0.15

-0.1

-0.05

0.05

0.1

0.15

0.2

0.25

0.3

-20 -10 0 10 20 30 40



(



)



-0.3

-0.2

-0.1

0.1

0.2

0.3

-20 -10 0 10 20 30 40



(



)



Fig. 22: The eigenfunctions to the two largest eigenvalues of the kernel for m= 2 perturba-

tions in the doubly cut-out disk. The solid line is the real part, the dashed line the imaginary

part, and the dotted line is the absolute value.

The marginal m= 2 mode is also a good example to discuss the general appearance of the

modes. As in the m= 1 case, the eigenvalues come in pairs, with the difference that in the

two-armed case this only is so after the two branches of the spectrum unite (cf. Fig. 19a). The

obvious question is whether these pairs really are separate eigenvalues or rather a single eigen-

value with a two-dimensional eigenspace that are split up by numerical inaccuracies. While

a thorough treatment of this issue would require some hard-core numerical mathematics, a

“poor-man’s approach” might yield at least an idea of what is going on.

An eigenvector vito the eigenvalue λiof a matrix Mcertainly has δi:= ||Mvi−λivi||/||vi|| =

0. In the present case there is some numerical noise, but with n= 300 and λ1= 1.00207 −

0.0012 and λ2= 0.9889+0.0024, one has δ1= 2.3×10−9and δ2= 4.1×10−9for the marginal

mode determined above, close enough to zero (I have used the Euclidean norm here).

The Classic Modes 63

-0.15

-0.1

-0.05

0.05

0.1

0.15

-0.15 -0.1 -0.05 0 0.05 0.1 0.15

-0.08

-0.06

-0.04

-0.02

0.02

0.04

0.06

0.08

-0.08 -0.06 -0.04 -0.02 0 0.02 0.04 0.06 0.08

Fig. 23: Parametric plots of

δ1(ϕ)(solid line) and ¯

δ2(ϕ)(bro-

ken line). On the left, n= 50, on

the right, n= 300. See text.

If λ1=λ2(i.e., λ1is an eigenvalue f algebraic multiplicity ≥2), clearly the eigenvectors of

λ1would be the eigenvectors of λ2too, and indeed all linear combinations of these would be

eigenvectors as well. If numerical errors have separated the eigenvalues, one would expect

to see a remnant of this behaviour. Since ideally δi(v) := ||Mv −λjv||/||v|| would vanish

when vis a linear combination of the eigenvectors of λj, one should see that ¯

δ1,2(ϕ) :=

δ1,2sin(ϕ)v1+ cos(ϕ)v2is “small”, and it should furthermore not depend too heavily on

ϕ, except where the linear combination suppresses one of the eigenvectors entirely. Given

that the eigenvalues are different, the “small” in the last sentence will not mean 10−9, but it

should be significantly smaller than the δfor an unrelated eigenvalue; for the eigenvector the

the third largest eigenvalue one has, for example, δ1(v3) = 0.15.

Referring to Fig. 23, one sees that the eigenvalues quite probably are distinct single eigen-

values, since with increasing numerical resolution the parametric plot of ¯

δ(ϕ) becomes in-

creasingly elongated, indicating that the mismatch of the sum of the eigenvectors is quite

insensitive to the size of the grid. Note that this result most likely is not due to the fact

that the two eigenvectors are “almost” linear dependent, v1≈ −v2, which, when too severe,

would make the numerics tricky – 

v1/kv1k+v2/kv2k

still is about 0.14, thus leaving enough

headroom for numerics.

On the other hand, this issue is not too important for an assessment of the physics. The

eigenvectors belonging to each of a pairs of eigenvalues will have about the same growth rate

and will mix to some extent regardless if they really belong to a common eigenspace. And

from Fig. 22 and Fig. 24 it is clear that these two modes are very similar anyway.

The basic features of these modes are quite as expected: There is no perturbation within the

inner Lindblad resonance (ILR), the perturbation is quite strong between the ILR and the

corotation (CR), and outside the CR the wave fades away, only barely reaching the outer

Lindblad resonance. It is more or less coincidental that the position of the ILR quite nicely

matches the inner edge of the spiral pattern, although it is tempting to specualte that the

marginal mode occurs at that point in the (σ, ω) space right because the Q-barrier starts to

shield the ILR. In general, though, one can move the Q-barrier out or the ILR in, so that

there would be no waves that can reach the ILR in the first place.

Of course, this discussion is somewhat futile since a M=N= 2 cut-out disk will be badly

instable when the m= 2 set in, and thus one would see everything but structures like those

in Fig. 24.

64 The Classic Modes

Fig. 24: The modes belonging to the two eigenvalues close to one of the m= 2 kernel of

the cut-out disk marginally stable against m= 2 perturbations. The circles mark the outer

Lindblad resonance, the corotation, and the inner Lindblad resonance.

5.4. Three-Armed Modes

Spiral galaxies with a pronounced three-armed pattern are not very common, and for some

cases in which a galaxy shows one, Block and Wainscoat (1991) demonstrated that the stellar

backbone of the galaxy still has a more or less bifold symmetry. Inspired by modal theory,

Block et al (1994) suggested that the fact that the three-armed patterns are not found in

the stellar backbone might indicate an ISM that is dynamically decoupled from the stellar

backbone at least to some extent. On the other hand, it is as well conceivable that nonlinear

coupling between modes of different mcould take place and destabilise an m= 3 mode that

is stable by itself. Thus, although I will show that three-armed modes are yet more stable

than m= 2 modes, there might be some merit in their examination in anticipation of a more

complete theory taking into account nonlinear effects.

5.4.1. Neutral Mo des in the Self-Consistent Disk

As shown in Fig. 25, the kernel for the neutral m= 3 modes starts out below zero for all σ

and finally levels out to zero, never approaching unity in the α-range shown here. However,

the curves for σ= 0.05 and σ= 0.1 are still rising at α= 70, and it seems possible that they

might eventually reach unity at some very large α. However, even for a velocity dispersion as

low as σ= 0.05, this is not the case; the σ= 0.05 curve, for example, reaches a wide maximum

of 0.51 at α≈380. An even cooler disk can be expected to admit neutral modes, but this is

difficult to check because at large α, the Gamma functions in the kernel (120) become really

large, and although they could safely be replaced by their asymptotic expansions, it did not

seem worth the effort to further pursue this issue. A σ= 0.025 kernel is still growing from

0.6 at α≈480, and it seems likely that it will eventually cross the K= 1 line.

The Classic Modes 65

-16

-14

-12

-10

-8

-6

-4

-2

0 10 20 30 40 50 60 70

K(α′)

α′

σ= 0.05

σ= 0.1

Q= 1

Fig. 25: The kernel for neutral m= 3 perturbations in the self-consistent disk.

This is exactly the behaviour Read found for her β= 0 disks.

5.4.2. Mo des in the Cut-Out Disk

The m= 3 modes can be investigated quite analogously to the other non-axisymmetric

modes. The kernel in Fig. 28 continues the pattern of off-diagonal contributions increasingly

dominating over the diagonal as mincreases. Also, the width of the off-diagonal humps has

further increased. Fig. 26a shows the spectra of some kernels. Again, increasing ωbasically

rotates a pattern counterclockwisely around the origin. The pattern itself is now a twofold

arc with something like a tongue in between. It is striking that while the m= 1 spectra had

one branch, the m= 2 spectra had two of them and the m= 3 spectra have three.

-30

-20

-10

-30

-20

-10

α’

-0.8

-0.6

-0.4

-0.2

0.2

0.4

0.6

0.8

Fig. 28: The real part of the

kernel for m= 3 perturbations

of the cut-out disk. As always for

the kernel plots, σ= 0.3.

The curves of largest eigenvalues versus σshows that cut out disks are more stable against

m= 3 perturbations than they are against m= 2 perturbations. The σ= 0.2 disks near

the margin of stability of the two-armed case are safely cross the real axis safely inside unity.

66 The Classic Modes

-0.2

-0.15

-0.1

-0.05

0.05

0.1

0.15

0.2

0.25

0.3

-0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7

Re(λ)

Im(λ)

(a)

-0.4

-0.2

0.2

0.4

0 0.5 1 1.5 2 2.5

Re(λ)

Im(λ)

σ= 0.100

σ= 0.200

σ= 0.386

(b)

()

-0.3

-0.2

-0.1

0.1

0.2

-0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7

Re(λ)

Im(λ)

ℑ(ω) = 0.0

ℑ(ω) = 0.1

ℑ(ω) = 0.2

ℑ(ω) = 0.3

ℑ(ω) = 0.4

(a)

-0.4

-0.2

0.2

0.4

0 0.2 0.4 0.6 0.8 1

Re(λ)

Im(λ)

Im(ω) = 0.0

Im(ω) = 0.1

Im(ω) = 0.2

Im(ω) = 0.3

Im(ω) = 0.4

(b)

Fig. 26: Eigenvalues of the cut-out m= 3 kernel. (a) shows the spectra for σ= 0.205 kernels

from ω= 0.1···2.0in 20 steps, the upmost points belonging to the ω= 2 kernel; (b) is plot

of the largest mathematical eigenvalue versus the indicated values of the velocity dispersion σ

and Re(ω) = 0.8···3.2in 5 steps; (c) are two plots of the largest mathematical eigenvalues

of a Q= 1 (left) and a Q= 0.5disk versus Re(ω) = 1.5···3.9and growth rate s= 0 ···0.4

for Read’s exact calculations, taken from Read (1997) with kind permission of the author; (d)

are two plots of the largest mathematical eigenvalue versus the indicated values of the growth

rate Im(ω)and Re(ω) = 0.5···3.9in 20 steps for a Q= 1 (left) and Q= 0.5disk.

The Classic Modes 67

The plots of largest mathematical eigenvalues versus ωat various growth rates show the

familiar picture except that at very low ωthe curve for the marginal mode bends towards

large Im(λ). This is an artifact of too low numerical resolution, the corresponding kernels

should be evaluated with a larger L.

Lowering the velocity dispersion further finally turns three-armed modes on at about σ=

0.141 with ω= 1.829, corresponding to a pattern speed of Ωp= 0.61.

5.5. Four-Armed Modes

Generally speaking, the four-armed modes continue the pattern of the non-axisymmetric

modes investigated so far in that Mestel disks are yet more stable against them than against

their m= 3 cousins. However, the neutral modes do not fit into the picture drawn so far.

-0.5

0.5

1.5

0 5 10 15 20 25 30 35 40 45 50

K(α′)

α′

σ= 0.2

σ= 0.3

Q= 1

σ= 0.5

Fig. 29: The kernel for neutral m= 4 perturbations in the self-consistent disk.

5.5.1. Neutral Mo des in the Self-Consistent Disk

In Fig. 29 one sees that the m= 4 neutral modes show a behaviour completely different from

all previous cases. The cooler the disk, the lower the value of K(0), but also the greater the

ascent afterwards. Even the σ= 0.5 disk can maintain a neutral mode. This reproduces what

Read found out; in her more general perspective, the emergence of neutral modes at high

temperatures seems to be linked to the existence or non-existence of closed orbits of m-fold

symmetry in the case of three- and four-armed spirals, which, depending on the disk index

β, may occur at m= 3 as well. Her hypothesis that stars on closed orbits “act as absorbers”

and damp away the high-temperature modes, however, cannot be literally true. For one,

neither for m= 1 nor for m= 2 do closed orbits of matching azimuthal symmetry exist, so

high-temperature modes should exist, which, at least for m= 2, they do not. And secondly,

in my disks no closed orbits of any azimuthal symmetry exist, since after the linearisation of

Ω and κthe two frequencies have the fixed ratio √2 and are thus strictly incommensurable.

Still, it seems likely that some kind of resonance does play a role here. Turning to (120), on

sees that the contribution from the terms except the one with the sum is fairly independent

68 The Classic Modes

of m, whereas the sum does significantly change its value with m. The term to watch here is

the m′2/(m′2−2k2). When m′≈√2k, the corresponding summand is much larger than its

neighbours since the ratio Ik(x) and Ik+1(x) basically is k(plus some constant) and everything

else does not depend on k; when mis not too large, that “resonant” summand will dominate

the value of the entire series. The sign of this resonant summand is determined by the side

relative to the nearest kon which m′−√2kreaches its minimum—for the m= 3 series, the

near-resonance is at k= 2 and thus one has a positive contribution with m′−√2k= 0.17

that enters into the kernel negatively (enhancing stability), whereas the resonance for the

m= 4 occurs for m′−√2k=−0.24 and thus destabilises the disk. The next azimuthal

harmonic that has an mclose to a √2kis 7, and it is again negative; not surprisingly, the

m= 7 harmonic is quite unstable as well.

Unfortunately, khas no immediate meaning in physical terms, although it does have a loose

connection to Read’s radial harmonic number l. This latter quantity arises via a discrete

Fourier transformation, whereas the source of the kis (85) that can be viewed as a sort

of modified discrete Fourier transform of some oscillatory function of the auxiliary angle w′

However, it is hard to attach a tangible meaning to w′

1since it only enters the the distribution

function through an integral, and the presence of the Bessel function in the sum does not

exactly clarify the situation. But although kis quite far from being anything like the index

to a Fourier series of as lis in Read’s calculation, it seems that it should eventually end

up imprinting itself as some kind of symmetry in the epicyclic variable w1of the perturbed

distribution function (cf. Eq. 44).

If this is true, the interpretation of the above findings is quite obvious. Let us assume that m

is the azimuthal symmetry of the forcing and ksomehow dictates, say, the number of nodes

in the k-th member of the series describing the response as a function of the epicyclic angle

w′

1. The factor of √2 is simply the epicyclic frequency, so that the two frequencies in question

are the frequency with which a body in the disk encounters a forcing, mΩ, and the frequency

with it would like to respond, kκ. If the forcing is near resonant with with the response times

the epicyclic frequency but has a slightly higher angular symmetry, the response is damped,

whereas in the reverse case it is amplified. This coincides nicely with the well-known fact

that a harmonic oscillator responds out of phase when the forcing is slightly faster than its

natural frequency (the epicyclic frequency times the angular symmetry), whereas it responds

in phase in the reverse case. This situation is quite akin to a Lindblad resonance, with the

difference that for one, the resonance is not exact (if it were, the series would evaluate to

infinity), and, for second, it occurs on the entire disk.

All this does not really answer the question whether the bewildering stability behaviour is an

artifact of the linearisation or not, although it does suggest that in general m= 4 modes in

the self-consistent disk should be more unstable than their m= 3 counterparts. For rotating

modes, on the other hand, this quasi-resonance should be broken for most part of the disk,

since m(Ω+Ωp) and kκ no longer have a constant ratio. Of course, for a perfectly self-similar

disk, the phrase “most part of the disk” is quite meaningless, and thus there is little to keep

rotating modes from showing a similar behaviour. Again, the unusual self-similarity of the

Mestel disk makes for some quite bizarre features.

The Classic Modes 69

-0.6

-0.4

-0.2

0.2

-0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7

Re(λ)

Im(λ)

(a)

-0.2

-0.15

-0.1

-0.05

0.05

0.1

0.15

0.2

-0.1 0 0.1 0.2 0.3 0.4 0.5 0.6

Re(λ)

Im(λ)

σ= 0.00

σ= 0.25

σ= 0.50

(b)

()

-0.2

-0.15

-0.1

-0.05

0.05

0.1

0.15

0.2

0.25

-0.1 0 0.1 0.2 0.3 0.4 0.5 0.6

Re(λ)

Im(λ)

Im(ω) = 0.0

Im(ω) = 0.1

Im(ω) = 0.2

Im(ω) = 0.3

Im(ω) = 0.4

(a)

-0.4

-0.2

0.2

0.4

0 0.2 0.4 0.6 0.8 1

Re(λ)

Im(λ)

Im(ω) = 0.0

Im(ω) = 0.1

Im(ω) = 0.2

Im(ω) = 0.3

Im(ω) = 0.4

(b)

Fig. 30: Eigenvalues of the cut-out m= 4 kernel. (a) shows the spectra for σ= 0.205 kernels

from ω= 0.1···2.0in 20 steps, the upmost points belonging to the ω= 2 kernel; (b) is plot

of the largest mathematical eigenvalue versus the indicated values of the velocity dispersion σ

and Re(ω) = 0.8···3.2in 5 steps; (c) are two plots of the largest mathematical eigenvalues

of a Q= 1 (left) and a Q= 0.5disk versus Re(ω) = 1.2···5.2and growth rate s= 0 ···0.4

for Read’s exact calculations, taken from Read (1997) with kind permission of the author; (d)

are two plots of the largest mathematical eigenvalue versus the indicated values of the growth

rate Im(ω)and Re(ω) = 1 ···5.2in 20 steps for a Q= 1 (left) and Q= 0.5disk.

70 The Classic Modes

-30

-20

-10

-30

-20

-10

α’

-0.4

-0.2

0.2

0.4

Fig. 32: The real part of the

kernel for m= 4 perturbations

of the cut-out disk. As always,

σ= 0.3.

5.5.2. Mo des in the Cut-Out Disk

In contrast to the neutral modes in the singular disk, almost nothing special happens in the

cut-out disk. The kernel in Fig. 32 has yet wider and more dominant humps, and the spectra

in Fig. 15a are still somewhat more complex than they were with m= 3, and the disk is yet

a little more stable against m= 4 perturbations than it is against three-armed ones.

The quite complex spectra also show in the plots of the largest eigenvalues, particularly in the

Q= 0.5 case where quite a few kinks in both Read’s and my diagram witness that eigenvalues

corresponding to different eigenvectors overtake one another rather frequently.

With this kind of stability behaviour, it is not surprising that m= 4 set in at even lower

temperatures; the limit below which modes set in here is σ= 0.124 at ω= 2.35, corresponding

to a pattern speed of Ωp= 0.588.

The qualification “almost” with the above statement that nothing special happens in the

cut-out disk was necessary nevertheless. Just as in Read’s formalism, for very low ωthe

eigenvalues start to return through the third quadrant to the real axis, and given that for

N= 200, σ= 0.2 and ω= 0.01 the largest mathematical eigenvalue is λ= 1 = 1.354 −0.001,

one has to expect that modes will show up in this region. Clearly these modes feel the same

resonance that caused the non-quenchable neutral modes, which had to be expected after the

above discussion. Again it is hard to tell if these modes have much physical significance.

More Modes? 71

6. More Mo des?

One major point that remained somewhat open in the works of both Zang and Read was

the issue whether or not the self-consistent disk does indeed admit other modes than the

rather exotic neutral modes investigated above. The behaviour in some limits in which the

cut-out disk should approach the self-consistent disk (most notably, ω→0) suggests that no

such modes exist, but a conclusive answer could not be given due to the singular nature of

the kernel. The main discomfort with the notion of modes in the singular disk is that in its

governing equation ωfactors out, so that without further devices a disk is either stable or

becomes unstable at all growth rates and pattern speeds at once. This seemed weird enough

to dismiss the possibility of such modes. On the other hand, this would imply that even

a completely cold disk would be stable to all nonaxisymmetric perturbations, which again

would be quite an exotic situation.

Quite recently, Goodman and Evans (1999) have tackled the problem from a different per-

spective, employing the Jeans equations to revisit the problem of modes in the self-consistent

Mestel disk. The basic results of their work are that the Mestel disk’s weirdness starts in its

centre, and once one fixes boundary conditions there, the Mestel disk becomes quite tame.

By the requirement that the centre does not absorb or emit energy, they curb down the

two-dimensional continuum of solutions to a spectrum of one-dimensional continua, in which

the ratio between growth rate and pattern speed is fixed and only the phase of ωremains

unknown. Even this ambiguity can be removed by fixing the phase shift during the reflection

of a wave at the disk’s centre, thus reducing the spectrum to a countable set of discrete

points; however, there currently are no good indications as to how this phase shift should

fixed. Of course, this only raises new questions. A self-similar disk cannot just admit a single

mode without also admitting its scaled brethren; so, somewhere Goodman and Evans must

have broken the self-similarity simply by requiring a centre that conserves energy. How this

happens is completely unknow at this point.

With a closed expression of the kernel of the Mestel disk, one is in a position to reconsider

the stellar-dynamical equivalent of the Goodman and Evans’ problem using the tools of the

theory of singular integral equations: Does the integral equation (111) admit continuous

solutions, and if, for what disk parameters?

6.1. The Singular Disk revisited

The kernel (118) of the equation governing the self-consistent disk has a singularity on the

diagonal, and thus the familiar Fredholm theory does not apply to it since the operator

RK(α, α′)Φ(α)dα is not compact. The theory that replaces Fredholm theory in this case is

developed in Muskhelishvili (1953) and Gakhov (1990), and is nicely summarised in the more

modern work of Polyanin and Manzhirov (1998).

72 More Modes?

The integral equation (111) has the form

0 = a(α′)Φ(α′) + 1

iπZL

K(α, α′)

α−α′Φ(α)dα, (146)

where ˆ

Kis at least H¨older continuous, Lis the real axis (integration over the real axis is

henceforth implied for integral signs not otherwise marked), I wrote Φ(α′) instead of Φα′,m

to accentuate that Φ is a (complex) function of the real variable α′, and, of course, the

integral is to be understood as a Cauchy principal value. This equation is a homogenous

equation of the second kind (or, in Pipkin’s (1991) terminology, of the third kind). In my

case, a(α′) = F(α′)−1 and ˆ

K= (α−α′)K(α, α′)/iπ. The multiplication with (α−α′) ensures

that the kernel is bounded on the diagonal, since the single zero in the denominator stemming

from exp2π(α−α′)−1 now cancels out. The only part of the kernel that could jeopardise

the condition of H¨older continuity is the square root, but since α′only enters squared in the

arguments of the square roots, ˆ

Kdoes indeed satisfy a H¨older condition.

It is convenient to rewrite (146) to the equivalent form

0 = a(α′)Φ(α′) + 1

iπb(α′)ZΦ(α)

α−α′dα +1

iπZˆ

Kr(α, α′)Φ(α)dα, (147)

where

b(α) = ˆ

K(α′, α′) and

Kr(α, α′) = ˆ

K(α, α′)−ˆ

K(α′, α′)

α−α′.

(148)

It is quite easy to see that for the integral equation (111), b(α′) = F(α′). After this manipu-

lation, ˆ

Kris a compact operator, whereas

Φ := a(α′)Φ(α′) + 1

iπb(α′)ZΦ(α)

α−α′dα = 0,(149)

known as the dominant or characteristic part of the integral equation, decides if the integral

equation has a nontrivial solution.

To find a solution of (149), let us introduce an auxiliary function

Ψ(z) = 1

2πiZΦ(α)

(α−z)dα. (150)

If Φ(α) is continuous, this function is piecewise analytic in the upper and lower half planes

but may be discontinuous across the real axis. The merit of this auxiliary function is that it

allows to transform the problem of the solution of the integral equation (149) into a Riemann

boundary value problem.

This can be done by applying the Sokhotski-Plemelj formulae, that, in their formulation for

the real axis, state that when Φ and Ψ are linked by (150),

Ψ+(α′) + Ψ−(α′) = 1

iπZΦ(α)

α−α′dα

Ψ+(α′)−Ψ−(α′) = Φ(α′)

α′∈R(151)

More Modes? 73

hold. Here, Ψ+(α′) and Ψ−(α′) are the limiting values of Ψ(z) as z→α′from above and

below, respectively.

Inserting (151) into (149) and collecting the terms yields

Ψ+(α′) = a(α′)−b(α′)

a(α′) + b(α′)Ψ−(α′) =: D(α′)Ψ−(α′).(152)

Going back to the integral equation for the singular disk, one has

D(α′) = −1/(2F(α′)−1).

This represents a Riemann boundary value problem, asking (in the special case used here) for

two functions Ψ+and Ψ−that are analytic on the upper and lower half planes, respectively,

and that satisfy Ψ+(α′) = D(α′)Ψ−(α′) for all real α′.

For reasons that will become clear below, I for the time being demand that Dhas neither

zeroes nor poles on the real axis. When will this be the case? To answer this question, consider

the plots of Fover αin Fig. 33 for m= 1 ···4 and various values of the velocity dispersion σ

(as shown above, the axisymmetric case does not lead to a singular integral equation). The

condition of a bounded and nonzero Dwill be satisfied if F(α)6= 1/2 (α∈R), where F(α)

clearly tends to zero for large α.

For m= 1, Freaches the 0.5 line even for very large values of σ, although this might happen

at very large αwhen σis small. In this case, Dwould have two poles of order one over

the entire range of σfor which this analysis could even remotely be applicable. This nicely

corresponds to the stability behaviour for neutral modes.

The bisymmetric case, on the other hand is “well-behaved” as long as σ > 0.155, showing

four poles of order one below this threshold, reproducing the stability limit of the neutral

modes. After this, it is probably not too surprising that since this stability limit was very

low indeed in the three-armed case, Dis well-behaved for all but very low σfor m= 3. At

m= 4, the quite pathologic behaviour found for the neutral modes crops up again and leads

to four poles of order one for all but ridiculously high σ.

Calling this suggestive probably is a gross understatement. On physical grounds, one would

thus expect (111) to have no solution in the regular case of a bounded and nonzero D, and at

least one marginal solution (given that the neutral modes have two or four branches, it would

not be too surprising at least the trailing branches would show up here again) otherwise.

With this concept mind, let us return to the analysis of the boundary value problem, still

restricting the analysis to the regular case.

The D(α′) treated here is the limiting value of an analytic function (in particular, according

to my assumptions it has no poles on the real axis), and thus is guaranteed to be H¨older

continuous. Furthermore, D(α′) is a positive real function for real α′. This ensures that its

index, the increment of its argument over the integration contour,

Ind(D) = 1

2πZdarg D(α)dα, (154)

is zero. This in turn implies that ln Dis a well-defined (single-valued) function. Then one

can take the logarithm of both sides of (152), to arrive at

ln Ψ+(α′)−ln Ψ−(α′) = ln D(α′).(155)

74 More Modes?

-0.5

0.5

1.5

2.5

3.5

-50 -40 -30 -20 -10 0 10 20 30 40 50

(



)





= 0



= 0



= 0

-1.2

-1

-0.8

-0.6

-0.4

-0.2

0.2

0.4

0.6

0.8

-50 -40 -30 -20 -10 0 10 20 30 40 50

(



)





= 0



= 0

155



= 0

-8.0

-7.0

-6.0

-5.0

-4.0

-3.0

-2.0

-1.0

0.0

1.0

-200 -150 -100 -50 0 50 100 150 200

(



)





= 0



= 0



= 0

-0.5

0.5

1.5

-50 -40 -30 -20 -10 0 10 20 30 40 50

(



)





= 0



= 0



= 0

Fig. 33: Plots of the function F(α)for m= 1,2,3,4(from top to bottom) and the indicated

values of σ. The solid horizontal line marks the F(α) = 1/2line. When Fcrosses this line,

Dwill have singularities.

More Modes? 75

Using the second of the Sokhotski-Plemelj formulae (151), it follows that with

G(z) = 1

2πiZln D(α)

α−zdα (156)

one can write down a solution to the boundary value problem (152)

X+(α′) = eG+(α′)X−(α′) = eG−(α′),(157)

provided that the problem does have a solution.

Thus one has D=X+/X−, and the boundary value problem (152) can be written as

Ψ+(α′)

X+(α′)=Ψ−(α′)

X−(α′).(158)

Now, the functions G±by construction have no pole in their respective domains, and therefore

X±has no zeroes. If Φ±are defined via an integral as in (150), this implies that they are

analytic in their respective half-planes, with the possible exception of ∞. Thus, the functions

on the two sides of (158) are analytic in their respective half-planes, and they are identical

on the real axis. By Morera’s theorem this equation implies that the right-hand side and the

left hand side represent the same analytic function, S(z), say. Now, what are the properties

of S(z)?

I already pointed out that S(z) is analytic in the entire plane with the possible exception of

the point at infinity, where Ψ±(∞) = ±1

2Φ(∞). Since I require any solution of the original

integral equation (111) to vanish at infinity, one has Ψ±(∞) = 0 as well. Since D→1 for

large α,Gapproaches 0 and X±tends to 1. This implies that Svanishes at infinity. Since

Sis an entire function, Liouville’s theorem then states that Smust be identically zero. By

the first Sokhotski-Plemelj formula one now has Φ(α′) = S(α)(X++X−), and hence the

homogenous integral equation has no nontrivial solution. This is a special case of the general

result that is proven somewhat more rigorously in the literature cited above: A homogenous

Cauchy integral equation with an index of zero is unsolvable.

Loosening the conditions on D, I now allow it to have “poles” at some points αi, which, for

the sake of simplicity, are supposed to be poles of order one. The term “pole” has been put

in quotes since it is not necessary that Dis analytic in C\{αi}, only that the D(αi) behaves

like (α′−αi)−1at the exceptional points.

Writing

D(α′) = D′(α′)

i=1(α′−αi),(159)

the boundary condition (152) takes the form

Ψ+(α′) = D′(α′)

i=1(α′−αi)Ψ−(α′),(160)

where D′now satisfies the conditions I required of the Din (152). Solving the boundary

value problem for this D′and substituting the result back into (160), the equivalent of (158)

now is Ψ+(α′)

X+(α′)=Ψ−(α′)

X−(α′)Qn

i=1(α′−αi),(161)

76 More Modes?

where this time X= eGand G(α′) = 1

2πiRln D′(α′)/(α−α′)dα. Let us again call the

analytic function defined in the left and right sides of (161) S(z).

By the right side of (161), Snow may have poles of order 1 at αi. By the generalzed Liouville

theorem, a solution of the integral equation would then have the form

Φ = 1

Qci(α′−αi)X++X−Y(α′−αi).(162)

Evidently, solutions of this type cannot be bounded at αi, since that would require the bracket

to vanish at αi, implying X+(αi) = 0 which is impossible due to the definition of Xvia an

exponential of a function bounded for finite arguments. Thus, no solutions exist even in the

exceptional case.

Certainly, this is not what was to be expected after the above analysis. The reason for this

surprise is of course the requirement of boundedness for Φ and thus for Ψ±. Unfortunately,

the requirement that Φ does not possess poles is necessary for (150) to be even defined. If

Φ(α) has as much as one singularity of order one, the integrad in (150) will have a singularity

of order two at that point and will thus no longer be integrable. Indeed, the integral equation

(147) itself cannot admit any Φ(α) with a pole, since it would render the integral in the second

summand undefined. It is still conceivable that Φ(α) might have an integrable singularity of

some kind, but then of course the basic assumption that any singularity Smight have is a

pole breaks down and the entire theory does not work any more. Hence, for want of a theory

for nonbounded solutions of Cauchy integral equations, once again the existence of solutions

cannot be conclusively proven, although I can rule out bounded solutions.

Reviewing the solutions put forward by Goodman and Evans (1999), one finds that they do

show singularities right at the exceptional points of (161), which is where one would expect

singularities in the present analysis. This might be sufficient encouragement to actually try

and assume solutions that may be found numerically. To do this, one has to solve the full

equation (147). If both Fand the solution were bounded, this would be quite an easy

task, since dealing with the singularity in the kernel appears to be feasible. In preparatory

experiments, the midpoint quadrature formula (constructed along the ideas laid down in

Press et al, 1992)

−L

f(α)

(α−α′)dα =

i=0

j=0

χj(αi+j, α′)Φ(αi+j),(163)

where αi=L((1+2i)

4n−1) and χi(α, α′) is the (straightforward but lengthy) solution of

Za+7/2h

a−h/2

f(α)

α−α′dα =f(a)χ0(a, α′)+f(a+h)χ1(a, α′)+f(a+2h)χ2(a, α′)+f(a+3h)χ3(a, α′),

(164)

f(x) = xi, i = 0 ···3, with h=L/(2n) worked quite well for the kernel and steers clear of

evaluating the kernel on the diagonal. However, the only cases that can have solutions are

those with a nonbounded Fand a solution that will be in some way “bad” at the points

where Fdiverges. Handling this appears to be at least difficult.

Once one had a solution, after the findings of Goodman and Evans one would certainly want

to weed out solutions not satisfying some appropriate boundary conditions. This again does

More Modes? 77

not seem trivial in that it is hard to separate incoming and outgoing components in solutions

like Fig. 22. All in all, the scope of such a project seems more like another thesis than

a footnote to this one, and consequently I have little choice but to leave this issue open,

although I strongly suspect that the neutral modes do indeed have families of rotating and

growing modes behind them. Given that the stability of the neutral modes is governed by the

fact that for them an approximate resonance between the frequency of the forcing and some

harmonic of the epicyclic frequence stretches over the entire disk and for true rotating modes

that is not the case, one can expect this procedure to lead to quite different stability limits for

rotating modes. On the other hand, any pattern speed can be scaled to be arbitrarily small

as long as the disk’s self-similarity is not broken, and thus without “tricks” to introduce some

scale into the disk the stability limits suggested here (if indeed the exceptional case admits

modes) are valid. The most interesting case in an analysis employing such “tricks” should

be the three-armed one. If the above analysis is correct, the inhibition of these modes by

the resonance should pertain to rotating modes as well, and only when the self-similarity is

broken in some way should one expect to see the stability limit rise.

78 Summary and Conclusions

7. Summary and Conlusions

In this work, I have investigated the modes admitted by Mestel disks. This problem has

already been investigated by Zang (1976) and Read (1997), who chose a semianalytic approach

that essentially relied on two approximations:

•The treatment of the disk as razor-thin.

•The linearisation of the Boltzmann equation (i.e., the perturbation amplitude is a small

quantity)

While their results are as exact as one wishes them to be, their method does not deliver

the kernel in a closed form—to compute it, some numerical integration is required. I traded

in the exactness for (almost) closed expressions for the kernel by employing two further

approximations:

•I approximate the orbits of bodies in the Mestel disk by epicyclic orbits, first deriving

the two periods of the orbits by the method of harmonic balance, then linearising those

simplified expressions.

•I (sort of) linearise the logarithmic spirals used as the potential basis

While the second of these approximations should not have much impact since slightly chang-

ing basis functions for an expansion is rather safe, the first one is quite serious, in particular

because I still pretend that the coordinates I derive from these approximations are canonic

action-angle coordinates, which, of course, they are only approximately. In particular the

second linearisation of the expressions for the actions and angles can be expected to have a

quite drastic impact on the results, basically reducing my approach to an epicyclic approxi-

mation of order 1.1, as it were. A strict first-order epicyclic approximation is unusable in the

Mestel disk, since the resulting Hamiltonian does not contain the second integral of motion.

Following Zang, in addition to the self-consistent Mestel disk I also investigated a disk in which

a part of the matter does not take part in the perturbation. While the main purpose of doing

this is to remove some properties of the self-consistent disk that appear quite unphysical and

also to introduce a Q-barrier, it is more than a pure mathematical device in that the cut-outs

have real-world counterparts in dark haloes and bulges that influence the rotation curve of

the disks without being dynamically active enough to take part in spiral perturbations. The

cut-out function treated here is one of the cut-out functions investigated by Read in order

to facilitate comparisons with her results. For further investigations, I suggest a different

(one-parameter) cut-out that results in much more tractable expressions.

Despite the cautioning words above, in both the self-consistent disk and the cut-out disk the

precision of the results is good enough to justify the approximations—the physics come out

right to within 10% of the results of Read even in the more critical cases, and to within 2%

when the problems are tame. Table 1 summarises the stability properties of the self-consistent

and the cut-out disks with respect to the various azimuthal harmonics.

Summary and Conclusions 79

self-consistent cut-out

m σmin α σmin Ωp

0 0.386 3.51 0.349 −1

1−2−20.5730.003

2 0.155 25.18 0.214 0.554

3<0.05 >600 0.141 0.61

4−4−40.124 0.588

1m= 0 modes have no orientation and therefore no

pattern speed.

2See the discussion in chapter 5.2.

3Even slightly rotating modes have more realistic stabil-

ity properties. See chapter 5.2.

4See chapter 5.5.

Table 1: A summary of the stability properties of the Mestel disk. “Cut-out” refers to the

M=N= 2 cut-out disk (see chapter 4.1), self-consistent refers to the plain Mestel disk. The

modes in the self-consistent disk become unstable for all Ωpat the same time, at least as far

as one can tell from my investigation. The parameter αis the logarithmic wave number of

the neutral mode. The modes in the cut-out disks have contributions from a continuum of α

but have a definite pattern speed Ωp.

So, the picture of the Mestel disk is mainly like the one painted by Zang and Read: The

axisymmetric stability is quite like predicted by the local theory. The disk readily admits

growing m= 1 modes, whereas modes of higher azimuthal symmetry become unstable only

when the disk already is axisymmetrically unstable and should thus not occur in disks suf-

ficiently similar to a Mestel disk. The stability against neutral modes in the self-consistent

disk in retrospect appears to be closely linked to the proximity of a harmonic of the epicyclic

frequency κthis the frequency of the forcing mΩ. For m= 2, mΩ lies quite in between

two harmonics of κ, and the disk becomes unstable at σ= 0.155; for m= 3, the driving is

somewhat faster than than the closest kκ, inducing a response out of phase, so that the disk

is remarkably stable; and for m= 4, the driving is somewhat slower, and the disk is very

unstable. This becomes a global property of the disk for the (non-rotating) neutral modes

because the frequency ratio between forcing and epicyclic frequency is constant over the disk,

which for rotating modes it is not. Basically, the neutral modes can exhibit a near Lindblad

resonance over the entire disk. The m= 1 case would fit into the m= 4 picture, although

certain differences suggest that its instability might be linked to, e.g., an instability of the

position of the cusp.

The most valuable asset this work adds to the body of knowledge piled up by Zang and Read

is that it gives a closed form of the kernel. This may lead a way for further investigating

Mestel disks or even more general disks, in particular with respect to their properties in the

regime of mildly nonlinear perturbations. I did not follow those possible ways but restricted

myself to a treatment of the issues raised by Zang and Read, where I shed some light on them

from a slightly different perspective.

The main issue that remains to be solved to facilitate a completely analytic treatment of

the problem is a closer investigation of the series involving Bessel functions in the kernels.

80 Summary and Conclusions

Although I have been unsuccessful in trying to find a way to find at least approximative

analytic solutions, it is quite conceivable that the kernel is close enough to a degenerate one,

so that the eigenvalues and eigenfunctions governing the problem might be writable in a

closed form.

One point in which this work goes beyond Zang and Read is that the closed form of the

kernel enabled me to explore the stability properties of the self-consistent disk with a broader

set of tools. Although I can show that no modes with bounded Mellin transforms exist in

the self-consistent disk, for velocity dispersions below the stability limit of the neutral modes

this does not rule out the existence of rotating modes with some integrable singularities.

Indeed, the structure of the problem strongly suggests that they should exist, even more so

since the modes Goodman and Evans (1999) put forward from their gas-dynamical analysis

show singularities right at the points where one would expect them to be from my analysis.

However, without a solution theory for Cauchy integral equations on the space of functions

with integrable singularities, their existence cannot be proven. Thus, once more the existence

of ω6= 0-modes in the unmodified Mestel disk remains open.

Summary and Conclusions 81

Aknowledgement

A significant part of this work was done while the author was supported by a grant under

the Landesgraduiertenf¨orderungsgesetz des Landes Baden-W¨urttemberg.

I want to express my gratitude to the students and the staff at the Astronomisches Rechen-

Institut Heidelberg, in particular my supervisor Burkhard Fuchs, who not only provided

advice and encouragement but also kindly arranged the defence while the author was not in

Heidelberg, and to Holger Baumgart, Michael Biermann, Sabine Frink, and my office-mate

Heiko Reffert, who provided a friendly environment.

Also, I want to thank Jenny Read for useful discussions, sending me a copy of her thesis, and

the permission to reprint some figures from that work.

Furthermore, I want to acknowledge Guido von Rossum, the author of the wonderful Python

language, Linus Thorvalds, the author of the Linux operating system, Donald E. Knuth, the

author of the T

EX typesetting system, Richard Stallman, spiritus rector of the GNU project,

and countless other contributors to tools like Gnuplot, netpbm, NumPy, meschach and others

that were instrumental to the present work.

Finally, I feel obliged to mention Klaus von Trotha, Peter Ulmer and J¨urgen Siebke. Their

unrelenting quest for the University of the Future provided endless distraction and entertain-

ment.

82 Appendix

8. App endix

In chapter 4.4 I did not write down the results of the final quadrature over Lzfor the cut-out

disk. In this appendix I give the results of their computation by means of computer algebra.

These terms could be simplified to result in expressions about 1/4 the size of the expressions

below, if one may extrapolate the behaviour of the analogous expressions for the simpler

cut-out (122). However, even then the kernel would be extremely unwieldy, and thus I chose

to simply dump the expressions actually used to compute the kernel in the C++ programme

to a T

EX-source.

Maple sometimes abbreviates common subexpressions with %n. The definitions of the subex-

pressions can be found at the end of the formulae.

I(Lz,1)

2=1

4i ( L0( i ( α−α′) ) −Lc( i ( α−α′) ) ) %2 

Γ1

4+1

2m+1

2iα′

Lc21−2σ2(1/2−1/21

σ2)e(1/21

σ2)

e−1/42+σ4α2−2σ4α α′+σ4α′2

σ2(−1 + e(−π(α−α′) ) ) ( L0+Lc) ( L0−Lc)



Γ3

4+1

2m+1

2iα′

eσ2−1

2i1−2σ2(1/2−1/21

σ2)e(1/21

σ2)Lc2



Γ1

4+1

2m+1

2iα′

e−1/42+4 σ4m2+σ4α′2+σ4α2

σ2

BesselI 0,1

2qσ42m2α2+ 2 m2α′2+α2α′2+ 4 m4

−1

2L0( i ( α−α′) ) 

−%2 L04ω4−2 %2 L02ω2m2−%2 L03ω3σ2m(α−α′)

+ i %2 L02ω2σ2m2(α−α′)−%2 L0ω σ2m3(α−α′)−%2 m4

+ i %2 ( α−α′)m4σ2−%2 σ2m2L02ω2+ 2 i %2 σ2m3ω L0+ %2 σ2m4

−%1 L04ω4−2 %1 L02ω2m2+ %1 L03ω3σ2m(α−α′)

+ i %1 L02ω2σ2m2(α−α′) + %1 L0ω σ2m3(α−α′)−%1 m4

+ i %1 ( α−α′)m4σ2−%1 σ2m2L02ω2−2 i %1 σ2m3ω L0+ %1 σ2m4.

(L0+Lc) ( L0−Lc) ( −1 + e(−2π(α−α′) ) )−L02ω2−2 i m ω L0+m2

−L02ω2+ 2 i m ω L0+m2+1

2Lc( i ( α−α′) ) −%2 m4−2 %2 m2Lc2ω2

−%2 Lc4ω4+ %2 σ2m4+ 2 i %2 σ2m3ω Lc−%2 σ2m2Lc2ω2

Appendix 83

−%2 Lcω m3σ2(α−α′) + i %2 Lc2ω2m2σ2(α−α′)−%2 Lc3ω3m σ2(α−α′)

+ i %2 ( α−α′)m4σ2−%1 m4−2 %1 m2Lc2ω2−%1 Lc4ω4+ %1 σ2m4

−2 i %1 σ2m3ω Lc−%1 σ2m2Lc2ω2+ %1 Lcω m3σ2(α−α′)

+ i %1 Lc2ω2m2σ2(α−α′) + %1 Lc3ω3m σ2(α−α′) + i %1 ( α−α′)m4σ2.

(−1 + e(−2π(α−α′) ) ) ( L0−Lc) ( L0+Lc)m2−2 i m ω Lc−Lc2ω2

m2+ 2 i m ω Lc−Lc2ω2

σ2m2ω2−ω4L02Lc2+m2Lc2ω2+ω2m2L02+ 3 m4ω

m(−i ( α−α′) )

Lc2ω2+m22m2+L02ω22(−1 + e( 2 π(α−α′) ) )

e

σ2

Γ3

4+1

2m+1

2iα′

2

%1 := e(−3/2π(α−α′) )

%2 := e(−1/2π(α−α′) )

I(Lz,2)

2=−1

8i1−2σ2(1/2−1/21

σ2)√2 e(1/21

σ2)Lc2

Γ1

4+1

2m+1

2iα′

e−1/42+4 σ4m2+σ4α′2+σ4α2

σ2

BesselI k, 1

2qσ42m2+α′2(2 m2+α2)−√%1 L0( i ( α−α′) ) π

−4σ2m4cos( %2 ) √2 %1 ξ82−iL0σ2m ω √2 cos( %2 ) ξ73

−2m2cos( %2 ) √2 %1 ξ82ξ7 + 8 i L0σ2m ω √2 cos( %2 ) k2ξ72

−2 i σ2m k sin( %2 ) ξ73−2 i σ2m k sin( %2 ) %1 ξ82ξ7

+ 2 i σ2m2cos( %2 ) √2 ( α−α′) %1 ξ82ξ7−4L0ω k sin( %2 ) %1 ξ82ξ7

+ 2 L02ω2√2 cos( %2 ) %1 ξ82ξ7−σ2m2cos( %2 ) √2ξ73

−2 i σ2m k sin( %2 ) ξ83%1 −2m2cos( %2 ) √2ξ8ξ72

+ 4 L0ω k sin( %2 ) ξ8ξ72−8 i σ2m3ksin( %2 ) ξ72

−4σ2m k sin( %2 ) ( α−α′)ξ8ξ72+ 16 i σ2m k3sin( %2 ) %1 ξ82

+ i L0σ2m ω √2 cos( %2 ) ξ83%1 −4σ2m4cos( %2 ) √2ξ72

−iL0σ2m ω √2 cos( %2 ) %1 ξ82ξ7−4 i L0ωcos( %2 ) √2m ξ8ξ72

+ 2 L0σ2m ω √2 cos( %2 ) ( α−α′) %1 ξ82ξ7

−8 i L0σ2m ω √2 cos( %2 ) k2%1 ξ82−σ2m2cos( %2 ) √2 %1 ξ82ξ7

+ 8 σ2m2cos( %2 ) √2k2%1 ξ82−2 i σ2m k sin( %2 ) ξ8ξ72

−4 i m k sin( %2 ) ξ8ξ72−8 i σ2m3ksin( %2 ) %1 ξ82

+ 4 i L0ωcos( %2 ) √2m%1 ξ82ξ7−2L0 σ2m ω √2 cos( %2 ) ( α−α′)ξ8ξ72

−σ2m2cos( %2 ) √2ξ8ξ72+ 2 L02ω2√2 cos( %2 ) ξ8ξ72

+ 8 σ2m2cos( %2 ) √2k2ξ72+ 16 i σ2m k3sin( %2 ) ξ72

84 Appendix

−4σ2m k sin( %2 ) ( α−α′) %1 ξ82ξ7 + 2 i σ2m2cos( %2 ) √2 ( α−α′)ξ8ξ72

+ i L0σ2m ω √2 cos( %2 ) ξ8ξ72+ 4 i L0σ2m3ω√2 cos( %2 ) %1 ξ82

−4 i m k sin( %2 ) %1ξ82ξ7−4 i L0σ2m3ω√2 cos( %2 ) ξ72

−σ2m2cos( %2 ) √2ξ83%1.ξ72ξ82( %1 −1 ) ( %1 + 1 ) ( L0 −Lc) ( L0+Lc)

+√%1 Lc( i ( α−α′) ) πim σ2Lcω√2 cos( %2 ) ξ103%1

+ i m σ2Lcω√2 cos( %2 ) ξ92ξ10 + 4 i ωcos( %2 ) √2mLc %1 ξ9ξ102

−8 i m σ2Lcω√2 cos( %2 ) k2%1 ξ102+ 8 i m σ2Lcω√2 cos( %2 ) k2ξ92

+ 4 i m3σ2Lcω√2 cos( %2 ) %1 ξ102−im σ2Lcω√2 cos( %2 ) %1 ξ9ξ102

+ 2 i σ2m2cos( %2 ) √2 %1 ( α−α′)ξ9ξ102

+ 2 m σ2Lcω√2 cos( %2 ) %1 ( α−α′)ξ9ξ102−im σ2Lcω√2 cos( %2 ) ξ93

−2m2cos( %2 ) √2ξ92ξ10 −4 i ωcos( %2 ) √2m Lcξ92ξ10

−2 i σ2m k sin( %2 ) %1 ξ9ξ102+ 2 i σ2m2cos( %2 ) √2 ( α−α′)ξ92ξ10

−4 i m3σ2Lcω√2 cos( %2 ) ξ92−4σ2m k sin( %2 ) %1 ( α−α′)ξ9ξ102

−σ2m2cos( %2 ) √2ξ93−4σ2m4cos( %2 ) √2ξ92

−4 i m k sin( %2 ) %1ξ9ξ102−2m σ2Lcω√2 cos( %2 ) ( α−α′)ξ92ξ10

+ 8 σ2m2cos( %2 ) √2 %1 k2ξ102+ 2 Lc2ω2√2 cos( %2 ) %1 ξ9ξ102

−2 i σ2m k sin( %2 ) ξ103%1 + 16 i σ2m k3sin( %2 ) %1 ξ102

−4ω k sin( %2 ) Lc %1 ξ9ξ102−2 i σ2m k sin( %2 ) ξ92ξ10

−σ2m2cos( %2 ) √2 %1 ξ9ξ102−8 i σ2m3ksin( %2 ) %1 ξ102

+ 2 Lc2ω2√2 cos( %2 ) ξ92ξ10 −4σ2m k sin( %2 ) ( α−α′)ξ92ξ10

−σ2m2cos( %2 ) √2ξ92ξ10 −8 i σ2m3ksin( %2 ) ξ92

+ 8 σ2m2cos( %2 ) √2k2ξ92−σ2m2cos( %2 ) √2ξ103%1

−4σ2m4cos( %2 ) √2 %1 ξ102+ 16 i σ2m k3sin( %2 ) ξ92

−2m2cos( %2 ) √2 %1 ξ9ξ102+ 4 ω k sin( %2 ) Lcξ92ξ10

−2 i σ2m k sin( %2 ) ξ93−4 i m k sin( %2 ) ξ92ξ10.ξ102( %1 −1 ) ( %1 + 1 )

ξ92(L0−Lc) ( L0+Lc )−2ω6πei lnm−√2k

ω(α−α′)m−√2k−m−√2k5

ω4

+mm−√2k4

ω4−Lc2m−√2k3

ω2−3σ2mm−√2k4

ω4−Lc2L02m−√2k

+Lc2mm−√2k2

ω2+L02mm−√2k2

ω2+Lc2L02m−L02σ2mm−√2k2

ω2

−m σ2Lc2m−√2k2

ω2+m σ2L02Lc2−L02m−√2k3

ω2

Appendix 85

cos( %2 ) m−cos( %2 ) m−√2k+ i sin( %2 ) √2k.( %1 −1 ) ( %1 + 1 ) k

2k2−2√2k m +m2+L02ω222k2−2√2k m +m2+Lc2ω22+ 2ω6π

ei lnm+√2k

ω(α−α′)m+√2kLc2L02m+m σ2L02Lc2+Lc2mm+√2k2

ω2

−L02m+√2k3

ω2−Lc2L02m+√2k−3σ2mm+√2k4

ω4−m+√2k5

ω4

−m σ2Lc2m+√2k2

ω2−L02σ2mm+√2k2

ω2+mm+√2k4

ω4

−Lc2m+√2k3

ω2+L02mm+√2k2

ω2

cos( %2 ) m+ i sin( %2 ) √2k−cos( %2 ) m+√2k.( %1 −1 ) ( %1 + 1 ) k

2k2+ 2 √2k m +m2+L02ω222k2+ 2 √2k m +m2+Lc2ω22πeσ2



Γ3

4+1

2m+1

2iα′

2

%1 := e(π(α−α′) )

%2 := karctan 1

2σ2√2m(α−α′),1

2σ2α′2+ 2 m2+ ( α−α′)α′

where

ξ7 = L02ω2+ 2 i m ω L0−m2+ 2 k2

ξ8 = L02ω2−2 i m ω L0−m2+ 2 k2

ξ9 = Lc2ω2+ 2 i m ω Lc−m2+ 2 k2

ξ10 = Lc2ω2−2 i m ω Lc−m2+ 2 k2

86 References

9. Referenes

Abramowitz, M, and Stegun, I.A., 1968: Handbook of Mathematical Functions, New

York: Dover Publishing [1968hmfw.book.....A]

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

338

, 78

[1989ApJ...338...78B]

Bertin, G., and Lin, C.C., 1996: Spiral Structure in Galaxies: A Density Wave theory,

Cambridge, Mass.: The MIT Press [1996ssgd.conf.....B]

Binney J., Tremaine S., 1987: Galactic Dynamics, Princeton: Princeton Univ. Press

[1987gady.book.....B]

Biviano, A., Girardi, M., Giuricin, G., Mardirossian, F., and Mezzetti, M., 1991: ApJ

376

, 458 [1991ApJ...376..458B]

Block, D.L., and Wainscoat, R.J., 1991: Nature

353

, 48 [1991Natur.353...48B]

Block, D.L., Bertin, G., Stockton, A., Grosbøl, P., Moorwood, A.F.M., Peletier, R.F. 1994:

A&A

288

, 365 [1994A&A...288B]

Bogoliubov N.N., Mitropolski Y.A., 1961: Asymptotic Methods in the Theory of Non-

Linear Oscillations, New York: Gordon and Breach

Burton, W.B., 1971: A&A

, 76 [1971A&A....10...76B]

Clutton-Brock, 1972: Ap&SS

, 292 [1972Ap&SS..17..292C]

Evans, N.W., and Read, J.C.A., 1998a: MNRAS

300

, 83 [1998MNRAS.300...83E]

Evans, N.W., and Read, J.C.A., 1998b: MNRAS

300

, 106 [1998MNRAS.300..106E]

Goodman, J., and Evans, N.W., 1999: MNRAS,

309

, 599 [1999MNRAS.309..599G]

Fuchs, B, 1991: In: Sundelius, B. (ed.), Dynamics of Disc Galaxies, 359

[1991ddg..conf..359F]

Gakhov, F.D., 1990: Boundary value Problems, Mineola: Dover Publications (based on

the the second Russian edition published in 1963)

Goldreich, P., and Lynden-Bell, D., 1965: MNRAS

130

, 125 [1965MNRAS.130..125G]

Gradshteyn, I.S., and Ryzhik, I.M., 1994: Table of Integrals, Series, and Products, New

York: Academic Press [1994tisp.book.....G]

Hockney, R.W., and Hohl, F., 1969: ApJ

, 1102 [1969AJ.....74.1102H]

Jeans, J. H., 1919, Problems of Cosmology and Stellar Dynamics, Cambridge: Cambridge

University Press [1919QB981.J35......]

Julian, W.H., and Toomre, A., 1966: ApJ

146

, 810 [1966ApJ...146..810J]

Jungwiert, B., and Palouˇs, J., 1994: A&A

287

, 55 [1994A&A...287...55J]

Kalnajs, A.J., 1971: ApJ

166

, 275 [1971ApJ...166..275K]

Kalnajs, A.J., 1976a: ApJ

205

, 745 [1976ApJ...205..745K]

Kalnajs, A.J., 1976b: ApJ

205

, 751 [1976ApJ...205..751K]

Kalnajs, A.J., 1977: ApJ

212

, 637 [1977ApJ...212..637K]

Kuno, N, Tosaki, T, Nakai, N., and Nishiyama, K., 1997: PASJ

, 275

[1997PASJ...49..275K]

References 87

Kuypers F., 1990: Klassische Mechanik, Weinheim: VCH

Landau, L.D., 1946: J. Phys. USSR

, 25 [reprinted in: Ter Haar, D., 1969: Men of

Physics: L.D. Landau II, Oxford: Pergamon Press]

Lemos, J.P.S., Kalnajs, A.J., and Lynden-Bell, D., 1991: ApJ

375

, 484

[1991ApJ...375..484L]

Lin, C.C., and Shu, F., 1964: ApJ

140

, 646 [1964ApJ...140..646L]

Lindblad, B, 1926: Uppsala Medd. No. 13 [1926MeUpp..13....1L]

Lindblad, B, 1940: ApJ

, 1 [1940ApJ....92....1L]

Lynden-Bell, D., and Ostriker, J.P., 1967: MNRAS

136

, 87P [1967MNRAS.136P..87L]

Lynden-Bell, D., and Kalnajs, A., 1972: MNRAS

157

, 1 [1972MNRAS.157....1L]

Mark, J. W.-K., 1976: ApJ

205

, 363 [1976ApJ...205..363M]

Mestel, L., 1963: MNRAS

126

, 553 [1963MNRAS.126..553M]

Miller, R.H., and Smith, B.F., 1992: ApJ

393

, 508 [1992ApJ...393..508M]

Muskhelishvili, N. I., 1953: Singular integral equations, Groningen: Noordhoff (based on

the second Russian edition published in 1946)

Naim, A., Lahav, O., Sodre, L. jr., Storrie-Lombardi, M.C., 1995: MNRAS

275

, 567

[1995MNRAS.275..567N]

Persic, M. and Salucci, P., 1991: ApJ bf 368, 60 [1991ApJ...368...60P]

Pipkin, A.C., 1991: A course on integral equations, New York: Springer

Press, W.H, Teukolsky, S.A., Vetterling, W.T., and Flannery, B.P., 1992: Numerical Re-

cipies in C, Cambridge, New York: Cambridge University Press [1992nrca.book.....P]

Polyanin, A.D., and Manzhirov, A.V., 1998: Handbook of integral equations. Boca Raton:

CRC Press

Read, J.C.A., 1997: The Stability of Model Disk Galaxies, Ph.D. Thesis, University of

Oxford

Read, J.C.A., 1999: private communication

Rix, H.-W., and Rieke, M.J., 1993: ApJ

418

, 123 [1993ApJ...418..123R]

Rohlfs, K., 1977: Lectures on Density Wave Theory, Berlin, Heidelberg, New York:

Springer [1977ldwt.book.....R]

Rosse, William Parsons Earl of, 1852: Drawings to illustrate Recent Observations on

Nebulæ, in: Parsons, Charles, 1926: The scientific Papers of William Parsons, Third

Earl of Rosse 1800-1867, Bradford and London: Percy Lund, Humphreys&Co

Rots A.H., 1975: A&A

, 43 [1975A&A....45..43R]

Rubin, V.C., Thonnard, N., Ford, W.K., and Burstein, D., 1982: ApJ

261

, 439

[1982ApJ...261..439R]

Sellwood, J.A., and Carlberg, R.G., 1984: ApJ,

282

, 61 [1984ApJ...282...61S]

Sellwood, J.A., 1985: MNRAS

217

, 127 [1985MNRAS.217..127S]

Taga, M., 1998: In: ASP Conf. Ser.

182

, 71 [1998ASPC..182...71T]

Thornley, M.D., and Mundy, L.G., 1997: ApJ

490

, 682 [1997ApJ...490..682T]

Toomre, A., 1964: ApJ

139

, 1218 [1964ApJ...139.1218T]

Toomre, A., 1969: ApJ

158

, 899 [1969ApJ...158..899T]

Toomre, A., 1981: In: The structure and evolution of normal galaxies, Cambridge, UK,

New York: Cambridge University Press, 111 [1981seng.proc..111T]

Toomre, A., and Toomre, J., 1972: ApJ

178

, 623 [1972ApJ...178..623T]

von Hoerner,S., 1963: Zs. f. Astrophysik

, 47 [1963ZA.....57...47V]

88 References

Wladimirow, W.S., 1972: Gleichungen der mathematischen Physik, Berlin: VEB

Deutscher Verlag der Wissenschaften (Hochschulb¨ucher f¨ur Mathematik, Band 74)

Zang, T.A., 1976: The Stability of a Model Galaxy, Ph.D. Thesis, Massachussetts Insti-

tute of Technology, Cambrigde, MA