On stationary oscillations of galaxies

Abstract

This paper describes a general investigation of stationary oscillations of galaxies. It begins with a linear analysis of modes of oscillation with continuous spectra of real frequencies. Such modes are gravitational analogues of the van Kampen modes of oscillation in plasmas. The characteristic value problem governing modes of the van Kampen type in a galaxy is solved with the aid of a modified version of the matrix method of Kalnajs in which the perturbation of the distribution function is expressed in terms of generalized functions. In general, there is no characteristic equation governing the frequencies in the continuous spectrum. However, isolated frequencies in the continuous spectrum do satisfy a characteristic equation, which, for stellar systems, is a counterpart of the dispersion relation proposed by Vlasov for plasma oscillations. The linear analysis also provides a characteristic equation for modes with a discrete spectrum of real and/or complex frequencies. The second part of the paper describes a perturbation theory for a stationary oscillation of a galaxy with a small but finite amplitude. Integrals of the stellar motion are constructed with the aid of canonical perturbation theory and used in conjunction with Jeans' theorem to specify the density of stars in the six-dimensional phase space. These oscillations are slightly non-linear counterparts of the modes of the van Kampen type, and they are stellar-dynamical counterparts of the non-linear plasma waves described by Bernstein, Greene and Kruskal in the 1950s. Fully non-linear models of stationary oscillations of galaxies can be constructed with the aid of Schwarzschild's numerical method for the solution of the fundamental integral equation describing the self-consistency of a stellar system.

instabilities, plasmas, stellar dynamics, galaxies: kinematics and dynamics, galaxies: structure

1 Introduction

Stationary oscillations of galaxies are oscillations with constant amplitudes. Relatively little attention has been given to stationary modes of oscillation in conventional investigations of small perturbations of galaxies. A major reason for this apparent neglect is that the principal investigations of the subject have been based on methods that have their origin in Landau's initial-value formulation of the theory of plasma oscillations (Landau 1946). A particular feature of such methods is the adoption of the so-called ‘Landau contour’ or ‘Landau prescription’ for the integration of singular quantities over the momentum space (Fridman & Polyachenko 1984; Palmer 1994).

Early investigations of the oscillations and the stability of infinite, homogeneous stellar systems (Lynden-Bell 1962; Sweet 1963) closely follow Landau's formulation and solution of the initial-value problem for plasma oscillations. Landau's prescription also underlies the matrix method formulated by Kalnajs (1977) and used by him and by many others (e.g. Polyachenko & Shukhman 1981; Fridman & Polyachenko 1984; Palmer & Papaloizou 1987; Weinberg 1989, 1991a,b; Saha 1991; Bertin et al. 1994; Palmer 1994; Robijn 1995) for investigations of instabilities in galaxies. The investigations by Zang (1976) and by Evans & Read (1998a,b) of coplanar oscillations and instabilities of thin discs in terms of a homogeneous integral equation of the Fredholm type also require Landau's prescription (see also Kalnajs 1971). Although there is considerable variation in the details of the methods employed, all of these investigations are based, at least implicitly, on initial-value formulations for the solution of the governing equations. In addition, they all reduce to the investigation of oscillations of galaxies with complex frequencies or, where methods of Laplace transforms are used (e.g. Kalnajs 1971; Weinberg 1989, 1991a,b), to the investigation of superpositions of such oscillations. Initial conditions are imposed in such a way that the reduction and solution of the governing perturbation equations are valid only in the case that the imaginary parts of the frequencies have algebraic signs such that the amplitudes of the oscillations grow exponentially with time. Landau's prescription for the interpretation and evaluation of integrals of singular quantities over the phase space follows from the construction of such solutions for exponentially growing oscillations and the analytic continuation of those solutions into other regions of the complex-frequency plane (see e.g. Jackson 1960; Stix 1992; Palmer 1994).

Such methods are suitable for the investigation of instabilities of galaxies. However, in place of such initial-value formulations, one needs a normal-mode analysis, in the sense of van Kampen (1955, 1957), in order to study stationary oscillations of galaxies within the framework of linear perturbation theory.

The understanding and expectations that one has regarding small perturbations in stellar systems rest heavily on what is known about plasma oscillations. In particular, it is well known (see e.g. Jackson 1960; Clemmow & Dougherty 1969; Ecker 1972; Stix 1992) that the initial-value formulation of Landau (1946) and the normal-mode analysis of van Kampen (1955, 1957) provide mutually complementary descriptions of electrostatic waves in a homogeneous plasma. The van Kampen modes with continuous spectra of real frequencies are stationary modes of oscillation. However, those modes are singular in the sense that perturbations of the distribution of particles in the six-dimensional phase space of a single particle are expressed in terms of generalized functions. It is frequently suggested that such distributions are artificial and would be difficult to establish as initial conditions in a plasma. On the other hand, the van Kampen modes are complete (Case 1959), and they may be superposed in order to represent more physical initial conditions. As results of dispersion and phase mixing, such superpositions evolve consistently with Landau's solution of the initial-value problem and exhibit forms of apparently irreversible behaviour such as Landau damping. Conversely, when a van Kampen mode is adopted as an initial condition, the solution of the initial-value problem in accordance with Landau's prescription is that van Kampen mode. Thus, the study of the van Kampen modes, the linear modes of stationary oscillation in a plasma, contributes significantly to our understanding of plasma oscillations, notwithstanding that such modes might be difficult to excite individually and superpositions of such modes would suffer Landau damping (see especially section 8-10 in Stix 1992).

The literature of plasma physics shows that van Kampen modes are significant and important in other respects as well. The van Kampen modes are the small-amplitude limits of non-linear plasma waves (‘BGK waves’) of the kind described by Bernstein, Greene & Kruskal (1957) (see also Bohm & Gross 1949; Jackson 1960; Stix 1992). The non-linearity of BGK waves precludes a principle of superposition. Thus, BGK waves are stationary oscillations of a homogeneous plasma, and they do not suffer Landau damping. As the linear limit of a BGK wave, each van Kampen mode must be considered individually. The interpretation in terms of a BGK wave shows how a van Kampen mode owes its self-consistency to the way in which orbits trapped in the electrostatic potential of the wave are populated with particles. The connection also establishes the van Kampen modes as a starting point in linear perturbation theory for the study of a certain class of non-linear waves.

For stellar systems, this paper describes an investigation of stationary oscillations, which are the gravitational counterparts of van Kampen modes and BGK waves.Section 2 describes the construction of modes of oscillation and instability of a stellar system within the framework of linear perturbation theory and with the aid of a modified version of the matrix method of Kalnajs (1977) (see Fridman & Polyachenko 1984; Palmer 1994). In particular, this formulation describes singular modes with continuous spectra of real frequencies. The treatment unites the methods of Kalnajs and van Kampen and generalizes both. Section 3 describes the construction of slightly non-linear oscillations of stellar systems. In order to solve the governing equations, we make use of canonical perturbation theory (Goldstein 1980) in order to construct suitable integrals of the motion for the stellar orbits, and we apply Jeans' theorem in order to specify the distribution of stars in the six-dimensional phase space. The oscillations derived in this way are slightly non-linear counterparts of the singular modes of the van Kampen type that are described in Section 2. As in the case of BGK waves, such stationary oscillations of a galaxy owe their self-consistency to the way in which resonant orbits are populated with stars. The non-linear treatment imposes certain constraints on the frequencies of modes of the van Kampen type which seem not to arise within the framework of linear perturbation theory. Section 4 contains a very brief discussion of fully non-linear oscillations of stellar systems, and the paper concludes in Section 5 with a summary and discussion of the important results of this work. For the sake of clarity and specificity, we include Appendix A in which the general analysis described in Sections 2 and 3 is worked out explicitly for plane waves in an infinite, homogeneous system of stars. Although an important point of this paper is to emphasize the distinction between initial-value formulations and normal-mode formulations of methods for the study of the oscillations and the stability of galaxies, an equally important point is to recognize the complementarity of the two approaches and to identify features that they have in common. Accordingly, Appendix B contains a comparison of the representations of perturbations in the conventional initial-value formulation and the present normal-mode formulation of the matrix method of Kalnajs.

Remarkably, many of the elements of what is fully worked out in the present investigation have been anticipated by Louis & Gerhard (1988), who asked the question ‘Can galaxies oscillate?’ and then constructed a self-consistent model of a spherically symmetric stellar system in a state of stationary oscillation. They constructed that model with the aid of a suitably generalized version of the Schwarzschild (1979) method for the solution of the fundamental integral equation describing the self-consistency of a stellar system. That paper provides an important example of a fully non-linear (albeit numerical) model of a galaxy in a state of stationary oscillation. The work of Louis & Gerhard and its relationship to the present investigation are the main subjects of Section 4. It is relevant to the motivation for the present investigation to note that Louis & Gerhard (1988) begin their paper with a brief review of observational evidence that states of equilibrium may not be the only ‘normal’ states of galaxies. They suggest that some observed galaxies could be in states of stationary oscillation.

The generality of the present investigation should be emphasized. In the first place, what follows establishes for finite, inhomogeneous stellar systems many results that have been well established for infinite, homogeneous plasmas and stellar systems. Moreover, these results apply to a wide class of stellar systems. The only significant restriction is that we consider the stationary oscillations of a system in which the motions of stars in the unperturbed gravitational potential are integrable.

2 Modes Of Oscillation Of The Van Kampen Type

Infinitesimal perturbations in a galaxy are governed by the collisionless Boltzmann equation in the linearized form  

1

where the distribution function ƒ₀(x, v) represents the unperturbed density of stars in the six-dimensional phase space of a single star, V₀ (x) denotes the unperturbed gravitational potential of the system, and ƒ₁(x, v, t) and V₁(x, t) denote the perturbations of the distribution function and the potential, respectively. Here, x and v denote the position and velocity of a star, respectively, and t denotes the time. The density of the perturbation  

2

where m_* is the mass of a star, is the source of the potential V₁(x, t). In other words,  

3

where the integration extends over the volume V of the system.

The normal modes of the system are described by solutions of equations (1)–(3) for f₁(x, v, t), ρ₁(x, t) and V₁(x, t) with a time dependence exp(−iσt), where the characteristic frequency σ is a constant.

2.1 Solution for \mathbf{\mathit{f}}_{1}(\mathbf{\mathit{x}},<?TeX\bv?>,\mathbf{\mathit{t}}) in terms of action-angle variables

We assume that the motion of a star in the unperturbed potential V₀(x) is integrable, i.e. that there exist three mutually independent isolating integrals of the unperturbed motion. Accordingly, with the aid of a suitable canonical transformation (x, v) → (w, J), we introduce a set of angle variables w = (w₁, w₂, w₃) and conjugate action variables J = (J₁, J₂, J₃) which describe the motion of a star in the unperturbed system. When expressed as functions of the action-angle variables, functions of x and v are quasi-periodic functions of the angle variables w with periods 2π.

The actions J are isolating integrals of the motion of a star in the gravitational field of the unperturbed system. Therefore, the unperturbed distribution function can be expressed as a function ƒ₀(x, v) = ƒ₀(J) of the actions in accordance with Jeans' theorem.

Transformed to action-angle variables, the linearized, collisionless Boltzmann equation now reads  

4

where the quantities ω(J) = (ω₁, ω₂, ω₃) are the frequencies of the unperturbed motion of a star. In other words, ω(J) = ∂H₀/∂J, where H₀ (J) is the Hamiltonian that governs the unperturbed motion.

In order to solve equation (4), we express the perturbations of the distribution function and the gravitational potential as trigonometric series in the angle variables of the forms  

5

and  

6

respectively, where the sums extend over sets of integers n = (n₁, n₂, n₃). In this representation, trigonometric components of equation (4) separate, and the coefficients f₁(n, J, σ) and V₁(n, J, σ) satisfy the equation  

7

The solution of equation (7) for f₁(n, J, σ) must allow the integration of quantities involving f₁(x, v, t) over the entire region of the phase space that is accessible to stars in the unperturbed system. In particular, f₁(n, J, σ) must be defined for this purpose at values of J such that σ−n·ω(J) = 0. For each set of integers n, that condition defines a surface of two dimensions in the three-dimensional space of the actions J. Let two variables (ϖ₁, ϖ₂) = ϖ, say, label points on that surface. Then the values of the actions at a given point ϖ on the surface are given by the functions J = J_R (n, ϖ, σ) (say).

The solution for f₁(n, J, σ) is a generalized function of the form  

8

where P signifies that the integral of a quantity must be interpreted as the Cauchy principal value, λ (n, ϖ, σ) is an arbitrary function, and δ denotes Dirac's delta function. The integral on the right-hand side of equation (8) extends over the entire surface in the action space on which σ−n·ω(J) = 0. Equation (8) expresses f₁(n, J, σ) as a sum of a particular solution of equation (7) and a general solution of the reduced equation [σ−n·ω(J)] f₁(n, J, σ) = 0. In particular, λ (n, ϖ, σ) δ[ J−J_R (n, ϖ, σ) ] d ϖ is the contribution to the solution of the reduced equation (for fixed values of n and σ) at values of the actions in the neighbourhood of the values J = J_R(n, ϖ, σ). The integration over ϖ adds such contributions over all points J on the surface σ−n·ω(J) = 0.

Equation (8) clearly satisfies equation (7) where σ−n·ω(J) ≠ 0 and J ≠ J_R(n, ϖ, σ). Where σ−n·ω(J) = 0, equation (8) also satisfies equation (7) in the sense that  

Here the integral over ϖ vanishes in the limit, because  

10

Substituting from equations (8) for the coefficients ƒ₁(n, J, σ) in equation (5), we obtain the solution for the perturbation of the distribution function in the form  

11

It is in the solution of equation (7) that the present formulation of the matrix method differs from the formulation by Kalnajs (1977). Where the solution by Kalnajs is based on Landau's initial-value analysis of plasma oscillations (Landau 1946), the present solution is based on van Kampen's normal-mode analysis (van Kampen 1955). The difference is manifest in the interpretation of equations (8) and (11) in terms of Cauchy principal values of integrals and in the appearance in those equations of Dirac's delta function. Thus, for example, equation (11) is to be contrasted with equation (10) in section 4 of the appendix of Fridman & Polyachenko (1984) or with equations (9.5) and (9.6) in Palmer (1994). The contrast is considered further in Appendix B of this paper.

2.2 A matrix method for the solution of the characteristic value problem governing normal modes

Substituting from equation (11) for f₁(x, v, t) in equation (2) and rearranging terms in the resulting equation, we obtain  

12

Equation (12) is a condition of self-consistency to be satisfied by perturbations of the system, inasmuch as the density ρ₁(x, t) must be the source of the gravitational potential V₁(x, t). In other words, the characteristic value problem governing the normal modes of the system now reduces to equations (3) and (12). We can solve the characteristic value problem with the aid of a modified version of the matrix method of Kalnajs (1977).

Let ρ₁(x, β) and V₁(x, β) be members of a biorthonormal set of densities and potentials, in the sense of Clutton-Brock (1972) and Kalnajs (1976), where the constants β are sets of parameters which label the members of that basis set, and each density ρ₁(x, β) is the source of the corresponding potential V₁(x, β). A biorthonormal set of densities and potentials satisfies the conditions  

13

where the integration extends over the volume V of the configuration, and δ(α, β) is the Kronecker delta. The minus sign on the right-hand side of equation (13) is required by the physical interpretation of the integral on the left-hand side as the gravitational potential energy of the mass distribution described by the density ρ₁(x, β) in the gravitational potential V₁^*(x, α). When β = α, the integral represents the self-energy of the perturbation, which is negative.

The matrix formulation of the characteristic value problem is based on representations of the perturbations of the density and the potential of the forms  

14

respectively, where the coefficients μ(β, σ) are constants to be determined. We shall need representations of the potentials V₁(x, β) as trigonometric series in the angle variables of the form  

15

It follows from the second of equations (14) that the coefficients in the trigonometric series on the right-hand sides of equations (6) and (15) must be related in the manner  

16

We obtain a matrix representation of the condition of self-consistency by multiplying each term in equation (12) by V₁^*(x, α) and integrating the product over the volume of the configuration space accessible to stars in the unperturbed system. In the reduction of the resulting equation, we transform integrals over the phase space from integrals over x and v to integrals over w and J. Inasmuch as this is a canonical transformation of coordinates and momenta, the Jacobian determinant of the transformation is unity. Making use of 13141516, and noting that integrals of terms of the form exp[i(n−n′) ·w] over the angle variables vanish unless n′ = n, where n and n′ are sets of integers, we reduce the resulting equation to  

17

where we have suppressed a common factor exp(−iσt).

2.3 Normal modes

Equation (17) represents a system of algebraic equations of the form  

18

for the determination of the coefficients μ(β, σ), where  

19

and  

20

2.3.1 The continuous spectrum of modes

As they stand, equations (18)19(20) provide a matrix representation of the characteristic value problem governing modes with a continuous spectrum of real frequencies. For an assigned value of σ in the continuous spectrum, equations (18) are an inhomogeneous system of linear equations governing the coefficients μ(β, σ). In general, those equations admit of a solution only if det∣M(σ, α, β)∣≠ 0.

Without loss of generality, we can normalize the perturbation in such a way that  

21

inasmuch as the gravitational potential energy of the perturbation must be negative. In virtue of equations (13) and (14), this normalization reduces to  

22

For the assigned value of σ, equation (22) imposes one constraint on the quantities Λ(σ, α) and, hence, on the functions λ(n, ϖ, σ) (see equations 20). Inasmuch as the functions λ(n, ϖ, σ) are arbitrary, apart from that constraint, there is no characteristic equation for the determination of the value of σ. Moreover, the modes belonging to an assigned value of σ in the continuous spectrum are highly degenerate.

The physical meaning of the degeneracy of the modes in the continuous spectrum is the following. The functions λ(n, ϖ, σ) determine the perturbation f₁(x, v, t) on the surfaces in the action space on which σ−n·ω(J) = 0 and the unperturbed stellar orbits are resonant. The arbitrariness of the functions λ(n, ϖ, σ) represents a freedom to populate the resonant orbits in different ways in order to maintain the self-consistency of the perturbation.

One particular construction of modes of the van Kampen type can be achieved by choosing λ(n, ϖ, σ) in the manner  

23

where n₀ and ϖ₀ are fixed constants, δ(n, n₀) denotes a Kronecker delta, and δ(ϖ−ϖ₀) denotes Dirac's delta function. With this choice, we populate only the resonance n = n₀ with stars, and, on the resonant surface σ−n₀·ω(J) = 0 in the action space, we include only stars whose actions have the values J = J_R(n₀, ϖ₀, σ). In this case, equation (20) reduces to  

24

Therefore, all of the coefficients μ(β, σ) in the solution of equation (18) are proportional to λ(n₀, ϖ₀, σ). It follows that the normalization of the coefficients specified by equation (22) determines the value of |λ(n₀, ϖ₀, σ)|².

In Section 3, we shall find that non-linear effects impose a constraint on the frequencies of the modes of the van Kampen type. Specifically, we shall find that the spectrum of frequencies is bounded by the frequency of a mode of the ‘Vlasov type’.

2.3.2 Modes of the Vlasov type

At exceptional frequencies in the continuous spectrum, it can happen that det ∣M(σ, α, β)∣ = 0. This is essentially a characteristic equation for the determination of such exceptional values of σ. Under these conditions, a solution for the coefficients μ(β, σ) exists if and only if Λ(σ, α) = 0 in equation (18). In other words the quantities λ(n, ϖ, σ) must all vanish. Nevertheless, the integrals that appear in the definition of M (σ, α, β) in equation (19) must be interpreted as Cauchy principal values. In this case, the characteristic equation det|M(σ, α, β)|= 0 is the counterpart of the dispersion relation proposed by Vlasov (1945) for plasma oscillations. This is a mode in the continuous spectrum in which stars on resonant orbits do not contribute to the self-consistency of the perturbation. The interpretation of Vlasov's dispersion relation in terms of an absence of particles trapped in a plasma wave has been discussed by Bohm & Gross (1949) and Jackson (1960). The interpretation of modes of the Vlasov type is illuminated by the discussion of their slightly non-linear counterparts at the end of Section 3 below.

2.3.3 The discrete spectrum of modes

In the cases of modes with complex frequencies and of certain modes with real frequencies, σ−n·ω(J) ≠ 0 for all values of n and all values J accessible to stars in the unperturbed system. In such cases, there are no terms in the sums over integrals involving the functions λ(n, ϖ, σ) on the right-hand sides of equations (12), (17) and (20). In other words, the quantities Λ(σ, α) vanish identically, and equations (18) reduce to  

25

In this case, it is not necessary to interpret the integrals on the right-hand side of equation (19) as Cauchy principal values. Solutions of equations (25) for the coefficients μ(β, σ) exist if and only if det ∣M(σ, α, β)∣ = 0. This condition provides a characteristic equation for the determination of the frequencies σ of the discrete modes.

In Appendix B, we show that the modes of instability derived from equation (25) are the unstable perturbations that would be found in applications of conventional versions of the matrix method of Kalnajs.

3 Slightly Non-Linear Oscillations Of A Galaxy

In the theory of modes of the van Kampen type described in the preceding section, we have attributed singularities in the perturbation of the distribution function to contributions of stars on resonant orbits. In the case of plasma oscillations, such an interpretation of the singularities in the van Kampen modes has been verified by Bernstein et al. (1957) in their investigation of an exact model of non-linear plasma waves (see also Stix 1992).

In this section, we investigate stationary oscillations of galaxies in cases in which the amplitudes of the oscillations are small but finite. As in the investigation of Bernstein et al. (1957), we find that such oscillations can be represented in terms of series in half-integral powers of the amplitude of the perturbation of the gravitational potential. We solve the governing equations through terms of the first order in that amplitude. The solution requires an explicit treatment of the contributions of resonant orbits to the density of stars in the six-dimensional phase space. This treatment of stationary oscillations confirms and clarifies the interpretation described above of the singularities in the modes of the van Kampen type.

3.1 Formulation of the problem and reduction in terms of action-angle variables

Consider an oscillating galaxy in which the density and the gravitational potential are of the forms  

26

and  

27

respectively, where σ is a constant, real frequency and the small, dimensionless parameter ϵ is a positive, real number. We require that the oscillation be self-consistent in the sense that the density ρ(x, t) is the source of the potential V(x, t). Evidently, the density ρ₀ (x) must be the source of the gravitational potential V₀(x). We shall find that ρ₀(x) and V₀(x) characterize an unperturbed state of the galaxy. Let f₀(x, v) denote the distribution function of the galaxy in that unperturbed state. In equations (26) and (27), we represent ρ(x, t) and V(x, t), respectively, as real quantities, expressed, for the sake of convenience in what follows, in terms of a complex density ρ₁(x) and a complex potential V₁(x), respectively. The requirement of self-consistency also implies that the density ρ₁(x) must be the source of the potential V₁(x).

The goal in what follows is to construct the distribution function ƒ(x, v, t) of the galaxy in the state of oscillation described by equations (26) and (27). In order to accomplish this, we investigate the motions of stars in the gravitational potential described by equation (27), we construct a set of integrals of the motion, and, with the aid of Jeans' theorem, we express ƒ(x, v, t) as a function of those integrals. We carry out this programme under conditions in which the oscillatory perturbation of the system is small but finite, and, in the construction of the distribution function, we evaluate all quantities involved consistently through terms of the first order in the parameter ϵ. We have written the expressions for ρ(x, t) and V(x, t) as we have in terms of ϵ in order to facilitate the ordering of terms in the series that arise in this process.

As in Section 2, we consider a system in which the motion of a star in the unperturbed potential V₀(x) is integrable, and we introduce the actions J and angles w of the unperturbed motion as canonical momenta and coordinates, respectively. The Hamiltonian is expressible as a function H = H₀(J) of the actions alone. It follows from the canonical equations of motion in that case that the actions J are integrals of the unperturbed motion of a star and the frequencies of the unperturbed motion are given by the relations ω(J) = ∂H₀/∂J.

We now express the function V₁(x) as a function of the action-angle variables and particularly as a trigonometric series in the angle variables. Thus, we write  

28

where the sum extends over sets of integers n = (n₁, n₂, n₃). The motion of a star in the potential V(x, t) is governed by the Hamiltonian  

29

3.2 Construction of the integrals of the motion for a resonant orbit

Consider now the motion of a star in the neighbourhood of the resonance at which σ−m·ω(J) = 0, where m = (m₁, m₂, m₃) is a particular set of three integers. We can construct the integrals of the motion through order ϵ in this case with the aid of two successive canonical transformations. This procedure is essentially von Zeipel's method for the construction of a canonical perturbation theory (von Zeipel 1916; Contopoulos 1975; Goldstein 1980).

The first of these transformations replaces the actions J and angles w with new momenta I = (I₁, I₂, I₃) and new coordinates ψ = (ψ₁, ψ₂, ψ₃) with the aid of the generating function  

30

a function of the new momenta and the old coordinates. The transformation equations are  

31

and  

32

and the new Hamiltonian is  

33

(Goldstein 1980). In equation (33), we have eliminated the old angles w in favour of the new angles ψ with the aid of equation (32). This canonical transformation removes the non-resonant terms of order ϵ from the perturbed Hamiltonian.

The second canonical transformation introduces final momenta G = (G_R, G₂, G₃) and final coordinates θ = (θ_R, θ₂, θ₃), and it is designed to make the quantity m·ψ−σt a canonical coordinate. The generating function is  

34

again a function of the new momenta and the old coordinates. The transformation equations are  

35

and  

36

and the final Hamiltonian is  

37

(Goldstein 1980), where we have defined G_⊥ = (0, G₂, G₃) and we have abbreviated equations (36) to I = mG_R+G_⊥.

If we neglect terms of order ϵ² and higher on the right-hand side of equation (37), then the Hamiltonian described there has the following properties. The new momenta G₂ and G₃ are integrals of the motion, inasmuch as the Hamiltonian does not depend on the coordinates θ₂ and θ₃. Thus H_R(θ_R, G_R, G_⊥) is effectively a Hamiltonian governing a motion in one degree of freedom in the canonical variables (θ_R, G_R). The Hamiltonian H_R(θ_R, G_R, G_⊥) itself is a third integral of the motion inasmuch as it does not depend explicitly on the time.

3.3 Reduction of the perturbed Hamiltonian to a pendulum model

In this subsection, we show that the Hamiltonian of a physical pendulum is a well-defined approximation to the Hamiltonian described by equation (37). As is explained and emphasized by Chirikov (1979) and Henrard (1993), the pendulum model is one of the fundamental paradigms in the study of resonances in non-linear systems. Louis & Gerhard (1988) have made use of a pendulum model of stellar orbits in their interpretation of their investigation of non-linear, radial oscillations of spherical galaxies. Likewise, Lynden-Bell (1973) has made use of a pendulum model in order to represent the behaviour of stellar orbits near the Lindblad resonances of spiral galaxies.

Let G₀ be the ‘resonance value’ of G_R satisfying the condition  

38

for assigned values of the components G₂ and G₃ of G_⊥. We expect, and we shall find, that G_R−G₀ = O(ϵ^1/2) for stars on resonant orbits. By virtue of equation (38), an expansion of the Hamiltonian H_R(θ_R, G_R, G_⊥) in powers of G_R−G₀ accordingly yields, in place of equation (37),  

39

Let  

40

and write the function V₁(m, mG₀+G_⊥) in terms of its amplitude and phase in the manner  

41

Then, equation (39) reduces to  

42

This Hamiltonian has a stable fixed point at (G_R, θ_R) = (G₀, θ₀), where the value of θ₀ is such that cos(θ₀+ϕ) =−η. Inasmuch as the value of η is either 1 or −1 (see equations 40), it follows that an alternative form of equation (42) is  

43

This is essentially the Hamiltonian of a physical pendulum in the canonical momentum ΔG = G_R−G₀ and the conjugate coordinate Δθ = θ_R−θ₀. The separatrix of the system, for an assigned value of G_⊥, satisfies the equation  

44

where we are now ignoring terms of order ϵ^3/2 and higher in the Hamiltonian. In other words, the region of the phase space in which the angle θ_R librates about the value θ₀ and the motion is resonant is defined by the inequality G₋ ⩽ G_R ⩽ G₊, where  

45

This definition of the resonance region confirms the expectation described above that G_R−G₀ = O(ϵ^1/2) for stars in the resonance region. In other words, the size of the resonance region is of order ϵ^1/2. Outside the resonance region, the angle θ_R precesses, and the motion is essentially non-resonant.

The resonance region described by equation (43) for the Hamiltonian and equation (44) for the separatrix consists of one or more islands in a closed chain in the phase space. The island or chain of islands circulates in the phase space. The guiding centre for the circulation of each island is a periodic orbit whose motion is given by J = mG_R+G_⊥+O(ϵ) and m·w−σt+O(ϵ) = θ₀.

The foregoing analysis of the pendulum model described by equation (43) ignores the possibility that the parameter μ vanishes (see equation 40). If, for assigned values of m and G_⊥, the quantity m·ω(mG₀+G_⊥) were not a monotonic function of G₀, then μ would vanish at isolated values of G₀. At such points, the resonance region described by equations (44) and (45) would appear to blow up, but such behaviour would be an artefact of the approximation leading to equation (43) for the pendulum model of the Hamiltonian. A proper treatment of resonances near values of G₀ at which μ vanishes would require a more accurate representation of the Hamiltonian, for example, one in which we explicitly include terms of order ϵ^3/2 and higher on the right-hand side of equation (43). The analysis of that case is beyond the scope of this paper. It should be noted, however, that the solution for the distribution function for the resonant stars which is presented in equation (65) below exhibits no pathologies where μ = 0, provided that the resonance region remains finite (as it should) where μ = 0.

3.4 A matrix representation of a self-consistent oscillation of a galaxy

We are now ready to make use of Jeans' theorem in order to construct the distribution function ƒ(x, v, t) underlying the density distribution and the gravitational potential described by equations (26) and (27), respectively. We shall construct the distribution function in the non-resonance and resonance regions of the phase space separately. Moreover, we shall consider a system in which only a single resonance, defined by the set of integers m = (m₁, m₂, m₃), is populated with stars. Finally, we shall develop the distribution function as a series in powers of ϵ, working consistently to an order of approximation such that the fundamental relation  

46

is satisfied through terms of order ϵ.

In this connection, we recall that we have identified the functions ρ₀(x) and V₀(x) in equations (26) and (27) as the density and gravitational potential, respectively, that characterize an unperturbed state of the galaxy. According to Jeans' theorem, the distribution function ƒ₀(x, v) for the galaxy in that state can be expressed as ƒ₀(J), a function of the unperturbed action integrals. The self-consistency of the equilibrium of the galaxy in that state requires that  

47

3.4.1 Non-resonant stars

We first construct the distribution function for non-resonant stars, that is, for stars outside the region of the phase space defined by the inequality G₋ ⩽ G_R ⩽ G₊, where G_± is defined by equation (45). For the construction of the integrals of the motion required for this purpose, we modify slightly the canonical transformation described by equations (30)31(32). In the trigonometric series in those equations, we now include the terms n = m that were previously omitted. In place of equation (33), the equation for the new Hamiltonian is now K(ψ, I, t) = H₀(I) +O(ϵ²). In this case, the canonical equations of motion imply that the new momenta I are integrals of the motion through terms of order ϵ. According to Jeans' theorem, the distribution function (for the non-resonant stars) can be expressed as a function of the integrals I. Inasmuch as the unperturbed distribution function of the system is of the form f₀(x, v) = f₀(J) described above, we write the perturbed distribution function for the non-resonant stars as  

48

[see equation (31), but recall that we now include the terms n = m in the trigonometric series there].

The physical picture underlying this choice for the non-resonant distribution function is that the non-resonant orbits are orbits that are not altered qualitatively by the perturbation. In particular, the perturbed angle variables ψ precess as do the unperturbed angle variables w. Thus, although the non-resonant orbits are perturbed, the dependence of the distribution function on the integrals of the motion is not altered for the non-resonant orbits.

The contribution of the non-resonant stars to the density of the system is  

49

where the integrals include only regions of the phase space in which the stellar orbits are non-resonant, and P signifies that integrals of functions containing resonance denominators are to be interpreted as Cauchy principal values. We reduce equation (49) to the more useful form  

50

where we have rewritten the integral of f₀(J) with the aid of equation (47) and we have made use of the label R in order to signify that an integral includes only the region of the phase space in which the orbits are resonant. In equation (50), we have also enlarged the region of integration of the last integral in equation (49) to include the region of the phase space containing resonant orbits. Inasmuch as the size of the resonant region is of order ϵ^1/2, the result is to change the integral by a quantity of order ϵ^3/2. That is beyond the order of approximation required here.

The interpretation of integrals in equations (49) and (50) as Cauchy principal values now applies at the resonance n = m as well as at all other resonances. The reduction of equation (49) to equation (50) shows that, in general, the interpretation of the integrals as Cauchy principal values neglects the finite sizes of the resonance regions in the phase space and consequently neglects corrections of order ϵ^3/2 to the density.

3.4.2 Resonant stars

We found in Section 3.2 that the Hamiltonian H_R(θ_R, G_R, G_⊥) and the components G₂ and G₃ of G_⊥ are integrals of the motion in the resonance region of the phase space. According to Jeans' theorem, the distribution function for the resonant stars can be written as  

51

where the function f_R(H_R, G_⊥) remains to the determined (see Section 3.5).

The contribution of the resonant stars to the density is  

52

3.4. A matrix representation of self-consistency

The condition that the oscillation of the galaxy be self-consistent requires that the sum of the densities ρ_NR (x, t) and ρ_R (x, t) be equal to the ‘imposed’ density ρ(x, t). Combining equations (26), (50) and (52), cancelling the term ρ₀(x) on both sides of the resulting equation, and rearranging terms in what is left of the resulting equation, we obtain  

53

Equation (53), which expresses the condition of self-consistency for a slightly non-linear oscillation of the system, is to be compared (apart from terms of order ϵ^3/2 and higher) with equation (12), which expresses the condition of self-consistency for an infinitesimal perturbation of the system. In particular, the right-hand side of equation (53) is essentially the real part of the left-hand side of equation (12). Thus, the essential difference between the two equations lies in the different representations of the contributions of the resonant stars on the right-hand side of equation (12) and on the left-hand side of equation (53).

Beginning with equation (53), we can now formulate the matrix method for the determination of the functions ρ₁(x), V₁(x) and ƒss_R(H_R, G_⊥). In this formulation, we represent ρ₁(x) and V₁(x) in terms of the biorthonormal set of densities and potentials introduced in Section 2.2. In particular, we now let equations (13)14 (16) apply to ρ₁(x, t) = ρ₁(x) exp (− iσt) and V₁(x, t) = V₁(x) exp(−iσt).

We next multiply each term in equation (53) by V₁^*(x, α) exp(iσt), integrate the result over the configuration space, and then average the result over one period 2 π/σ of the oscillation. The reduction here of the terms appearing on the right-hand side of equation (53) follows closely the similar reduction in Section 2.2 of terms appearing on the left-hand side of equation (12) and leading to the terms on the left-hand side of equation (17). In other words, the immediate result of these operations on the terms in equation (53) is  

54

On the left-hand side of equation (54), we transform the integral over the resonance region of the phase space to an integral over the action-angle variables J and w and we recall that the Jacobian determinant of the transformation is unity. Writing V₁^*(x, α) as a trigonometric series in the angle variables in accordance with equation (15), we obtain  

55

We next transform the integral over the actions J on the right-hand side of equation (55) to an integral over the momenta G_R, G₂ and G₃. For this purpose, we make use of the transformation equations (31), (32), (35) and (36), and we observe accordingly that  

In the transformed integral, the integration over G_R is restricted to the resonance region G₋ ⩽ G_R ⩽ G₊, where G_± is given by equation (45). Inasmuch as J₁ = m₁G_R+O(ϵ), the limits of integration over G_R depend on the algebraic sign of m₁. Thus, if m₁ > 0, then the upper and lower limits of integration are G₊ and G₋, respectively. However, if m₁ < 0, then the upper and lower limits of integration are G₋ and G₊, respectively. In order to include both cases in equation (56) below, we replace m₁ dG_R with ∣m₁∣ dG_R and consistently adopt G₊ and G₋ as the upper and lower limits of integration, respectively.

With the further observation that  

and  

we finally reduce equation (55) to  

56

Contributions of order ϵ^1/2 vanish on the right-hand side of equation (56), because they involve integrals over G_R of odd functions of (G_R−G₀) with limits of integration that are symmetric with respect to the value G_R = G₀.

In the next subsection, we shall seek a solution for f_R(H_R, G_⊥) such that  

57

where λ(m, G_⊥, G₀) is a function to be determined. In the remainder of this subsection, we shall complete the formulation of the matrix method. Accordingly, we now combine equations (56) and (57) in order to rewrite the left-hand side of equation (54). Inasmuch as exp(iθ_R) = exp(im·w− iσt)[1 +O(ϵ)] (see equations 32 and 35), we bring equation (54) to the form  

58

Equation (58), with the terms of order ϵ^3/2 and higher suppressed, is an inhomogeneous, linear equation for the determination of the coefficients μ(β, σ). To order ϵ, equation (58) is essentially the matrix representation of self-consistency obtained previously in equation (17) for modes of the van Kampen type. The most substantial difference between equations (17) and (58) is that the inhomogeneous terms on the right-hand side of equation (17) represent a sum of contributions of resonant stars, whereas the inhomogeneous term on the left-hand side of equation (58) represents the contribution of a single resonance. The remaining differences between the two equations are entirely superficial. The components G₂ and G₃ of G_⊥, which label points on the surface σ−m·ω(mG₀+G_⊥) = 0 in the action space, are equivalent to the components ϖ₁ and ϖ₂ of ϖ, which label points on the surfaces σ−n·ω(J) = 0. Likewise, the quantity mG₀ + G_⊥ is equivalent to the function J_R (n, ϖ, σ). The dependence of λ(m, G_⊥, G₀) on G₀ is implicitly a dependence on m, σ and G_⊥ by virtue of equation (38).

Equation (58) is readily generalized in order to include contributions of a number of resonances on the left-hand side. For this purpose, the analysis in Sections 3.2, 3.3 and 3.4.2 would be performed within each of the resonance regions that contain stars, and the results incorporated in equations (53)54555657(58). In each of the resonance regions included, the solution for the perturbation of the distribution function would be determined as described in Section 3.5.

The solution of equation (58) for the coefficients μ(β, σ) can be obtained and the constraint on λ(m, G_⊥, G₀) imposed along the lines described in Sections 2.2 and 2.3 in the case of equation (17). However, in the case of slightly non-linear oscillations, we shall find in the next subsection that the solution of equation (57) for f_R (H_R, G_⊥) imposes additional constraints on λ(m, G_⊥, G₀) which do not arise in the linear theory of modes of the van Kampen type.

3.5 Construction of the distribution function for the resonant stars

We shall now solve equation (57) for f_R(H_R, G_⊥). Let  

59

where F_R(E_R;m, G_⊥, G₀) is the component of f_R(H_R,G_⊥) which remains to be determined. We have defined an ‘energy’  

60

(see equation 43), and we have introduced new canonical variables ΔG = G_R−G₀ and Δθ = θ_R−θ₀. Then we can rewrite equation (57) in the form  

61

It follows from equations (45) and (60) that the integral on the left-hand side of equation (61) must be a function of cos(Δθ) and, hence, an even function of Δθ. Therefore, equation (61) can be satisfied only if  

62

In other words, the quantity λ(m, G_⊥, G₀) exp(i θ₀) must be real, and equation (61) reduces to  

63

Equation (63) can be transformed into an integral equation of the Abelian form. Making use of equation (60) in the inversion of equation (63), we obtain the solution  

64

Upon combining this result with equation (59), we obtain the distribution function for the resonant stars in the form  

65

The term f₀(mG₀ + G_⊥) on the right-hand side of equation (65) is of order unity, and it describes a uniform density of stars in the resonance region of the phase space. Moreover, the contribution of f₀(mG₀+G_⊥) to the density of stars in the resonance region and the contribution of f₀(J) to the density in the non-resonance region (see equation 48) combine in equation (46) to give the unperturbed density ρ₀(x) on the right-hand side of equation (26).

The second term on the right-hand side of equation (65) is formally of order ϵ^1/2, inasmuch as E_R is of order ϵ. That contribution to the distribution function in the resonance region and the contribution of the terms of order ϵ on the right-hand side of equation (48) to the distribution function in the non-resonance region combine in equation (46) to give the sinusoidal contributions to the density on the right-hand side of equation (26).

The second term on the right-hand side of equation (65) is mildly singular as E_R→ 0. Therefore, the requirement that the total density of stars in the six-dimensional phase space must be non-negative implies that  

66

where it will be remembered that λ(m, G_⊥, G₀) exp (i θ₀) is a real quantity (see equation 62).

The singularity of f_R (H_R, G_⊥) at E_R = 0 implies that the resonant stars tend to concentrate near the separatrix defined by equation (44) (cf. equation 60). In other words, the density of stars is low on resonant orbits near the stable periodic orbit, whereas the density is high on resonant orbits near the unstable periodic orbit.

If λ(m, G_⊥, G₀) = 0, then the oscillation described here is a slightly non-linear counterpart of a mode of the Vlasov type. According to equation (58), resonant stars are not required to support the oscillation of the system. Moreover, the mode satisfies the characteristic equation det ∣M(σ, α, β) ∣ = 0 as described in Section 2.3.2. According to equation (65), the density of stars in the resonance region of the phase space is uniform in that case.

3.6 A constraint on modes of the van Kampen type

Inequality (66) constrains the possible frequencies of a stationary oscillation described by equation (58). In order to derive the constraint, we multiply each term in equation (58) by μ*(α, σ), sum over the values of α, and reduce the resulting equation with the aid of equations (16), (22) and (41). Finally, we make use of the condition cos(θ₀ + ϕ) = −η governing the value of θ_R = θ₀ at a stable fixed point of the Hamiltonian represented by equation (42) (see also equation 40). We obtain  

67

It follows from inequality (66) that the product of η and the left-hand side of equation (67) must be negative or zero. Therefore,  

68

The equality holds in condition (68) in the case that λ(m, G_⊥, G₀) = 0. In other words, condition (68) implies that the continuous spectrum of frequencies of the modes of the van Kampen type is bounded by the frequency of a mode of the Vlasov type. The continuous spectrum might be bounded either from above or from below, depending on the sign of η (see equations 40) and on the structure of the integrand in condition (68).

The constraint imposed by condition (68) on the frequencies of plane waves in an infinite homogeneous stellar system is described in Appendix A3.

4 Non-Linear Oscillations Of A Galaxy

In principle, the series in powers of ϵ constructed in Section 3 for the representation of slightly non-linear oscillations could be extended to include terms of higher order in ϵ. However, there does not appear to be any outstanding issue of principle that would require the investigation of terms of higher order. Moreover, the construction of terms of higher order in the series representation of slightly non-linear oscillations of a galaxy would rapidly become complicated. Therefore, it seems preferable to find alternative ways to satisfy equations (26), (27) and (46) for fully non-linear oscillations without resorting to the series solutions started in Section 3.

For a particular model of a radial oscillation of a spherically symmetric galaxy, Louis & Gerhard (1988) have solved equations (26), (27) and (46) with the aid of a generalized version of the Schwarzschild (1979) numerical method for the construction of models of galaxies. Specifically, they have made use of Lucy's method (Lucy 1974) in order to solve the discrete set of equations with which Schwarzschild replaces the integral equation (46).

In their discussion and interpretation of their results, Louis & Gerhard (1988) touch upon certain features of their model which are reproduced more generally in the present work. By investigating the arrangement of resonant orbits in the phase space and exhibiting their solution of Schwarzschild's equations for the distribution of stars in their model, they determine the distribution of stars on resonant orbits that is required in order to sustain the self-consistency of the oscillation. Stars on resonant orbits near the unstable periodic orbit tend to support the oscillation of the model, and those orbits are populated with a relatively high density of stars. Stars on resonant orbits near the stable periodic orbit tend to damp the oscillation of the model, and those orbits are populated with a relatively low density of stars. These are precisely the characteristics of the distribution function described by equation (65), in general, for resonant stars. In the model of Louis & Gerhard, the low density of stars on resonant orbits near the stable periodic orbits produces distinctive ‘gaps’ in the distribution of stars in the phase space. Elsewhere in their paper, Louis & Gerhard present and analyse the pendulum model described in Section 3.3 above in the form appropriate to conditions of spherical symmetry. Their surfaces of section exhibit the arrangement of the resonance regions of the phase space as closed chains of islands as described at the end of Section 3.3 above.

5 Discussion And Concluding Remarks

The work described in this paper generalizes the theory of van Kampen modes in plasmas to the gravitational case of modes of oscillation in inhomogeneous stellar systems, and it does so with the aid of a modified formulation of the matrix method of Kalnajs. The present construction of normal modes, particularly modes with continuous spectra of real frequencies, complements more conventional, initial-value treatments of small perturbations in galaxies.

We have found that the connection between van Kampen modes and non-linear plasma waves carries over to a connection in stellar dynamics between linear modes of stationary oscillation and slightly non-linear, stationary oscillations. That connection also bridges the gap between the linear theory of small perturbations in stellar systems and non-linear models of oscillating systems of the kind considered by Louis & Gerhard (1988). The slightly non-linear theory of stationary oscillations described in Section 3 should provide an analytical foundation for more general and more extensive numerical investigations along the lines described by Louis & Gerhard.

The treatment of stationary oscillations described in this paper is quite general, and it has wide applications. The principal restriction underlying the construction of linear oscillations in Section 2 and slightly non-linear oscillations in Section 3 is that the stellar motions in the unperturbed gravitational field of the system must be integrable. The KAM theorems imply that the series in powers of the small parameter ϵ developed in Section 3 are a valid approximation, for sufficiently small values of ϵ, in the sense that the perturbed motions of stars remain integrable for most initial conditions. Thus, the theories of stationary oscillations of stellar systems described in this paper apply to spherically symmetric systems and to axisymmetric and triaxial Stäckel systems. The KAM theorems also imply that the present representation of stationary oscillations in galaxies would apply, at least approximately, to systems in which the unperturbed stellar motions are only approximately integrable. This would include applications to systems such as slowly rotating, triaxial systems.

On the other hand, for sufficiently large values of ϵ, the overlap of resonances will lead to chaotic behaviour of the relevant stellar orbits (Chirikov 1979). Such chaotic behaviour would be expected to inhibit stationary oscillations at sufficiently large amplitudes. Inasmuch as orbits near the separatrices that form the boundaries of resonance regions would be the first to become chaotic, the onset of chaos in the oscillating systems considered here would probably be enhanced by the tendency of resonant stars to concentrate near the separatrices.

Although the oscillations described in this paper are sinusoidal, it appears that more general periodic oscillations would arise in a more fully non-linear treatment.Louis & Gerhard (1988) showed this to be the case for their model, but they argued that the relevant non-linearities were small enough that they could approximate the oscillation of the model as sinusoidal.

It is an important result of the present investigation that non-linear considerations impose the constraint represented by inequality (68) on the continuous spectrum of modes of the van Kampen type and that a mode of the Vlasov type plays a special role in this connection. The mode of the Vlasov type is the limiting case of modes of the van Kampen type in which resonant stars do not contribute to the self-consistency of the oscillation. Moreover, according to inequality (68), the continuous spectrum of frequencies of modes of the van Kampen type is bounded by the frequency of the mode of the Vlasov type (see particularly the example described in Appendix A3 of this paper). These results, which seem not to be recognized in the literature of plasma physics, raise an important issue regarding the connection between Landau's and van Kampen's treatments of plasma oscillations. The point of inequality (68) is that the excluded modes of the van Kampen type are unphysical in the sense that they involve negative densities of resonant stars in the phase space. One might ask if unphysical van Kampen modes can be included in superpositions of modes that would represent solutions of the initial-value problem of the Landau type. There is certainly no mathematical objection to the use of unphysical modes in a superposition, and there can be no physical objection if the superposition preserves positive densities of stars everywhere in the phase space. Thus, inequality (68) serves mainly to identify those modes of the van Kampen type which, in isolation, may be regarded as physical perturbations.

Stix (1992, see particularly chapters 7 and 8) has emphasized that the connection between van Kampen modes and BGK waves validates the consideration of the van Kampen modes individually as distinct modes of oscillation in a plasma. The results described in this paper serve likewise to validate the consideration of stationary oscillations of galaxies. It remains, however, to identify mechanisms that could excite stationary oscillations of galaxies in nature. That is a question for future research. Nevertheless, one might speculate that, in the formation or growth of a galaxy through mergers, a stationary oscillation could be excited as an accreted fragment becomes trapped in a resonance. In addition, one might ask if dynamical friction acting on an accreted fragment could contribute to such trapping.

This paper has been typeset from a TEX/LATEX file prepared by the author.

Acknowledgments

I warmly thank the referee, whose comments and suggestions significantly increased my understanding of this subject and resulted in a substantial improvement of the paper.

Appendix

Appendix A: Plane Waves In An Infinite, Homogeneous Stellar System

Plane waves in a system in which the density is uniform and the distribution function f₀(v) depends only on the stellar velocities provide concrete illustrations of the linear and slightly non-linear treatments of small perturbations described in the body of this paper. The model is the gravitational counterpart of the theory of plane waves in a homogeneous plasma (Stix 1992).

A Van Kampen Modes

The evolution of small perturbations in an infinite, homogeneous stellar system has been treated as an initial-value problem by Lynden-Bell (1962) and Sweet (1963), and the complementary normal-mode analysis, following van Kampen's treatment of plasma oscillations (van Kampen 1955), is discussed briefly by Binney & Tremaine (1987).

In order to construct a suitable set of action-angle variables and a suitable biorthonormal set of densities and potentials, we employ a ‘box normalization’ and concentrate on the unit cell 0 ⩽x_i⩽L(i = 1, 2, 3). The action-angle variables of the unperturbed motion of a star are  

(A1)

The Hamiltonian and the frequencies of the unperturbed motion are  

(A2)

respectively. The members of the basis set are of the form  

(A3)

where κ = (κ₁, κ₂, κ₃) is a set of three integers and k = 2 πκ/L. Thus the basis functions satisfy periodic boundary conditions on the unit cell. In the orthogonality relations, equations (13), the unit cell is the volume of integration. In particular, the normalization of these bases gives  

(A4)

We must reduce equations (11) and (17) with the aid of equations (A1)A2A3(A4). Upon writing equation (15) explicitly for this purpose, we find that  

(A5)

where δ(n, β) is a Kronecker delta. When we simplify equation (17) with the aid of equation (A5), we find that the modes must be plane waves exp(ik·x− i σt) = exp(iκ·w− i σt). Thus, in equations (11) and (17), we have μ(α, σ) = 0 and λ(α, ϖ, σ) = 0 (α≠κ) and, by virtue of equation (22), μ(κ, σ) = 1. This reduces each of the trigonometric series in equation (11) to a single term in exp(iκ·w).

The condition σ−n·ω(J) = 0 in the case that n = κ reduces to σ−k·v = 0. By resolving v into a component v_k parallel to k and a component v_⊥ perpendicular to k, we can reduce the conditions J = J_R(n, ϖ, σ) to v_k =σ/|k| and v_⊥ = ϖ.

Equation (11) reduces to  

(A6)

Similarly, equation (17) reduces to  

(A7)

Equations (A6) and (A7) describe the gravitational analogue of van Kampen's solution for the electrostatic modes of oscillation of a homogeneous plasma (van Kampen 1955).

For a mode in which the frequency σ is real, equation (A7) fixes the value of the integral of λ(κ, v_⊥, σ) over v_⊥ on the right-hand side, and there is no dispersion relation. Of course, the exceptional case is one in which λ(κ, v_⊥, σ) vanishes, and equation (A7) reduces to the gravitational counterpart of Vlasov's dispersion relation (Vlasov 1945).

When |k|⩽k_J, where k_J is the Jeans wavenumber, and σ−k·v≠ 0, equation (A7) reduces to  

(A8)

a dispersion relation that admits unstable Jeans modes (Binney & Tremaine 1987). If the velocity distribution f₀(v) is Maxwellian, then k²_J = 12πGρ/〈v²〉, where 〈v²〉 is the mean-square velocity of a star in three dimensions.

8.2 A BGK waves

We now construct a gravitational counterpart of a BGK wave (Bernstein et al. 1957) in an infinite, homogeneous stellar system as an example of a slightly non-linear oscillation of the kind described in Section 3.

As in the case of a van Kampen mode of the kind described in Appendix A1, the perturbation of the system consists of a single plane wave. Symbols not defined in what follows have the meanings assigned in the body of the paper or in Appendix A1. Thus the density of the system is of the form  

(A9)

where ρ₀ is the constant density of the system in the absence of the wave, and the gravitational force per unit mass acting on a body is derived from the potential  

(A10)

Here, the wave is described in terms of a single member of the biorthonormal basis set of potentials and densities represented in equations (A3). The Hamiltonian is now of the form  

(A11)

By virtue of (A1) and (A2), the resonance condition reduces to σ−k·v = 0 in this case. The Hamiltonian contains no non-resonant terms. Therefore, the canonical perturbation theory described in Section 3.2 simplifies so that the transformation equations (A1), (31), (32), (35) and (36) reduce to  

(A12)

and we are left with the Hamiltonian  

(A13)

[see equations (37) and (A5) and note that we are now replacing m with κ].

In the reduction of the right-hand side of equation (A13) to the Hamiltonian of a pendulum model, we find that equations (40) reduce to η = +1 and μ = ∣k∣² by virtue of the second of equations (A2) and the first of equations (A12). As in Appendix A1, we introduce components v_k and v_⊥ of v parallel and perpendicular, respectively, to the wavevector k. With the aid of the first of equations (A12), we find that  

(A14)

Thus for the Hamiltonian of the pendulum model in the case at hand, we replace equation (43) with  

(A15)

where v₀ =σ/k is the phase velocity of the wave. Here, we have also defined x_R and x₀ by writing θ_R = kx_R−σt and θ₀ = kx₀−σt, respectively, where it will be remembered that θ₀ is the value of the angle variable θ_R at the stable fixed point of the Hamiltonian.

The distribution function for the stars trapped in the wave, i.e. for the resonant stars, is now  

(A16)

where v₀ = (kv₀/k) +v_⊥ and  

(A17)

(see equations 60 and 65).

The distribution function for the untrapped stars, i.e. for the non-resonant stars, is  

(A18)

(see equation 48).

Finally, inasmuch as the perturbation of the potential consists of a single basis function V₁(κ) exp(ik·x− iσt), equation (58) expressing the self-consistency of the wave reduces to  

(A19)

by virtue of equations (A1) and (A4).Equation (A19) is to be compared with equation (A7).

For a plane wave described by equation (A19), inequality (68) reads  

(A20)

Where the equality holds, this condition gives the gravitational counterpart of Vlasov's dispersion relation for a plasma wave. The manner in which the mode of the Vlasov type bounds the allowed, continuous spectrum of modes of the van Kampen type is illustrated in what follows.

8.3 A An example of an allowed spectrum of van Kampen modes

We define a one-dimensional velocity distribution  

(A21)

for the unperturbed system, and we let  

(A22)

where j is a constant.

For this model, the integral in inequality (A20) is readily evaluated and the inequality brought to the form  

(A23)

where y = k/(σj) and, without loss of generality, we let σ੾ 0.

For an assigned value of y, the right-hand side of inequality (A23) sets an upper bound σ_max on the magnitude of σ. Thus, σ = σ_max at a wavenumber k =σjy. We have σ = σ_max when the equality holds in conditions (A20) and (A23). In other words, σ_max satisfies the gravitational counterpart of the Vlasov (1945) dispersion relation.

Fig.A1 is a plot of σ_max against k for a plane wave in the model that we are considering. Here, we measure frequencies in the unit (6πGρ₀)^1/2 and wavenumbers in the unit (6πGρ₀j²)^1/2. The frequencies in the continuous spectrum of the van Kampen modes which are allowed by inequality (A23) fill the region under the curve in Fig.A1. Note that σ_max→ 0 as k→k_J = (12πGρ₀j²)^1/2, where k_J (=2^1/2 in the adopted system of units) is the Jeans wavenumber for the model. Evidently, the allowed modes of the van Kampen type occur in this model only at wavenumbers such that k⩽k_J and the system is gravitationally unstable.

Figure A1.

The allowed real frequencies in the continuous spectrum of modes of the van Kampen type in an infinite, homogeneous stellar system. For a plane wave, the curve represents the solution of a dispersion relation of the Vlasov (1945) type for the dependence of the frequency σ, measured in the unit (6πGρ0)^1/2, on the wavenumber k, measured in the unit (6πGρ0 j 2)^1/2. The allowed frequencies of the van Kampen modes fill the region under the curve representing the Vlasov mode.

Appendix

Appendix B: Complementarity Of The Initial-Value Formulation And The Normal-Mode Formulation Of The Matrix Method

In this appendix, we compare the representations of perturbations in the conventional, initial-value formulation and the present, normal-mode formulation of the matrix method of Kalnajs. This comparison serves to illuminate the distinction between an initial-value formulation and a normal-mode analysis in the study of small perturbations in galaxies and to demonstrate the complementarity of the two approaches. Although the present comparison is made within the framework of the matrix methods, the general results would appear to have a wider validity. It should be noted that what is derived below for stellar systems is essentially what one would infer from the literature of plasma physics. Moreover, such results are already well understood by some authors in stellar dynamics (see e.g.section 2.3 in Weinberg 1991b). On the other hand, the previous literature on stationary oscillations and modes of the van Kampen type in galaxies is sparse, and it does not appear that the results that follow have been described explicitly elsewhere.

We begin the comparison by considering modes of instability, i.e. exponentially growing modes. For such perturbations, we have Im(σ) > 0, and it follows that σ−n·ω(J) ≠ 0 for all values of n and all values J accessible to stars in the unperturbed system. Accordingly, the solution for the trigonometric components of the perturbation of the distribution function given in equation (8) reduces to  

(B1)

As we have explained in Section 1, an initial-value formulation of the matrix method leads promptly to the investigation of perturbations of the form of exponentially growing oscillations. The reduction of the governing equations follows an analysis like that in Section 2 through equation (7). However, in place of equation (8), one adopts the solution of equation (7) given by equation (B1), as explained in the preceding paragraph. With that choice, the initial-value formulation leads, with appropriate modifications of the reductions in the remainder of Section 2, to a system of matrix equations of the form of equations (25). Thus, for growing normal modes, the initial-value formulation and the normal-mode analysis are equivalent.

For stationary and exponentially damped oscillations the two approaches differ. Whereas a normal-mode analysis makes use of equation (8) in all cases, an initial-value formulation makes use of an analytic continuation of equation (B1) into regions of the complex-frequency plane in which Im(σ) ⩽ 0.

The differences are illustrated by the following question of practical interest in applications of the matrix method and related methods to the study of small perturbations in galaxies. Consider a linear sequence of models of a stellar system that are distinguished by the values of a parameter τ, say, which can be varied continuously. Let there be a mode which is unstable for models described by values of τ in some range. In general, the characteristic frequency σ of the mode will be a function of τ. Let the mode be marginally unstable when τ takes the value τ = τ_crit, say. In other words, let  

(B2)

Now, inasmuch as Im(σ) = 0, the initial-value formulation requires that the distribution function for the marginal mode be given by an analytic continuation of equation (B1). Evidently, the marginal mode would be a stationary oscillation. The question arises as to whether or not the oscillation is a normal mode of the van Kampen type. In other words, is the analytic continuation of equation (B1) a distribution function of the form given by equation (8)?

According to the reduction described, for example, in section 4 of Jackson (1960), sections 8-7 to 8-9 of Stix (1992), and section 5.4 of Palmer (1994), the analytic continuation of equation (B1) to the real-frequency axis reduces to an application of Plemelj's identity  

(B3)

Thus, the initial-value formulation gives  

(B4)

as the distribution function for the mode of marginal instability.

It remains to show that equation (B4) is a special case of equation (8). For this purpose, we write the second term on the right-hand side of equation (8) in terms of the canonical momenta G_R, G₂ and G₃ generated, with ϵ set equal to zero, in the sequence of canonical transformations represented by equations (30)313233343536(37). The required transformation of variables is J = nG_R+G_⊥, where G_⊥ = (0, G₂, G₃). We recall that the variables ϖ₁ and ϖ₂ label points on the surface in the action space on which σ−n·ω(J) = 0. As we have seen in Section 3, the canonical momenta G_R, G₂ and G₃ were constructed so that the variables G₂ and G₃ also label points on that surface. Accordingly, we let (ϖ₁, ϖ₂) = (G₂, G₃) and we define a function G₀(n, ϖ, σ) = G₀(n, G_⊥, σ) by requiring that σ−n·ω[nG₀(n, G_⊥, σ) +G_⊥]= 0. Here, and in what follows, it is to be understood that a function of ϖ is a function of the variables ϖ₁ and ϖ₂, whereas a function of G_⊥ is essentially a function of the variables G₂ and G₃. It follows from the definition of G₀(n, G_⊥, σ) that J_R(n, ϖ, σ) = nG₀(n, G_⊥, σ) +G_⊥ and that we may write  

(B5)

on the right-hand side of equation (8), where n⁻¹₁ is the Jacobian determinant of the transformation J = nG_R+G_⊥. After making that substitution, we can reduce equation (8) to  

(B6)

The delta functions on the right-hand sides of equations (B4) and (B6) are essentially the same delta function inasmuch as G_R = G₀(n, G_⊥, σ) whenever σ−n·ω(J) =σ−n·ω(nG_R+G_⊥) = 0. The precise relationship is  

(B7)

It follows that equation (B4) is a special case of equation (B6) in which  

(B8)

This verifies that the mode of marginal instability derived in the initial-value formulation of the matrix method is a normal mode of the van Kampen type.

Although modes of the van Kampen type with real frequencies are not constrained by a characteristic equation, in general, the frequencies of modes of marginal instability do satisfy a characteristic equation. For, if we now replace (G₂, G₃) with (ϖ₁, ϖ₂) in equation (B8), substitute from the resulting equation for the function λ(n, ϖ, σ) in equations (20), and make use of equations (16) in the result, then we find that equations (18) reduce to a homogeneous system of linear equations governing the coefficients μ(β, σ). A solution for the coefficients exists if and only if the secular determinant of the system vanishes. That condition is the characteristic equation governing the frequencies of the modes of marginal instability.

Thus, the initial-value formulation and the normal-mode analysis of small perturbations in galaxies are equivalent for both unstable and marginally unstable modes of oscillation. The equivalence is of practical interest for investigations in which unstable modes are delineated in order to determine criteria for instability (Weinberg 1991a; Evans & Read 1998b). The foregoing analysis shows that a mode of marginal instability is a stationary oscillation of the van Kampen type and indicates that, numerically, the marginal modes can be located (and criteria for instability determined) by solving the characteristic value problem for unstable modes with arbitrarily small growth rates.

The analytic continuation to the real-frequency axis of solutions of the governing equations in the initial-value formulation of the matrix method does not represent the continuous spectrum of normal modes of the van Kampen type. This might appear to be a paradox, because a van Kampen mode would be a solution of an initial-value problem in which the initial perturbation, at a finite time, has the structure of a van Kampen mode. The paradox is only apparent, inasmuch as the version of the matrix method that we are considering incorporates an initial condition in which the perturbation vanishes in the distant past. A van Kampen mode does not vanish in the distant past and cannot satisfy the required initial condition.

As it stands, equation (B1) also describes normal modes in which Im(σ) < 0. It was found in Section 2.3.3 that the characteristic equation governing σ is det|M(σ, α, β)|= 0, where M(σ, α, β) is the matrix defined in equation (19). It is readily verified that det|M(σ*, α, β)|= det|[M(σ, β, α)]*|= det|[M(σ, α, β)]|*= 0, inasmuch as we can swap the rows and columns of a determinant without altering its value. Thus, if σ is the characteristic frequency of a mode, then σ* is also the characteristic frequency of a mode. In other words, exponentially damped normal modes and exponentially growing normal modes occur in pairs with complex-conjugate frequencies.

In accordance with Landau's prescription, the analytic continuation of equation (B1) into the region of the complex-frequency plane in which Im(σ) < 0 gives  

(B9)

(Jackson 1960; Stix 1992; Palmer 1994). Evidently, equation (B9) does not represent a normal mode. The exponential decay of a perturbation described by equation (B9) is a manifestation of Landau damping.