text:
# Introduction to relativistic astrophysics and cosmology through Maple

## 1 Introduction

A rapid progress of the observational astrophysics, which resulted from the active use of orbital telescopes, essentially intensifies the astrophysical researches at the last decade and allows to choose the more definite directions of further investigations. At the same time, the development of high-performance computers advances in the numerical astrophysics and cosmology. Against a background of these achievements, there is the renascence of analytical and semi-analytical approaches, which is induced by new generation of high-efficient computer algebra systems.

Here we present the pedagogical introduction to relativistic astrophysics and cosmology, which is based on computational and graphical resources of Maple 6. The pedagogical aims define the use only standard functions despite the fact that there are the powerful General Relativity (GR) oriented extensions like GRTensor . The knowledge of basics of GR and differential geometry is supposed. It should be noted, that our choice of metric signature (+2) governs the definitions of Lagrangians and energy-momentum tensors.

The computations in this worksheet take about of 6 min of CPU time (PIII-500) and 9 Mb of memory.

## 2 Relativistic celestial mechanics in weak gravitational field

### 2.1 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 ).

### 2.2 Schwarzschild metric

Let us consider the centrally symmetric gravitational field, which is produced by mass M. The small cell $`K_{\mathrm{}}`$ 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 $`K_{\mathrm{}}`$. This results in the locally Lorenzian metric with linear element (c is the velocity of light):

d $`s^2`$ = d $`x_{\mathrm{}}^{}{}_{}{}^{2}`$ \+ d $`y_{\mathrm{}}^{}{}_{}{}^{2}`$ \+ d $`z_{\mathrm{}}^{}{}_{}{}^{2}`$ \- $`c^2`$ d $`t_{\mathrm{}}^{}{}_{}{}^{2}`$

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 $`x_{\mathrm{}}`$ = $`\frac{\mathrm{𝑑𝑟}}{\sqrt{1\beta ^2}}`$  ( $`\beta `$ = $`\frac{v}{c}`$ )

d $`t_{\mathrm{}}`$ = $`\sqrt{1\beta ^2}`$ dt

d $`y_{\mathrm{}}`$ =rd $`\theta `$

d $`z_{\mathrm{}}`$ =r $`\mathrm{sin}(\theta )`$ d $`\varphi `$

The first and second relations are the Lorentzian length shortening and time slowing down in the moving system. As result, $`K_{\mathrm{}}`$ from K looks as:

d $`s^2`$ = $`(1\beta ^2)^{(1)}`$ d $`r^2`$ \+ $`r^2`$ (d $`\theta ^2`$ \+ $`\mathrm{sin}(\theta )^2`$ d $`\varphi ^2`$ ) \- $`c^2`$ (1 \- $`\beta ^2`$ )d $`t^2`$

The sense of the additional terms in metric has to connect with the characteristics of gravitational field. What is the energy of $`K_{\mathrm{}}`$ in K? If the mass of $`K_{\mathrm{}}`$ is m, and $`m_0`$ 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

$$(m\mathit{m0})c^2\frac{GMm}{r}=0$$

$$1\sqrt{1\beta ^2}\frac{GM}{c^2r}=0$$

$$\sqrt{1\beta ^2}=1\frac{GM}{c^2r}$$

$$1\beta ^2=12\frac{1}{r}\alpha +\mathrm{O}(\alpha ^2)$$

The last results in:

d $`s^2`$ = $`\frac{\mathrm{𝑑𝑟}^2}{1\frac{2\alpha }{r}}`$ \+ $`r^2`$ (d $`\theta ^2`$  \+ $`\mathrm{sin}(\theta )^2`$ d $`\varphi ^2`$ ) - $`c^2`$ (1 - $`\frac{2\alpha }{r}`$ )d $`t^2`$

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} );

$`\mathrm{𝑚𝑒𝑡𝑟𝑖𝑐}:=\mathrm{d}(s)^2=`$

$`(c^2+{\displaystyle \frac{2c^2\alpha }{r}})\mathrm{d}(t)^2+(1+{\displaystyle \frac{2\alpha }{r}})\mathrm{d}(r)^2+r^2(\mathrm{d}(\theta )^2+\mathrm{sin}(\theta )^2\mathrm{d}(\varphi )^2)`$

### 2.3 Equations of motion

Let us consider a motion of the small unit mass in Newtonian gravitational potential $`\mathrm{\Phi }`$ = \- $`\frac{GM}{r}`$ . 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);

$$L:=\frac{1}{2}(\frac{}{t}\mathrm{r}(t))^2+\frac{1}{2}\mathrm{r}(t)^2(\frac{}{t}\theta (t))^2+\frac{GM}{\mathrm{r}(t)}$$

The transition to Schwarzschild metric transforms the time and radial differentials of coordinates (see, for example, ): dr –$`>`$ (1 \+ $`\frac{\alpha }{r}`$ )dr and dt –$`>`$ $`(1+\frac{\alpha }{r})^{(1)}`$ 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)

$$L\mathrm{\_}n:=\frac{1}{2}\gamma (\mathrm{r}(t))^4(\frac{}{t}\mathrm{r}(t))^2+\frac{1}{2}\mathrm{r}(t)^2\gamma (\mathrm{r}(t))^2(\frac{}{t}\theta (t))^2+\frac{GM}{\mathrm{r}(t)}$$

Next step for the obtaining of the equations of motion from Lagrangian is the calculation of the force $`\frac{}{x}\mathrm{L}(x,y)`$ and the momentum $`\frac{}{y}\mathrm{L}(x,y)`$ , where y= $`\frac{}{t}x`$ :

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

$`\mathit{e1}:=\frac{}{r}\mathrm{Lagrangian}(r,\frac{}{t}r)=`$ $`2\gamma (r)^3(\frac{}{t}r)^2(\frac{}{r}\gamma (r))+r\gamma (r)^2(\frac{}{t}\theta (t))^2+r^2\gamma (r)(\frac{}{t}\theta (t))^2(\frac{}{r}\gamma (r))`$  $`{\displaystyle \frac{GM}{r^2}}`$

$$\mathit{e2}:=\mathrm{Diff}(\mathrm{Lagrangian}(r,\frac{}{t}r),\frac{}{t}r)=\gamma (\mathrm{r}(t))^4(\frac{}{t}\mathrm{r}(t))$$

$$\mathit{e3}:=\frac{}{\theta }\mathrm{Lagrangian}(\theta ,\frac{}{t}\theta )=\mathrm{r}(t)^2\gamma (\mathrm{r}(t))^2(\frac{}{t}\theta )(\frac{^2}{\theta t}\theta )$$

$$\mathit{e4}:=\mathrm{Diff}(\mathrm{Lagrangian}(\theta ,\frac{}{t}\theta ),\frac{}{t}\theta )=\mathrm{r}(t)^2\gamma (\mathrm{r}(t))^2(\frac{}{t}\theta (t))$$

The equations of motion are the so-called Euler-Lagrange equations
$`\frac{^2}{ty}\mathrm{L}(x,y)\frac{}{x}\mathrm{L}(x,y)`$ = 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 $`\frac{^2}{ty_2}\mathrm{L}(r,\theta ,y_1,y_2)`$  \- $`\frac{}{\theta }\mathrm{L}(r,\theta ,y_1,y_2)`$  = 0 ( $`y_1`$ = $`\frac{}{t}r`$ , $`y_2`$ = $`\frac{}{t}\theta `$ ) in the form:

Eu\_Lagr\_1 := Diff(rhs(e4),t) = 0;

$$\mathrm{𝐸𝑢}\mathrm{\_}\mathrm{𝐿𝑎𝑔𝑟}\mathrm{\_}\mathit{1}:=\frac{}{t}\mathrm{r}(t)^2\gamma (\mathrm{r}(t))^2(\frac{}{t}\theta (t))=0$$

Hence

$`\mathrm{r}(t)^2\gamma (\mathrm{r}(t))^2(\frac{}{t}\theta (t))=\text{H}=\text{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

$$\mathrm{𝑠𝑜𝑙}\mathrm{\_}\mathit{1}:=\frac{}{t}\theta (t)=\frac{H\mathrm{u}(\theta )^2}{\gamma ^2}$$

The introduced replacement u= $`\frac{1}{r}`$  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

$$\frac{}{t}\mathrm{r}(t)=\frac{\frac{}{t}\mathrm{u}(t)}{\mathrm{u}(t)^2}$$

$$\frac{}{t}\mathrm{r}(t)=\frac{(\frac{}{\theta }\mathrm{u}(\theta ))(\frac{}{t}\theta (t))}{\mathrm{u}(\theta )^2}$$

$$\mathrm{𝑠𝑜𝑙}\mathrm{\_}\mathit{2}:=\frac{}{t}\mathrm{r}(t)=\frac{(\frac{}{\theta }\mathrm{u}(\theta ))H}{\gamma ^2}$$

The last result will be used for the manipulations with second Euler-Lagrange equation $`\frac{^2}{ty_1}\mathrm{L}(r,\theta ,y_1,y_2)`$  \- $`\frac{}{r}\mathrm{L}(r,\theta ,y_1,y_2)`$  = 0. We have for the right-hand side of e2:

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

$$\frac{}{t}(\gamma ^2(\frac{}{\theta }\mathrm{u}(\theta ))H)$$

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

$$2\frac{\gamma H\alpha (\frac{}{t}\mathrm{r}(t))(\frac{}{\theta }\mathrm{u}(\theta ))}{\mathrm{r}(t)^2}H\gamma ^2(\frac{^2}{\theta ^2}\mathrm{u}(\theta ))(\frac{}{t}\theta (t))$$

$$\mathrm{𝑓𝑖𝑟𝑠𝑡}\mathrm{\_}\mathrm{𝑡𝑒𝑟𝑚}:=2\frac{H^2\alpha \mathrm{u}(\theta )^2(\frac{}{\theta }\mathrm{u}(\theta ))^2}{\gamma }H^2(\frac{^2}{\theta ^2}\mathrm{u}(\theta ))\mathrm{u}(\theta )^2$$

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), %);

$`2\gamma (r)^3(\frac{}{t}\mathrm{r}(t))^2(\frac{}{r}\gamma (r))+r\gamma (r)^2(\frac{}{t}\theta (t))^2+r^2\gamma (r)(\frac{}{t}\theta (t))^2(\frac{}{r}\gamma (r))`$ $`{\displaystyle \frac{GM}{r^2}}`$

$`\mathrm{𝑠𝑒𝑐𝑜𝑛𝑑}\mathrm{\_}\mathrm{𝑡𝑒𝑟𝑚}:=`$ $`2{\displaystyle \frac{H^2\alpha \mathrm{u}(\theta )^2(\frac{}{\theta }\mathrm{u}(\theta ))^2}{\gamma }}+{\displaystyle \frac{\mathrm{u}(\theta )^3H^2}{\gamma ^2}}{\displaystyle \frac{H^2\mathrm{u}(\theta )^4\alpha }{\gamma ^3}}GM\mathrm{u}(\theta )^2`$

And finally:

Eu\_Lagr\_2 := expand( simplify(first\_term -\
second\_term)/u(theta)^2/H^2);

$$\mathrm{𝐸𝑢}\mathrm{\_}\mathrm{𝐿𝑎𝑔𝑟}\mathrm{\_}\mathit{2}:=(\frac{^2}{\theta ^2}\mathrm{u}(\theta ))\frac{\mathrm{u}(\theta )}{\gamma ^2}+\frac{\mathrm{u}(\theta )^2\alpha }{\gamma ^3}+\frac{GM}{H^2}$$

In the first-order approximation:

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

$$\gamma ^n=1+nu\alpha +\mathrm{O}(\alpha ^2)$$

So

taylor( subs(gamma = 1+alpha\*u(theta),Eu\_Lagr\_2),\
alpha=0,2 ):
basic\_equation := convert(%, polynom) = 0;

$$\mathrm{𝑏𝑎𝑠𝑖𝑐}\mathrm{\_}\mathrm{𝑒𝑞𝑢𝑎𝑡𝑖𝑜𝑛}:=(\frac{^2}{\theta ^2}\mathrm{u}(\theta ))\mathrm{u}(\theta )+\frac{GM}{H^2}+3\mathrm{u}(\theta )^2\alpha =0$$

### 2.4 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 $`s^2`$ =0 for light, but we simplify a problem and consider the trajectory of the particle moving from the infinity. In this case H= $`\mathrm{}`$ .

eq\_def := subs(G\*M/(H^2)=0,basic\_equation);

$$\mathrm{𝑒𝑞}\mathrm{\_}\mathrm{𝑑𝑒𝑓}:=(\frac{^2}{\theta ^2}\mathrm{u}(\theta ))\mathrm{u}(\theta )+3\mathrm{u}(\theta )^2\alpha =0$$

The free propagation ( $`\alpha `$ =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

$$(\frac{^2}{\theta ^2}\mathrm{u}(\theta ))\mathrm{u}(\theta )=0$$

$$\mathrm{𝑠𝑜𝑙}:=\mathrm{u}(\theta )=\frac{\mathrm{cos}(\theta )}{R}$$

That is r= $`\frac{R}{\mathrm{cos}(\theta )}`$  . The last expression corresponds to the straight ray passing through point $`\theta `$ =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

$$(\frac{^2}{\theta ^2}\mathrm{u}(\theta ))+\mathrm{u}(\theta )=3\frac{\mathrm{cos}(\theta )^2\alpha }{R^2}$$

$$\mathrm{𝑠𝑜𝑙}:=\mathrm{u}(\theta )=\frac{(\alpha +R)\mathrm{cos}(\theta )}{R^2}\frac{1}{2}\frac{\alpha \mathrm{cos}(2\theta )}{R^2}+\frac{{\displaystyle \frac{3}{2}}\alpha }{R^2}$$

The equation for asymptote $`lim_{\theta (\frac{\pi }{2})}\mathrm{u}(\theta )`$ describes the observed ray. Then angle of ray deflection is $`\frac{2R}{r}`$ (symmetrical deflection before and after perihelion). Hence

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

$$4\frac{GM}{Rc^2}$$

This is a 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"

$$1.754485794$$

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

### 2.5 Planet’s perihelion motion

Now we return to basic\_equation. Without relativistic correction to the metric ( $`\alpha `$ =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));

$$\mathrm{𝑠𝑜𝑙}:=\mathrm{u}(\theta )=k\frac{(kR1)\mathrm{cos}(\theta )}{R}$$

This equation describes an elliptical orbit:

u=k(1+e\* $`\mathrm{cos}(\theta )`$ ),

where e = ( $`\frac{1}{kR}`$  \- 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));

$$(\frac{^2}{\theta ^2}\mathrm{u}(\theta ))\mathrm{u}(\theta )+k+3\alpha k^2(1+e\mathrm{cos}(\theta ))^2=0$$

$`\mathrm{u}(\theta )=`$ $`3\alpha k^2e\mathrm{cos}(\theta )+3\mathrm{sin}(\theta )\alpha k^2e\theta +{\displaystyle \frac{3}{2}}\alpha k^2e^2{\displaystyle \frac{1}{2}}\alpha k^2e^2\mathrm{cos}(2\theta )+k+3\alpha k^2`$  $`{\displaystyle \frac{(3\alpha k^2eR+\alpha k^2e^2R+kR+3\alpha k^2R1)\mathrm{cos}(\theta )}{R}}`$

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 $`\omega `$ : u( $`\theta `$ ) = k(1 + e $`\mathrm{cos}(\theta \omega \theta )`$ ).

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

$$6\frac{\pi GM}{c^2R(1+e)}$$

$$6\frac{\pi GM}{c^2a(1e)(1+e)}$$

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(%);#\["\]

$$43.08704513$$

Now we will consider the basic\_equation in detail . The implicit solutions of this equation are:

dsolve( subs(G\*M/H^2=k, basic\_equation), u(theta));

$`{\displaystyle _{}^{\mathrm{u}(\theta )}}{\displaystyle \frac{1}{\sqrt{\mathrm{\_}\mathit{C1}\mathrm{\_}a^2+2k\mathrm{\_}a+2\mathrm{\_}a^3\alpha }}}𝑑\mathrm{\_}a\theta \mathrm{\_}\mathit{C2}=0,`$ $`{\displaystyle _{}^{\mathrm{u}(\theta )}}{\displaystyle \frac{1}{\sqrt{\mathrm{\_}\mathit{C1}\mathrm{\_}a^2+2k\mathrm{\_}a+2\mathrm{\_}a^3\alpha }}}d\mathrm{\_}a\theta \mathrm{\_}\mathit{C2}=0`$

Hence, these implicit solutions result from the following equation ( $`\beta `$ 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 );

$$(\frac{}{\theta }\mathrm{u}(\theta ))^2=2M\mathrm{u}(\theta )^3\mathrm{u}(\theta )^2+\frac{2\mathrm{u}(\theta )M}{H^2}+\beta $$

$$2M\mathrm{u}(\theta )^3\mathrm{u}(\theta )^2+\frac{2\mathrm{u}(\theta )M}{H^2}+\beta =2M(\mathrm{u}(\theta )\mathit{u1})(\mathrm{u}(\theta )\mathit{u2})(\mathrm{u}(\theta )\mathit{u3})$$

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

fun := rhs(%);

$$\mathrm{𝑓𝑢𝑛}:=2M(\mathrm{u}(\theta )\mathit{u1})(\mathrm{u}(\theta )\mathit{u2})(\mathrm{u}(\theta )\mathit{u3})$$

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–$`>`$ $`\mathrm{}`$ ) 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;

$$2M\mathrm{u}(\theta )^3\mathrm{u}(\theta )^2+\frac{2\mathrm{u}(\theta )M}{H^2}+\beta $$

In this expression, one can eliminate the second term by the substitution u( $`\theta `$ )= y( $`\theta `$ ) $`(\frac{2}{M})^{(\frac{1}{3})}`$ \+ $`\frac{1}{6M}`$ :

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

$`4\mathrm{y}(\theta )^3+`$ $`\left({\displaystyle \frac{1}{6}}{\displaystyle \frac{2^{(1/3)}({\displaystyle \frac{1}{M}})^{(1/3)}}{M}}+{\displaystyle \frac{\mathrm{2\hspace{0.17em}2}^{(1/3)}({\displaystyle \frac{1}{M}})^{(1/3)}M}{H^2}}\right)\mathrm{y}(\theta ){\displaystyle \frac{1}{54}}{\displaystyle \frac{1}{M^2}}+\beta +{\displaystyle \frac{{\displaystyle \frac{1}{3}}}{H^2}}`$

This substitution reduced our equation to canonical form for the Weierstrass P function :

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

$$\mathit{g2}=\frac{1}{6}\frac{2^{(1/3)}({\displaystyle \frac{1}{M}})^{(1/3)}(H^2+12M^2)}{MH^2}$$

$$\mathit{g3}=\frac{1}{54}\frac{H^2+54\beta M^2H^2+18M^2}{M^2H^2}$$

$$(\frac{}{\theta }\mathrm{y}(\theta ))^2=4\mathrm{y}(\theta )^3\mathit{g2}\mathrm{y}(\theta )\mathit{g3}$$

that results in:

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

$$\mathrm{y}(\theta )=\mathrm{WeierstrassP}(\theta ,\mathit{g2},\mathit{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\*u^3+a\*u^2+a\*u+a‘):
degree(f,x):
if(% = 3) then
a := coeff(f, x^3):# coefficients of polynomial
a := coeff(f, x^2):
a := coeff(f, x):
a := coeff(f, x, 0):
sol := solve(f = 0, x):
print(‘Roots of polynomial u \< u \< u:‘):
print(sol, sol, sol):
solution := u + (u-u)\*JacobiSN(theta\*sqrt(\
2\*M\*(u-u) )/2 + delta,\
sqrt((u-u)/(u-u)))^2:
print(‘Result through Jacobi sn - function‘):
print(u(theta) = solution):
else
print(‘The polynomial degree is not 3‘)
fi
end:

$`\mathrm{𝐸𝑞𝑢𝑎𝑡𝑖𝑜𝑛}\mathrm{𝑖𝑛}\mathrm{𝑡ℎ𝑒}\mathrm{𝑓𝑜𝑟𝑚}:`$ $`u^{}(\mathrm{𝑡ℎ𝑒𝑡𝑎})\mathrm{^}\mathit{2}=a[\mathit{0}]u\mathrm{^}\mathit{3}+a[\mathit{1}]u\mathrm{^}\mathit{2}+a[\mathit{2}]u+a[\mathit{3}]`$

$$\mathrm{𝑅𝑜𝑜𝑡𝑠}\mathrm{𝑜𝑓}\mathrm{𝑝𝑜𝑙𝑦𝑛𝑜𝑚𝑖𝑎𝑙}u[\mathit{1}]<u[\mathit{2}]<u[\mathit{3}]:$$

$`{\displaystyle \frac{1}{6}}{\displaystyle \frac{\%1^{(1/3)}}{MH}}{\displaystyle \frac{1}{6}}{\displaystyle \frac{H^2+12M^2}{MH\%1^{(1/3)}}}+{\displaystyle \frac{{\displaystyle \frac{1}{6}}}{M}},`$ $`{\displaystyle \frac{1}{12}}{\displaystyle \frac{\%1^{(1/3)}}{MH}}+{\displaystyle \frac{{\displaystyle \frac{1}{12}}(H^2+12M^2)}{MH\%1^{(1/3)}}}+{\displaystyle \frac{{\displaystyle \frac{1}{6}}}{M}}+`$  $`{\displaystyle \frac{1}{2}}I\sqrt{3}({\displaystyle \frac{1}{6}}{\displaystyle \frac{\%1^{(1/3)}}{MH}}+{\displaystyle \frac{{\displaystyle \frac{1}{6}}(H^2+12M^2)}{MH\%1^{(1/3)}}}),`$  $`{\displaystyle \frac{1}{12}}{\displaystyle \frac{\%1^{(1/3)}}{MH}}+{\displaystyle \frac{{\displaystyle \frac{1}{12}}(H^2+12M^2)}{MH\%1^{(1/3)}}}+{\displaystyle \frac{{\displaystyle \frac{1}{6}}}{M}}`$  $`{\displaystyle \frac{1}{2}}I\sqrt{3}({\displaystyle \frac{1}{6}}{\displaystyle \frac{\%1^{(1/3)}}{MH}}+{\displaystyle \frac{{\displaystyle \frac{1}{6}}(H^2+12M^2)}{MH\%1^{(1/3)}}})`$  $`\%1:=H^354\beta M^2H^318HM^2`$  $`+6\sqrt{3}\sqrt{M^2H^2+16M^4H^6\beta +27H^6\beta ^2M^2+18H^4\beta M^2}M`$

$$\mathrm{𝑅𝑒𝑠𝑢𝑙𝑡}\mathrm{𝑡ℎ𝑟𝑜𝑢𝑔ℎ}\mathrm{𝐽𝑎𝑐𝑜𝑏𝑖}\mathrm{𝑠𝑛}\mathrm{𝑓𝑢𝑛𝑐𝑡𝑖𝑜𝑛}$$

$$\mathrm{u}(\theta )=u_1+(u_2u_1)\mathrm{JacobiSN}(\frac{1}{2}\theta \sqrt{2}\sqrt{M(u_3u_1)}+\delta ,\sqrt{\frac{u_2u_1}{u_3u_1}})^2$$

As $`u_1`$ = $`\frac{1}{r_1}`$ and $`u_2`$ = $`\frac{1}{r_2}`$ are the small values for the planets ( $`r_1`$ and $`r_2`$ are the perihelion and aphelion points) and $`u_1`$ \+ $`u_2`$ \+ $`u_3`$ = $`\frac{1}{2M}`$ , we have 2M $`u_3`$ =1 and

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

$$\mathrm{u}(\theta )u_1=(u_2u_1)\mathrm{sin}(\frac{1}{2}\theta +\delta )^2$$

that is the equation of orbital motion (see above) with e = $`\frac{u_2u_1}{u_2+u_1}`$ . $`u_1`$ $`>`$ 0 corresponds to the elliptical motion, $`u_1`$ $`<`$ 0 corresponds to hyperbolical motion. The period of the orbital motion in the general case:

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

$$2\mathrm{EllipticK}(\sqrt{\frac{u_2u_1}{u_3+u_1}})$$

and can be found approximately for small $`u_1`$ and $`u_2`$ as result of expansion:

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

$$\pi +\frac{1}{2}\pi M(u_2u_1)$$

From the expression for u( $`\theta `$ ) and obtained expression for the period of $`\mathrm{𝑠𝑛}^2`$ we have the change of angular coordinate over period:

%/(1/2\*sqrt( 2\*M\*(u-u) )):
simplify(%);

$$\frac{1}{2}\frac{\pi (2+Mu_2Mu_1)\sqrt{2}}{\sqrt{M(u_3+u_1)}}$$

But $`\sqrt{2M(u_3u_1)}`$ = $`\sqrt{12M(u_1u_2u_1)}`$ $``$ 1 + $`M(2u_1+u_2)`$ . And the result is

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

$$\pi (2+3Mu_1+3Mu_2)$$

The deviation of the period from 2 $`\pi `$ causes the shift of the perihelion over one rotation of the planet around massive star.

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

$$6\frac{G\pi M}{a(1+e)(1+e)c^2}$$

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(%);

$$24\frac{G\pi ^3a^2}{T^2(1+e)(1+e)c^2}$$

### 2.6 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.

## 3 Relativistic stars and black holes

### 3.1 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 .

### 3.2 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 $`\mathrm{𝑐𝑚}^{(2)}`$ ), the right-hand side is to have same dimension. Let’s the gravitational constant is G = $`\mathrm{6.673\hspace{0.17em}10}^{(8)}`$ $`\frac{\mathrm{𝑐𝑚}^3}{gs^2}`$ = 1 and the light velocity is c = $`\mathrm{2.998\hspace{0.17em}10}^{10}`$ $`\frac{\mathrm{𝑐𝑚}}{s}`$ = 1,

then

G/ $`c^2`$ = $`\mathrm{.7425\hspace{0.17em}10}^{(28)}`$  cm/g = 1

$`c^5`$ / G = $`\mathrm{3.63\hspace{0.17em}10}^{59}`$ erg/s = 1 (power unit)

G/c = $`\mathrm{2.23\hspace{0.17em}10}^{(18)}`$ Hz\* $`\mathrm{𝑐𝑚}^2`$ /g = 1 (characteristic of absorption)

$`c^2`$ / $`\sqrt{\text{G}}`$ = $`\mathrm{3.48\hspace{0.17em}10}^{24}`$ CGSE units (field strength)

h/2 $`\pi `$ = $`\mathrm{1.054\hspace{0.17em}10}^{(27)}`$ g\* $`\mathrm{𝑐𝑚}^2`$ /s = $`\mathrm{2.612\hspace{0.17em}10}^{(66)}`$  $`\mathrm{𝑐𝑚}^2`$

elementary charge e = $`\mathrm{1.381\hspace{0.17em}10}^{(34)}`$  cm

1 ps = $`\mathrm{3.0856\hspace{0.17em}10}^{18}`$ cm

1 eV = $`\mathrm{1.324\hspace{0.17em}10}^{(61)}`$  cm

There are the following extremal values of length, time, mass and density, which are useful in the context of the consederation of GR validity:

$`\sqrt{\frac{\text{G}\text{h}}{2\pi c^3}}`$ = $`\mathrm{1.616\hspace{0.17em}10}^{(33)}`$ cm (Planck length)

$`\sqrt{\frac{\text{G}\text{h}}{2\pi c^5}}`$ = $`\mathrm{5.391\hspace{0.17em}10}^{(44)}`$ s (Planck time)

$`\sqrt{\frac{hc}{2\pi \text{G}}}`$ = $`\mathrm{2.177\hspace{0.17em}10}^{(5)}`$ g (Planck mass)

$`\frac{2\pi c^5}{h\text{G}^2}`$ = $`\mathrm{5.157\hspace{0.17em}10}^{93}`$ g/ $`\mathrm{𝑐𝑚}^3`$ (Planck density)

### 3.3 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
g\_compts :=\
array(symmetric,sparse,1..4,1..4):# metric components
g\_compts := -exp(2\*Phi(r)):# component\
of interval attached to d(t)^2
g\_compts := exp(2\*Lambda(r)):# component\
of interval attached to d(r)^2
g\_compts := r^2:# component of interval\
attached to d(theta)^2
g\_compts := 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

$`g:=`$ $`\mathrm{table}([\mathrm{𝑐𝑜𝑚𝑝𝑡𝑠}=\left[\begin{array}{cccc}e^{(2\mathrm{\Phi }(r))}& 0& 0& 0\\ 0& e^{(2\mathrm{\Lambda }(r))}& 0& 0\\ 0& 0& r^2& 0\\ 0& 0& 0& r^2\mathrm{sin}(\theta )^2\end{array}\right],`$  $`\mathrm{𝑖𝑛𝑑𝑒𝑥}\mathrm{\_}\mathrm{𝑐ℎ𝑎𝑟}=[1,1]])`$

$`\mathrm{𝑔𝑖𝑛𝑣}:=`$ $`\mathrm{table}([\mathrm{𝑐𝑜𝑚𝑝𝑡𝑠}=\left[\begin{array}{cccc}{\displaystyle \frac{1}{e^{(2\mathrm{\Phi }(r))}}}& 0& 0& 0\\ 0& {\displaystyle \frac{1}{e^{(2\mathrm{\Lambda }(r))}}}& 0& 0\\ 0& 0& {\displaystyle \frac{1}{r^2}}& 0\\ 0& 0& 0& {\displaystyle \frac{1}{r^2\mathrm{sin}(\theta )^2}}\end{array}\right],`$  $`\mathrm{𝑖𝑛𝑑𝑒𝑥}\mathrm{\_}\mathrm{𝑐ℎ𝑎𝑟}=[1,\mathrm{\hspace{0.17em}1}]])`$

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

$$\mathrm{𝑇ℎ𝑒}\mathrm{𝐸𝑖𝑛𝑠𝑡𝑒𝑖𝑛}\mathrm{𝑇𝑒𝑛𝑠𝑜𝑟}$$

$$\mathrm{𝑛𝑜𝑛}\mathrm{𝑧𝑒𝑟𝑜}\mathrm{𝑐𝑜𝑚𝑝𝑜𝑛𝑒𝑛𝑡𝑠}:$$

$$\mathit{G11}=\frac{e^{(2\mathrm{\Phi }(r))}(2(\frac{}{r}\mathrm{\Lambda }(r))r+e^{(2\mathrm{\Lambda }(r))}1)}{r^2e^{(2\mathrm{\Lambda }(r))}}$$

$$\mathit{G22}=\frac{2(\frac{}{r}\mathrm{\Phi }(r))re^{(2\mathrm{\Lambda }(r))}+1}{r^2}$$

$`\mathit{G33}=\frac{r}{e^{(2\mathrm{\Lambda }(r))}}`$ $`\left[\frac{}{r}\mathrm{\Phi }(r)\frac{}{r}\mathrm{\Lambda }(r)+r\frac{^2}{r^2}\mathrm{\Phi }(r)+r(\frac{}{r}\mathrm{\Phi }(r))^2r\frac{}{r}\mathrm{\Phi }(r)\frac{}{r}\mathrm{\Lambda }(r)\right]`$

$`\mathit{G44}=\frac{\mathrm{sin}(\theta )^2r}{e^{(2\mathrm{\Lambda }(r))}}`$ $`\left[\frac{}{r}\mathrm{\Phi }(r)\frac{}{r}\mathrm{\Lambda }(r)+r\frac{^2}{r^2}\mathrm{\Phi }(r)+r(\frac{}{r}\mathrm{\Phi }(r))^2r\frac{}{r}\mathrm{\Phi }(r)\frac{}{r}\mathrm{\Lambda }(r)\right]`$

$$\mathrm{𝑐ℎ𝑎𝑟𝑎𝑐𝑡𝑒𝑟}:[\mathit{1},\mathit{1}]$$

In the beginning we will consider the star in the form of drop of liquid. In this case the energy-momentum tensor is $`T_{\mu ,\nu }`$ = ( $`p+\rho `$ ) $`u_\mu `$ $`u_\nu `$ \+ p $`g_{\mu ,\nu }`$ (all components of u except for $`u_0`$ are equal to zero, and -1 = $`g^{(0,\mathrm{\hspace{0.17em}0})}`$ $`u_0`$ $`u_0`$ ; the signature (-2) results in $`u_\alpha u^\alpha `$ =1 and $`T_{\mu ,\nu }`$ = ( $`p+\rho `$ ) $`u_\mu `$ $`u_\nu `$ \- p $`g_{\mu ,\nu }`$ ):

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

$`T:=`$ $`\mathrm{table}([\mathrm{𝑐𝑜𝑚𝑝𝑡𝑠}=\left[\begin{array}{cccc}e^{(2\mathrm{\Phi }(r))}\rho (r)& 0& 0& 0\\ 0& e^{(2\mathrm{\Lambda }(r))}\mathrm{p}(r)& 0& 0\\ 0& 0& \mathrm{p}(r)r^2& 0\\ 0& 0& 0& \mathrm{p}(r)r^2\mathrm{sin}(\theta )^2\end{array}\right],`$  $`\mathrm{𝑖𝑛𝑑𝑒𝑥}\mathrm{\_}\mathrm{𝑐ℎ𝑎𝑟}=[1,1]])`$

To write the Einstein equations (sign corresponds to )

$`G_{\mu ,\nu }`$= -8 $`\pi `$ $`T_{\mu ,\nu }`$

let’s extract the tensor components:

Energy\_momentum := get\_compts(T);
Einstein := get\_compts(Estn);

$`\mathrm{𝐸𝑛𝑒𝑟𝑔𝑦}\mathrm{\_}\mathrm{𝑚𝑜𝑚𝑒𝑛𝑡𝑢𝑚}:=`$ $`\left[\begin{array}{cccc}e^{(2\mathrm{\Phi }(r))}\rho (r)& 0& 0& 0\\ 0& e^{(2\mathrm{\Lambda }(r))}\mathrm{p}(r)& 0& 0\\ 0& 0& \mathrm{p}(r)r^2& 0\\ 0& 0& 0& \mathrm{p}(r)r^2\mathrm{sin}(\theta )^2\end{array}\right]`$

$`\mathrm{𝐸𝑖𝑛𝑠𝑡𝑒𝑖𝑛}:=`$ $`[{\displaystyle \frac{e^{(2\mathrm{\Phi }(r))}(2(\frac{}{r}\mathrm{\Lambda }(r))r+\%11)}{r^2\%1}},\mathrm{\hspace{0.17em}0},\mathrm{\hspace{0.17em}0},\mathrm{\hspace{0.17em}0}]`$  $`[0,{\displaystyle \frac{2(\frac{}{r}\mathrm{\Phi }(r))r\%1+1}{r^2}},\mathrm{\hspace{0.17em}0},\mathrm{\hspace{0.17em}0}]`$  $`[0,\mathrm{\hspace{0.17em}0},`$  $`{\displaystyle \frac{r((\frac{}{r}\mathrm{\Phi }(r))(\frac{}{r}\mathrm{\Lambda }(r))+r(\frac{^2}{r^2}\mathrm{\Phi }(r))+r(\frac{}{r}\mathrm{\Phi }(r))^2r(\frac{}{r}\mathrm{\Phi }(r))(\frac{}{r}\mathrm{\Lambda }(r)))}{\%1}},\mathrm{\hspace{0.17em}0}]`$  $`[0,\mathrm{\hspace{0.17em}0},\mathrm{\hspace{0.17em}0},`$  $`{\displaystyle \frac{\mathrm{sin}(\theta )^2r((\frac{}{r}\mathrm{\Phi }(r))(\frac{}{r}\mathrm{\Lambda }(r))+r(\frac{^2}{r^2}\mathrm{\Phi }(r))+r(\frac{}{r}\mathrm{\Phi }(r))^2r(\frac{}{r}\mathrm{\Phi }(r))(\frac{}{r}\mathrm{\Lambda }(r)))}{\%1}}`$  $`]`$  $`\%1:=e^{(2\mathrm{\Lambda }(r))}`$

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

8\*Pi\*Energy\_momentum + Einstein:
expand(%/exp(Phi(r))^2):
eq1 := simplify(%) = 0;#first Einstein equation

$$\mathit{eq1}:=\frac{8\pi \rho (r)r^2+2e^{(2\mathrm{\Lambda }(r))}(\frac{}{r}\mathrm{\Lambda }(r))r+1e^{(2\mathrm{\Lambda }(r))}}{r^2}=0$$

This equation can be rewritten as:

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

$$\mathit{eq1}:=8\pi \rho (r)r^2+(\frac{}{r}r(1e^{(2\mathrm{\Lambda }(r))}))=0$$

The formal substitution results in:

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

$$\mathit{eq1}\mathrm{\_}n:=8\pi \rho (r)r^2+(\frac{}{r}(2\mathrm{m}(r)))=0$$

dsolve(eq1\_n,m(r));

$$\mathrm{m}(r)=4\pi \rho (r)r^2𝑑r+\mathrm{\_}\mathit{C1}$$

So, m(r) is the mass inside sphere with radius r. To estimate \_C1 we have to express $`\mathrm{\Lambda }`$ from m:

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

$$1\frac{2\mathrm{m}(r)}{r}$$

Hence (see expression for $`g_{\mu ,\nu }`$ through $`\mathrm{\Lambda }`$ ) 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 +\
Einstein ) = 0;#second Einstein equation

$$\mathit{eq2}:=\frac{8\pi e^{(2\mathrm{\Lambda }(r))}\mathrm{p}(r)r^22(\frac{}{r}\mathrm{\Phi }(r))r+e^{(2\mathrm{\Lambda }(r))}1}{r^2}=0$$

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

$$\mathit{eq2}\mathrm{\_}\mathit{2}:=8\pi r^3\mathrm{p}(r)+2(\frac{}{r}\mathrm{\Phi }(r))r^24(\frac{}{r}\mathrm{\Phi }(r))r\mathrm{m}(r)2\mathrm{m}(r)=0$$

As result we have:

eq2\_3 := Diff(Phi(r),r) = solve(eq2\_2, diff(Phi(r),r) );

$$\mathit{eq2}\mathrm{\_}\mathit{3}:=\frac{}{r}\mathrm{\Phi }(r)=\frac{4\pi r^3\mathrm{p}(r)+\mathrm{m}(r)}{r(r+2\mathrm{m}(r))}$$

We can see, that the gradient of the gravitational potential $`\mathrm{\Phi }`$ is greater than in the Newtonian case $`\frac{}{r}\mathrm{\Phi }`$ = $`\frac{m}{r^2}`$ , 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(, compts):# 4-velocity
Cf2 := Christoffel2 ( ginv, Cf1 ):
(rho(r)+p(r))\*get\_compts(\
cov\_diff( u, coord, Cf2 ))/(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) );

$$(\rho (r)+\mathrm{p}(r))(\frac{}{r}\mathrm{\Phi }(r))=(\frac{}{r}\mathrm{p}(r))$$

$$\mathit{eq3}:=\frac{}{r}\mathrm{\Phi }(r)=\frac{\frac{}{r}\mathrm{p}(r)}{\rho (r)+\mathrm{p}(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)));

$$\frac{}{r}\mathrm{p}(r)=\frac{(4\pi r^3\mathrm{p}(r)+\mathrm{m}(r))(\rho (r)+\mathrm{p}(r))}{r(r+2\mathrm{m}(r))}$$

One can see that the pressure gradient is greater than in the classical limit ( $`\frac{}{r}p`$ = \- $`\frac{\rho m}{r^2}`$ ) 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));

$$\frac{}{r}\mathrm{\Phi }(r)=\frac{M}{r(r+2M)}$$

$$\mathit{eq4}:=\mathrm{\Phi }(r)=\frac{1}{2}\mathrm{ln}(r)+\frac{1}{2}\mathrm{ln}(r2M)+\mathrm{\_}\mathit{C1}$$

The boundary condition

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

$$0=\mathrm{\_}\mathit{C1}$$

results in

subs(\_C1=0,eq4);

$$\mathrm{\Phi }(r)=\frac{1}{2}\mathrm{ln}(r)+\frac{1}{2}\mathrm{ln}(r2M)$$

So, Schwarzschild metric out of star is:

g\_matrix := get\_compts(g);

$$g\mathrm{\_}\mathrm{𝑚𝑎𝑡𝑟𝑖𝑥}:=\left[\begin{array}{cccc}e^{(2\mathrm{\Phi }(r))}& 0& 0& 0\\ 0& e^{(2\mathrm{\Lambda }(r))}& 0& 0\\ 0& 0& r^2& 0\\ 0& 0& 0& r^2\mathrm{sin}(\theta )^2\end{array}\right]$$

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

$`\mathrm{𝑠𝑐ℎ}:=`$ $`\mathrm{table}([\mathrm{𝑐𝑜𝑚𝑝𝑡𝑠}=\left[\begin{array}{cccc}1+{\displaystyle \frac{2M}{r}}& 0& 0& 0\\ 0& {\displaystyle \frac{1}{1{\displaystyle \frac{2M}{r}}}}& 0& 0\\ 0& 0& r^2& 0\\ 0& 0& 0& r^2\mathrm{sin}(\theta )^2\end{array}\right],`$  $`\mathrm{𝑖𝑛𝑑𝑒𝑥}\mathrm{\_}\mathrm{𝑐ℎ𝑎𝑟}=[1,1]])`$

Now we consider the star, which is composed of an incompressible substance $`\rho `$ = $`\rho `$ 0 = const (later we will use also the following approximation: p=( $`\gamma `$ -1) $`\rho `$ , where $`\gamma `$ =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

$$\mathrm{m}(r)=\frac{4}{3}\pi \rho 0r^3$$

$$M=\frac{4}{3}\pi \rho 0R^3$$

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));

$$\frac{}{r}\mathrm{p}(r)=\frac{(4\pi r^3\mathrm{p}(r)+{\displaystyle \frac{4}{3}}\pi \rho 0r^3)(\rho 0+\mathrm{p}(r))}{r(r{\displaystyle \frac{8}{3}}\pi \rho 0r^3)}$$

$`\mathit{eq5}:=\mathrm{p}(r)=\rho 0({\displaystyle \frac{6\%1+2\sqrt{3\%1+8\%1\pi \rho 0r^2}}{3+8\pi \rho 0r^29\%1}}1),`$ $`\mathrm{p}(r)=\rho 0({\displaystyle \frac{6\%12\sqrt{3\%1+8\%1\pi \rho 0r^2}}{3+8\pi \rho 0r^29\%1}}1)`$  $`\%1:=e^{(16\mathrm{\_}\mathit{C1}\pi \rho 0)}`$

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

$$3+8\pi \rho 0R^2$$

$$3+8\pi \rho 0R^2$$

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

$`\mathit{sol1}:={\displaystyle \frac{\rho 0}{4}}`$ $`{\displaystyle \frac{3+12\pi \rho 0R^2+\sqrt{924\pi \rho 0R^224\pi \rho 0r^2+64\pi ^2\rho 0^2r^2R^2}4\pi \rho 0r^2}{3\pi \rho 0r^2+9\pi \rho 0R^2}}`$

$`\mathit{sol2}:={\displaystyle \frac{\rho 0}{4}}`$ $`{\displaystyle \frac{312\pi \rho 0R^2+\sqrt{924\pi \rho 0R^224\pi \rho 0r^2+64\pi ^2\rho 0^2r^2R^2}+4\pi \rho 0r^2}{3\pi \rho 0r^2+9\pi \rho 0R^2}}`$

In the center of star we have

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

$$\frac{1}{12}\frac{\rho 0(3+12\pi \rho 0R^2+\sqrt{924\pi \rho 0R^2})}{1+3\pi \rho 0R^2}$$

$$\frac{1}{12}\frac{\rho 0(3+12\pi \rho 0R^2\sqrt{924\pi \rho 0R^2})}{1+3\pi \rho 0R^2}$$

The pressure is infinity when the radius is equal to

R\_crit1 =\
simplify(solve(expand(denom(%)/12)=0,R),\
radical,symbolic);
R\_crit2 =\
simplify(solve(expand(denom(%%)/12)=0,R),\
radical,symbolic);

$$R\mathrm{\_}\mathit{crit1}=\frac{1}{3}\frac{\sqrt{3}\sqrt{\pi \rho 0}}{\pi \rho 0}$$

$$R\mathrm{\_}\mathit{crit2}=\frac{1}{3}\frac{\sqrt{3}\sqrt{\pi \rho 0}}{\pi \rho 0}$$

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 ( $`\theta `$ = $`\pi `$ /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 $`\mathrm{𝑑𝑠}^2`$ = $`\frac{\mathrm{𝑑𝑟}^2}{1\frac{2m}{r}}`$ \+ $`r^2`$ $`\mathrm{𝑑𝑝ℎ𝑖}^2`$ and for Euclidean 3-space we have $`\mathrm{𝑑𝑠}^2`$ = $`\mathrm{𝑑𝑧}^2`$ \+ $`\mathrm{𝑑𝑟}^2`$ \+ $`r^2`$ $`\mathrm{𝑑𝑝ℎ𝑖}^2`$ . We will investigate the 2-surface z=z(r). As dz= $`\frac{\mathrm{𝑑𝑧}}{\mathrm{𝑑𝑟}}`$ dr, one can obtain Euclidian line element: $`\mathrm{𝑑𝑠}^2`$ = \[1+ $`(\frac{\mathrm{𝑑𝑧}}{\mathrm{𝑑𝑟}})^2`$ \] $`\mathrm{𝑑𝑟}^2`$ \+ $`r^2`$ $`\mathrm{𝑑𝑝ℎ𝑖}^2`$

subs(theta=Pi/2,get\_compts(sch));#Schwarzschild metric\
dr^2\*% + dphi^2\*% = \
(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));
dsolve(%,z(r));# embedding

$$\left[\begin{array}{cccc}1+\frac{2M}{r}& 0& 0& 0\\ 0& \frac{1}{1{\displaystyle \frac{2M}{r}}}& 0& 0\\ 0& 0& r^2& 0\\ 0& 0& 0& r^2\mathrm{sin}(\frac{1}{2}\pi )^2\end{array}\right]$$

$$\frac{\mathrm{𝑑𝑟}^2}{1{\displaystyle \frac{2M}{r}}}+r^2\mathrm{𝑑𝑝ℎ𝑖}^2=(1+(\frac{}{r}\mathrm{z}(r))^2)\mathrm{𝑑𝑟}^2+r^2\mathrm{𝑑𝑝ℎ𝑖}^2$$

$$\frac{}{r}\mathrm{z}(r)=\frac{\sqrt{2}\sqrt{(r2M)M}}{r2M}$$

$$\mathrm{z}(r)=2\sqrt{2rM4M^2}+\mathrm{\_}\mathit{C1}$$

Now let’s take the penultimate equation and to express M through $`\rho _0`$ = 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

$$\mathrm{𝑟ℎ𝑜}\mathrm{\_}\mathrm{𝑠𝑜𝑙}:=\frac{3}{4}\frac{M}{\pi R^3}$$

$$\mathit{fun1}:=2\frac{1}{\sqrt{3{\displaystyle \frac{1}{r^2\pi \rho 0}}4}}𝑑r$$

$$\mathit{fun2}:=2\frac{1}{\sqrt{3{\displaystyle \frac{r}{\pi \rho 0}}4}}𝑑r$$

fun3 := value(subs(rho0=rho\_sol,fun1));# inside
fun4 := value(subs(rho0=rho\_sol,fun2));# outside

$$\mathit{fun3}:=\frac{R^3+r^2M}{\sqrt{{\displaystyle \frac{R^3+r^2M}{r^2M}}}rM}$$

$$\mathit{fun4}:=\frac{\sqrt{4{\displaystyle \frac{rR^3}{M}}4}M}{R^3}$$

The resulting embedding for equatorial and vertical sections is presented below (Newtonian case corresponds to horizontal surface, that is the asymptote for r–$`>`$ $`\mathrm{}`$ ). 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=,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= $`\frac{\mathrm{𝑑𝑟}}{\sqrt{1\frac{2M}{r}}}`$  (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=2M? 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: $`g^{(\alpha ,\beta )}`$ $`p_\alpha `$ $`p_\beta `$ \+ $`\mu ^2`$ = 0 (p is 4-momentum, $`\mu `$ is the rest mass). Then (if p =\[ $`E,\frac{}{\lambda }r,\mathrm{\hspace{0.17em}0},L`$ \], where L corresponds to angle momentum, $`\lambda `$ 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

$$\frac{E^2}{1{\displaystyle \frac{2M}{\mathrm{r}(\tau )}}}+\frac{(\frac{}{\tau }\mathrm{r}(\tau ))^2}{1{\displaystyle \frac{2M}{\mathrm{r}(\tau )}}}+\frac{L^2}{\mathrm{r}(\tau )^2}+1=0$$

$$\mathit{eq6}:=(\frac{}{\tau }\mathrm{r}(\tau ))^2=E^2+\frac{2M}{\mathrm{r}(\tau )}\frac{L^2}{\mathrm{r}(\tau )^2}+\frac{2L^2M}{\mathrm{r}(\tau )^3}1$$

$$V:=\sqrt{\frac{(\mathrm{r}(\tau )^2+L^2)(\mathrm{r}(\tau )2M)}{\mathrm{r}(\tau )^3}}$$

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 a 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 $`V_{\mathrm{𝑚𝑎𝑥}}`$ with the following potential decrease as result of radial coordinate decrease. If $`V_{\mathrm{𝑚𝑎𝑥}}`$ $`>`$E $`>`$ 1 the motion is infinite (it is analogue of the hyperbolical motion in Newtonian case). $`V_{\mathrm{𝑚𝑎𝑥}}`$ $`>`$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);

$$\frac{1}{2}\frac{(L+\sqrt{L^212M^2})L}{M},\frac{1}{2}\frac{(L\sqrt{L^212M^2})L}{M}$$

As consequence, the circle motion is possible if L $``$ 2M $`\sqrt{3}`$ , when r $``$ 3(2M). 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);

$$3M$$

Lack of extremes (red curve) or transition to r $`<`$ 3M results in the falling motion (the similar case is E $`>`$ $`V_{\mathrm{𝑚𝑎𝑥}}`$ ). 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

$$\tau =\frac{1}{\sqrt{E^2+{\displaystyle \frac{2M}{r}}1}}𝑑r$$

$$\tau =\frac{1}{\sqrt{2{\displaystyle \frac{M}{r}}{\displaystyle \frac{2M}{R}}}}𝑑r$$

$`{\displaystyle \frac{1}{4}}\sqrt{{\displaystyle \frac{2M+R}{R}}}R(2\sqrt{2}\sqrt{(2M+R)M}+R\mathrm{ln}(2)`$ $`R\mathrm{ln}(R+4M+2\sqrt{2}\sqrt{(2M+R)M}))\sqrt{2}`$  $`/\sqrt{(2M+R)M}`$

For the remote exterior observer we have $`\frac{}{\tau }r`$ = $`\frac{}{t}r`$ $`\frac{}{\tau }t`$ = $`\frac{}{t}r`$ $`\frac{E}{1\frac{2M}{r}}`$ =
E $`\frac{}{t}r_o`$  (here $`r_0`$ 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

$$r_o:=\frac{1}{1{\displaystyle \frac{2M}{r}}}𝑑r=r+2M\mathrm{ln}(r2M)$$

$$\mathrm{signum}(M)\mathrm{}$$

These equations will be utilized below.

### 3.4 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 2M ( $`\frac{2GM}{c^2}`$ in dimensional case, that is the so-called ”gravitational radius$`R_g`$ ). The finite proper time $`\tau `$ of fall corresponds to the infinite ”remote” time $`r_o`$ , when r–$`>`$2M. This means, that for the remote observer the fall does not finish never. It results from the relativistic 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: $`\frac{\omega _1}{\omega _2}`$ = $`\frac{\mathrm{\Delta }\tau _2}{\mathrm{\Delta }\tau _1}`$ = $`\sqrt{\frac{g_{0,\mathrm{\hspace{0.17em}0}}(r_2)}{g_{0,\mathrm{\hspace{0.17em}0}}(r_1)}}`$ = $`\frac{\sqrt{1\frac{2M}{r_2}}}{\sqrt{1\frac{2M}{r_1}}}`$ (here $`\mathrm{\Delta }`$ $`\tau `$ are the proper time intervals between light flash in different radial points). Simultaneously, the escape velocity is equal to $`\sqrt{\frac{2GM}{R}}`$ that is the velocity of light for sphere with R= $`R_g`$ . So, we cannot receive any information from interior of black hole.

How can appear the object with the radius, which is less than $`R_g`$ ? 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 =\
%\*d(t)^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"

$$\mathrm{d}(s)^2=(1+\frac{2M}{\mathrm{r}(t)})\mathrm{d}(t)^2+\frac{\mathrm{d}(r)^2}{1{\displaystyle \frac{2M}{\mathrm{r}(t)}}}$$

$$\mathrm{d}(\tau )^2=\left(1+\frac{2M}{\mathrm{r}(t)}+\frac{(\frac{}{t}\mathrm{r}(t))^2}{1{\displaystyle \frac{2M}{\mathrm{r}(t)}}}\right)\mathrm{d}(t)^2$$

$$1=\frac{\left(1+{\displaystyle \frac{2M}{\mathrm{r}(t)}}+{\displaystyle \frac{(\frac{}{t}\mathrm{r}(t))^2}{1{\displaystyle \frac{2M}{\mathrm{r}(t)}}}}\right)\mathrm{d}(t)^2}{\mathrm{d}(\tau )^2}$$

$$1=\frac{\left(1+{\displaystyle \frac{2M}{\mathrm{r}(t)}}+{\displaystyle \frac{(\frac{}{t}\mathrm{r}(t))^2}{1{\displaystyle \frac{2M}{\mathrm{r}(t)}}}}\right)E^2}{(1{\displaystyle \frac{2M}{\mathrm{r}(t)}})^2}$$

$$\mathrm{𝑝𝑜𝑡}\mathrm{\_}\mathit{1}:=\frac{(\mathrm{r}(t)2M)^2(\mathrm{r}(t)+E^2\mathrm{r}(t)+2M)}{E^2\mathrm{r}(t)^3}$$

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 $`\frac{}{t}r`$ (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);

$$2M,\mathrm{\hspace{0.17em}2}M,2\frac{M}{1+E^2}$$

that are the apoastr (see above) and the gravitational radius. Hence (from the value for apoastr and pot\_1)
$`(\frac{}{t}r)^2`$ = $`\frac{(1\frac{2M}{r})^2(\frac{2M}{r}1+E^2)}{E^2}`$ = $`\frac{(1\frac{2M}{r})^2(\frac{2M}{r}\frac{2M}{R})}{1\frac{2M}{R}}`$ = $`(1\frac{R_g}{r})^2(1\frac{1\frac{R_g}{r}}{1\frac{R_g}{R}})`$.

The slowing-down of the collapse for remote observer in vicinity of $`R_g`$ results from the time slowing-down (see above). But the collapsing observer has the proper velocity $`\frac{}{\tau }r`$ :

pot\_2 := simplify(pot\_1\*(E/(1-2\*M/r(t)))^2);#
d/d(t)=(d/d(tau))\*(1-2\*M/r)/E

$$\mathrm{𝑝𝑜𝑡}\mathrm{\_}\mathit{2}:=\frac{\mathrm{r}(t)+E^2\mathrm{r}(t)+2M}{\mathrm{r}(t)}$$

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 $`lim_{rR_g}\frac{}{t}r`$ =1 (that is the velocity of light), when E = 1 (R= $`\mathrm{}`$ ). 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);

$$\frac{{\displaystyle \frac{1}{\sqrt{{\displaystyle \frac{R}{r}}1}}𝑑r}}{\sqrt{1E^2}}$$

$$\frac{1}{2}\frac{2\sqrt{r(Rr)}+R\mathrm{arctan}({\displaystyle \frac{1}{2}}{\displaystyle \frac{R2r}{\sqrt{r(Rr)}}})}{\sqrt{1E^2}}$$

It is obvious, that the obtained integral is $`\frac{1}{\sqrt{1E^2}}`$ $`_r^R\frac{1}{\sqrt{\frac{R}{r}1}}𝑑r`$ = $`\frac{\pi R}{\sqrt{1E^2}}`$ = $`\frac{\pi M}{(1E^2)^{(\frac{3}{2})}}`$ .

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 $`\frac{GM^2}{R}`$ 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 ~ $`\lambda ^3`$ , where $`\lambda `$ = $`\frac{h}{2\pi p_e}`$ is the Compton wavelength ( $`p_e`$ 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

$$E_e:=\frac{p_{e}^{}{}_{}{}^{2}}{m_e}$$

$$E_k:=\frac{1}{4}\frac{n_{e}^{}{}_{}{}^{(5/3)}R^3h^2}{\pi ^2m_e}$$

$$E_k:=\frac{1}{4}\frac{({\displaystyle \frac{M}{R^3m_p}})^{(5/3)}R^3h^2}{\pi ^2m_e}$$

$$E:=\frac{GM^2}{R}+\frac{{\displaystyle \frac{1}{4}}({\displaystyle \frac{M}{R^3m_p}})^{(5/3)}R^3h^2}{\pi ^2m_e}$$

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

$`{\displaystyle \frac{1}{2}}{\displaystyle \frac{\%1h^2}{Mm_{p}^{}{}_{}{}^{2}\pi ^2m_eG}},{\displaystyle \frac{1}{2}}{\displaystyle \frac{\left({\displaystyle \frac{1}{2}}{\displaystyle \frac{\%1}{Mm_p}}+{\displaystyle \frac{{\displaystyle \frac{1}{2}}I\sqrt{3}\%1}{Mm_p}}\right)h^2}{\pi ^2m_pm_eG}},`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \frac{\left({\displaystyle \frac{1}{2}}{\displaystyle \frac{\%1}{Mm_p}}+{\displaystyle \frac{{\displaystyle \frac{1}{2}}I\sqrt{3}\%1}{Mm_p}}\right)h^2}{\pi ^2m_pm_eG}}`$  $`\%1:=(M^2m_p)^{(1/3)}`$

So, we have the equilibrium state with $`R_{\mathrm{𝑚𝑖𝑛}}`$ ~ $`\frac{h^2}{GM^{(\frac{1}{3})}m_em_{p}^{}{}_{}{}^{(\frac{5}{3})}}`$ 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

$$n_e=\frac{M^2G^3m_{e}^{}{}_{}{}^{3}m_{p}^{}{}_{}{}^{4}}{h^6}$$

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 solid or liquid defined by atomic packing is ~20 g/ $`\mathrm{𝑐𝑚}^3`$ , for white dwarf this value is ~ $`10^7`$ g/ $`\mathrm{𝑐𝑚}^3`$ !). But, in compliance with principle of uncertainty, the last causes the electron momentum increase ( $`p_en_{e}^{}{}_{}{}^{(\frac{1}{3})}`$ ~ h). Our nonrelativistic approximation implies $`p_e`$ $`<`$$`<`$ $`m_ec`$ or

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

$$(\frac{M^2G^3m_{e}^{}{}_{}{}^{3}m_{p}^{}{}_{}{}^{4}}{h^6})^{(1/3)}hcm_e$$

$$\frac{\sqrt{Gch}hc}{G^2m_{p}^{}{}_{}{}^{2}},\frac{\sqrt{Gch}hc}{G^2m_{p}^{}{}_{}{}^{2}}$$

The last result gives the nonrelativistic criterion M $`<`$$`<`$ $`\frac{(hc)^{(\frac{3}{2})}}{m_{p}^{}{}_{}{}^{2}G^{(\frac{3}{2})}}`$ . 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 );

$$E_e:=hcn_{e}^{}{}_{}{}^{(1/3)}$$

$$E_k:=n_{e}^{}{}_{}{}^{(4/3)}R^3hc$$

$$E_k:=(\frac{M}{R^3m_p})^{(4/3)}R^3hc$$

$$E:=\frac{GM^2}{R}+(\frac{M}{R^3m_p})^{(4/3)}R^3hc$$

$$(GMm_p+(\frac{M}{R^3m_p})^{(1/3)}hcR)M$$

$$0,\frac{\sqrt{Gch}hc}{G^2m_{p}^{}{}_{}{}^{2}},\frac{\sqrt{Gch}hc}{G^2m_{p}^{}{}_{}{}^{2}}$$

We obtained the estimation for so-called critical mass $`M_c`$ ~ $`\frac{(hc)^{(\frac{3}{2})}}{m_{p}^{}{}_{}{}^{2}G^{(\frac{3}{2})}}`$ =1.4 $`M_{\mathrm{𝑠𝑢𝑛}}`$ (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_{\mathrm{𝑠𝑢𝑛}}`$ $`<`$ M $`<`$ 3 $`M_{\mathrm{𝑠𝑢𝑛}}`$. Collapse for larger masses has to result in the black hole formation.

### 3.5 Schwarzschild black hole

Now we return to Schwarzschild metric.

get\_compts(sch);

$$\left[\begin{array}{cccc}1+\frac{2M}{r}& 0& 0& 0\\ 0& \frac{1}{1{\displaystyle \frac{2M}{r}}}& 0& 0\\ 0& 0& r^2& 0\\ 0& 0& 0& r^2\mathrm{sin}(\theta )^2\end{array}\right]$$

One can see two singularities: r=2M and r=0. What is a sense of first singularity? When we cross the horizon, $`\mathrm{𝑠ℎ}_{1,\mathrm{\hspace{0.17em}1}}`$ and $`\mathrm{𝑠ℎ}_{2,\mathrm{\hspace{0.17em}2}}`$ (i.e. $`g_{0,\mathrm{\hspace{0.17em}0}}`$ and $`g_{1,\mathrm{\hspace{0.17em}1}}`$ ) change signs. The space and time exchange the roles! The fall gets inevitable as the time flowing. As consequence, when particle or signal cross the gravitational radius, they cannot escape the falling on r=0. This fact can be illustrated by infinite value of acceleration on r=2M, which is \- $`\frac{(\mathrm{\Gamma }^\alpha )_{0,\mathrm{\hspace{0.17em}0}}}{g_{0,\mathrm{\hspace{0.17em}0}}}`$  ( $`\mathrm{\Gamma }`$ are the Christoffel symbols):

D1sch := d1metric( sch, coord ):
Cf1 := Christoffel1 ( D1sch ):
displayGR(Christoffel1,%);

$$\mathrm{𝑇ℎ𝑒}\mathrm{𝐶ℎ𝑟𝑖𝑠𝑡𝑜𝑓𝑓𝑒𝑙}\mathrm{𝑆𝑦𝑚𝑏𝑜𝑙𝑠}\mathrm{𝑜𝑓}\mathrm{𝑡ℎ𝑒}\mathrm{𝐹𝑖𝑟𝑠𝑡}\mathrm{𝐾𝑖𝑛𝑑}$$

$$\mathrm{𝑛𝑜𝑛}\mathrm{𝑧𝑒𝑟𝑜}\mathrm{𝑐𝑜𝑚𝑝𝑜𝑛𝑒𝑛𝑡𝑠}:$$

$$[\mathit{11},\mathit{2}]=\frac{M}{r^2}$$

$$[\mathit{12},\mathit{1}]=\frac{M}{r^2}$$

$$[\mathit{22},\mathit{2}]=\frac{M}{(r2M)^2}$$

$$[\mathit{23},\mathit{3}]=r$$

$$[\mathit{24},\mathit{4}]=r\mathrm{sin}(\theta )^2$$

$$[\mathit{33},\mathit{2}]=r$$

$$[\mathit{34},\mathit{4}]=r^2\mathrm{sin}(\theta )\mathrm{cos}(\theta )$$

$$[\mathit{44},\mathit{2}]=r\mathrm{sin}(\theta )^2$$

$$[\mathit{44},\mathit{3}]=r^2\mathrm{sin}(\theta )\mathrm{cos}(\theta )$$

-get\_compts(Cf1)/get\_compts(sch);#radial\
component of acceleration

$$\frac{M}{r^2(1+{\displaystyle \frac{2M}{r}})}$$

Such particles and signals will be expelled from the cause-effect chain of universe. Therefore an imaginary surface r= $`R_g`$ is named ”event horizon”.

The absence of true physical singularity for r= $`R_g`$ can be illustrated in the following way. The invariant of Riemann tensor $`R_{\alpha ,\beta ,\gamma ,\delta }R^{(\alpha ,\beta ,\gamma ,\delta )}`$ is

schinv := invert( sch, ’detg’ ):
D2sch := d2metric( D1sch, coord ):
Cf1 := Christoffel1 ( D1sch ):
RMN := Riemann( schinv, D2sch, Cf1 ):
raise(schinv,RMN,1):#raise of indexes\
in Riemann tensor
raise(schinv,%,2):
raise(schinv,%,3):
RMNinv := raise(schinv,%,4):
prod(RMN,RMNinv,,,,);

$$\mathrm{table}([\mathrm{𝑐𝑜𝑚𝑝𝑡𝑠}=48\frac{M^2}{r^6},\mathrm{𝑖𝑛𝑑𝑒𝑥}\mathrm{\_}\mathrm{𝑐ℎ𝑎𝑟}=[]])$$

and has no singularity on horizon. This is only coordinate singularity and its sense is the lack of rigid coordinates inside horizon. True physical singularity is r=0 and has the character of the so-called space-like singularity (the inavitable singularity for the observer crossing horizon).

Now try to embed the instant equator section of curved space into flat space (see above):

z(r) = int(sqrt(2\*r\*M-4\*M^2)/(-r+2\*M),r);
z(r) = int(-sqrt(2\*r\*M-4\*M^2)/(-r+2\*M),r);

$$\mathrm{z}(r)_1=2\sqrt{2rM4M^2}$$

$$\mathrm{z}(r)_2=2\sqrt{2rM4M^2}$$

plot3d({subs({M=1,r=sqrt(x^2+y^2)},rhs(%%)),\
subs({M=1,r=sqrt(x^2+y^2)},rhs(%))},\
x=-10..10,y=-10..10,axes=boxed,style=PATCHCONTOUR,\
grid=,title=‘Schwarzschild black hole‘);

This is the black hole, which looks as a neck of battle (Einstein-Rosen bridge): the inside way to horizon and the outside way to horizon are the ways between different but identically asymptotically flat universes (”wormhole” through 2-dimensional sphere with minimal radius 2M). But since the static geometry is not valid upon horizon (note, that here t–$`>`$t+dt is not time translation), this scheme is not stable.

Let’s try to exclude the above-mentioned coordinate singularity on the horizon by means of coordinate change. For example, we consider the null geodesics ( $`\mathrm{𝑑𝑠}^2`$ =0) in Schwarzschild space-time corresponding to radial motion of photons: $`\mathrm{𝑑𝑡}^2`$ = $`\frac{\mathrm{𝑑𝑟}^2}{(1\frac{R_g}{r})^2}`$ . Hence (see above considered expression for $`r_0`$ ):

t = combine( int(-1/(1-R\[g\]/r),r) );

$$t=rR_g\mathrm{ln}(rR_g)$$

Then if v is the constant, which defines the radial coordinate for fixed t, we have

t = -r - R\[g\]\*ln( abs(-r/R\[g\]+1) ) + v;#module\
allows to extend the expression for r\<R\[g\]

$$t=rR_g\mathrm{ln}(\left|\frac{r}{R_g}1\right|)+v$$

Differentiation of this equation with subsequent substitution of $`\mathrm{𝑑𝑡}^2`$ in expression for interval in Schwarzschild metric produce

defform(f=0,w1=1,w2=1,w3=1,v=1,R\[g\]=0,r=0,v=0);
d(t)^2 = expand(\
subs( d(R\[g\])=0,\
d( -r-R\[g\]\*ln( r/R\[g\]-1 )+v ) )^2);# differentials\
in new coordinates
subs( d(t)^2=rhs(%),sch\_compts\*d(t)^2 ) +\
sch\_compts\*d(r)^2 + sch\_compts\*d(theta)^2 +\
sch\_compts\*d(phi)^2:
collect( simplify( subs(R\[g\]=2\*M,%) ),\
{d(r)^2,d(v)^2});#new metric

$$\mathrm{d}(t)^2=\mathrm{d}(v)^2\frac{2\mathrm{d}(v)r\mathrm{d}(r)}{rR_g}+\frac{r^2\mathrm{d}(r)^2}{(rR_g)^2}$$

$$\frac{(r2M)\mathrm{d}(v)^2}{r}+2\mathrm{d}(r)\mathrm{d}(v)\frac{r^3\mathrm{d}(\theta )^2r^3\mathrm{d}(\varphi )^2+r^3\mathrm{d}(\varphi )^2\mathrm{cos}(\theta )^2}{r}$$

So, we have a new linear element (in the so-called Eddington - Finkelstein coordinates) $`\mathrm{𝑑𝑠}^2`$ = - ( $`1\frac{R_g}{r}`$ ) $`\mathrm{𝑑𝑉}^2`$ \+ 2dvdr + $`r^2`$ $`\mathrm{d}(\mathrm{\Omega })^2`$  ( $`\mathrm{d}(\mathrm{\Omega })`$ is the spherical part). The corresponding metric has the regular character in all region of r (except for r=0). It should be noted, that the regularization was made by transition to ”light coordinate” v therefore such coordinates can not be realized physically, but formally we continued analytically the coordinates to all r$`>`$0. We can see, that for future directed (i.e. $`\mathrm{𝑑𝑣}^2`$ $`>`$0) null ( $`\mathrm{𝑑𝑠}^2`$ =0) or time like ( $`\mathrm{𝑑𝑠}^2`$ $`<`$0) worldlines dr$`<`$0 for r$`<`$ $`R_g`$ that corresponds to above-mentioned conclusion about inevitable fall on singularity.

In the conclusion, we consider the problem of the deviation from spherical symmetry for static black hole. Such deviation can be described by the characteristic of quadrupole momentum q. Erez and Rosen found the corresponding static metric with axial symmetry:

coord := \[t, lambda, mu, phi\]:
er\_compts :=\
array(symmetric,sparse,1..4,1..4):# metric components
er\_compts :=\
-exp(2\*psi):# coefficient of d(t)^2 in interval
er\_compts :=\
M^2\*exp(2\*gamma-2\*psi)\*(lambda^2-mu^2)/\
(lambda^2-1):# coefficient of d(lambda)^2 in interval
er\_compts :=\
M^2\*exp(2\*gamma-2\*psi)\*(lambda^2-mu^2)/\
(1-mu^2):# coefficient of d(mu)^2 in interval
er\_compts :=\
M^2\*exp(-2\*psi)\*(lambda^2-1)\*(1-mu^2):#coefficient\
of d(phi)^2 in interval
er := create(\[-1,-1\], eval(er\_compts));# axially\
symmetric metric

$`\mathrm{𝑒𝑟}:=\mathrm{table}([\mathrm{𝑐𝑜𝑚𝑝𝑡𝑠}=`$ $`\left[\begin{array}{cccc}e^{(2\psi )}& 0& 0& 0\\ 0& {\displaystyle \frac{M^2e^{(2\gamma 2\psi )}(\lambda ^2\mu ^2)}{\lambda ^21}}& 0& 0\\ 0& 0& {\displaystyle \frac{M^2e^{(2\gamma 2\psi )}(\lambda ^2\mu ^2)}{1\mu ^2}}& 0\\ 0& 0& 0& M^2e^{(2\psi )}(\lambda ^21)(1\mu ^2)\end{array}\right]`$  $`,\mathrm{𝑖𝑛𝑑𝑒𝑥}\mathrm{\_}\mathrm{𝑐ℎ𝑎𝑟}=[1,1]`$  $`])`$

where

f1 := psi = 1/2\*(
(1+q\*(3\*lambda^2-1)\*(3\*mu^2-1)/4)\*ln((lambda-1)/\
(lambda+1))+3/2\*q\*lambda\*(3\*mu^2-1) );
f2 := gamma =\
1/2\*(1+q+q^2)\*ln((lambda^2-1)/(lambda^2-mu^2))\
-3/2\*q\*(1-mu^2)\*(lambda\*ln((lambda-1)/(lambda+1))+2)\
+9/4\*q^2\*(1-mu^2)\*((lambda^2+mu^2-1-9\*lambda^2\*mu^2)\*\
(lambda^2-1)/16\*ln((lambda-1)/(lambda+1))^2+\
(lambda^2+7\*mu^2-5/3-9\*mu^2\*lambda^2)\*lambda\*\
ln((lambda-1)/(lambda+1))/4+1/4\*lambda^2\*(1-9\*mu^2)+\
(mu^2-1/3) );
f3 := lambda = r/M-1;
f4 := mu = cos(theta);

$`\mathit{f1}:=`$ $`\psi ={\displaystyle \frac{1}{2}}(1+{\displaystyle \frac{1}{4}}q(3\lambda ^21)(3\mu ^21))\mathrm{ln}({\displaystyle \frac{\lambda 1}{\lambda +1}})+{\displaystyle \frac{3}{4}}q\lambda (3\mu ^21)`$

$`\mathit{f2}:=\gamma ={\displaystyle \frac{1}{2}}(1+q+q^2)\mathrm{ln}({\displaystyle \frac{\lambda ^21}{\lambda ^2\mu ^2}})`$ $`{\displaystyle \frac{3}{2}}q(1\mu ^2)(\lambda \mathrm{ln}({\displaystyle \frac{\lambda 1}{\lambda +1}})+2)+{\displaystyle \frac{9}{4}}q^2(1\mu ^2)(`$  $`{\displaystyle \frac{1}{16}}(\lambda ^2+\mu ^219\lambda ^2\mu ^2)(\lambda ^21)\mathrm{ln}({\displaystyle \frac{\lambda 1}{\lambda +1}})^2`$  $`+{\displaystyle \frac{1}{4}}(\lambda ^2+7\mu ^2{\displaystyle \frac{5}{3}}9\lambda ^2\mu ^2)\lambda \mathrm{ln}({\displaystyle \frac{\lambda 1}{\lambda +1}})+`$  $`{\displaystyle \frac{1}{4}}\lambda ^2(19\mu ^2)+\mu ^2{\displaystyle \frac{1}{3}})`$

$$\mathit{f3}:=\lambda =\frac{r}{M}1$$

$$\mathit{f4}:=\mu =\mathrm{cos}(\theta )$$

In the case q=0 we have

get\_compts(er):
map2(subs,{psi=rhs(f1),gamma=rhs(f2)},%):
map2(subs,q=0,%):
map2(subs,{lambda=rhs(f3),mu=rhs(f4)},%):
map(simplify,%);

$$\left[\begin{array}{cccc}\frac{r2M}{r}& 0& 0& 0\\ 0& \frac{rM^2}{r2M}& 0& 0\\ 0& 0& \frac{r^2}{1+\mathrm{cos}(\theta )^2}& 0\\ 0& 0& 0& r^2(1+\mathrm{cos}(\theta )^2)\end{array}\right]$$

That is the Schwarzschild metric with regard to f3 and f4.

Now let’s find the horizon in the general case of nonzero quadrupole moment. For static field ( $`\frac{}{t}g_{\alpha ,\beta }`$ =0) the corresponding condition is $`g_{0,\mathrm{\hspace{0.17em}0}}`$ =0. Then

map2(subs,psi=rhs(f1),er):
er2 := map2(subs,gamma=rhs(f2),%):
get\_compts(er2);

$$e^{((1+1/4q(3\lambda ^21)(3\mu ^21))\mathrm{ln}(\frac{\lambda 1}{\lambda +1})+3/2q\lambda (3\mu ^21))}$$

That is r=2M ( $`\lambda `$ =1). But the result for invariant of curvature is:

erinv := invert( er2, ’detg’ ):
D1er := d1metric( er2, coord ):
D2er := d2metric( D1er, coord ):
Cf1 := Christoffel1 ( D1er ):
RMN := Riemann( erinv, D2er, Cf1 ):
raise(erinv,RMN,1):#raise of indexes\
in Riemann tensor
raise(erinv,%,2):
raise(erinv,%,3):
RMNinv := raise(erinv,%,4):
prod(RMN,RMNinv,,,,):
get\_compts(%):
series(%,q=0,2):#expansion on q
convert(%,polynom):
res := simplify(%):

In the spherical case the result corresponds to above obtained:

factor( subs(q=0,res) ):#spherical symmetry
subs(lambda=1,%);# this is 48\*M^2/r^6

$$\frac{3}{4}\frac{1}{M^4}$$

There is no singularity. But for nonzero q (let’s choose $`\mu `$ = 0 for sake of simplification):

subs(mu=0,res):
#L’Hospital’s rule for calculation of limit
limit(diff(numer(%),lambda),lambda=1);
limit(diff(denom(%),lambda),lambda=1);

$$\mathrm{signum}(q)\mathrm{}$$

$$0$$

Hence, there is a true singularity on horizon that can be regarded as the demonstration of the impossibility of static axially symmetric black hole with nonzero quadrupole momentum. Such momentum will be ”taken away” by the gravitational waves in the process of the black hole formation.

### 3.6 Reissner-Nordström black hole (charged black hole)

The generalization of Schwarzschild metric on the case of spherically symmetric vacuum solution of bounded Einstein-Maxwell equations results in

coord := \[t, r, theta, phi\]:
rn\_compts := \
array(symmetric,sparse,1..4,1..4):# metric components
rn\_compts :=\
-(1-2\*M/r+Q^2/r^2):# coefficient of d(t)^2 in interval
rn\_compts :=\
1/(1-2\*M/r+Q^2/r^2):# coefficient of d(r)^2 in interval
rn\_compts := \
g\_matrix:# coefficient of d(theta)^2 in interval
rn\_compts :=\
g\_matrix:# coefficient of d(phi)^2 in interval
rn :=\
create(\[-1,-1\], eval(rn\_compts));# Reissner-Nordstrom\
(RN) metric

$`\mathrm{𝑟𝑛}:=`$ $`\mathrm{table}([\mathrm{𝑐𝑜𝑚𝑝𝑡𝑠}=\left[\begin{array}{cccc}1+{\displaystyle \frac{2M}{r}}{\displaystyle \frac{Q^2}{r^2}}& 0& 0& 0\\ 0& {\displaystyle \frac{1}{1{\displaystyle \frac{2M}{r}}+{\displaystyle \frac{Q^2}{r^2}}}}& 0& 0\\ 0& 0& r^2& 0\\ 0& 0& 0& r^2\mathrm{sin}(\theta )^2\end{array}\right],`$  $`\mathrm{𝑖𝑛𝑑𝑒𝑥}\mathrm{\_}\mathrm{𝑐ℎ𝑎𝑟}=[1,1]`$  $`])`$

Here Q is the electric charge. The metric has three singularities: r=0 and

denom(get\_compts(rn)) = 0;
r\_p := solve(%, r);
r\_n := solve(%%, r);

$$r^22rM+Q^2=0$$

$$r\mathrm{\_}p:=M+\sqrt{M^2Q^2}$$

$$r\mathrm{\_}n:=M\sqrt{M^2Q^2}$$

Let’s calculate the invariant of curvature:

rninv := invert( rn, ’detg’ ):
D1rn := d1metric( rn, coord ):
D2rn := d2metric( D1rn, coord ):
Cf1 := Christoffel1 ( D1rn ):
RMN := Riemann( rninv, D2rn, Cf1 ):
raise(rninv,RMN,1):#raise of indexes\
in Riemann tensor
raise(rninv,%,2):
raise(rninv,%,3):
RMNinv := raise(rninv,%,4):
prod(RMN,RMNinv,,,,);

$$\mathrm{table}([\mathrm{𝑐𝑜𝑚𝑝𝑡𝑠}=8\frac{6r^2M^212rMQ^2+7Q^4}{r^8},\mathrm{𝑖𝑛𝑑𝑒𝑥}\mathrm{\_}\mathrm{𝑐ℎ𝑎𝑟}=[]])$$

solve(get\_compts(%),r);#nonphysical roots for\
numerator with nonzero Q

$$\frac{(1+{\displaystyle \frac{1}{6}}I\sqrt{6})Q^2}{M},\frac{(1{\displaystyle \frac{1}{6}}I\sqrt{6})Q^2}{M}$$

That is the situation like to one in Schwarzschild metric and two last singularities have a coordinate character. Now let us plot the signs of two first terms in linear element (we plot inverse value for second term in order to escape the divergence due to coordinate singularities).

plot({subs({M=1,Q=1/2},get\_compts(rn)),\
subs({M=1,Q=1/2},1/get\_compts(rn))},\
r=0.1..2, title=‘signs of first and second terms of\
linear element‘);

One can see the radical difference from the Schwarzschild black hole. The space and time terms exchange the roles in region r\_n $`<`$ r $`<`$ r\_p (between so-called inner and outer horizons with $`g_{0,\mathrm{\hspace{0.17em}0}}`$ =0). But there are the usual signs in the vicinity of physical singularity, i. e. it has a time-like character and the falling observer can avoid this singularity.

Next difference is the lack of coordinate singularities for $`M^2`$ $`<`$ $`Q^2`$ . In this case we have the so-called naked singularity. One can demonstrate that there is no such singularity as result of usual collapse of charged shell with mass M and charge Q. The total energy (in Newtonian limit but with correction in framework of special relativity, M\_0 is the rest mass) is

en := M(r) = M\_0 + Q^2/r - M(r)^2/r;#we use\
the geometric units of charge so that the Coulomb low is\
G\*Q\_1\*Q\_2/r^2
sol := solve(en,M(r));

$$\mathrm{𝑒𝑛}:=\mathrm{M}(r)=M\mathrm{\_}\mathit{0}+\frac{Q^2}{r}\frac{\mathrm{M}(r)^2}{r}$$

$$\mathrm{𝑠𝑜𝑙}:=\frac{1}{2}r+\frac{1}{2}\sqrt{r^2+4M\mathrm{\_}\mathit{0}r+4Q^2},\frac{1}{2}r\frac{1}{2}\sqrt{r^2+4M\mathrm{\_}\mathit{0}r+4Q^2}$$

The choice of the solution is defined by the correct asymptotic $`lim_r\mathrm{}\mathrm{𝑠𝑜𝑙}`$ = M\_0:

if limit(sol,r=infinity)=M\_0 then true\_sol :=\
sol fi:
if limit(sol,r=infinity)=M\_0 then true\_sol :=\
sol fi:
true\_sol;

$$\frac{1}{2}r+\frac{1}{2}\sqrt{r^2+4M\mathrm{\_}\mathit{0}r+4Q^2}$$

diff(true\_sol,r):
subs(M\_0=solve(en,M\_0),%):
simplify(%,radical,symbolic);

$$\frac{Q^2+\mathrm{M}(r)^2}{(r+2\mathrm{M}(r))r}$$

So, $`\frac{}{r}\mathrm{M}(r)`$ = $`\frac{\mathrm{M}(r)^2Q^2}{r(r+2\mathrm{M}(r))}`$ . The collapse is possible if M decreases with decreasing R (domination of gravity over the Coulomb interaction in the process of collapse) that is possible, when $`M^2`$ $`>`$ $`Q^2`$ . It should be noted, that the limit

limit(true\_sol,r=0);

$$\sqrt{Q^2}$$

resolves the problem of the infinite proper energy of charged particle.

Now we will consider the pressure free collapse of charged sphere of dust by analogy with Schwarzschild metric.

r := ’r’:
E := ’E’:
subs( r=r(t),get\_compts(rn) ):
d(s)^2 = %\*d(t)^2 + %\*d(r)^2;#RN 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)+Q^2/r^2),\
d(tau)=1},%);#we used d(t)/d(tau)=E/(1-2\*M/r+Q^2/r^2)
pot\_1 :=\
factor( solve(%,(diff(r(t),t))^2) ):#"potential"\
for remote observer
pot\_2 :=\
simplify(pot\_1\*(E/(1-2\*M/r(t)+Q^2/r^2))^2):#"potential"\
for collapsing observer\
d/d(t)=(d/d(tau))\*(1-2\*M/r+Q^2/r^2)/E

$$\mathrm{d}(s)^2=(1+\frac{2M}{\mathrm{r}(t)}\frac{Q^2}{\mathrm{r}(t)^2})\mathrm{d}(t)^2+\frac{\mathrm{d}(r)^2}{1{\displaystyle \frac{2M}{\mathrm{r}(t)}}+{\displaystyle \frac{Q^2}{\mathrm{r}(t)^2}}}$$

$$\mathrm{d}(\tau )^2=\left(1+\frac{2M}{\mathrm{r}(t)}\frac{Q^2}{\mathrm{r}(t)^2}+\frac{(\frac{}{t}\mathrm{r}(t))^2}{1{\displaystyle \frac{2M}{\mathrm{r}(t)}}+{\displaystyle \frac{Q^2}{\mathrm{r}(t)^2}}}\right)\mathrm{d}(t)^2$$

$$1=\frac{\left(1+{\displaystyle \frac{2M}{\mathrm{r}(t)}}{\displaystyle \frac{Q^2}{\mathrm{r}(t)^2}}+{\displaystyle \frac{(\frac{}{t}\mathrm{r}(t))^2}{1{\displaystyle \frac{2M}{\mathrm{r}(t)}}+{\displaystyle \frac{Q^2}{\mathrm{r}(t)^2}}}}\right)\mathrm{d}(t)^2}{\mathrm{d}(\tau )^2}$$

$$1=\frac{\left(1+{\displaystyle \frac{2M}{\mathrm{r}(t)}}{\displaystyle \frac{Q^2}{\mathrm{r}(t)^2}}+{\displaystyle \frac{(\frac{}{t}\mathrm{r}(t))^2}{1{\displaystyle \frac{2M}{\mathrm{r}(t)}}+{\displaystyle \frac{Q^2}{\mathrm{r}(t)^2}}}}\right)E^2}{(1{\displaystyle \frac{2M}{\mathrm{r}(t)}}+{\displaystyle \frac{Q^2}{r^2}})^2}$$

plot({subs({E=0.5,M=1,Q=1/2,r(t)=r},pot\_2),\
0\*r},r=0.11..3,axes=boxed,title=‘(dr/dtau)^2 vs r for\
collapsing observer‘);

The difference from Schwarzschild collapse is obvious: the observer crosses the outer and inner horizons but does not reach the singularity because of the collapsar explodes as white hole due to repulsion with consequent recollapse and so on.

And at last, we consider the ”extreme” case $`M^2`$ = $`Q^2`$ .

subs(Q^2=M^2, get\_compts(rn));
factor(%);

$$\left[\begin{array}{cccc}1+\frac{2M}{r}\frac{M^2}{r^2}& 0& 0& 0\\ 0& \frac{1}{1{\displaystyle \frac{2M}{r}}+{\displaystyle \frac{M^2}{r^2}}}& 0& 0\\ 0& 0& r^2& 0\\ 0& 0& 0& r^2\mathrm{sin}(\theta )^2\end{array}\right]$$

$$\frac{(rM)^2}{r^2}$$

So, we have one coordinate singularity in r=M. What happen with second horizon? Let’s find the distance between horizons for fixed t and angular coordinates for RN-metric:

get\_compts(rn):
d(s)^2 = %\*d(r)^2;#RN metric
\# or
d(s)^2 = d(r)^2/expand( (1-r\_p/r)\*(1-r\_n/r) );#second\
representation of expression

$$\mathrm{d}(s)^2=\frac{\mathrm{d}(r)^2}{1{\displaystyle \frac{2M}{r}}+{\displaystyle \frac{Q^2}{r^2}}}$$

$$\mathrm{d}(s)^2=\frac{\mathrm{d}(r)^2}{1{\displaystyle \frac{2M}{r}}+{\displaystyle \frac{Q^2}{r^2}}}$$

Hence, when r\_p–$`>`$r\_n ( $`M^2Q^2`$ $`>`$0)

r\_p := ’r\_p’:
r\_n := ’r\_n’:
s = Int(1/sqrt((1-r\_p/r)\*(1-r\_n/r)),r=r\_n..r\_p) ;
simplify( value(rhs(%)),radical,symbolic );

$$s=_{r\mathrm{\_}n}^{r\mathrm{\_}p}\frac{1}{\sqrt{(1{\displaystyle \frac{r\mathrm{\_}p}{r}})(1{\displaystyle \frac{r\mathrm{\_}n}{r}})}}𝑑r$$

$`{\displaystyle \frac{1}{2}}r\mathrm{\_}n\mathrm{ln}(r\mathrm{\_}pr\mathrm{\_}n)+{\displaystyle \frac{1}{2}}r\mathrm{\_}p\mathrm{ln}(r\mathrm{\_}pr\mathrm{\_}n)`$ $`{\displaystyle \frac{1}{2}}r\mathrm{\_}n\mathrm{ln}(r\mathrm{\_}p+r\mathrm{\_}n){\displaystyle \frac{1}{2}}r\mathrm{\_}p\mathrm{ln}(r\mathrm{\_}p+r\mathrm{\_}n)`$

we have s–$`>`$ $`\mathrm{}`$ . So, there is the infinitely long Einstein-Rosen bridge (charged string) between horizons that means a lack of wormhole between asymptotically flat universes. This fact can be illustrated by means of embedding of equatorial section of static RN space in flat Euclidian space (see above).

d(r)^2/(1-1/r)^2 =\
(1+diff(z(r),r)^2)\*d(r)^2;#equality of radial elements\
of intervals, M=1
diff(z(r),r) = solve(%,diff(z(r),r));
dsolve(%,z(r));# embedding

$$\frac{\mathrm{d}(r)^2}{(1{\displaystyle \frac{1}{r}})^2}=(1+(\frac{}{r}\mathrm{z}(r))^2)\mathrm{d}(r)^2$$

$$\frac{}{r}\mathrm{z}(r)=\frac{\sqrt{2r1}}{r1}$$

$$\mathrm{z}(r)=2\sqrt{2r1}+\mathrm{ln}(\sqrt{2r1}1)\mathrm{ln}(\sqrt{2r1}+1)+\mathrm{\_}\mathit{C1}$$

plot3d(subs(r=sqrt(x^2+y^2),\
2\*sqrt(2\*r-1)+ln(sqrt(2\*r-1)-1)-ln(sqrt(2\*r-1)+1)),\
x=-10..10,y=-10..10,axes=boxed,style=PATCHCONTOUR,\
grid=,title=‘extreme RN black hole‘);\
#we expressed arctanh through ln

The asymptotic behavior of RN-metric as r–$`>`$ $`\mathrm{}`$  is Minkowski. For investigation of the situation r–$`>`$M let introduce the new coordinate (see, for example, ):

rn\_assym := subs(\
{r=M\*(1+lambda),Q^2=M^2},get\_compts(rn) );

$`\mathrm{𝑟𝑛}\mathrm{\_}\mathrm{𝑎𝑠𝑠𝑦𝑚}:=`$ $`\left[\begin{array}{cccc}1+{\displaystyle \frac{2}{\lambda +1}}{\displaystyle \frac{1}{(\lambda +1)^2}}& 0& 0& 0\\ 0& {\displaystyle \frac{1}{1{\displaystyle \frac{2}{\lambda +1}}+{\displaystyle \frac{1}{(\lambda +1)^2}}}}& 0& 0\\ 0& 0& M^2(\lambda +1)^2& 0\\ 0& 0& 0& M^2(\lambda +1)^2\mathrm{sin}(\theta )^2\end{array}\right]`$

m1 :=\
series(rn\_assym,lambda=0,3);# we keep only\
leading term in lambda
m2 :=\
series(rn\_assym,lambda=0,3);# we keep only\
leading term in lambda
d(s)^2 = convert(m1,polynom)\*d(t)^2\
\+ M^2\*convert(m2,polynom)\*d(lambda)^2 +\
M^2\*d(Omega)^2;#d(Omega) is spherical part

$$\mathit{m1}:=\lambda ^2+\mathrm{O}(\lambda ^3)$$

$$\mathit{m2}:=\lambda ^2+\mathrm{O}(\lambda ^1)$$

$$\mathrm{d}(s)^2=\lambda ^2\mathrm{d}(t)^2+\frac{M^2\mathrm{d}(\lambda )^2}{\lambda ^2}+M^2\mathrm{d}(\mathrm{\Omega })^2$$

This is the Robinson-Bertotti metric. The last term describes two - dimensional sphere with radius M (these dimensions are compactified in the vicinity of horizon) and the first terms corresponds to anti-de Sitter space-time (see below) with constant negative curvature.

### 3.7 Kerr black hole (rotating black hole)

In the general form the stationary rotating black hole is described by so-called Kerr-Newman metric, which in the Boyer-Linquist coordinates can be presented as:

coord := \[t, r, theta, phi\]:
kn\_compts :=\
array(sparse,1..4,1..4):# metric components,\
a=J/M, J is the angular momentum
kn\_compts :=\
-(Delta-a^2\*sin(theta)^2)/Sigma:#coefficient of d(t)^2
kn\_compts :=\
-2\*a\*sin(theta)^2\*(r^2+a^2-Delta)/\
Sigma:# coefficient of d(t)\*dphi
kn\_compts :=\
Sigma/Delta:# coefficient of d(r)^2
kn\_compts :=\
Sigma:# coefficient of d(theta)^2
kn\_compts :=\
(((r^2+a^2)^2-Delta\*a^2\*sin(theta)^2)/\
Sigma)\*sin(theta)^2:#coefficient of d(phi)^2
kn := create(\[-1,-1\], eval(kn\_compts));# Kerr-Newman\
(KN) metric
#where
sub\_1 := Sigma = r^2+a^2\*cos(theta)^2;
sub\_2 := Delta = r^2-2\*M\*r+a^2+sqrt(Q^2+P^2);# P is\
the magnetic (monopole) charge

$`\mathrm{𝑘𝑛}:=\mathrm{table}([\mathrm{𝑐𝑜𝑚𝑝𝑡𝑠}=`$ $`\left[\begin{array}{cccc}{\displaystyle \frac{\mathrm{\Delta }a^2\mathrm{sin}(\theta )^2}{\mathrm{\Sigma }}}& 0& 0& 2{\displaystyle \frac{a\mathrm{sin}(\theta )^2(r^2+a^2\mathrm{\Delta })}{\mathrm{\Sigma }}}\\ 0& {\displaystyle \frac{\mathrm{\Sigma }}{\mathrm{\Delta }}}& 0& 0\\ 0& 0& \mathrm{\Sigma }& 0\\ 0& 0& 0& {\displaystyle \frac{((r^2+a^2)^2\mathrm{\Delta }a^2\mathrm{sin}(\theta )^2)\mathrm{sin}(\theta )^2}{\mathrm{\Sigma }}}\end{array}\right],`$  $`\mathrm{𝑖𝑛𝑑𝑒𝑥}\mathrm{\_}\mathrm{𝑐ℎ𝑎𝑟}=[1,1]`$  $`])`$

$$\mathrm{𝑠𝑢𝑏}\mathrm{\_}\mathit{1}:=\mathrm{\Sigma }=r^2+a^2\mathrm{cos}(\theta )^2$$

$$\mathrm{𝑠𝑢𝑏}\mathrm{\_}\mathit{2}:=\mathrm{\Delta }=r^22rM+a^2+\sqrt{Q^2+P^2}$$

In the absence of charges, this results in Kerr metric. The obvious singularities are (except for an usual singularity of spherical coordinates $`\theta `$ =0):

r\_p :=\
solve( subs({Q=0,P=0},\
rhs(sub\_2))=0,r );#outer horizon
r\_n :=\
solve( subs({Q=0,P=0},\
rhs(sub\_2))=0,r );#inner horizon
solve( subs({Q=0,P=0},rhs(sub\_1))=0,theta );

$$r\mathrm{\_}p:=M+\sqrt{M^2a^2}$$

$$r\mathrm{\_}n:=M\sqrt{M^2a^2}$$

$$\frac{1}{2}\pi I\mathrm{arcsinh}(\frac{r}{a}),\frac{1}{2}\pi +I\mathrm{arcsinh}(\frac{r}{a})$$

The last produces r=0, $`\theta `$ = $`\frac{\pi }{2}`$ .

As it was in the case of charged static black hole, there are three different situations: $`M^2`$ $`<`$ $`a^2`$ , $`M^2`$ = $`a^2`$ , $`M^2`$ $`>`$ $`a^2`$ .

Let us consider the signs of $`g_{0,\mathrm{\hspace{0.17em}0}}`$ , $`g_{1,\mathrm{\hspace{0.17em}1}}`$ , $`g_{4,\mathrm{\hspace{0.17em}4}}`$ ( $`M^2`$ $`>`$ $`a^2`$ ).

plot3d(subs({M=1,a=1/2,P=0,Q=0},subs(\
{Sigma=rhs(sub\_1),Delta=rhs(sub\_2)},\
get\_compts(kn))),r=0.1..4,theta=0..Pi,color=red):
plot3d(subs({M=1,a=1/2,P=0,Q=0,theta=2\*Pi/3},\
subs({Sigma=rhs(sub\_1),Delta=rhs(sub\_2)},\
1/get\_compts(kn))),\
r=0.1..4,theta=0..Pi,color=green):
display(%,%%,\
title=‘signs of first and second diagonal elements\
of metric‘,axes=boxed);
plot3d(subs({M=1,a=1/2,P=0,Q=0},subs(\
{Sigma=rhs(sub\_1),Delta=rhs(sub\_2)},\
get\_compts(kn))),r=-0.01..0.01,\
theta=Pi/2-0.01..Pi/2+0.01,color=blue,\
title=‘four diagonal element of metric in vicinity of\
singularity‘,axes=boxed);

One can see, that the approach to r=0 in the line, which differs from $`\theta `$ = $`\frac{\pi }{2}`$, corresponds to usual signs of diagonal elements of metric, i.e. the observer crosses r=0 and comes into region r$`<`$0 without collision with singularity. But from the second picture we can see the change of the $`g_{3,\mathrm{\hspace{0.17em}3}}`$ -sign for r$`<`$0. Now $`\varphi `$ is the time like coordinate. But this coordinate has circle character and, as consequence, we find oneself in the world with closed time lines. The approach to r=0 in the line of $`\theta `$ = $`\frac{\pi }{2}`$ produces the change of $`g_{0,\mathrm{\hspace{0.17em}0}}`$ -sign, i.e. we find a true singularity in this direction. These facts demonstrate that the singularity in Kerr black hole has a more complicated character than in above considered black holes.

The more careful consideration gives the following results:

1) $`M^2`$ $`<`$ $`a^2`$ . There exist no horizon (r\_p and r\_n are complex), but the singularity in r=0, $`\theta `$ = $`\frac{\pi }{2}`$ keeps. To remove the coordinate singularity in $`\theta `$ =0 we introduce the Kerr-Schild coordinates with linear element

macro(ts=‘t\*‘):
ks\_le := d(s)^2 = -d(ts)^2 + d(x)^2 + d(y)^2\
\+ d(z)^2 + (2\*M\*r^3/(r^4+a^2\*z^2))\*\
((r\*(x\*d(x)+y\*d(y))-a\*(x\*d(y)-y\*d(x)))/\
(r^2+a^2)+z\*d(z)/r+d(ts) )^2;
x + I\*y = (r + I\*a)\*sin(theta)\*\
exp(I\*(Int(1,phi)+Int(a/Delta,r)));
z = r\*cos(theta);
ts = Int(1,t) + Int((r^2+a^2/Delta),r) - r;

$`\mathrm{𝑘𝑠}\mathrm{\_}\mathrm{𝑙𝑒}:=\mathrm{d}(s)^2=\mathrm{d}(t)^2+\mathrm{d}(x)^2+\mathrm{d}(y)^2+\mathrm{d}(z)^2`$ $`+{\displaystyle \frac{2Mr^3({\displaystyle \frac{r(x\mathrm{d}(x)+y\mathrm{d}(y))a(x\mathrm{d}(y)y\mathrm{d}(x))}{r^2+a^2}}+{\displaystyle \frac{z\mathrm{d}(z)}{r}}+\mathrm{d}(t))^2}{r^4+a^2z^2}}`$

$$x+Iy=(r+Ia)\mathrm{sin}(\theta )e^{(I({\scriptscriptstyle 1𝑑\varphi }+{\scriptscriptstyle {\scriptscriptstyle \frac{a}{\mathrm{\Delta }}}𝑑r}))}$$

$$z=r\mathrm{cos}(\theta )$$

$$t=1dt+r^2+\frac{a^2}{\mathrm{\Delta }}drr$$

which is reduced to Minkowski metric by M–$`>`$0.

int(subs( {P=0,Q=0},\
subs(Delta=rhs(sub\_2),a/Delta) ),r):
x+I\*y=(r+I\*a)\*sin(theta)\*exp(I\*(phi+%));

$$x+Iy=(r+Ia)\mathrm{sin}(\theta )e^{(I(\varphi +\frac{a\mathrm{arctan}(1/2\frac{2r2M}{\sqrt{a^2M^2}})}{\sqrt{a^2M^2}}))}$$

When r=0, $`\theta `$ = $`\frac{\pi }{2}`$ the singularity is the ring $`x^2+y^2`$ = $`a^2`$ , z = 0.

2) $`M^2`$ $`>`$ $`a^2`$ . As before we have the ring singularity, but there are the horizons r\_p and r\_n. As an additional feature of Kerr metric we note here the existence of coordinate singularity:

get\_compts(kn)=0;
subs( Delta=rhs( sub\_2 ),numer( lhs(%) ) ):
solve( subs({Q=0,P=0},%) = 0,r );

$$\frac{\mathrm{\Delta }a^2\mathrm{sin}(\theta )^2}{\mathrm{\Sigma }}=0$$

$$M+\sqrt{M^2a^2+a^2\mathrm{sin}(\theta )^2},M\sqrt{M^2a^2+a^2\mathrm{sin}(\theta )^2}$$

The crossing of ellipsoid r_1= $`M+\sqrt{M^2a^2\mathrm{cos}(\theta )^2}`$ produces the change of $`g_{0,\mathrm{\hspace{0.17em}0}}`$ -sign. As it was for Schwarzschild black hole this fact demonstrates the lack of static coordinates under this surface, which is called ”ergosphere” and the region between r\_p and r\_1 is the ergoregion. The absence of singularity for $`r_{r,r}`$ suggests that the nonstatic behavior results from entrainment of observer by black hole rotation, but not from fall on singularity, as it takes place for observer under horizon in Schwarzschild black hole.

### 3.8 Conclusion

So, the elementary analysis by means of basic Maple 6 functions allows to obtain the main results of black hole physics including conditions of collapsar formation, the space-time structure of static spherically symmetric charged and uncharged black holes and stationary axially symmetric black hole.

## 4 Cosmological models

### 4.1 Introduction

I suppose that the most impressive achievements of the GR-theory belong to the cosmology. The description of the global structure of the universe changes our opinion about reality and radically widens the horizon of knowledge. The considerable progress in the observations and measures allows to say that we are living in a gold age of cosmology. Here we will consider some basic conceptions of relativistic cosmology by means of analytical capabilities of Maple 6.

### 4.2 Robertson-Walker metric

We will base on the conception of foliation of 4-dimensional manifold on 3-dimensional space-like hyperplanes with a time-dependent geometry. Now let’s restrict oneself to the isotropic and homogeneous evolution models, i.e. the models without dependence of curvature on the observer’s location or orientation. We will suppose also, that such foliation with isotropic and homogeneous geometry and energy-momentum distribution exists at each time moment. These space-like foliations are the so-called hyperplanes of homogeneity. The time vectors are defined by the proper time course in each space point.

Let us consider 3-geometry at fixed time moment. As it is known, in 3-dimensional space the curvature tensor has the following form:

$`R_{\mathrm{𝑖𝑘𝑙𝑚}}`$ =

$`R_{\mathrm{𝑖𝑙}}`$ $`g_{\mathrm{𝑘𝑚}}`$  \- $`R_{\mathrm{𝑖𝑚}}`$ $`g_{\mathrm{𝑘𝑙}}`$  \+ $`R_{\mathrm{𝑘𝑚}}`$ $`g_{\mathrm{𝑖𝑙}}`$  \- $`R_{\mathrm{𝑘𝑙}}`$ $`g_{\mathrm{𝑖𝑚}}`$  \+ $`\frac{R}{2}`$ ( $`g_{\mathrm{𝑖𝑚}}`$ $`g_{\mathrm{𝑘𝑙}}`$ - $`g_{\mathrm{𝑖𝑙}}`$ $`g_{\mathrm{𝑘𝑚}}`$ ).

But if (as it takes a place in our case) there is not the dependence of curvature on space-direction in arbitrary point, our space has a constant curvature (Schur’s theorem) that results in

$`R_{\mathrm{𝑖𝑘𝑙𝑚}}`$ = $`\frac{R}{2}`$ ( $`g_{\mathrm{𝑖𝑚}}`$ $`g_{\mathrm{𝑘𝑙}}`$ - $`g_{\mathrm{𝑖𝑙}}`$ $`g_{\mathrm{𝑘𝑚}}`$ ), (1)

where R is the constant Ricci scalar.

Now we will consider the so-called Robertson-Walker metric:

$`g_{\mathrm{𝑖𝑗}}`$ = $`\frac{\delta _{\mathrm{𝑖𝑗}}}{(1+\frac{R\delta _{\mathrm{𝑖𝑗}}x^ix^j}{8})^2}`$ (2)

For this metric the spatial curvature tensor is:

restart:
with(tensor):
with(linalg):
with(difforms):
coord := \[x, y, z\]:# coordinates
g\_compts :=\
array(symmetric,sparse,1..3,1..3):# metric components
g\_compts := 1/\
(1+R\*(x^2+y^2+z^2)/8)^2:# component of interval\
attached to d(x)^2
g\_compts :=\
1/(1+R\*(x^2+y^2+z^2)/8)^2:# component of interval\
attached to d(y)^2
g\_compts :=\
1/(1+R\*(x^2+y^2+z^2)/8)^2:# component of\
interval attached to d(z)^2
g := create(\[-1,-1\], eval(g\_compts));# covariant\
metric tensor
ginv := invert( g, ’detg’ ):# contravariant\
metric tensor
D1g := d1metric( g, coord ):# calculation of \
curvature tensor
D2g := d2metric( D1g, coord ):
Cf1 := Christoffel1 ( D1g ):
RMN := Riemann( ginv, D2g, Cf1 ):
displayGR(Riemann, RMN);

$`g:=\mathrm{table}([\mathrm{𝑖𝑛𝑑𝑒𝑥}\mathrm{\_}\mathrm{𝑐ℎ𝑎𝑟}=[1,1],\mathrm{𝑐𝑜𝑚𝑝𝑡𝑠}=\left[\begin{array}{ccc}\%1& 0& 0\\ 0& \%1& 0\\ 0& 0& \%1\end{array}\right]])`$ $`\%1:={\displaystyle \frac{1}{(1+{\displaystyle \frac{1}{8}}R(x^2+y^2+z^2))^2}}`$

$$\mathrm{𝑇ℎ𝑒}\mathrm{𝑅𝑖𝑒𝑚𝑎𝑛𝑛}\mathrm{𝑇𝑒𝑛𝑠𝑜𝑟}$$

$$\mathrm{𝑛𝑜𝑛}\mathrm{𝑧𝑒𝑟𝑜}\mathrm{𝑐𝑜𝑚𝑝𝑜𝑛𝑒𝑛𝑡𝑠}:$$

$$\mathit{R1212}=2048\frac{R}{(8+Rx^2+Ry^2+Rz^2)^4}$$

$$\mathit{R1313}=2048\frac{R}{(8+Rx^2+Ry^2+Rz^2)^4}$$

$$\mathit{R2323}=2048\frac{R}{(8+Rx^2+Ry^2+Rz^2)^4}$$

$$\mathrm{𝑐ℎ𝑎𝑟𝑎𝑐𝑡𝑒𝑟}:[\mathit{1},\mathit{1},\mathit{1},\mathit{1}]$$

Now we try to perform the direct calculations from Eq. (1):

A := get\_compts(g):
B := array(1..3, 1..3, 1..3, 1..3):
for i from 1 to 3 do
for j from 1 to 3 do
for k from 1 to 3 do
for l from 1 to 3 do
B\[i,j,k,l\] :=\
eval(R\*(A\[i,k\]\*A\[j,l\] - A\[i,l\]\*A\[j,k\])/2):# Eq. 1.\
These calculations can be performed also by means of\
tensor\[prod\] and tensor\[permute\_indices\]
if B\[i,j,k,l\]\<\>0 then \
print(\[i,j,k,l\],B\[i,j,k,l\]);# non-zero components
else
fi:
od: od: od: od:

$$[1,\mathrm{\hspace{0.17em}2},\mathrm{\hspace{0.17em}1},\mathrm{\hspace{0.17em}2}],\frac{1}{2}\frac{R}{(1+{\displaystyle \frac{1}{8}}R(x^2+y^2+z^2))^4}$$

$$[1,\mathrm{\hspace{0.17em}2},\mathrm{\hspace{0.17em}2},\mathrm{\hspace{0.17em}1}],\frac{1}{2}\frac{R}{(1+{\displaystyle \frac{1}{8}}R(x^2+y^2+z^2))^4}$$

$$[1,\mathrm{\hspace{0.17em}3},\mathrm{\hspace{0.17em}1},\mathrm{\hspace{0.17em}3}],\frac{1}{2}\frac{R}{(1+{\displaystyle \frac{1}{8}}R(x^2+y^2+z^2))^4}$$

$$[1,\mathrm{\hspace{0.17em}3},\mathrm{\hspace{0.17em}3},\mathrm{\hspace{0.17em}1}],\frac{1}{2}\frac{R}{(1+{\displaystyle \frac{1}{8}}R(x^2+y^2+z^2))^4}$$

$$[2,\mathrm{\hspace{0.17em}1},\mathrm{\hspace{0.17em}1},\mathrm{\hspace{0.17em}2}],\frac{1}{2}\frac{R}{(1+{\displaystyle \frac{1}{8}}R(x^2+y^2+z^2))^4}$$

$$[2,\mathrm{\hspace{0.17em}1},\mathrm{\hspace{0.17em}2},\mathrm{\hspace{0.17em}1}],\frac{1}{2}\frac{R}{(1+{\displaystyle \frac{1}{8}}R(x^2+y^2+z^2))^4}$$

$$[2,\mathrm{\hspace{0.17em}3},\mathrm{\hspace{0.17em}2},\mathrm{\hspace{0.17em}3}],\frac{1}{2}\frac{R}{(1+{\displaystyle \frac{1}{8}}R(x^2+y^2+z^2))^4}$$

$$[2,\mathrm{\hspace{0.17em}3},\mathrm{\hspace{0.17em}3},\mathrm{\hspace{0.17em}2}],\frac{1}{2}\frac{R}{(1+{\displaystyle \frac{1}{8}}R(x^2+y^2+z^2))^4}$$

$$[3,\mathrm{\hspace{0.17em}1},\mathrm{\hspace{0.17em}1},\mathrm{\hspace{0.17em}3}],\frac{1}{2}\frac{R}{(1+{\displaystyle \frac{1}{8}}R(x^2+y^2+z^2))^4}$$

$$[3,\mathrm{\hspace{0.17em}1},\mathrm{\hspace{0.17em}3},\mathrm{\hspace{0.17em}1}],\frac{1}{2}\frac{R}{(1+{\displaystyle \frac{1}{8}}R(x^2+y^2+z^2))^4}$$

$$[3,\mathrm{\hspace{0.17em}2},\mathrm{\hspace{0.17em}2},\mathrm{\hspace{0.17em}3}],\frac{1}{2}\frac{R}{(1+{\displaystyle \frac{1}{8}}R(x^2+y^2+z^2))^4}$$

$$[3,\mathrm{\hspace{0.17em}2},\mathrm{\hspace{0.17em}3},\mathrm{\hspace{0.17em}2}],\frac{1}{2}\frac{R}{(1+{\displaystyle \frac{1}{8}}R(x^2+y^2+z^2))^4}$$

The direct calculation on the basis of Eq. (1) and metric (2) results in the tensor with symmetries \[i,j,k,l\] = \[k,l,i,j\] = \[j,i,k,l\] = -\[ i,j,l,k\] and \[i,j,k,l\] + \[i,l,j,k\] + \[i,k,l,j\] = 0. These symmetries and values of nonzero components correspond to curvature tensor, which was calculated previously from (2). As consequence, metric (2) satisfies to condition of constant curvature (Schur’s theorem, Eq. (1)). We will base our further calculations on this metric.

The next step is the choice of the appropriate coordinates. We consider two types of coordinates: spherical and pseudospherical.

In the case of the spherical coordinates the transition of Eq. (2) to a new basis results in:

g\_sp :=\
simplify(subs(\
{x=r\*cos(phi)\*sin(theta),y=r\*sin(phi)\*sin(theta),\
z=r\*cos(theta)},get\_compts(g))):# This is transformation\
of coordinates in metric tensor
Sp\_compts := array(1..3,1..3):# matrix of\
conversion from (NB!) spherical coordinates
Sp\_compts := cos(phi)\*sin(theta):
Sp\_compts := sin(phi)\*sin(theta):
Sp\_compts := cos(theta):
Sp\_compts := r\*cos(phi)\*cos(theta):
Sp\_compts := r\*sin(phi)\*cos(theta):
Sp\_compts := -r\*sin(theta):
Sp\_compts := -r\*sin(phi)\*sin(theta):
Sp\_compts := r\*cos(phi)\*sin(theta):
Sp\_compts := 0:
Sp := eval(Sp\_compts):
simplify(\
multiply(transpose(Sp), g\_sp, Sp));#transition\
to new coordinates in metric

$`\left[\begin{array}{ccc}64{\displaystyle \frac{1}{\%1}}& 0& 0\\ 0& 64{\displaystyle \frac{r^2}{\%1}}& 0\\ 0& 0& 64{\displaystyle \frac{r^2(1+\mathrm{cos}(\theta )^2)}{\%1}}\end{array}\right]`$ $`\%1:=(8+Rr^2\mathrm{cos}(\varphi )^2\mathrm{sin}(\theta )^2+Rr^2\mathrm{sin}(\varphi )^2\mathrm{sin}(\theta )^2+Rr^2\mathrm{cos}(\theta )^2)^2`$

The denominator is

denom(%):
expand(%):
simplify(%,trig):
factor(%);

$$(Rr^2+8)^2$$

After substitution r\* = r $`\sqrt{\frac{R}{2}}`$  (we will omit asterix and suppose R $`>`$ 0), the metric is converted as:

coord := \[r, theta, phi\]:# spherical coordinates
g\_compts\_sp :=\
array(symmetric,sparse,1..3,1..3):# metric components
g\_compts\_sp :=\
(2/R)/(1+r^2/4)^2:# component of interval\
attached to d(r)^2
g\_compts\_sp :=\
(2/R)\*r^2/(1+r^2/4)^2:# component of interval\
attached to d(theta)^2
g\_compts\_sp :=\
(2/R)\*r^2\*sin(theta)^2/(1+r^2/4)^2:# component of\
interval attached to d(phi)^2
g\_sp := create(\[-1,-1\], eval(g\_compts\_sp));# covariant\
metric tensor

$`g\mathrm{\_}\mathrm{𝑠𝑝}:=\mathrm{table}([\mathrm{𝑖𝑛𝑑𝑒𝑥}\mathrm{\_}\mathrm{𝑐ℎ𝑎𝑟}=[1,1],`$ $`\mathrm{𝑐𝑜𝑚𝑝𝑡𝑠}=\left[\begin{array}{ccc}2{\displaystyle \frac{1}{R(1+{\displaystyle \frac{1}{4}}r^2)^2}}& 0& 0\\ 0& 2{\displaystyle \frac{r^2}{R(1+{\displaystyle \frac{1}{4}}r^2)^2}}& 0\\ 0& 0& 2{\displaystyle \frac{r^2\mathrm{sin}(\theta )^2}{R(1+{\displaystyle \frac{1}{4}}r^2)^2}}\end{array}\right]`$  $`])`$

The denominator’s form suggests the transformation of radial coordinate: sin( $`\chi `$ ) = r/(1 + $`r^2`$ /4). Then for $`\mathrm{𝑑𝑟}^2`$ -component we have

diff(sin(chi),chi)\*d(chi) =\
simplify(diff(r/(1+r^2/4),r))\*d(r);
lhs(%)^2 = rhs(%)^2;
subs(cos(chi)^2=1-r^2/(1+r^2/4)^2,lhs(%))\
= rhs(%):
factor(%);

$$\mathrm{cos}(\chi )\mathrm{d}(\chi )=4\frac{(4+r^2)\mathrm{d}(r)}{(4+r^2)^2}$$

$$\mathrm{cos}(\chi )^2\mathrm{d}(\chi )^2=16\frac{(4+r^2)^2\mathrm{d}(r)^2}{(4+r^2)^4}$$

$$\frac{(r2)^2(r+2)^2\mathrm{d}(\chi )^2}{(4+r^2)^2}=16\frac{(r2)^2(r+2)^2\mathrm{d}(r)^2}{(4+r^2)^4}$$

that leads to a new metric (d $`\chi `$  = $`\frac{\mathrm{𝑑𝑟}}{(1+\frac{r^2}{4})^2}`$):

coord := \[r, theta, phi\]:# spherical coordinates
g\_compts\_sp :=\
array(symmetric,sparse,1..3,1..3):# metric components
g\_compts\_sp := \
2/R:# component of interval attached to d(chi)^2
g\_compts\_sp :=\
(2/R)\*sin(chi)^2:# component of interval\
attached to d(theta)^2
g\_compts\_sp :=\
(2/R)\*sin(chi)^2\*sin(theta)^2:# component of\
interval attached to d(phi)^2
g\_sp :=\
create(\[-1,-1\], eval(g\_compts\_sp));# covariant\
metric tensor

$`g\mathrm{\_}\mathrm{𝑠𝑝}:=\mathrm{table}([\mathrm{𝑖𝑛𝑑𝑒𝑥}\mathrm{\_}\mathrm{𝑐ℎ𝑎𝑟}=[1,1],\mathrm{𝑐𝑜𝑚𝑝𝑡𝑠}=`$ $`\left[\begin{array}{ccc}2{\displaystyle \frac{1}{R}}& 0& 0\\ 0& 2{\displaystyle \frac{\mathrm{sin}(\chi )^2}{R}}& 0\\ 0& 0& 2{\displaystyle \frac{\mathrm{sin}(\chi )^2\mathrm{sin}(\theta )^2}{R}}\end{array}\right]])`$

Here R/2 is the Gaussian curvature K (K$`>`$0) so that the linear element is:

l\_sp := d(s)^2 = a^2\*(d(chi)^2 +\
sin(chi)^2\*(d(theta)^2+sin(theta)^2\*d(phi)^2));

$$l\mathrm{\_}\mathrm{𝑠𝑝}:=\mathrm{d}(s)^2=a^2(\mathrm{d}(\chi )^2+\mathrm{sin}(\chi )^2(\mathrm{d}(\theta )^2+\mathrm{sin}(\theta )^2\mathrm{d}(\varphi )^2))$$

Here a = $`\sqrt{\frac{1}{K}}`$ is the scaling factor. Note, that substitution r = a $`\mathrm{sin}(\chi )`$  produces the alternative form of metric:

subs(\
{sin(chi)=r/a,d(chi)^2=d(r)^2/(a^2\*(1-r^2/a^2))},\
l\_sp):#we use d(r)=a\*cos(chi)\*d(chi)
expand(%);

$$\mathrm{d}(s)^2=\frac{\mathrm{d}(r)^2}{1{\displaystyle \frac{r^2}{a^2}}}+r^2\mathrm{d}(\theta )^2+r^2\mathrm{sin}(\theta )^2\mathrm{d}(\varphi )^2$$

For negative curvature:

coord := \[r, theta, phi\]:# spherical coordinates
g\_compts\_hyp :=\
array(symmetric,sparse,1..3,1..3):# metric components
g\_compts\_hyp :=\
(-2/R)/(1-r^2/4)^2:# component of interval\
attached to d(r)^2
g\_compts\_hyp :=\
(-2/R)\*r^2/(1-r^2/4)^2:# component of interval\
attached to d(theta)^2
g\_compts\_hyp :=\
(-2/R)\*r^2\*sin(theta)^2/(1-r^2/4)^2:#component of\
interval attached to d(phi)^2
g\_hyp :=\
create(\[-1,-1\], eval(g\_compts\_hyp));# covariant metric\
tensor

$`g\mathrm{\_}\mathrm{ℎ𝑦𝑝}:=\mathrm{table}([\mathrm{𝑖𝑛𝑑𝑒𝑥}\mathrm{\_}\mathrm{𝑐ℎ𝑎𝑟}=[1,1],\mathrm{𝑐𝑜𝑚𝑝𝑡𝑠}=`$ $`\left[\begin{array}{ccc}2{\displaystyle \frac{1}{R(1{\displaystyle \frac{1}{4}}r^2)^2}}& 0& 0\\ 0& 2{\displaystyle \frac{r^2}{R(1{\displaystyle \frac{1}{4}}r^2)^2}}& 0\\ 0& 0& 2{\displaystyle \frac{r^2\mathrm{sin}(\theta )^2}{R(1{\displaystyle \frac{1}{4}}r^2)^2}}\end{array}\right]`$  $`])`$

where R is the module of Ricci scalar and new radial coordinate is I r $`\sqrt{\frac{R}{2}.}`$  Let’s try to use the substitution sh( $`\chi `$ ) = r/(1- $`r^2`$ /4). Then for $`\mathrm{𝑑𝑟}^2`$ -component we have

diff(sinh(chi),chi)\*d(chi) =\
simplify(diff(r/(1-r^2/4),r))\*d(r);
lhs(%)^2 = rhs(%)^2;
subs(cosh(chi)^2=1+r^2/(1-r^2/4)^2,lhs(%))\
= rhs(%):
factor(%);

$$\mathrm{cosh}(\chi )\mathrm{d}(\chi )=4\frac{(4+r^2)\mathrm{d}(r)}{(4+r^2)^2}$$

$$\mathrm{cosh}(\chi )^2\mathrm{d}(\chi )^2=16\frac{(4+r^2)^2\mathrm{d}(r)^2}{(4+r^2)^4}$$

$$\frac{(4+r^2)^2\mathrm{d}(\chi )^2}{(r2)^2(r+2)^2}=16\frac{(4+r^2)^2\mathrm{d}(r)^2}{(r2)^4(r+2)^4}$$

And finally

coord := \[r, theta, phi\]:# spherical coordinates
g\_compts\_hyp :=\
array(symmetric,sparse,1..3,1..3):# metric components
g\_compts\_hyp :=\
-2/R:# component of interval\
attached to d(chi)^2
g\_compts\_hyp :=\
(-2/R)\*sinh(chi)^2:# component of interval\
attached to d(theta)^2
g\_compts\_hyp :=\
(-2/R)\*sinh(chi)^2\*sin(theta)^2:# component\
of interval attached to d(phi)^2
g\_hyp :=\
create(\[-1,-1\], eval(g\_compts\_hyp));# covariant\
metric tensor

$`g\mathrm{\_}\mathrm{ℎ𝑦𝑝}:=`$ $`\mathrm{table}([\mathrm{𝑖𝑛𝑑𝑒𝑥}\mathrm{\_}\mathrm{𝑐ℎ𝑎𝑟}=[1,1],\mathrm{𝑐𝑜𝑚𝑝𝑡𝑠}=`$  $`\left[\begin{array}{ccc}2{\displaystyle \frac{1}{R}}& 0& 0\\ 0& 2{\displaystyle \frac{\mathrm{sinh}(\chi )^2}{R}}& 0\\ 0& 0& 2{\displaystyle \frac{\mathrm{sinh}(\chi )^2\mathrm{sin}(\theta )^2}{R}}\end{array}\right]])`$

Here R/2 is the Gaussian curvature K  (K$`<`$0). So, the linear element is

l\_hyp := d(s)^2 = a^2\*(d(chi)^2 +\
sinh(chi)^2\*(d(theta)^2+sin(theta)^2\*d(phi)^2));

$$l\mathrm{\_}\mathrm{ℎ𝑦𝑝}:=\mathrm{d}(s)^2=a^2(\mathrm{d}(\chi )^2+\mathrm{sinh}(\chi )^2(\mathrm{d}(\theta )^2+\mathrm{sin}(\theta )^2\mathrm{d}(\varphi )^2))$$

where a = $`\sqrt{\frac{1}{K}}`$ is the imaginary scale factor. The substitution r = a sinh( $`\chi `$ ) results in:

subs(\
{sinh(chi)=r/a,d(chi)^2=d(r)^2/(a^2\*(1+r^2/a^2))},\
l\_hyp):#we use d(r)=a\*cosh(chi)\*d(chi)
expand(%);

$$\mathrm{d}(s)^2=\frac{\mathrm{d}(r)^2}{1+{\displaystyle \frac{r^2}{a^2}}}+r^2\mathrm{d}(\theta )^2+r^2\mathrm{sin}(\theta )^2\mathrm{d}(\varphi )^2$$

At last, for flat space:

l\_flat := d(s)^2 = a^2\*(d(chi)^2 +\
chi^2\*(d(theta)^2+sin(theta)^2\*d(phi)^2));

$$l\mathrm{\_}\mathrm{𝑓𝑙𝑎𝑡}:=\mathrm{d}(s)^2=a^2(\mathrm{d}(\chi )^2+\chi ^2(\mathrm{d}(\theta )^2+\mathrm{sin}(\theta )^2\mathrm{d}(\varphi )^2))$$

Here a is the arbitrary scaling factor.

To imagine the considered geometries we can use the usual technique of embedding. In the general case the imbedding of 3-space is not possible if the dimension of enveloping flat space is only 4 (we have 6 functions $`g_{i,j}`$ , but the number of free parameters in 4-dimensional space is only 4). However, the constant nonzero curvature allows to embed our curved space in 4-dimensional flat space as result of next transformation of coordinates (for spherical geometry):

defform(f=0,w1=1,w2=1,w3=1,v=1,chi=0,theta=0,phi=0);
e\_1 := Dw^2 = a^2\*d(cos(chi))^2;
e\_2 := Dx^2 =\
(a\*d(sin(chi)\*sin(theta)\*cos(phi)))^2;
e\_3 := Dy^2 =\
(a\*d(sin(chi)\*sin(theta)\*sin(phi)))^2;
e\_4 := Dz^2 =\
(a\*d(sin(chi)\*cos(theta)))^2;

$$e\mathrm{\_}\mathit{1}:=\mathrm{𝐷𝑤}^2=a^2\mathrm{sin}(\chi )^2\mathrm{d}(\chi )^2$$

$`e\mathrm{\_}\mathit{2}:=\mathrm{𝐷𝑥}^2=a^2`$ $`(\mathrm{sin}(\theta )\mathrm{cos}(\varphi )\mathrm{cos}(\chi )\mathrm{d}(\chi )+\mathrm{sin}(\chi )\mathrm{cos}(\varphi )\mathrm{cos}(\theta )\mathrm{d}(\theta )`$  $`\mathrm{sin}(\chi )\mathrm{sin}(\theta )\mathrm{sin}(\varphi )\mathrm{d}(\varphi ))^2`$

$`e\mathrm{\_}\mathit{3}:=\mathrm{𝐷𝑦}^2=a^2`$ $`(\mathrm{sin}(\theta )\mathrm{sin}(\varphi )\mathrm{cos}(\chi )\mathrm{d}(\chi )+\mathrm{sin}(\chi )\mathrm{sin}(\varphi )\mathrm{cos}(\theta )\mathrm{d}(\theta )+`$  $`\mathrm{sin}(\chi )\mathrm{sin}(\theta )\mathrm{cos}(\varphi )\mathrm{d}(\varphi ))^2`$

$$e\mathrm{\_}\mathit{4}:=\mathrm{𝐷𝑧}^2=a^2(\mathrm{cos}(\theta )\mathrm{cos}(\chi )\mathrm{d}(\chi )\mathrm{sin}(\chi )\mathrm{sin}(\theta )\mathrm{d}(\theta ))^2$$

#so we have for Euclidian 4-metric:
simplify(\
rhs(e\_1) + rhs(e\_2) + rhs(e\_3) + rhs(e\_4), trig);

$`a^2(\mathrm{d}(\varphi )^2\mathrm{cos}(\chi )^2\mathrm{cos}(\theta )^2+\mathrm{d}(\theta )^2+\mathrm{d}(\varphi )^2\mathrm{d}(\theta )^2\mathrm{cos}(\chi )^2`$ $`\mathrm{d}(\varphi )^2\mathrm{cos}(\theta )^2\mathrm{d}(\varphi )^2\mathrm{cos}(\chi )^2+\mathrm{d}(\chi )^2)`$

collect(%,d(phi)^2);

$`a^2(1\mathrm{cos}(\theta )^2\mathrm{cos}(\chi )^2+\mathrm{cos}(\chi )^2\mathrm{cos}(\theta )^2)\mathrm{d}(\varphi )^2+`$ $`a^2(\mathrm{d}(\theta )^2+\mathrm{d}(\chi )^2\mathrm{d}(\theta )^2\mathrm{cos}(\chi )^2)`$

This expression can be easily converted in l\_sp, i. e. the imbedding is correct.

Since

w^2 + x^2 + y^2 + z^2 =\
simplify(\
(a\*cos(chi))^2 +(a\*sin(chi)\*sin(theta)\*cos(phi))^2 +\
(a\*sin(chi)\*sin(theta)\*sin(phi))^2 +\
(a\*sin(chi)\*cos(theta))^2);

$$w^2+x^2+y^2+z^2=a^2$$

the considered hypersurface is the 3-sphere in flat 4-space. The scale factor is the radius of the closed spherical universe.

In the case l\_hyp there is not the imbedding in flat Euclidian 4-space, but there is the imbedding in flat hyperbolical flat space with interval

-d $`w^2`$ \+ d $`x^2`$ \+ d $`y^2`$ \+ d $`z^2`$

As result we have

w^2 - x^2 - y^2 - z^2 =\
simplify(\
(a\*cosh(chi))^2 - (a\*sinh(chi)\*sin(theta)\*cos(phi))^2 -\
(a\*sinh(chi)\*sin(theta)\*sin(phi))^2 -\
(a\*sinh(chi)\*cos(theta))^2);

$$w^2x^2y^2z^2=a^2$$

This is expression for 3-hyperboloid in flat 4-space.

### 4.3 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:#dt^2
g\_compts := a(t)^2:#chi^2
g\_compts:=a(t)^2\*sin(chi)^2:#dtheta^2
g\_compts:=a(t)^2\*sin(chi)^2\*sin(theta)^2:#dphi^2
g := create(\[-1,-1\], eval(g\_compts));
ginv := invert( g, ’detg’ ):

$`g:=\mathrm{table}([\mathrm{𝑖𝑛𝑑𝑒𝑥}\mathrm{\_}\mathrm{𝑐ℎ𝑎𝑟}=[1,1],\mathrm{𝑐𝑜𝑚𝑝𝑡𝑠}=`$ $`\left[\begin{array}{cccc}\hfill 1& 0& 0& 0\\ \hfill 0& \mathrm{a}(t)^2& 0& 0\\ \hfill 0& 0& \mathrm{a}(t)^2\mathrm{sin}(\chi )^2& 0\\ \hfill 0& 0& 0& \mathrm{a}(t)^2\mathrm{sin}(\chi )^2\mathrm{sin}(\theta )^2\end{array}\right]])`$

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 );

$`G:=\mathrm{table}([\mathrm{𝑖𝑛𝑑𝑒𝑥}\mathrm{\_}\mathrm{𝑐ℎ𝑎𝑟}=[1,1],`$ $`\mathrm{𝑐𝑜𝑚𝑝𝑡𝑠}=`$  $`\left[\begin{array}{c}{\displaystyle \frac{2\mathrm{sin}(\chi )^23(\frac{}{t}\mathrm{a}(t))^2\mathrm{sin}(\chi )^21+\mathrm{cos}(\chi )^2}{\mathrm{a}(t)^2\mathrm{sin}(\chi )^2}},\mathrm{\hspace{0.17em}0},\mathrm{\hspace{0.17em}0},\mathrm{\hspace{0.17em}0}\\ 0,{\displaystyle \frac{2\mathrm{a}(t)\mathrm{sin}(\chi )^2(\frac{^2}{t^2}\mathrm{a}(t))(\frac{}{t}\mathrm{a}(t))^2\mathrm{sin}(\chi )^21+\mathrm{cos}(\chi )^2}{\mathrm{sin}(\chi )^2}},\mathrm{\hspace{0.17em}0},\mathrm{\hspace{0.17em}0}\\ 0,\mathrm{\hspace{0.17em}0},\mathrm{\hspace{0.17em}2}\mathrm{a}(t)\mathrm{sin}(\chi )^2(\frac{^2}{t^2}\mathrm{a}(t))+\mathrm{sin}(\chi )^2+(\frac{}{t}\mathrm{a}(t))^2\mathrm{sin}(\chi )^2,\mathrm{\hspace{0.17em}0}\\ 0,\mathrm{\hspace{0.17em}0},\mathrm{\hspace{0.17em}0},\mathrm{\hspace{0.17em}2}\mathrm{a}(t)\mathrm{sin}(\chi )^2\mathrm{sin}(\theta )^2(\frac{^2}{t^2}\mathrm{a}(t))+\mathrm{sin}(\theta )^2\mathrm{sin}(\chi )^2+(\frac{}{t}\mathrm{a}(t))^2\mathrm{sin}(\chi )^2\mathrm{sin}(\theta )^2\end{array}\right]`$  $`])`$

Let’s the matter is defined by the energy-momentum tensor $`T_{\mu ,\nu }`$ = (p+ $`\rho `$ ) $`u_\mu `$ $`u_\nu `$  +p $`g_{\mu ,\nu }`$ :

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

$`T:=\mathrm{table}([\mathrm{𝑖𝑛𝑑𝑒𝑥}\mathrm{\_}\mathrm{𝑐ℎ𝑎𝑟}=[1,1],\mathrm{𝑐𝑜𝑚𝑝𝑡𝑠}=`$ $`\left[\begin{array}{cccc}\rho (t)& 0& 0& 0\\ 0& \mathrm{a}(t)^2\mathrm{p}(t)& 0& 0\\ 0& 0& \mathrm{a}(t)^2\mathrm{sin}(\chi )^2\mathrm{p}(t)& 0\\ 0& 0& 0& \mathrm{a}(t)^2\mathrm{sin}(\chi )^2\mathrm{sin}(\theta )^2\mathrm{p}(t)\end{array}\right]`$  $`])`$

Then the first Einstein equation $`G_{0,\mathrm{\hspace{0.17em}0}}`$ \- $`\mathrm{\Lambda }`$ $`g_{0,\mathrm{\hspace{0.17em}0}}`$ = - 8 $`\pi `$ $`T_{0,\mathrm{\hspace{0.17em}0}}`$ (we add the so-called cosmological constant $`\mathrm{\Lambda }`$ , which can be considered as the energy density of vacuum):

get\_compts(G):
%:
e1 := simplify(%):
get\_compts(g):
%:
e2 := simplify(%):
get\_compts(T):
%:
e3 := simplify(%):
E := expand(e1/(-3)) = -8\*Pi\*expand(e3/(-3)) +\
expand(e2\*Lambda/(-3));#first Einstein equation

$$E_1:=\frac{1}{\mathrm{a}(t)^2}+\frac{(\frac{}{t}\mathrm{a}(t))^2}{\mathrm{a}(t)^2}=\frac{8}{3}\pi \rho (t)+\frac{1}{3}\mathrm{\Lambda }$$

Second equation is:

get\_compts(G):
%:
e1 := simplify(%):
get\_compts(g):
%:
e2 := simplify(%):
get\_compts(T):
%:
e3 := simplify(%):
E :=\
expand(e1/a(t)^2) = -8\*Pi\*expand(e3/a(t)^2) +\
expand(e2\*Lambda/a(t)^2);#second Einstein equation

$$E_2:=2\frac{\frac{^2}{t^2}\mathrm{a}(t)}{\mathrm{a}(t)}+\frac{(\frac{}{t}\mathrm{a}(t))^2}{\mathrm{a}(t)^2}+\frac{1}{\mathrm{a}(t)^2}=8\pi \mathrm{p}(t)+\mathrm{\Lambda }$$

Two other equations are tautological.

After elementary transformation we have for the second equation:

expand(simplify(E-E)/2);

$$\frac{\frac{^2}{t^2}\mathrm{a}(t)}{\mathrm{a}(t)}=4\pi \mathrm{p}(t)+\frac{1}{3}\mathrm{\Lambda }\frac{4}{3}\pi \rho (t)$$

Similarly, for the hyperbolical geometry we have:

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

$`g:=\mathrm{table}([\mathrm{𝑖𝑛𝑑𝑒𝑥}\mathrm{\_}\mathrm{𝑐ℎ𝑎𝑟}=[1,1],\mathrm{𝑐𝑜𝑚𝑝𝑡𝑠}=`$ $`\left[\begin{array}{cccc}\hfill 1& 0& 0& 0\\ \hfill 0& \mathrm{a}(t)^2& 0& 0\\ \hfill 0& 0& \mathrm{a}(t)^2\mathrm{sinh}(\chi )^2& 0\\ \hfill 0& 0& 0& \mathrm{a}(t)^2\mathrm{sinh}(\chi )^2\mathrm{sin}(\theta )^2\end{array}\right]`$  $`])`$

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 );

$`G:=\mathrm{table}([\mathrm{𝑖𝑛𝑑𝑒𝑥}\mathrm{\_}\mathrm{𝑐ℎ𝑎𝑟}=[1,1],\mathrm{𝑐𝑜𝑚𝑝𝑡𝑠}=`$ $`\left[\begin{array}{c}{\displaystyle \frac{2\mathrm{sinh}(\chi )^2+3(\frac{}{t}\mathrm{a}(t))^2\mathrm{sinh}(\chi )^2+1\mathrm{cosh}(\chi )^2}{\mathrm{a}(t)^2\mathrm{sinh}(\chi )^2}},\mathrm{\hspace{0.17em}0},\mathrm{\hspace{0.17em}0},\mathrm{\hspace{0.17em}0}\\ 0,{\displaystyle \frac{2\mathrm{a}(t)\mathrm{sinh}(\chi )^2(\frac{^2}{t^2}\mathrm{a}(t))+(\frac{}{t}\mathrm{a}(t))^2\mathrm{sinh}(\chi )^2+1\mathrm{cosh}(\chi )^2}{\mathrm{sinh}(\chi )^2}},\mathrm{\hspace{0.17em}0},\mathrm{\hspace{0.17em}0}\\ 0,\mathrm{\hspace{0.17em}0},\mathrm{\hspace{0.17em}2}\mathrm{a}(t)\mathrm{sinh}(\chi )^2(\frac{^2}{t^2}\mathrm{a}(t))\mathrm{sinh}(\chi )^2+(\frac{}{t}\mathrm{a}(t))^2\mathrm{sinh}(\chi )^2,\mathrm{\hspace{0.17em}0}\\ 0,\mathrm{\hspace{0.17em}0},\mathrm{\hspace{0.17em}0},\mathrm{\hspace{0.17em}2}\mathrm{a}(t)\mathrm{sinh}(\chi )^2\mathrm{sin}(\theta )^2(\frac{^2}{t^2}\mathrm{a}(t))\mathrm{sin}(\theta )^2\mathrm{sinh}(\chi )^2+(\frac{}{t}\mathrm{a}(t))^2\mathrm{sinh}(\chi )^2\mathrm{sin}(\theta )^2\end{array}\right]`$  $`])`$

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

$`T:=\mathrm{table}([\mathrm{𝑖𝑛𝑑𝑒𝑥}\mathrm{\_}\mathrm{𝑐ℎ𝑎𝑟}=[1,1],\mathrm{𝑐𝑜𝑚𝑝𝑡𝑠}=`$ $`\left[\begin{array}{cccc}\rho (t)& 0& 0& 0\\ 0& \mathrm{a}(t)^2\mathrm{p}(t)& 0& 0\\ 0& 0& \mathrm{a}(t)^2\mathrm{sinh}(\chi )^2\mathrm{p}(t)& 0\\ 0& 0& 0& \mathrm{a}(t)^2\mathrm{sinh}(\chi )^2\mathrm{sin}(\theta )^2\mathrm{p}(t)\end{array}\right]`$  $`])`$

get\_compts(G):
%:
e1 := simplify(%):
get\_compts(g):
%:
e2 := simplify(%):
get\_compts(T):
%:
e3 := simplify(%):
E :=\
expand(e1/(-3)) = -8\*Pi\*expand(e3/(-3)) +\
expand(e2\*Lambda/(-3));#first Einstein equation

$$E_1:=\frac{1}{\mathrm{a}(t)^2}+\frac{(\frac{}{t}\mathrm{a}(t))^2}{\mathrm{a}(t)^2}=\frac{8}{3}\pi \rho (t)+\frac{1}{3}\mathrm{\Lambda }$$

get\_compts(G):
%:
e1 := simplify(%):
get\_compts(g):
%:
e2 := simplify(%):
get\_compts(T):
%:
e3 := simplify(%):
E :=\
expand(e1/a(t)^2) = -8\*Pi\*expand(e3/a(t)^2) +\
expand(e2\*Lambda/a(t)^2):#second Einstein equation
expand(simplify(E-E)/2);

$$\frac{\frac{^2}{t^2}\mathrm{a}(t)}{\mathrm{a}(t)}=4\pi \mathrm{p}(t)+\frac{1}{3}\mathrm{\Lambda }\frac{4}{3}\pi \rho (t)$$

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

$$\mathrm{𝑏𝑎𝑠𝑖𝑐}\mathrm{\_}\mathit{1}:=\frac{(\frac{}{t}\mathrm{a}(t))^2}{\mathrm{a}(t)^2}=\frac{K}{\mathrm{a}(t)^2}+\frac{1}{3}\mathrm{\Lambda }+\frac{8}{3}\pi \rho (t)$$

$$\mathrm{𝑏𝑎𝑠𝑖𝑐}\mathrm{\_}\mathit{2}:=\frac{\frac{^2}{t^2}\mathrm{a}(t)}{\mathrm{a}(t)}=4\pi (\mathrm{p}(t)+\frac{1}{3}\rho (t))+\frac{1}{3}\mathrm{\Lambda }$$

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 recent value of this parameter as $`H_0`$ , then $`H_0`$ = 65 $`\frac{\mathrm{𝑘𝑚}}{s\mathrm{𝑀𝑝𝑐}}`$ (ps is the parallax second corresponding to distance about of 3.26 light-years or 3\* $`10^{13}`$ km) .

When K= $`\mathrm{\Lambda }`$ =0, the universe is defined by the so-called critical density $`\rho _c`$ = $`\frac{3H_{0}^{}{}_{}{}^{2}}{8\pi }`$ =7.9\* $`10^{(30)}`$ $`\frac{g}{\mathrm{𝑐𝑚}^3}`$ . Let’s denote the recent relative values $`\mathrm{\Omega }_i`$ = $`\frac{\rho _i}{\rho _c}`$ of the matter, radiation, cosmological constant and curvature as:

$`\mathrm{\Omega }_M`$ = $`\frac{8\pi \rho _M}{3H_{0}^{}{}_{}{}^{2}}`$ , $`\mathrm{\Omega }_R`$ = $`\frac{8\pi \rho _R}{3H_{0}^{}{}_{}{}^{2}}`$ , $`\mathrm{\Omega }_\mathrm{\Lambda }`$ = $`\frac{\mathrm{\Lambda }}{3H_{0}^{}{}_{}{}^{2}}`$ , $`\mathrm{\Omega }_K`$ = \- $`\frac{K}{a_{0}^{}{}_{}{}^{2}H_{0}^{}{}_{}{}^{2}}`$ .

Here $`a_0`$ is the recent scaling factor. At this moment the contribution of the radiation to density is negligible $`\mathrm{\Omega }_R`$ $`<`$$`<`$ $`\mathrm{\Omega }_M`$ . To rewrite the evolutionary equation in terms of $`\mathrm{\Omega }_i`$ we have to take into account the time dependence of these parameters: $`\mathrm{\Lambda }`$ does not depend on a(t) (in framework of standard model), the dependence of curvature ~ $`\frac{1}{a^2}`$ , density of matter ~ $`\frac{1}{a^3}`$ , but for radiation ~ $`\frac{1}{a^4}`$ . The last results from the change of density of quanta ~ $`\frac{1}{a^3}`$ and their energy ~ $`\frac{1}{a}`$ . The resulting equation can be written as:

basic\_3 := H(a)^2 =\
H^2\*(Omega\[R\]\*a^4/a^4 + Omega\[M\]\*a^3/a^3\
\+ Omega\[Lambda\] + Omega\[K\]\*a^2/a^2);

$$\mathrm{𝑏𝑎𝑠𝑖𝑐}\mathrm{\_}\mathit{3}:=\mathrm{H}(a)^2=H_{0}^{}{}_{}{}^{2}(\frac{\mathrm{\Omega }_Ra_{0}^{}{}_{}{}^{4}}{a^4}+\frac{\mathrm{\Omega }_Ma_{0}^{}{}_{}{}^{3}}{a^3}+\mathrm{\Omega }_\mathrm{\Lambda }+\frac{\mathrm{\Omega }_Ka_{0}^{}{}_{}{}^{2}}{a^2})$$

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

1 = $`\mathrm{\Omega }_M+\mathrm{\Omega }_R`$ \+ $`\mathrm{\Omega }_\mathrm{\Lambda }`$ \+ $`\mathrm{\Omega }_K`$ .

Let us transit to the new variables: y= $`\frac{a}{a_0}`$ , $`\tau `$ = $`H_0`$ ( $`tt_0`$ ), where $`t_0`$ 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\*d(t) --\> a\*H\*(d(y)/d(tau))/a =\
(d(a)/d(t))/a = H

$$\mathrm{𝑏𝑎𝑠𝑖𝑐}\mathrm{\_}\mathit{4}:=\frac{}{\tau }\mathrm{y}(\tau )=\sqrt{1+(\frac{1}{\mathrm{y}(\tau )}1)\mathrm{\Omega }_M+(\mathrm{y}(\tau )^21)\mathrm{\Omega }_\mathrm{\Lambda }}$$

Now we will consider some typical standard models. For review see, for example, , or on-line survey .

#### 4.3.1 Einstein static

rhs(basic\_1) = 0:
rhs(basic\_2) = 0:
allvalues(\
solve({%, %%},{a(t), Lambda}) );# static universe

$`\{\mathrm{\Lambda }=12\pi \mathrm{p}(t)+4\pi \rho (t),\mathrm{a}(t)=\sqrt{{\displaystyle \frac{K}{4\pi \mathrm{p}(t)+4\pi \rho (t)}}}\},`$ $`\{\mathrm{\Lambda }=12\pi \mathrm{p}(t)+4\pi \rho (t),\mathrm{a}(t)=\sqrt{{\displaystyle \frac{K}{4\pi \mathrm{p}(t)+4\pi \rho (t)}}}\}`$

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

$`a^2`$ = $`\frac{1}{4\pi \rho },\mathrm{\Lambda }=4\pi \rho `$

This solution has a very simple form and requires K = +1 (the closed spherical universe) and positive $`\mathrm{\Lambda }`$ (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).

#### 4.3.2 Einstein-de Sitter

$`\mathrm{\Omega }_M`$ = 1 and, as consequence $`\mathrm{\Omega }_\mathrm{\Lambda }`$ = $`\mathrm{\Omega }_K`$ = 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(%);

$$\sqrt{\mathrm{y}(\tau )}(\frac{}{\tau }\mathrm{y}(\tau ))=1$$

$`\mathrm{y}(\tau )=`$ $`{\displaystyle \frac{1}{4}}(12\tau +8)^{(2/3)},\mathrm{y}(\tau )=({\displaystyle \frac{1}{4}}(12\tau +8)^{(1/3)}+{\displaystyle \frac{1}{4}}I\sqrt{3}(12\tau +8)^{(1/3)})^2,`$  $`\mathrm{y}(\tau )=`$  $`({\displaystyle \frac{1}{4}}(12\tau +8)^{(1/3)}{\displaystyle \frac{1}{4}}I\sqrt{3}(12\tau +8)^{(1/3)})^2`$

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 $`\mathrm{y}(\tau )=\frac{(12\tau +8)^{(\frac{2}{3})}}{4}`$ it is possible to find the universe’s age:

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

$$\frac{2}{3}\frac{1}{H_0}$$

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

$$\mathrm{.9756859072\hspace{0.17em}10}^{10}$$

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

#### 4.3.3 de Sitter and anti-de Sitter

The empty universes ( $`\mathrm{\Omega }_M`$ =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));

$$(\frac{}{t}\mathrm{a}(t))^2=1+\mathrm{a}(t)^2\lambda \mathrm{~}$$

$`\mathrm{a}(t)={\displaystyle \frac{1}{\sqrt{\lambda \mathrm{~}}}},\mathrm{a}(t)={\displaystyle \frac{1}{\sqrt{\lambda \mathrm{~}}}},`$ $`\mathrm{a}(t)={\displaystyle \frac{1}{2}}{\displaystyle \frac{(1+{\displaystyle \frac{\%2^2}{\%1^2}})\%1}{\%2\sqrt{\lambda \mathrm{~}}}},\mathrm{a}(t)={\displaystyle \frac{1}{2}}{\displaystyle \frac{(1+{\displaystyle \frac{\%1^2}{\%2^2}})\%2}{\%1\sqrt{\lambda \mathrm{~}}}}`$  $`\%1:=e^{(\mathrm{\_}\mathit{C1}\sqrt{\lambda \mathrm{~}})}`$  $`\%2:=e^{(t\sqrt{\lambda \mathrm{~}})}`$

It is obvious, that the first two solutions don’t satisfy the basic\_2. Another solutions give:

a(t) = $`\frac{\mathrm{cosh}(\sqrt{\lambda }(tC))}{2\sqrt{\lambda }}`$

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));

$$(\frac{}{t}\mathrm{a}(t))^2=1+\mathrm{a}(t)^2\lambda \mathrm{~}$$

$`\mathrm{a}(t)={\displaystyle \frac{\sqrt{\lambda \mathrm{~}}}{\lambda \mathrm{~}}},\mathrm{a}(t)={\displaystyle \frac{\sqrt{\lambda \mathrm{~}}}{\lambda \mathrm{~}}},`$ $`\mathrm{a}(t)={\displaystyle \frac{\mathrm{sin}(t\sqrt{\lambda \mathrm{~}}\mathrm{\_}\mathit{C1}\sqrt{\lambda \mathrm{~}})}{\sqrt{\lambda \mathrm{~}}}},`$  $`\mathrm{a}(t)={\displaystyle \frac{\mathrm{sin}(t\sqrt{\lambda \mathrm{~}}\mathrm{\_}\mathit{C1}\sqrt{\lambda \mathrm{~}})}{\sqrt{\lambda \mathrm{~}}}}`$

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.

#### 4.3.4 Closed Friedmann-Lemaitre

$`\mathrm{\Omega }_M`$ $`>`$ 1 $`\mathrm{\Lambda }`$ = 0 and, as consequence $`\mathrm{\Omega }_K`$ $`<`$ 0 (K=+1). This is the closed spherical universe.

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

$$\mathrm{y}(\tau )((\frac{}{\tau }\mathrm{y}(\tau ))^2+b)=\mathrm{\Omega }_M$$

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

$$\frac{\sqrt{\mathrm{y}(\tau )(\mathrm{y}(\tau )b+\mathrm{\Omega }_M)}}{\mathrm{y}(\tau )},\frac{\sqrt{\mathrm{y}(\tau )(\mathrm{y}(\tau )b+\mathrm{\Omega }_M)}}{\mathrm{y}(\tau )}$$

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

$$\mathrm{d}(t)=\sqrt{\frac{y}{b({\displaystyle \frac{\mathrm{\Omega }_M}{b}}y)}}\mathrm{d}(y)$$

Introducing a new variable $`\varphi `$ we have :

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

$$\sqrt{\frac{y}{{\displaystyle \frac{\mathrm{\Omega }_M}{b}}y}}=\mathrm{tan}(\varphi )$$

$$y=\frac{\mathrm{tan}(\varphi )^2\mathrm{\Omega }_M}{b(1+\mathrm{tan}(\varphi )^2)}$$

$$y=\frac{\mathrm{sin}(\varphi )^2\mathrm{\Omega }_M}{\mathrm{cos}(\varphi )^2b(1+{\displaystyle \frac{\mathrm{sin}(\varphi )^2}{\mathrm{cos}(\varphi )^2}})}$$

But this is y = $`\frac{\mathrm{\Omega }_M\mathrm{sin}(\varphi )^2}{b}`$  = $`\frac{\mathrm{\Omega }_M(1\mathrm{cos}(2\varphi ))}{2b}`$ . And

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

$$\mathrm{d}(y)=2\mathrm{sin}(\varphi )\mathrm{cos}(\varphi )(\mathrm{d}(\varphi )\&\mathrm{^}\frac{\mathrm{\Omega }_M}{b})+\mathrm{sin}(\varphi )^2\mathrm{d}(\frac{\mathrm{\Omega }_M}{b})$$

Then our equation results in

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

$$\sqrt{\frac{1+\mathrm{cos}(\varphi )^2}{b\mathrm{cos}(\varphi )^2}}$$

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

$$\mathrm{d}(t)=2\frac{\mathrm{\Omega }_M(1+\mathrm{cos}(\varphi )^2)\mathrm{d}(\varphi )}{b^{(3/2)}}$$

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

$$\mathrm{𝑠𝑜𝑙}\mathrm{\_}\mathit{1}:=t=\frac{1}{2}\frac{\mathrm{\Omega }_M(2\varphi +\mathrm{sin}(2\varphi ))}{(\mathrm{\Omega }_M1)^{(3/2)}}$$

This equation in the combination with y( $`\varphi `$ ) 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,rhs(sol\_1));#age \
of universe vs Omega\[M\]
plot(%, Omega\[M\]=1.01..10, axes=BOXED, \
title=‘age of universe \[1/H \]‘);

$$\mathrm{𝑠𝑜𝑙}\mathrm{\_}\mathit{2}:=\mathrm{arcsin}(\frac{\sqrt{\mathrm{\Omega }_M(\mathrm{\Omega }_M1)}}{\mathrm{\Omega }_M}),\mathrm{arcsin}(\frac{\sqrt{\mathrm{\Omega }_M(\mathrm{\Omega }_M1)}}{\mathrm{\Omega }_M})$$

$$\frac{1}{2}\frac{\mathrm{\Omega }_M(2\mathrm{arcsin}({\displaystyle \frac{\sqrt{\mathrm{\Omega }_M(\mathrm{\Omega }_M1)}}{\mathrm{\Omega }_M}})+\mathrm{sin}(2\mathrm{arcsin}({\displaystyle \frac{\sqrt{\mathrm{\Omega }_M(\mathrm{\Omega }_M1)}}{\mathrm{\Omega }_M}})))}{(\mathrm{\Omega }_M1)^{(3/2)}}$$

So, the increase of $`\mathrm{\Omega }_M`$ decreases the age of universe, which is less than in Einstein-de Sitter model. The last and the modern observations of microwave background suggesting $`\mathrm{\Omega }_K`$ ~0 do not agree with this model.

#### 4.3.5 Open Friedmann-Lemaitre

$`\mathrm{\Omega }_M`$ $`<`$ 1 $`\mathrm{\Lambda }`$ = 0 and, as consequence $`\mathrm{\Omega }_K`$ $`>`$ 0 (K= -1). This is the 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

$$\mathrm{y}(\tau )((\frac{}{\tau }\mathrm{y}(\tau ))^2a)=\mathrm{\Omega }_M$$

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

$$\frac{\sqrt{\mathrm{y}(\tau )(\mathrm{y}(\tau )a+\mathrm{\Omega }_M)}}{\mathrm{y}(\tau )},\frac{\sqrt{\mathrm{y}(\tau )(\mathrm{y}(\tau )a+\mathrm{\Omega }_M)}}{\mathrm{y}(\tau )}$$

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

$$\frac{}{\tau }\mathrm{y}(\tau )=\sqrt{\frac{\mathrm{y}(\tau )a+\mathrm{\Omega }_M}{\mathrm{y}(\tau )}}$$

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

$`\mathrm{𝑠𝑜𝑙}\mathrm{\_}\mathit{1}:={\displaystyle \frac{1}{2}}(2\tau \%1\mathrm{\Omega }_M+\mathrm{\Omega }_M\mathrm{ln}({\displaystyle \frac{1}{2}}`$ $`{\displaystyle \frac{\mathrm{\Omega }_M2\mathrm{y}(\tau )+2\mathrm{y}(\tau )\mathrm{\Omega }_M2\sqrt{\mathrm{y}(\tau )(\mathrm{y}(\tau )+\mathrm{y}(\tau )\mathrm{\Omega }_M\mathrm{\Omega }_M)}\%1}{\sqrt{1\mathrm{\Omega }_M}}})`$  $`2\mathrm{\_}\mathit{C1}\%1\mathrm{\Omega }_M+2\tau \%12\sqrt{\mathrm{y}(\tau )(\mathrm{y}(\tau )+\mathrm{y}(\tau )\mathrm{\Omega }_M\mathrm{\Omega }_M)}\%1+`$  $`2\mathrm{\_}\mathit{C1}\%1)/(\sqrt{1\mathrm{\Omega }_M}(\mathrm{\Omega }_M1))=0`$  $`\%1:=\sqrt{1\mathrm{\Omega }_M}`$

#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

$`\mathrm{𝑠𝑜𝑙}\mathrm{\_}\mathit{2}:={\displaystyle \frac{1}{4}}(4\%1\mathrm{\Omega }_M+\mathrm{\Omega }_M\mathrm{ln}({\displaystyle \frac{(\mathrm{\Omega }_M22\%1)^2}{\mathrm{\Omega }_M1}})`$ $`\mathrm{ln}\left({\displaystyle \frac{(\mathrm{\Omega }_M+2y2y\mathrm{\Omega }_M+2\sqrt{y(y+y\mathrm{\Omega }_M\mathrm{\Omega }_M)}\%1)^2}{\mathrm{\Omega }_M1}}\right)\mathrm{\Omega }_M`$  $`+4\sqrt{y(y+y\mathrm{\Omega }_M\mathrm{\Omega }_M)}\%1)/(\sqrt{1\mathrm{\Omega }_M}(\mathrm{\Omega }_M1))`$  $`\%1:=\sqrt{1\mathrm{\Omega }_M}`$

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 $`\mathrm{\Omega }_M`$ $`>`$0 (almost empty world)

limit(sol\_2,Omega\[M\]=0);

$$\sqrt{y^2}$$

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\]

$$\mathrm{.1463528861\hspace{0.17em}10}^{11}$$

The dependence on $`\mathrm{\Omega }_M`$ looks as:

-(subs(y=0,sol\_2)-subs(y=1,sol\_2)):#age\
of universe \[ 1/H \]
plot(%, Omega\[M\]=0..1,axes=boxed,\
title=‘age of universe \[ 1/H\]‘);

It is natural, the growth of matter density decreases the universe age. If estimations of $`\mathrm{\Omega }_M`$ 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).

#### 4.3.6 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 ($`\mathrm{\Omega }_M`$ =1) but the cosmological term slightly dominates $`\mathrm{\Omega }_\mathrm{\Lambda }`$ $`>`$ $`2\mathrm{\Omega }_M`$ . There is the delay of expansion, which looks like transition to collapse, but the domination of $`\mathrm{\Omega }_\mathrm{\Lambda }`$ causes the further expansion. This model can explain the large age of universe for large $`H_0`$ and to explain the misalignment in the ”red shift-distance” distribution for distant quasars.

The next example demonstrate the collapsing behavior of an 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‘);

#### 4.3.7 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‘);

#### 4.3.8 Our universe (?)

The modern observations of microwave background suggest : $`\mathrm{\Omega }_K`$ $``$ 0 (flat), $`\mathrm{\Omega }_\mathrm{\Lambda }`$ $``$ 0.7, $`\mathrm{\Omega }_M`$ $``$ 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”).

### 4.4 Beginning

#### 4.4.1 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 ). 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:
g\_compts := a(t)^2:
g\_compts := b(t)^2:
g\_compts := c(t)^2:
g := create(\[-1,-1\], eval(g\_compts));
ginv := invert( g, ’detg’ ):

$`g:=\mathrm{table}([\mathrm{𝑖𝑛𝑑𝑒𝑥}\mathrm{\_}\mathrm{𝑐ℎ𝑎𝑟}=[1,1],\mathrm{𝑐𝑜𝑚𝑝𝑡𝑠}=`$ $`\left[\begin{array}{cccc}\hfill 1& 0& 0& 0\\ \hfill 0& \mathrm{a}(t)^2& 0& 0\\ \hfill 0& 0& \mathrm{b}(t)^2& 0\\ \hfill 0& 0& 0& \mathrm{c}(t)^2\end{array}\right]])`$

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

$$\mathrm{𝑇ℎ𝑒}\mathrm{𝐸𝑖𝑛𝑠𝑡𝑒𝑖𝑛}\mathrm{𝑇𝑒𝑛𝑠𝑜𝑟}$$

$$\mathrm{𝑛𝑜𝑛}\mathrm{𝑧𝑒𝑟𝑜}\mathrm{𝑐𝑜𝑚𝑝𝑜𝑛𝑒𝑛𝑡𝑠}:$$

$`\mathit{G11}=`$ $`{\displaystyle \frac{(\frac{}{t}\mathrm{b}(t))(\frac{}{t}\mathrm{a}(t))\mathrm{c}(t)+(\frac{}{t}\mathrm{c}(t))(\frac{}{t}\mathrm{a}(t))\mathrm{b}(t)+(\frac{}{t}\mathrm{c}(t))(\frac{}{t}\mathrm{b}(t))\mathrm{a}(t)}{\mathrm{a}(t)\mathrm{b}(t)\mathrm{c}(t)}}`$

$`\mathit{G22}=`$ $`{\displaystyle \frac{\mathrm{a}(t)^2((\frac{^2}{t^2}\mathrm{b}(t))\mathrm{c}(t)+(\frac{^2}{t^2}\mathrm{c}(t))\mathrm{b}(t)+(\frac{}{t}\mathrm{c}(t))(\frac{}{t}\mathrm{b}(t)))}{\mathrm{b}(t)\mathrm{c}(t)}}`$

$`\mathit{G33}=`$ $`{\displaystyle \frac{\mathrm{b}(t)^2((\frac{^2}{t^2}\mathrm{a}(t))\mathrm{c}(t)+(\frac{^2}{t^2}\mathrm{c}(t))\mathrm{a}(t)+(\frac{}{t}\mathrm{c}(t))(\frac{}{t}\mathrm{a}(t)))}{\mathrm{a}(t)\mathrm{c}(t)}}`$

$`\mathit{G44}=`$ $`{\displaystyle \frac{\mathrm{c}(t)^2((\frac{^2}{t^2}\mathrm{a}(t))\mathrm{b}(t)+(\frac{^2}{t^2}\mathrm{b}(t))\mathrm{a}(t)+(\frac{}{t}\mathrm{b}(t))(\frac{}{t}\mathrm{a}(t)))}{\mathrm{a}(t)\mathrm{b}(t)}}`$

$$\mathrm{𝑐ℎ𝑎𝑟𝑎𝑐𝑡𝑒𝑟}:[\mathit{1},\mathit{1}]$$

In the vacuum state the right-hand sides of Einstein equations are 0. Then

E\_eqn := get\_compts(Estn):
e1 := numer( E\_eqn ) = 0;
e2 :=\
expand( numer( E\_eqn )/a(t)^2 ) = 0;
e3 :=\
expand( numer( E\_eqn )/b(t)^2 ) = 0;
e4 :=\
expand( numer( E\_eqn )/c(t)^2 ) = 0;

$`\mathit{e1}:=(\frac{}{t}\mathrm{b}(t))(\frac{}{t}\mathrm{a}(t))\mathrm{c}(t)(\frac{}{t}\mathrm{c}(t))(\frac{}{t}\mathrm{a}(t))\mathrm{b}(t)`$ $`(\frac{}{t}\mathrm{c}(t))(\frac{}{t}\mathrm{b}(t))\mathrm{a}(t)=0`$

$$\mathit{e2}:=(\frac{^2}{t^2}\mathrm{b}(t))\mathrm{c}(t)+(\frac{^2}{t^2}\mathrm{c}(t))\mathrm{b}(t)+(\frac{}{t}\mathrm{c}(t))(\frac{}{t}\mathrm{b}(t))=0$$

$$\mathit{e3}:=(\frac{^2}{t^2}\mathrm{a}(t))\mathrm{c}(t)+(\frac{^2}{t^2}\mathrm{c}(t))\mathrm{a}(t)+(\frac{}{t}\mathrm{c}(t))(\frac{}{t}\mathrm{a}(t))=0$$

$$\mathit{e4}:=(\frac{^2}{t^2}\mathrm{a}(t))\mathrm{b}(t)+(\frac{^2}{t^2}\mathrm{b}(t))\mathrm{a}(t)+(\frac{}{t}\mathrm{b}(t))(\frac{}{t}\mathrm{a}(t))=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) );

$`\mathit{e5}:=b\mathrm{\_}\mathit{0}t^{(p\mathrm{\_}\mathit{2}2+p\mathrm{\_}\mathit{1}+p\mathrm{\_}\mathit{3})}p\mathrm{\_}\mathit{2}a\mathrm{\_}\mathit{0}p\mathrm{\_}\mathit{1}c\mathrm{\_}\mathit{0}`$ $`c\mathrm{\_}\mathit{0}t^{(p\mathrm{\_}\mathit{2}2+p\mathrm{\_}\mathit{1}+p\mathrm{\_}\mathit{3})}p\mathrm{\_}\mathit{3}a\mathrm{\_}\mathit{0}p\mathrm{\_}\mathit{1}b\mathrm{\_}\mathit{0}`$  $`c\mathrm{\_}\mathit{0}t^{(p\mathrm{\_}\mathit{2}2+p\mathrm{\_}\mathit{1}+p\mathrm{\_}\mathit{3})}p\mathrm{\_}\mathit{3}b\mathrm{\_}\mathit{0}p\mathrm{\_}\mathit{2}a\mathrm{\_}\mathit{0}=0`$

$`\mathit{e6}:=b\mathrm{\_}\mathit{0}t^{(p\mathrm{\_}\mathit{3}2+p\mathrm{\_}\mathit{2})}p\mathrm{\_}\mathit{2}^2c\mathrm{\_}\mathit{0}b\mathrm{\_}\mathit{0}t^{(p\mathrm{\_}\mathit{3}2+p\mathrm{\_}\mathit{2})}p\mathrm{\_}\mathit{2}c\mathrm{\_}\mathit{0}`$ $`+c\mathrm{\_}\mathit{0}t^{(p\mathrm{\_}\mathit{3}2+p\mathrm{\_}\mathit{2})}p\mathrm{\_}\mathit{3}^2b\mathrm{\_}\mathit{0}c\mathrm{\_}\mathit{0}t^{(p\mathrm{\_}\mathit{3}2+p\mathrm{\_}\mathit{2})}p\mathrm{\_}\mathit{3}b\mathrm{\_}\mathit{0}`$  $`+c\mathrm{\_}\mathit{0}t^{(p\mathrm{\_}\mathit{3}2+p\mathrm{\_}\mathit{2})}p\mathrm{\_}\mathit{3}b\mathrm{\_}\mathit{0}p\mathrm{\_}\mathit{2}=0`$

$`\mathit{e7}:=a\mathrm{\_}\mathit{0}t^{(p\mathrm{\_}\mathit{3}2+p\mathrm{\_}\mathit{1})}p\mathrm{\_}\mathit{1}^2c\mathrm{\_}\mathit{0}a\mathrm{\_}\mathit{0}t^{(p\mathrm{\_}\mathit{3}2+p\mathrm{\_}\mathit{1})}p\mathrm{\_}\mathit{1}c\mathrm{\_}\mathit{0}`$ $`+c\mathrm{\_}\mathit{0}t^{(p\mathrm{\_}\mathit{3}2+p\mathrm{\_}\mathit{1})}p\mathrm{\_}\mathit{3}^2a\mathrm{\_}\mathit{0}c\mathrm{\_}\mathit{0}t^{(p\mathrm{\_}\mathit{3}2+p\mathrm{\_}\mathit{1})}p\mathrm{\_}\mathit{3}a\mathrm{\_}\mathit{0}`$  $`+c\mathrm{\_}\mathit{0}t^{(p\mathrm{\_}\mathit{3}2+p\mathrm{\_}\mathit{1})}p\mathrm{\_}\mathit{3}a\mathrm{\_}\mathit{0}p\mathrm{\_}\mathit{1}=0`$

$`\mathit{e8}:=a\mathrm{\_}\mathit{0}t^{(p\mathrm{\_}\mathit{1}2+p\mathrm{\_}\mathit{2})}p\mathrm{\_}\mathit{1}^2b\mathrm{\_}\mathit{0}a\mathrm{\_}\mathit{0}t^{(p\mathrm{\_}\mathit{1}2+p\mathrm{\_}\mathit{2})}p\mathrm{\_}\mathit{1}b\mathrm{\_}\mathit{0}`$ $`+b\mathrm{\_}\mathit{0}t^{(p\mathrm{\_}\mathit{1}2+p\mathrm{\_}\mathit{2})}p\mathrm{\_}\mathit{2}^2a\mathrm{\_}\mathit{0}b\mathrm{\_}\mathit{0}t^{(p\mathrm{\_}\mathit{1}2+p\mathrm{\_}\mathit{2})}p\mathrm{\_}\mathit{2}a\mathrm{\_}\mathit{0}`$  $`+b\mathrm{\_}\mathit{0}t^{(p\mathrm{\_}\mathit{1}2+p\mathrm{\_}\mathit{2})}p\mathrm{\_}\mathit{2}a\mathrm{\_}\mathit{0}p\mathrm{\_}\mathit{1}=0`$

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) );

$$\mathit{e9}:=b\mathrm{\_}\mathit{0}a\mathrm{\_}\mathit{0}c\mathrm{\_}\mathit{0}(p\mathrm{\_}\mathit{2}p\mathrm{\_}\mathit{1}+p\mathrm{\_}\mathit{3}p\mathrm{\_}\mathit{1}+p\mathrm{\_}\mathit{3}p\mathrm{\_}\mathit{2})=0$$

$$\mathit{e10}:=b\mathrm{\_}\mathit{0}c\mathrm{\_}\mathit{0}(p\mathrm{\_}\mathit{2}^2p\mathrm{\_}\mathit{2}+p\mathrm{\_}\mathit{3}^2p\mathrm{\_}\mathit{3}+p\mathrm{\_}\mathit{3}p\mathrm{\_}\mathit{2})=0$$

$$\mathit{e11}:=a\mathrm{\_}\mathit{0}c\mathrm{\_}\mathit{0}(p\mathrm{\_}\mathit{1}^2p\mathrm{\_}\mathit{1}+p\mathrm{\_}\mathit{3}^2p\mathrm{\_}\mathit{3}+p\mathrm{\_}\mathit{3}p\mathrm{\_}\mathit{1})=0$$

$$\mathit{e12}:=a\mathrm{\_}\mathit{0}b\mathrm{\_}\mathit{0}(p\mathrm{\_}\mathit{1}^2p\mathrm{\_}\mathit{1}+p\mathrm{\_}\mathit{2}^2p\mathrm{\_}\mathit{2}+p\mathrm{\_}\mathit{2}p\mathrm{\_}\mathit{1})=0$$

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);

$$\mathit{e13}:=p\mathrm{\_}\mathit{2}p\mathrm{\_}\mathit{1}+p\mathrm{\_}\mathit{3}p\mathrm{\_}\mathit{1}+p\mathrm{\_}\mathit{3}p\mathrm{\_}\mathit{2}=0$$

$$\mathit{e14}:=p\mathrm{\_}\mathit{2}^2p\mathrm{\_}\mathit{2}+p\mathrm{\_}\mathit{3}^2p\mathrm{\_}\mathit{3}+p\mathrm{\_}\mathit{3}p\mathrm{\_}\mathit{2}=0$$

$$\mathit{e15}:=p\mathrm{\_}\mathit{1}^2p\mathrm{\_}\mathit{1}+p\mathrm{\_}\mathit{3}^2p\mathrm{\_}\mathit{3}+p\mathrm{\_}\mathit{3}p\mathrm{\_}\mathit{1}=0$$

$$\mathit{e16}:=p\mathrm{\_}\mathit{1}^2p\mathrm{\_}\mathit{1}+p\mathrm{\_}\mathit{2}^2p\mathrm{\_}\mathit{2}+p\mathrm{\_}\mathit{2}p\mathrm{\_}\mathit{1}=0$$

The following manipulation gives a connection between parameters:

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

$$p\mathrm{\_}\mathit{2}^2p\mathrm{\_}\mathit{2}+p\mathrm{\_}\mathit{3}^2p\mathrm{\_}\mathit{3}+p\mathrm{\_}\mathit{1}^2p\mathrm{\_}\mathit{1}=0$$

That is $`p\mathrm{\_}\mathit{1}^2`$ \+ $`p\mathrm{\_}\mathit{2}^2`$ \+ $`p\mathrm{\_}\mathit{3}^2`$ = p\_1 \+ p\_2 \+ p\_3 (3).

collect(e15-e16,{p\_2,p\_3});
collect(e14-e15,{p\_1,p\_2});

$$p\mathrm{\_}\mathit{3}^2+(p\mathrm{\_}\mathit{1}1)p\mathrm{\_}\mathit{3}p\mathrm{\_}\mathit{2}^2+(1p\mathrm{\_}\mathit{1})p\mathrm{\_}\mathit{2}=0$$

$$p\mathrm{\_}\mathit{2}^2+(p\mathrm{\_}\mathit{3}1)p\mathrm{\_}\mathit{2}p\mathrm{\_}\mathit{1}^2+(1p\mathrm{\_}\mathit{3})p\mathrm{\_}\mathit{1}=0$$

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);

$$p\mathrm{\_}\mathit{3}+p\mathrm{\_}\mathit{1}3+p\mathrm{\_}\mathit{2}$$

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;

$$p\mathrm{\_}\mathit{3}+p\mathrm{\_}\mathit{1}1+p\mathrm{\_}\mathit{2}$$

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

$`\mathrm{𝑠𝑜𝑙}\mathrm{\_}\mathit{1}:=\{p\mathrm{\_}\mathit{1}={\displaystyle \frac{1}{2}}p+{\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{2}}\%1,p\mathrm{\_}\mathit{2}={\displaystyle \frac{1}{2}}p+{\displaystyle \frac{1}{2}}+{\displaystyle \frac{1}{2}}\%1\},`$ $`\{p\mathrm{\_}\mathit{1}={\displaystyle \frac{1}{2}}p+{\displaystyle \frac{1}{2}}+{\displaystyle \frac{1}{2}}\%1,p\mathrm{\_}\mathit{2}={\displaystyle \frac{1}{2}}p+{\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{2}}\%1\}`$  $`\%1:=\sqrt{3p^2+2p+1}`$

$`\mathrm{𝑠𝑜𝑙}\mathrm{\_}\mathit{2}:=`$ $`\{p\mathrm{\_}\mathit{2}={\displaystyle \frac{1}{2}}p+{\displaystyle \frac{3}{2}}+{\displaystyle \frac{1}{2}}I(p1)\sqrt{3},p\mathrm{\_}\mathit{1}={\displaystyle \frac{1}{2}}p+{\displaystyle \frac{3}{2}}{\displaystyle \frac{1}{2}}I(p1)\sqrt{3}\},`$  $`\{p\mathrm{\_}\mathit{2}={\displaystyle \frac{1}{2}}p+{\displaystyle \frac{3}{2}}{\displaystyle \frac{1}{2}}I(p1)\sqrt{3},p\mathrm{\_}\mathit{1}={\displaystyle \frac{1}{2}}p+{\displaystyle \frac{3}{2}}+{\displaystyle \frac{1}{2}}I(p1)\sqrt{3}\}`$

So, if the parameters p\_1, p\_2, p\_3 are real, we have $`p\mathrm{\_}\mathit{1}^2`$ \+ $`p\mathrm{\_}\mathit{2}^2`$ \+ $`p\mathrm{\_}\mathit{3}^2`$ = 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 - %% - %;

$$\mathrm{𝑠𝑜𝑙}\mathrm{\_}p\mathrm{\_}\mathit{1}:=\frac{1}{2}p+\frac{1}{2}\frac{1}{2}\sqrt{(3p+1)(p1)}$$

$$\mathrm{𝑠𝑜𝑙}\mathrm{\_}p\mathrm{\_}\mathit{2}:=\frac{1}{2}p+\frac{1}{2}+\frac{1}{2}\sqrt{(3p+1)(p1)}$$

$$\mathrm{𝑠𝑜𝑙}\mathrm{\_}p\mathrm{\_}\mathit{3}:=p$$

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–$`>`$ $`\mathrm{}`$ (so-called ”pancake” singularity) or one direction of expansion and two direction of contraction for t–$`>`$0 (”sigar” singularity).

Now let’s 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) :

$`g_{0,\mathrm{\hspace{0.17em}0}}`$ = -1, $`g_{\alpha ,\beta }`$ = $`\eta _{a,b}`$ $`(e^a)_\alpha `$ $`(e^b)_\beta `$ ,

where $`\eta _{a,b}`$ is the time dependent matrix, $`(e^a)_\alpha `$ is the frame vector (a changes from 1 to 3, that is the number of space triad). We will use the following representation for $`\eta `$ :

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

$$\eta :=\left[\begin{array}{ccc}\mathrm{a}(t)^2& 0& 0\\ 0& \mathrm{b}(t)^2& 0\\ 0& 0& \mathrm{c}(t)^2\end{array}\right]$$

$$\kappa 2:=\left[\begin{array}{ccc}2\frac{\frac{}{t}\mathrm{a}(t)}{\mathrm{a}(t)}& 0& 0\\ 0& 2\frac{\frac{}{t}\mathrm{b}(t)}{\mathrm{b}(t)}& 0\\ 0& 0& 2\frac{\frac{}{t}\mathrm{c}(t)}{\mathrm{c}(t)}\end{array}\right]$$

The first vacuum Einstein equation has a form

$`R_{0}^{}{}_{}{}^{0}`$ = \- $`\frac{1}{2}`$ $`\frac{}{t}\kappa _{a}^{}{}_{}{}^{a}`$  - $`\frac{1}{4}`$ $`\kappa _{a}^{}{}_{}{}^{b}`$ $`\kappa _{b}^{}{}_{}{}^{a}`$ = 0

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

$$\mathrm{𝑒𝑞}\mathrm{\_}\mathit{1}:=\frac{\frac{^2}{t^2}\mathrm{a}(t)}{\mathrm{a}(t)}\frac{\frac{^2}{t^2}\mathrm{b}(t)}{\mathrm{b}(t)}\frac{\frac{^2}{t^2}\mathrm{c}(t)}{\mathrm{c}(t)}=0$$

The next Einstein equations demand the definition of Bianchi structured constants: $`C^{\mathrm{𝑎𝑏}}`$ = $`n^{\mathrm{𝑎𝑏}}+e^{\mathrm{𝑎𝑏𝑐}}a_c`$ , $`(C^c)_{\mathrm{𝑎𝑏}}`$ = $`e_{\mathrm{𝑎𝑏𝑑}}`$ $`C^{\mathrm{𝑑𝑐}}`$ . Here $`e_{\mathrm{𝑎𝑏𝑐}}`$ = $`e^{\mathrm{𝑎𝑏𝑐}}`$ is the unit antisymmetric symbol, $`n^{\mathrm{𝑎𝑏}}`$ is the symmetric ”tensor”, which can be expressed through its principal values $`n_i`$ , $`a_c`$ is the vector. The Bianchi types are:

We will analyze type IX model, where $`C^{11}`$ = $`C^{22}`$ = $`C^{33}`$ = $`C_{23}^{}{}_{}{}^{1}`$ = $`C_{31}^{}{}_{}{}^{2}`$ = $`C_{12}^{}{}_{}{}^{3}`$ = 1. The curvature $`P_{a}^{}{}_{}{}^{b}`$ = $`\frac{1}{2\eta }`$ {2 $`C^{\mathrm{𝑏𝑑}}`$ $`C_{\mathrm{𝑎𝑑}}`$ \+ $`C^{\mathrm{𝑑𝑏}}`$ $`C_{\mathrm{𝑎𝑑}}`$ \+ $`C^{\mathrm{𝑏𝑑}}`$ $`C_{\mathrm{𝑑𝑎}}`$  \- $`C_{d}^{}{}_{}{}^{d}`$ ( $`(C^b)_a+C_{a}^{}{}_{}{}^{b}`$ ) + $`\delta _{a}^{}{}_{}{}^{b}`$ \[ $`(C^d)_{d}^{}{}_{}{}^{2}`$ -2 $`C^{\mathrm{𝑑𝑓}}`$ $`C_{\mathrm{𝑑𝑓}}`$ \]} is

Cab :=\
array(\[,,\]);#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(\[,,\]):
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,%);

$$\mathrm{𝐶𝑎𝑏}:=\left[\begin{array}{ccc}\hfill 1& \hfill 0& \hfill 0\\ \hfill 0& \hfill 1& \hfill 0\\ \hfill 0& \hfill 0& \hfill 1\end{array}\right]$$

$$C\mathrm{\_}\mathrm{𝑎𝑏}:=\left[\begin{array}{ccc}\mathrm{a}(t)^4& 0& 0\\ 0& \mathrm{b}(t)^4& 0\\ 0& 0& \mathrm{c}(t)^4\end{array}\right]$$

$$\mathrm{𝐶𝑎}\mathrm{\_}b:=\left[\begin{array}{ccc}\mathrm{a}(t)^2& 0& 0\\ 0& \mathrm{b}(t)^2& 0\\ 0& 0& \mathrm{c}(t)^2\end{array}\right]$$

$$C\mathrm{\_}a\mathrm{\_}b:=\left[\begin{array}{ccc}\mathrm{a}(t)^2& 0& 0\\ 0& \mathrm{b}(t)^2& 0\\ 0& 0& \mathrm{c}(t)^2\end{array}\right]$$

$$P:=\left[\begin{array}{c}\frac{1}{2}\frac{\mathrm{a}(t)^4+\mathrm{b}(t)^42\mathrm{b}(t)^2\mathrm{c}(t)^2+\mathrm{c}(t)^4}{\mathrm{a}(t)^2\mathrm{b}(t)^2\mathrm{c}(t)^2},\mathrm{\hspace{0.17em}0},\mathrm{\hspace{0.17em}0}\\ 0,\frac{1}{2}\frac{\mathrm{b}(t)^4+\mathrm{a}(t)^42\mathrm{a}(t)^2\mathrm{c}(t)^2+\mathrm{c}(t)^4}{\mathrm{a}(t)^2\mathrm{b}(t)^2\mathrm{c}(t)^2},\mathrm{\hspace{0.17em}0}\\ 0,\mathrm{\hspace{0.17em}0},\frac{1}{2}\frac{\mathrm{c}(t)^4\mathrm{a}(t)^4+2\mathrm{a}(t)^2\mathrm{b}(t)^2\mathrm{b}(t)^4}{\mathrm{a}(t)^2\mathrm{b}(t)^2\mathrm{c}(t)^2}\end{array}\right]$$

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:

$`R_{a}^{}{}_{}{}^{b}`$ = \- $`\frac{1}{2\sqrt{\eta }}`$  \[( $`\frac{}{t}\sqrt{\eta }\kappa _{a}^{}{}_{}{}^{b}`$ ) - $`P_{a}^{}{}_{}{}^{b}`$ \] = 0

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

$`[{\displaystyle \frac{(\frac{^2}{t^2}\mathrm{a}(t))\mathrm{b}(t)\mathrm{c}(t)+(\frac{}{t}\mathrm{b}(t))(\frac{}{t}\mathrm{a}(t))\mathrm{c}(t)+(\frac{}{t}\mathrm{c}(t))(\frac{}{t}\mathrm{a}(t))\mathrm{b}(t)}{\mathrm{a}(t)\mathrm{b}(t)\mathrm{c}(t)}}`$ $`+{\displaystyle \frac{{\displaystyle \frac{1}{2}}(\mathrm{a}(t)^4+\mathrm{b}(t)^42\mathrm{b}(t)^2\mathrm{c}(t)^2+\mathrm{c}(t)^4)}{\mathrm{a}(t)^2\mathrm{b}(t)^2\mathrm{c}(t)^2}},\mathrm{\hspace{0.17em}0},\mathrm{\hspace{0.17em}0}]`$  $`[0,{\displaystyle \frac{(\frac{^2}{t^2}\mathrm{b}(t))\mathrm{a}(t)\mathrm{c}(t)+(\frac{}{t}\mathrm{b}(t))(\frac{}{t}\mathrm{a}(t))\mathrm{c}(t)+(\frac{}{t}\mathrm{c}(t))(\frac{}{t}\mathrm{b}(t))\mathrm{a}(t)}{\mathrm{a}(t)\mathrm{b}(t)\mathrm{c}(t)}}`$  $`+{\displaystyle \frac{{\displaystyle \frac{1}{2}}(\mathrm{b}(t)^4+\mathrm{a}(t)^42\mathrm{a}(t)^2\mathrm{c}(t)^2+\mathrm{c}(t)^4)}{\mathrm{a}(t)^2\mathrm{b}(t)^2\mathrm{c}(t)^2}},\mathrm{\hspace{0.17em}0}]`$  $`[0,\mathrm{\hspace{0.17em}0},`$  $`{\displaystyle \frac{(\frac{^2}{t^2}\mathrm{c}(t))\mathrm{a}(t)\mathrm{b}(t)+(\frac{}{t}\mathrm{c}(t))(\frac{}{t}\mathrm{a}(t))\mathrm{b}(t)+(\frac{}{t}\mathrm{c}(t))(\frac{}{t}\mathrm{b}(t))\mathrm{a}(t)}{\mathrm{a}(t)\mathrm{b}(t)\mathrm{c}(t)}}`$  $`{\displaystyle \frac{1}{2}}{\displaystyle \frac{\mathrm{c}(t)^4\mathrm{a}(t)^4+2\mathrm{a}(t)^2\mathrm{b}(t)^2\mathrm{b}(t)^4}{\mathrm{a}(t)^2\mathrm{b}(t)^2\mathrm{c}(t)^2}}]`$

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);

$$\mathrm{𝑒𝑞}\mathrm{\_}\mathit{2}:=\frac{\frac{}{t}(\frac{}{t}\mathrm{a}(t))\mathrm{b}(t)\mathrm{c}(t)}{\mathrm{a}(t)\mathrm{b}(t)\mathrm{c}(t)}=\frac{1}{2}\frac{(\mathrm{b}(t)^2\mathrm{c}(t)^2)^2\mathrm{a}(t)^4}{\mathrm{a}(t)^2\mathrm{b}(t)^2\mathrm{c}(t)^2}$$

$$\mathrm{𝑒𝑞}\mathrm{\_}\mathit{3}:=\frac{\frac{}{t}(\frac{}{t}\mathrm{b}(t))\mathrm{a}(t)\mathrm{c}(t)}{\mathrm{a}(t)\mathrm{b}(t)\mathrm{c}(t)}=\frac{1}{2}\frac{(\mathrm{a}(t)^2\mathrm{c}(t)^2)^2\mathrm{b}(t)^4}{\mathrm{a}(t)^2\mathrm{b}(t)^2\mathrm{c}(t)^2}$$

$$\mathrm{𝑒𝑞}\mathrm{\_}\mathit{4}:=\frac{\frac{}{t}(\frac{}{t}\mathrm{c}(t))\mathrm{b}(t)\mathrm{a}(t)}{\mathrm{a}(t)\mathrm{b}(t)\mathrm{c}(t)}=\frac{1}{2}\frac{(\mathrm{a}(t)^2\mathrm{b}(t)^2)^2\mathrm{c}(t)^4}{\mathrm{a}(t)^2\mathrm{b}(t)^2\mathrm{c}(t)^2}$$

The last group of Einstein equations results from

$`R_{a}^{}{}_{}{}^{0}`$ = \- $`\frac{1}{2}`$ $`\kappa _{b}^{}{}_{}{}^{c}`$ ( $`(C^b)_{\mathrm{𝑐𝑎}}`$ - $`\delta _{a}^{}{}_{}{}^{b}`$ $`(C^d)_{\mathrm{𝑑𝑐}}`$ ) = 0

It is obvious, that in our case this equation is identically zero (see expression for $`(C^b)_{\mathrm{𝑐𝑎}}`$ ). The substitutions of a= $`e^\alpha `$ , b= $`e^\beta `$ , c= $`e^\gamma `$ (do not miss these exponents with frame vectors!), dt=a b c (d $`\tau `$ ) and the expressions $`\frac{}{t}a`$ = $`\frac{}{\tau }a`$ $`\frac{}{t}\tau `$ = $`\frac{\frac{}{\tau }a}{abc}`$ = $`\frac{\frac{}{\tau }\alpha }{bc}`$ and $`\frac{^2}{t^2}a`$ = $`\frac{1}{abc}`$ $`\frac{}{\tau }\frac{\frac{}{\tau }\alpha }{bc}`$ = $`\frac{1}{abc}`$ \[ $`\frac{\frac{^2}{\tau ^2}\alpha }{bc}`$ - ( $`\frac{(\frac{}{\tau }\gamma )+(\frac{}{\tau }\beta )}{bc}`$ ) $`\frac{}{\tau }\alpha `$ \] 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;

$`\mathrm{𝑒𝑞}\mathrm{\_}\mathit{1}:={\displaystyle \frac{1}{2}}(\frac{^2}{\tau ^2}\alpha (\tau ))+{\displaystyle \frac{1}{2}}(\frac{^2}{\tau ^2}\beta (\tau ))+{\displaystyle \frac{1}{2}}(\frac{^2}{\tau ^2}\gamma (\tau ))=`$ $`(\frac{}{\tau }\alpha (\tau ))(\frac{}{\tau }\beta (\tau ))+(\frac{}{\tau }\alpha (\tau ))(\frac{}{\tau }\gamma (\tau ))+(\frac{}{\tau }\beta (\tau ))(\frac{}{\tau }\gamma (\tau ))`$

$$\mathrm{𝑒𝑞}\mathrm{\_}\mathit{2}:=2(\frac{^2}{\tau ^2}\alpha (\tau ))=(\mathrm{b}(t)^2\mathrm{c}(t)^2)^2\mathrm{a}(t)^4$$

$$\mathrm{𝑒𝑞}\mathrm{\_}\mathit{3}:=2(\frac{^2}{\tau ^2}\beta (\tau ))=(\mathrm{a}(t)^2\mathrm{c}(t)^2)^2\mathrm{b}(t)^4$$

$$\mathrm{𝑒𝑞}\mathrm{\_}\mathit{4}:=2(\frac{^2}{\tau ^2}\gamma (\tau ))=(\mathrm{a}(t)^2\mathrm{b}(t)^2)^2\mathrm{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 $`\tau `$ = $`\mathrm{ln}(t)`$ \+ const (dt=a b c d $`\tau `$ = $`t^{(p_1+p_2+p_3)}`$ d $`\tau `$ , as result d $`\tau `$  = $`\frac{\mathrm{𝑑𝑡}}{t}`$ ) and $`\frac{}{\tau }\alpha `$ = $`p_1`$ , $`\frac{}{\tau }\beta `$ = $`p_2`$ , $`\frac{}{\tau }\gamma `$ = $`p_3`$ . We assume that the right-hand sides are close to zero, but contribution of $`e^{(4\alpha )}`$ \- 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;

$$\mathrm{𝑒𝑞}\mathrm{\_}\mathit{5}:=\frac{^2}{\tau ^2}\alpha (\tau )=\frac{1}{2}e^{(4\alpha (\tau ))}$$

$$\mathrm{𝑒𝑞}\mathrm{\_}\mathit{6}:=\frac{^2}{\tau ^2}\beta (\tau )=\frac{1}{2}e^{(4\alpha (\tau ))}$$

$$\mathrm{𝑒𝑞}\mathrm{\_}\mathit{7}:=\frac{^2}{\tau ^2}\gamma (\tau )=\frac{1}{2}e^{(4\alpha (\tau ))}$$

The first equation is analog of equation describing the one-dimensional motion ( $`\alpha `$ 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 $`p_1`$ = -0.5 (green curve) is reflected on barrier (red curve) with change of velocity sign. But eq\_5, eq\_6, eq\_7 results in $`p_1+p_2`$ = const, $`p_1+p_3`$ = const therefore

$`(p^{\mathrm{𝑛𝑒𝑤}})_2`$ = $`p_2`$ \+ $`p_1`$  - $`(p^{\mathrm{𝑛𝑒𝑤}})_1`$ = $`p_2`$ \+ 2 $`p_1`$ , $`(p^{\mathrm{𝑛𝑒𝑤}})_3`$ = $`p_3`$ \+ $`p_1`$  - $`(p^{\mathrm{𝑛𝑒𝑤}})_1`$ = $`p_3`$ \+ 2 $`p_1`$ ,

t = a b c = $`e^{((1+2p_1)\tau )}`$  ,

and the result is:

a = $`t^{(\frac{p_1}{1+2p_1})}`$ , b = $`t^{(\frac{p_2+2p_1}{1+2p_1})}`$ , c = $`t^{(\frac{p_3+2p_1}{1+2p_1})}`$

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 $`p\mathrm{\_}\mathit{1}^2`$ \+ $`p\mathrm{\_}\mathit{2}^2`$ \+ $`p\mathrm{\_}\mathit{3}^2`$ =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 the 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 ).

#### 4.4.2 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 $`\mathrm{𝑑𝑠}^2`$ = 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);

$$r=_0^{\mathit{t0}}\frac{1}{\mathrm{a}(t)}𝑑t$$

$$r=\frac{{\displaystyle _{\mathit{t0}\mathit{H0}}^0}{\displaystyle \frac{1}{\mathrm{y}(\tau )}}𝑑\tau }{\mathit{a0}\mathit{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 );

$$2\frac{1}{\mathit{H0}}$$

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)

$$2\frac{\varphi }{\sqrt{b}\mathit{a0}\mathit{H0}}$$

$$2\varphi $$

The last example is of interest. The increase of $`\varphi `$ from 0 to $`\pi `$ /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 $`\varphi `$ $`<`$ $`\pi `$ /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, , on-line reviews ).

For escape the horizon problem we have to suppose the accelerated expansion with $`\frac{^2}{t^2}\mathrm{a}(t)`$ $`>`$ 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:

$`\frac{\frac{^2}{t^2}\mathrm{a}(t)}{\mathrm{a}(t)}=4\pi (\mathrm{p}(t)+\frac{\rho (t)}{3})+\frac{\mathrm{\Lambda }}{3}`$

This equation demonstrates that in the usual conditions (p $`>`$ 0) the pressure is a source of gravitation, but if p $`<`$ \- $`\frac{\rho }{3}`$  ( $`\mathrm{\Lambda }`$ = 0) or $`\mathrm{\Lambda }`$ $`>`$ 4 $`\pi `$ $`\rho `$ (p=0) the repulsion dominates and the expansion is accelerated. As the source of the repulsion, the scalar ”inflaton” field $`\varphi `$ is considered.

Let us introduce the simple Lagrangian for this homogeneous field:

L=- $`\frac{1}{2}`$ $`\frac{}{x_i}\varphi `$ $`\frac{}{x_j}\varphi `$ $`g^{\mathrm{𝑖𝑗}}`$ - V( $`\varphi `$ ),

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

$`\frac{1}{\sqrt{g}}`$ $`\frac{}{x_i}\sqrt{g}g^{\mathrm{𝑖𝑗}}(\frac{}{x_j}\varphi )`$ \+ $`\frac{}{\varphi }V`$ = 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:
g\_compts := a(t)^2:
g\_compts:= a(t)^2\*r^2:
g\_compts:= 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,):
Tt := prod(g\_det\_sq,%):
Ttt := cov\_diff( Tt, coord, Cf2 ):
get\_compts(%)/det\_scal:
expand(-%):
field\_eq := % + diff(V(phi),phi) = 0;

$$\mathrm{𝑓𝑖𝑒𝑙𝑑}\mathrm{\_}\mathrm{𝑒𝑞}:=3\frac{(\frac{}{t}\varphi (t))(\frac{}{t}\mathrm{a}(t))}{\mathrm{a}(t)}+(\frac{^2}{t^2}\varphi (t))+(\frac{}{\varphi }\mathrm{V}(\varphi ))=0$$

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

$`T_{\mathrm{𝑖𝑗}}`$ = $`g_{\mathrm{𝑖𝑗}}`$ $`L`$ \- $`g_{\mathrm{𝑗𝑙}}`$ $`f_{x\mathrm{\_}i}`$ $`\frac{}{f_{x\mathrm{\_}l}}L`$ ,

where $`f_{x\mathrm{\_}i}`$ = $`\frac{}{x_i}\varphi `$ . As consequence, $`\rho `$ = \[- $`\frac{1}{2}`$ $`(\frac{}{t}\varphi )^2`$ \+ V( $`\varphi `$ )\] + $`(\frac{}{t}\varphi )^2`$ = $`\frac{1}{2}`$ $`(\frac{}{t}\varphi )^2`$ \+ V( $`\varphi `$ ), and (for i, j = 1,2,3) $`T_{\mathrm{𝑖𝑗}}`$ = $`g_{\mathrm{𝑖𝑗}}`$ ( $`\frac{1}{2}`$ $`(\frac{}{t}\varphi )^2`$ \- V( $`\varphi `$ )) $``$ p = $`\frac{1}{2}`$ $`(\frac{}{t}\varphi )^2`$ \- V( $`\varphi `$ ). Last expression emplies that the necessary condition p $`<`$ \- $`\frac{\rho }{3}`$  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
basic\_inf\_2 := field\_eq;

$$\mathrm{𝑏𝑎𝑠𝑖𝑐}\mathrm{\_}\mathrm{𝑖𝑛𝑓}\mathrm{\_}\mathit{1}:=\frac{K}{\mathrm{a}(t)^2}+\frac{(\frac{}{t}\mathrm{a}(t))^2}{\mathrm{a}(t)^2}=\frac{8}{3}\pi (\frac{1}{2}(\frac{}{t}\varphi (t))^2+\mathrm{V}(\varphi ))$$

$$\mathrm{𝑏𝑎𝑠𝑖𝑐}\mathrm{\_}\mathrm{𝑖𝑛𝑓}\mathrm{\_}\mathit{2}:=3\frac{(\frac{}{t}\varphi (t))(\frac{}{t}\mathrm{a}(t))}{\mathrm{a}(t)}+(\frac{^2}{t^2}\varphi (t))+(\frac{}{\varphi }\mathrm{V}(\varphi ))=0$$

The meaning of K is considered earlier, here we will suppose K=0 (local flatness). In the slow-rolling approximation $`\frac{^2}{t^2}\varphi (t)`$ , $`(\frac{}{t}\varphi (t))^2`$ $`>`$0 and for potential V( $`\varphi `$ ) = $`\frac{m\varphi ^2}{2}`$ 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
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;

$$\mathrm{𝑏𝑎𝑠𝑖𝑐}\mathrm{\_}\mathrm{𝑖𝑛𝑓}\mathrm{\_}\mathit{1}\mathrm{\_}n:=\frac{(\frac{}{t}\mathrm{a}(t))^2}{\mathrm{a}(t)^2}=\frac{4}{3}\pi m^2\varphi (t)^2$$

$$\mathrm{𝑏𝑎𝑠𝑖𝑐}\mathrm{\_}\mathrm{𝑖𝑛𝑓}\mathrm{\_}\mathit{2}\mathrm{\_}n:=3\frac{(\frac{}{t}\varphi (t))(\frac{}{t}\mathrm{a}(t))}{\mathrm{a}(t)}+m^2\varphi (t)=0$$

dsolve({basic\_inf\_1\_n,basic\_inf\_2\_n,\
phi(0)=phi0,a(0)=a0},{phi(t),a(t)}):
expand(%);

$$\{\mathrm{a}(t)=\frac{e^{(1/6m^2t^2)}\mathit{a0}}{e^{(2/3\varphi 0\sqrt{\pi }m\sqrt{3}t)}},\varphi (t)=\frac{1}{6}\frac{\sqrt{3}mt}{\sqrt{\pi }}+\varphi 0\}$$

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

a(t)= $`a_0`$ $`e^{(Ht)}`$ , where H = 2 $`\sqrt{\frac{\pi }{3}}`$ $`\varphi _0`$ m

When t ~ $`\frac{2\sqrt{3\pi }}{m}`$ $`\varphi _0`$ 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\*m,t=2\*sqrt(3\*Pi)/m\*phi},\
a\*exp(H\*t));

$$\mathrm{a}(t)=a_0e^{(4\pi \varphi _{0}^{}{}_{}{}^{2})}$$

If at the beginning of inflation the energy and size of universe are defined by Plank density and length, then $`m^2`$ $`(\varphi ^2)_0`$ ~1, $`a_0`$ $`\mathrm{~}\mathit{10}^{(33)}`$ cm. The necessary value of m $`\mathrm{~}\mathit{10}^{(6)}`$  is constrained from the value of fluctuation of background radiation ( $`\mathrm{~}\mathit{10}^{(5)}`$ ). As result a ~ $`10^{(33)}`$ $`e^{(\frac{4\pi }{m^2})}`$ = $`10^{(33)}`$ $`e^{(4\pi \mathrm{\hspace{0.17em}10}^{12})}`$ 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);

$$\mathrm{𝑟𝑒ℎ}\mathrm{\_}\mathrm{𝑖𝑛𝑓}:=3(\frac{}{t}\varphi (t))H+\mathrm{\Gamma }(\frac{}{t}\varphi (t))+(\frac{^2}{t^2}\varphi (t))=m^2\varphi (t)$$

Here H is the Hubble constant, $`\mathrm{\Gamma }`$ is the decay rate, which is defined by coupling with material boson field. When H $`>`$ $`\mathrm{\Gamma }`$ , there are the coherent oscillations:

dsolve(subs(Gamma=0,reh\_inf),phi(t));
diff(%,t);

$`\varphi (t)=\mathrm{\_}\mathit{C1}e^{(1/2(3H+\sqrt{(3H2m)(3H+2m)})t)}+`$ $`\mathrm{\_}\mathit{C2}e^{(1/2(3H\sqrt{(3H2m)(3H+2m)})t)}`$

$`\frac{}{t}\varphi (t)=\mathrm{\_}\mathit{C1}({\displaystyle \frac{3}{2}}H{\displaystyle \frac{1}{2}}\%1)e^{(1/2(3H+\%1)t)}+`$ $`\mathrm{\_}\mathit{C2}({\displaystyle \frac{3}{2}}H+{\displaystyle \frac{1}{2}}\%1)e^{(1/2(3H\%1)t)}`$  $`\%1:=\sqrt{(3H2m)(3H+2m)}`$

The last expression demonstrates the decrease of energy density via time increase ( $`\rho `$ ~ $`(\frac{}{t}\varphi )^2`$ ) and, as result, the decrease of H (see basic\_inf\_1). This leads to H $`<`$ $`\mathrm{\Gamma }`$ , 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.

## 5 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.