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 celestial mechanics in weak gravitational field

Introduction

The first results in the GR-theory were obtained without exact knowledge of the field equations (Einstein's equations for space-time geometry). The leading idea was the equivalence principle based on the equality of inertial and gravitational masses. We will demonstrate here that the natural consequences of this principle are the Schwarzschild metric and the basic experimental effects of GR-theory in weak gravitational field, i. e. planet's orbit precession and light ray deflection [see Sommerfeld A. Lectures on theoretical physics, v. 3, Academic Press, NY, 313-315 (1952) ].

Schwarzschild metric

Let us consider the centrally symmetric gravitational field, which is produced by mass M . The small cell falls along x -axis on the central mass. In the agreement with the equivalence principle, the uniformly accelerated motion locally compensates the gravitational force hence there is no a gravitational field in the free falling system . This results in the locally Lorenzian metric with linear element ( c is the velocity of light):

d = d + d + d - d

The velocity v and radial coordinate r are measured in the spherical system K, which are connected with central mass. It is natural, the observer in this motionless system "feels" the gravitational field. Since the first system moves relatively second one there are the following relations between coordinates:

d = (b = )

d = dt

d =rdq

d =r df

The first and second relations are the Lorentzian length shortening and time slowing down in the moving system. As result, from K looks as:

d = d + ( d + d ) - ( 1 - ) d

The sense of the additional terms in metric has to connect with the characteristics of gravitational field. What is the energy of in K ? If the mass of is m , and is the rest mass, the sum of kinetic and potential energies is:

> restart:
with(plots):
(m - m0)*c^2 - G*M*m/r=0;# energy conservation law (we suppose that the Newtonian law of gravitation is correct in the first approximation), G is the gravitational constant
%/(m*c^2):
subs( m=m0/sqrt(1-beta^2),% ):# relativistic mass
expand(%);
solve( %,sqrt(1-beta^2) ):
sqrt(1-beta^2) = expand(%);
1-beta^2 = taylor((1-subs(op(2,%)=alpha/r,rhs(%)))^2,alpha=0,2);# alpha=G*M/c^2, we use the first-order approximation on alpha

Warning, the name changecoords has been redefined

The last results in:

d = + ( d + d ) - (1 - ) d

This linear element describes the so-called Schwarzschild metric . In the first-order approximation:

> taylor( d(r)^2/(1-2*alpha/r)+r^2*(d(theta)^2+sin(theta)^2*d(phi)^2)-c^2*(1-2*alpha/r)*d(t)^2, alpha=0, 2 ):
convert( % , polynom ):
metric := d(s)^2 = collect( % , {d(r)^2, d(t)^2} );

Equations of motion

Let us consider a motion of the small unit mass in Newtonian gravitational potential F = - . Lagrangian describing the motion in centrally symmetric field is:

> L := (diff(r(t),t)^2 + r(t)^2*diff(theta(t),t)^2)/2 + G*M/r(t);

The transition to Schwarzschild metric transforms the time and radial differentials of coordinates (see, for example, [ Kai-Chia Cheng. A simple calculation of the perihelion Mercury from the principle of equivalence. Nature (Lond.), v. 155, 574 (1945); N.T. Roseveare, Mercury's perihelion from Le Verrier to Einstein, Clarendon Press - Oxford, 1982 ]):

dr --> (1 + ) dr and dt --> dt (it should be noted, that the replacement will be performed for differentials, not coordinates, and we use the weak field approximation for square-rooting):

> L_n := collect(\
expand(\
subs(\
{diff(r(t),t)=gamma(r(t))^2*diff(r(t),t),\ diff(theta(t),t)=gamma(r(t))*diff(theta(t),t)}, L ) ),\
{diff(r(t),t)^2, diff(theta(t),t)^2});# modified Lagrangian, gamma(r) = 1+alpha/r(t)

Next step for the obtaining of the equations of motion from Lagrangian is the calculation of the force and the momentum , where y = :

> e1 := Diff(Lagrangian(r, Diff(r,t)), r) = diff(subs(r(t) = r, L_n),r);# first component of force
e2 := Diff(Lagrangian(r, Diff(r,t)), Diff(r,t)) = subs( x=diff(r(t),t), diff(subs(diff(r(t),t)=x, L_n), x));# first component of momentum
e3 := Diff(Lagrangian(theta, Diff(theta,t)), theta) = diff(subs(theta(t) = theta, L_n), theta);# second component of force
e4 := Diff(Lagrangian(theta, Diff(theta,t)), Diff(theta,t)) = subs(y=diff(theta(t),t), diff(subs(diff(theta(t),t)=y, L_n), y));# second component of momentum

The equations of motion are the so-called Euler-Lagrange equations - = 0 (in fact, these equations are the second Newton's law and result from the law of least action). Now let us write the equations of motion in angular coordinates. Since e3 = 0 due to an equality to zero of the mixed derivative, we have from e4 the equation of motion - = 0 ( = , = ) in the form:

> Eu_Lagr_1 := Diff(rhs(e4),t) = 0;

Hence

= H = const

> sol_1 := Diff(theta(t),t) = solve( gamma^2*diff(theta(t),t)/u(theta)^2 = H, diff(theta(t),t) );# u=1/r is the new variable

The introduced replacement u= leads to the next relations:

> Diff(r(t), t) = diff(1/u(t),t);
Diff(r(t), t) = diff(1/u(theta),theta)*diff(theta(t),t);# change of variables
sol_2 := Diff(r(t), t) = subs(diff(theta(t),t) = rhs(sol_1), rhs(%));# the result of substitution of above obtained Euler-Lagrange equation

The last result will be used for the manipulations with second Euler-Lagrange equation - = 0. We have for the right-hand side of e2:

> Diff( subs( {diff(r(t),t)=rhs(%),gamma(r(t))=gamma}, rhs(e2)), t );

and can rewrite this expression:

> -2*gamma*H*diff(1+alpha/r(t),t)*diff(u(theta),theta) - H*gamma^2*diff(u(theta),theta$2)*diff(theta(t),t);# from the previous expression, definition of gamma and definition of derivative of product
first_term := subs({r(t)=1/u(theta), diff(theta(t),t)=rhs(sol_1)},\
subs(diff(r(t),t) = rhs(sol_2),%));# this is a first term in second Euler-Lagrange equation

e1 results in:

> 2*gamma(r)^3*diff(r(t),t)^2*diff(gamma(r),r) + \
r*gamma(r)^2*diff(theta(t),t)^2 + \
r^2*gamma(r)*diff(theta(t),t)^2*diff(gamma(r),r) - \
G*M/(r^2);# This is e1

subs(\
{diff(r(t),t) = rhs(sol_2), diff(gamma(r),r) = diff(1+alpha/r,r), diff(theta(t),t) = rhs(sol_1)},\
%):# we used the expressions for diff(r(t),t), gamma(r) and the first equation of motion

subs(gamma(r) = gamma, %):

second_term := subs(r = 1/u(theta), %);

And finally:

> Eu_Lagr_2 := expand( simplify(first_term - second_term)/u(theta)^2/H^2);

In the first-order approximation:

> gamma^n = taylor((1+alpha*u)^n, alpha=0,2);

So

> taylor( subs(gamma = 1+alpha*u(theta),Eu_Lagr_2), alpha=0,2 ):
basic_equation := convert(%, polynom) = 0;

Light ray deflection

In the beginning the obtained equation will be used for the search of the light ray deflection in the vicinity of a star. The fundamental consideration has to be based on the condition of null geodesic line d =0 for light, but we simplify a problem and consider the trajectory of the particle moving from the infinity. In this case H= .

> eq_def := subs(G*M/(H^2)=0,basic_equation);

The free propagation ( a=0) results in:

> subs(alpha=0, lhs(eq_def)) = 0;
sol := dsolve({%, u(0) = 1/R, D(u)(0) = 0}, u(theta));# theta is measured from perihelion, where r = R

That is r= . The last expression corresponds to the straight ray passing through point q=0, r = R . To find the corrected solution in the gravitational field let substitute the obtained solution into eliminated term in eq_def :

> -op(1,lhs(eq_def)) - op(2,lhs(eq_def)) = subs(u(theta)=rhs(sol),op(3,lhs(eq_def)));
sol := dsolve({%, u(0) = 1/R, D(u)(0) = 0}, u(theta));#corrected solution in the presence of field

The equation for asymptote describes the observed ray. Then angle of ray deflection is (symmetrical deflection before and after perihelion). Hence

> simplify( 2*subs({theta=Pi/2, alpha=G*M/c^2},rhs(sol))*R );

This is correct expression for the light ray deflection in the gravitational field. For sun we have:

> subs({kappa=0.74e-28, M=2e33, R=6.96e10},4*kappa*M/R/4.848e-6);
# in ["],where kappa = G/c^2 [cm/g], 4.848e-6 rad corresponds to 1"

The ray trajectory within Pluto's orbit distance is presented below:

> K := (theta, alpha, R) -> 1 / (-(alpha-R)*cos(theta)/(R^2)-1/2*alpha*cos(2*theta)/(R^2)+3/2*alpha/(R^2)):
S := theta -> K(theta, 2.125e-6, 1):#deflected ray
SS := theta -> 1/cos(theta):#ray without deflection
p1 := polarplot([S,theta->theta,Pi/2..-Pi/2],axes=boxed):
p2 := polarplot([SS,theta->theta,Pi/2..-Pi/2],axes=boxed,color=black):
display(p1,p2,view=-10700..10700,title=`deflection of light ray`);#distance of propagation corresponds to 100 AU=1.5e10 km, distance is normalyzed to Sun radius

Planet's perihelion motion

Now we return to basic_equation . Without relativistic correction to the metric ( a =0) the solution is

> sol := dsolve({subs({G*M/(H^2)=k,alpha=0},basic_equation),u(0)=1/R,D(u)(0)=0},u(theta));

This equation describes an elliptical orbit:

u=k ( 1 + e ),

where e = ( - 1) is the eccentricity. For Mercury k = 0.01, e = 0.2056, R = 83.3

> K := (theta, k, R) -> 1/(k-(k*R-1)*cos(theta)/R):
S := theta -> K(theta, .01, 83.3):
polarplot([S,theta->theta,0..2*Pi],title=`Orbit of Mercury`);

The correction to this expression results from the substitution of obtained solution in basic_equation .

> subs({G*M/(H^2)=k,3*u(theta)^2*alpha=3*alpha*(k*(1+e*cos(theta)))^2},\
basic_equation);
dsolve({%,u(0)=1/R,D(u)(0)=0},u(theta));

Now it is possible to plot the corrected orbit (we choose the exaggerated parameters for demonstration of orbit rotation):

> K := (theta, k, R, alpha, e) -> 1 / (3*alpha*k^2*e*cos(theta)+3*sin(theta)*alpha*k^2*e*theta+3/2*alpha*k^2*e^2-1/2*alpha*k^2*e^2*cos(2*theta)+k+3*alpha*k^2-(3*alpha*k^2*e*R+alpha*k^2*e^2*R+k*R+3*alpha*k^2*R-1)*cos(theta)/R):
S := theta -> K(theta, .42, 1.5, 0.01, 0.6):
polarplot([S,theta->theta,0..4.8*Pi],title=`rotation of orbit`);

Now we try to estimate the perihelion shift due to orbit rotation. Let's suppose that the searched solution of basic_equation differs from one in the plane space-time only due to ellipse rotation. The parameter describing this rotation is w: u ( q ) = k (1 + e ).

> subs({G*M/(H^2)=k,u(theta)=k*(1+e*cos(theta-omega*theta))},\
basic_equation):#substitution of approximate solution
simplify(%):
lhs(%):
collect(%, cos(-theta+omega*theta)):
coeff(%, cos(-theta+omega*theta)):#the coefficient before this term gives algebraic equation for omega
subs(omega^2=0,%):#omega is the small value and we don't consider the quadratic term
solve(% = 0, omega):
subs({alpha = G*M/c^2, k = 1/R/(1+e)},2*Pi*%);#result is expressed through the minimal distance R between planet and sun; 2*Pi corresponds to the transition to circle frequency of rotation
subs(R=a*(1-e),%);#result expressed through larger semiaxis a of an ellipse

Hence for Mercury we have the perihelion shift during 100 years:

> subs({kappa=0.74e-28,M=2e33,a=57.9e11,e=0.2056},\
6*Pi*kappa*M/a/(1-e^2)*(100*365.26/87.97)/4.848e-6):#here we took into account the periods of Earth's and Mercury's rotations
evalf(%);#["]

Now we will consider the basic_equation in detail [ J. L. Synge, Relativity: the general theory, Amsterdam (1960) ]. The implicit solutions of this equation are:

> dsolve( subs(G*M/H^2=k, basic_equation), u(theta));

Hence, these implicit solutions result from the following equation (b is the constant depending on the initial conditions):

> diff( u(theta), theta)^2 = 2*M*u(theta)^3 - u(theta)^2 + 2*u(theta)*M/H^2 + beta;# we use the units, where c=1, G=1 (see definition of the geometrical units in the next part)
f := rhs(%):
f = 2*M*( u(theta) - u1 )*( u(theta) - u2 )*( u(theta) - u3 );

Here u1 , u2 , u3 are the roots of cubic equation, which describes the "potential" defining the orbital motion:

> fun := rhs(%);

The dependence of this function on u leads to the different types of motion. The confined motion corresponds to an elliptical orbit

> plot(subs({u3=2, u2=1, u1=0.5, M=1/2},fun),u=0.4..1.1, title=`elliptical motion`, axes=boxed, view=0..0.08);

The next situation with u-- >0 ( r-- > ) corresponds to an infinite motion:

> plot(subs({u3=2, u2=1, u1=-0.1, M=1/2},fun),u=0..1.1, title=`hyperbolical motion`, axes=boxed, view=0..0.6);

Now we return to the right-hand side of the modified basic equation.

> f;

In this expression, one can eliminate the second term by the substitution u (q )= y ( q ) + :

> collect( subs(u(theta) = (2/M)^(1/3)*y(theta) + 1/(6*M), f),y(theta) );

This substitution reduced our equation to canonical form for the Weierstrass P function [ E. T. Whittaker, G. N. Watson, A course of modern analysis, Cambridge (1927) ]:

> g2 = simplify(coeff(%, y(theta)));
g3 = -simplify(coeff(%%, y(theta), 0));
diff( y(theta), theta)^2 = 4*y(theta)^3 - g2*y(theta) - g3;

that results in:

> y(theta) = WeierstrassP(theta, g2, g3);

In the general case, the potential in the form of three-order polynomial produces the solution in the form of Jacobi sn-function:

> Orbit := proc(f, x)
print(`Equation in the form: u'(theta)^2 = a[0]*u^3+a[1]*u^2+a[2]*u+a[3]`):
degree(f,x):
if(% = 3) then
a[0] := coeff(f, x^3):# coefficients of polynomial
a[1] := coeff(f, x^2):
a[2] := coeff(f, x):
a[3] := coeff(f, x, 0):
sol := solve(f = 0, x):
print(`Roots of polynomial u[1] < u[2] < u[3]:`):
print(sol[1], sol[2], sol[3]):
solution := u[1] + (u[2]-u[1])*JacobiSN(theta*sqrt( 2*M*(u[3]-u[1]) )/2 + delta,sqrt((u[2]-u[1])/(u[3]-u[1])))^2:
print(`Result through Jacobi sn - function`):
print(u(theta) = solution):
else
print(`The polynomial degree is not 3`)
fi
end:

> Orbit(subs(u(theta)=x,f),x);

As = and = are the small values for the planets ( and are the perihelion and aphelion points) and + + = , we have 2 M =1 and

> u(theta) - u[1] = (u[2] - u[1])*JacobiSN(1/2*theta+delta,0)^2;

that is the equation of orbital motion (see above) with e = . > 0 corresponds to the elliptical motion, < 0 corresponds to hyperbolical motion. The period of the orbital motion in the general case:

> (u[2]-u[1])/(u[3]-u[1]):
kernel := 2/sqrt((1-t^2)*(1-t^2*%)):#2 in the numerator corresponds to sn^2-period
int(kernel, t=0..1);

and can be found approximately for small and as result of expansion:

> series(2*EllipticK(x), x=0,4):
convert(%,polynom):
subs(x = sqrt( 2*M*(u[2]-u[1]) ), %);

From the expression for u ( q ) and obtained expression for the period of we have the change of angular coordinate over period:

> %/(1/2*sqrt( 2*M*(u[3]-u[1]) )):
simplify(%);

But = ~ 1 + . And the result is

> 2*Pi*(1 + M*(u[2]-u[1])/2)*(1 + M*(2*u[1]+u[2])):
series(%, u[1]=0,2):
convert(%,polynom):
series(%, u[2]=0,2):
convert(%,polynom):
expand(%):
expand(%-op(4,%)):
factor(%);

The deviation of the period from 2p causes the shift of the perihelion over one rotation of the planet around massive star.

> simplify(%-2*Pi):
subs({u[1]=1/r[1], u[2]=1/r[2]}, %):
subs({r[1]=a*(1+e), r[2]=a*(1-e)},G*%/c^2):# we returned the constants G and c
simplify(%);

This result coincides with the expression, which was obtained on the basis of approximate solution of basic_equation. From the Kepler's low we can express this result through orbital period:

> subs(M=4*Pi^2*a^3/T^2, %):# T is the orbital period
simplify(%);

Conclusion

So, we found the expression for the Schwarzschild metric from the equivalence principle without introducing of the Einstein's equations. On this basis and from Euler-Lagrange equations we obtained the main experimental consequences of GR-theory: the light ray deflection and planet's perihelion motion.