Introduction to relativistic astrophysics and cosmology through Maple

Vladimir L. Kalashnikov ,

Belarussian Polytechnical Academy,

[email protected]

www.geocities.com/optomaplev

Abstract: The basics of the relativistic astrophysics including the celestial mechanics in weak field, black holes and cosmological models are illustrated and analyzed by means of Maple 6

Application Areas/Subjects: Science, Astrophysics, General Relativity, Tensor Analysis, Differential geometry, Differential equations

Relativistic stars and black holes

Introduction

The most wonderful prediction of GR-theory is the existence of black holes, which are the objects with extremely strong gravitational field. The investigation of these objects is the test of our understanding of space-time structure. We will base our consideration on the analytical approach that demands to consider only symmetrical space-times. But this restriction does not decrease the significance of the obtained data because of the rich structure of analytical results and possibilities of clear interpretation clarify the physical basis of the phenomenon in the strongly curved space-time. The basic principles can be found in C. W. Misner, K. S. Thorne, J. A. Wheeler, Gravitation, San Francisco, 1973 , V. Frolov, I. Novikov, Physics of Black Holes, Kluwer, Dordrecht, 1998 .

Geometric units

The very useful normalization in GR utilizes the so-called geometric units. Since the left-hand side of the Einstein equations describes the curvature tensor (its dimension is ), the right-hand side is to have same dimension. Let's the gravitational constant is G = = 1 and the light velocity is c = = 1,

then

G / = cm/g = 1

/G = erg/s = 1 (power unit)

G / c = Hz* / g = 1 (characteristic of absorption)

/ = CGSE units (field strength)

h /2p = g* /s =

elementary charge e = cm

1 ps = cm

1 eV = cm

There are the following extremal values of length, time, mass and density (the so-called Planck units), which are useful in the context of the consideration of GR validity:

= cm (Planck length)

= s (Planck time)

= g (Planck mass)

= g / (Planck density)

Relativistic star

As stated above the first realistic metric was introduced by Schwarzschild for description of the spherically symmetric and static curved space. Let us introduce the spherically symmetric metric in the following form:

> restart:
with(tensor):
with(plots):
with(linalg):
with(difforms):

coord := [t, r, theta, phi]:# spherical coordinates, which will be designated in text as [0,1,2,3]
g_compts := array(symmetric,sparse,1..4,1..4):# metric components
g_compts[1,1] := -exp(2*Phi(r)):# component of interval attached to d(t)^2
g_compts[2,2] := exp(2*Lambda(r)):# component of interval attached to d(r)^2
g_compts[3,3] := r^2:# component of interval attached to d(theta)^2
g_compts[4,4] := r^2*sin(theta)^2:# component of interval attached to d(phi)^2

g := create([-1,-1], eval(g_compts));# covariant metric tensor

ginv := invert( g, 'detg' );# contravariant metric tensor

Now we can use the standard Maple procedure for Einstein tensor definition

> # intermediate values
D1g := d1metric( g, coord ):
D2g := d2metric( D1g, coord ):
Cf1 := Christoffel1 ( D1g ):
RMN := Riemann( ginv, D2g, Cf1 ):
RICCI := Ricci( ginv, RMN ):
RS := Ricciscalar( ginv, RICCI ):

> Estn := Einstein( g, RICCI, RS ):# Einstein tensor
displayGR(Einstein,Estn);# Its nonzero components

In the beginning we will consider the star in the form of drop of liquid. In this case the energy-momentum tensor is = (p+r ) + p (all components of u except for are equal to zero, and - 1 = ; the signature (-2) results in =1 and = (p+r ) - p ):

> T_compts := array(symmetric,sparse,1..4,1..4):# energy-momentum tensor for drop of liquid
T_compts[1,1] := exp(2*Phi(r))*rho(r):
T_compts[2,2] := exp(2*Lambda(r))*p(r):
T_compts[3,3] := p(r)*r^2:
T_compts[4,4] := p(r)*r^2*sin(theta)^2:
T := create([-1,-1], eval(T_compts));

To write the Einstein equations (sign corresponds to W.Pauli, Theory of relativity, Pergamon Press (1958) )

= -8 p

let's extract the tensor components:

> Energy_momentum := get_compts(T);
Einstein := get_compts(Estn);

First Einstein equation for (0,0) - component is:

> 8*Pi*Energy_momentum[1,1] + Einstein[1,1]:
expand(%/exp(Phi(r))^2):
eq1 := simplify(%) = 0;#first Einstein equation

This equation can be rewritten as:

> eq1 := -8*Pi*rho(r)*r^2+Diff(r*(1-exp(-2*Lambda(r))),r) = 0;

The formal substitution results in:

> eq1_n := subs(\
r*(1-exp(-2*Lambda(r))) = 2*m(r),\
lhs(eq1)) = 0;

> dsolve(eq1_n,m(r));

So, m(r) is the mass inside sphere with radius r. To estimate _C1 we have to express L from m :

> r*(1-exp(-2*Lambda(r))) = 2*m(r):
expand(solve(%, exp(-2*Lambda(r)) ));

Hence (see expression for through L) the Lorenzian metric in the absence of matter is possible only if _C1=0 .

Second Einstein equation for (1,1)-component is:

> eq2 := simplify( 8*Pi*Energy_momentum[2,2] + Einstein[2,2] ) = 0;#second Einstein equation

> eq2_2 := numer(\
lhs(\
simplify(\
subs(exp(2*Lambda(r)) = 1/(1-2*m(r)/r),eq2)))) = 0;

As result we have:

> eq2_3 := Diff(Phi(r),r) = solve(eq2_2, diff(Phi(r),r) );

We can see, that the gradient of the gravitational potential F is greater than in the Newtonian case = , that is the pressure in GR is the source of gravitation.

For further analysis we have to define the relativist equation of hydrodynamics (relativist Euler's equation):

> compts := array([u_t,u_r,u_th,u_ph]):
u := create([1], compts):# 4-velocity
Cf2 := Christoffel2 ( ginv, Cf1 ):
(rho(r)+p(r))*get_compts(cov_diff( u, coord, Cf2 ))[1,2]/(u_t) = -diff(p(r),r);# radial component of Euler equation, u_r=u_th=u_ph=0
eq3 := Diff(Phi(r),r) = solve(%, diff(Phi(r),r) );

As result we obtain the so-called Oppenheimer-Volkoff equation for hydrostatic equilibrium of star:

> Diff(p(r),r) = factor( solve(rhs(eq3) = rhs(eq2_3),diff(p(r),r)) );

One can see that the pressure gradient is greater than in the classical limit ( = - ) and this gradient is increased by pressure growth (numerator) and r decrease (denominator) due to approach to star center. So, one can conclude that in our model the gravitation is stronger than in Newtonian case.

Out of star m ( r ) =M , p= 0 ( M is the full mass of star). Then

> diff(Phi(r),r) = subs({m(r)=M,p(r)=0},rhs(eq2_3));
eq4 := dsolve(%, Phi(r));

The boundary condition

> 0 = limit(rhs(eq4),r=infinity);

results in

> subs(_C1=0,eq4);

So, Schwarzschild metric out of star is:

> g_matrix := get_compts(g);

> coord := [t, r, theta, phi]:
sch_compts := array(symmetric,sparse,1..4,1..4):# metric components
sch_compts[1,1] := expand( subs(Phi(r)=-1/2*ln(r)+1/2*ln(r-2*M),g_matrix[1,1]) ):# coefficient of d(t)^2 in interval
sch_compts[2,2] := expand( subs(Lambda(r)=-ln(1-2*M/r)/2,g_matrix[2,2]) ):# coefficient of d(r)^2 in interval
sch_compts[3,3]:=g_matrix[3,3]:# coefficient of d(theta)^2 in interval
sch_compts[4,4]:=g_matrix[4,4]:# coefficient of d(phi)^2 in interval
sch := create([-1,-1], eval(sch_compts));# Schwarzschild metric

Now we consider the star, which is composed of an incompressible substance r = r0 = const (later we will use also the following approximation: p =( g -1) r, where g =1 (dust), 4/3 (noncoherent radiation), 2 (hard matter) ).

Then the mass is

> m(r) = int(4*Pi*rho0*r^2,r);
M = subs(r=R,rhs(%));# full mass

and the pressure is

> diff(p(r),r) = subs( {rho(r)=rho0,m(r)=4/3*Pi*rho0*r^3},\
-(4*Pi*r^3*p(r)+m(r))*(rho(r)+p(r))/(r*(r-2*m(r))) );# Oppenheimer-Volkoff equation
eq5 := dsolve(%,p(r));

> solve(\
simplify(subs(r=R,rhs(eq5[1])))=0,\
exp(-16*_C1*Pi*rho0));#boundary condition
solve(\
simplify(subs(r=R,rhs(eq5[2])))=0,\
exp(-16*_C1*Pi*rho0));#boundary condition

> sol1 := simplify(subs(exp(-16*_C1*Pi*rho0) = -3+8*Pi*rho0*R^2,rhs(eq5[1])));
sol2 := simplify(subs(exp(-16*_C1*Pi*rho0) = -3+8*Pi*rho0*R^2,rhs(eq5[2])));

In the center of star we have

> simplify(subs(r=0,sol1));
simplify(subs(r=0,sol2));

The pressure is infinity when the radius is equal to

> R_crit1 = simplify(solve(expand(denom(%)/12)=0,R)[1],radical,symbolic);
R_crit2 = simplify(solve(expand(denom(%%)/12)=0,R)[2],radical,symbolic);

> plot3d({subs(R=1,sol1),subs(R=1,sol2)},r=0..1,rho0=0..0.1,axes=boxed,title=`pressure`);#only positive solution has a physical sense

To imagine the space-time geometry it is necessary to consider an equator section of the star (q = p/2) at fixed time moment in 3-dimensional flat space. The corresponding procedure is named as " embedding ". From the metric tensor, the 2-dimensional line element is = + and for Euclidean 3-space we have = + + . We will investigate the 2-surface z=z ( r ). As dz= dr, one can obtain Euclidian line element: = [1+ ] +

> subs(theta=Pi/2,get_compts(sch));#Schwarzschild metric
dr^2*%[2,2] + dphi^2*%[3,3] = (1+diff(z(r),r)^2)*dr^2 +r^2*dphi^2;#equality of intervals of flat and embedded spaces
diff(z(r),r) = solve(%,diff(z(r),r))[1];
dsolve(%,z(r));# embedding

Now let's take the penultimate equation and to express M through = const:

> rho_sol := solve(M = 4/3*Pi*rho0*R^3,rho0);# density from full mass
fun1 := Int(1/sqrt(r/((4/3)*Pi*rho0*r^3)-1),r);# inside space
fun2 := Int(1/sqrt(r/((4/3)*Pi*rho0)-1),r);# outside space

> fun3 := value(subs(rho0=rho_sol,fun1));# inside
fun4 := value(subs(rho0=rho_sol,fun2));# outside

The resulting embedding for equatorial and vertical sections is presented below (Newtonian case corresponds to horizontal surface, that is the asymptote for r -- > ). The outside space lies out of outside ring representing star border.

> fig1 := plot3d(subs(M=1,subs(R=2.66*M,subs(r=sqrt(x^2+y^2),fun3))),x=-5..5,y=-5..5,grid=[100,100],style=PATCHCONTOUR):# the inside of star [r in units of M]

fig2 := plot3d(subs(M=1,subs(R=2.66*M,subs(r=sqrt(x^2+y^2),fun4))),x=-5..5,y=-5..5,style=PATCHCONTOUR):# the outside of star

display(fig1,fig2,axes=boxed);

We can see, that the change of radial coordinate dr versus the variation of interval dl= ( dl is the length of segment of curve on the depicted surface) in vicinity of the star is smaller in compare with Newtonian case (plane z =0, where dr=dl ). The star "rolls" the space so that an observer moves away from the star, but the radial coordinate increases slowly. But what is a hole in the vicinity of the center of surface, which illustrates the outside space? What happens if the star radius becomes equal or less than the radius of this hole R= 2 M ? Such unique object is the so-called black hole (see below).

But before consideration of the main features of black hole let us consider the motion of probe particle in space with Schwarzschild metric. We will base our consideration on the relativistic low of motion: + = 0 ( p is 4-momentum, m is the rest mass). Then (if p =[ ], where L corresponds to angle momentum, l is the affine parameter, second term describes the radial velocity):

> -E^2/(1-2*M/r(tau)) + diff(r(tau),tau)^2/(1-2*M/r(tau)) + L^2/r(tau)^2+1 = 0;# tau=lambda*mu, E=E/mu, L=L/mu
eq6 := diff(r(tau),tau)^2 = expand(solve(%,diff(r(tau),tau)^2));# integral of motion
V := sqrt(-factor(op(2,rhs(eq6)) + op(3,rhs(eq6)) + op(4,rhs(eq6)) + op(5,rhs(eq6))));# effective potential

The effective potential determining the orbital motion is (we vary the angular momentum L and r ):

> subs({M=1,r(tau)=x,L=y},V):
plot({subs(y=3,%),subs(y=4,%),subs(y=6,%),1},x=0..30,axes=boxed,color=[red,green,blue,brown],title=`potentials`,view=0.6..1.4);#level E=1 corresponds to particle, which was initially motionless

We considered the motion in the relativistic potential in the first part, but only in the framework of linear approximation (weak field). As result of such motion the perihelion rotation and light ray deflection appear. But now, in the more general case, one can see a very interesting feature of the relativistic orbital motion: unlike Newtonian case there is the maximum with the following potential decrease as result of radial coordinate decrease. If > E > 1 the motion is infinite (it is analogue of the hyperbolical motion in Newtonian case). > E =1 corresponds to parabolic motion. When E lies in the potential hole or E = V , we have the finite motion. The motion with energy, which is equal to extremal values of potential, corresponds to circle orbits (the stable motion for potential minimum, unstable one for the maximum). The existence of the extremes is defined by expressions:

> numer( simplify( diff( subs(r(tau)=r,V), r) ) ) = 0:# zeros of potential's derivative
solve(%,r);

As consequence, the circle motion is possible if L > 2 M , when r > 3(2 M ). The minimal distance, which allows the circular unstable motion, is defined by

> limit(1/2*(L-sqrt(L^2-12*M^2))*L/M, L=infinity);

Lack of extremes (red curve) or transition to r < 3 M results in the falling motion (the similar case is E > ). This is the so-called gravitational capture and has no analogy in Newtonian mechanics of point-like masses. The finite orbital motion in the case of the right local minimum existence (blue curve) is the analog of one in the Newtonian case, but has no elliptical character (see first part).

As the particular case, let's consider the radial ( L =0) fall of particle. The proper time of falling particle is defined by next expression

> tau = Int( 1/sqrt(subs({L=0,r(tau)=r},rhs(eq6))),r);
tau = Int(1/(sqrt(2*M/r-2*M/R)),r);#R=2*M/(1-E^2) - apoastr, i.e. radius of zero velocity
limit( value( rhs(%) ),r=2*M ):
simplify(%);#So, this is convergent integral for r-->2*M

For the remote exterior observer we have = = = E

(here is the time coordinate relatively infinitely remote observer):

> r[o] := Int(1/(1-2*M/r),r) = int(1/(1-2*M/r),r);
limit( value( rhs(%) ),r=2*M );#this is divergent integral for r-->2*M

These equations will be utilized below.

Degeneracy stars and gravitational collapse

The previously obtained expressions for the time of radial fall have an interesting peculiarity if the radius of star is less or equal to 2 M ( in dimensional case, that is the so-called " gravitational radius " ). The finite proper time t of fall corresponds to the infinite "remote" time , when r -- >2 M. This means, that for the remote observer the fall does not finish never. It results from the relativist time's slowing-down. The consequence of this phenomena is the infinite red shift of "falling" photons because of the red shift in Schwarzschild metric is defined by the next frequencies relation: = = = (here Dt are the proper time intervals between light flash in different radial points). Simultaneously, the escape velocity is equal to that is the velocity of light for sphere with R = . So, we cannot receive any information from interior of black hole.

How can appear the object with the radius, which is less than ? In the beginning we consider the pressure free radial (i.e. angular part of metric is equal to 0) collapse of dust sphere with mass M .

> r := 'r':
E := 'E':
subs( r=r(t),get_compts(sch) ):
d(s)^2 = %[1,1]*d(t)^2 + %[2,2]*d(r)^2;#Schwarzschild metric
-d(tau)^2 = collect(subs( d(r)=diff(r(t),t)*d(t),rhs(%) ),d(t));#tau is the proper time for the observer on the surface of sphere
%/d(tau)^2;
subs({d(t)=E/(1-2*M/r(t)),d(tau)=1},%);#we used d(t)/d(tau) = E/(1-2*M/r)
pot_1 := factor( solve(%,(diff(r(t),t))^2) );#"potential"

> plot({subs({E=0.5,M=1,r(t)=r},pot_1),0*r},r=1.7..3,axes=boxed,title=`(dr/dt)^2 vs r`);

The collapse is the motion from the right point of zero velocity (this is the velocity for remote observer) to the left point of zero velocity. These points are

> solve(subs(r(t)=r,pot_1)=0,r);

that are the apoastr (see above) and the gravitational radius. Hence (from the value for apoastr and pot_1 ) = = = .

The slowing-down of the collapse for remote observer in vicinity of results from the time slowing-down (see above). But the collapsing observer has the proper velocity :

> pot_2 := simplify(pot_1*(E/(1-2*M/r(t)))^2);# d/d(t)=(d/d(tau))*(1-2*M/r)/E

> plot({subs({E=0.5,M=1,r(t)=r},pot_2),0*r},r=1.7..3,axes=boxed,title=`(dr/dtau)^2 vs r`);

As result, the collapsing surface crosses the gravitational radius at finite time moment and =1 (that is the velocity of light), when E = 1 ( R = ). Time of collapse for falling observer is

> Int(1/sqrt(R/r-1),r)/sqrt(1-E^2);#from pot_2, R is apoastr
simplify( value(%), radical, symbolic);

It is obvious, that the obtained integral is = = .

It is natural, that the real stars are not clouds of dust and the equilibrium is supported by nuclear reactions. But when the nuclear reactions in the star finish, the "fall" of the star's surface can be prevented only by pressure of electron or baryons. The equilibrium state corresponds to the minimum of the net-energy composed of gravitational energy and thermal kinetic energy of composition. Let consider the hydrogen ball. When the temperature tends to zero, the kinetic energy of electrons is not equal to zero due to quantum degeneracy. In this case one electron occupies the cell with volume ~ , where l = is the Compton wavelength ( is the electron's momentum). For non-relativistic electron:

> E[e] := p[e]^2/m[e];#kinetic energy of electron
E[k] := simplify( n[e]*R^3*subs( p[e]=h*n[e]^(1/3)/(2*pi),E[e] ) );# full kinetic energy, n[e] is the number density of electrons
E[k] := subs( n[e]=M/R^3/m[p],%);# m[p] is the mass of proton. We supposed m[e]<<m[p] and n[p]=n[e]
E[g] := -G*M^2/R:# gravitational energy
E := E[g] + E[k];# full energy

> pot := -1/R+1/R^2:# the dependence of energy on R
plot(pot,R=0.5..10, title=`energy of degenerated star`);

One can see that the dependence of energy on radius has the minimum and, as consequence, there is the equilibrium state.

> solve( diff(E, R) = 0, R);# minimum of energy corresponding to equilibrium state

So, we have the equilibrium state with ~ and the star, which exists in such equilibrium state as result of the termination of nuclear reactions with the following radius decrease, is the white dwarf . The value of number density in this state is

> simplify( subs( R[min]=h^2/(G*M^(1/3)*m[e]*m[p]^(5/3)),n[e]=M/R[min]^3/m[p]) );#equilibrium number density

So, the mass's increase decreases the equilibrium radius (unlike model of liquid drop, see above) and to increase the number density (the maximal density of solids or liquid defined by atomic packing is ~20 g/ , for white dwarf this value is ~ g/ !). But, in compliance with principle of uncertainty, the last causes the electron momentum increase ( ~ h ). Our nonrelativistic approximation implies << or

> simplify( rhs(%)^(1/3)*h ) - c*m[e];#this must be large negative value
solve(%=0,M);

The last result gives the nonrelativistic criterion M << . Therefore in the massive white dwarf the electrons have to be the relativistic particles:

> E[e] := h*c*n[e]^(1/3);#E=p*c for relativistic particle
E[k] := simplify( n[e]*R^3*E[e] );
E[k] := subs( n[e]=M/R^3/m[p],%);
E[g] := -G*M^2/R:
E := E[g] + E[k];#this is dependence -a/R+b/R
#equilibrium state:
simplify( diff(E, R), radical ):
numer(%);#there is not dependence on R therefore there is not stable configuration with energy minimum
solve( %=0, M );

We obtained the estimation for so-called critical mass ~ =1.4 (so-called Chandrasekhar limit for white dwarfs). The smaller masses produce the white dwarf with non-relativistic electrons but larger masses causes collapse of star, which can not be prevented by pressure of degeneracy electrons. Such collapse can form a neutron star for 1.4 < M < 3 . Collapse for larger masses has to result in the black hole formation.