The plan is: 1. To have a generalized form of Ricci tensor and Ricci scalar. 2. Substitute into Einstein’s equation to find solution.
To begin, Schwarzschild solution is spherical, we can have the spherical spacetime metric \[g_{ij}=\left| {\begin{array} *{1}&{0}&{0}&{0}\\ {0}&{-1}&{0}&{0}\\ {0}&{0}&{-r^2}&{0}\\ {0}&{0}&{0}&{-r^2\sin^2\theta} \end{array}} \right|\] Then we can make some modifications by adding some new information \(U,V,W,X\) \[g_{ij}=\left| {\begin{array} *{U}&{0}&{0}&{0}\\ {0}&{-V}&{0}&{0}\\ {0}&{0}&{-Wr^2}&{0}\\ {0}&{0}&{0}&{-Xr^2\sin^2\theta} \end{array}} \right|\] Schwarzschild solution is non-angular. So we do not need to consider angular parts of the metric, we can just treating \(W=X=1\) to get \[g_{ij}=\left| {\begin{array} *{U}&{0}&{0}&{0}\\ {0}&{-V}&{0}&{0}\\ {0}&{0}&{-r^2}&{0}\\ {0}&{0}&{0}&{-r^2\sin^2\theta} \end{array}} \right|\] Schwarzschild solution is static. We only have 4 variables, we said it is non-angular so \(\theta, \phi\) are out of picture, Now it is static, so we do not need to consider time \(t\). If we take derivative with respect to \(t, \theta, \phi\) we only can get 0, vanishes. \(U,V\) are only functions of \(r\), \(U=U(r),V=V(r)\). \[g_{ij}=\left| {\begin{array} *{U(r)}&{0}&{0}&{0}\\ {0}&{-V(r)}&{0}&{0}\\ {0}&{0}&{-r^2}&{0}\\ {0}&{0}&{0}&{-r^2\sin^2\theta} \end{array}} \right|\] Schwarzschild solution is vacuum. This means there is no matter and energy in the space, the stress-energy tensor is all of 0s, \(T_{ij}=0\).
The 4 non-zero metric components are \[g_{00}=U(r),g_{11}=-V(r),g_{22}=-r^2,g_{33}=-r^2\sin^2\theta\] The metric inverse or the dual metric 4 non-zero components are \[g^{00}=\frac{1}{U(r)},g^{11}=-\frac{1}{V(r)},g^{22}=-\frac{1}{r^2},g^{33}=-\frac{1}{r^2\sin^2\theta}\]
library(Deriv)## Warning: package 'Deriv' was built under R version 4.0.5g=matrix(c("U(r)","0","0","0","0","-V(r)","0","0","0","0","-r^2","0","0","0","0","-r^2*sin(theta)^2"),nrow=4)
inv_g=matrix(c("1/U(r)","0","0","0","0","-1/V(r)","0","0","0","0","-1/r^2","0","0","0","0","-1/(r^2*sin(theta)^2)"),nrow=4)
g## [,1] [,2] [,3] [,4]
## [1,] "U(r)" "0" "0" "0"
## [2,] "0" "-V(r)" "0" "0"
## [3,] "0" "0" "-r^2" "0"
## [4,] "0" "0" "0" "-r^2*sin(theta)^2"inv_g## [,1] [,2] [,3] [,4]
## [1,] "1/U(r)" "0" "0" "0"
## [2,] "0" "-1/V(r)" "0" "0"
## [3,] "0" "0" "-1/r^2" "0"
## [4,] "0" "0" "0" "-1/(r^2*sin(theta)^2)"Christoffel symbols of second kind, use a previous function expand to input \(4\times4\) metric tensor. I wrote a Comma() function for comma derivative, to provide more rules of differentiation. \[\Gamma^{k}_{\text{ }\text{ }ij} = \frac{1}{2}g^{k\alpha} \left(\frac{\partial g_{i\alpha}}{\partial x^{j}}+ \frac{\partial g_{\alpha j}}{\partial x^{i}} - \frac{\partial g_{ij}}{\partial x^{\alpha}}\right)\]
#comma derivative
Comma<-function(x,y){
#print(x)
if(x=="U(r)" & y=="r") return("U1(r)")
if(x=="-V(r)" & y=="r") return("-V1(r)")
if(x=="1/U(r)" & y=="r") return("-U1(r)/U(r)^2")
if(x=="-1/V(r)" & y=="r") return("V1(r)/V(r)^2")
if(x=="U1(r)/(2 * U(r))" & y=="r") return("U2(r)/(2*U(r))-U1(r)^2/(2*U(r)^2)")
if(x=="V1(r)/(2 * V(r))" & y=="r") return("V2(r)/(2*V(r))-V1(r)^2/(2*V(r)^2)")
if(x=="U1(r)/(2 * V(r))" & y=="r") return("U2(r)/(2*V(r))-U1(r)*V1(r)/(2*V(r)^2)")
if(x=="-(V1(r)/(2 * V(r)))" & y=="r") return("V2(r)/(2*V(r))-V1(r)^2/(2*V(r)^2)")
if(x=="-(r/V(r))" & y=="r") return("r*V1(r)/V(r)^2-1/V(r)")
if(x=="-(r * sin(theta)^2/V(r))" & y=="r") return("r * sin(theta)^2*V1(r)/V(r)^2-sin(theta)^2/V(r)")
return(Deriv(x,y))
}dot_prod<-function(a,b){
#paste0("(",c("1","2","3"),")")
a=paste0("(",a,")")
b=paste0("(",b,")")
res=Simplify(paste(paste(a,b,sep="*"),collapse="+") )
return(paste0("(",res,")"))
}
Gamma2nd<-function(k,i,j){
para=c("t","r","theta","phi")
g1=inv_g[k,]
g2=sapply(g[i,],function(m) Comma(m,para[j]))
g3=sapply(g[j,],function(m) Comma(m,para[i]))
g4=sapply(para,function(m) Comma(g[i,j],m))
g5=sapply(paste(paste(g2,g3,sep="+"),g4,sep="-"),Simplify)
res=Simplify(paste0(dot_prod(g1,g5),"/2"))
return(res)
}
data=expand.grid(1:4,1:4,1:4)
test=apply(data,1,function(x) Gamma2nd(x[3],x[2],x[1]))
Gamma=array(test,dim=c(4,4,4))
Gamma## , , 1
##
## [,1] [,2] [,3] [,4]
## [1,] "0" "U1(r)/(2 * U(r))" "0" "0"
## [2,] "U1(r)/(2 * U(r))" "0" "0" "0"
## [3,] "0" "0" "0" "0"
## [4,] "0" "0" "0" "0"
##
## , , 2
##
## [,1] [,2] [,3]
## [1,] "U1(r)/(2 * V(r))" "0" "0"
## [2,] "0" "V1(r)/(2 * V(r))" "0"
## [3,] "0" "0" "-(r/V(r))"
## [4,] "0" "0" "0"
## [,4]
## [1,] "0"
## [2,] "0"
## [3,] "0"
## [4,] "-(r * sin(theta)^2/V(r))"
##
## , , 3
##
## [,1] [,2] [,3] [,4]
## [1,] "0" "0" "0" "0"
## [2,] "0" "0" "1/r" "0"
## [3,] "0" "1/r" "0" "0"
## [4,] "0" "0" "0" "-(cos(theta) * sin(theta))"
##
## , , 4
##
## [,1] [,2] [,3] [,4]
## [1,] "0" "0" "0" "0"
## [2,] "0" "0" "0" "1/r"
## [3,] "0" "0" "0" "cos(theta)/sin(theta)"
## [4,] "0" "1/r" "cos(theta)/sin(theta)" "0"Use previous Ricci tensor function and generalize it to \(4\times4\).
Ricci<-function(i,j){
para=c("t","r","theta","phi")
term1=paste(Comma(Gamma[i,1,1],para[j]),Comma(Gamma[i,2,2],para[j]),Comma(Gamma[i,3,3],para[j]),Comma(Gamma[i,4,4],para[j]),sep="+")
term1=Simplify(term1)
term2=paste(Comma(Gamma[i,j,1],para[1]),Comma(Gamma[i,j,2],para[2]),Comma(Gamma[i,j,3],para[3]),Comma(Gamma[i,j,4],para[4]),sep="+")
term2=Simplify(paste0("-(",term2,")"))
mydata=expand.grid(1:4,1:4)
g3=apply(mydata,1,function(x) Gamma[i,x[1],x[2]] )
g4=apply(mydata,1,function(x) Gamma[x[2],j,x[1]] )
term3=Simplify(dot_prod(g3,g4))
g5=apply(mydata,1,function(x) Gamma[i,j,x[1]] )
g6=apply(mydata,1,function(x) Gamma[x[1],x[2],x[2]] )
term4=dot_prod(g5,g6)
term4=paste0("-(",term4,")")
res=paste(term1,term2,term3,term4,sep="+")
return(Simplify(res))
}
data=expand.grid(1:4,1:4)
Ricci=array(apply(data,1,function(x) Ricci(x[1],x[2])),dim = c(4, 4))
Ricci## [,1]
## [1,] "((0.25 * (U1(r)/U(r)) + 0.25 * (V1(r)/V(r)) - 1/r) * U1(r) - U2(r)/2)/V(r)"
## [2,] "0"
## [3,] "0"
## [4,] "0"
## [,2]
## [1,] "0"
## [2,] "(U2(r)/2 - 0.25 * (U1(r)^2/U(r)))/U(r) - (1/r + U1(r)/(4 * U(r))) * V1(r)/V(r)"
## [3,] "0"
## [4,] "0"
## [,3]
## [1,] "0"
## [2,] "0"
## [3,] "(1 + r * (U1(r)/(2 * U(r)) - 0.5 * (V1(r)/V(r))))/V(r) - 1"
## [4,] "0"
## [,4]
## [1,] "0"
## [2,] "0"
## [3,] "0"
## [4,] "((1 + r * (U1(r)/(2 * U(r)) - 0.5 * (V1(r)/V(r))))/V(r) - 1) * sin(theta)^2"Ricci_scalar=sapply(paste0("(",inv_g,")*(",Ricci,")"),Simplify)
Ricci_scalar=Simplify(paste(Ricci_scalar,collapse="+"))
Ricci_scalar## [1] "(((0.25 * (V1(r)/V(r)) + 0.5 * (U1(r)/U(r)) - 1/r) * U1(r) - U2(r))/U(r) + (1/r + U1(r)/(4 * U(r))) * V1(r)/V(r))/V(r) - 2 * (((1 + r * (U1(r)/(2 * U(r)) - 0.5 * (V1(r)/V(r))))/V(r) - 1)/r^2)"uvtoab<-function(x){
x=gsub("U(r)", "U", x,fixed=TRUE)
x=gsub("V(r)", "V", x,fixed=TRUE)
x=gsub("U1(r)", "a", x,fixed=TRUE)
x=gsub("U2(r)", "b", x,fixed=TRUE)
x=gsub("V1(r)", "m", x,fixed=TRUE)
x=gsub("V2(r)", "n", x,fixed=TRUE)
x=gsub("sin(theta)^2", "z", x,fixed=TRUE)
x=gsub(" ", "", x,fixed=TRUE)
return(x)
}
abtouv<-function(x){
x=gsub("U", "U(r)", x,fixed=TRUE)
x=gsub("V", "V(r)", x,fixed=TRUE)
x=gsub("a", "U1(r)", x,fixed=TRUE)
x=gsub("b", "U2(r)", x,fixed=TRUE)
x=gsub("m", "V1(r)", x,fixed=TRUE)
x=gsub("n", "V2(r)", x,fixed=TRUE)
x=gsub("z", "sin(theta)^2", x,fixed=TRUE)
x=gsub(" ", "", x,fixed=TRUE)
return(x)
}
U=1.2
V=1.4
a=1.5
b=1.6
m=1.7
n=1.8
z=1.9
r=1.1
eval_str<-function (x){eval(parse(text=uvtoab(x)))}I computed the Ricci tensor and Ricci scalar values numerically, and checked with published results got matching.
eval(parse(text=uvtoab(Ricci[1,1])))## [1] -0.885378e1="-U2(r)/2/V(r)+U1(r)/4*V1(r)/V(r)^2 +U1(r)^2/4/U(r)/V(r)-U1(r)/r/V(r)"
eval(parse(text=uvtoab(e1)))## [1] -0.885378eval(parse(text=uvtoab(Ricci[2,2])))## [1] -1.207319e2="U2(r)/2/U(r)-U1(r)^2/4/U(r)^2 -U1(r)/2/U(r)*V1(r)/2/V(r)-V1(r)/V(r)/r"
eval(parse(text=uvtoab(e2)))## [1] -1.207319eval(parse(text=uvtoab(Ricci[3,3])))## [1] -0.2716837e3="r*U1(r)/2/U(r)/V(r)+1/V(r)-r*V1(r)/2/V(r)^2 - 1"
eval(parse(text=uvtoab(e3)))## [1] -0.2716837eval(parse(text=uvtoab(Ricci[4,4])))## [1] -0.516199e4=paste0("(",e3,")*sin(theta)^2")
eval(parse(text=uvtoab(e4)))## [1] -0.516199eval(parse(text=uvtoab(Ricci_scalar)))## [1] 0.5736194eR="-U2(r)/U(r)/V(r)+U1(r)*V1(r)/2/U(r)/V(r)^2 +U1(r)^2/2/U(r)^2/V(r)-2*U1(r)/r/U(r)/V(r)+2*V1(r)/r/V(r)^2 +2/r^2*(1 -1/V(r))"
eval(parse(text=uvtoab(eR)))## [1] 0.5736194Now we can have the 3 famous equations, the 4th equation is a duplicate. To help simplifying, equation 1 has common factor \(-1/U\), equation 2 has common factor \(1/V\), equation 3 has common factor \(2V/r\), can be canceled. Here are some numerical checks with published results. As you see, I don’t have anything to make simplification. \[R_{11}-\frac{1}{2}g_{11}R=0\\ R_{22}-\frac{1}{2}g_{22}R=0\\ R_{33}-\frac{1}{2}g_{33}R=0\]
eq1=paste0(Ricci[1,1],"-1/2*(",g[1,1],")*(",Ricci_scalar,")")
eq2=paste0(Ricci[2,2],"-1/2*(",g[2,2],")*(",Ricci_scalar,")")
eq3=paste0(Ricci[3,3],"-1/2*(",g[3,3],")*(",Ricci_scalar,")")
eval_str(eq1)/(-U)## [1] 1.024625eq1.="1/r*V1(r)/V(r)^2+1/r^2*(1-1/V(r))"
eval(parse(text=uvtoab(eq1.)))## [1] 1.024625eval_str(eq2)/V## [1] -0.5755608eq2.="-U1(r)/r/U(r)/V(r)+1/r^2 *(1 -1/V(r))"
eval(parse(text=uvtoab(eq2.)))## [1] -0.5755608eval_str(eq3)*2*V/r## [1] 0.1918155eq3.="-U1(r)/U(r)+V1(r)/V(r)-r*U2(r)/U(r)+r*U1(r)*V1(r)/2/U(r)/V(r)+r*U1(r)^2/2/U(r)^2"
eval(parse(text=uvtoab(eq3.)))## [1] 0.1918155