A numerical investigation of the stability of steady states and critical phenomena for the spherically symmetric Einstein-Vlasov system

0.26 MB·English

by  Hakan Andreasson

A numerical investigation of the stability of steady
states and critical phenomena for the spherically
6 symmetric Einstein-Vlasov system
0
0
2 H˚akan Andr´easson

n Department of Mathematics, Chalmers,
a
S-41296 G¨oteborg, Sweden
J
6 email: [email protected]
2

Gerhard Rein
1
v Department of Mathematics, University of Bayreuth,
2 D-95440 Bayreuth, Germany
1
email: [email protected]
1
1
0 February 7, 2008
6
0
/
c Abstract
q
- The stability features of steady states of the spherically symmetric
r
g Einstein-Vlasov system are investigated numerically. We find support
: for the conjecture by Zeldovich and Novikov that the binding energy
v
i maximumalongasteadystatesequencesignalstheonsetofinstability,
X
a conjecture which we extend to and confirm for non-isotropic states.
r The sign of the binding energy of a solution turns out to be relevant
a
for its time evolution in general. We relate the stability properties to
thequestionofuniversalityincriticalcollapseandfindthatforVlasov
matter universality does not seem to hold.
1 Introduction
For any given dynamical system the existence and stability of steady states
is essential both from a mathematics and from a physics point of view. In
this paper we investigate these questions by numerical means for the spher-
ically symmetric, asymptotically flat Einstein-Vlasov system. This system
describes a self gravitating collisionless gas in the framework of general rela-
tivity. Here the matter is thought of as a large ensemble of particles, which
1
is described by a density function on phase space, and the individual par-
ticles move along geodesics. The precise formulation of this system will be
given in the next section; for further information on the Einstein-Vlasov
system we refer to [1]. In astrophysics the system is used to model compact
star clusters and galaxies. In this context the stability question was first
studied by Zeldovich e. a. in the sixties [38, 37]. These authors character-
ize a steady state by its central redshift and binding energy and conjecture
that the binding energy maximum along a steady state sequence signals the
onset of instability. For isotropic steady states Ipser [17] and Shapiro and
Teukolsky [33] found numerical support for this conjecture. In our investi-
gation we find that this conjecture holds for non-isotropic steady states as
well, where the density on phase space depends on the particle energy and
angular momentum. Moreover, we usethreedifferent kindsof perturbations
(fordetailsseeSection3)whileinthepreviousstudiesmentionedaboveonly
one type of perturbation was implemented. We also find that the sign of
the binding energy is crucial for the evolution of the perturbed solutions. A
positive value implies that the solution is bound in the sense that not all
matter can disperse to infinity, and in this case the perturbation of a stable
state seems to lead to periodic oscillations. In the case of a negative value
of the binding energy we observe that the solution disperses to infinity.
Our initial motivation for studying the stability of steady states was its
role in critical collapse. This topic started with the work of Choptuik [7]
where he studied the Einstein-scalar-field system. He took a fixed initial
profile for the scalar field and scaled it by an arbitrary constant factor. This
gives rise to a family of initial data depending on a real parameter A. It
turned out that there exists a critical parameter A such that for A < A
∗ ∗
the corresponding solutions disperses as predicted by the theoretical results
of Christodoulou [9] for small data, while for A > A the corresponding
∗
solutions collapse and produce a black hole in accordance with [10]. The
surprising result was that the limit of the mass M(A) of the black hole
tends to zero for A → A so that in such a one-parameter family there are
∗
black holes with arbitrarily small mass. Choptuik found that this fact is
related to the existence of self-similar solutions of the Einstein-scalar-field
system, in particular, the critical solution is self-similar and universal, i.e.,
independentoftheinitialprofilewhichdeterminestheone-parameterfamily.
Later on Choptuik, Chmaj, and Bizon´ performed a similar investigation for
the Einstein-Yang-Mills system [8]. Here both cases lim M(A) = 0
A→A∗
and lim M(A) > 0 were found and called type II and type I
A→A∗,A>A∗
respectively. In the latter case there is a mass gap in the M(A)-curve. The
possibility of type I behavior in this system was related to the existence
2
of the Bartnik-McKinnon solutions [5, 34], which are static. Again, the
conjecture is that the critical solution is universal.
As opposed to the field theoretic matter models mentioned above the
Vlasov model is phenomenological, and in contrast to fluid models several
globalresultshavebeenobtainedforVlasovmatter. Forthesphericallysym-
metric, asymptotically flat case it was shown in [23] that sufficiently small
initial data launch global, geodesically complete solutions which dispersefor
large times. It is also known that there do exist initial data, necessarily
large, which develop singularities [29]. The proof relies on the Penrose sin-
gularity theorem. There are no general results on the behavior of large data
solutions yet, except for the following: If data on a hypersurface of constant
Schwarzschild time give rise to a solution which develops a singularity after
a finite amount of Schwarzschild time, then the first singularity occurs at
the center of symmetry [26]. An analogous result where Schwarzschild time
is replaced by maximal slicing has also been proved [29]. In [2] further re-
sults on the global behavior of solutions for large data can be found. The
transition between dispersion and gravitational collapse was numerically in-
vestigated by Rein, Rendall, and Schaeffer [27], and it was found that there
is amassgap intheM(A)curve. Thisresultwas later confirmedby Olabar-
rieta and Choptuik [19]. In addition, the latter authors reported evidence
that the mass gap is due to the presence of static solutions, and they found
support that the critical solution is universal.
Inthepresentinvestigation weaddresstheroleofsteadystatesincritical
phenomena for the Einstein-Vlasov system and the question of universality
by explicitly exploiting the fact that for this system the existence of an
abundance of steady states is well-established [25], and that these steady
states can easily be computed numerically. Computing a steady state f we
s
consider the family Af of initial data. Within every family of steady states
s
givenbyaspecificdependenceonparticleenergyandangularmomentumwe
find unstable ones which act as critical solutions: If they are perturbedwith
A > 1 they collapse to a black hole, if they are perturbed with A < 1 they
either disperse or oscillate in a neighborhood of the steady state, depending
on the sign of the binding energy. Due to the abundance of possible such
dependences on particle energy and angular momentum there cannot be a
universal critical solution in spherically symmetric collapse for the Einstein-
Vlasov system.
The paper proceeds as follows: In the next section we formulate the
Einstein-Vlasov system, first in general coordinates and then in coordinates
adaptedtosphericalsymmetry. In[27]Schwarzschildcoordinateswereused.
Here we use maximal areal coordinates, as was done in [19]. This has the
3
advantage that regions of spacetime containing trapped surfaces can becov-
ered. In Section 3 we discuss the numerical scheme which we use; it is a
particle-in-cell scheme of the typeused for kinetic models in plasma physics.
It should be noted that for the analogous scheme used in [27] a rigorous
convergence proof has been established in [28], and we conjecture that this
can be done for the present scheme as well. In Section 4 we investigate the
stability properties of certain steady states and compare our findings with
the earlier work mentioned above. Section 5 is devoted to critical phenom-
ena and non-universality, and in a final section we discuss the reliability of
our code.
Our main motivation for this numerical analysis is that it may lead to
conjectures onthebehaviorofsolutionsoftheEinstein-Vlasov systemwhich
may eventually be proven rigorously.
2 The Einstein-Vlasov system
We first write down the Einstein-Vlasov system in general coordinates on
the tangent bundle TM of the spacetime manifold M. Following standard
practice we normalize the physical constants to one. The system then reads
as follows:
pa∂xaf −Γabcpbpc∂paf = 0,
Gab = 8πTab,
d4p
Tab = papbf |g|1/2 .
m
Z
Heref isthenumberdensityoftheparticles onphasespace, Γa andGab de-
bc
notetheChristoffelsymbolsandtheEinsteintensorobtainedfromthespace-
time metric g , |g| denotes its determinant, Tab is the energy-momentum
ab
tensor generated by f, xa are coordinates on M, (xa,pb) the corresponding
coordinates on the tangent bundle TM, Latin indices run from 0 to 3, and
m = |g papb|1/2
ab
is the rest mass of a particle at the corresponding point in phase space. We
assume that all particles have rest mass 1 and move forward in time so that
the distribution function f lives on the mass shell
PM = g papb = −1, p0 > 0 .
ab
n o
4
We consider this system in the spherically symmetric, asymptotically flat
case and write the metric in the following form:
ds2 = −(α2+a2β2)dt2+2a2βdtdr+a2dr2+r2 dθ2+sin2θdφ2 .
Here the metric coefficients α,β, and a depend on t(cid:0)∈ R and r ≥ 0,(cid:1)α and
a are positive, and the polar angles θ ∈ [0,π] and φ ∈ [0,2π] parameterize
the unit sphere. The radial coordinate r is thus the area radius. Let Ka be
b
the second fundamental form and define
β
κ:= Kθ = .
θ rα
By imposing the maximal gauge condition, which means that each hyper-
surface of constant t has vanishing mean curvature, we obtain the following
field equations:
3 a
a = a3rκ2+4πra3ρ+ (1−a2), (1)
r
2 2r
3
κ = − κ−4πa, (2)
r
r
a = 2αaκ+(aβ) , (3)
t r
a 2 2α a
α = α r − + 2r r +a2−1 +4πa2α(S −3ρ). (4)
rr r a r r2 a
(cid:18) (cid:19)
(cid:16) (cid:17)
The Vlasov equation takes the form
αw α ǫ αL
r
∂ f + −β ∂ f + − −2ακw+ ∂ f = 0, (5)
t aǫ r a ar3ǫ w
(cid:18) (cid:19)
(cid:16) (cid:17)
where
ǫ =ǫ(r,w,L) = 1+w2 +L/r2. (6)
The variables w and L can be thougpht of as the momentum in the radial
direction and the square of the angular momentum respectively, see [20] for
more details. The matter quantities are defined by
π ∞ ∞
ρ(t,r) = ǫf(t,r,w,L)dLdw, (7)
r2
Z−∞Z0
π ∞ ∞
(t,r) = wf(t,r,w,L)dLdw, (8)
r2
Z−∞Z0
π ∞ ∞ w2+L/r2
S(t,r) = f(t,r,w,L)dLdw. (9)
r2 ǫ
Z−∞Z0
5
We impose the following boundary conditions which ensure asymptotic flat-
ness and a regular center:
a(t,0) = a(t,∞) =α(t,∞) = 1. (10)
Theequations(1)–(9)together withtheboundaryconditions(10)constitute
the Einstein-Vlasov system for a spherically symmetric, asymptotically flat
spacetime in maximal areal coordinates.
Let us consider some simple properties of this system, in particular such
as will berelevant for its numerical simulation; a more careful mathematical
analysis of the system will be performed elsewhere [3]. First we note that
the phase space density f is constant along solutions of the characteristic
system
α(τ,r)w
r˙ = −β(τ,r), (11)
a(τ,r)ǫ
α (τ,r)ǫ α(τ,r)L
r
w˙ =− −2α(τ,r)κ(τ,r)w+ , (12)
a(τ,r) a(τ,r)r3ǫ
L˙ =0 (13)
of the Vlasov equation. If τ 7→ (R,W,L)(τ,t,r,w,L) =: Z(τ,t,z) de-
notes the solution of the characteristic system with (R,W,L)(t,t,r,w,L) =
(r,w,L) then
◦
f(t,r,w,L) = f((R,W,L)(0,t,r,w,L)),
◦
with f = f(0,·) the initial datum for f. Due to our choice of coordinates
w,L in momentum space the characteristic flow is not measure preserving.
Indeed, ifbyD wedenotedifferentiation along acharacteristic of theVlasov
equation, then
αw 2β
D(f dLdwdr)= − a +β + fdLdwdr; (14)
a2ǫ r r r
(cid:18) (cid:19)
the factor on the right hand sideis the (r,w,L)-divergence of the right hand
side of the characteristic system. Notice that by Eqn. (3) this factor equals
−Dlna. If A ⊂ R × R × R is measurable and A(t) := Z(t,0,A) =
+ +
{(R,W,L)(t,0,r,w,L) |(r,w,L) ∈ A} then
◦ a(0,r)
f(t,z)dz = f(z) dz,
a(t,R(t,0,z))
ZA(t) ZA
while
◦
a(t,r)f(t,z)dz = a(0,r)f(z)dz.
ZA(t) ZA
6
In particular, the total number of particles, which, since all particles have
rest mass one, equals the rest mass of the system, is conserved:
∞ ∞ ∞
4π2 a(t,r)f(t,r,w,L)dLdwdr = M .
0
Z0 Z−∞Z0
Next we notice that Eqn. (2) can be rewritten as
r3κ = −4πr3a. (15)
r
(cid:0) (cid:1)
Using Eqn. (1) the second order equation (4) can be rewritten as
r2
α = 4πr2aα(ρ+S)+6r2aακ2. (16)
r
a
(cid:18) (cid:19)r
The Hawking mass m is given by
r β2 1
m = 1+ − . (17)
2 α2 a2
(cid:18) (cid:19)
We also introduce the quantity
r 1
µ := 1− ,
2 a2
(cid:18) (cid:19)
and note that by Eqn. (1), µ can be written in the form
r 3
µ(t,r) = 4πρ(t,s)+ κ2(t,s) s2ds.
2
Z0 (cid:18) (cid:19)
Assuming that the matter is compactly supported initially and hence also
for later times Eqn. (15) implies that κ(t,r) ∼ r−3 for r large. Hence the
limits as r tends to ∞ of m and µ are equal so that the ADM mass M can
be written as
∞ 3
M = 4πρ(t,r)+ κ2(t,r) r2dr.
2
Z0 (cid:18) (cid:19)
The ADM mass is conserved.
It is of interest to find out when trapped surfaces form, which means
that 2m/r > 1, and in view of (17) this condition can be written as
α2
a2 > .
β2
7
Note also that 2m/r < 1/2 implies that 2µ/r < 1/2. From [18] we obtain
the following inequalities
1 β 1 β
− ≤ 1, + ≤ 1. (18)
a α a α
(cid:12) (cid:12) (cid:12) (cid:12)
(cid:12) (cid:12) (cid:12) (cid:12)
(cid:12) (cid:12) (cid:12) (cid:12)
Finally we notice that th(cid:12)e secon(cid:12)d order(cid:12) equati(cid:12)on for α in the form (16) can
be integrated:
a(t,r) r
α (t,r) = 4παa(ρ+S)+6aακ2 s2ds. (19)
r r2
Z0
(cid:0) (cid:1)
Thus α is monotonically increasing outwards, and from (10) and (18) it
follows that
0 < α ≤ 1, |β| ≤ α≤ 1.
3 The numerical method
◦ ◦
Let usconsider an initial condition f= f(r,w,L), defineadditional variables
u ≥ 0 and φ ∈[0,π] by the relations
L
u2 = w2+ , w = u cosφ, L = r2u2sin2φ
r2
◦
and assume that f vanishes outside the set (r,u,φ) ∈ [R ,R ]×[U ,U ]×
0 1 0 1
[Φ ,Φ ]. We will approximate the solution using a particle method. For a
0 1
thorough treatment of particle methods in the context of plasma physics we
refer to [6]. To initialize the particles we take integers N ,N ,N and define
r u φ
R −R U −U Φ −Φ
1 0 1 0 1 0
∆r = , ∆u= , ∆φ = ,
N N N
r u φ
1 1 1
r = R + i− ∆r, u = U + j − ∆u, φ = Φ + k− ∆φ,
i 0 j 0 k 0
2 2 2
(cid:18) (cid:19) (cid:18) (cid:19) (cid:18) (cid:19)
◦
f0 = f(r ,u ,φ ) 4πr2∆r 2πu2∆u sinφ ∆φ,
ijk i j k i j k
r0 = r , w0 = u cosφ , L = (r u sinφ )2.
ijk i ijk j k ijk i j k
At this point it should be noted that f0 contains the phase space volume
ijk
element, afactwhichisconvenientwhencomputingtheinducedcomponents
of the energy momentum tensor, but which has to be observed when this
quantity is time-stepped. At each point in phase space with coordinates
8
r0 ,w0 ,L we imagine a particle with weight f0 which is smeared out
ijk ijk ijk ijk
in the radial direction, i.e., with respect to r it is represented by a hat
function of width 2∆r. Once the particles are initialized the grid covering
the supportof the initial datum plays nofurther role. From these numerical
particles, approximations are made of the quantities ρ, , S defined in (7),
(8), (9) at the spatial grid points
r := j∆r, j = 0,...,N +1;
j
note that this is now a new grid, which we extend from r = 0 to the radius
of the spatial support of the matter quantities. The field equations (1) and
(15), which do not contain the metric coefficient α, can now be integrated
on the support of the matter from the origin outward, using the boundary
conditions (10) and (r3κ) = 0; notice that κ is then also known outside
|r=0
the support of the matter.
To obtain an approximation of the metric coefficient α at the grid points
r we discretize Eqn. (16) in the following way: We write
j
1
r 1 := j + ∆r,
j+2 2
(cid:18) (cid:19)
and from what is already computed we obtain approximations for a, κ, ρ, S
both at the grid points r and by interpolation at the intermediate points
j
r 1, which we denote by
j+
2
aj, aj+1, ...
2
Using the approximations
r2α (r ) ≈ 1 rj2+21αr(rj+21) − rj2−21αr(rj−21)
r j
a ∆r  a 1 a 1 
(cid:18) (cid:19)r j+2 j−2
 
and
α −α α −α
j+1 j j j−1
αr(rj+12) ≈ ∆r , αr(rj−21) ≈ ∆r
we obtain the following linear system for the approximations α for α(r )
j j
on the grid points: For j = 0,...,N −1,
r2 r2
1 1
j+ j−
2 α + 2 α −
j+1 j−1
a 1 a 1
j+ j−
2 2
r2 r2
1 1
− j+2 + j−2 + 6κ2+4π(ρ +S ) r2(∆r)2a α = 0,
a 1 a 1 j j j j j j
j+ j−
2 2
(cid:2) (cid:3)
 
9
and
2M∗
α −α =0, α = 1− .
1 0 N+1
s rN+1
The last two equations arise from the boundary condition
r2
α (t,0) = 0
r
a
(cid:18) (cid:19)r
and the approximation
2µ(t,r)
α(t,r) ≈ 1−
r
r
forr large, inparticular,r welloutsidethesupportofthematter. Thelatter
approximation is motivated by the known representation of the asymptoti-
cally flat Schwarzschild solution in maximal areal coordinates. Theapproxi-
mation M∗ for µ(t,r) at large values of r can be computed from the already
approximated quantities ρ and κ.
The linear system for α is tridiagonal, obviously diagonally dominated,
j
and it can easily be solved. Thus we now have approximations for all the
field quantities on the spatial grid points r . In passing we note that in [19]
j
thefieldequation (4)was discretized directly resultingagain inatridiagonal
system for α , but our approach of discretizing the rewritten version (16)
j
instead seems to perform better numerically.
Toperformthetimestep wepropagate thenumericalparticles according
to the characteristic system of the Vlasov equation (11), (12), (13). To do
so we interpolate the field quantities to particle locations and use a sim-
ple Euler time stepping method to define the new phase space coordinates
r1 ,w1 ,L1 ofthenumericalparticlewithlabelijk. Westillneedtoprop-
ijk ijk ijk
agate thephasespacevolumeelement alongthecharacteristics. Discretizing
(14) in time we obtain the relation
1 αw 2β
f1 −f0 = − a +β + f0
∆t ijk ijk a2ǫ r r r ijk
(cid:18) (cid:19)
(cid:0) (cid:1)
which we use to compute f1 , and one time step is complete.
ijk
4 Stability issues for steady states
For static solutions, β = κ =  = 0 so that
r
m(r)= µ(r)= 4π ρ(s)s2ds,
Z0
10