Unstable Galaxy Models

The dynamics of collisionless galaxy can be described by the Vlasov-Poisson system. By the Jean's theorem, all the spherically symmetric steady galaxy models are given by a distribution of {\Phi}(E,L), where E is the particle energy and L the angular momentum. In a celebrated Doremus-Feix-Baumann Theorem, the galaxy model {\Phi}(E,L) is stable if the distribution {\Phi} is monotonically decreasing with respect to the particle energy E. On the other hand, the stability of {\Phi}(E,L) remains largely open otherwise. Based on a recent abstract instability criterion of Guo-Lin, we constuct examples of unstable galaxy models of f(E,L) and f(E) in which f fails to be monotone in E.


Introduction
A galaxy is an ensemble of billions of stars, which interact by the gravitational field they create collectively. For galaxies, the collisional relaxation time is much longer than the age of the universe ( [6]). The collisions can therefore be ignored and the galactic dynamics is well described by the Vlasov -Poisson system (collisionless Boltzmann equation) where (x, v) ∈ R 3 × R 3 , f (t, x, v) is the distribution function and U f (t, x) is its self-consistent gravitational potential. The Vlasov-Poisson system can also be used to describe the dynamics of globular clusters over their period of orbital revolutions ( [9]). One of the central questions in such galactic problems, which has attracted considerable attention in the astrophysics literature, of [5], [6], [9], [17] and the references there, is to determine dynamical stability of steady galaxy models. Stability study can be used to test a proposed configuration as a model for a real stellar system. On the other hand, instabilities of steady galaxy models can be used to explain some of the striking irregularities of galaxies, such as spiral arms as arising from the instability of an initially featureless galaxy disk ( [5]), ( [18]). In this article, we consider stability of spherical galaxies, which are the simplest elliptical galaxy models. Though most elliptical galaxies are known to be nonspherical, the study of instability and dynamical evolution of spherical galaxies could be useful to understand more complicated and practical galaxy models. By Jeans's Theorem, a steady spherical galaxy is of the form f (x, v) ≡ µ(E, L), where the particle energy and total momentum are E = 1 2 |v| 2 + U (x), L = |x × v| , and U µ (x) = U (|x|) satisfies the self-consistent nonlinear Poisson equation The isotropic models take the form f (x, v) ≡ µ(E). The case when µ ′ (E) < 0 (on the support of µ(E)) has been widely studied and these models are known to be 1 stable ( [1], [2], [7], [2], [8], [20], [14], [19], [13]). To understand such a stability, we expand the well-known Casimir-Energy functional (as a Liapunov functional) which is conserved for all time with (this is possible only if µ ′ < 0!). Upon using Taylor expansion around a steady galaxy model of [f (E, L), U ], the first order variation vanishes due to the choice of (3), and the second order variation takes the form of The remarkable feature for stability lies in the fact that H ′′ [2], [7], [14]). This crucial observation leads to the conclusion that galaxy models with monotonically decreasing energy µ ′ < 0 are linearly stable. On the other hand, for the case µ is not monotone in E, (3) breaks down, and H ′′ µ is not well-defined, which indicates possible formation of instability. It is important to note that a negative direction of the quadratic form H ′′ f [g, g] < 0 does not imply instability. In [12], an oscillatory instability was found for certain generalized 2 , m < 0, by using N-body code. This instability was later reanalyzed by more sophisticated N-body code in [3]. Despite progresses made over the years (e.g., [3], [10], [17], [18]), no explicit example of isotropic galaxy model µ(E) are known to be unstable.
The difficulty of finding instability lies in the complexity of the linearized Vlasov-Poisson system around a spatially non-homogeneous µ(E, L) : for which the construction for dispersion relation for a growing mode is mathematically very challenging. In a recent paper ( [11]), a sufficient condition was derived rigorously as follows: be a steady galaxy model. Assume that f (E, L) has a compact support in x and v, and µ ′ is bounded. Define auxiliary quadratic form for a spherically symmetric function φ(|x|) as where r 1 (E, L) and r 2 (E, L) are two distinct roots to the equation and the averageφ is defined as .
The quadratic form [A 0 φ, φ] in the above instability criterion involves delicate integrals along particle paths. The purpose of the current paper is to use numerical computations to construct two explicit examples of galaxy models for which a test function φ satisfying [A 0 φ, φ] < 0 exists. This ensures their radial instability by Guo-Lin's Theorem 1. The first example is an anisotropic model with radial instability. There are two differences to the oscillatory instability found in [12] and [3]: here the distribution function is non-singular and the instability is non-oscillatory. The second example is a non-singular isotropic model with radial instability. To our knowledge, this provides the first example of unstable isotropic model.
Even though the galaxy models studied here are not actual ones observed, our results demonstrates how to apply Guo-Lin's Theorem 1 to detect possible instability for a given galaxy model. It is our hope to foster interactions between mathematical and astronomical communities and to advance the study of instability of real galaxy models.

Examples of Unstable galaxy models
The graph of µ(E) is showed in Figure 1 below. Choosing U (0) = 3, we numerically calculate U (r), and the graph of U (r) is showed in Figure 2. By choosing the test function φ = e −r , then numerical computation below shows and hence f 0 is unstable by Theorem 1.  choose the test function with v 1 , v 2 , ...v 7 given by (16). Then [A 0 φ, φ] = −24. 016 < 0, and instability of the profile (9) follows from Theorem 1.

Numerical implementation
3.1. Numerical computation for Example 1. By the results in [21], for a steady state with f 0 (E, L) = µ (E) L 2l of Vlasov-Poisson system, the steady potential U (r) satisfies the equation where Γ denotes the gamma function. The boundary condition is lim r→∞ U (r) = 0. The computation of finding φ with [A 0 φ, φ] < 0 is carried out in the following steps.
Step (8), so l = 2 and Letting M = 2 l+7/2 π 2 c l,−1/2 , we obtain from (10) which is equivalent to We use the boundary conditions U ′ (0) = 0 and U (0) = 3. Let u 1 (r) = U (r), u 2 (r) = u ′ 1 (r), we transform the above equation to the following system (11) u ′ 1 (r) = u 2 u ′ 2 (r) = − 2 r u 2 + M r 2l g l+1/2 (u 1 ) with u 1 (r) = 3 and u 2 (0) = 0. We remark that we can subtract the finite limit of U at infinity from U and redefine E 0 accordingly. We apply Runge-Kutta method to solve equations (11), as follows First let h = 0.1 to get the values of U and U ′ at points x n = nh, then we use piecewise cubic Hermite interpolation to get an approximation of U (r).
Step 2. Computation of roots of the equation For fixed E, this equation has two solutions r 1 < r 2 when L < L max , one solution r * when L = L max and no solution when L > L max . Here, L max (E) = r * 3 U ′ (r * ) and r * (E) satisfies the equation which comes from the combination of the equations −U ′ (r) + L 2 /r 3 = 0 and (12). We employ Newton method to find the unique root r * of (13) by the Newton iteration Chose an initial point r 0 = 1, we get r * (with the stopping criterion |r n − r n+1 | < 10 −10 ). Table 1 and Figure 5 show the relation of L max to E. For L < L max , we apply Newton's iteration   Figures 8 and 9 show how r 1 and r 2 (r 1 < r 2 ) change with respect to L, when E is given. (Note that the equation (12) has a unique root r * when L = L max , thus when L approaches L max , the distance r 2 − r 1 tends to zero.) Step 3. Computing A(φ, φ) in (6). Denote the integrand Table 2 shows some of the values of I(E, L).
For given E, the integration interval of I(E, L) in (6) is [0, L max ] for L. When L = 0 and L = L max , (7) has only one solution instead of two. To avoid this problem, we integrate I(E, L) for L in the truncated interval [a 1 , L max − a 2 ], where a 1 , a 2 > 0. Choose a 2 = 10 −5 , and compare the value of Lmax−a2 a1 I(E, L)dL when a 1 = 0.01 and a 1 = 0.001 in Table 3. We see that the error is approximately 10 −7 .
In the above, we use the estimate where the effective potential due to the facts that lim and T (E, L) is mono- Numerically computing the coefficient before a 2 in the right hand side of (15), we get 0.605275913474830. Choose a 2 = 10 −5 , the right side of (15) if of order 10 −5 , which is good enough for what we need later on.
Step 5. Conclusion. We choose φ(r) = e −r . The first term of (6) is |∇φ| 2 = π, and the second term is computed to be −31.733535998660550, by choosing h = 0.1 in Runga-Kutta Method when solving U (r). Therefore (Aφ, φ) = π − 31.733535998660550 < 0. We gradually decrease h, and repeat the whole process to calculate the second term of (6). The result is shown in Table 4. Thus when h becomes smaller, the second term of (6) is less than −45. Combining with the error estimates of (14) and (15) in Step 4, we can ensure (A 0 φ, φ) < 0.  (9), we choose φ(r) to be a linear combination of several functions e −αir , 1 ≤ i ≤ n, φ(r) = a i e −αir , with a i to be determined. In the quadratic form (6), the first term is computed to be and the second term is The corresponding matrix for this quadratic form is S = (s ij ), where s ii = b i + c i , s ij = b ij + c ij . Denote λ min to be the minimum eigenvalue of the matrix S, then λ min < 0 guarantees that there exists φ such that (A 0 φ, φ) < 0. We choose In (9), there are five free parameters in µ: a, k, t 1 , E 0 , C 1 and E 1 is determined by E 1 = (E 0 + t 1 a)/(t 1 + 1). Using U 0 (the initial value of the potential) as an additional parameter, there are totally six parameters in our calculations. The idea is to view λ min as a function of these parameters λ min = λ min (a, k, t 1 , E 0 , C 1 , U 0 ). Our goal now is to find proper parameters such that λ min < 0. The numerical methods are the same as in Step 1 -Step 3 in the computation of Example 1, except that now we use self-adaptive Runge-Kutta Method to solve (11). We next verify that the following choice of parameters C 1 = 2, E 0 = 5.1, a = 1.9, k = 2.01, t 1 = 1.5, U 0 = −15.1, would lead to λ min < 0 and (A 0 φ, φ) < 0. In Table 5 we show the results of λ min under different accuracies. where a 1 and a 2 are the parameters used to truncate [0, L max ] to [a 1 , L max − a 2 ]. We also calculate the eigenvector corresponding to   Let We draw a picture of this function φ in Figure 10. Using (6), we obtain (Aφ, φ) = 35.231484866282720 − 59.247837621248109 < 0, when the integration accuracy is 10 −13 , a 1 = 10 −6 , a 2 = 10 −7 , and the root accuracy is 10 −11 .

Summary
We study the instability of spherical galaxy models in the Vlasov theory for collisionless stars. Based on the instability criterion of Theorem 1 and careful numerical computations, we have constructed two explicit unstable galaxy models f 0 (E, L) and f 0 (E) with distributions (8) and (9) respectively. In particular, f 0 (E) in Example 2 provides the first example of unstable isotropic galaxy which has not been found in literature. The instability in these examples are radial and non-oscillatory. Compared with the usual N-body codes of finding instability of galaxy models, our method only requires numerical evaluation of certain explicit integrals. Therefore, it is much more reliable and easier to implement. It is hoped that Theorem 1 can be employed to detest instability for other galaxy models in the future.