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

Cosmological models

Standard models

Now we turn to some specific isotropic and homogeneous models. The assumption of time-dependent character of the geometry results in the time-dependence of scaling factor a ( t ). Then 4-metric (spherical geometry) is

> coord := [t, chi, theta, phi]:
g_compts := array(symmetric,sparse,1..4,1..4):
g_compts[1,1] := -1:#dt^2
g_compts[2,2] := a(t)^2:#chi^2
g_compts[3,3]:=a(t)^2*sin(chi)^2:#dtheta^2
g_compts[4,4]:=a(t)^2*sin(chi)^2*sin(theta)^2:#dphi^2
g := create([-1,-1], eval(g_compts));
ginv := invert( g, 'detg' ):

The Einstein tensor is

> D1g := d1metric( g, coord ):
D2g := d2metric( D1g, coord ):
Cf1 := Christoffel1 ( D1g ):
RMN := Riemann( ginv, D2g, Cf1 ):
RICCI := Ricci( ginv, RMN ):
RS := Ricciscalar( ginv, RICCI ):
G := Einstein( g, RICCI, RS );

Let's the matter is defined by the energy-momentum tensor = ( p+ r ) + p :

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

Then the first Einstein equation - L = - 8 p (we add the so-called cosmological constant L, which can be considered as the energy density of vacuum):

> get_compts(G):
%[1,1]:
e1 := simplify(%):
get_compts(g):
%[1,1]:
e2 := simplify(%):
get_compts(T):
%[1,1]:
e3 := simplify(%):
E[1] := expand(e1/(-3)) = -8*Pi*expand(e3/(-3)) + expand(e2*Lambda/(-3));#first Einstein equation

Second equation is:

> get_compts(G):
%[2,2]:
e1 := simplify(%):
get_compts(g):
%[2,2]:
e2 := simplify(%):
get_compts(T):
%[2,2]:
e3 := simplify(%):
E[2] := expand(e1/a(t)^2) = -8*Pi*expand(e3/a(t)^2) + expand(e2*Lambda/a(t)^2);#second Einstein equation

Two other equations are tautological.

After elementary transformation we have for the second equation:

> expand(simplify(E[2]-E[1])/2);

Similarly, for the hyperbolical geometry we have:

> coord := [t, chi, theta, phi]:
g_compts := array(symmetric,sparse,1..4,1..4):
g_compts[1,1] := -1:#dt^2
g_compts[2,2] := a(t)^2:#chi^2
g_compts[3,3]:=a(t)^2*sinh(chi)^2:#dtheta^2
g_compts[4,4]:=a(t)^2*sinh(chi)^2*sin(theta)^2:#dphi^2
g := create([-1,-1], eval(g_compts));
ginv := invert( g, 'detg' ):

> D1g := d1metric( g, coord ):
D2g := d2metric( D1g, coord ):
Cf1 := Christoffel1 ( D1g ):
RMN := Riemann( ginv, D2g, Cf1 ):
RICCI := Ricci( ginv, RMN ):
RS := Ricciscalar( ginv, RICCI ):
G := Einstein( g, RICCI, RS );

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

> get_compts(G):
%[1,1]:
e1 := simplify(%):
get_compts(g):
%[1,1]:
e2 := simplify(%):
get_compts(T):
%[1,1]:
e3 := simplify(%):
E[1] := expand(e1/(-3)) = -8*Pi*expand(e3/(-3)) + expand(e2*Lambda/(-3));#first Einstein equation

> get_compts(G):
%[2,2]:
e1 := simplify(%):
get_compts(g):
%[2,2]:
e2 := simplify(%):
get_compts(T):
%[2,2]:
e3 := simplify(%):
E[2] := expand(e1/a(t)^2) = -8*Pi*expand(e3/a(t)^2) + expand(e2*Lambda/a(t)^2):#second Einstein equation
expand(simplify(E[2]-E[1])/2);

The second equation is identical with one in spherical geometry.

In general, the first and second differential equations for scaling factor in our ideal universe:

> basic_1 := (diff(a(t),t)/a(t))^2 = -K/a(t)^2 + Lambda/3 + 8*Pi*rho(t)/3;#first basic equation for homogeneous universe
basic_2 := diff(a(t),`$`(t,2))/a(t) = -4*Pi*(p(t)+rho(t)/3)+1/3*Lambda;# second basic equation for homogeneous universe

Here K =+1 for spherical, -1 for hyperbolical and 0 for flat geometries. The left-hand side of the first equation is the square of Hubble constant H . If we denote the recent value of this parameter as , then = 65 ( ps is the parallax second corresponding to distance about of 3.26 light-years or 3* km ) [ J. Garcia-Bellido, Astrophysics and cosmology, arXiv: hep-ph/0004188].

When K= L = 0, the universe is defined by the so-called critical density = =7.9* . Let's denote the recent relative values = of the matter, radiation, cosmological constant and curvature as:

= , = , = , = - .

Here is the recent scaling factor. At this moment the contribution of the radiation to density is negligible << . To rewrite the evolutionary equation in terms of we have to take into account the time dependence of these parameters: L does not depend on a ( t ) (in framework of standard model), the dependence of curvature ~ , density of matter ~ , but for radiation ~ . The last results from the change of density of quanta ~ and their energy ~ . The resulting equation can be written as:

> basic_3 := H(a)^2 = H[0]^2*(Omega[R]*a[0]^4/a^4+Omega[M]*a[0]^3/a^3+Omega[Lambda]+Omega[K]*a[0]^2/a^2);

As result of these definitions we have the cosmic sum rule :

1 = + + .

Let us transit to the new variables: y = , t = ( ), where is the present time moment. Then

> basic_4 := diff(y(tau),tau) = sqrt(1 + (1/y(tau)-1)*Omega[M] + (y(tau)^2-1)*Omega[Lambda]);#we neglect the radiation contribution and use the cosmic sum rule. d(tau)=H[0]*d(t) --> a[0]*H[0]*(d(y)/d(tau))/a = (d(a)/d(t))/a = H

Now we will consider some typical standard models. For review see, for example, P.J.E. Peebles, Principles of physical cosmology, Prinston Univ. Press, 1993 , or on-line surveys P. B. Pal, Determination of cosmological parameters, arXiv: hep-ph/9906447, J. R. Bondono, The cosmological models.

Einstein static

> rhs(basic_1) = 0:
rhs(basic_2) = 0:
allvalues( solve({%, %%},{a(t), Lambda}) );# static universe

At this moment for the matter p =0 and the contribution of the radiation is small, therefore:

=

This solution has a very simple form and requires K= +1 (the closed spherical universe) and positive L (repulsive action of vacuum), but does not satisfy the observations of red shift in the spectra of extragalactic sources. The last is the well-known fact demonstrating the expansion of universe (the increase of scaling factor).

Einstein-de Sitter

= 1 and, as consequence = = 0. This is the flat universe described by equation

> y(tau)^(1/2)*diff(y(tau),tau) = 1;
dsolve({%, y(0)=1}, y(tau)):
allvalues(%);

> plot(1/4*(12*tau+8)^(2/3), tau= -2/3..2, axes=BOXED, title=`Einstein-de Sitter universe`);

We have the infinite expansion from the initial singularity (Big Bang).

From the value of scaling factor it is possible to find the universe's age:

> 0 = (12*H[0]*(0-t[0]) + 8)^(2/3)/4:# y(0)=0
solve(%,t[0]);

> evalf((2/3)*3*10^19/65/60/60/24/365);#[yr]

The obtained universe's age does not satisfy the observations of star's age in older globular clusters.

de Sitter and anti-de Sitter

The empty universes ( =0) with R > 0 (de Sitter) and R < 0 (anti-de Sitter). Turning to basic_1 we have for de Sitter universe:

> assume(lambda,positive):
subs({rho(t)=0,K=1,Lambda=lambda*3},basic_1):#lambda=Lambda/3
expand(%*a(t)^2);
dsolve(%,a(t));

It is obvious, that the first two solutions don't satisfy the basic_2 . Another solutions give:

a ( t ) =

where C is the positive or negative constant defining the time moment corresponding to the minimal scaling factor.

> plot(subs(lambda=0.5,\
cosh(sqrt(lambda)*tau)/2/sqrt(lambda)),\
tau=-3..3,title=`de Sitter universe`,axes=boxed);

For anti-de Sitter model we have:

> assume(lambda,negative):
subs({rho(t)=0,K=-1,Lambda=lambda*3},basic_1):#lambda=Lambda/3
expand(%*a(t)^2);
dsolve(%,a(t));

De Sitter and anti-de Sitter universes play a crucial role in many modern issues in GR and cosmology, in particular, in inflationary scenarios of the beginning of universe evolution.

Closed Friedmann-Lemaitre

> 1, L = 0 and, as consequence < 0 ( K =+1). This is the closed spherical universe.

> y(tau)*(diff(y(tau),tau)^2+b) =Omega[M];# b=Omega[M]-1 > 0

> solve( %,diff(y(tau),tau) );

> d(t)=sqrt(y/b/(Omega[M]/b - y))*d(y);

Introducing a new variable f we have:

> sqrt(y/(Omega[M]/b - y)) = tan(phi);
y = solve(%, y);
y = convert(rhs(%),sincos);

But this is y = = . And

> defform(f=0,w1=0,w2=0,v=1,phi=0,y=0);
d(y) = d( (Omega[M]/b)*sin(phi)^2 );

Then our equation results in

> subs( y=Omega[M]*sin(phi)^2/b,sqrt(y/(b*(Omega[M]/b-y))) ):
simplify(%);

> d(t) = simplify( (Omega[M]/b)*tan(phi)*2*sin(phi)*cos(phi) )*d(phi)/sqrt(b);

Hence, the age of universe is:

> #initial conditions: y(0)=0 --> phi=0
int( Omega[M]*(1-cos(2*x))/(Omega[M]-1)^(3/2), x=0..phi):
sol_1 := t = simplify(%);#age of universe vs phi

This equation in the combination with y ( f ) is the parametric representation of cycloid:

> y = Omega[M]*(1-cos(2*phi))/(2*(Omega[M]-1)):
plot([subs(Omega[M]=5,rhs(sol_1)),subs(Omega[M]=5,rhs(%)),phi=0..Pi], axes=boxed, title=`closed Friedmann-Lemaitre universe`);

In this model the expansion from the singularity changes into contraction and terminates in singularity (Big Crunch)

The age of universe is:

> #initial conditions: y(t0)=1 --> phi=?
1 = sin(phi)^2*Omega[M]/(Omega[M]-1):
sol_2 := solve(%,phi);
subs(phi=sol_2[1],rhs(sol_1));#age of universe vs Omega[M]
plot(%, Omega[M]=1.01..10, axes=BOXED, title=`age of universe [ 1/H[0] ]`);

So, the increase of decreases the age of universe, which is less than in Einstein-de Sitter model. The last and the modern observations of microwave background suggesting ~0 do not agree with this model.

Open Friedmann-Lemaitre

< 1, L = 0 and, as consequence > 0 ( K = -1). This is an open hyperbolical spherical universe. From the above introduced equation we have

> y(tau)*(diff(y(tau),tau)^2 - a) = Omega[M];# a = -b = 1-Omega[M]>0

> solve( %,diff(y(tau),tau) );

> diff(y(tau), tau)=sqrt((Omega[M] + a*y(tau))/y(tau));

> dsolve(%, y(tau)):
sol_1 := simplify( subs(a=1-Omega[M],%),radical,symbolic);#implicit solution for y(tau)

> #definition of _C1
subs(tau=1,sol_1):
subs(y(1)=1,%):
solve(%,_C1):
simplify( subs(_C1=%,sol_1) ):
subs(y(tau)=y,%):
sol_2 := solve(%,tau);#implicit solution with defined _C1

This solution can be plotted in the following form:

> subs(Omega[M]=0.3,sol_2):
plot(%, y=0.01..2, axes=BOXED, colour=RED, title=`hyperbolical universe (time vs scaling factor)` );

The extreme case -- >0 (almost empty world)

> limit(sol_2,Omega[M]=0);

results in the linear low of universe expansion and gives the estimation for maximal age of universe in this model:

> evalf(3*10^19/65/60/60/24/365);#[yr]

The dependence on looks as:

> -(subs(y=0,sol_2)-subs(y=1,sol_2)):#age of universe [ 1/H[0] ]
plot(%, Omega[M]=0..1,axes=boxed, title=`age of universe [ 1/H[0] ]`);

It is natural, the growth of matter density decreases the universe age. If estimations of give 0.3 (that results in the age ~12 gigayears for considered model), then this model does not satisfies the observations. Moreover, it is possible, that the nonzero curvature is not appropriate (see below).

Expanding spherical and recollapsing hyperbolical universes

This model demonstrates the possibility of the infinite expansion of universe with spherical geometry in the presence of nonzero cosmological constant.

> with(DEtools):
DEplot(subs({Omega[M]=1,Omega[Lambda]=2.59},basic_4),[y(tau)],tau=-2.5..1,[[y(0)=1]],stepsize=0.01,axes=boxed,linecolor=red,arrows=none,title=`loitering expansion`);

We selected the parameters, which are close to ones for static universe ( =1) but the cosmological term slightly dominates > . There is the delay of expansion, which looks like transition to collapse, but the domination of causes the further expansion. This model can explain the large age of universe for large and to explain the misalignment in the "red shift-distance" distribution for distant quasars.

The next example demonstrate the collapsing behavior of open (hyperbolical) universe:

> p1 := DEplot(subs({Omega[M]=0.5,Omega[Lambda]=-1},basic_4),[y(tau)],tau=-0.65..0.7,[[y(0)=1]],stepsize=0.01,axes=boxed,linecolor=red,arrows=none):#expansion

> p2 := DEplot(subs({Omega[M]=0.5,Omega[Lambda]=-1},lhs(basic_4)=-rhs(basic_4)),[y(tau)],tau=0.71..1.98,[[y(0.71)=1.36]],stepsize=0.01,axes=boxed,linecolor=red,arrows=none):#we have to change the sign in right-hand side of basic_4 and to change the initial condition

> with(plots):
display(p1,p2,title=`recollapsing open universe`);

Bouncing model

Main feature of the above considered nonstatic models (except for de Sitter model) is the presence of singularity at the beginning of expansion and, it is possible, at the end of evolution. On the one hand, the singularity is not desirable, but, on the other hand, the standard model of Big Bang demands a high-density and hot beginning of expansion. The last corresponds to the conditions in the vicinity of the singularity. The choice of the parameters allows the behavior without singularity:

> p1 := DEplot(subs({Omega[M]=0.1,Omega[Lambda]=1.5},basic_4),[y(tau)],tau=-1.13..1,[[y(0)=1]],stepsize=0.01,axes=boxed,linecolor=red,arrows=none):#expansion
p2 := DEplot(subs({Omega[M]=0.1,Omega[Lambda]=1.5},lhs(basic_4)=-rhs(basic_4)),[y(tau)],tau=-3..-1.14,[[y(-3)=2.3]],stepsize=0.01,axes=boxed,linecolor=red,arrows=none):#contraction
display(p1,p2,title=`bouncing universe`);

Our universe (?)

The modern observations of microwave background suggest [ J. Garcia-Bellido, Astrophysics and cosmology, arXiv: hep-ph/0004188]: ~ 0 (flat universe), ~ 0.7, ~ 0.3. The corresponding behavior of scaling factor looks as:

> DEplot(subs({Omega[M]=0.3,Omega[Lambda]=0.7},basic_4),[y(tau)],tau=-.96...2,[[y(0)=1]],stepsize=0.01,axes=boxed,linecolor=red,arrows=none,title=`Our universe?`);

We are to note, that our comments about appropriate (or not appropriate) character of the considered models are not rigid. The picture of topologically nontrivial space-time can combine the universes with the very different local geometries and dynamics in the framework of a hypothetical Big Universe ("Tree of Universes").

Beginning

Bianchi models and Mixmaster universe

A good agreement of the isotropic homogeneous models with observations does not provide for their validity nearby singularity. In particular, in framework of analytical approach, we can refuse the isotropy of universe in the beginning of expansion [for review see Ya.B. Zel'dovich, I.D. Novikov, Relativistic Astrophysics, V. 2: The structure and evolution of the universe, The university of Chicago Press, Chicago (1983) ]. Let us consider the anisotropic homogeneous flat universe without rotation:

> coord := [t, x, y, z]:
g_compts := array(symmetric,sparse,1..4,1..4):# metric components
g_compts[1,1] := -1:
g_compts[2,2] := a(t)^2:
g_compts[3,3] := b(t)^2:
g_compts[4,4] := c(t)^2:

g := create([-1,-1], eval(g_compts));

ginv := invert( g, 'detg' ):

The Einstein tensor is:

> # 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 vacuum state the right-hand sides of Einstein equations are 0. Then

> E_eqn := get_compts(Estn):
e1 := numer( E_eqn[1,1] ) = 0;
e2 := expand( numer( E_eqn[2,2] )/a(t)^2 ) = 0;
e3 := expand( numer( E_eqn[3,3] )/b(t)^2 ) = 0;
e4 := expand( numer( E_eqn[4,4] )/c(t)^2 ) = 0;

We will search the power-behaved solutions of this system ( Kasner metric ):

> e5 := simplify( subs({a(t)=a_0*t^p_1,b(t)=b_0*t^p_2,c(t)=c_0*t^p_3},e1) );
e6 := simplify( subs({a(t)=a_0*t^p_1,b(t)=b_0*t^p_2,c(t)=c_0*t^p_3},e2) );
e7 := simplify( subs({a(t)=a_0*t^p_1,b(t)=b_0*t^p_2,c(t)=c_0*t^p_3},e3) );
e8 := simplify( subs({a(t)=a_0*t^p_1,b(t)=b_0*t^p_2,c(t)=c_0*t^p_3},e4) );

If t is not equal to zero, we have:

> e9 := factor( e5/t^(p_2-2+p_1+p_3) );
e10 := factor( e6/t^(p_3-2+p_2) );
e11 := factor( e7/t^(p_3-2+p_1) );
e12 := factor( e8/t^(p_2-2+p_1) );

> e13 := e9/(-b_0*a_0*c_0);
e14 := e10/(b_0*c_0);
e15 := e11/(a_0*c_0);
e16 := e12/(a_0*b_0);

The following manipulation gives a connection between parameters:

> ((e14 + e15 + e16) - e13)/2;

That is + + = p_1 + p_2 + p_3 (3) .

> collect(e15-e16,{p_2,p_3});
collect(e14-e15,{p_1,p_2});

Hence

p_3 ( p_3 + p_1 - 1) = p_2 ( p_2 + p_1 - 1),

p_2 ( p_3 + p_2 - 1) = p_1 ( p_3 + p_1 - 1).

The last expressions suggest the simple solution:

( p_3 + p_1 - 1) = p_2 , ( p_2 + p_1 - 1) = p_3 and

( p_3 + p_2 - 1) = p_1 , ( p_3 + p_1 - 1) = p_2 ,

or

( p_3 + p_1 - 1) = - p_2 , ( p_2 + p_1 - 1) = - p_3 and

( p_3 + p_2 - 1) = - p_1 , ( p_3 + p_1 - 1) = - p_2

The first results in

> p_3 + p_1 - 1 + p_2 + p_1 - 1 + p_3 + p_2 - 1 - (p_1+p_2+p_3);

The second produces

> (p_3 + p_1 - 1 + p_2 + p_1 - 1 + p_3 + p_2 - 1 + (p_1+p_2+p_3))/3;

These expressions in the combination with Eq. (3) suggest that there is only one independent parameter p :

> solve({p_1 + p_2 + p = 1, p_1^2 + p_2^2 + p^2 = 1},{p_1,p_2}):
sol_1 := allvalues(%);#first choice
solve({p_1 + p_2 + p = 3, p_1^2 + p_2^2 + p^2 = 3},{p_1,p_2}):
sol_2 := allvalues(%);#second choice

So, if the parameters p_1 , p_2 , p_3 are real, we have + + = p_1 + p_2 + p_3 =1 and

> sol_p_1 := factor(subs(sol_1,p_1));
sol_p_2 := factor(subs(sol_1,p_2));
sol_p_3 := 1 - %% - %;

The values of the analyzed parameters are shown in the next figure:

> plot({sol_p_1,sol_p_2,sol_p_3},p=-1/3..1, axes=boxed, title=`powers of t`);

Hence, there are two directions of expansion and one direction of contraction for t-- > (so-called "pancake" singularity ) or one direction of expansion and two direction of contraction for t-- >0 ( "sigar" singularity ).

Now let consider more complicated situation with homogeneous anisotropic metric.

We will use the tetrad notation of Einstein equations that allows to avoid the explicit manipulations with coordinates. In the synchronous frame we have for the homogeneous geometry (here Greek indexes change from 1 to 3) [ L.D. Landau, E.M. Lifshitz, The classical theory of fields, Pergamon Press, Oxford (1962) ]:

= - 1, = ,

where is the time dependent matrix, is the frame vector ( a changes from 1 to 3, that is the number of space triad). We will use the following representation for h:

> eta := array([[a(t)^2,0,0],[0,b(t)^2,0],[0,0,c(t)^2]]);#a, b, c are the time-dependent coefficients of anisotropic deformation
eta_inv := inverse(eta):
kappa1 := map(diff,eta,t):#kappa1[a,b]=diff(eta[a,b],t)
kappa2 := multiply( map(diff,eta,t),eta_inv );#kappa2[a]^b=diff(eta[a,b],t)*eta_inv, i.e. we raise the index by eta_inv

The first vacuum Einstein equation has a form

= - - = 0

> eq_1 := evalm( -trace(map(diff,kappa2,t))/2 - trace( multiply(kappa2,transpose(kappa2))/4) ) = 0;#first Einstein equation

The next Einstein equations demand the definition of Bianchi structured constants: = , = . Here = is the unit antisymmetric symbol, is the symmetric "tensor", which can be expressed through its principal values , is the vector. The Bianchi types are:

Type	a
I	0	0	0	0
II	0	1	0	0
VII	0	1	1	0
VI	0	1	-1	0
IX	0	1	1	1
VIII	0	1	1	-1
V	1	0	0	0
IV	1	0	0	1
VII	z	0	1	1
III (z=1), VI (z >< 1)	z	0	1	-1

We will analyze type IX model, where = = = [1] = [2]= [3] = 1. The curvature = {2 + + - ( ) + [ - 2 ]} is

> Cab := array([[1,0,0],[0,1,0],[0,0,1]]);#C^(ab)
C_ab := multiply(eta,eta,Cab);#C[ab]
Ca_b := multiply(eta,Cab);#C[a]^b
C_a_b := multiply(Cab,eta);#C^b[a];
delta := array([[1,0,0],[0,1,0],[0,0,1]]):
evalm( (4*multiply(Cab,C_ab) - trace(C_a_b)*evalm(C_a_b+Ca_b) + delta*(trace(C_a_b)^2-2*trace(multiply(Cab,C_ab))))/(2*det(eta)) ):#here we use the diagonal form of eta and symmetry of structured constants
P := map(simplify,%);

So, in the degenerate case a = b = c , t = const this model describes the closed spherical universe with a positive scalar curvature.

Hence we can obtained the next group of vacuum Einstein equations:

= - [( ) - ]= 0

> evalm( -map(diff,evalm(sqrt(det(eta))*kappa2),t)/(2*sqrt(det(eta)))):
evalm( map(simplify,%) - P );

The obtained expression gives three Einstein equations:

> eq_2 := Diff( (diff(a(t),t)*b(t)*c(t)),t )/(a(t)*b(t)*c(t)) = ( (b(t)^2-c(t)^2)^2 - a(t)^4 )/(2*a(t)^2*b(t)^2*c(t)^2);
eq_3 := Diff( (diff(b(t),t)*a(t)*c(t)),t )/(a(t)*b(t)*c(t)) = ( (a(t)^2-c(t)^2)^2 - b(t)^4 )/(2*a(t)^2*b(t)^2*c(t)^2);
eq_4 := Diff( (diff(c(t),t)*b(t)*a(t)),t )/(a(t)*b(t)*c(t)) = ( (a(t)^2-b(t)^2)^2 - c(t)^4 )/(2*a(t)^2*b(t)^2*c(t)^2);

The last group of Einstein equations results from

= - ( - ) = 0

It is obvious, that in our case this equation is identically zero (see expression for ). The substitutions of a = , b = , c = (do not miss these exponents with frame vectors!), dt = a b c ( d t ) and the expressions = = = and = = [ - ( ) ] result in the simplification of the equations form:

> eq_1 := diff((alpha(tau)+beta(tau)+gamma(tau)), tau$2)/2 = diff(alpha(tau), tau)*diff(beta(tau),tau) + diff(alpha(tau),tau)*diff(gamma(tau),tau) + diff(beta(tau),tau)*diff(gamma(tau),tau);
eq_2 := 2*diff(alpha(tau),tau$2) = (b(t)^2-c(t)^2)^2-a(t)^4;
eq_3 := 2*diff(beta(tau),tau$2) = (a(t)^2-c(t)^2)^2-b(t)^4;
eq_4 := 2*diff(gamma(tau),tau$2) = (a(t)^2-b(t)^2)^2-c(t)^4;

Comparison of obtained equations with above considered Kasner's type equations suggests that these systems are identical if the right-hand sides of eq_2 , eq_3 and eq_4 are equal to zero. Therefore we will consider the system's evolution as dynamics of Kasner's solution perturbation. In this case t = ln(t) + const ( dt = a b c dt = d t, as result dt = ) and = , = , = . We assume that the right-hand sides are close to zero, but contribution of -terms dominates. Hence:

> eq_5 := diff(alpha(tau),`$`(tau,2)) = -exp(4*alpha(tau))/2;
eq_6 := diff(beta(tau),`$`(tau,2)) = exp(4*alpha(tau))/2;
eq_7 := diff(gamma(tau),`$`(tau,2)) = exp(4*alpha(tau))/2;

The first equation is analog of equation describing the one-dimensional motion (a is the coordinate) in the presence of exponential barrier:

> p:= dsolve({eq_5,alpha(5)=-2,D(alpha)(5)=-0.5},alpha(tau),type=numeric):

> with(plots):
odeplot(p,[tau,diff(alpha(tau),tau)],-5..5,color=green):
plot(exp(4*(-0.5)*tau)/2,tau=-5..5,color=red):
display(%,%%,view=-1..1,axes=boxed, title=`reflection on barrier`);

We can see, that the"particle" with initial velocity = -0.5 (green curve) is reflected on barrier (red curve) with change of velocity sign. But eq_5 , eq_6 , eq_7 results in = const , = const therefore

= + - = + 2 , = + - = + 2 ,

t = a b c = ,

and the result is:

a = , b = , c =

The following procedure gives a numerical algorithm for iteration of powers of t .

> iterations := proc(iter,p)
p_1_old := evalhf( -1/2*p+1/2-1/2*sqrt(-(3*p+1)*(p-1)) );
p_2_old := evalhf( -1/2*p+1/2+1/2*sqrt(-(3*p+1)*(p-1)) );
p_3_old := evalhf( p ):

for m from 1 to iter do
if(p_1_old<0) then
pp := evalhf(1+2*p_1_old):
p_1_n := evalhf(-p_1_old/pp):
p_2_n := evalhf((p_2_old+2*p_1_old)/pp):
p_3_n := evalhf((p_3_old+2*p_1_old)/pp):
else fi:

if(p_2_old<0) then
pp := evalhf(1+2*p_2_old):
p_1_n := evalhf((p_1_old+2*p_2_old)/pp):
p_2_n := evalhf(-p_2_old/pp):
p_3_n := evalhf((p_3_old+2*p_2_old)/pp):
else fi:

if(p_3_old<0) then
pp := evalhf(1+2*p_3_old):
p_1_n := evalhf((p_1_old+2*p_3_old)/pp):
p_2_n := evalhf((p_2_old+2*p_3_old)/pp):
p_3_n := evalhf(-p_3_old/pp):
else fi:

p_1_old := p_1_n:
p_2_old := p_2_n:
p_3_old := p_3_n:
if m = iter then pts := [p_1_old, p_2_old,p_3_old] fi;
od:
pts
end:

> with(plots):
pointplot3d({seq(iterations(i,0.6), i=1 .. 8)},symbol=diamond,axes=BOXED,labels=[p_1,p_2,p_3]);
pointplot3d({seq(iterations(i,0.6), i=7 .. 12)},symbol=diamond,axes=BOXED,labels=[p_1,p_2,p_3]);

One can see, that the evolution has character of switching between different Kasner's epochs (all points lie on the section of sphere + + =1 by plane p_1 + p_2 + p_3 =1). The negative power of t increases the corresponding scaling factor if t-- >0. In the upper figure the periodical change of signs takes place between y - and z - directions. The universe oscillates in these directions and monotone contracts in x- direction ( t-- >0). Next the evolution changes (lower figure): x- z oscillations are accompanied by the contraction in y -direction. The whole picture is the switching between different oscillations as the result of logarithmical approach to singularity. It is important, that there is a strong dependence on initial conditions that can cause the chaotical scenario of oscillations:

> pointplot3d({seq(iterations(i,0.601), i=1..100)},symbol=diamond,axes=BOXED,labels=[p_1,p_2,p_3],title=`Kasner's epochs (I)`);
pointplot3d({seq(iterations(i,0.605), i=1..100)},symbol=diamond,axes=BOXED,labels=[p_1,p_2,p_3],title=`Kasner's epochs (II)`);

So, we can conclude that the deflection from isotropy changes the scenario of the universe's evolution in the vicinity of singularity, so that the dynamics of the geometry looks like chaotical behavior of the deterministic nonlinear systems (so-called Mixmaster universe , see J. Wainwright, C.F.R. Ellis, eds., Dynamical systems in cosmology, Cambridge University press, Cambridg (1997) ).

Inflation

The above-considered standard homogeneous models face the challenge of so-called horizon . The homogeneity of the universe results from the causal connection between its different regions. But the distance between these regions is defined by the time of expansion. Let's find the maximal distance of the light signal propagation in the expanding region. As = 0 the accompanying radial coordinate of horizon is

> r = Int(1/a(t),t=0..t0);# t0 is the age of the universe
#or
r = Int(1/y(tau),tau=-t0*H0..0)/(a0*H0);

For instance, in the Einstein-de Sitter model we have

> int( subs(y(tau)=1/4*(12*tau+8)^(2/3),1/y(tau)), tau=-2/3..0)/(a0*H0):# we use the above obtained expressions for this model
simplify( subs(a0=1/4*(12*0+8)^(2/3),%),radical );

and in the Friedmann-Lemaitre:

> int( subs(y(phi)=Omega[M]*(1-cos(2*phi))/(2*b),1/2*Omega[M]*(2-2*cos(2*phi))/(b^(3/2))/y(phi)),phi )/a0/H0;# we use the above obtained expressions for this model
simplify( subs(K=-1,subs(b=-K/a0^2/H0^2,%)),radical,symbolic );#b=Omega[M]-1=Omega[K] and K=-1 (spherical universe)

The last example is of interest. The increase of f from 0 to p/2 corresponds to expansion of universe up to its maximal radius (see above). At the moment of maximal expansion the observer can see the photons from antipole of universe that corresponds to formation of full causal connection:

> plot([2*phi-sin(2*phi),2*phi],phi=0..Pi, title=`age of universe and time of photon travel (a.u.)`);

This figure shows the age of Friedmann-Lemaitre universe in comparison with the time of photon propagation from outermost point of universe. For f < p/2 there are the regions, which don't connect with us. After this points there is the causal connection. At the moment of recollapse the photons complete revolutions around universe.

The different approaches were developed in order to overcome the problem of horizon (see above considered loitering and Mixmaster models). But the most popular models are based on the so-called inflation scenario , which can explain also the global flatness of universe (see, for example, A.D. Linde, Particle physics and inflationary cosmology, Harwood academic press, Switzerland (1990) , on-line reviews A. Linde, Lectures on inflationary cosmology, arXiv: hep-th/9410082, R.H. Brandenberger, Inflationary cosmology: progress and problems, arXiv: hep-ph/9910410).

For escape the horizon problem we have to suppose the accelerated expansion with > 0. As result of this condition, the originally causally connected regions will expand faster, than the horizon expands. So, we observe the light from the regions, which was connected at the early stage of universe's evolution. In the subsection "Standard models" we derived from Einstein equations the energy conservation low basic_2 :

This equation demonstrates that in the usual conditions ( p > 0) the pressure is a source of gravitation, but if p < - ( L =0) or L > 4pr ( p =0) the repulsion dominates and the expansion is accelerated. As the source of the repulsion, the scalar "inflaton" field f is considered.

Let us introduce the simple Lagrangian for this homogeneous field:

L =- - V ( f ),

where V is the potential energy. The procedure, which was considered in first section, allows to derive from this Lagrangian the field equation:

+ = 0.

In the case of flat RW metric this equation results in:

> coord := [t, r, theta, phi]:
g_compts := array(symmetric,sparse,1..4,1..4):
g_compts[1,1] := -1:
g_compts[2,2] := a(t)^2:
g_compts[3,3]:= a(t)^2*r^2:
g_compts[4,4]:= a(t)^2*r^2*sin(theta)^2:
g := create([-1,-1], eval(g_compts)):#definition of flat RW metric
d1g:= d1metric(g, coord):
g_inverse:= invert( g, 'detg'):
Cf1:= Christoffel1 ( d1g ):
Cf2:= Christoffel2 ( g_inverse, Cf1 ):
det_scal := simplify( sqrt( -det( get_compts(g) ) ),radical,symbolic ):
g_det_sq := create([], det_scal):#sqrt(-g)
field := create([], phi(t)):#field
#calculation of first term in the field equation
cov_diff( field, coord, Cf2 ):
T := prod(g_inverse,%):
contract(T,[2,3]):
Tt := prod(g_det_sq,%):
Ttt := cov_diff( Tt, coord, Cf2 ):
get_compts(%)[1,1]/det_scal:
expand(-%):
field_eq := % + diff(V(phi),phi) = 0;

Now let's define the energy-momentum tensor for first Einstein equation:

= L - ,

where = . As consequence, r = [- + V (f )] + = + V ( f ) , and (for i , j = 1,2,3) = ( - V (f )) --> p = - V (f ). Last expression emplies that the necessary condition p < - may result from large V .

The dynamics of our system is guided by coupled system of equations:

> basic_inf_1 := K/(a(t)^2)+diff(a(t),t)^2/(a(t)^2) = 8/3*Pi*(1/2*diff(phi(t),t)^2+V(phi));#see above E[1]
basic_inf_2 := field_eq;

The meaning of K is considered earlier, here we will suppose K =0 (local flatness). In the slow-rolling approximation , -- >0 and for potential V (f)= we have:

> basic_inf_1_n := subs(K=0,K/(a(t)^2)+diff(a(t),t)^2/(a(t)^2)) = 8/3*Pi*m^2*phi(t)^2/2;#see above E[1]
basic_inf_2_n := op(1,lhs(field_eq))+
subs(phi=phi(t),
expand(subs(V(phi)=m^2*phi^2/2,op(3,lhs(field_eq)))))=0;

> dsolve({basic_inf_1_n,basic_inf_2_n,phi(0)=phi0,a(0)=a0},{phi(t),a(t)}):
expand(%);

One can see, that for >>1 there is an quasi-exponential expansion (inflation) of universe:

a ( t )= , where H = 2 m

When t ~ the regime of slow rolling ends and the potential energy of inflaton field converts into kinetic form that causes the avalanche-like creation of other particles from vacuum (so-called " reheating "). Only at this moment the scenario become similar to Big Bang picture. What is the size of universe after inflation?

> a(t)=subs( {H=2*sqrt(Pi/3)*phi[0]*m,t=2*sqrt(3*Pi)/m*phi[0]},a[0]*exp(H*t) );

If at the beginning of inflation the energy and size of universe are defined by Plank density and length, then ~1, cm . The necessary value of m is constrained from the value of fluctuation of background radiation ( ). As result a ~ = cm , that is much larger of the observable universe. This is the reason why observable part of universe is flat, homogeneous and isotropic.

Let's return to basic_inf_2 . The reheating results from oscillation of inflaton field around potential minimum that creates new boson particles. The last damps the oscillations that can be taking into consideration by introducing of dumping term in basic_inf_2 :

> reh_inf := 3*diff(phi(t),t)*H+Gamma*diff(phi(t),t)+diff(phi(t),`$`(t,2)) = -m^2*phi(t);

Here H is the Hubble constant, G is the decay rate, which is defined by coupling with material boson field. When H > G, there are the coherent oscillations:

> dsolve(subs(Gamma=0,reh_inf),phi(t));
diff(%,t);

The last expression demonstrates the decrease of energy density via time increase ( r ~ ) and, as result, the decrease of H (see basic_inf_1 ). This leads to H < G, that denotes the reheating beginning.

In the conclusion of this short and elementary description it should be noted, that the modern inflation scenarios modify the standard model essentially. The Universe here looks like fractal tree of bubbles-universes, inflating and recollapsing foams, topological defects etc. The main features of this picture wait for thorough investigations.

Conclusion

Our consideration was brief and I omitted some important topics like perturbation theory for black holes and cosmological models, Penrose diagrams and conformal mappings for investigation of causal structure, gravitational waves etc. But I suppose to advance this worksheet hereafter.

LaTeX - version of this worksheet can be found on arXiv.org.

2000-2001�Kalashnikov