A numerical stability analysis for the Einstein-Vlasov system
Sebastian Günther, Jacob Körner, Timo Lebeda, Bastian Pötzl, Gerhard Rein, Christopher Straub, Jörg Weber
aa r X i v : . [ g r- q c ] S e p A numerical stability analysis for theEinstein-Vlasov system
Sebastian G¨unther ∗ , Jacob K¨orner † , Timo Lebeda ‡ , Bastian P¨otzl ∗ ,Gerhard Rein ∗ , Christopher Straub ∗ , J¨org Weber § September 18, 2020
Abstract
We investigate stability issues for steady states of the sphericallysymmetric Einstein-Vlasov system numerically in Schwarzschild, max-imal areal, and Eddington-Finkelstein coordinates. Across all coordi-nate systems we confirm the conjecture that the first binding energymaximum along a one-parameter family of steady states signals theonset of instability. Beyond this maximum perturbed solutions eithercollapse to a black hole, form heteroclinic orbits, or eventually fullydisperse. Contrary to earlier research, we find that a negative bindingenergy does not necessarily correspond to fully dispersing solutions.We also comment on the so-called turning point principle from theviewpoint of our numerical results. The physical reliability of the lat-ter is strengthened by obtaining consistent results in the three differentcoordinate systems and by the systematic use of dynamically accessibleperturbations.
We consider in the context of general relativity a large ensemble of masspoints which interact only through the gravitational field which they createcollectively. Such a self-gravitating collisionless gas is used in astrophysicsto model galaxies or globular clusters. Gravity is described by the Einsteinequations G αβ = 8 πT αβ , (1.1) ∗ Department of Mathematics, University of Bayreuth, Germany † Institute of Mathematics, Julius-Maximilians-Universit¨at W¨urzburg, Germany ‡ Department of Physics, University of Bayreuth, Germany § Centre for Mathematical Sciences, Lund University, Sweden G αβ is the Einstein tensor induced by the Lorentzian metric g αβ withsignature ( − + + +) on the smooth spacetime manifold M , and T αβ isthe 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 oflight and the gravitational constant are equal to 1. The evolution equationfor a collisionless gas is the collisionless Boltzmann or Vlasov equation sothat we obtain the Einstein-Vlasov system. We study this system under theassumption that the spacetime is spherically symmetric and asymptoticallyflat, but we first formulate it in general.The world line of a test particle on M obeys the geodesic equation˙ x α = p α , ˙ p α = − Γ αβγ p β p γ , where x α denote general coordinates on M , p α are the corresponding canon-ical momenta, Γ αβγ are the Christoffel symbols induced by the metric g αβ ,the dot indicates differentiation with respect to proper time along the worldline of the particle, and the Einstein summation convention is applied. Weassume that all the particles in the ensemble have the same rest mass, nor-malized to 1, and move forward in time, i.e., their number density f is anon-negative function supported on the mass shell P M := n g αβ p α p β = − , p α future pointing o , a submanifold of the tangent bundle T M of the spacetime manifold M whichis invariant under the geodesic flow. Letting Latin indices range from 1 to 3we choose coordinates ( t, x a ) such that on the mass shell P M the variable p becomes a function of the remaining variables ( t, x a , p b ); t should bethought of as a time-like variable. Since the particles in the ensemble movelike test particles, their number density f = f ( t, x a , p b ) is constant along thegeodesics and hence satisfies the Vlasov equation ∂ t f + p a p ∂ x a f − p Γ aβγ p β p γ ∂ p a f = 0 . (1.2)The energy-momentum tensor is given by T αβ = Z p α p β f | g | / dp dp dp − p , (1.3)where | g | denotes the modulus of the determinant of the metric, and indicesare 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 to2odel isolated systems and therefore require that the spacetime is asymp-totically flat. In order to simplify the system we only consider sphericallysymmetric solutions. In the next section we formulate the Einstein-Vlasovsystem in coordinates adapted to this symmetry. For background on theEinstein-Vlasov system we refer to [1] and the references there.The Einstein-Vlasov system possesses a plethora of steady state solu-tions. For a stationary metric the Killing vector ∂/∂t gives rise to the quan-tity E = − g ( ∂/∂t, p α ) which represents the particle energy and is constantalong geodesics. Hence the ansatz f ( x a , p b ) = φ ( E ) (1.4)satisfies the stationary Vlasov equation and reduces the system to the fieldequations. In Section 2.5 we recall how any such “microscopic equation ofstate φ ” gives rise to a one-parameter family of steady states, where theparameter can be identified with the central redshift of the configurationand is therefore a measure of how relativistic it is; the steady states actuallyconsidered below can also depend on the angular momentum of the particles.A natural question is which of these steady states are stable or unstable. Forthe Vlasov-Poisson system, which is the non-relativistic limit of the Einstein-Vlasov system, such steady states essentially are stable if the microscopicequation of state φ is a decreasing function of the particle energy, cf. [23]and the references there. For the relativistic case of the Einstein-Vlasovsystem the situation is quite different. At least on the linearized level ithas been shown in [10, 11, 12] that such steady states are stable if theircentral redshift is sufficiently small, but for the same microscopic equationof state they become unstable if their central redshift is large. The questionwhether sufficiently relativistic matter distributions become unstable playedan important role in the discovery and subsequent discussion of quasars, cf.[7, 15, 30]. But unstable steady states of the Einstein-Vlasov system arealso important for conceptual reasons, since they can possibly explain theso-called type I behavior in critical collapse observed for the Einstein-Vlasovsystem, the latter being related to the cosmic censorship hypothesis, cf.[3, 18, 25]. All this motivates the present investigation where we analyze thetransition from stability to instability along one-parameter families of steadystates of the Einstein-Vlasov system by numerical means. In particular, weinvestigate where this transition takes place and what happens to weaklyperturbed steady states which lie in the unstable regime.Concerning the former question there are various possibilities. One canfor example consider the so-called binding energy as a function of the central3edshift, cf. Figure 1. It has been conjectured [29, 30] that the transitionfrom stability to instability happens at the first (local) maximum of thiscurve. For all microscopic equations of state which we consider we confirmthis conjecture.Alternatively, one can plot for a fixed microscopic equation of state andeach value of the central redshift the ADM mass and radius of the supportof the corresponding steady state. This results in a so-called mass-radiuscurve, cf. Figure 2. It would be conceivable that the stability propertieschange at the turning points of this curve since a precise version of this so-called turning point principle has recently been proven both for the Euler-Poisson and the Einstein-Euler system, cf. [9, 16]. In these models, matteris described as an ideal, compressible fluid. The macroscopic quantitiesinduced by an isotropic steady state of the Einstein-Vlasov system of theform (1.4) yield a steady state of the Einstein-Euler system with a suitable,induced macroscopic equation of state, cf. [12]. In particular, the mass-radius curves are then the same for both systems. However, we clearlydisprove this turning point principle for the Einstein-Vlasov system.Another issue is to understand the behavior of solutions which arelaunched by small perturbations of a stable or an unstable steady state.In the former case we find that the system starts to oscillate, i.e., to expandand contract in a seemingly time-periodic fashion. This behavior was ob-served in [20] for the Vlasov-Poisson system and in [3] for shell-like solutionsof the Einstein-Vlasov system; the code employed in [3] was not able toproperly handle steady states which have matter at the center instead of avacuum region. The behavior of an unstable steady state after perturbationis more interesting. The solution either collapses and forms a black hole,or it seems to follow a heteroclinic orbit to a different, stable steady stateabout which (the bulk of) it starts to oscillate. The terminology “hetero-clinic orbit” may not be quite appropriate here, but it captures the observedbehavior. In [28] a similar observation is claimed without further comment.For steady states with a large central redshift the perturbed state may alsodisperse instead of following a heteroclinic orbit as explained above.In all the numerical simulations we used dynamically accessible pertur-bations which in particular preserve all the so-called Casimir functionals(3.1) of the system; perturbing the steady state by some external forceresults in such dynamically accessible states. All the simulations were per-formed in three different coordinate systems, namely Schwarzschild, maxi-mal areal, and Eddington-Finkelstein coordinates, and our observations werecompletely consistent across these. This is a priori not obvious. In a stabilityanalysis for the Einstein-Vlasov system one necessarily must compare func-4ions like metric components or mass-energy densities which are defined ontwo different spacetimes, the stationary one and the perturbed one. Thereis no canonical way of identifying points on these two spacetimes so thatone could compare the values of certain functions at those identified points.What we do is to simply identify points which have the same coordinates inthe coordinate system at hand, and it is therefore not a priori clear that thestability findings in different coordinate systems must be consistent. Thispoint has also been made in the astrophysics literature, cf. [15].The paper proceeds as follows. In the next section we formulate theEinstein-Vlasov system in the coordinate systems mentioned above, com-pare these coordinate systems, and recall how steady states of the Einstein-Vlasov system are obtained. In Section 3 we explain the numerical methodwhich we employ, which is a particle-in-cell scheme and lends itself well toparallelization. The numerical results are presented and discussed in Sec-tion 4, and in the last section we comment on their precision and reliability. In this section we formulate the Einstein-Vlasov system in Schwarzschild,maximal areal, and Eddington-Finkelstein coordinates and compare variousproperties of these coordinate systems. We also recall how steady statesolutions can be obtained.
In Schwarzschild coordinates the metric reads ds = − e µ ( t,r ) dt + e λ ( t,r ) dr + r (cid:0) dθ + sin θdϕ (cid:1) (2.1)with ( t, r, θ, ϕ ) ∈ R × [0 , ∞ [ × [0 , π ] × [0 , π ] and metric coefficients µ = µ ( t, r )and λ = λ ( t, r ). Here, t corresponds to the proper time of an observerlocated at spatial infinity and r denotes the areal radius. For numerical andanalytical reasons it is convenient to introduce Cartesian coordinates x = ( x , x , x ) = r (sin θ cos ϕ, sin θ sin ϕ, cos θ ) ∈ R and corresponding non-canonical momentum variables v i = p i + (cid:16) e λ − (cid:17) x · pr x i r . r →∞ λ ( t, r ) = lim r →∞ µ ( t, r ) = λ ( t,
0) = 0 , t ∈ R . (2.2)Inserting the metric into the Einstein equations yields the following fieldequations: e − λ (cid:0) rλ ′ − (cid:1) + 1 = 8 πr ρ, (2.3) e − λ (cid:0) rµ ′ + 1 (cid:1) − πr p, (2.4)˙ λ = − πre µ + λ j. (2.5)These are the 00, 11, and 01 components of the general equation (1.1).Equation (2.5) is not independent, but follows from (2.3) and (2.4) togetherwith the Vlasov equation. However, it is useful for the numerics. The alsonon-trivial 22 and 33 components of (1.1) follow as well, but they are notused in the numerics. In the above, ˙ and ′ denote the derivative with respectto t or r respectively. The Vlasov equation takes the form ∂ t f + e µ − λ vε · ∂ x f − (cid:16) ˙ λ x · vr + µ ′ e µ − λ ε (cid:17) xr · ∂ v f = 0 , (2.6)where we introduce ε = p | v | = r w + Lr . Here | v | denotes the Euclidean length and x · v the Euclidean scalar product.The variables w = x · vr and L = | x × v | can be thought of as the momentum inthe radial direction and the square of the angular momentum respectively.We assume f to be spherically symmetric, i.e., f ( t, x, v ) = f ( t, Ax, Av )for A ∈ SO(3), and under abuse of notation we may write f ( t, x, v ) = f ( t, r, w, L ). The source terms in the field equations are defined by ρ ( t, r ) = πr Z ∞ Z ∞−∞ εf ( t, r, w, L ) dw dL, (2.7) p ( t, r ) = πr Z ∞ Z ∞−∞ w ε f ( t, r, w, L ) dw dL, (2.8) j ( t, r ) = πr Z ∞ Z ∞−∞ wf ( t, r, w, L ) dw dL. (2.9)6ere, ρ is the energy density, p the radial pressure, and j the particlecurrent. Equations (2.2)–(2.9) constitute the Einstein-Vlasov system inSchwarzschild coordinates. Unless stated otherwise, we employ the abovenotation for the other coordinate systems as well.For later use, we derive some formulas for the metric coefficients. To thisend, we introduce the Hawking mass defined as m ( t, r ) = 4 π Z r ρ ( t, s ) s ds. (2.10)Integrating (2.3) yields e − λ = 1 − mr ; (2.11)the right hand side of this equation remains positive as long as the solutionto (2.2)–(2.9) exists. Solving the field equation (2.4) for µ ′ and using (2.11),we obtain µ ′ = e λ (cid:16) πrp + mr (cid:17) . (2.12) In maximal areal coordinates the line element can be written as ds = ( − α + a β ) dt + 2 a βdtdr + a dr + r ( dθ + sin θdϕ )with positive metric coefficients a and α . We note that in the present case t is not introduced as the proper time of any physical observer, but is fixed byimposing the maximal gauge condition, i.e., each hypersurface of constant t has vanishing mean curvature. The non-canonical momentum variablesintroduced in the Schwarzschild case translate to v i = p i + (cid:18) a − (cid:19) x · pr x i r . In analogy to (2.2) the metric coefficients satisfy the boundary conditions a ( t,
0) = lim r →∞ a ( t, r ) = lim r →∞ α ( t, r ) = 1 , β ( t,
0) = 0 . (2.13)7e obtain the field equations κ = βαr , (2.14) a ′ = 4 πrρa + 32 rκ a + a r (cid:0) − a (cid:1) , (2.15) κ ′ = − κr − πaj, (2.16) α ′′ = α ′ (cid:18) a ′ a − r (cid:19) + 6 αa κ + 4 παa ( S + ρ ) , (2.17)which are coupled to the Vlasov equation ∂ t f + h αa vε − β xr i · ∂ x f + (cid:20) − ε α ′ a xr + ακ (cid:16) v − x · vr xr (cid:17)(cid:21) · ∂ v f = 0 (2.18)via the source terms (2.7), (2.9), and S ( t, r ) = πr Z ∞ Z ∞−∞ ε − ε f ( t, r, w, L ) dw dL, (2.19)the trace of the spatial part of the energy-momentum tensor.As in Schwarzschild coordinates we introduce the Hawking mass, whichin maximal areal coordinates becomes m = r (cid:18) − a + r κ (cid:19) . In order to solve the field equations numerically, it is useful to consider thequantity η = r (cid:18) − a (cid:19) , which immediately implies η ( t, r ) = Z r (cid:18) πρ ( t, s ) + 32 κ ( t, s ) (cid:19) s ds. (2.20)Furthermore, the field equation (2.16) yields the implicit formula κ ( t, r ) = − πr Z r a ( t, s ) j ( t, s ) s ds, while α ′ ( t, r ) = a ( t, r ) r Z r (cid:0) πaα ( ρ + S ) + 6 aακ (cid:1) s ds holds because of the second order equation (2.17). Note that κ ( t, r ) ∼ r − for large r provided that the matter is compactly supported.8 .3 Eddington-Finkelstein coordinates In Eddington-Finkelstein coordinates the metric takes the form ds = − a ( t, r ) b ( t, r ) dt + 2 b ( t, r ) dtdr + r ( dθ + sin θ dϕ ) . Similar to (2.2) and (2.13) the metric coefficients a and b satisfy the bound-ary conditions a ( t,
0) = lim r →∞ a ( t, r ) = lim r →∞ b ( t, r ) = 1 . (2.21)Notice that the metric coefficient a here is not the same as the coeffi-cient a appearing in maximal areal coordinates, but we nevertheless usethis quite common, albeit equivocal notation. As opposed to Schwarzschildand maximal areal coordinates we use the canonical momentum coordinates( p , p , p , p ). The angular momentum is given by L = ( p ) + 1sin θ ( p ) . The particle density f can be written as a function of ( t, r, p , L ) and theVlasov equation reads ∂ t f + b (cid:16) a − L/r ( p ) (cid:17) ∂ r f + 12 bLr p − ∂ r ( ab ) p − ∂ r b L/r p ! ∂ p f = 0 . Here, the Hawking mass m is given by m = r − a ) . The metric coefficients a and b as well as the Hawking mass can be computeddirectly from f via b ( t, r ) = exp (cid:18) − π Z ∞ r ηT ( t, η ) dη (cid:19) ,m ( t, r ) = 2 πb ( t, r ) Z r η ( T + S )( t, η ) b ( t, η ) dη,a ( t, r ) = 1 − m ( t, r ) r , T ( t, r ) = πr Z ∞ Z ∞ p f ( t, r, p , L ) dL dp ,S ( t, r ) = πr Z ∞ Z ∞ Lr p f ( t, r, p , L ) dL dp . We have to keep in mind that the physical interpretation of the timelikevariable t differs across the three coordinate systems, and we should men-tion that in Eddington-Finkelstein coordinates, what we called t is usuallydenoted as v . Before we investigate the system numerically, we briefly discuss and compareselected properties of the coordinate systems. Two conserved quantities ofthe system are the total number of particles or total rest mass and the ADMmass. Across all coordinate systems the latter is given by M = lim r →∞ m ( t, r ) . (2.22)For compactly supported matter this equals the Hawking mass evaluated atthe outer boundary of the radial support.The total number of particles is computed differently across the coordi-nate systems. In the Schwarzschild case we have N = 4 π Z ∞ Z ∞−∞ Z ∞ e λ ( t,r ) f ( t, r, w, L ) dr dw dL, (2.23)while in maximal areal coordinates, N = 4 π Z ∞ Z ∞−∞ Z ∞ a ( t, r ) f ( t, r, w, L ) dr dw dL, (2.24)and in Eddington-Finkelstein coordinates, N = 4 π Z ∞ Z ∞ Z ∞ f ( t, r, p , L ) dr dp dL. (2.25)The weights e λ and a appear because the characteristic flow is not measurepreserving. This is due to the use of non-canonical momentum variables. Infact, if by D we denote differentiation along a characteristic of the Vlasovequation, then in the Schwarzschild case D ( f dx dv ) = − (cid:16) ˙ λ + λ ′ e µ − λ wε (cid:17) f dx dv D ( f dx dv ) = − (cid:18) αa ′ a wε + β ′ + 2 βr (cid:19) f dx dv. In Eddington-Finkelstein coordinates canonical momentum variables areused which implies D ( f dx dv ) = 0.Another important difference between the three coordinate systems isthe existence of a criterion for the formation of trapped surfaces. Thereexists no such criterion in Schwarzschild coordinates since these coordinatescannot cover an open region which contains a trapped surface. In maximalareal coordinates a trapped surface is present when1 a ( t, r ) − rκ ( t, r ) < . (2.26)In this case, the expansion of both outgoing and ingoing null geodesics isnegative at the time t on the sphere of radius r . This signals the developmentof a spacetime singularity, cf. [19]. In Eddington-Finkelstein coordinates thecondition a ( t, r ) < Despite the differences of the above coordinate systems, there exists anexplicit coordinate transformation which maps stationary solutions in onecoordinate system into stationary solutions in the other two. For this reasonwe discuss steady states in Schwarzschild coordinates first.11 simple calculation shows that for a time-independent metric the par-ticle energy E = εe µ for Schwarzschild, εα for maximal areal, b (cid:16) ap + L/r p (cid:17) for Eddington-Finkelsteinis conserved along characteristics of the Vlasov equation, and due to spher-ical symmetry the same is true for the square of the angular momentum L .Thus, every ansatz of the form f = φ ( E, L ) solves the Vlasov equation inthe time-independent metric and reduces the system to the field equationsfor that metric, where the source terms now depend on the latter throughthe given ansatz. In the present paper we consider two different types ofansatz functions, namely the polytropic ansatz f = (cid:18) − EE (cid:19) k + ( L − L ) l + (2.28)and the King model f = (cid:16) e − EE − (cid:17) + . (2.29)Here z + denotes the positive part of some number z . The constants E > L ≥ k ≥ l ≥ k < l + , but also consider the borderline case of k = and l = 0. The cut-off energy is necessary to ensure finite extensionand finite ADM mass of the corresponding stationary solution. Prescribing L > L = 0 implies that the support ofthe steady state contains the origin. Note that the factor ( L − L ) l + can alsobe multiplied to the King ansatz, with the analogous effect on the supportof the resulting steady state.Even though multiple radially separated shells may arise from a poly-tropic ansatz with L >
0, we only consider the innermost shell as our steadystate, cf. [4]. Formally, these steady states depend not only on E and L butalso on the choice of the shell. However, we still refer to them as polytropes.In order to solve for the metric coefficients, it turns out to be moreconvenient to consider y = ln E − µ . For an ansatz as above we can write12he spatial mass density and pressure as functions of y , i.e., ρ ( r ) = g ( y ( r )) , p ( r ) = h ( y ( r )) . For the definition of g and h and further details we refer to [20]. The sta-tionary Einstein-Vlasov system is then reduced to the differential equation y ′ = − − πr R r g ( y ( s )) s ds (cid:18) πr Z r g ( y ( s )) s ds + 4 πrh ( y ( r )) (cid:19) . (2.30)This equation is (2.12) with (2.11) substituted in. For every choice of y (0) = y >
0, [20] guarantees a unique solution of (2.30) with finite ADMmass and compact support. Hence one given ansatz of the form (2.28) or(2.29) yields a one-parameter family of such steady states. The parameteris closely connected to the central redshift z c , measuring the redshift of aphoton which is emitted at the center r = 0 and received at the boundaryof the steady state: z c = e y − . (2.31)This is not the standard definition of the central redshift where the photonis received at infinity, but our definition is more suitable here. If y is asolution to (2.30) and the cut-off energy is defined as E = lim r →∞ e y ( r ) ,the metric coefficient µ = ln E − y satisfies the proper boundary conditionat infinity, and λ is given by λ ( r ) = −
12 ln (cid:18) − πr Z r ρ ( s ) s ds (cid:19) . We now discuss how to obtain the stationary solutions in the other two co-ordinate systems. In maximal areal coordinates stationary solutions satisfy κ = β = j = 0. The metric is then equivalent to the Schwarzschild metricby simply setting µ = ln α and λ = ln a , cf. [3]. In the case of Eddington-Finkelstein coordinates the metric coefficients a and b can be obtained bythe change of variables t t + Z r e λ ( s ) − µ ( s ) ds, and the metric coefficients are related via a = e − λ , b = e λ + µ . y >
0, we computethe corresponding solution y to the equation (2.30). We use the explicitEuler method enhanced with a leap frog scheme to solve (2.30) and applySimpson’s rule to calculate the integrals appearing in g , h , and (2.30). Notethat the whole steady state computation has to be done only once, whichallows us to use very high accuracy for this part. The distribution function f of the steady state is then given via (2.28) or (2.29). The algorithms used to investigate the stability of steady states of the abovesystem are based on the particle-in-cell scheme. This scheme has also beenused in [3, 21, 25]. For the spherically symmetric Vlasov-Poisson system itsconvergence has been shown in [27], the analogous result for the Einstein-Vlasov system in Schwarzschild coordinates is shown in [26]. In a particle-in-cell scheme the support of the distribution function f is initially splitinto distinct cells. Into each cell a numerical particle is placed to representthe contribution to f of this cell, and these particles are then propagatedaccording to the Einstein-Vlasov system. Most of the following steps differfor the three coordinate systems under investigation. We focus on the algo-rithm for Schwarzschild coordinates, but always highlight difficulties arisingin the other two coordinate systems and how to overcome them.To initialize the numerical particles we use variables adapted to thespherical symmetry of the steady state f , i.e., we write f = f ( r, w, L ) = f ( r, u, ψ ) , where we use the additional variables u ≥ ψ ∈ [0 , π ] given by u = w + Lr , w = u cos( ψ ) . Assuming that f vanishes outside of the set [ R − , R + ] × [ U − , U + ] × [Ψ − , Ψ + ],we prescribe a radial step length ∆ r > N u and N ψ todefine the step lengths∆ u = U + − U − N u , ∆ ψ = Ψ + − Ψ − N ψ r i = (cid:18) i − (cid:19) ∆ r, u j = (cid:18) j − (cid:19) ∆ u, ψ k = (cid:18) k − (cid:19) ∆ ψ. At each point ( r i , u j , ψ k ) we generate a numerical particle carrying theweight f i,j,k = f ( r i , u j , ψ k ) 4 πr i ∆ r πu j ∆ u sin( ψ k )∆ ψ, where f ( r i , u j , ψ k ) is calculated using the steady state from above. Weuse ( r, u, ψ )-variables in Schwarzschild and maximal areal coordinates togenerate the numerical particles. In the Eddington-Finkelstein case a similarinitialization scheme based on ( r, p , L )-variables is used. For the followingsteps it is convenient for Schwarzschild and maximal areal coordinates towrite the particle positions in ( r, w, L )-variables, since L is conserved alongcharacteristics.We then compute the matter quantities ρ , p , and j on the fixed radialgrid given by r j = j ∆ r by integrating f according to (2.7), (2.8), and (2.9).This is implemented by adding up f i,j,k with the appropriate weight givenby (2.7)-(2.9) and by a linear interpolation of f in the radial direction; notethat f i,j,k contains the phase space volume element.Next, we compute the Hawking mass m on the fixed radial grid by using(2.10) and a quadratic interpolation, taking into account the possible orderof the integrand near the origin. Afterwards, the metric quantity µ is calcu-lated by applying (2.11), (2.12), and the boundary condition µ ( R ) = λ ( R ),where R denotes the outer boundary of the radial support.In maximal areal coordinates we instead solve for η and κ simultaneouslyby employing a fourth-order Runge-Kutta method in order to guarantee suf-ficient numerical precision and to appropriately handle the implicit structureof the field equations. Determining α requires more effort due to the elliptic-ity of equation (2.17) and the boundary condition given at spatial infinity.We choose a sufficiently large grid, approximate the various derivatives andarrive at a tridiagonal system which can be solved explicitly. For moredetails consider the scheme used in [3].The time step is then performed by propagating the numerical particlesaccording to the characteristic system corresponding to the Vlasov equation.In order to avoid numerical errors at the spatial origin caused particularlyby particles coming close to the origin, it is advantageous to use Cartesiancoordinates for the propagation. Here, all functions involved are interpolatedaccording to their order, especially near the origin. However, we do not use15artesian coordinates in the case of Eddington-Finkelstein coordinates sincethe metric written in these coordinates is not continuous at the origin, cf.[6]. The new position of each particle is then computed by a proper steppingmethod with a prescribed time step size ∆ t >
0. Note that we also have toupdate the weight f i,j,k during each time step in the case of Schwarzschildand maximal areal coordinates, since the characteristic flow of the Vlasovequation system is not measure preserving in these coordinates.At this point the particle coordinates and weights are known at the time∆ t and we can repeat the iteration until we reach some final time T whichis prescribed from the beginning.So far we have numerically evolved a steady state f itself, but in orderto analyze its stability properties we have to perturb it. In [3] this is doneby taking A f for some parameter A ≈ C ( f ( t )) = Z Z e λ ( t,x ) χ ( f ( t, x, v )) dx dv (3.1)where χ ∈ C ( R ) with χ (0) = 0; a physically viable perturbation for exampleby some external force should preserve these. We provide a perturbationprocedure preserving these invariants. An analogous type of perturbationhas been used in [21] for the Vlasov-Poisson system, but as far as we knowthis is the first implementation of dynamically accessible perturbations tosteady states of the Einstein-Vlasov system. During an initial time interval[0 , T pert ] we propagate the numerical particles according to the modifiedcharacteristic system˙ x = e µ − λ vε , ˙ v = − (cid:16) ˙ λ x · vr + µ ′ e µ − λ ε (cid:17) xr + γ xr (3.2)in the Schwarzschild case for prescribed γ ≈
0. Compared to the origi-nal characteristic system, we add the term γ xr to the right hand side of theequation for ˙ v . Since Casimir functionals are preserved along solutions of theVlasov equation and the above perturbation does not contribute to the di-vergence of the right hand side of the characteristic system, this perturbationis indeed dynamically accessible. In maximal areal coordinates we also addthe divergence free term (0 , γ xr ) to the right hand side of the respective char-acteristic system during an initial time interval. In Eddington-Finkelsteincoordinates, where E = E ( r, p ) is the Hamiltonian governing the motion ofthe particles, we add γrp as a perturbation of the Hamiltonian.16espite the fact that v does not denote the canonical momentum, allthese perturbations can be interpreted as particles being accelerated eitherradially outwards if γ > γ < t = T pert as the perturbed state. Todetermine the intensity of the perturbation we consider e µ ( t, , α ( t, b ( t,
0) respectively. We prescribe a small number ǫ pert and then choose γ such that the relative error of e µ ( t, between t = 0 and t = T pert is closeto ǫ pert . Obviously, numerical and evolutionary effects contribute to theperturbation of the steady state as well. To this end, we choose ǫ pert suchthat it dominates the relative error of the above quantities resulting fromthe numerical initialization. In addition, we choose T pert sufficiently smallsuch that evolutionary effects are negligible.We emphasize that we have tested several other perturbations. On theone hand, we used rather simple perturbations where we scale the steadystate by some factor or shift it slightly in one direction. On the other hand,we employed further dynamically accessible perturbations where we prop-agate the numerical particles using modified metric quantities during aninitial time interval. For every kind of perturbation there seem to exist pre-cisely two distinct characteristic behaviors: pushing the steady state eithertowards collapse or towards dispersion. However, all the effects described inthe next section seem to only depend on the direction of the perturbationbut not on its specific type.Let us end this section with a remark concerning the numerical imple-mentation. In order to obtain reliable physical results, we have to workwith tens of millions of numerical particles and sufficiently small time steps.For these simulations to run within a reasonable time-frame, the programshave to be parallelized, which however fits very well with the particle-in-cellscheme. We refer to [17] for a detailed discussion. In the following we investigate the stability behavior of the one-parameterfamilies of steady states which we obtain by choosing the King model or thepolytropic ansatz with fixed parameters k , l , and L . Due to its physicalrelevance we mainly present the results for the King model. The parameterthat determines the steady state is the central value y where larger valuesof y generally indicate more relativistic scenarios, cf. (2.31). For a detailedoverview of properties of the steady states we refer to [4].17 physically meaningful quantity of a steady state is its binding energy E b = N − MN , where N is the number of particles and M the ADM mass. In Figure 1the binding energy is plotted against y for different models, namely forthe King model and polytropes with ( k, l, L ) = (0 . , , . , , . . , . , . − . − . − . − . .
05 0 0 . . E b y King model k = 0 . l = 0, L = 0 k = 0 . l = 0, L = 0 . k = 0 . l = 0 . L = 0 . Figure 1: Binding energy for different models. L = 0 .
001 may seem very small, it still suffices to change the correspondingone-parameter family of steady states significantly. Table 1 gives the values y max and y zero where the binding energy attains its first maximum or its firstzero respectively for the different models used in Figure 1. Another way tovisualize a one-parameter family of steady states originating from a givenansatz function is the so-called mass-radius diagram where for each valueof y > y max y zero King 0.334 0.784 k = 0 . l = 0, L = 0 0.298 0.693 k = 0 . l = 0, L = 0 .
001 0.267 0.583 k = 0 . l = 0 . L = 0 .
001 0.265 0.588Table 1: First binding energy maximizer and zero.The spiral structure is a general feature of these mass-radius diagrams forthe Einstein-Vlasov system, cf. [4]. According to the so-called turning pointprinciple the steady state should pass from being stable to being unstableas it crosses the first maximum point along the curve, but we find that thisis not true. . . . . . . . . . .
13 1 2 3 4 5 6first maximum of E b M R Figure 2: Mass-radius spiral for the King model; as y increases the corre-sponding ( R, M ) moves into the spiral.As mentioned at the end of Section 3, perturbations can be categorizedas either promoting collapse or promoting dispersion. For example, the dy-namically accessible perturbation of f with parameter γ described in (3.2)19ushes the steady state towards collapse if γ <
0, and towards dispersion if γ > For small values of y we find that the steady states are stable with respectto every reasonable perturbation. However, perturbations do not leave thesteady state unchanged but make it oscillate in a pulsating manner. Sim-ilar oscillations have been numerically observed as perturbations of stablesteady states of the Vlasov-Poisson system, cf. [21]. Figure 3 illustrates thisbehavior on the level of the mass distribution. Alternatively, the quantity r t Figure 3: The oscillation of the weighted mass density 4 πr ρ for the Kingmodel in Schwarzschild coordinates for y = 0 . γ >
0. The color represents the valueof 4 πr ρ , increasing from black to yellow. e µ ( t, , α ( t, b ( t, γ when perturbing a sta-ble steady state consider Figure 4. On the one hand, for a perturbation . . . . . .
885 0 50 100 150 200 e µ ( t , ) t γ < γ = 0 γ > Figure 4: The King model in Schwarzschild coordinates for y = 0 . γ < γ = 0,and γ > γ <
0) the quantity e µ ( t, initially decreases. Onthe other hand, e µ ( t, initially increases for a perturbation that promotesdispersion ( γ > y . At some threshold value the stabilitybehavior changes. 21 .2 The first binding energy maximum—onset of instability One main goal of our numerical investigation of the time dependent sys-tem is to detect when instability of steady states first occurs along a one-parameter family. We find convincing evidence that the so-called “bindingenergy maximum hypothesis” holds which has also been analyzed and con-firmed in [3, 13, 14, 28] for different models and perturbations. It was firstproposed by Zel’dovich et al. in [29, 30]. The hypothesis states that the firstbinding energy maximum along a steady state sequence signals the onset ofinstability.The starting point of our investigation is [3] where the binding energymaximum hypothesis has been verified numerically in maximal areal coor-dinates for the case of polytropic shell steady states, i.e., L > L = 0, and in particular to isotropic models, i.e., L = 0 and l = 0, which allows particles to pass through the origin. Our findings clearlysupport the binding energy maximum hypothesis in these cases.We present the results of simulations for four cases across the threecoordinates systems: the isotropic case for k = 0 . , , . y max , with sufficient accuracy. We then consider steady stateswith values of y close to y max . As mentioned in Section 4.1 we find that thesteady states are stable for both types of perturbations for y < y max acrossall models and coordinate systems. For y > y max the stability behaviorchanges. In fact, in all coordinate systems and models we observe that aperturbation promoting collapse leads to the actual collapse of the steadystate to a black hole, which is illustrated in Figure 5. Notice that for y
314 and y = 0 .
324 in Figure 6, i.e., y < y max , the oscillation startsimmediately and the first local maximum is at about t = 15. On the otherhand, for y = 0 .
344 and y = 0 .
354 there exists an initial phase over thecourse of which the increase of e µ ( t, seems to occur independently from theoscillation. Therefore, the first local maximum appears much later at about t = 50. The general behavior for y > y max will be analyzed more closely inSection 4.4.We now come back to the question whether the turning point principle isvalid for the Einstein-Vlasov system. As explained above, our numerics showthat in a one-parameter family steady states are stable as long as y < y max ,and they become unstable as soon as y > y max ; the change from stabilityto instability occurs at the first maximum of the binding energy. In themass-radius spiral in Figure 2 we have marked the point which corresponds23 . . . .
68 0 40 80 120 160 e µ ( t , ) t y = 0 . y = 0 . y = 0 . y = 0 .
354 0 . . . .
75 0 40 80 120 160 e µ ( t , ) t y = 0 . y = 0 . y = 0 . y = 0 . Figure 6: The King model in Schwarzschild coordinates for y close to y max perturbed by a dynamically accessible perturbation with direction γ > y max , i.e., to the first maximum of the binding energy. But this point isclearly to the left of the first turning point of the mass-radius spiral, i.e.,steady states for the Einstein-Vlasov system remain stable even after theturning point principle would predict their instability and after the verysame steady states have become unstable as steady states of the Einstein-Euler system. We have observed this discrepancy or failure of the turningpoint principle in all cases which we investigated numerically, and it seemsa very interesting mathematical problem to properly understand this issue. One of the main advantages of using three coordinate systems is the possibil-ity to observe the properties of collapsing matter from different perspectives.In Schwarzschild coordinates we have no analytical criterion for the detec-tion of a trapped surface. However, the ratio m ( t,r ) r approaching 1 at someradius signals the development of a black hole. On the contrary, in max-imal areal and Eddington-Finkelstein coordinates trapped surfaces can bedetected analytically, cf. (2.26) and (2.27). Note that rapidly increasingmetric coefficients eventually cause the program to break down in the eventof a collapse.First, let us explain our observations when a collapse to a black hole oc-curs. As far as the numerical particles within the simulation are concerned,we notice that the matter focuses towards the origin and the matter density ρ increases at the center. In fact, in our simulations no particles remain unaf-fected by the collapse, and all particles are eventually and irreversibly suckedtowards the interior of the trapped surface. This eventually leads to the ra-24ius R of the support of the solution tending to the Schwarzschild radius2 M . This general behavior should be compared to the analytical findings in[5]. The observations on the particle level can be linked to the characteristicbehavior of metric coefficients. In Schwarzschild and maximal areal coor-dinates the collapse is reflected in the behavior of the functions e µ and α ,respectively. These quantities rapidly decrease to zero in close proximity ofthe origin while still satisfying the boundary conditions at infinity, cf. (2.2)and (2.13). In the literature, this phenomenon is commonly known as the“collapse of the lapse” since e µ and α determine the amount of proper timewhich elapses from one hypersurface of constant t to the next. The collapseof the lapse is depicted in Figure 7 where we can see an exponential decay of α for late times and r close to zero. In Eddington-Finkelstein coordinates,the metric coefficient b shows a similar behavior. Furthermore, the criteria − − − − − − − .
05 0 . .
15 0 . t = 0 t = 8 . t = 8 . t = 8 . t = 8 . t = 8 . t = 9 . t = 9 . t = 9 . t = 9 . t = 9 . α ( t , r ) r Figure 7: Evolution of α ( t, r ) for the King model in maximal areal coordi-nates for y = 0 . γ <
0. Notice the logarithmic scale.for the formation of a trapped surface mentioned above are satisfied fromsome moment in time onwards. Figure 8 shows the evolution of the metriccoefficient a in Eddington-Finkelstein coordinates where a < r . This supports the conjecture that the weakcensorship hypothesis holds for the Einstein-Vlasov system in our setting.For an overview of this topic we refer to [1]. Independently of the choice − − . .
51 0 0 .
05 0 . .
15 0 . t = 0 t = 4 . t = 4 . t = 4 . t = 4 . a ( t , r ) r Figure 8: Evolution of a ( t, r ) for the King model in Eddington-Finkelsteincoordinates for y = 0 . γ < y is enlarged. This can beexplained by the fact that larger values of y correspond to more relativisticsteady states. However, one should recall that the time coordinate t has adifferent meaning across the three coordinate systems. When perturbing an unstable steady state with a perturbation promotingdispersion, i.e., y > y max and γ >
0, the resulting solution disperses at first,which means that the matter distributes more evenly in space and spacetimebecomes flatter. This corresponds to an initial increase of e µ ( t, , α ( t, ( t, y is not too large, we can clearly see thatthe solution starts reimploding at some point in time, which means that e µ ( t, decreases again, giving rise to an oscillating behavior. However, thisoscillation occurs at values of e µ ( t, larger than e µ (0 , . We call this differencethe initial elevation. As indicated in Section 4.2, the initial elevation can beobserved for all y > y max . These effects occur across all classes of steadystates, types of perturbations promoting dispersion, and coordinate systemsin a similar way and are illustrated in Figure 9 in the case of the King modelin Schwarzschild coordinates for various y > y max . . . . . . . .
91 0 50 100 150 200 e µ ( t , ) t y = 0 . y = 0 . y = 0 . y = 0 . y = 0 . y = 0 . Figure 9: The King model in Schwarzschild coordinates for y > y max per-turbed by a dynamically accessible perturbation with direction γ > e µ ( t, , completingthe first oscillation. This behavior of particles getting expelled and return-ing repeats several times until the system eventually oscillates around aseemingly unchanging configuration. In particular, the oscillating behaviorcannot only be seen in e µ ( t, , but also on the particle level as well as in27 . . r t Figure 10: The weighted mass density 4 πr ρ for the King model inSchwarzschild coordinates for y = 0 . > y max perturbed by a dynami-cally accessible perturbation with direction γ >
0. The color represents thevalue of 4 πr ρ , increasing from black to yellow.all metric and matter quantities. The oscillation of the mass density, i.e.,the radial derivative of the Hawking mass, in the case of the King model isdepicted in Figure 10.Another illustration of the effect described above is given in Figure 11where we plot the energy E of individual particles against their distancefrom the origin. Thereby, we can investigate trajectories of particles. Notethat the scales of both the energy E and the radius r change from the plotfor t = 0 to the plot for t = 15 and stay the same for t >
15. Plotting thesespecific quantities allows us to identify a portion of low energy particlesthat seems to form some dense structure which is independent of the lessdense particles with higher energy. In addition, the dispersion of clusters ofparticles can be observed in both Figure 10 and 11.A possible interpretation of these observations is that after perturbationa part of the solution corresponding to high energy particles disperses whileafter some transition period the remainder starts oscillating around a new28 . . . . . . . . . .
95 0 0 . . . . E r t = 0 . . . . .
05 0 1 2 3 4 5 E r t = 15 . . . . .
05 0 1 2 3 4 5 E r t = 30 . . . . .
05 0 1 2 3 4 5 E r t = 45 . . . . .
05 0 1 2 3 4 5 E r t = 120 . . . . .
05 0 1 2 3 4 5 E r t = 180 Figure 11: Numerical particles for the King model in Schwarzschild coordi-nates for y = 0 . > y max perturbed by a dynamically accessible perturba-tion with direction γ >
0. The movie corresponding to these snapshots isavailable at [31]. 29teady state. This new state is less relativistic than the original one andseems to depend solely on the original equilibrium, i.e., it is independent ofnumerical parameters and the type and strength of the perturbation used toperturb the original unstable steady state, at least for sufficiently weak per-turbations. In the context of dynamical systems such a new state is usuallya stable steady state of the system, and the above behavior is reminiscent ofa heteroclinic orbit, as the original steady state after perturbation migratesto a different one. On the other hand, it does not really seem to convergetowards this new state, as would be the case for a genuine heteroclinic orbit,so we use this terminology very loosely.Since the Einstein-Vlasov system possesses a plethora of steady states,the explicit identification of the target state seems to be very delicate.We developed a method which acts on the level of particles and basically“deletes” all the “departed” particles from our solution before fitting theremainder by some stable steady state. However, a weak point of our pro-cedure is that we have to decide manually whether or not a particle belongsto the “remainder”. In our simulations, we always made this decision basedon the radius or the particle energy. All this being said the observationswith the above procedure seem to support the “heteroclinic orbit picture”.Developing appropriate criteria for the latter question would definitely openup new possibilities for the search of these target states. In passing we notethat Figure 11 actually presents snapshots from a movie which illustratesthe above behavior and which, together with similar simulations, can beviewed via the link [31].As can be seen in Figure 9, the difference between the value of e µ ( t, around which the solution oscillates at later times and the initial value e µ (0 , increases in y . In particular, the initial elevation also increases in y . Thisis caused by a relative increase of the number and mass of particles whichinitially get expelled from the configuration as well as a decrease of thenumber of particles which return. This difference seems to change smoothlyin y which causes the initial elevation to be very subtle for y ≈ y max ,cf. Figure 6, consistent with the fact that no such elevation occurs in thestable regime. Furthermore, the time-frame of the initial reimplosion of thesolution increases rapidly in y when choosing y large enough, which meansthat the existence of fully dispersing solutions induced by large y can notbe ruled out numerically.It has been suggested in [28] that perturbations promoting dispersion ofcertain isotropic steady states with negative binding energy lead to fully dis-persing solutions. In [3], this suggestion has been extended to non-isotropicsteady states with an inner vacuum region, i.e., L >
0. However, our results30isprove this conjecture in all three coordinate systems under consideration.In the case of the King model, the binding energy of a steady state is nega-tive for y > .
784 which means that Figure 9 clearly shows the reimplosionand oscillation of solutions emerging from weakly perturbed steady stateswith negative binding energy. The same phenomenon can also be observedfor other isotropic and non-isotropic models, at least if L is not too large.In the case of a polytropic ansatz with k = 0 . l = 0 . L = 0 .
001 this isdepicted in Figure 12. . . . . . .
91 0 100 200 300 400 500 e µ ( t , ) t y = 0 . y = 0 . y = 0 . y = 0 . y = 0 . y = 0 . Figure 12: The polytropic ansatz with k = 0 . l = 0 . L = 0 .
001 inSchwarzschild coordinates with y > y max perturbed by a dynamically ac-cessible perturbation with direction γ > L is closely connected to the size of the inner vacuum regionof steady states. Increasing this parameter in the ansatz function seems toslow down all occurring effects, in particular causing the reimplosion timeto increase. This effect is the main reason why we choose L rather small.In fact, when setting L significantly larger, we cannot decide whether somesteady states with y max < y < y zero reimplode or fully disperse. A similareffect can be observed for large k . Furthermore, large y pose a problemwhen examining dispersing steady states for every model since we cannot31umerically distinguish solutions with large reimplosion times from fullydispersing ones. It therefore remains an open question whether or not thereexist fully dispersing solutions emerging from weakly perturbed unstablesteady states and if this behavior is connected to the binding energy of theoriginal steady states in some way. We conclude our paper with an overview of the parameter settings whichinfluence the numerical accuracy. Compared to former numerical projects inthis field of research, we benefit from more computational power and froma parallelized code. This allows us to perform computations with a largenumber of numerical particles within a reasonable computation time. Theactual computation time further depends on the discretization of the timevariable t and the spatial variable r . We typically use ∆ r and ∆ t in theorder of magnitude of 10 − . Moreover, we ensure that at least 15 millionnumerical particles are used for our computations by choosing the number ofsteps N u and N ψ suitably after fixing ∆ r and ∆ t . In order to guarantee thatperturbations resulting of errors due to our initialization are small comparedto our applied perturbation and that the applied perturbation is reasonablysmall we choose T pert = 0 . ǫ pert of the order 10 − and determine asuitable γ as described in Section 3.To monitor the validity of our simulation, we keep track of the ADMmass M ( t ) and the analytical number of particles N ( t ) which are conservedquantities along solutions of the Einstein-Vlasov system. We define therelative errors as e M ( t ) = | M ( t ) − M | M , e N ( t ) = | N ( t ) − N | N where M = M ( T pert ) and N = N ( T pert ) are given by the perturbed steadystate, i.e., the evolved state at t = T pert . In the case of oscillating or hetero-clinic solutions the errors stay very small—of the order 10 − until t = 100—when choosing the numerical parameters as mentioned above. When con-sidering the formation of a trapped surface and the formation of a blackhole, it turns out that the errors become larger and—not surprisingly—thesimulation eventually breaks down.A further test of our codes is to evolve an unperturbed steady state.When choosing a stable steady state it is tracked faithfully for very longtimes. Obviously, for an unstable steady state the errors due to the initial-ization can eventually cause a deviation from the steady state. In conclusion,32t seems fair to say that our simulations provide conclusive results at leastif we are not considering the long time behavior of collapsing solutions aftertrapped surfaces have formed. References [1]
Andr´easson, H. , The Einstein-Vlasov System/Kinetic Theory.
LivingReviews in Relativity , 4 (2011).[2] Andr´easson, H., Kunze, M., Rein, G. , Global existence for thespherically symmetric Einstein-Vlasov system with outgoing matter,
Comm. Partial Differential Eqns. , 656–668 (2008).[3] Andr´easson, H., Rein, G. , A numerical investigation of the stabilityof steady states and critical phenomena for the spherically symmetricEinstein-Vlasov system.
Class. Quantum Gravity , 3659–3677 (2006).[4] Andr´easson, H., Rein, G.
On the steady states of the sphericallysymmetric Einstein-Vlasov system.
Class. Quantum Gravity , 1809–1832 (2007).[5] Andr´easson, H., Rein, G. , The asymptotic behaviour inSchwarzschild time of Vlasov matter in spherically symmetric gravita-tional collapse.
Math. Proc. Camb. Phil. Soc. , 173–188 (2010).[6]
Andr´easson, H., Rein, G. , Formation of trapped surfaces for thespherically symmetric Einstein-Vlasov system.
J. Hyperbolic Differ. Equ. , 707–731 (2010).[7]
Bisnovatyi-Kogan, G. S., Zel’dovich, Ya. B. , Models of clusters ofpoint masses with great central red shift.
Astrofizika , , 223–234 (1969).[8] G¨unther, S. , The Einstein-Vlasov system in maximal areal coordi-nates. Master thesis, Bayreuth 2019.[9]
Hadˇzi´c, M., Lin, Z. , Turning point principle for relativistic stars.Preprint, https://arxiv.org/abs/2006.09749 (2020).[10]
Hadˇzi´c, M., Rein, G. , Stability for the spherically symmetricEinstein-Vlasov system—a coercivity estimate.
Math. Proc. Camb. Phil.Soc. , 529–556 (2013). 3311]
Hadˇzi´c, M., Rein, G. , On the small redshift limit of steady statesof the spherically symmetric Einstein-Vlasov system and their stability.
Math. Proc. Camb. Phil. Soc. , 529–546 (2015).[12]
Hadˇzi´c, M., Lin, Z., Rein, G. , Stability and instabil-ity of self-gravitating relativistic matter distributions. Preprint,https://arxiv.org/abs/1810.00809 (2018).[13]
Ipser, J. R. , Relativistic, Spherically Symmetric Star Clusters. III.Stability of Compact Isotropic Model.
Astrophys. J. , 17–43 (1969).[14]
Ipser, J. , A binding-energy criterion for the dynamical stability ofspherical stellar systems in general relativity.
Astrophys. J. , 1101–1110 (1980).[15]
Ipser, J., Thorne, K. S. , Relativistic, spherically symmetric starclusters I. Stability theory for radial perturbations.
Astrophys. J. ,251–270 (1968).[16]
Lin, Z., Zeng, C: , Separable Hamiltonian PDEs and Turn-ing point principle for stability of gaseous stars. Preprint,https://arxiv.org/abs/2005.00973 (2020).[17]
Korch, M., Ramming, R., Rein, G. , Parallelization of Particle-in-Cell Codes for Nonlinear Kinetic Models from Mathematical Physics.
ICPP 13: Proceedings of the 2013 42nd International Conference onParallel Processing , 523–529 (2013).[18]
I. Olabarrieta, M. W. Choptuik , Critical phenomena at thethreshold of black hole formation for collisionless matter in sphericalsymmetry,
Phys. Rev. D. , 024007 (2002).[19] Penrose, R. , Gravitational Collapse and Space-Time Singularities
Physical Review Letters , 57–59 (1965).[20]
Ramming, R., Rein, G. , Spherically symmetric equilibria for self-gravitating kinetic or fluid models in the nonrelativistic and relativisticcase—a simple proof for finite extension.
SIAM J. Math. Analysis ,900–914 (2013).[21] Ramming, R., Rein, G. , Oscillating solutions of the Vlasov-Poissonsystem—A numerical investigation.
Phys. D , 72–79 (2018).3422]
Rein, G. , The Vlasov-Einstein System with Surface Symmetry , Habili-tationsschrift, M¨unchen 1995.[23]
Rein, G. , Collisionless kinetic equations from astrophysics—TheVlasov-Poisson system. In
Handbook of Differential Equations, Evolu-tionary Equations, vol. 3 , edited by C. M. Dafermos and E. Feireisl,Elsevier (2007).[24]
Rein, G., Rendall, A. , Global existence of solutions of the spher-ically symmetric Vlasov-Einstein system with small initial data.
Com-mun. Math. Phys. , 561–583 (1992). Erratum:
Commun. Math. Phys. , 475–478 (1996).[25]
Rein, G., Rendall, A. D., Schaeffer, J. , Critical collapse of col-lisionless matter: A numerical investigation.
Phys. Rev. D , 044007(1998).[26] Rein, G., Rodewis, T. , Convergence of a particle-in-cell scheme forthe spherically symmetric Vlasov-Einstein system,
Indiana UniversityMath. J. , 821–862 (2003).[27] Schaeffer, J.
Discrete approximation of the Poisson-Vlasov system.
Q. Appl. Math. , 59–73 (1987).[28] Shapiro, S. L., Teukolsky, S. A. , Relativistic Stellar Dynamics onthe Computer—Part Two—Physical Applications.
Astrophysical Journal , 58 (1985).[29]
Zel’dovich, Y. B., Podurets, M. A. , The Evolution of a System ofGravitationally Interacting Point Masses.
Soviet Astronomy , 742–749(1966).[30] Zel’dovich, Y. B., Novikov, I. D. , Relativistic astrophysics. Vol.1:Stars and relativity.
Chicago: University of Chicago Press