The Robertson-Walker metric

In the following the Einstein field equations will be solved for the Robertson-Walker standard model cosmology  metric. See Weinberg, page 471. Note: Weinberg conventions for Riemann differ in sign from definition used here.

IndepVar = {t, r, θ,  φ}

The function   metric[]calculates the coefficients of the metric tensor from a given line element.   metric[] takes the line element  ds^2  and a list of coordinates as input variables. The result is the metric tensor of the underlying space.

metric[lineelement_, independentvars_List] := 
        Block[{lenindependent, differentials, diffmatrix, metricform, 
            varmetric, gh, sum, equation, rule, varhelp, zeros, zerorule}, 
          lenindependent = Length[independentvars]; 
          differentials = Dt /@ independentvars; 
          diffmatrix = Outer[Times, differentials, differentials]; 
          metricform = Array[gh, {lenindependent, lenindependent}]; 
          varmetric = Variables[metricform]; 
          If[Length[metricform] == Length[diffmatrix], sum = 0; 
            Do[Do[sum = 
                  sum + metricform[[i, 
                        j]] diffmatrix[[
                        i, j]], {j, 1, 
                  lenindependent}], {i, 1, lenindependent}], sum = 0]; 
          If[sum === 0, Return[sum], sum = sum - lineelement; 
            equation = CoefficientList[sum, differentials] == 0; 
            rule = Solve[equation, varmetric]; 
            metricform = metricform /. rule; 
            varmetric = Variables[metricform]; varhelp = {}; 
            Do[If[! (FreeQ[
                      varmetric[[i]], 
                      gh]), 
                AppendTo[varhelp, 
                  varmetric[[i]]]], {i, 
                1, Length[varmetric]}]; zeros = Table[0, {Length[varhelp]}]; 
            SubstRule[x_, y_] := x -> y; 
            zerorule = Thread[SubstRule[varhelp, zeros]]; 
            metricform = Flatten[metricform /. zerorule, 1]; 
            metricform = 
              Expand[(1/2) (metricform + Transpose[metricform])]]; 
          metricform];

The Robertson-Walker  ds^2  in the coordinates  t ,r, th  and  phi   is

ds =  - 
          Dt[t]^2 +  ((R[t] )^2) ((1/(1 - k (r^2))) 
                Dt[r]^2 + 
              (r^2) (Dt[θ]^2 + 
                    Sin[θ]^2 Dt[φ]^2))

The corresponding metric is

g = metric[ds, {t, r, θ, φ}];

MatrixForm[g]

ing  =  Inverse[g];

The energy-momentum tensor for a perfect fluid is given by

en  =  {{p[t] + ρ[t], 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 
              0}, {0, 0, 0, 0}}  +  p[t] g ;

MatrixForm[%]

EnMomT[a_, b_, g_]  :=  Block[{}, en[[a, b]] ]

Table[EnMomT[a, b, g],  {a, 1, 4}, {b, 1, 4}]

The Einstein quations are then R_mu_nu - (1/2) g_mu_nu R = - 8 π T_mu_nu .

Up to now we have only discussed the line element of the metric and its related metric tensor. To derive the field equations for the vacuum case in GR, we have to introduce other tensors. One of the essential quantities determining the field equations are the Christoffel symbols. These symbols are related to the metric tensor in a straightforward way.

Christoffel[m_, a_, b_, g_, ing_] := 
      Block[{n}, 
        Expand[(1/2) (∑_(n = 1)^(Length[g]))ing[[
                  m, n]] ((∂_(IndepVar[[b]]))g[[a, 
                        n]] + (∂_(IndepVar[[a]]))g[[b, 
                        n]] - (∂_(IndepVar[[n]]))g[[a, 
                        b]])]]

Riemann[a_, b_, c_, d_, g_, ing_] := 
      Block[{}, 
        Expand[(∂_(IndepVar[[c]]))Christoffel[a, b, d, g, 
                ing] - (∂_(IndepVar[[d]]))Christoffel[a, b, c, g, 
                ing] + (∑_(e = 1)^(Length[g]))Christoffel[e, b, d, 
                  g, ing] Christoffel[a, e, c, g, 
                  ing] - (∑_(e = 1)^(Length[g]))Christoffel[e, b, c, 
                  g, ing] Christoffel[a, e, d, g, ing]]]

Ricci[m_, q_, g_, ing_] := 
      Block[{a}, 
        Expand[(∑_(a = 1)^(Length[g]))Riemann[a, m, a, q, g, 
              ing]]]

RicciScalar[g_, ing_] := 
        Block[{}, 
          Expand[(∑_(a = 1)^(Length[g]))(∑_(b = 1)^(Length[g]))ing[[a, b]] Ricci[a, 
                    b, g, ing]]];

Einstein[m_, n_, g_, ing_] := 
      Ricci[m, n, g, ing] - 
        (1/2) RicciScalar[g, ing] g[[m, 
            n]]

Bianchi[a_, g_, ing_] := 
        Block[{}, 
          Expand[(∑_(n = 1)^(Length[g]))(∂_(IndepVar[[n]]))((∑_(m = 1)^(Length[g]))ing[[n, m]] Einstein[m, a, g, 
                          ing]) + (∑_(n = 1)^(Length[g]))(∑_(m = 1)^(Length[g]))Christoffel[n, m, n, g, 
                      ing] (∑_(l = 1)^(Length[g]))ing[[m, l]] Einstein[l, a, g, 
                          ing] - (∑_(n = 1)^(Length[g]))(∑_(m = 1)^(Length[g]))Christoffel[n, m, a, g, 
                      ing] (∑_(l = 1)^(Length[g]))ing[[m, l]] Einstein[l, n, g, 
                          ing]]];

\(Simplify[
	
      Table[Christoffel[i, j, k, g, ing], {i, 1, 4}, {j, 1, 4}, {k, 1, 
          4}]]\)

\(Simplify[
	
      Table[Ricci[i, j, g, ing], {i, 1, 4}, {j, 1, 4}]]\)

Simplify[Ricci[1, 1, g, ing]]

Simplify[Ricci[2, 2, g, ing]]

Simplify[Ricci[3, 3, g, ing]]

Simplify[Ricci[4, 4, g, ing]]

Simplify[RicciScalar[g,  ing]]

Simplify[
      Table[Einstein[a, b, g, ing], {a, 1, 4}, {b, 1, 4}]]

The Einstein equations are given by

Simplify[
      Table[Einstein[a, b, g, ing]  -  8 π  G EnMomT[a, b, g], {a, 
          1, 4}, {b, 1, 4}]]

Simplify[
      Einstein[1, 1, g, ing]  -  
        8 π  G EnMomT[1, 1, g] ]

Simplify[
      Einstein[2, 2, g, ing]  -  
        8 π  G EnMomT[2, 2, g] ]

Simplify[
      Einstein[3, 3, g, ing]  -  
        8 π  G EnMomT[3, 3, g] ]

Simplify[
      Einstein[4, 4, g, ing]  -  
        8 π  G EnMomT[4, 4, g] ]

Simplify[
      Einstein[1, 2, g, ing]  -  
        8 π  G EnMomT[1, 2, g] ]

Check Bianchi identity:

Simplify[Table[Bianchi[a, g, ing], {a, 1, 4} ]]

In the case of a cosmological constant (λ not 0) one has for the Einstein equations instead

Simplify[
      Table[Einstein[a, b, g, ing]  +  λ g[[a, b]]  -  
          8 π G EnMomT[a, b, g], {a, 1, 4}, {b, 1, 4}]]

A solution of the above equations is given by the pair of equations

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

eq2 = k + R'[t]^2+ 2 R[t] R''[t] == - 8 π G p[t] R[t]^2

For non-relativistic matter with negligible pressure

p[t]  =  0

eq2 /. { R[t]  ->  (t^α), R'[t]  ->  D[(t^α), t] , R''[t] -> D[(t^α), {t, 2}]}

Case k=0  (otherwise α=1, k=-1) :

Solve[ 
      2  α (α - 1) +   (α^2)  == 
        0,  α]

eq1  /.   ρ[
          t]  ->  ρ0 R[t]^(γ )

%  /.  { 
        R[t]  ->  (t^α),  R'[t]  ->  
          D[(t^α), t]  }

Solve[ -2  +  
          2 α ==  α (2 + γ) ,  γ]

This then gives R(t) ~ (t^(2/3)) (α=2/3) and ρ(t) ~ R(t)\(TraditionalForm\`\^\(-3\)\) or ~ t\(TraditionalForm\`\^\(-2\)\) .

The Riemann tensor itself is singular at t=0

Simplify[ Riemann[1, 2, 1, 2, g, ing] ]

Simplify[ Riemann[2, 3, 2, 3, g, ing]]