Cosmology

Section 1: Setup the metric

(* Defining coordinate system *)
coord={t,r,θ,ϕ};
(* Set these variable for coordinate, hopefully mathematica will be smarter when we Simplify *)
(* Noted that e is an expansions parameter *)
$Assumptions=t∈Reals&&r∈Reals&&θ∈Reals&&ϕ∈Reals&&e∈Reals;
(* Defining the metric*)
metric = DiagonalMatrix[{1,-(((R[t])^2)/(1-k (r^2)+D/r)),-((R[t])^2)r^2,-((R[t])^2)(r Sin[θ])^2}];
metricsign=-1;

metric //MatrixForm

(* Load the package from the folder which you put diffgeo.m *)
SetDirectory["C:\Users\15-P240TX\Documents"]
<<diffgeo.m

Section 2 : The equation of motion, typical Einstein equation + cosmological constant

(* axion can be straightforwardly added *)

?Einstein

display[Riemann]

RicciScalar

display[Einstein]

display[RicciTensor]

?lower

u={u0,u1,u2,u3};
ul=lower[u]//FullSimplify
(*{ud0,ud1,ud2,ud3}=ul;

usub={u0->1,u1->0,u2->0,u3->0};(*Comoving observer*)

Tmunu=(ρ[t]+P[t])TensorProduct[ul,ul]-(P[t])*metric/.usub;
display[Tmunu]

Einstein[[1,1]]==8 π Tmunu[[1,1]]
(Einstein[[2,2]]==8 π Tmunu[[2,2]]/.D->0)//Simplify//Expand

((-1+k (r^2))/((R[t])^2))((k/(-1+k (r^2)))+((8 π P[t] ((R[t])^2))/(-1+k (r^2)))+(((R'[t])^2)/(-1+k (r^2)))+((2 R[t] R''[t])/(-1+k (r^2))))==0//Simplify//Expand

Tdu1=raise[Tmunu,{2}];
display[Tdu1]

EOM=RicciTensor-8 π (Tmunu-metric (Tdu1[[1,1]]+Tdu1[[2,2]]+Tdu1[[3,3]]+Tdu1[[4,4]])/2)//Simplify

display[EOM]

(EOM[[3,3]]/.{R''[t]->-((4 π)/3)R[t] (ρ[t]+3P[t])}) (3/(2 (r^2)))//Simplify

(*Only way [[1,1]] consistent with [[3,3]] is D=0*)

(D/(4 (r^3)))+(8/3) π ((R[t])^2) ρ[t]==(k+((R'[t])^2))

(*ALSO Friedmann Eqn for D=0*)

(8/3) π ((R[t])^2) ρ[t]==(k+((R'[t])^2))

?covariant

?contract

contract[covariant[Tmunu],{1,2}]
contract[covariant[Tmunu],{1,3}]

contract[Tmunu]

CON=contract[covariant[Tmunu],{1,2}];

CON[[1]]

((8/3) π ((R[t])^2) ρ[t]==(k+((R'[t])^2))/.{P[t]->w ρ[t]}/.{ρ[t]->ρ0 ((((R0/(R[t]))))^(3 (1+w)))})

DSolve[(8/3) π ρ0 (((R0))^(3 (1+w))) ((R[t])^(-1-3 w))==((R'[t])^2),R[t],t]//Simplify

DSolve[(8/3) π ((a[t])^(-1-3 w)) (((a0))^(3 (1+w))) ρ0==((a'[t])^2),a[t],t]