Coordinate approach to General Relativity

Whether to perturb

Clear[Perturb];
Perturb=False

Algorithm for simplifying anything; this one involves taking all equations to linear order in the small parameter ϵ

Clear[S];
If[Perturb&&(ϵ∈Integers)=!=True,
S[x__]:=Simplify[Normal[Series[Together[x, Trig->True],{ϵ,0,1}]],TimeConstraint->Infinity],
S[x__]:=Simplify[Together[x, Trig->True],TimeConstraint->Infinity],
];

Choose number of spacetime dimensions

(* Number of time and space dimensions, usually 1 and 3 *)
Clear[TimeDimension,SpaceDimension,Dimension];
TimeDimension=1;
SpaceDimension=3;

(* Number of spacetime dimensions *)
Clear[Dimension];
Dimension=TimeDimension+SpaceDimension

Choose coordinates

(* Choose time coordinate *)
Clear[TimeOption,Time,Conformal];
TimeOption=Time

(* Choose spatial coordinates *)
Clear[SpaceOption,Cartesian,Polar,ComovingPolar,Cylindrical];
SpaceOption=ComovingPolar

Clear[FirstIndex,LastIndex];
(* FirstIndex=0 for a time dimension, or FirstIndex=1 for only spatial dimensions *)
FirstIndex=If[TimeDimension>0,0,1];
LastIndex=SpaceDimension;

Clear[t,r,x,y,z,η,θ,ϕ];
Clear[Coordinate];
Which[
TimeOption===Time,
Coordinate[0]=t;
,TimeOption===Conformal,
Coordinate[0]=η;
];
Which[
SpaceOption===Cartesian,
Coordinate[1]=x;
Coordinate[2]=y;
If[LastIndex≥3,Coordinate[3]=z];
If[LastIndex≥4,Coordinate[4]=w];
If[LastIndex≥5,Coordinate[5]=v];
If[LastIndex≥6,Coordinate[6]=u];
,SpaceOption===Polar,
Coordinate[1]=r;
Coordinate[2]=θ;
If[LastIndex≥3,Coordinate[3]=ϕ];
If[LastIndex≥4,Coordinate[4]=χ];
If[LastIndex≥5,Coordinate[5]=ξ];
If[LastIndex≥6,Coordinate[6]=ψ];
If[LastIndex≥7,Coordinate[7]=ζ];
,SpaceOption===ComovingPolar,
Coordinate[1]=x;
Coordinate[2]=θ;
Coordinate[3]=ϕ;,
SpaceOption===Cylindrical,
Coordinate[1]=r;
Coordinate[2]=ϕ;
Coordinate[3]=z;
];
Table[Coordinate[i],{i,FirstIndex,LastIndex}]

(* Cunning way to define a function f[(x^μ)] of the coordinates (x^μ), whatever they may be. *)Clear[FunctionCoordinates,FunctionCoordinatesFirstIndex,FunctionCoordinatesLastIndex];
FunctionCoordinatesFirstIndex=1;
FunctionCoordinatesLastIndex=LastIndex;
FunctionCoordinates[func_]:=(func[Table[Coordinate[i],{i,FunctionCoordinatesFirstIndex,FunctionCoordinatesLastIndex}]]/.{x__}->x)

Choose metric sign (-1 for +---, or +1 for -+++)

Clear[MetricSign];
MetricSign=1

Imaginary

(* Define our own imaginary 𝕀 because Mathematica's I does not allow constructs such as I -> -I *)
Clear[𝕀];
Unprotect[Power];
Power[𝕀,n_]:=(((-1))^(Quotient[n,2]))(𝕀^(Mod[n,2]))/;n∈Integers&&Quotient[n,2]≠0;
Protect[Power];

Simplification rules for printing out quantities

(* Define simplification rules for printing out quantities *)
Clear[SPrint];
Which[
SpaceOption===Cartesian,SPrint[f__]:=Simplify[f//.{(x^2)+(y^2)+(z^2)->(r^2),((((x^2)+(y^2)+(z^2)))^(n_Integer))/;EvenQ[n]->(r^n)}/.{FunctionCoordinates[func_]->func}/.{𝕀->I},Assumptions->r≥0];
,True,
SPrint[f__]:=Simplify[f/.{FunctionCoordinates[func_]->func}/.{𝕀->I}];
]

Define the metric tensor as a function of the coordinates

Clear[MetricOption,Newton,Rindler,UniformlyRotating,FriedmannLemaitreRobertsonWalker,Schwarzschild,ReissnerNordstrom,GullstrandPainleve,StaticSpherical,BoyerLindquist,Doran,ADM,DupinCyclide];
MetricOption=FriedmannLemaitreRobertsonWalker

Clear[Metric];
Do[Metric[i,j]=0,{i,FirstIndex,LastIndex},{j,FirstIndex,LastIndex}];

Which[
(* Newtonian (weak field, pressureless matter) metric *)
MetricOption===Newton,
Which[
TimeOption===Time&&SpaceOption===Cartesian,
Metric[0,0]=-1;
Do[Metric[i,i]=1,{i,1,3}];
,TimeOption===Time&&SpaceOption===Polar,
Metric[0,0]=-1;
Metric[1,1]=1;
Metric[2,2]=(r^2);
Metric[3,3]=(r^2)((Sin[θ])^2);
];

Clear[Φ,Ψ];
Clear[Φf,Ψf];
If[Perturb,
(*Ψf=FunctionCoordinates[Ψ];*)
Φf=Φ[Coordinate[1],Coordinate[2],Coordinate[3]];
(*Wf[i_Integer]:=FunctionCoordinates[W[i]];*)
Metric[0,0]*=1+MetricSign 2ϵ Φf;
Do[Metric[i,i]*=1-MetricSign 2ϵ Φf,{i,1,3}]
(*Do[Metric[0,i]=Metric[i,0]=-MetricSign ϵ Wf[i],{i,1,3}];*)
];

(* Rindler metric *)
,MetricOption===Rindler,
Which[
TimeOption===Time&&SpaceOption===Cartesian,
Metric[0,0]=-(x^2);
Do[Metric[i,i]=1,{i,1,3}];
];

(* Uniformly rotating metric *)
,MetricOption===UniformlyRotating,
Which[
TimeOption===Time&&SpaceOption===Cylindrical,
Metric[0,0]=-1+(r^2)(ω^2);
Metric[1,1]=Metric[3,3]=1;
Metric[0,2]=Metric[2,0]=(r^2)ω;
Metric[2,2]=(r^2);
];

(* Friedmann-Robertson-Walter metric *)
,MetricOption===FriedmannLemaitreRobertsonWalker,
Clear[a,κ];
Which[
TimeOption===Time&&SpaceOption===Cartesian,
Metric[0,0]=-1;
Do[Metric[i,i]=((a[t])^2),{i,1,3}];

,TimeOption===Time&&SpaceOption===ComovingPolar,
Metric[0,0]=-1;
Metric[1,1]=((a[t])^2)/(1-κ (x^2));
Metric[2,2]=((a[t])^2)(x^2);
Metric[3,3]=((a[t])^2)(x^2)((Sin[θ])^2);

,TimeOption===Conformal&&SpaceOption===Cartesian,
Metric[0,0]=-((a[η])^2);
Do[Metric[i,i]=((a[η])^2),{i,1,3}];

,TimeOption===Conformal&&SpaceOption===ComovingPolar,
(*a[η]=a;*)
Metric[0,0]=-((a[η])^2);
Metric[1,1]=((a[η])^2)/(1-κ (x^2));
Metric[2,2]=((a[η])^2)(x^2);
Metric[3,3]=((a[η])^2)(x^2)((Sin[θ])^2);
];

Clear[Φ,Ψ];
Clear[Φf,Ψf];
If[Perturb,
Ψf=FunctionCoordinates[Ψ];
Φf=FunctionCoordinates[Φ];
(*Wf[i_Integer]:=FunctionCoordinates[W[i]];*)
Metric[0,0]*=1+MetricSign 2ϵ Ψf;
Do[Metric[i,i]*=1-MetricSign 2ϵ Φf,{i,1,3}]
(*Do[Metric[0,i]=Metric[i,0]=-MetricSign ϵ Wf[i],{i,1,3}];*)
];

(* static spherical metric *)
,MetricOption===StaticSpherical,
Clear[A,B];
Which[
TimeOption===Time&&SpaceOption===Polar,
Metric[0,0]=-A[r];
Metric[1,1]=B[r];
Metric[2,2]=(r^2);
Metric[3,3]=(r^2)((Sin[θ])^2);
];

(* Schwarzschild metric *)
,MetricOption===Schwarzschild,
Clear[M];
Which[
TimeOption===Time&&SpaceOption===Polar,
Δ[r_]:=1-((2M)/r);
Metric[0,0]=-Δ[r];
Metric[1,1]=(1/(Δ[r]));
Metric[2,2]=(r^2);
Metric[3,3]=(r^2)((Sin[θ])^2);
];

(* Reissner-Nordström metric *)
,MetricOption===ReissnerNordstrom,
Clear[M,Q];
Which[
TimeOption===Time&&SpaceOption===Polar,
Δ[r_]:=1-((2M)/r)+((Q^2)/(r^2));
Metric[0,0]=-Δ[r];
Metric[1,1]=(1/(Δ[r]));
Metric[2,2]=(r^2);
Metric[3,3]=(r^2)((Sin[θ])^2);
];

(* Gullstrand-Painlevé metric *)
,MetricOption===GullstrandPainleve,
Clear[β,M,Q];
Q=0;
β[r_]:=-(Sqrt[((2M)/r)-((Q^2)/(r^2))]);
Which[
TimeOption===Time&&SpaceOption===Polar,
Metric[0,0]=-(1-((β[r])^2));
Metric[0,1]=Metric[1,0]=β[r];
Metric[1,1]=1;
Metric[2,2]=(r^2);
Metric[3,3]=(r^2)((Sin[θ])^2);
];

(* Boyler-Lindquist metric *)
,MetricOption===BoyerLindquist,
Clear[a,M,Q];
Which[
TimeOption===Time&&SpaceOption===Polar,
Clear[β,M,Q,R,ρ];
(*Q=0;*)
β[r_]:=-((Sqrt[2M r-(Q^2)])/(R[r]));
R[r]=(Sqrt[(r^2)+(a^2)]);
ρ[r,θ]=(Sqrt[(r^2)+(a^2)((Cos[θ])^2)]);
Metric[0,0]=-1+((((R[r])^2)((β[r])^2))/((ρ[r,θ])^2));
Metric[1,1]=(((ρ[r,θ])^2)/(((R[r])^2)(1-((β[r])^2))));
Metric[2,2]=((ρ[r,θ])^2);
Metric[0,3]=Metric[3,0]=-((a ((Sin[θ])^2)((R[r])^2)((β[r])^2))/((ρ[r,θ])^2));
Metric[3,3]=((((Sin[θ])^2)((R[r])^2)(((ρ[r,θ])^2)+(a^2)((Sin[θ])^2)((β[r])^2)))/((ρ[r,θ])^2));
];

(* Doran metric *)
,MetricOption===Doran,
Clear[β,M,Q,R,ρ];
(*Q=0;*)
β[r_]:=-((Sqrt[2M r-(Q^2)])/(R[r]));
Which[
TimeOption===Time&&SpaceOption===Polar,
Clear[M,Q,R,ρ];
Q=0;
R[r]=(Sqrt[(r^2)+(a^2)]);
ρ[r,θ]=(Sqrt[(r^2)+(a^2)((Cos[θ])^2)]);
Metric[0,0]=-((((ρ[r,θ])^2)-2M r+(Q^2))/((ρ[r,θ])^2));
Metric[0,1]=Metric[1,0]=-β[r];
Metric[1,1]=(((ρ[r,θ])^2)/((R[r])^2));
Metric[2,2]=((ρ[r,θ])^2);
Metric[0,3]=Metric[3,0]=-((a ((Sin[θ])^2)((R[r])^2)((β[r])^2))/((ρ[r,θ])^2));
Metric[1,3]=Metric[3,1]=a ((Sin[θ])^2)β[r];
Metric[3,3]=((((Sin[θ])^2)((R[r])^2)(((ρ[r,θ])^2)+(a^2)((Sin[θ])^2)((β[r])^2)))/((ρ[r,θ])^2));
];

(* ADM metric *)
,MetricOption===ADM,
Which[
TimeOption===Time&&SpaceOption===Cartesian,
Clear[α,αf,β,βf,g,gf,Metric];
αf:=FunctionCoordinates[α];
βf[i_]:=FunctionCoordinates[β[i]];
g[i_Integer,j_Integer]:=g[j,i]/;i>j;
gf[i_,j_]:=FunctionCoordinates[g[i,j]];
Metric[0,0]:=Module[{i,j},-(α^2)+(∑_(i=1)^LastIndex)(∑_(j=1)^LastIndex)Metric[i,j]βf[i]βf[j]];
Metric[i_Integer,0]:=Metric[0,i]/;i>0;
Metric[0,i_Integer]:=Module[{j},-(∑_(j=1)^LastIndex)Metric[i,j]βf[j]]/;1≤i≤LastIndex;
Metric[i_Integer,j_Integer]:=gf[i,j]/;1≤i≤LastIndex&&1≤j≤LastIndex;
];

(* Dupin Cyclide *)
,MetricOption===DupinCyclide,
Which[
TimeOption===Time&&SpaceOption===Cartesian,
Metric[0,0]=-1;
Metric[1,1]=(((((2(1+2(y^2)))/(1+(x^2)+(y^2)))))^2);
Metric[2,2]=(((((2(1+2(x^2)))/(1+(x^2)+(y^2)))))^2);
];
];

(* change sign of metric to have desired signature *)If[MetricSign===-1,Do[Metric[i,j]=-Metric[i,j],{i,FirstIndex,LastIndex},{j,FirstIndex,LastIndex}]]

Clear[MetricMatrix];
MetricMatrix=Table[S[Metric[i,j]],{i,FirstIndex,LastIndex},{j,FirstIndex,LastIndex}];

SPrint[MetricMatrix]//MatrixForm

(* metric determinant *)
Clear[DeterminantMetricMatrix];
DeterminantMetricMatrix=S[Det[MetricMatrix]];
SPrint[DeterminantMetricMatrix]

Clear[SqrtAbsDeterminantMetricMatrix];
SqrtAbsDeterminantMetricMatrix=PowerExpand[S[Sqrt[Sqrt[(DeterminantMetricMatrix^2)]]]];
SPrint[SqrtAbsDeterminantMetricMatrix]

(* Totally antisymmetric (Levi-Civita) tensor (Subscript[ε, ijkl]) *)
Clear[TotallyAntisymmetric];
(* sign of (Subscript[ε, 0123]) in locally inertial frame *)
SignTotallyAntisymmetric=-1;
TotallyAntisymmetric[i_Integer,j_Integer,k_Integer,l_Integer]:=SignTotallyAntisymmetric SqrtAbsDeterminantMetricMatrix Signature[{i,j,k,l}]

Compute the inverse metric

Clear[InverseMetricMatrix];
InverseMetricMatrix=S[Inverse[Table[Metric[i,j],{i,FirstIndex,LastIndex},{j,FirstIndex,LastIndex}]]];

Clear[InverseMetric];
InverseMetric[i_,j_]:=InverseMetric[i,j]=InverseMetricMatrix[[i+1-FirstIndex,j+1-FirstIndex]];

Table[SPrint[InverseMetric[i,j]],{i,FirstIndex,LastIndex},{j,FirstIndex,LastIndex}]//MatrixForm

Compute the connection coefficients (all indices down)

Clear[CovariantConnection];
CovariantConnection[i_Integer,j_Integer,k_Integer]:=CovariantConnection[i,j,k]=S[(1/2)((∂_(Coordinate[k]))Metric[i,j]+(∂_(Coordinate[j]))Metric[i,k]-(∂_(Coordinate[i]))Metric[j,k])];
Do[Print["CovariantConnection[",i,",",j,",",k,"] = ",SPrint[CovariantConnection[i,j,k]]],{i,FirstIndex,LastIndex},{j,FirstIndex,LastIndex},{k,j,LastIndex}]

The connection coefficients (first index up)

Clear[Connection];
Connection[i_Integer,j_Integer,k_Integer]:=Connection[i,j,k]=Module[{l},S[(∑_(l=FirstIndex)^LastIndex)InverseMetric[i,l]CovariantConnection[l,j,k]]]

Do[Print["Connection[",i,",",j,",",k,"] = ",SPrint[Connection[i,j,k]]],{i,FirstIndex,LastIndex},{j,FirstIndex,LastIndex},{k,j,LastIndex}]

Compute the components of the completely covariant Riemann tensor (all indices down).  This uses all the symmetries.

(* Covariant Riemann tensor *)
Clear[CovariantRiemannf];
CovariantRiemannf[i_Integer,j_Integer,k_Integer,l_Integer]:=Module[{m},S[(∂_(Coordinate[i]))CovariantConnection[k,l,j]-(∂_(Coordinate[j]))CovariantConnection[k,l,i]+(∑_(m=FirstIndex)^LastIndex)(Connection[m,k,j]CovariantConnection[m,l,i]-Connection[m,k,i]CovariantConnection[m,l,j])]];

(* Covariant Riemann tensor invoking symmetries *)
Clear[CovariantRiemann];
CovariantRiemann[i_,i_,k_,l_]:=0;
CovariantRiemann[i_,j_,k_,k_]:=0;
CovariantRiemann[i_Integer,j_Integer,k_,l_]:=-CovariantRiemann[j,i,k,l]/;i>j;
CovariantRiemann[i_,j_,k_Integer,l_Integer]:=-CovariantRiemann[i,j,l,k]/;k>l;
CovariantRiemann[i_Integer,j_Integer,k_Integer,l_Integer]:=CovariantRiemann[k,l,i,j]/;i<j&&k<l&&(i>k||i==k&&j>l);
CovariantRiemann[i_Integer,j_Integer,k_Integer,l_Integer]:=-CovariantRiemann[i,l,j,k]-CovariantRiemann[i,k,l,j]/;i<k<l<j;
CovariantRiemann[i_Integer,j_Integer,k_Integer,l_Integer]:=CovariantRiemann[i,j,k,l]=CovariantRiemannf[i,j,k,l];

Do[If[i<j&&k<l&&(i<k||(i==k&&j≤l)),Print["CovariantRiemann[",i,",",j,",",k,",",l,"] = ",SPrint[CovariantRiemann[i,j,k,l]]]],{i,FirstIndex,LastIndex},{j,i,LastIndex},{k,FirstIndex,LastIndex},{l,k,LastIndex}]

The Riemann tensor with first two indices up

Clear[UpRiemann];
UpRiemann[i_,i_,k_,l_]:=0;
UpRiemann[i_,j_,k_,k_]:=0;
UpRiemann[i_Integer,j_Integer,k_,l_]:=-UpRiemann[j,i,k,l]/;i>j;
UpRiemann[i_,j_,k_Integer,l_Integer]:=-UpRiemann[i,j,l,k]/;k>l;
UpRiemann[i_Integer,j_Integer,k_Integer,l_Integer]:=UpRiemann[i,j,k,l]=Module[{m,n},S[(∑_(m=FirstIndex)^LastIndex)InverseMetric[i,m](∑_(n=FirstIndex)^LastIndex)InverseMetric[j,n]CovariantRiemann[m,n,k,l]]]

(*Do[If[i<j&&k<l,temp=SPrint[UpRiemann[i,j,k,l]];If[temp=!=0,Print["UpRiemann[",i,",",j,",",k,",",l,"] = ",temp]]],{i,FirstIndex,LastIndex},{j,i,LastIndex},{k,FirstIndex,LastIndex},{l,k,LastIndex}]*)

The Ricci tensor

Clear[CovariantRicci];
CovariantRicci[i_Integer,j_Integer]:=CovariantRicci[j,i]/;i>j;
CovariantRicci[i_Integer,j_Integer]:=CovariantRicci[i,j]=Module[{k,l},S[(∑_(k=FirstIndex)^LastIndex)(∑_(l=FirstIndex)^LastIndex)InverseMetric[k,l]CovariantRiemann[k,i,l,j]]];

Do[Print["CovariantRicci[",i,",",j,"] = ",SPrint[CovariantRicci[i,j]]],{i,FirstIndex,LastIndex},{j,i,LastIndex}]

The Ricci scalar

Clear[RicciScalar];
RicciScalar=Module[{i,j},S[(∑_(i=FirstIndex)^LastIndex)(∑_(j=FirstIndex)^LastIndex)InverseMetric[i,j]CovariantRicci[i,j]]];
SPrint[RicciScalar]

Compute the covariant Einstein tensor (both indices down)

Clear[CovariantEinstein];
CovariantEinstein[i_Integer,j_Integer]:=CovariantEinstein[j,i]/;i>j;
CovariantEinstein[i_Integer,j_Integer]:=CovariantEinstein[i,j]=S[CovariantRicci[i,j]-(1/2)Metric[i,j]RicciScalar]

Module[{temp},Do[temp=SPrint[CovariantEinstein[i,j]];If[temp=!=0,Print["CovariantEinstein[",i,",",j,"] = ",temp]],{i,FirstIndex,LastIndex},{j,i,LastIndex}]]

Compute the Einstein tensor (first index up, the second down)

Clear[Einstein];
Einstein[i_Integer,j_Integer]:=Einstein[i,j]=S[(∑_(k=FirstIndex)^LastIndex)InverseMetric[i,k]CovariantEinstein[k,j]]

Module[{temp},Do[temp=SPrint[Einstein[i,j]];If[temp=!=0,Print["Einstein[",i,",",j,"] = ",temp]],{i,FirstIndex,LastIndex},{j,FirstIndex,LastIndex}]]

The zeroth order (unperturbed) part of the Einstein tensor

If[Perturb,
Module[{temp},Do[temp=SPrint[S[Einstein[i,j]/.ϵ->0]];If[temp=!=0,Print["(Einstein^((0)))[",i,",",j,"] = ",temp]],{i,FirstIndex,LastIndex},{j,FirstIndex,LastIndex}]]]

The first order (perturbed) part of the Einstein tensor

If[Perturb,
Module[{temp},Do[temp=SPrint[Einstein[i,j]];If[temp=!=0,temp=SeriesCoefficient[Series[Einstein[i,j],{ϵ,0,1}],1]];If[temp=!=0,Print["(Einstein^((1)))[",i,",",j,"] = ",temp]],{i,FirstIndex,LastIndex},{j,FirstIndex,LastIndex}]]]

Compute the completely covariant Weyl tensor (all indices down)

Clear[CovariantRicciCombo];
CovariantRicciCombo[k_,k_,m_,n_]:=0;
CovariantRicciCombo[k_,l_,m_,m_]:=0;
CovariantRicciCombo[k_Integer,l_Integer,m_,n_]:=-CovariantRicciCombo[l,k,m,n]/;k>l;
CovariantRicciCombo[k_,l_,m_Integer,n_Integer]:=-CovariantRicciCombo[k,l,n,m]/;m>n;
CovariantRicciCombo[k_Integer,l_Integer,m_Integer,n_Integer]:=CovariantRicciCombo[m,n,k,l]/;k<l&&m<n&&(k>m||(k==m&&l>n));
CovariantRicciCombo[k_Integer,l_Integer,m_Integer,n_Integer]:=-CovariantRicciCombo[k,n,l,m]-CovariantRicciCombo[k,m,n,l]/;k<m<n<l;
CovariantRicciCombo[k_Integer,l_Integer,m_Integer,n_Integer]:=CovariantRicciCombo[k,l,m,n]=S[-(1/2) (Metric[k,m]CovariantRicci[l,n]-Metric[k,n]CovariantRicci[l,m]+Metric[l,n] CovariantRicci[k,m]-Metric[l,m]CovariantRicci[k,n])+(1/6)(Metric[k,m]Metric[l,n]-Metric[k,n]Metric[l,m])RicciScalar];

Clear[CovariantWeyl];
CovariantWeyl[k_,k_,m_,n_]:=0;
CovariantWeyl[k_,l_,m_,m_]:=0;
CovariantWeyl[k_Integer,l_Integer,m_,n_]:=-CovariantWeyl[l,k,m,n]/;k>l;
CovariantWeyl[k_,l_,m_Integer,n_Integer]:=-CovariantWeyl[k,l,n,m]/;m>n;
CovariantWeyl[k_Integer,l_Integer,m_Integer,n_Integer]:=CovariantWeyl[m,n,k,l]/;k<l&&m<n&&(k>m||(k==m&&l>n));
CovariantWeyl[k_Integer,l_Integer,m_Integer,n_Integer]:=-CovariantWeyl[k,n,l,m]-CovariantWeyl[k,m,n,l]/;k<m<n<l;
CovariantWeyl[k_Integer,l_Integer,m_Integer,n_Integer]:=CovariantWeyl[k,l,m,n]=S[CovariantRiemann[k,l,m,n]+CovariantRicciCombo[k,l,m,n]];

Do[If[i<j&&k<l&&(i<k||(i==k&&j≤l)),Print["CovariantWeyl[",i,",",j,",",k,",",l,"] = ",SPrint[CovariantWeyl[i,j,k,l]]]],{i,FirstIndex,LastIndex},{j,i,LastIndex},{k,FirstIndex,LastIndex},{l,k,LastIndex}]

The Weyl tensor with first two indices up

Clear[UpWeyl];
UpWeyl[i_,i_,k_,l_]:=0;
UpWeyl[i_,j_,k_,k_]:=0;
UpWeyl[i_Integer,j_Integer,k_,l_]:=-UpWeyl[j,i,k,l]/;i>j;
UpWeyl[i_,j_,k_Integer,l_Integer]:=-UpWeyl[i,j,l,k]/;k>l;
UpWeyl[i_Integer,j_Integer,k_Integer,l_Integer]:=UpWeyl[i,j,k,l]=Module[{m,n},S[(∑_(m=FirstIndex)^LastIndex)InverseMetric[i,m](∑_(n=FirstIndex)^LastIndex)InverseMetric[j,n]CovariantWeyl[m,n,k,l]]]

Do[If[i<j&&k<l,temp=SPrint[UpWeyl[i,j,k,l]];If[temp=!=0,Print["UpWeyl[",i,",",j,",",k,",",l,"] = ",temp]]],{i,FirstIndex,LastIndex},{j,i,LastIndex},{k,FirstIndex,LastIndex},{l,k,LastIndex}]

Define complexification projection operator

(* (1/2)(((Subscript[δ, i])^k)((Subscript[δ, j])^l) + (I/2)((Subscript[ε, ij])^kl)) *)
Clear[ComplexificationProjectionOperator];
ComplexificationProjectionOperator[i_Integer,j_Integer,k_Integer,l_Integer]:=ComplexificationProjectionOperator[i,j,k,l]=Module[{m,n},S[(1/2)(KroneckerDelta[i,k]KroneckerDelta[j,l]+(I/2)(∑_(m=FirstIndex)^LastIndex)InverseMetric[k,m](∑_(n=FirstIndex)^LastIndex)InverseMetric[l,n]TotallyAntisymmetric[i,j,m,n])]]

(*Do[temp=SPrint[ComplexificationProjectionOperator[i,j,k,l]];If[temp=!=0,Print["ComplexificationProjectionOperator[",i,",",j,",",k,",",l,"] = ",temp]],{i,FirstIndex,LastIndex},{j,FirstIndex,LastIndex},{k,FirstIndex,LastIndex},{l,FirstIndex,LastIndex}]*)

The complexified Weyl tensor

ComplexifiedCovariantWeyl[i_,i_,k_,l_]:=0; 
ComplexifiedCovariantWeyl[i_,j_,k_,k_]:=0;
ComplexifiedCovariantWeyl[i_Integer,j_Integer,k_,l_]:=-ComplexifiedCovariantWeyl[j,i,k,l]/;i>j;
ComplexifiedCovariantWeyl[i_,j_,k_Integer,l_Integer]:=-ComplexifiedCovariantWeyl[i,j,l,k]/;k>l;
ComplexifiedCovariantWeyl[i_Integer,j_Integer,k_Integer,l_Integer]:=ComplexifiedCovariantWeyl[k,l,i,j]/;i<j&&k<l&&(i>k||i==k&&j>l);ComplexifiedCovariantWeyl[i_,j_,k_,l_]:=Module[{m,n,p,q},S[(∑_(m=FirstIndex)^LastIndex)(∑_(n=FirstIndex)^LastIndex)ComplexificationProjectionOperator[i,j,m,n](∑_(p=FirstIndex)^LastIndex)(∑_(q=FirstIndex)^LastIndex)ComplexificationProjectionOperator[k,l,p,q]CovariantWeyl[m,n,p,q]]];

Do[If[i<j&&k<l&&(i<k||(i==k&&j≤l)),temp=SPrint[ComplexifiedCovariantWeyl[i,j,k,l]];If[temp=!=0,Print["ComplexifiedCovariantWeyl[",i,",",j,",",k,",",l,"] = ",temp]]],{i,FirstIndex,LastIndex},{j,i,LastIndex},{k,FirstIndex,LastIndex},{l,k,LastIndex}]

The complexified Weyl tensor with first two indices up

Clear[UpComplexifiedCovariantWeyl];
UpComplexifiedCovariantWeyl[i_,i_,k_,l_]:=0; 
UpComplexifiedCovariantWeyl[i_,j_,k_,k_]:=0;
UpComplexifiedCovariantWeyl[i_Integer,j_Integer,k_,l_]:=-UpComplexifiedCovariantWeyl[j,i,k,l]/;i>j;
UpComplexifiedCovariantWeyl[i_,j_,k_Integer,l_Integer]:=-UpComplexifiedCovariantWeyl[i,j,l,k]/;k>l;
UpComplexifiedCovariantWeyl[i_Integer,j_Integer,k_Integer,l_Integer]:=UpComplexifiedCovariantWeyl[i,j,k,l]=Module[{m,n},S[(∑_(m=FirstIndex)^LastIndex)InverseMetric[i,m](∑_(n=FirstIndex)^LastIndex)InverseMetric[j,n]ComplexifiedCovariantWeyl[m,n,k,l]]]

(*Do[If[i<j&&k<l,temp=SPrint[UpComplexifiedCovariantWeyl[i,j,k,l]];If[temp=!=0,Print["UpComplexifiedCovariantWeyl[",i,",",j,",",k,",",l,"] = ",temp]]],{i,FirstIndex,LastIndex},{j,i,LastIndex},{k,FirstIndex,LastIndex},{l,k,LastIndex}]*)