Legal Notice, Disclaimer, and Book Ordering Information:

(* Legal Notice, Disclaimer, and Book Ordering Information:

Last Modified: FIO June 8, 2002

Mathematica For Physics: 2nd Edition

Robert L. Zimmerman and Fredrick I. Olness

WebSite:  
  www.physics.smu.edu/~olness
  darkwing.uoregon.edu/~phys600/

MathSource Number: 0206-862

ISBN 0-201-53796-6
For ordering information, call 1-800-282-0693
Addison-Wesley Publishing Company

Communication with the authors:
Fredrick I. Olness:  olness@mail.physics.smu.edu
Robert L. Zimmerman: bob@zim.uoregon.edu

Copyright 2002, Addison-Wesley Publishing Company, Inc. The material in
this file may be distributed freely  so long as the content remains
unchanged, and the copyright and  reference notices are included. The
publisher grants permission for the noncommercial use of these programs
and program segments.  All other uses require the prior written consent
of the publisher.    

All rights are reserved with regard to the  material in the
"Mathematica For Physics" text, and may not be reproduced, stored in a
retrieval system, or transmitted,  in any form or by any means,
electronic, mechanical, photocopying,  recording, or otherwise, without
the prior written permission  of the publisher. 

The programs and applications presented in this book have been included
for  their instructional value. They have been tested with care but are
not  guaranteed for any particular purpose. The publisher does not
offer any  warranties or representations, nor does it accept any
liabilities with  respect to the programs or applications. 

 *)

10. Relativity and Cosmology

Introduction

Off[General::spell ];  
Off[General::spell1];

10.1 Special Relativity

Overview of Special Relativity

Problem 1:  Decay of a particle

Remarks and outline

Solution

Clear["Global`*"];

momRule = 
{p1  -> 
            m1 {1, 0, 0, 0},  
 
          p2  -> 
            m2 γ2 {1, 0, 0, -v2}, 
p3  -> 
            m3  γ3 {1, 0, 0, v3}};

conRule = {p3 -> p1 - p2} ;

η = DiagonalMatrix[{1, -1, -1, -1}];

eq1= (m3^2) == p3 . η . p3 /. conRule /. momRule/.{v2->(Sqrt[1 - (1/γ)(2^2)])} //ExpandAll

sol =  Solve[eq1, γ2  ] // Simplify // Flatten

FullSimplify[m2  γ2   - m2 /. sol]

Problem 2:  Two-particle collision

Solution

Clear["Global`*"];

Part a

pL1 = {E1, p1, 0, 0}; 
pL2 = {m2, 0, 0, 0};

pLT = pL1 + pL2

pCM1 = {ECM1, kCM1, 0, 0};
pCM2 = {ECM2, -kCM1, 0, 0};

pCMT = pCM1 + pCM2

Part b

η = DiagonalMatrix[{1, -1, -1, -1}];

eq1 = pLT . η . pLT == pCMT . η . pCMT

eq2 = ECM1 + ECM2 == m3 + m4;

{eq3, eq4}={E1 == T1 + m1,p1==(Sqrt[(E1^2) - (m1^2)])}

sol = Solve[{eq1, eq2, eq3, eq4}, {T1}, {E1, p1, ECM1}] // Flatten // 
        Simplify

Problem 3:  Compton scattering

Solution

Clear["Global`*"];

Part a

p1=(h/λ)0 {1, 0, 0, 1}; 
p2 = m {1, 0, 0, 0};

p3=(h/λ)1{1, Sin[θ1], 0, 
            Cos[θ1]};
p4 = 
        m γf {1, -vf Sin[θ2], 0, 
            vf Cos[θ2]};

η = DiagonalMatrix[{1, -1, -1, -1}];

eq1 = ((m^2) == (p1 + p2 - p3) . η . (p1 + p2 - 
                    p3) // Simplify // ExpandAll)

lamRule = Solve[ eq1 ,  λ1 ][[1, 
            1]] //ExpandAll

λ1  - λ0 /. 
        lamRule//Collect[#, {m, h}] &

Part b

(eq2 = p1 + p2 == p3 + p4 // Thread) // 
      ColumnForm

eq3 = Solve[{eq2[[1]]}, γf  ] // 
          Flatten//Simplify

((  m (γf - 1)  /. lamRule /. eq3 // 
                  Factor)  /. lamRule  /. λ0 -> 
              h/m temp  // Factor) /. 
      temp ->  λ0 m/h

Part c

eq4 = eq2[[{2, 4}]]

eq5 = Eliminate[ 
          eq2[[{2, 
                4}]],   γf ]//Simplify

eq6 =  
      eq5 /. lamRule  // ExpandAll[#, Trig -> True] &

Solve[
              eq6 /. { 
                  Sin[θ2] -> 
                    tan[θ2] Cos[θ2]} ,   tan[θ2]][[1, 1]]  /. tan  -> Tan /.   Rule -> 
          Equal //Simplify

Problem 4:  Moving mirror and generalized Snell's law

Solution

Clear["Global`*"];

Part a

Ki = νi {1, Cos[θi], Sin[θi], 0}; 
    Kr = νr {1, -Cos[θr], Sin[θr], 0};

Kmi = νm {1, Cos[θm], Sin[θm], 0}; 
    Kmr = νm {1, -Cos[θm], Sin[θm], 0};

boostx=(({{γ, -v γ, 0, 0}, {-v γ, γ, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1}}));

(eq1 = (Ki == boostx . Kmi) // 
            Thread)//ColumnForm

(eq2 = (Kr == boostx . Kmr) // Thread)  // 
      ColumnForm

sol1 = Solve[{eq1, eq2} // Flatten,  νi , {Cos[θr], 
            Sin[θr], Cos[θm], 
            Sin[θm], νm}]//Flatten

Part b

eq3 = 
        Eliminate[{eq1, 
                  eq2  ,   νr ≠ 0, νi ≠ 
                    0  ,  v ≠ 0, v ≠ -1  , 
                  v ≠ 1 } /. sol1 // Flatten, {Cos[θm], 
              Sin[θm], γ  }]  // 
          Map[Simplify , #] & ;

    
 eq3 //ColumnForm

eq4 = Solve[ eq3[[2]] , v ]  // Simplify // 
          PowerExpand  // Simplify

((1 - v)/(v + 1))==(((1 - v)/(v + 1))/.eq4[[2]])//Simplify

Problem 5:  One-dimensional motion of a relativistic particle with constant acceleration

Solution

Clear["Global`*"];

Part a

eqMotion=D[((x'[t])/(Sqrt[1-((x'(t))^2)])),t]==g

eq1  =  DSolve[{eqMotion } // Flatten , x[t], t] // 
              Simplify // PowerExpand // 
          Simplify //Flatten

eqMotion /. D[eq1[[#]], t] /. 
                D[eq1[[#]], {t, 2}] & /@ {1, 2}  // Simplify // 
        PowerExpand //Simplify

eq2 = x[t] /. eq1[[1]]

initialConditions = {( eq2) == x0, (D[eq2, t]) == 
              v0} /. t -> 0 // PowerExpand

gamRule={1 - (v0^2) -> (1/(γ^2)), 
            (v0^2) - 1 -> -(1/(γ^2))};

eq3 = ((Solve[initialConditions, {C[1], C[2]}]   // 
                  Simplify) // PowerExpand) /. 
          gamRule //PowerExpand

crule = eq3[[1]]

xsol =   eq2 /. crule//Simplify

nonRelResult =Series[( xsol/.γ->(1/(Sqrt[1 - (v0^2)]))) ,   {g, 0, 1},{v0, 0, 1}]//Normal

Part b

values={x0 -> 0,v0 -> 1/2,g -> -0.1,γ->(1/(Sqrt[1 - (v0^2)]))};

text = {Text[ 
            StyleForm["Relativistic", FontFamily -> "Times-Bold", 
              FontSize -> 10 ], {10, 1.2}], 
          Text[ StyleForm["non-relaltivistic ", 
              FontFamily -> "Times-Bold", FontSize -> 10 ], {8,  
              0.2}] };

\(\(\(Plot[{\ xsol, nonRelResult} //. values // 
          Evaluate, 
{t, 0, 
          12}, 
GridLines -> 
          Automatic, 
AxesLabel -> {"t", "x"}, 
        PlotStyle -> {{Thickness[0.01]}, {Dashing[{0.06, 0.1}], 
              Thickness[
                0.001]}}, 
Epilog -> \ \ text\ ];\)\(\ \ \)\)\)

Problem 6: Two-dimensional motion of a relativistic particle in a uniform electric field

Solution

Clear["Global`*"];

Part a

pVec[t_] = {px[t], py[t], pz[t]};

eqMotion = pVec'[t] == {0,  e E0, 
            0}//Thread

initial = {px[0] == px0, py[0] == 0, 
          pz[0] == 0};

momRule = DSolve[eqMotion~Join~initial, {px[t], py[t], pz[t]}, 
          t][[1]]

γRule = m γ == 
            Sqrt[pVec[t] . pVec[t] + (m^2) /. momRule] // 
          Solve[#, γ] & // Flatten

Part b

eq1 =          pVec[t] == 
              m γ {x'[t], y'[t], z'[t]} /. 
            momRule /. γRule // Thread

initial = {x[0] == 0, y[0] == 0, 
          z[0] == 0};

path = 
        DSolve[{eq1, initial} // Flatten, {x[t], y[t], z[t]}, t] // 
          Flatten;

    
path // ColumnForm

values = {m -> 1, e -> 1, E0 -> 0.1, 
          px0 -> 0.1};

Plot[{x[t], y[t]} /. path //. values // 
          Evaluate, 
{t, 0, 
          5}, 
GridLines -> Automatic, 
        PlotStyle -> { Thickness[0.01] ,  Dashing[{0.06, 0.1}] }, 
        AxesLabel -> {"t", "y or x"} ];

relXY = {x[t], y[t]} /. path;

nonrelXY = (Map[Series[#, {t, 0, 2}] &, 
                  relXY] /. {px0^2 -> 0} // PowerExpand // Simplify // 
          Normal)

Plot[{relXY[[1]],  nonrelXY[[1]]} /. values // 
          Evaluate, 
{t, 0, 
          30}, 
GridLines -> Automatic, 
        PlotStyle -> { Thickness[0.01] , Dashing[{0.06, 0.1}] }, 
        AxesLabel -> {"t", "x"} ];

Plot[{relXY[[2]],  nonrelXY[[2]]} /. values // 
          Evaluate, 
{t, 0, 
          30}, 
GridLines -> Automatic, 
        PlotStyle -> { Thickness[0.01] ,  Dashing[{0.06, 0.1} ] }, 
        AxesLabel -> {"t", "y"} ];

10.2 General Relativity

Overview of General Relativity

Problem 1: Killing vectors and Spherical Symmetry

Clear["Global`*"];

gμν=DiagonalMatrix[{g00(r), 
              g11(r), -(r^2), -(r^2) (sin^2)(θ)}];
gμν // MatrixForm

Remarks and outline

killingEq[ContraV_, g_, Γ_, 
          x_: {t, r, θ, ϕ}] := 
Module[{temp1, 
            temp2, temp3, temp4}, 
temp1 = 
            g . ContraV; 
temp2 = 
            Outer[D, temp1, x] - 
              temp1 . Γ; 
temp3 = 
            temp2 + Transpose[
                temp2]; 
temp4 = (temp3 // 
                      Flatten) == 0 // Thread  // 
              Simplify; 
Complement[
            Select[temp4, 
              FreeQ[#, True & ] ] ]
];

christ[met_, 
          var_: {t, r, θ, ϕ}]  := 
Module[{ 
            temp1, imet, result}, 
temp1 = 
            Outer[D, met, var] ; 
imet = 
            Inverse[met] // Simplify; 
result = 
            1/2 ( Transpose[imet . temp1, {1, 3, 2}]  +    Transpose[
                    imet . temp1, {1, 2, 3}]   -     imet . 
                    Transpose[
                      temp1, {3, 2, 1}] ); 
Return[
            result]];

Solution

Part a

kt = {1, 0, 0, 0};
kz = {0, 0, 0, -1};
kx = {0, 
          0, -Sin[ϕ], -Cos[ϕ] 
            Cot[θ]};  
ky = {0, 0, 
          Cos[ϕ], -Sin[ϕ] Cot[θ]};

Γ=christ[gμν]

Γ

(killingEq[#, gμν, Γ] & /@ {kt, kz, kx, 
            ky}  // Simplify)

Part b

velocity = {t'[s], r'[s], θ'[s], ϕ'[
            s]};

kconst = {
en == 
            gμν . velocity . kt, 
jx == 
            gμν . velocity . kx, 
jy == 
            gμν . velocity . ky, 
jz == 
            gμν . velocity . kz};

    
kconst // ColumnForm

Solve[kconst[[{ 1, 4}]], {t'[s], ϕ'[s]}] // 
        Flatten // ColumnForm

Solve[kconst[[{ 2, 3}]], θ'[s], ϕ'[s] ] // 
        Flatten//Simplify

Problem 2: Schwarzschild solution

Clear["Global`*"];

gμν=DiagonalMatrix[{g00(r), 
                g11(r), -(r^2), -(r^2) (sin^2)(θ)}];

Remarks and outline

christ[met_, 
          var_: {t, r, θ, ϕ}]  := 
Module[{ 
            temp1, imet},  
temp1 = 
            Outer[D, met, var] ; 
imet = 
            Inverse[met] // Simplify; 
1/2 ( 
              Transpose[imet . temp1, {1, 3, 2}]  +    Transpose[
                  imet . temp1, {1, 2, 3}]   -     imet . 
                  Transpose[
                    temp1, {3, 2, 1}] )
];

curvR[christ_, 
        var_: {t, 
            r, θ, ϕ}] := 
Module[{temp1}, 
temp1 = 
          Outer[D, christ, 
            var]; 
(Transpose[
              Outer[D, christ, var], {1, 2, 4, 3}] - 
            Transpose[Outer[D, christ, var], {1, 2, 3, 4}] + 
            Transpose[Inner[Times, christ, christ, Plus, 2], {1, 3, 4, 2}] - 
            Transpose[
              Inner[Times, christ, christ, Plus, 2], {1, 4, 3, 
                2}])
]

einstein[curv_, 
        met_] := 
Module[{dem, temp1, temp2, i, j, 
          s}, 
dem = 
          Length[met]; 
temp1 = 
          Simplify[
            ExpandAll[
              Table[Sum[curv[[s, i, s, j]], {s, dem}],    {i, 
                  dem}, {j, dem}]]]; 
temp2 = 
          Simplify[
            ExpandAll[
              Sum[Inverse[met][[j, i]] curv[[s, i, s, j]], {s, 
                  dem}, {i, dem}, {j, dem}]]]; 
Simplify[
          ExpandAll[
            Table[temp1[[i, j]] - (1/2) met[[i, j]] temp2, {i, 
                dem}, {j, dem}]]]]

Solution

Part a

Γ =christ[
        gμν, {t, r, θ, ϕ}]

R = 
        curvR[Γ, {t, r, θ, ϕ}] // 
          Simplify ;

Gμν = einstein[R, gμν];

{eqGtt, eqGrr, eqGθθ, 
          eqGϕϕ} =  MapAt[Flatten, Gμν == 0, 1] // 
              Thread // Select[#1, FreeQ[#, True] &] &  // 
          Simplify ;

    
{eqGtt, eqGrr, eqGθθ, eqGϕϕ}  // 
      ColumnForm

vacuumeq = ( { eqGtt, eqGrr} // 
                Solve[#, {g11'[r], g00'[r]}] & // Flatten) /. 
          Rule -> Equal;

    
vacuumeq // ColumnForm

sol=DSolve[vacuumeq,{g11[r],g00[r]},r]// Simplify//Flatten

bdcond={C[1]->Log[((c^2)/(2 G M))],C[2]-> ((2 G M)/(c^2))}

(gμν /. sol /. bdcond /. G -> temp (c^2) r  // 
          Simplify) /. temp -> (G/(r (c^2)))

Part b

data = {{7.35 *(10^25)  g, 
              1.738 *(10^8)   cm},  
           {5.977 *(10^27)  g, 
              6.378 *(10^8)  cm}, 
           {Mo, 6.9598 *(10^10)    cm}, 
           {Mo , 1.738 *(10^10)  cm}, 
           {Mo , 
              (10^6)  cm }} /. 
Mo ->  
            1.989 *(10^33) g  ;

names = {Moon, Earth, Sun, "white dwarf" , 
          neutron star};
 TableForm[data, TableHeadings -> {names,  { Mass,  Radius}}] /. 
      cm -> (10^(-5)) km

constants = {c -> 3.0* (10^10) cm/sec ,  
          G -> 7.425 *(10^(-29)) (((c^2) cm)/g)};

Ratio[{M_, r_}] :=((2 G M)/((c^2) r))//.constants

Transpose[{names, Ratio[#] & /@ data }]  // TableForm

Problem 3: Geodesics for the Schwarzschild metric

Clear["Global`*"];

Remarks and outline

geodesic[met_, var_, 
        parameter_:  s] := 
Module[{met1, temp1, vel1, 
          var1}, 
vel1 = 
          D[Through[var[parameter]], parameter]; 
var1 = 
          Through[var[parameter]]; 
met1 = 
          met /. Thread[
              var -> 
                Through[
                  var[parameter]]]; 
temp1 = (D[
                    met1 . vel1, 
                    parameter] - (1/2 D[met1, #] . vel1 . vel1 & /@ 
                      var1) == {0, 0, 0, 0}) // Thread]

Solution

Part a

gμν=DiagonalMatrix[{g00(r), 
            g11(r), -(r^2), -(r^2) (sin^2)(θ)}]

comp=\({g00 -> \((1 - \(2\ m\)\/#1 &)\), 
            g11 -> \((\(-\(1\/\(1 - \(2\ m\)\/#1\)\)\) &)\)}\);

eq1 = 
        geodesic[gμν , {t, r, θ, ϕ}] // 
          Simplify ;

eq2 = (Solve[
                      eq1, {t''[s], r''[s], θ''[
                          s], ϕ''[s]}] // Flatten) /. 
                comp /. θ -> (π/2 &) /. 
            Rule -> Equal  // Simplify;

    
 eq2//ColumnForm

Part b

eq3 =  DSolve[eq2[[1]] , r[s], s][[1, 
            1]]  /. Rule->Equal

trule= Solve[eq3, t'[s]][[1, 1]]/.C[1]->-(1/(2 m))  Log[En]//Flatten

ϕrule=(DSolve[eq2[[4]] , r[s], s][[1, 1]]  /. 
              Rule -> 
                Equal // Solve[#, ϕ'[s]][[1, 1]] & )/.C[1]->(Sqrt[h])

eq4  = eq2[[2]] /. trule /. ϕrule // Simplify

Part c

vel = {t'[s], r'[s], 0, ϕ'[s]};

met = (gμν /. comp /. θ -> (π/2) /. 
              r -> r[s] /. θ -> θ[
                s]) ;

eq5 = 1 == vel . met . vel

eq6 = eq5 /. trule /. ϕrule  // Simplify

r'[s]^2 /. Solve[eq6, r'[s]][[1]] // 
      Simplify

Problem 4: Time it takes to fall into a black hole

Solution

Clear["Global`*"];

Part a

eq1=((r'(s))^2)==(En^2) + ((2 m)/(r(s))) - 
              ((h^2)/((r(s))^2)) + ((2 (h^2) m)/((r(s))^3)) - 1;

eq2=Solve[eq1, r'[s]]  /. h -> 0  /. 
            Rule -> Equal /.(Sqrt[x])_ :> (Sqrt[Simplify])[x]//Flatten

eq3 = eq2[[1]]

eq4 = 0 == eq3[[2]] /. s -> 0 /. 
        r[0] -> r0

rprimerule =  Solve[{eq3, eq4}, {r'[s], En}][[2]]    // 
          Simplify // ExpandAll

radialeq=rprimerule[[1]]/.(Sqrt[x])_ :> (Sqrt[Factor])[x]//PowerExpand

Stime = Integrate[ 
          (1/(r'))[s] /. rprimerule, {r[s], r0 , 
            r }, 
Assumptions -> {  r0  > r  > 
                0}]//Apart

Stime /. r -> 0/.m-> ((G M)/(c^2))//PowerExpand

Part b

trule=t'[s]->(En/(1 - ((2 m)/(r(s)))));

eq5=(((r'(s))/(t'(s)))/.trule/.rprimerule/.r[s] -> r[t])  /. \(x_\^n_ :> Factor[Together[x]]\^n\) //PowerExpand //Expand //Simplify

change= {r(t) -> r0 (cos^2)(x),r'(x)->-2 r0 cos(x) sin(x)};

eq6 =  ((1/eq5)  r'[x] /. change) // 
          Simplify//PowerExpand

Ttime=(Integrate[  eq6 ,{x,0,(cos^(-1))((Sqrt[(r/r0)])) }  ]/.{2 m - r0 -> -temp1 }//PowerExpand)/.temp1 -> r0 - 2 m //Simplify

Part c

value = {m -> 1, r0 -> 8};

pt1 = 
        Plot[Evaluate[(Ttime /. value)], {r, 2, 8}, 
          DisplayFunction -> Identity];

    
pt2 = 
        Plot[Evaluate[Stime /. value], {r, 0, 8}, 
          DisplayFunction -> Identity];

text = {Line[{{2, 0}, {2, 
                  40}}], 

 
            Text[ StyleForm["Proper Time", 
                FontFamily -> "Times-Bold", FontSize -> 10 ], {1.2, 
                22}], 
 
Text[ 
              StyleForm["Coordinate Time", 
                FontFamily -> "Times-Bold", FontSize -> 10 ], {4.5, 
                30}], 
 
Text[ 
              StyleForm["r=2m", FontFamily -> "Times-Bold", 
                FontSize -> 10 ], {2.6, 
                10}] 

} // 
          Graphics;

\(\(Show[{pt1, pt2, \ 
          text}, 
AxesLabel -> {"Distance", "Time"}, 
        DisplayFunction -> $DisplayFunction, 
PlotRange -> {0, 40}];\)\)

Problem 5: Circular geodesics for the Schwarzschild metric

Clear["Global`*"];

Solution

teq=t'(s)==(En/(1 - ((2 m)/(r(s)))))/.r[s]  -> r0

ϕeq= ϕ'(s)==(h/((r(s))^2))/.r[s]  -> r0

acc=r''[s] == ((m ((En^2)-((r'(s))^2)) ((r(s))^3)-(h^2) ((r(s))^2)+4 (h^2) m (r(s) - m))/((2 m - 
                      r(s)) ((r(s))^4)))/.r  -> (r0 &)

vel= ((r'(s))^2)==(En^2) + ((2 m)/(r(s))) - 
              ((h^2)/((r(s))^2)) + ((2 (h^2) m)/((r(s))^3)) - 1 /.r  -> (r0 &)

eq2= (Solve[{acc, vel}, {En,  h} ] // Simplify  // 
            PowerExpand)/.{(1/(Sqrt[3 m - r0])) -> -(I/(Sqrt[r0 - 3 m]))}

Enhrule = eq2[[4]]

sol = DSolve[{ teq, ϕeq, 
                t[0] == 0, ϕ[0] == 0}, {t[s], ϕ [s]}, 
              s] /. Enhrule // Flatten // Simplify

s/.Solve[(ϕ[s] /. sol)  == 
            2 π, s][[1, 1]]

Period/.Solve[ {Period == 
                t[s ], ϕ[s] == 2 π } /. sol, Period, 
          s][[1, 1]]

Problem 6: Precession of the Perihelion

Clear["Global`*"];

Solution

Part a

eq1=((r'(s))^2)==(En^2) + ((2 m)/(r(s))) - 
              ((h^2)/((r(s))^2)) + ((2 (h^2) m)/((r(s))^3)) - 1;

change={r'(s)->((h r'(ϕ))/(((r(ϕ))^2))),r(s) -> r(ϕ)};

eq2 = eq1 /. change

eq3 = eq2 /. {r -> (1/u[#] &)}

uEq = (D[
                eq3, ϕ] // Solve[#, u''[ϕ]][[1, 
                    1]] &) /. Rule -> Equal // ExpandAll

Part b

pertu = {u  -> (u0 + du[#] &)};

uEq /. pertu//ExpandAll

eq4 =  (uEq /. pertu // ExpandAll)  /. 
        du[ϕ]^2 -> 0

eq5 = eq4 /. du''[ϕ] -> 0 /. 
          du[ϕ] -> 0 // Simplify

eq6 = Solve[{eq5, eq4}, du''[ϕ]][[1, 
                  1]] /. {6 m u0 -> 1 - (γ^2)} /. Rule -> Equal // 
          ExpandAll // Simplify

γrule=  γ -> (Sqrt[1 - 6 m u0]);

dusol=DSolve[{eq6  , du[0] == e u0, du'[0] == 0}, 
          du[ϕ], ϕ][[1, 1]]

uSol = u0 + du[ϕ]/.dusol

u0rule= u0 -> (1/(a (1 - (e^2))));

r == (1/ uSol  ) /. 
        u0rule  //Simplify

eq7 =  2 π/γ  -  2 π;

((eq7 /. γrule /. u0rule) + O[m]^2 // Normal)/.m -> ((G M)/(c^2))

Δ=((6 G M π)/(a (c^2) (1 - (e^2))));

Part c

const={M  -> 1.989 *(10^33)  gm, G -> 7.42426 * (10^(-29)) (((c^2) cm)/gm)};

mercury  = { a -> 0.3871 AU,  
          AU -> 1.495985* (10^13) cm,   e -> 
            0.20561421 };

radTosec =  360* 60 *60 /(2 π) ;

Mnumber =    (100/0.24085);

Δ (radTosec) ( Mnumber) /. const //. 
        mercury//NumberForm[#, 4] &

venus = { a -> 108.2  *(10^11) cm, 
          e -> 0.007};

Vnumber = (100/0.62)  ;

Δ (radTosec) ( Vnumber) /. 
          const //. venus // NumberForm[#, 2] &

earth = {  a -> 149.6 * (10^11) cm , 
          e -> 0.017};

Δ (radTosec ) (100) /. const //. earth // 
      NumberForm[#, 2] &

10.3 Cosmology

Overview of Friedmann, Robertson, and Walker Cosmology

Problem  1: Field equations for Friedmann-Robertson-Walker cosmology

Clear["Global`*"];

Remarks and outline

christ[gμν_, 
          var_: {t, r, θ, ϕ}]  := 
Module[{ 
            temp1, imet}, 
temp1 =  
            Outer[D, gμν, var] ; 
imet =  
            Inverse[gμν] // Simplify; 
  1/2 ( 
              Transpose[imet . temp1, {1, 3, 2}]  + 
                Transpose[imet . temp1, {1, 2, 3}]   -  
                imet . Transpose[
                    temp1, {3, 2, 1}] )
];

curvR[Γ_, 
        var_: {t, r, θ, ϕ}] := 
Module[{temp1}, 
        temp1 = Outer[D, Γ, 
            var]; 
(Transpose[
              Outer[D, Γ, var], {1, 2, 4, 
                3}] - 
Transpose[
              Outer[D, Γ, var], {1, 2, 3, 
                4}] + 
Transpose[
              Inner[Times, Γ, Γ, Plus, 2], {1, 3, 
                4, 2}] - 
Transpose[
              Inner[Times, Γ, Γ, Plus, 2], {1, 4, 
                3, 2}])
]

einstein[curv_, 
        gμν_] := 
Module[{dem, temp1, temp2, i, j, 
          s}, dem = Length[gμν]; 
temp1 = 
          Simplify[
            ExpandAll[
              Table[Sum[curv[[s, i, s, j]], {s, dem}], {i, dem}, {j, 
                  dem}]]]; 
temp2 = 
          Simplify[
            ExpandAll[
              Sum[Inverse[gμν][[j, i]] 
                  curv[[s, i, s, j]], {s, dem}, {i, dem}, {j, 
                  dem}]]]; 
Simplify[
          ExpandAll[
            Table[temp1[[i, j]] - (1/2) gμν[[i, j]] 
                  temp2, {i, dem}, {j, dem}]]]]

Solution

Part a

gμν=DiagonalMatrix[{1, -(((R(t))^2)/(1 - 
                      k (r^2))), -(r^2) ((R(t))^2), -(r^2) ((R(t))^2) (sin^2)(θ)}];

Γ = christ[gμν];

Rαβγδ = 
        curvR[Γ] // Simplify;

Gμν = 
      einstein[Rαβγδ, gμν]

vel = {1, 0, 0, 
          0};                 
covel = gμν . {1, 0, 0, 0};

T = 
        Table[(ρ[t] + P[t]) covel[[i]] covel[[j]] - 
              gμν[[i, j]] P[t], {i, 4}, {j, 4}] // 
          Simplify;

(eq1 = Map[Flatten, 
                  Gμν == 
                    8 π G T + λ gμν] // Thread // 
              Select[#1, FreeQ[#, True] &] &  // Simplify ) // 
      ColumnForm

(fieldeq = Solve[
                    eq1 /. R'[t] -> Sqrt[x], {R''[t], 
                      x}][[1]] /. x -> R'[t]^2 /. 
              Rule -> Equal // ExpandAll) // ColumnForm

Part b

conT = 
      Table[(ρ[t] + P[t]) 
            vel[[i]] vel[[j]] - Inverse[gμν][[i, 
                j]] P[t], {i, 4}, {j, 4}]

\(coder[{2, U, U}, conT_, christ_, var_]\  := 
	
      Outer[D, conT, var] + \ 	
	
        Transpose[
          Inner[Times, christ, conT, Plus, 2], {1, 3, 2}] + \ \ Transpose[
          Inner[Times, christ, Transpose[conT, {2, 1}], Plus, 2], {2, 3, 
            1}]\)

eq2 = 
        coder[{2, U, U}, conT, Γ, {t, r, θ, ϕ}] // 
          Simplify ;

eq3=Table[
        Sum[eq2[[i, s, s]], {s, 1, 4}], {i, 1, 4}]

eq4 = eq3[[1]] == 0

eqstate = {P[t] -> α ρ[t]};

eq5 = DSolve[{eq4 /. eqstate, R[to] == Ro}, R[t], t][[1, 
                1]] /. {ρ[to] -> ρo, Rule -> Equal} // 
          Simplify // PowerExpand

ρrule = Solve[eq5, ρ[t]][[1, 1]] // PowerExpand // 
        Simplify

Part c

fieldeq/.ρrule/.eqstate/.α->0//
                           ExpandAll  //ColumnForm

Problem  2: Zero-pressure (Dust) Cosmological Models

Remarks and outline

Solution

Clear["Global`*"];

Part a

eq1=R''(t)==(1/3) λ R(
                  t) - ((4 G π (Ro^3) ρo)/(3 ((R(t))^2)));
eq2=((R'(t))^2)==((8 G π ρo (Ro^3))/(3 R(t))) + 
              (1/3) λ ((R(t))^2) - k;

parameterrule={R(to) -> Ro,R'(to)->Ho Ro,R''(to)->-(Ho^2) qo Ro,ρo -> ((3 (Ho^2) σ)/(4 G π)),λ -> 3 (Ho^2) Γ};

qrule = Solve[{eq1, eq2} /. t -> to /. parameterrule, {qo, 
              k}] // Flatten // Simplify

Part b

change = {R'[t] -> X'[t] Ro, 
          R[t] -> X[t] Ro};

(eq2 /. parameterrule /. qrule /. change // 
        Solve[#, X'[t]] &)

eq3 = {X'[t] -> 
          Ho (1 - Γ - 
                  2 σ + (2 σ)/
                    X[t] + Γ X[t]^2)^(1/2)}

integrand=(1/X'[t] /. 
            eq3 /. {Γ -> 0, 
            Ho -> 1, σ -> qo}  )

eq4 = 
      Integrate[ integrand, {X[t], 0, 1} , 
          GenerateConditions -> False ] // Simplify

eq4 = Integrate[ integrand, {X[t], 0, 1}  ] // Simplify

Clear[ageo]

ageo[q1_]  :=  
        eq4[[2]] /. qo -> q1 // Chop;
ageo[1/2] = 2/3;  ageo[ .5] = 2/3;
ageo[0] = 1; geo[ .0] = 1;

label = {PointSize[ .03], Point[{ .5, ageo[ .5]}], 
          Text[ StyleForm["k=0, Flat", FontFamily -> "Times-Bold", 
              FontSize -> 10 ], { .9,  .67}], 
          Text[ StyleForm["k=1, Closed", FontFamily -> "Times-Bold", 
              FontSize -> 10 ], {1.4,  .58}], 
          Text[ StyleForm["k=-1, Open", FontFamily -> "Times-Bold", 
              FontSize -> 10 ], { .5,  .82}]};

\(\(\(\ \)\(ageplot[Ho_] := 
        Plot[ageo[qo], {qo,  .00001, 2}, Frame -> True, 
          PlotLabel -> "Age for Γ=0", 
          FrameLabel -> {"10^-29 gm/cc", "  ", "ϱo", " "}, 
          FrameTicks -> \((1/10 
                   Round[\ 
                    10\ {{0, {1/2,  .000196\ Ho^2}, {1, 
                            2* .000196\ Ho^2}, {1.5, 
                            3* .000196\ Ho^2}}, {{ .5,  .5*9.78*\((100/
                                  Ho)\)}, { .6,  .6*9.78*\((100/
                                  Ho)\)}, { .7,  .7*9.78*\((100/
                                  Ho)\)}, { .8,  .8*9.78*\((100/
                                  Ho)\)}, { .9,  .9*9.78*\((100/Ho)\)}, {1, 
                            1*9.78*\((100/Ho)\)}}}\ \ ] // N)\), 
          GridLines -> Automatic, Epilog -> label];\)\)\)

ageplot[70];

Part c

ageL[
        q1_, σ1_]:=NIntegrate[1/X'[t] /. eq3 /. 
                  Ho -> 1 /. Γ -> σ1 - q1 /. 
              qo -> q1 /. σ -> σ1 /. 
          X[t] -> x, {x, 0, 1}]

agetable[Ho_] := 
        1/100  Round[
              100  Table[ 
                  9.78*(100/Ho) ageL[q, sigma] , {q, 0, 
                    2,  .5}, {sigma, 0, 2.5,  .5}]] // N;

\(TableForm[agetable[70], TableSpacing -> {1}, 
      TableHeadings -> {{"qo=0", "qo=.5", "qo=1", "qo=1.5", "qo=2"}, {"σo=0", ".5", "1", "1.5", "2", "2.5"}}]\)

Problem  3: The Expansion and Age for the Friedmann-Robertson-Walker models

Remarks and outline

Solution

Clear["Global`*"];

Part a

eq1=((R'(t))^2)==(((Ho^2) (Ro^3))/(R(t)));

eq2 = 
      t == Integrate[Sqrt[1/eq1[[2]]], R[t]] // 
        PowerExpand

eq2[[2]] /. R[t] -> Ro

scale0 = Solve[eq2, R[t]][[2, 1]] /. 
        Rule -> Equal

Plot[
        Evaluate[scale0[[2]] /. {Ho -> 1, Ro -> 1}], {t, 0, 1}, 
        Frame -> True, 
        FrameTicks -> {{{0, "0"}, {1/3, "1/(3Ho)"}, {2/
                  3, "Now "}, {1, "1/Ho"}}, {0, {1/
                  2, " Ro/2 "}, {1, "Ro"}}}, 
        Epilog -> Graphics[{PointSize[ .03], 
                Point[{2/3, 1}]}][[1]]];

Plot[
        Evaluate[scale0[[2]] /. {Ho -> 1, Ro -> 1}], {t, 0, 1}, 
        Frame -> True, PlotLabel -> scale0, 
        FrameTicks -> {{{0, "0"}, {1/3, "1/(3Ho)"}, {2/
                  3, "Now "}, {1, "1/Ho"}}, {0, {1/
                  2, " Ro/2 "}, {1, "Ro"}}}, 
        Epilog -> Graphics[{PointSize[ .03], 
                Point[{2/3, 1}]}][[1]]];

Part b

eq3 = R'[t]^2 == -1 + 
            2 qo Ho^2 Ro^3/(R[t]);

Rorule = Ro -> 1/Ho Sqrt[1/(2 qo - 1)];

eq4 = eq3 /. Rorule // PowerExpand

change = 
        R -> ((2 qo/(Ho (2  qo - 1)^(3/2)) 
                  Sin[x[#]]^2) &);

scale1 = R[x] /. change /. x[x] -> x // PowerExpand

eq5 = Solve[eq4 /. change, x'[t]][[2]] // Simplify // 
        PowerExpand

time1=Integrate[(1/x'[t] /. eq5), {x[t], 0, 
          x}]

xage = Solve[Ro == scale1, x][[2]] /. Rorule // 
          ExpandAll // PowerExpand

age1 = time1 /. xage // ExpandAll // Simplify

value = {qo -> 1, Ho -> 1};

label = {PointSize[ .03],  
          Text[ StyleForm[
              1/100  Round[100 (time1 /. xage /. value  ) ] // 
                N  , FontFamily -> "Times-Bold", 
              FontSize -> 10 ], {2.6, 1}], 
 
          Text[ StyleForm["age/Ho ="  , 
              FontFamily -> "Times-Bold", FontSize -> 10 ], {1.4, 
              1}], 
Point[{time1 /. x -> π/2 /. 
                value, scale1 /. x -> π/2 /. value}], 
          Point[{time1 /. xage /. value, scale1 /. xage /. 
                value}]};

ParametricPlot[Evaluate[{time1, scale1} /. value], {x, 0, π}, 
        Frame -> True, PlotLabel -> "Closed Model, qo=1 ", 
        FrameTicks -> {{{0, "0"}, {age1 /. value // 
                  N, "age"}, {3 age1 /. value // 
                  N, "3* age"}, {time1 /. value /. 
                  x -> π/2, "Max"}, {9 age1 /. value // 
                  N, "9* age"}}, {0, 1/2, {1, "Ro"}, 1.5, 2}}, 
        Epilog -> label];

Part c

eq6=R'[t]^2 == 1+ 2 qo Ho^2 Ro^3/(R[t]);

Rorule = Ro -> 1/Ho Sqrt[1/(1 - 2 qo)];

eq7 = eq6 /. Rorule // PowerExpand

change = 
        R -> ((2 qo/(Ho (1 - 2 qo)^(3/2)) 
                  Sinh[x[#]]^2) &);

scale2 = R[x] /. change /. x[x] -> x // PowerExpand

eq8 = Solve[eq7 /. change, x'[t]][[2]] // Simplify // 
        PowerExpand

time2 = Integrate[(1/x'[t] /. eq8), {x[t], 0, x}]

xage = Solve[Ro == scale2, x][[2]] /. Rorule // 
          PowerExpand // Simplify

age2 = time2 /. xage

value = {qo ->  .02, Ho -> 1};

label = {PointSize[ .03], 
 
          Text[ StyleForm["age/Ho = "  , 
              FontFamily -> "Times-Bold", 
              FontSize -> 10 ], {1.5 time2 /. xage /. 
                value, scale2 /. xage /. value}], 
 
          Text[ StyleForm[
              1/100  Round[100 (time2 /. xage /. value  ) ] // 
                N  , FontFamily -> "Times-Bold", 
              FontSize -> 10 ], {1.9 time2 /. xage /. 
                value, scale2 /. xage /. value}], 
          Point[{time2 /. xage /. value, scale2 /. xage /. 
                value}]};

ParametricPlot[Evaluate[{time2, scale2} /. value], {x, 0, 3}, 
        Frame -> True, PlotLabel -> "Open Model, qo=.02", 
        FrameTicks -> {{{0, " 0"}, { .5 age1 /. value // 
                  N, "age/2"}, {age1 /. value // 
                  N, "age"}, {2  age1 /. value // 
                  N, "2age"}, {3/2 age1 /. value // 
                  N, "3/2 age"}}, {0, 1/2, {1, "Ro"}, 1.5, 2}}, 
        Epilog -> label];

10.4 Exercises

Exercise 1 Lorentz boosts

Exercise 2 Radioactive decay of moving Nucleus

Exercise 3 Schwarzschild solution

Exercise 4 Potential analysis for timelike geodesics

Exercise 5 Schwarzschild solution in null coordinates

Exercise 6 Null geodesics for the Schwarzschild metric

Exercise 7 Mercury's Anomalous Precession and Modification of Newton's Gravity

Exercise 8 Cosmological Force and Mercury