Stability and instability of self-gravitating relativistic matter distributions

We consider steady state solutions of the massive, asymptotically flat, spherically symmetric Einstein-Vlasov system, i.e., relativistic models of galaxies or globular clusters, and steady state solutions of the Einstein-Euler system, i.e., relativistic models of stars. Such steady states are embedded into one-parameter families parameterized by their central redshift $\kappa>0$. We prove their linear instability when $\kappa$ is sufficiently large, i.e., when they are strongly relativistic, and that the instability is driven by a growing mode. Our work confirms the scenario of dynamic instability proposed in the 1960s by Zel'dovich \&Podurets (for the Einstein-Vlasov system) and by Harrison, Thorne, Wakano, \&Wheeler (for the Einstein-Euler system). Our results are in sharp contrast to the corresponding non-relativistic, Newtonian setting. We carry out a careful analysis of the linearized dynamics around the above steady states and prove an exponential trichotomy result and the corresponding index theorems for the stable/unstable invariant spaces. Finally, in the case of the Einstein-Euler system we prove a rigorous version of the turning point principle which relates the stability of steady states along the one-parameter family to the winding points of the so-called mass-radius curve.


Introduction and main results
We consider a smooth spacetime manifold M equipped with a Lorentzian metric g αβ with signature (− + + +). The Einstein equations read where G αβ is the Einstein tensor induced by the metric, and T αβ is the energy-momentum tensor given by the matter content of the spacetime; Greek indices run from 0 to 3, and we choose units in which the speed of light and the gravitational constant are equal to 1.
To obtain a closed system, the field equations (1.1) have to be supplemented by evolution equations for the matter and by the definition of T αβ in terms of the matter and the metric. We consider two matter models, namely a collisionless gas as described by the collisionless Boltzmann or Vlasov equation and an ideal fluid as described by the Euler equations. This results in the Einstein-Vlasov and Einstein-Euler system respectively. We study these systems under the assumption that the spacetime is spherically symmetric and asymptotically flat, but we first formulate them in general, together with our main results.

The Einstein-Vlasov system
In the case of a collisionless gas matter is described by the number density f of the particles on phase space. The world line of a test particle on M obeys the geodesic equatioṅ x α = p α ,ṗ α = −Γ α βγ p β p γ , where x α denote general coordinates on M , p α are the corresponding canonical momenta, Γ α βγ are the Christoffel symbols induced by the metric g αβ , the dot indicates differentiation with respect to proper time along the world line of the particle, and the Einstein summation convention is applied. We assume that all the particles have the same rest mass m 0 ≥ 0 and move forward in time, i.e., their number density f is a non-negative function supported on the mass shell P M := g αβ p α p β = −m 2 0 , p α future pointing , a submanifold of the tangent bundle T M of the spacetime manifold M which is invariant under the geodesic flow. Since we are interested in the massive case m 0 > 0, we normalize m 0 = 1. But at a crucial point in our analysis the massless case m 0 = 0 will play an important role. Letting Latin indices range from 1 to 3 we use coordinates (t, x a ) with zero shift which implies that g 0a = 0. On the mass shell P M the variable p 0 then becomes a function of the remaining variables (t, x a , p b ): p 0 = −g 00 1 + g ab p a p b .
Since the particles move like test particles in the given metric, their number density f = f (t, x a , p b ) is constant along the geodesics and hence satisfies the Vlasov equation The energy-momentum tensor is given by T αβ = p α p β f |g| 1/2 dp 1 dp 2 dp 3 −p 0 , (1.3) where |g| denotes the modulus of the determinant of the metric, and indices are raised and lowered using the metric, i.e., p α = g αβ p β . The system (1.1), (1.2), (1.3) is the Einstein-Vlasov system in general coordinates. We want to describe isolated systems, and therefore we require that the spacetime is asymptotically flat. In Section 3 we will see that steady states of this system can be obtained via an ansatz where E = −p 0 is the local or particle energy and φ is a prescribed ansatz function; we refer to (1.4) as the microscopic equation of state. We can now state our main result for the Einstein-Vlasov system in an informal way; the precise meaning of the parameter κ > 0 is explained in Section 3, in particular (3.10) and (3.11).
Theorem 1. Let (f κ ) κ>0 be a one-parameter family of spherically symmetric steady states to the Einstein-Vlasov system with a fixed, decreasing microscopic equation of state φ and κ the central redshift of the corresponding steady state.
(a) For κ sufficiently large, the associated steady state is dynamically unstable in the sense that the linearized Einstein-Vlasov system possesses an exponentially growing mode.
(b) For any κ > 0 the phase space of the linearized system splits into three invariant subspaces: the stable, unstable, and center space. The dimension of the stable/unstable space is equal to the negative Morse index of a suitably defined, macroscopic Schrödinger-type operator.
The Einstein-Vlasov system models large stellar systems such as galaxies or globular clusters, and the central redshift is a measure of how relativistic the corresponding state is. Part (a) of the theorem says that strongly relativistic galaxies or globular clusters are unstable. This is in sharp contrast to the corresponding Newtonian situation, i.e., to the Vlasov-Poisson system. For this system any steady state induced by a strictly decreasing microscopic equation of state is non-linearly stable, cf. [28,29,46,54]. The instability result in Theorem 1 is therefore a genuinely relativistic phenomenon.
A proof of instability which does not rely on the refined spectral information about the existence of growing modes can be found in Theorem 4.8, while the existence of a growing mode is given by Theorem 4.25. A rigorous statement of exponential trichotomy is provided in Theorem 4.28.

The Einstein-Euler system
In this case the matter is described by its mass-energy density ρ, its four velocity u α which is a future-pointing, time-like unit vector field, and its pressure p. These quantities are defined on the spacetime manifold M and induce the energy-momentum tensor T αβ = (ρ + p)u α u β + p g αβ . (1.5) The Bianchi identity applied to the field equations yields the evolution equations for the fluid where ∇ α denotes the covariant derivative associated with the metric. More explicitly, (ρ + p)u α ∇ α u β + (g αβ + u α u β )∇ α p = 0. (1.7) In order to close the system we need to prescribe a (macroscopic) equation of state which relates p and ρ, p = P (ρ), (1.8) with a prescribed function P . Notice that p will always denote the pressure as a function on spacetime, while P will denote its functional relation to the mass-energy density ρ. We state our main result for the Einstein-Euler system, which consists of (1.1), (1.5), (1.6), (1.7), (1.8), in an informal way.
Theorem 2. Let (ρ κ , p κ , u α κ ≡ 0) κ>0 be a one-parameter family of spherically symmetric steady states to the Einstein-Euler system with a fixed, strictly increasing equation of state (satisfying suitable assumptions) and κ the central redshift of the corresponding steady state.
(a) For κ sufficiently large, the associated steady state is dynamically unstable in the sense that the linearized Einstein-Euler system possesses an exponentially growing mode.
(b) For any κ > 0, the phase space of the linearized system splits into three invariant subspaces: the stable, unstable, and center space. The dimension of the stable/unstable space is equal to the negative Morse index of a suitably defined Schrödinger-type operator Σ κ .
(c) A version of the turning point principle holds, i.e., as κ ∈]0, ∞[ varies, the stability of the steady state can be inferred from the so-called mass-radius diagram and the knowledge of the negative Morse index of Σ κ .
The Einstein-Euler system models stars, and hence part (a) of the theorem says that strongly relativistic steady stars-those with very large central redshift-are unstable. The analogous comment as in the Vlasov case applies concerning the Newtonian situation. The rigorous statement and proof of the turning point principle is given in Theorem 5.17, the exponential trichotomy is shown in Theorem 5.20, and the instability for large κ can be found in Corollary 5.25.

Motivation and history of the problem
The Newtonian analogue of the Einstein-Vlasov system is the gravitational Vlasov-Poisson system Here f = f (t, x, v) ≥ 0, a function of time t, position x ∈ R 3 , and velocity v ∈ R 3 , is the phase-space density of the stars in a galaxy, and U = U (t, x) is the gravitational potential induced by the macroscopic, spatial mass density ρ = ρ(t, x); integrals without explicitly specified domain always extend over R 3 . A convenient way of constructing steady states of this system is to make an ansatz is the local particle energy. Modulo some technical assumptions the profile Φ ≥ 0, which we refer to as the microscopic equation of state, is arbitrary, but must vanish on ] − ∞, 0[. Hence E 0 < 0 represents the maximal particle energy allowed in the system. Any such choice of f satisfies (1.9) with the given potential U , and the problem of finding a steady state is reduced to solving (1.10), where the right hand side becomes a function of U obtained by substituting the ansatz into (1.11). Since such steady states necessarily are spherically symmetric, it seems convenient to parameterize them by U (0). But in view of the boundary condition in (1.10) and the need of a cut-off energy E 0 in the ansatz this seems one free parameter too many. Instead, one can use y = E 0 − U as the basic unknown function and use κ = y(0) > 0 as the free parameter. Under suitable assumptions on Φ the function y has a unique zero at some radius R > 0 which marks the boundary of the support of the induced steady state which then also has finite mass, cf. [56] and the references there. With E 0 = lim r→∞ y(r) < 0 the corresponding potential U vanishes at infinity, and the parameter κ = y(0) = U (R) − U (0) is the potential energy difference between the center and the boundary of the compactly supported steady state. Examples for ansatz functions Φ which yield such one-parameter families of steady states are the polytropes where Φ(η) = η k for η > 0 with −1 < k < 7/2, and the King model where Φ(η) = e η , both of which are used in astrophysics, cf. [9] . A central challenge in the qualitative description of the Vlasov-Poisson dynamics is to understand the behavior of solutions close to the above steady states. In the astrophysics literature formal arguments towards linearized stability were given for example by Antonov [7], Doremus, Feix, Baumann [18], and Kandrup, Sygnet [41]. The stability criterion proposed by these authors is that Φ is strictly increasing on [0, ∞[. This is physically reasonable, as it states that the number of stars in the galaxy is a decreasing function of their energy. Much work has been invested in a mathematically rigorous, nonlinear proof of this stability criterion, cf. [28,29,30,45,46] and the references there. Remarkably, the size of κ is irrelevant for these stability results.
Theorem 1 shows that the above stability criterion for Vlasov-Poisson is false for the Einstein-Vlasov system: When the central redshift κ is sufficiently large, there exists a growing mode despite the monotonicity assumption on Φ. This fundamental difference between the Vlasov-Poisson and the Einstein-Vlasov system is driven by relativistic effects which become dominant at sufficiently large values of the central redshift κ. Another difference between these two systems seems to be that in the non-relativistic case nonlinearly stable steady states can be obtained as minimizers of suitable energy-Casimir functionals, cf. [60] and the references there, but this does not work so well in the relativistic case, as the problem is supercritical; in [72] an attempt in this direction has been made, but as shown in [5] that paper is not correct. We also mention an interesting recent work of Andersson and Korzyński [1] on the variational derivation of the Einstein-Vlasov system.
The Newtonian analogue of the Einstein-Euler system is the Euler-Poisson system, where the Euler equations ∂ t ρ + div(ρ u) = 0, ρ (∂ t + u · ∇) u + ∇p + ρ∇U = 0 are coupled to the Poisson equation (1.10). Here ρ is the macroscopic fluid density, and u is the fluid 3-velocity. To close this system, one must prescribe an equation of state p = P (ρ). A well-investigated choice is again the polytropic one, i.e., P (ρ) = ρ γ , γ > 1. When γ > 6 5 there exist spherically symmetric stationary distributions of compact support, cf. [9,56]. A classical, linear stability analysis shows that polytropic steady stars are stable if γ > 4 3 . By analogy to the Vlasov-Poisson case this result holds independently of the size of the central density ρ(0) of the steady state. Theorem 2 shows that this changes drastically in the relativistic context: Strongly relativistic stars, i.e., stars with large central density or equivalently large central redshift, are unstable.
The subject of stability of isolated self-gravitating solutions of the Einstein-Euler system has a long history. The first linear stability result and a variational characterization of the stability question in spherical symmetry goes back to the pioneering work of Chandrasekhar [16] in 1964. The same stability criterion was then derived by Harrison, Thorne, Wakano, and Wheeler in 1965 [34], who were the first to show that Chandrasekhar's result is equivalent to the positive-definiteness of the second variation of the ADM mass, subject to the constant baryon number constraint, see also [8]. Moreover, a turning point principle is formulated in [34] and investigated numerically. It proposes that along the one-parameter family stability changes to instability or vice versa when the steady state passes through an extremal point of the so-called mass-radius curve. Recently, a proof of such a turning point principle was given in [31]. For a detailed survey of these and related results from the 1960s see [71]. A variational approach to stability and the turning point principle are elegantly formulated in [13,14,15]. In Theorem 5.17 we formulate and prove a rigorous version of the turning point principle for the Einstein-Euler system, i.e., of part (c) of Theorem 2. A general abstract discussion of turning point principle is astrophysics is given in [68]. For the different notion of thermodynamic instabilities and its relation to turning point principles, cf. the discussions in [43,66]. We refer the reader to [27,69] for a comprehensive overview of this topic from the physics point of view. By contrast to what is known in the vicinity of massive steady states, the vacuum solution of the Einstein-Vlasov system is asymptotically non-linearly stable against small perturbations [23,50,61].
A burst of interest in the stability of highly relativistic, self-gravitating bodies occurred in the mid 1960s with the discovery of quasars [10,75]. Since these objects are very powerful and very concentrated sources of energy, it was initially unclear whether their high redshifts were due to them being far away or due to them being very relativistic. Zel'dovich and Podurets gave an early contribution to the subject by showing numerically that certain spherically symmetric steady states turn unstable when their central redshift passes a threshold value [76]. This lead Ipser, Thorne, Fackerell and others to initiate a systematic study of the question of stability of steady relativistic galaxies [19,20,21,22,36,37,38,39]. In particular, a turning point principle for the Einstein-Vlasov system was formulated. It states that the transition from a stable to an unstable configuration occurs at a critical point of the binding energy, plotted as a function of the central redshift. For some numerical results in this direction see [57,67]. On the other hand, Bisnovatyi-Kogan and Zel'dovich pursued the question whether there exist stable self-gravitating configurations with arbitrarily high central redshifts [10,11]. The steady states studied in [10] are singular and have infinite mass and extent. However, in Section 3 we shall see that certain explicit steady states of the massless Einstein-Vlasov system have the same macroscopic properties as the solutions in [10]. These massless solutions, which we refer to as the BKZ solutions, capture the behavior of massive steady states close to the center of the galaxy when the central redshift is sufficiently large. An analogous assertion holds for the Einstein-Euler system, see Proposition 3.13. This observation appears to be new and of independent interest, but it also plays a fundamental role in our proof of instability in part (a) of Theorems 1 and 2.
In the context of relativistic stars, Meltzer and Thorne [53] discuss the stability properties of very dense relativistic stars. They predict that perturbations localized to the core of a very dense star will result in gravitational collapse, which is consistent with our methodology, see Sections 4.2 and 5.6.

Methodology
Under natural assumptions on a given equation of state f = φ(E) or p = P (ρ), both the Einstein-Vlasov and the Einstein-Euler systems posses a corresponding one-parameter family of compactly supported steady states with finite mass, parameterized by the value of the central density ρ(0), or equivalently, by the central redshift κ > 0, cf. Section 3, in particular (3.10) and (3.11). Even though the equations of equilibrium are classical [55], a rigorous proof of the existence and finiteness of the total mass and support of a steady star/galaxy is nontrivial and depends crucially on the assumptions on the micro-and macroscopic equations of state, see [56] and references therein. In what follows M κ denotes the ADM mass and R κ the radius of the support of the corresponding steady state.
In both cases, the formal linearization around the steady state leads to a Hamiltonian partial differential equation, which comes with a rich geometric structure. In the case of the Einstein-Vlasov system the steady state can be interpreted as a critical point of the socalled energy-Casimir functional [32,40]. To prove the instability for large values of κ > 0 it is therefore natural to investigate the sign of the second variation of this energy-Casimir. In Theorem 4.3 we construct an explicit test function that turns the second variation negative for κ > 0 sufficiently large. The key to the construction is a precise understanding of the limiting behavior of the sequence of steady states (f κ ) κ>0 in the singular limit κ → ∞. We show that in a suitably rescaled annulus around the center r = 0 the behavior of the steady states is at the leading order described by the BKZ solutions mentioned above, see Sections 3.3 and 3.4. In the physics literature such a limit is referred to as the ultrarelativistic limit. Since the BKZ solutions are known explicitly, we obtain sharp a priori control over the second variation in the large κ regime, and we can make a judicious choice of a test function whose support is spatially localized to the aforementioned annulus, close to the center of the galaxy. In Theorem 4.3 we carry out careful estimates showing that such a test function makes the second variation functional negative. A similar procedure applies to the Einstein-Euler system, see Theorem 5.24.
The existence of a negative energy direction is sufficient for a proof of linear exponential instability as shown in Theorem 4.8 for the Einstein-Vlasov system; a similar proof can be given for the Einstein-Euler system. To this end we adapt a strategy developed in the study of plasma instabilities by Laval, Mercier, and Pellat [44]. However, with nonlinear applications in mind it is important to show that the instability is driven by a growing mode, a statement not afforded by Theorem 4.8. In the context of the Vlasov theory, this is a highly nontrivial question, as the presence of the continuous spectrum is generally unavoidable, and a more refined analysis is necessary. Thus to get a more precise spectral information, we must carefully analyze the linearized operator. For the Euler-Einstein system Makino showed that the spectrum of the linearized operator is in fact discrete [51,52] by formulating it as a singular Sturm-Liouville type operator on a bounded domain.
Both the linearized Einstein-Euler and the linearized Einstein-Vlasov system can be written in the second order form where L κ : X ⊃ D(L κ ) → X a self-adjoint linear operator, and X a Hilbert space. In both cases, it is shown that the second variation of the energy corresponds to to the quadratic form (L κ f, f ) X . In the case of the Einstein-Vlasov system, a naive attempt to minimize the functional f → (L κ f, f ) X leads to difficulties due to the possible loss of compactness along minimizing sequences. A related loss of compactness has been well-known in the stability theory for various Vlasovmatter systems, most notably the Vlasov-Poisson and the Vlasov-Maxwell system [47]. We adopt a different approach.
The reduced operator. The key step in our analysis is to construct a suitable reduced operator which by definition is a self-dual, macroscopic, Schrödinger-type operator 12) and for any ψ ∈Ḣ 1 (R 3 ) there exists an f ∈ D(L κ ) such that where ∆ κ is an explicit, nondegenerate, linear, elliptic operator with variable coefficients, V κ a compactly supported potential, andḢ 1 (R 3 ) is a homogeneous Sobolev space. Since S κ has at most finitely many negative eigenvalues, we can use the bounds (1.12), (1.13) to conclude that the negative Morse index of L κ is finite. The derivation of S κ in both the kinetic and the fluid case is new. An attempt to construct a reduced operator was made earlier by Ipser [37], but in that work the bound (1.12) with a different choice of S κ is satisfied only under an additional hypothesis on the steady states, which appears to be hard to verify rigorously. We point out that a related construction of the reduced operators for the Newtonian analogues is nontrivial, yet considerably simpler; for the Vlasov-Poisson system see [30] and for the Euler-Poisson system [49].
While the reduced operator can be constructed for any κ > 0, see Theorems 4.24 and 5.15, only in the regime of large κ we know that there exists at least one negative direction. By proving that L κ is self-adjoint and not merely symmetric, we then infer that the unstable space consists of eigenfunctions of finite multiplicity. The proof of selfadjointness of the operator L κ in both the Einstein-Vlasov and the Einstein-Euler case is not obvious and is presented in Lemmas 4.21 and 5.23 respectively.
Exponential trichotomy. To obtain further information about the linearized dynamics, we are forced to work with the first order Hamiltonian formulation of the problem. This part of our analysis applies to all κ > 0. Abstractly, both problems can be recast in the form ∂ t u = J κ L κ u, (1.14) where L κ : X → X ′ is a self-dual operator and J κ : X ′ → X is anti-self-dual; X is an appropriate Hilbert space. The operators J κ , L κ for the Einstein-Vlasov system are given in Lemma 4.20, for the Einstein-Euler system in Lemma 5.6. Using crucially the existence of the reduced operator and (1.14) we are able to apply the general framework developed by Lin and Zeng [48] to obtain a quantitative exponential trichotomy result, see Theorems 4.28 and 5.20. Roughly speaking, we show that the phase space around the steady state naturally splits into a stable, unstable, and center invariant subspace. Under a non-degeneracy assumption on the reduced operator S κ , we can go a step further and prove a quantitative Lyapunov stability statement on the center space in the natural energy topology, see parts (v) of Theorems 4.28 and 5.20. The exponential trichotomy statement shall provide a foundation for the understanding of the nonlinear dynamics in the vicinity of steady states. For instance, one expects that in the presence of growing modes the phase space will split into perturbations leading to gravitational collapse, and those leading to global existence via dispersion, separated by an invariant (co-dimension 1 center stable) manifold. This is consistent with the dynamical picture proposed in the study of criticality phenomena [6,25].
Turning point principle. In the context of self-gravitating relativistic bodies, it is desirable to have a simple criterion for stability that depends on certain bulk properties of the system under study. Precisely such an idea was put forward by Zel'dovich [74] and Harrison, Thorne, Wakano, and Wheeler in [34], wherein the so-called turning point principle was formulated. If one plots the mass M κ and the steady star radius R κ along a curve parameterized by the central redshift κ, it is proposed in [34] that the stability can be lost to instability and vice versa only at the extremal points of the mass plotted as a function of κ , see also [66,69].
The steady states of the Einstein-Euler system can be interpreted as critical points of the ADM mass among all densities of the same total baryon mass. This observation goes back to [34]. By computing the second variation of the ADM mass we can derive the well-known Chandrasekhar stability criterion [16], which states that the static star is (spectrally) stable if the linearized operator is non-negative on a codimension-1 subspace of a certain Hilbert space, see Proposition 5.11. Using this characterization, we show in Theorem 5.17 that when d dκ M κ d dκ

Mκ Rκ
= 0 the number of growing modes associated with the linearized operator J κ L κ equals n − (Σ κ ) − i κ , where Σ κ is the reduced operator associated with the Einstein-Euler system, n − (Σ κ ) is its negative Morse index, and the index i κ ∈ {0, 1} depends on the mass-radius curve through the formula Finally, in line with Antonov's First Law for Newtonian self-gravitating systems, we prove a relativistic "micro-macro" stability principle in Theorem 5.26, which in a precise way states that a steady galaxy with a certain microscopic equation of state is spectrally stable if an individual star with the corresponding macroscopic equation of state is spectrally stable. Some open questions. A natural question for further study is the nonlinear instability of steady states when κ ≫ 1. For the Einstein-Vlasov case a local well-posedness result was established in [61] which can be used (or if necessary adapted) for initial data close to a steady state. A non-trivial problem which needs to be understood before attacking the non-linear regime is the regularity of eigenfunctions corresponding to the growing modes of the linearized system. For the Einstein-Euler case we are naturally led to the vacuum free boundary problem wherein the location and the regularity of the boundary separating the star from the vacuum is an unknown. This question comes with a number of mathematical difficulties, and even the local-in-time well-posedness theory is an open problem in this context.
When κ ≪ 1 it is known [32,33] that the steady states of the Einstein-Vlasov system are linearly stable. Since the Einstein equations are energy supercritical, this type of a priori control coming from the linearized problem is far from sufficient for proving any kind of nonlinear stability; note that the situation is different for the Vlasov-Poisson system. It is an open problem to show that nonlinear orbital stability is true for κ ≪ 1. A linear stability result for so-called hard stars is given in [24].
A further open question is to show, as conjectured in the physics literature [34,69], that for the Einstein-Euler system the number of unstable directions grows to infinity as κ → ∞. As pointed out in Remark 5.19 there is strong evidence for this to be true.
2 The Einstein-Vlasov and Einstein-Euler system in spherical symmetry For the above systems, questions like the stability or instability of steady states are at present out of reach of a rigorous mathematical treatment, unless simplifying symmetry assumptions are made. We assume spherical symmetry, use Schwarzschild coordinates (t, r, θ, ϕ), and write the metric in the form ds 2 = −e 2µ(t,r) dt 2 + e 2λ(t,r) dr 2 + r 2 (dθ 2 + sin 2 θ dϕ 2 ). (2.1) Here t ∈ R is a time coordinate, and the polar angles θ ∈ [0, π] and ϕ ∈ [0, 2π] coordinatize the surfaces of constant t and r > 0. The latter are the orbits of SO(3), which acts isometrically on this spacetime, and 4πr 2 is the area of these surfaces. The boundary condition lim It is convenient to introduce the corresponding Cartesian coordinates x = (x 1 , x 2 , x 3 ) = r(sin θ cos ϕ, sin θ sin ϕ, cos θ).
Instead of the corresponding canonical momenta p = (p 1 , p 2 , p 3 ) we use the non-canonical momentum variables v a = p a + (e λ − 1) x · p r x a r , a = 1, 2, 3.
In these variables, and f is spherically symmetric iff The spherically symmetric, asymptotically flat Einstein-Vlasov system takes the following form: Here˙and ′ denote the derivatives with respect to t and r respectively. For a detailed derivation of these equations we refer to [59]. It should be noted that in this formulation no raising and lowering of indices using the metric appears anywhere. It is a completely explicit system of partial differential equations where x, v ∈ R 3 , x · v denotes the Euclidean scalar product, |v| 2 = v · v, and v is defined in (2.4); integrals without explicitly specified domain extend over R 3 . We note that p and p T are the pressure in the radial and tangential directions respectively, which for the Vlasov matter model in general are not equal. They are equal for the isotropic steady states which we consider in the next section, and of course also for the Euler matter model.
3 Steady states

The basic set-up; one-parameter families
We first consider the Vlasov case. The characteristic system of the stationary, spherically symmetric Vlasov equation readṡ x r , and the particle energy is constant along characteristics. Hence the static Vlasov equation is satisfied if f is taken to be a function of the particle energy, and the following form of this ansatz is convenient: Here E 0 > 0 is a prescribed cut-off energy-notice that the particle energy (3.1) is always positive. Because of spherical symmetry, the quantity is conserved along characteristics as well so that we could include a dependence on L in the ansatz (3.2). This leads to anisotropic steady states which cannot be treated in parallel with the Euler case and are not pursued in this paper. We require that Φ has the following properties; for the stability analysis in Section 4 this will be strengthened by (Φ2): and there exists −1/2 < k < 3/2 and constants c 1 , c 2 > 0 such that for all sufficiently small η > 0, These assumptions are sufficient for the results of the present section, in particular for the analysis of the steady states in the limit of large central redshift, which we believe has some interest in itself. A typical guiding example which satisfies (Φ1) is the function Φ(η) = c η k for η ≥ 0, Φ(η) = 0 for η < 0 with constants c > 0 and −1/2 < k < 3/2. These are the so-called polytropes, named so by analogy to the well-known polytropic equations of state in compressible fluid dynamics.
Since only the metric quantity µ enters into the definition of the particle energy E in (3.1), the field equations can be reduced to a single equation for µ. It is tempting to prescribe µ(0), but since the ansatz (3.2) contains the cut-off energy E 0 as another, in principle free parameter and since µ must vanish at infinity due to (2.2) this approach is not feasible. Instead we define y := log E 0 − µ so that e µ = E 0 e −y . For the ansatz (3.2) the spatial mass density and pressure become functions of y, i.e., ρ(r) = g(y(r)), p(r) = h(y(r)) = p T (r), (3.4) where g(y) := 4πe 4y The functions g and h are continuously differentiable on R, cf. [62, Lemma 2.2], and they vanish for y < 0. The metric coefficient λ can be eliminated from the system, because the field equation (2.6) together with the boundary condition (2.3) at zero imply that where the mass function m is defined by Hence the static Einstein-Vlasov system is reduced to the equation which is equivalent to (2.7); here m, ρ, and p are given in terms of y by (3.4) and (3.8).
In [56] it is shown that for every central value there exists a unique smooth solution y = y κ to (3.9), which is defined on [0, ∞[ and which has a unique zero at some radius R > 0. In view of (3.4)-(3.6) this implies that the induced quantities ρ and p are supported on the interval [0, R], and a non-trivial steady state of the Einstein-Vlasov system with compact support and finite mass is obtained. We observe that the limit y(∞) := lim r→∞ y(r) < 0 exists, since r → y(r) is a decreasing function and has a unique zero as mentioned above. The metric quantity µ is given by e µ(r) = E 0 e −y(r) , and in order that µ has the correct boundary value at infinity we must define E 0 := e y(∞) . Since y(R) = 0 we also see that E 0 = e µ(R) . We want to relate the parameter κ to the redshift factor z of a photon which is emitted at the center r = 0 and received at the boundary R of the steady state; this is not the standard definition of the central redshift where the photon is received at infinity, but it is a more suitable parameter here: Hence κ is in one-to-one correspondence with the central redshift factor z with κ → ∞ iff z → ∞, and although this is not the standard terminology we refer to κ as the central redshift. For a fixed ansatz function Φ we therefore obtain a family (y κ ) κ>0 of solutions to (3.9), which induce steady states to the Einstein-Vlasov system parameterized by the central redshift κ, and each member of this family represents a galaxy in equilibrium, which has finite mass and compact support.
Remark 3.1. The central redshift κ and the central density ρ c = ρ(0) are related via ρ c = g(κ), see (3.5). It can also be seen from (3.5) that g is monotonically increasing so that κ and ρ c are in a 1-1 relationship. In the literature both κ and ρ c are used to parameterize the steady state solutions, and the two parametrizations are equivalent.
We now consider the Euler case. For stationary, spherically symmetric solutions of the Einstein-Euler system the velocity field necessarily vanishes, u = 0, cf. (2.16 We define Then (3.12) can be written as which means that on the support of the matter, Q(ρ(r)) + µ(r) = const.
In analogy to the Vlasov case we introduce y = const − µ and find that ρ is given in terms of y, Taking into account the equation of state (1.8) it follows that p = h(y) := P (g(y)). (3.14) Hence the stationary system takes exactly the same form as in the Vlasov case and can be reduced to the equation (3.9), the only difference being the definitions of the functions g and h in (3.5), (3.6) or (3.13), (3.14) respectively. The interpretation of κ = y(0) remains as explained for the Vlasov case above.
We now specify the assumptions on the function P which defines the equation of state.
Assumption (P1). P ∈ C 1 ([0, ∞[) with P ′ (ρ) > 0 for ρ > 0, P (0) = 0, there exist constants c 1 , c 2 > 0 and 0 < n < 3 such that P ′ (ρ) ≤ c 1 ρ 1/n for all sufficiently small ρ > 0 (3.15) and the inverse of P exists on [0, ∞[ and satisfies the estimate Before we discuss examples for such equations of state in the next subsection, we briefly mention that the condition (3.15) forces the equation of state to take the approximate form P (ρ) ∼ ρ→0 ρ 1+ 1 n in the vicinity of the vacuum boundary, while (3.16) makes the equation of state linear at the leading order for large values of ρ, i.e. P (ρ) ∼ ρ→∞ 1 3 ρ. The first assumption is the well-known polytropic law from the classical gas dynamics, while the second assumption is necessary to ascertain that the speed of sound remains smaller than the speed of light at high densities.
Assumption (3.15) together with the required regularity guarantees that the function g and h defined in (3.13) and (3.14) are C 1 on R, and for each κ = y(0) > 0 the equation (3.9) has a unique solution with the same properties as stated for the Vlasov case, cf. [56]. For easier reference we collect the findings of this subsection.
(b) Let P satisfy (P1). Then there exists a one-parameter family of steady states (ρ κ , λ κ , µ κ ) κ>0 of the spherically symmetric, asymptotically flat Einstein-Euler system. In both cases these steady states are compactly supported on some interval [0, The following two identities hold on [0, ∞[, the second being known as the Tolman-Oppenheimer-Volkov equation: Proof. Eqn. (3.17), which holds also in the time-dependent case, follows by adding the field equations (2.6), (2.7) or (2.14), (2.15) respectively. For the Einstein-Euler case (3.18) was already stated above as (3.12). That this equation also holds in the Einstein-Vlasov case is due to the fact that a steady state of this system is macroscopically also one of the Einstein-Euler system, as will be seen in the next section.

Microscopic and macroscopic equations of state
Let us consider a microscopic equation of state Φ which satisfies the assumption (Φ1) (see (3.3)), and a corresponding steady state of the Einstein-Vlasov system. Given the relations (3.4) it is tempting to write and interpret this steady state as a solution to the Einstein-Euler system with the macroscopic equation of state given by we fix some κ > 0 and drop the corresponding dependence since it plays no role here. Indeed, this maneuver is perfectly rigorous and provides a class of examples P = P Φ which satisfy (P1) (see (3.15)-(3.16)). To see this, we first observe that The functions h and j and hence also g are continuously differentiable with and j ′ (y) = 2j(y) + 4π 3 cf. [62,Lemma 2.2]. In particular, the functions g and h are strictly increasing on [0, ∞[ with lim y→∞ g(y) = lim y→∞ h(y) = ∞ so that these functions are one-to-one and onto on [0, ∞[.
Proof. Let P be given by (3.19). Then for ρ > 0 and with y = g −1 (ρ) > 0, We claim that for some constant C > 0; such constants may depend only on Φ and may change their value from line to line. For 0 ≤ ρ ≤ 1 this is obvious, since by (3.25), which needs to be shown for y ≥ g −1 (1) =: y 0 . Now Combining both estimates yields (3.27) for y ≥ y 0 . We rewrite (3.26) in terms of P −1 to find that By (3.26), P (ρ) ≥ ρ/6 for ρ large which means that P −1 (p) ≤ 6p for p large so that (3.28) implies (3.16) in that case. We prove that there exist constants c 1 , c 2 > 0 such that with γ = 1/(k + 3/2), This implies (3.15) with n = k + 3/2 ∈]1, 3[, and P −1 (p) ≤ Cp 1/(γ+1) for p small. Thus for p small which completes the proof of (3.16) so that (P1) holds; notice that γ < 1 since k > −1/2. It therefore remains to prove (3.29), and it suffices to prove the estimate for P ′ which implies the one for P . In view of (3.24) and with ρ = g(y) the estimate for P ′ is equivalent to and j ′ is given by (3.23). Hence all the relevant terms are of the form e ay with a ∈ {0, 2, 4} and m ∈ {−1/2, 1/2, 3/2}, and these terms need to be estimated from above and from below for y small. We can drop the factor e ay . For Φ(η) we use the assumption (3.3), and we observe that Hence the relevant terms can be estimated from above and below by This implies that h ′ (y) can be estimated from above and below by (1 − e −y ) k+3/2 , and g(y) γ g ′ (y) can be estimated from above and below by (1 − e −y ) γ(k+3/2)+(k+1/2) . But by the choice of γ the two exponents are equal and the estimate (3.30) holds. It remains to show the regularity assertion. So let in addition Φ ∈ C 1 ([0, ∞[). We observe that We differentiate these expressions under the integral and change variables to find that Arguing as in [62, Lemma 2.2] it follows that g, h ∈ C 2 (R), and hence by (3.24), Remark 3.4. (a) The above proof shows that any equation of state derived from a microscopic ansatz satisfying (Φ1) obeys the asymptotic relation which for technical reasons is expressed in terms of P −1 in (P1). The limiting equation of state P (ρ) = ρ/3 is known in astrophysics and cosmology as the equation of state for radiation. It is remarkable that this physically reasonable behavior is taken care of automatically if the equation of state derives from a microscopic one.
(b) For the fluid model this asymptotic behavior must be put in by hand. In particular, it excludes equations of state of the form P (ρ) = cρ γ with γ > 1. Such equations of state, which also violate the requirement that the speed of sound √ P ′ must be less than the speed of light, are common for the Euler-Poisson system.
The above discussion indicates that macroscopic quantities induced by an isotropic steady state of the Einstein-Vlasov system represent a steady state of the Einstein-Euler system with an equation of state given by (3.19). To complete this argument, it remains to check the stationary Euler equation (3.18). But using the relations (3.20) and (3.22), As a genuinely fluid dynamical equation of state we consider the one used for neutron stars in [35]. Here for y ≥ 0, and the relation between p and ρ is given by or equivalently, by the equation of state Proposition 3.5. The equation of state (3.32) for a neutron star is well defined, satisfies (P1), and in addition Also, both functions are one-to-one and onto on the This implies that for ρ > 0, and For y > 0 large it holds thatg(y) ≤ Cy 4 and henceg −1 (ρ) ≥ Cρ 1/4 and for ρ > 0 large, but due to (3.34) this estimate also holds for ρ > 0 small. Hence Next we observe that for y > 0 small it holds that y 3 ≤g(y) ≤ Cy 3 and hence for ρ > 0 small, cρ 1/3 ≤g −1 (ρ) ≤ ρ 1/3 which together with (3.33) implies that for ρ ≥ 0 small; constants like 0 < c < C may change their value from line to line. In particular, (3.15) in (P1) holds. The estimate (3.35) implies that for ρ > 0, and the proof is complete.

The Bisnovatyi-Kogan-Zel'dovich (BKZ) solution
A key ingredient in our instability analysis is rather precise information on the form of the steady states for large central redshift κ. In order to obtain this information it is useful to understand that the steady states approach steady states of a certain version of the Einstein-Vlasov or Einstein-Euler system as κ → ∞. In this subsection we introduce this limiting system and discuss a special, explicit solution to it. We first consider the Vlasov case. For κ → ∞ we expect that close to the center the corresponding solution of (3.9), (3.10) is large, so that formally, and for the sake of notational simplicity we normalized Since g * and h * are strictly positive, a corresponding steady state is never compactly supported. Let us for the moment assume that we have a solution y of (3.9), (3.10) where g and h are replaced by g * and h * , and let µ, λ, ρ, p be induced by y. Then these quantities satisfy the stationary Einstein equations together with the radiative equation of state We want to understand what happens with the Vlasov equation in the limit κ → ∞. First, we note that the equation of state (3.38) can not come from an isotropic steady state particle distribution of the form (3.2) (and of course not from an anisotropic one either): The physical meaning of this is simply that massive particles are not radiation. However, this observation carries its own cure: In the κ → ∞ limit the stationary Vlasov equation (2.5) must be replaced by its ultrarelativistic, or massless version An isotropic solution of the latter equation is obtained by the ansatz it is straightforward that this satisfies (3.39). Moreover, for such a particle distribution f , as required by (3.38). Finally, which is as expected from (3.37) so that the above f is a consistent solution to the massless Vlasov equation ( Proof. The proof follows using the arguments in [58,Theorem 3.4]. We want to find a special, explicit solution of (3.40). To this end, we make the ansatz Then Substituting this into (3.40) it turns out that the latter equation holds iff β = 1 2 and 56 3 πe 4γ = 1. Thus These macroscopic data of the solution are the same as for the massive solutions found by Bisnovatyi-Kogan and Zel'dovich in [10], and we refer to it as the BKZ solution.
We note that the corresponding y is singular at the origin, and the solution violates the condition (2.3) for a regular center. It has infinite mass and also violates the boundary condition (2.2) at infinity, i.e., it does not represent an isolated system.
Remark 3.7 (Geometric properties of the BKZ solution). To understand the BKZ solution better we compute its Ricci curvature R(r) and its Kretschmann scalar K(r). With the help of Maple it turns out that so that this solution has a spacetime singularity at r = 0. A radially outgoing null geodesic satisfies the relation describe the radially outgoing null geodesics. They start at the singularity and escape to r = ∞, i.e., the singularity is visible for observers away from the singularity, and hence it violates the strong cosmic censorship conjecture. The concept of weak cosmic censorship is not applicable to this solution, since it is not asymptotically flat. According to the cosmic censorship hypothesis "naked" singularities, i.e., those which are visible to distant observers, should be "non-generic" and/or "unstable". Our analysis turns out to be nicely consistent with this hypothesis: regular steady states, by being close to the BKZ solution for large central redshift, seem to inherit this instability and are indeed unstable themselves, as we will show below.
Remark 3.8. The massless Einstein-Vlasov system has been considered in the literature, for example in [4,17,70]. The results of the present section, more precisely of the next subsection, show that the massless system captures the behavior of solutions of the massive system in a strongly relativistic situation, at least in the context of spherically symmetric steady states. The massless system plays a similar role in [2].
The size of the central redshift is not the only possible measure of the strength of relativstic effects in a given steady state. Another possible choice is the compactness ratio 2m(r) r , which by the Buchdahl inequality [12,3] is always less than 8 9 . We see from (3.43) that the compactness ration for the BKZ solution is 3 7 and it is therefore "far" from saturating the Buchdahl bound 8 9 .
So far, the discussion of the present subsection was restricted to the Einstein-Vlasov system, but in view of the relation between the corresponding steady states explained in Subsection 3.2 it applies equally well to the Einstein-Euler system with the equation of state (3.31), which we would refer to as the massless or radiative Einstein-Euler system.

The ultrarelativistic limit κ → ∞
In this subsection we prove that in a specified region close to the center the steady states provided by Proposition 3.2 are approximated by the BKZ solution when κ is sufficiently large. This is done in two steps. We first prove that the former steady states are approximated by solutions of the massless problem (3.40) with the same central value (3.41). In a second step we show that the behavior of the latter solutions is captured by the BKZ one.
For the first step we start with the Vlasov case. It is fairly simple to make the asymptotic relation between the functions g and h defined by (3.5) and (3.6) and their massless counterparts g * and h * precise: It should be noted that this is indeed a good approximation for large y, since all the terms on the left are then of order e 4y . This approximation would suffice to estimate the difference between solutions of (3.9), (3.10) and solutions of (3.40), (3.41). However, for the functions g and h defined by (3.13) and (3.14) in the Euler case the above estimate seems hard or impossible to obtain. For the latter case a different set-up must be used to obtain the desired asymptotics, but that set-up can be chosen such that it also works for the Vlasov case.
So from now on we treat both cases simultaneously and consider a steady state as provided by Proposition 3.2. For the moment we drop the subscript κ and note that by (3.18), which as we showed holds for the Vlasov case as well, We want to read this as a differential equation for p, and so we observe that the equation of state p = P (ρ) can be inverted to define ρ in terms of p: we notice that under our general assumptions the function P is one-to-one and onto on the interval [0, ∞[ both in the Vlasov and in the Euler case and both for the massive and the massless equations. We also notice that p determines all the other components of the corresponding steady state. Hence we consider the equation and its massless counterpart where according to (3.38) we set S(p) = 3p: In what follows, p κ and p * κ denote the solutions to (3.45) and (3.46) respectively, satisfying the boundary condition Remark 3.9. Strictly speaking we reparametrize our steady state family here and use the central pressure as the new parameter. However, the quantities y and p are in a strictly increasing, one-to-one correspondence in such a way that y → ∞ iff p → ∞, i.e., p(0) → ∞ is equivalent to y(0) → ∞. In view of (3.37) the above way of writing p(0) guarantees that at least for κ large this quantity asymptotically coincides with its original definition, which may justify the misuse of notation.
We now show that close to the origin and for large κ the behavior of p κ is captured by p * κ .
Lemma 3.10. There exists a constant C > 0 which depends only on Φ or P respectively such that for all κ > 0 and r ≥ 0, Proof. In order to keep the notation simple we drop the subscript κ and write p and p * for the two solutions which we want to compare. We introduce rescaled variables as follows: with α > 0. A straightforward computation shows that and notice that (3.50) is actually the same equation as (3.46). We choose In order to estimate the difference σ ′ − σ * ′ and apply a Gronwall argument we first collect a number of estimates. In these estimates C denotes a positive constant which depends only on the ansatz function Φ or P and which may change its value from line to line. In particular, such constants are independent of τ and κ. First we note that σ and σ * are decreasing, and hence for τ ≥ 0, By (3.16), This implies that A non-trivial issue is to get uniform control on the denominators in (3.49) and (3.50). But the spherically symmetric steady states of the Einstein-Vlasov or Einstein-Euler system which we consider here satisfy the estimate 2m(r) r < 8 9 , r > 0, (3.53) which is known as Buchdahl's inequality, cf. [12]; a proof for very general matter models which covers the present situation and also anisotropic steady states is given in [3]. Hence, Then these estimates together imply that Integration of this estimate yields Gronwall's lemma implies that = (e 4κ r 2 + e 8κ r 4 ) exp C e 4κ r 2 + e 8κ r 4 e 2κ , and the proof is complete.
In order to understand the behavior of p * κ in the limit κ → ∞ the following observation is useful; we recall (3.47) and Remark 3.9 to avoid miss-interpretation of the notation p * κ .
In view of Lemma 3.11 we now need to understand the behavior of p * 0 (r) for large r.
Proof. The proof relies on turning the massless steady state equations into a planar, autonomous dynamical system. We recall the Tolman-Oppenheimer-Volkov equation (3.18), which as we showed above holds in both the Euler and the Vlasov case, and combine it with the massless equation of state p = ρ/3 and (3.9). This leads to We rewrite this in terms of u 1 (r) = r 2 ρ(r), u 2 (r) = m(r)/r to obtain Now we multiply both equations with r and introduce w 1 (τ ) = u 1 (r), w 2 (τ ) = u 2 (r) with τ = ln r. Then We denote the right hand side of the system (3.54), (3.55) by F (w), which is defined and smooth for w 1 ∈ R and w 2 ∈] − ∞, 1/2[. The system has two steady states: We aim to show that the solution which corresponds to p * 0 converges to Z exponentially fast as τ → ∞. To this end we first observe that 4π −1 has eigenvalues so that Z is an exponential sink. On the other hand, with eigenvalues −1 and 2, stable direction (0, 1) and unstable direction (3, 4π). For the solution induced by p * 0 it holds that w(τ ) → (0, 0) for τ → −∞ which corresponds to r → 0. Since the corresponding trajectory lies in the first quadrant {w 1 > 0, w 2 > 0}, it must coincide with the corresponding branch T of the unstable manifold of (0, 0). We want to show that the trajectory T approaches the point Z. To this end, let D denote the triangular region bounded by the lines We want to show that T ⊂ D. Since the unstable direction (3, 4π) points into D, this is at least true for the part of T close to the origin. The line {w 1 = 0} is invariant and can not be crossed by T ; notice that this line is also the stable manifold of the point (0, 0). Due to Buchdahl's inequality, T must lie below the line {w 2 = 4/9}. Finally, along the line {w 1 = w 2 } the vector (−1, 1) is normal to it and points into the domain D. Since along this line, so that the vector field F points into the domain D, the trajectory T cannot leave this domain. Hence according to Poincaré-Bendixson theory, the ω-limit set of T must either coincide with Z, or with a periodic orbit. But according to Dulac's negative criterion, the set ]0, ∞[×]0, 1/2[ does not contain a periodic orbit, since In view of the real part of the eigenvalues λ 1,2 it follows that for any 0 < γ < 3/2 and all τ sufficiently large, i.e., all r sufficiently large, When rewritten in terms of the original variables this is the assertion.
We now combine the previous three lemmas to show that for κ large the steady state of the massive system which corresponds to p κ is approximated well by the BKZ solution on some interval close to the origin; ρ κ , m κ , µ κ , λ κ denote quantities induced by p κ . Proposition 3.13. There exist parameters 0 < α 1 < α 2 < 1 4 , κ 0 > 0 sufficiently large, and constants δ > 0 and C > 0 such that on the interval and for any κ ≥ κ 0 the following estimates hold: Proof. First we note that by Lemma 3.11, the radiative equation of state (3.38), and Lemma 3.12, so that together with Lemma 3.10, ≤ Ce 6κ r 4 + e 4κ r 6 exp C e 4κ r 2 + e 8κ r 4 + Ce −2κγ r −γ . (3.57) Below we will simplify the above right hand side on a suitably chosen r interval close to the origin, but first we estimate various other quantities in the same fashion. In order to estimate the difference in ρ we recall the rescaled variables introduced in (3.48) and the estimate (3.52). This implies that (3.58) If we combine this with the definition of τ and (3.56) it follows that ≤ Ce 2κ r 2 + Ce 6κ r 4 + e 4κ r 6 exp C e 4κ r 2 + e 8κ r 4 + Ce −2κγ r −γ where notice that E κ also bounds the right hand side in (3.57). The estimate (3.58) implies that m κ (r) r − m * κ (r) r ≤ Ce 2κ r 2 + Ce 6κ r 4 + e 4κ r 6 exp C e 4κ r 2 + e 8κ r 4 , and combining this with the second estimate in Lemma 3.12 and the scaling property in Lemma 3.11 it follows that m κ (r) integrating the ρ estimate directly yields the same result, but E κ is only integrable at the origin if we choose γ < 1. Next we use (3.7) to find that and using Buchdahl's inequality (3.53), Using (2.15) for µ ′ κ and µ ′ BKZ it follows that where we have again used Buchdahl's inequality (3.53); note that µ ′ BKZ refers to the BKZ solution (3.43). Finally, with the established bounds on r 2 ρ κ and e 2λκ above, equation . We now analyze the error term E κ (r) for r > 0 and κ > 0 such that where 0 < r 1 (κ) < r 2 (κ) and r 2 (κ) > 1. Then We choose r 1 (κ) = κ α 1 with some 0 < α 1 < 1 4 and Then for κ large, r 1 (κ) < r 2 (κ) and 1 < r 2 (κ) as required, and with δ := γα 1 and any α 2 ∈]α 1 , 1 4 [ the proof is complete.
In the above proposition no information is provided for µ κ itself. Such an estimate is derived next.
Proposition 3.14. Under the assumptions of Proposition 3.13 it follows that on the r interval specified in Proposition 3.13. Here C κ > 0 does depend on κ, but C > 0 does not. The point in this estimate is that the exponential terms on both sides converge to 1 as κ → ∞.
Proof. Clearly, Using Proposition 3.13 it follows that as desired, and the lower estimate is completely analogous.

Stability analysis for the Einstein-Vlasov system
In this section we fix some ansatz function Φ, but we need to strengthen the assumptions on Φ as follows.
with constants c > 0 and 1 ≤ k < 3/2; notice that we require the right hand side derivative at 0 to exist, but it need not vanish.
We consider a corresponding steady state (f κ , λ κ , µ κ ) as obtained in Proposition 3.2 with some κ > 0. We aim to prove that this steady state is linearly unstable when κ is sufficiently large. In order to keep the notation short, we write where we recall the definition of the particle energy and note the fact that the cut-off energy depends on κ as well; E 0 κ := e yκ(∞) = e µκ(Rκ) where R κ is the radius of the spatial support of the steady state, cf. Section 3.1. By (Φ2), The basic strategy on which our instability result relies is the following. The Einstein-Vlasov system conserves the ADM mass. The second variation of the ADM mass at the given steady state is a conserved quantity for the linearized system, restricted to linearly dynamically accessible states. The key to linear instability then is to prove that there is a negative energy direction for the linearized system, i.e., a linearly dynamically accessible state on which the second variation of the ADM mass is negative. This is the content of Theorem 4.3. This fact rather directly implies a linear, exponential instability result, cf. Theorem 4.8. However, much more precise instability information, including the existence of a simple growing mode, can be derived from Theorem 4.3. In order to do so the Hamiltonian character of the linearized system must be exploited which leads to the desired spectral properties of the generator of the C 0 group induced by the linearized system, cf. Theorem 4.28. In order to emphasize the basic mechanism leading to our instability result we try to bring in the more functional analytic, abstract tools and terminology only when they are finally needed to derive Theorem 4.28.

Dynamic accessibility and the linearized Einstein-Vlasov system
Sufficiently regular solutions of the spherically symmetric Einstein-Vlasov system conserve the ADM mass In addition, the Einstein-Vlasov flow conserves the Casimir functionals Here χ ∈ C 1 (R) is an arbitrary, prescribed function with χ(0) = 0, and λ f is defined by which is the solution of (2.6) induced by f and satisfying (2.3), cf. (3.7). In a stability analysis it is natural to restrict the admissible perturbations of the given steady state to such as preserve all the Casimir invariants. These dynamically accessible perturbations form the symplectic leaf S fκ through the given steady state. Its formal tangent space at f κ is the set of linearly dynamically accessible perturbations [40]. In order to proceed we need to recall the usual Poisson bracket of two continuously differentiable functions f and g of x, v ∈ R 3 . The product rule for the Poisson bracket reads {f, gh} = {f, g}h + {f, h}g.
In addition we denote the radial component of a vector v ∈ R 3 by An important subclass of linearly dynamically accessible states are of the form the operator B κ is defined in Definition 4.11, the concept of linearly dynamically accessible states which we adopt later is made precise in Definition 4.14. The generating function h = h(x, v) of the perturbation should be spherically symmetric and C 1 . In this case, we denote the perturbed metric field λ by λ h and we have On linearly dynamically accessible states of the form (4.5) the second variation of the ADM mass H ADM takes the form see [32].
Remark 4.1. We emphasize that the above variational structure and the discussion of dynamic accessibility is limited to radially symmetric perturbations; a detailed derivation is given in [32]. However, in our stability analysis of the steady states of the EV-system variational methods do not play a role.
The relation of the quadratic form A κ to the non-linear Einstein-Vlasov system is not relevant for the present analysis, but it will become so in a possible extension of the linear instability result to the non-linear system. For the present purposes it is important that A κ is conserved along the linearized flow, which we need to recall next for the special class of perturbations of the form (4.5). The dynamics of the generating function h is determined by where λ h and µ h are determined by together with the boundary conditions λ h (t, 0) = 0 = µ h (t, ∞), and It can be checked from (4.9) that λ h is indeed given by the formula (4.6), cf. [32]. Let us denote so that D κ is the support of the steady state under consideration. A simple iteration argument the details of which are indicated in [32] shows that for anyh ∈ C 1 (D κ ) there exists a unique solution h ∈ C 1 ([0, ∞[; C(D κ ))∩C([0, ∞[; C 1 (D κ )) of the above system with h(0) =h; in passing we note that Eqn. (5.10) in [32], which corresponds to (4.8), contains a sign error. Moreover, we note that the characteristic flow of (4.8) is the one of the steady state solution and leaves the set D κ invariant. Under the present regularity assumptions on Φ this characteristic flow is C 3 , and the limiting factor concerning regularity of solutions to the above system is the factor φ ′ κ in (4.5). We recall from [32] that the energy A κ (h, h) is conserved along solutions of the linearized system stated above. In [33] it was shown that for κ sufficiently small the quantity A κ is positive definite, which leads to a linear stability result. Here we aim for the opposite.

A negative energy direction for κ sufficiently large
Our next aim is to show the existence of a linearly dynamically accessible perturbation for which the bilinear form A κ is negative. Theorem 4.3. There exists κ 0 > 0 such that for all κ > κ 0 there exists a spherically symmetric function h ∈ C 2 (R 6 ) which is odd in v, such that where f h is given by (4.5).
For the proof of this result we first establish two auxiliary results. It will be convenient to use, in addition to the radial momentum component w = x · v/r, also the modulus of the tangential momentum component, i.e., L = r 2 ω 2 is the modulus of angular momentum squared mentioned in the introduction.
Lemma 4.4. The following identities hold: Proof. First we observe that all the integrands depend only on w and ω so we can change variables using dv = 2π ω dω dw.
With E κ = e µκ √ 1 + w 2 + ω 2 it follows that Using the first identity for the integrals containing w 2 or w 4 and the second one for the others, the assertions follow by integration by parts with respect to w or ω and the definitions of ρ κ and p κ ; no boundary terms arise since φ κ (E κ ) vanishes on the boundary of its support and the powers of ω vanish at ω = 0.
Lemma 4.5. On the interval [r 1 κ , r 2 κ ] introduced in Proposition 3.13 and for κ sufficiently large, Proof. The same change of variables which yields the formulas (3.4) implies that the point here is that only the factor e 2yκ appears, while for p κ (and ρ κ ) the factor is e 4yκ ; recall also that φ the constant C is finite, since by Assumption (Φ2) the last integral converges. On the other hand, Proposition 3.13 implies that on the interval [r 1 κ , r 2 κ ] and for κ sufficiently large, with α 2 > 0 from Proposition 3.13. Comparing the two estimates yields the assertion.
Proof of Theorem 4.3. We shall look for a negative energy direction h of the form where g ∈ C 2 ([0, ∞[) is to be chosen. We recall the definition of the coordinate w (4.4). As a function of v, h is clearly odd. With (4.18) and (4.6) in mind we can simplify the expression for λ h : where we used (4.15) and (3.17). On the other hand, According to Proposition 3.13, on the spatial interval [r 1 κ , r 2 κ ] the steady state (f κ , λ κ , µ κ ) is well approximated by the BKZ solution of the massless Einstein-Vlasov system, provided κ is sufficiently large. Hence we localize the perturbation h given by (4.18) to this interval by setting where 0 ≤ χ κ ≤ 1 denotes a smooth cut-off function supported in the interval [r 1 κ , r 2 κ ] and identically equal to 1 on the interval [2r 1 κ , r 2 κ /2]; note that the latter interval is non-trivial for κ sufficiently large. In addition, we require that Then the perturbation f h takes the form We plug (4.20) into (4.7) and obtain the following identity: The idea now is as follows: If we replace the steady state (f κ , λ κ , µ κ ) in A 1 by the corresponding limiting quantities according to Proposition 3.13, then a strictly negative term arises, together with various error terms, which like A 2 and A 3 do not destroy the overall negative sign of A κ (h, h), provided κ is sufficiently large.
Since by Proposition 3.13 1 rµ ′ κ ≈ 2 on the support of χ κ , we split A 1 further: Recalling the definition of ω some algebra implies that

If we apply Lemma 4.4 and Lemma 4.5 and collect terms it follows that
in the last estimate we used the fact that 3p κ ≤ ρ κ . Using (3.17) it follows that In order to estimate the various parts into which A κ (h, h) has now been split we observe that by Proposition 3.13 the following estimates hold, provided κ is sufficiently large: Hence using the fact that χ κ = 1 on [2r 1 κ , r 2 κ /2] and Proposition 3.14, (4.25) Here and in what follows C κ always denotes the constant introduced in Proposition 3.14, while C denotes a generic constant which does not depend on κ, but which may change its value from line to line.
We now proceed to prove that all the remaining terms are smaller in modulus than the negative term just obtained, provided κ is sufficiently large. In order to estimate A 112 we note that by Proposition 3.13, Combining this with (4.24) and Proposition 3.14, In order to estimate A 12 we first observe that Using Lemma 4.4 we find that We insert the above two estimates into the definition of A 12 , and using (4.24) and Proposition 3.14 we conclude that To estimate A 2 we notice that we use the bounds (4.19) for χ ′ κ and the estimate cf. (4.15) and (3.17). Together with (4.24) this implies that In order to estimate A 3 we first estimate the v-integral contained in that term: where we used (4.15) and (3.17). Hence We add up (4.25), (4.26), (4.27), (4.28), and (4.29) to find that besides the constants C κ and C introduced in Proposition 3.14 there exist positive constants C 1 , C 2 , C 3 such that provided κ is sufficiently large, and the proof is complete. ✷ The perturbation f h is carefully engineered to produce such a gain, taking advantage of the exact κ → ∞ asymptotics provided by Proposition 3.13.

Linear exponential instability
Theorem 4.3 is sufficient for showing a linear exponential instability result, without proving the existence of an exponentially growing mode. We present this simple argument first, before we turn to the growing mode in the following sections. In order to exploit Theorem 4.3 we further analyze Eqn. (4.8), which governs the dynamics of the linearly dynamically accessible perturbations. We split h into even and odd parts with respect to v, i.e., h ). It is then easy to see that λ h = λ h − , µ h = µ h − , and (4.8) can be turned into the system where we define the operator T essentially represents the transport along the characteristic flow of the steady state under consideration. In what follows we use the weight W := e λκ |φ ′ κ (E κ )| and denote by L 2 W = L 2 W (D κ ) the corresponding weighted L 2 space on the set D κ , cf. (4.13); (·, ·) L 2 W denotes the corresponding scalar product. We obtain the following linear, exponential instability result: Theorem 4.8. There exist initial datah + ,h − and constants c 1 , c 2 > 0 such that for the corresponding solution to the system (4.30), (4.31), Proof. The proof relies on the fact that solutions to (4.30), (4.31) preserve energy, i.e., and on the virial identity 1 2 (4.34) As to (4.33) we recall from [32] that the energy A κ (h, h) is conserved, but if we substitute h = h + + h − and use the fact that λ h = λ h − it follows easily that here we used the anti-symmetry of T stated in Lemma 4.10 below. The identity If we substitute this into the above expression and use the linearized field equations (4.9) and (4.10) to replace the terms ρ h and p h it follows by a straightforward computation that If we integrate the µ ′ h λ h term by parts the r integral vanishes and (4.34) follows. To prove the linear exponential instability we follow the approach used by Laval, Mercier, Pellat [44] in the context of stability of plasma flows. Since Jeans' Theorem does not hold for the Einstein-Vlasov system, cf. [65], we must modify the argument of [44]. (4.35) and to obtain initial data for the system (4.30), (4.31) we supplement this bẙ

By Theorem 4.3 and Remark 4.2 there exists an odd-in-v functionh
where ε > 0 is chosen so small that Conservation of energy (4.33) and (4.36) imply that By the Cauchy-Schwarz inequality, Combining this with the virial identity (4.34) and (4.37) it follows that we see that y(0) = 0 anḋ and thus y(t) ≥ c ǫ t. By definition of y this is the assertion for h − (t) L 2 W . By Cauchy-Schwarz, which implies the estimate for T h + (t) L 2 W , and the proof is complete.
Remark 4.9. It should be noted thath + defined in this proof is indeed even so thath ± are the even and odd parts of some initial datah for the original linearized system (4.8)-(4.12).
Proof. Clearly, For fixed ǫ > 0 we now integrate by parts: all the other terms which appear actually drop out, because they contain either ∂ and we may temporarily replace χ by a smooth approximation. But the outer unit normal ν = (ν x , ν v ) which appears in the boundary integral is parallel to (∇ x E κ , ∇ v E κ ) so that the boundary integral vanishes as well, and the assertion follows.
It should be noted that the above integration-by-parts formula does not require any boundary conditions on either g, h or on the weight χ.

The Hamiltonian structure of the linearized Einstein-Vlasov system
We begin by defining a Hilbert space which is the suitable state space for the linearized dynamics, and a certain transport operator.
Functions f ∈ X are extended by 0 to all of R 6 .
(b) For a function f ∈ X the transport term T κ f := −e −λκ {f, E κ } exists weakly in X, iff there exists a function η ∈ X such that for all test functions ψ ∈ C ∞ c (D κ ), and T κ f := η for f ∈ D(T κ ).
(c) The operator B κ : and where we recall the definition of w (4.4). (b) The fact that the operator T κ is anti-self-adjoint (or skew-adjoint) is shown in detail in [63].
In what follows we denote perturbations of f κ , µ κ , λ κ by f, µ, λ respectively, i.e., f κ + f, µ κ + µ, λ κ + λ is the solution to the Einstein-Vlasov system (2.5)-(2.13), which is then linearized in f, µ, λ. In particular, λ = λ f depends on f through the non-local relationship Remark 4.13. A simple application of the Cauchy-Schwarz inequality implies that λ is well-defined for any f ∈ X, i.e. there exists a C > 0 such that r 0 s 2 f v dv ds ≤ C, r ∈ [0, ∞[, as f is supported on D κ . In particular |λ(r)| ≤ C r and thus ∞ 0 λ(s) 2 ds < ∞. As shown in [32] the formal variation of (4.2) yields the relation Therefore, the formal tangent space of linearly dynamically accessible perturbations can be interpreted as the set of perturbations f satisfying the orthogonality condition for all χ ∈ C 1 (R), χ(0) = 0. Here λ is defined through (4.41). By Proposition 3.2 in [32] for any h ∈ C 1 (D κ ) the orthogonality condition (4.42) is satisfied if f = B κ h. A somewhat lengthy, but standard density argument shows that (4.42) is true for any f ∈ R(B κ ). This motivates the following definition.  In the context of the Vlasov-Poisson system, the analogous statement is true and can be inferred from [29]. The above question is intimately tied to the validity of Jeans' theorem; if the former is true then one can indeed show that R(B κ ) = R(B κ ). However, Jeans' theorem is known to be false in general for the Einstein-Vlasov system [65].
(b) We shall see in Lemma 4.20 that the operator B κ and the Hilbert space X arise naturally in the study of the linearized Einstein-Vlasov system around the steady state (f κ , µ κ , λ κ ).
A key ingredient in the following analysis is a modified potential induced by a state f ∈ X. with λ defined by (4.41).
(b) The operatorL κ : X → X is defined bȳ As it shall be often used in the remainder of this paper, we observe that this follows since µ ′ κ ≥ 0. in the weak sense.
Proof. The properties of the steady state imply that the quantity e µκ+λκ (2rµ ′ κ + 1) is bounded; this is seen from the field equation (2.7) and the boundedness of p κ , µ κ , λ κ which is explained in Section 3.1. For f ∈ X the Cauchy-Schwarz inequality implies that which implies the estimate forμ. Its continuity is obvious, and since λ is continuous,μ is continuously differentiable for r > 0 with which follows after multiplying (4.43) by e µκ+λκ and then differentiating with respect to r. Eqn. (4.46) follows directly from (4.49). By Remark 4.13 and (4.45) we conclude from (4.46) that d dr e µκ+λκμ ∈ L 2 r . Since µ ′ κ + λ ′ κ = 0 outside D κ , it follows thatμ ∈Ḣ 1 r . Multiplying both sides of (4.46) by e −2λκ r and taking one more radial derivative we obtain (4.47); note that r 2 ρ ∈ L 1 ([0, ∞[).
Proof. Boundedness follows from Lemma 4.17 and self-adjointness is obvious from an integration-by-parts argument and (4.47). To prove compactness, assume that (f n ) n∈N ⊂ X converges weakly to some f ∈ X. From (4.41), If we denote C n (s) := e 2λκ f n − f, e −λκ |φ ′ κ (E κ )| v 1 Bs(0) X , it follows from the weak convergence of (f n ) n∈N that lim n→∞ C n (s) = 0 for any s ≥ 0. Moreover for any s ≥ R κ , C n (s) = C n (R κ ), where [0, R κ ] is the spatial support of the steady state. For any s, s ′ ∈ [0, R κ ], where the last bound follows from Hölder's inequality. In particular, the sequence of functions (C n ) n∈N is equicontinuous on [0, R κ ]. By the Arzela-Ascoli theorem it follows that along some subsequence (f n k ) k∈N , lim k→∞ c n k = 0, where c n := sup s≥0 |C n (s)| = max s∈[0,Rκ] |C n (s)|. From this we conclude that This shows the compactness of the map K.
Lemma 4.19. The operator B κ defined in Definition 4.11 is densely defined and anti self-adjoint on X. The operatorL κ is bounded and symmetric on X.
Proof. We write B κ = T κ + R κ where for f ∈ X, notice that f is extended by 0 to all of R 6 . Since the functions r, µ κ , λ κ , φ ′ (E κ ), v are all bounded on the set D κ a straightforward computation shows that R κ is bounded on X, and its anti-symmetry is just as easy to check. Hence the assertion on B κ follows from Remark 4.12. The estimate forμ in Lemma 4.17 implies thatL κ is bounded on X, and it is easy to check thatL κ is symmetric.
We now establish a first order Hamiltonian formulation of the linearized Einstein-Vlasov system.
Lemma 4.20. The linearized, spherically symmetric Einstein-Vlasov system takes the form where the operators B κ andL κ are defined in (4.40) and (4.44) respectively. Furthermore D(B κLκ ) = D(B κ ). The operatorL κ induces a quadratic form on X which takes the form where λ depends on f through (4.41). The flow of (4.50) preserves A κ (f, f ).
Proof. First we notice that for f ∈ X at least formally, where recall that the coordinate w is given by (4.4). Since λ f is given by (4.41) and the bound (4.48) holds, we conclude that |λ f (r)| ≤ Cr 1/2 . This bound together with (4.46) implies that |μ ′ (r)| ≤ Cr −1/2 , and the right hand side of (4.52) lies in X. HenceL κ maps D(B κ ) = D(T κ ) into itself, and the operator product B κLκ is defined on D(B κ ).
To linearize the equations we expand f = f κ +ǫδf, µ = µ κ +ǫδµ, λ = λ κ +ǫδλ into (2.5), use (2.8), and neglect terms of quadratic order and higher. By a slight abuse of notation we rename δf, δλ, δµ into f, λ, µ respectively, and obtain the linearized Vlasov equation From the field equation (2.7) and the definition (2.11) of p it follows that Plugging this back into (4.53) we obtain the equation where B κ is given by (4.40). Using (4.52), (4.15), and (4.49) it follows that Plugging this into (4.54) we obtain which proves (4.50). The conservation of (L κ f, f ) X is a direct consequence of the antisymmetry of B κ and (4.50). Finally, from the definition ofL κ we have r 2μ e µκ+λκ ρ dr.
The above first order Hamiltonian formulation appears to be new and not exploited elsewhere in the literature. We now derive a second order formulation of the linearized Einstein-Vlasov system, which was first given in 1968 by Ipser and Thorne [36]. Lemma 4.21. A formal linearization of the spherically symmetric Einstein-Vlasov system takes the form where f − denotes the odd part of f with respect to the variable v,

56)
and B κ is defined in (4.40). The operator L κ is self-adjoint on X with the latter being defined in analogy to Definition 4.11 (b). Therefore ∂ t f − 2 X +(L κ f − , f − ) X is formally conserved along the flow of (4.55), and Proof. We begin by examining the operator L κ . First we note that B κ B κ f ∈ X is defined in a pointwise, classical sense, provided f ∈ C ∞ c (D κ ); here the regularity assumption Φ ∈ C 2 (]0, ∞[), i. e., φ κ (E κ ) ∈ C 2 (D κ ) enters. Thus B κ B κ is densely defined and selfadjoint by arguments analogous to those in Remark 4.12. If we define then it follows easily that R κ is bounded on X and symmetric, which proves the assertion for L κ . Since φ ′ κ (E κ )E κμ is even in v it is immediate from (4.50) that We take the time derivative of the second equation and use the relation which is the linearization of the field equation (2.8). Thus Next, here we recall (4.57) and observe that ∂ t λ = λ ∂tf . Therefore by (4.51) and the fact that f + and f − are orthogonal to each other in the space X.
Note that by (4.50), ∂ t f ∈ R (B κ ), the range of the operator B κ .
(b) As noted at the beginning of Section 4.1 the quadratic form (4.51) arises as the second variation of a suitably chosen energy-Casimir functional [32,33,40], which is conserved along the flow of the nonlinear Einstein-Vlasov system.

The reduced Einstein-Vlasov operator
The main result of this section is Theorem 4.24, which states in a precise manner that the operator L κ is bounded from below and above by a simpler, Schrödinger-type operator. Note that this operator acts only on functions that depend on the spatial variable r and we therefore refer to it as a macroscopic operator. In Section 4.6 we shall use this to obtain refined spectral information about the operator L κ . An attempt to derive such a reduced operator was made by Ipser in [37], yet the operator derived there bounds the operator L κ from below only under an additional assumption on the steady state, which appears to be hard to verify. Our approach is new and the key ingredient is the use of the modified potentialμ introduced in (4.43).
The operator ∆ κ is given by The reduced operator is given by The non-local reduced operatorS κ is given bỹ where Π denotes the projection onto the orthogonal complement of R(B κ ) in X, and id is the identity. On a flat background, i.e., for λ κ = µ κ = 0 the operator 4π∆ κ is the Laplacian applied to spherically symmetric functions.
To make sense classically of the above operators it is natural to require ψ ∈ H 2 r = H 2 (R 3 ) ∩ L 2 r , where L 2 r denotes the space of all spherically symmetric L 2 functions. However, this is too strong for our purposes, since the modified potentialsμ f from Definition 4.16 must belong to the domain of definition of S κ andS κ , but they in general belong toḢ 1 r . Therefore we define S κ andS κ using duality: Definition 4.23 (The operators S κ andS κ as symmetric quadratic forms). Let (Ḣ 1 r ) ′ denote the dual space ofḢ 1 r , the radial subspace of the homogeneous Sobolev spaceḢ 1 (R 3 ). The dual pairing is denoted by ·, · . Then S κ : The operatorS κ :Ḣ 1 r → (Ḣ 1 r ) ′ is defined analogously.
We use a standard definition ofḢ 1 as given e.g. in [26], and we recall that D κ = supp f κ . As defined above the operators S κ andS κ are clearly self-dual. Moreover, taking the r derivative of ρ κ = φ κ (E κ ) v dv it follows that Thus the operator S κ can formally be written in the form The key result of this section is the following theorem.
For every f ∈ X andμ f as defined in (4.43), For every f ∈ X, Proof. We first observe that for µ ∈Ḣ 1 r , in particular, using (4.43) and (4.46), and on the other hand by (4.47), The definitions of A κ and S κ together with (4.68)-(4.69) imply that and the lower bound (4.64) is proven. To obtain the upper bound (4.63) we observe that by (4.67) for µ ∈ H 2 r and λ = λ f with f ∈ X, Together with (4.70) this implies that and (4.63) is proven. We now turn to part (b) and assume that f ∈ R(B κ ). By definition of the projection Π and (4.68)-(4.69), and Using (4.70) and (4.71) it follows that which proves (4.66). To obtain the upper bound (4.65 An argument analogous to the one used for (4.70) implies that and the proof is complete.

Main results for the Einstein-Vlasov system
We first need to comment on some notation. Let L be either a self-adjoint, linear operator L : H ⊃ D(L) → H on some Hilbert space H with scalar product (·, ·) H or a self-dual operator L : H → H ′ . In either case we we denote by n − (L) its negative Morse index, which by definition is the maximal dimension of subspaces of H on which (L·, ·) H < 0 or L·, · < 0 in the self-dual case with ·, · the dual pairing. For any decomposition H = H − ⊕ ker L ⊕ H + such that (L·, ·) H < 0 on H − \ {0} and (Lu, u) H ≥ δ u 2 H on H + it can be shown that dim H − = n − (L), see Remark 2.2 and Lemma 12.1 in [48]. For any subspace S ⊂ H, we also use n − (L| S ) to denote the Morse index of L restricted to S.
The self-adjointness of the operator L κ and Theorem 4.24 imply that the negative part of the spectrum has to be discrete and that part (a) of Theorem 1 holds, more precisely: Theorem 4.25. The negative part of the spectrum of L κ : X ⊃ D(L κ ) → X is either empty or consists of at most finitely many eigenvalues with finite multiplicities. For κ sufficiently large, there exists at least one negative eigenvalue and therefore a growing mode for the linearization of the Einstein-Vlasov system around (f κ , λ κ , µ κ ).
Proof. By Lemma 4.21 the operator L κ is self-adjoint. We note that for f ∈ D(L κ ), for the sake of completeness we recall the general argument behind the first two estimates in Lemma A.1 in the appendix. In order to show that n − (S κ ) < ∞ we first observe that S κ φ, φ ≥ C S ′ κ φ, φ for φ ∈Ḣ 1 r , where the self-dual operator S ′ κ :Ḣ 1 r → (Ḣ 1 r ) ′ is formally given as S ′ κ = −∆ − V κ with a non-negative, continuous potential V κ which has compact support. This follows from the corresponding bounds on λ κ and µ κ . Now we observe that the mapping (−∆) 1/2 :Ḣ 1 (R 3 ) → L 2 (R 3 ), φ → (2π|ξ|φ)ˇis an isomorphism which respects spherical symmetry. Passing to ψ = (−∆) 1/2 φ we arrive at the relation that is compact, since V κ is bounded and supported on the compact setB Rκ (0), and the mappinġ so thatĥ ∈ L ∞ ∩ L 2 (R 3 ), and hence 1 2π|ξ|ĥ and its inverse Fourier transform are in L 2 (R 3 ). The spectral properties of compact operators imply that n − (id − K) < ∞, and invoking Lemma A.1 again it follows that n − (S κ ) < ∞ as claimed.
The assertion on the spectrum of L κ follows from its spectral representation; for the sake of completeness we include the argument in Proposition A.2 in the appendix.
Each negative eigenvalue of L κ gives rise to a pair of stable and unstable eigenvalues for (4.55). For κ sufficiently large, n − (L κ ) ≥ 1 by Theorem 4.3, and the claim follows. Note that the operator L κ is non-negative when restricted to the subspace of X of all even-in-v functions in X, since for any f ∈ X, ) is the even part of f . In particular, since the negative eigenvalues are of finite multiplicity, the associated eigenspace necessarily contains only odd-in-v functions.
Remark 4.26. The finiteness of n − (L κ ) follows alternatively from the observation that the difference id −L κ = K, which is a compact operator by Lemma 4.18. SinceL κ is a compact perturbation of the identity, the only possible accumulation point is 1 and therefore n − (L κ ) < ∞.
Remark 4.27. By the same arguments as in the proof above we can conclude that which will be used below.
We now come to the existence of the linearized flow and the corresponding exponential trichotomy decomposition into a stable, unstable, and center space. Under a nondegeneracy assumption on the operatorS κ we in fact obtain a Lyapunov stability of the flow on the center space for any κ > 0. The key in our analysis will be the first order formulation (4.50).
Theorem 4.28. The operator B κLκ generates a C 0 group e tBκLκ of bounded linear operators on X, and there exists a decomposition with the following properties: is the linear subspace spanned by the eigenvectors corresponding to positive (negative) eigenvalues of B κLκ , and (ii) The quadratic form L κ ·, · X vanishes on E u,s , but is non-degenerate on E u ⊕ E s , and where (v) When kerS κ = {0}, then E c can be further decomposed as (b) We expect the exponential trichotomy estimates in (iii) to become useful for the possible construction of invariant (stable, unstable and center) manifolds, as well as for a possible proof of nonlinear instability in the future.
(c) In the Newtonian limit the operatorS κ formally converges to its Newtonian counterpart, which was proved to be positive e.g. in [30]. This can be shown rigorously and for κ ≪ 1 it can be then shown thatS κ ≥ 0. By (4.78) we therefore obtain linear stability against general initial data in X. This improves the result in [32], where linear stability was proved for dynamically accessible initial data, i.e., when f ∈ R (B κ ).
To prove Theorem 4.28 we first need to prove several preparatory statements.
Proposition 4.30. There exists a decomposition of X into the direct sum of three closed subspaces X = X − ⊕ kerL κ ⊕ X + with the following properties: Proof. By Lemma 4.18 and the definition (4.44) ofL κ the operatorL κ is a compact perturbation of the identity, and therefore the only possible accumulation point of its spectrum is 1. Therefore the negative Morse index and the dimension of kerL κ are finite and the claimed decomposition follows, as well as the bound (H2).
In the study of linear stability of (4.50), it is important to restrict to the dynamically accessible space R (B κ ), for which we have the following result.
Lemma 4.31. The quadratic form L κ ·, · X is non-degenerate on the quotient space and it implies that f ∈ ker L κ ·, · X | R(Bκ) if and only if f ∈ R(B κ ) and here Π is the projection onto ker B κ in X. In view of the the definition of µ by (4.47), we deduce from (4.83) thatS κ µ = 0. On the other hand, ifS κ µ = 0, then it is easy to check that which together with (4.81) implies (4.82). Now we are ready to prove Theorem 4.28.
Proof of Theorem 4.28. To prove (i) we note that each pair of stable and unstable eigenvalues for the operator L κ provided by Theorem 4.25 also gives a pair of stable and unstable eigenvalues for (4.50). Together with (4.82) we obtain the dimension formula (4.73).
Consider the Hamiltonian formulation (4.50). The operator B κ can be defined as a self-dual operator X ′ ⊃ D (B κ ) → X. Moreover, by Proposition 4.30, the quadratic form L κ ·, · X satisfies the assumptions (H1)-(H3) in [48]. Thus the conclusions (ii)-(iv) follow directly from [ Since ker B κLκ is the steady solution space of (4.50) and is therefore contained in E c , we get the decomposition (4.77) for E c . Since [48,Cor. 6.1] and E c is the orthogonal complement of E s ⊕ E u in X, it follows that ThusL κ | E c ∩R(Bκ) ≥ δ 1 > 0, and by using L κ f, f X as the Lyapunov functional, we get the . Then the stability estimate |e tBκLκ | E c | ≤ M follows by the decomposition (4.77) and the fact that e tBκLκ f = f when f ∈ ker B κLκ . WhenS κ > 0, (4.78) follows from the decomposition (4.84) by the same arguments as above. This finishes the proof of the theorem.
Remark 4.32. (a) The polynomial growth on E c can be shown ( [49]) to be at most quadratic by using the second order formulation (4.55).
(b) It is tempting to find the most unstable eigenvalue λ 0 < 0 of L κ by trying to minimize (L κ f, f ) X over the constraint set {f ∈ D (L κ ) | f X = 1}. However, in the current case it is difficult to study this variational problem directly due to the lack of compactness. It is the self-adjointness of L κ and the finiteness of the negative Morse index that are crucial to the proof of the existence of unstable eigenvalues. and in particular, dN dρ is a positive and strictly decreasing function on ]0, ∞[. By (P1), P (s) ≤ cs γ for s ∈ [0, 1], where γ > 1. Hence for 0 < ρ ≤ 1, Together with (5.1), (5.2), and the monotonicity of dN dρ this implies that lim ρց0 N (ρ) = 0 and lim ρց0 Since the functional relationship between the pressure and the energy density ρ is prescribed by (1.8), the expression p(n) simply means P (ρ(n)). Using (5.3), the Tolman-Oppenheimer-Volkov equation (3.18) and the fact that is the support of the steady state. We integrate this differential equation to find that for all 0 ≤ r, r 0 < R κ , dN dρ (ρ κ (r)) = dN dρ (ρ κ (r 0 ))e µκ(r)−µκ(r 0 ) .

Now we observe that the function
with c κ > 0 has the same properties as the function N stated above, and we can choose the normalization constant c κ such that dNκ dρ (0) = e µκ(Rκ) . For the given steady state and with n κ = N κ (ρ κ ) the identity holds on the interval [0, R κ ].
Remark 5.1. The normalization can not be achieved by taking r → ∞, since outside the support of the steady state the middle term is necessarily constant, but µ κ is not.
Before we linearize the Einstein-Euler system we collect a few properties of the steady state under consideration.

Linearization
Let us now linearize the Einstein-Euler system (2.14)-(2.19) around a given steady state (n κ , u κ ≡ 0, λ κ , µ κ ). We write n, u, ρ, p, λ, µ for the Eulerian perturbations, i.e., n κ +n, u κ + u, ρ κ + ρ etc. correspond to the solution of the original, non-linear system. Linearizing ρ κ + ρ = N −1 κ (n κ + n) and p κ + p = P (N −1 κ (n κ + n)) we have the relations where we have used (5.9). Multiplying by e λκ we find that where we have used the field equation (2.15) in the third line to conclude that (5.17) In order to proceed we define a suitable phase space and a modified potential in complete analogy to Definition 4.16.
Definition 5.3. (a) The Hilbert space X 1 is the space of all spherically symmetric functions in the weighted L 2 space on the set B κ = B Rκ (the ball with radius R κ which is the support of ρ κ ) with weight e 2µκ+λκ Ψ κ and the corresponding inner product, X 2 is the space of radial functions in L 2 (B κ ), and the phase space for the linearized Einstein-Euler system is X := X 1 × X 2 .
(b) For ρ ∈ X 1 the induced modified potentialμ is defined as (c) The operatorsL κ : X 1 → X ′ 1 and L κ : X 1 → X 1 are defined bỹ L κ ρ := e 2µκ+λκ Ψ κ ρ + e µκ+λκμ ρ . (5.20) Here the dual pairing is realized through the L 2 -inner product, so that To check that the operatorsL κ and L κ are well-defined and self-dual or self-adjoint respectively, it clearly suffices to show only the self-adjointness of L κ .
Lemma 5.4. The operator L κ defined by (5.21) is well-defined, self-adjoint, and L κ − id is compact.
Proof. The operator L κ is clearly symmetric, i.e., We define the operator K : X 1 → X 1 by Kρ = 1 e µκ Ψκμ . It is straightforward to check that the modified potentialμ has the properties stated in Lemma 4.17, in particular Evaluating the L 2 (R 3 )-inner product of (5.22) with e µκ+λκμ , we get Therefore and thus ∇ e µκ+λκμ L 2 ≤ C ρ X 1 . Here we made use of the compact support of 1 Ψκ and the embeddingḢ 1 R 3 ֒→ L 2 (B κ ). Since and the mappingḢ 1 r R 3 ∋ ρ → ρ |S ∈ L 2 (S) is compact for any compact set S ⊂ R 3 , the operator K is compact. Therefore, L κ = id + K is self-adjoint on X 1 .
With the definition ofL κ the linearized velocity equation (5.16) takes the forṁ It easy to check that formally the energy is conserved along solutions of (5.15), (5.23). The functional I maps X 1 × L 2 e 3λκ nκ (B κ ) to R. From the definition (5.20) one readily checks that for ρ ∈ X 1 , In order to exhibit the Hamiltonian structure of this system it is convenient to define the modified divergence operator A κ : κ v , (5.27) and its dual operator (the modified gradient) A ′ κ : Then for any v ∈ X 2 and ρ ∈ X ′ 1 it holds that where ·, · is the dual bracket. We also define the adjoint operator A * κ : Then for any v ∈ X 2 and ρ ∈ X 1 We claim that both A κ and A ′ κ are densely defined and closed. The density is clear, and to show that they are closed we consider the isometries Then to show that A κ : The adjoint of A κ is given by Clearly A * κ is also densely defined. To show that A κ is closed it suffices to show A κ = A κ * * . To justify the latter we need to show that the adjoint of A * κ is precisely A κ , which is a simple consequence of integration-by-parts since by (5.8) n κ Ψ κ = dP dρ (ρ κ ) ∼ ρ≪1 ρ γ−1 κ and thus the boundary term vanishes.
Remark 5.5 (Notational conventions). We adopt the notation from Yosida's book [73]: For an operator A : X → Y between two Hilbert spaces X and Y , we use A ′ : Y ′ → X ′ to denote the dual operator and A * : Y → X to denote the adjoint operator of A. The operators A ′ and A * are related by where I X : X ′ → X and I Y : Y ′ → Y are the isomorphisms defined by the Riesz representation theorem.
Lemma 5.6. The formal linearization of the spherically symmetric Einstein-Euler system takes the Hamiltonian form d dt where (ρ, v) ∈ X, and Then J κ : X ′ → X and L κ : X → X ′ are anti-self-dual and self-dual respectively. The conserved energy functional takes the form Alternatively, one can rewrite (5.29) in the second-order form Proof. Eqn. (5.29) is just a rephrasing of (5.25), (5.26). The formal conservation of energy t → I(ρ, v)(t) is a direct consequence of the Hamiltonian structure (5.29). That J κ is anti-self-dual follows easily from the definition of A κ and A ′ κ . The self-duality of L κ follows from the self-duality ofL κ , which in turn, is a trivial consequence of the self-adjointness of L κ . Finally, (5.31) follows formally by applying d dt to the v-component of (5.29) and then plugging in the equation satisfied by the ρ-component of (5.29). Equation (5.32) is just a restatement of (5.32).
We will show below that the operator A * κ L κ A κ is in fact self-adjoint, cf. Lemma 5.23.

Variational approach to stability
The variational picture explained in this section is naturally constrained to radially symmetric perturbations. However, by comparison to the corresponding variational picture for the EV-system in Section 4.1, we need to provide more details, as it will play an important role in our formulation of the main results for the Einstein-Euler system in Section 5.5. The fundamental conserved quantities associated with the Einstein-Euler system are the total ADM mass M (ρ) = 4π ∞ 0 r 2 ρ(r) dr (5.33) and the total particle (or baryon) number We define the open set Since by Cauchy-Schwarz X 1 ⊂ L 1 (B κ ), it follows that both M and N κ are well-defined functionals on U .
Since M is a linear map it is clearly twice Fréchet differentiable and by a standard argument so is N κ . In order to compute the first and the second Fréchet derivative of these conserved quantities we will as above replace ρ by ρ κ + ρ etc., i.e., ρ now stands for the corresponding perturbation. Since by definition all perturbations ρ belong to X 1 we have supp ρ ⊂ B κ .
In order to compute I 2 we recall (5.6) and note that ρ = 1 4πr 2 e −2λκ rλ ′ . Hence since λ(r) ∼ 1/r for large r there is no boundary term at infinity in the integration by parts. Finally, using (5.37) it follows that Summing I j , j = 1, 2, 3, we obtain the asserted formula for D 2 E κ , cf. (5.24).
In the context of stellar perturbations the natural notion of dynamically accessible perturbations parallels Definition 4.14 for the Einstein-Vlasov system. This definition is motivated by (5.25), which states that for solutions of the linearized system,ρ ∈ R (A κ ). In other words, the space of linearly dynamically accessible perturbations is identified with the tangent space at ρ κ of the leaf of all density configurations of a fixed total mass M (ρ) = M (ρ κ ). Since R (A κ ) = (ker A ′ κ ) ⊥ and A ′ κ = e − 3 2 λκ n 1 2 κ ∂ r we have R (A κ ) = (ker ∂ r ) ⊥ and therefore A weaker notion of linear stability is spectral stability.
Definition 5.10. The steady state (ρ κ , λ κ , µ κ ) is spectrally stable if the spectrum of the linearized operator A ′ κLκ A κ is non-negative. Combining Lemmas 5.6 and 5.8, we can provide a variational criterion for spectral stability of steady states of the Einstein-Euler system, which in the physics literature goes back to Chandrasekhar [16].
Proposition 5.11 (Chandrasekhar stability criterion). A steady state (ρ κ , λ κ , µ κ ) of the Einstein-Euler system given by Proposition 3.2 is spectrally stable if and only if where the constraint set C κ = R(A κ ) is the set of all linearly dynamically accessible perturbations.
Proof. We use the second-order formulation (5.31) and note that the spectral stability is equivalent to the non-negativity of the quadratic form The claim now follows from (5.38).
Remark 5.12. By Lemma 5.8 spectral stability is equivalent to the positive semidefiniteness of the second variation of the energy along the manifold of linearly dynamically accessible perturbations.

The reduced Einstein-Euler operator
In analogy to Section 4.5 our goal is to show that the operatorL κ is bounded from above and below by a Schrödinger-type operator. Our main estimate is contained in Theorem 5.15. It is of independent interest, but it is also a crucial ingredient in our proof of Theorem 2. We first recall the definition (5.18) of the modified potentialμ and the fact that the relation (4.46) also holds in the present situation.
Definition 5.13 (Reduced operator for the Einstein-Euler system). By analogy to S κ the reduced operator Σ κ :Ḣ 1 r → (Ḣ 1 r ) ′ is defined by withḢ 1 r and ∆ κ defined as in Definition 4.23.
Remark 5.14. (a) We recall from Section 3.2 that the macroscopic quantities related to a steady state of the Einstein-Vlasov system constitute a steady state of the Einstein-Euler system. Indeed, in this situation S κ = Σ κ , since by (5.6) and (5.9), Theorem 5.15. Let Σ κ be the reduced operator given in Definition 5.13. For every µ ∈Ḣ 1 r and ρ = ρ µ := − e −µκ Ψκ µ it holds that For every ρ ∈ X 1 andμ =μ ρ as defined in (5.18) we have Proof. The proof is analogous to the corresponding result, Theorem 4.24, in the Einstein-Vlasov case. We first observe that the identities (4.68) and (4.69) remain valid in the present situation. Using these together with the definitions ofL κ and Σ κ implies that
Therefore, L κ * (and equivalently the reduced operator Σ κ * ) has a zero eigenvalue at κ * . This is a strong indication that the number of negative eigenvalues of L κ increases or decreases by 1 as κ crosses κ * . A detailed analysis of this subtle point will be given in [31].
(b) Since M κ is the ADM mass and R κ the area-radius of the 2-sphere t = 0 and r = R κ it follows that i κ has an invariant geometric meaning.
We now state the main result giving a detailed description of the linearized dynamics for the Einstein-Euler system.
Theorem 5.20. The operator J κ L κ generates a C 0 group e tJ κ L κ of bounded linear operators on X and there exists a decomposition where i κ is defined by (5.43).
(ii) The quadratic form (L κ ·, ·) X vanishes on E u,s , but is non-degenerate on E u ⊕ E s , and (iii) E c , E u , E s are invariant under e tJ κ L κ .
(iv) Let λ u = min{λ | λ ∈ σ(J κ L κ | E u )} > 0. Then there exist M > 0 such that where k 0 ≤ 1 + 2i κ . (5.51) The proof of Theorem 5.20 is similar to the proof of Theorem 4.28. We need several auxiliary lemmas. By an argument analogous to the proof of Lemma 4.18, one can show that the linear operator X 1 ∋ ρ → e −µκ Ψ −1 κμ ρ ∈ X 1 is compact and thus the only possible accumulation point of the spectrum of L κ and thereforeL κ is 1. Analogously to the proof of Proposition 4.30, we obtain the following result.
By the definition (5.30) of the operator L κ , the above lemma readily implies the following.
Proof of Theorem 5.20. To prove (i), we use the second order formulation (5.31). We note that When d dκ M κ d dκ (µ κ (R κ )) = 0, by (5.38) and the proof of Theorem 5.17, we have Thus there exist exactly n − (Σ κ ) − i κ negative eigenvalues of A ′ κL κ A κ in X 2 , each of which corresponds to a pair of stable and unstable eigenvalues for (5.31) and consequently for (5.29). This finishes the proof of part (i). The conclusions in (ii)-(iv) follow from [48, Thm. 2.2].

A negative energy direction for κ sufficiently large
In order to prove that the negative Morse index ofL κ is at least 1 we need to introduce an additional assumption on the equation of state. Assumption (P2). There exists a constant C > 0 such that P ′ (ρ) − 1 3 ≤ Cρ −1/2 for all ρ > 0. (5.52) We note that the examples P Φ and P NS which were discussed in Section 3.2 satisfy this assumption, cf. Propositions 3.3 and 3.5.
To analyze the sign of the quadratic form L κ ρ, ρ L 2 it will be convenient to rewrite in a slightly different form. Motivated by the Lagrangian coordinates point of view we consider perturbations of the form ρ = −(r 2 n κ ξ) ′ /r 2 and the associated metric perturbation λ. The second variation then can be written in the form A κ (ξ, ξ) := L κ ρ, ρ where ρ = −A κ ξ, and therefore ρ ∈ C κ is linearly dynamically accessible in the sense of Definition 5.9.
Proof. Some of the main ideas of this proof are analogous to the ones used in the proof of Theorem 4.3, but due to the local character of the operator L κ , the proof is actually simpler than the one of Theorem 4.3. Our starting point is the formula (5.53). We localize the perturbation ξ to the interval [r 1 κ , r 2 κ ] by setting ξ(r) = r a χ κ (r) for an a ∈ R to be specified later; the smooth cut-off function χ κ is as in (4.19). We split the quadratic form (5.53) into two parts: A κ = A κ,1 + A κ,2 , where A κ,1 = 4π e 2µκ+λκ Ψ κ 1 r 2 2 r 2+a n κ ′ r 2+a n κ χ κ + r 4+2a n 2 κ χ ′ κ χ ′ κ dr.
Remark 5.28. The proof of the full turning point principle for the Einstein-Vlasov system is considerably harder than in the case of the Einstein-Euler system since the space of linearly dynamically accessible perturbations R(B κ ) is of infinite codimension, by contrast to C κ , see (5.38).

A Auxiliary results for self-adjoint operators
In this appendix we collect some auxiliary results on symmetric and self-adjoint operators; to make the present paper self-contained, we include the proofs. First we show how estimates between two operators translate to corresponding information on their negative Morse index. Proof. Let V ⊂ D(A) be a subspace such that (Ax, x) X < 0 for all x ∈ V \ {0}, and define W := P (V ), which is a subspace of D(B). If x ∈ V such that P x = 0, then by (A.1) and the choice of V it follows that x = 0. Hence P : V → W is one-to-one and onto so that dim V = dim W . Let y ∈ W \ {0} so that y = P x for some x ∈ V \ {0}. Then by (A.1) and the choice of V it holds that (By, y) Y < 0. Hence n − (B) ≥ dim V , and the assertion follows.
Next we consider the spectral properties of a self-adjoint operator with negative Morse index.
Proposition A.2. Let A : X ⊃ D(A) → X be a self-adjoint operator on a Hilbert space X with scalar product (·, ·) and assume that n − (A) < ∞. Then the negative part of its spectrum, σ(A) ∩ ] − ∞, 0[, consists of at most finitely many eigenvalues which have finite multiplicities. If n − (A) > 0, then A has a negative eigenvalue of finite multiplicity.
For any x ∈ X + , S + x 2 = (L + x, x) = (Lx, x) ≥ δ x 2 , which implies that for all x ∈ X + , This lower bound of S + implies that T + := S + P + A is also closed with the dense domain D(T + ) = D(A) and thus T * + T + is self-adjoint. Then A * LA =A * P * + L + P + A − A * P * 1 L 1 P 1 A = (A * P * + S + )(S + P + A) − A * P * 1 L 1 P 1 A :=T * where B 1 = A * P * 1 L 1 P 1 A is bounded and symmetric. Therefore, A * LA is self-adjoint by the Kato-Rellich Theorem.