<<xAct`xTensor`
<<xAct`TexAct`
<<xAct`xPert`

DefManifold[M, 3, {i, j, k}];
(* the metric g works as 3D Kronecker delta in this code and PD is the spatial-partial derivative *)
 DefMetric[1, g[-i, -j], PD, SymbolOfCovD->",",PrintAs -> "g", FlatMetric-> True]; 

(* define possible indices *)
GetIndicesOfVBundle[TangentM, 25];
IndicesOfVBundle[TangentM];

(* A time can be treated as a parameter in 3+1 formalism *)
DefParameter[t];

(* Perturbations *)
(* the metric perturbations *)
DefTensor[Psi[],{M, t}, PrintAs-> "ψ" ];
DefTensor[Phi[],{M, t}, PrintAs -> "ϕ" ];
DefTensor[Bt[],{M, t}, PrintAs -> "β" ];

(* Background quantities *)
(* the Hubble parameter *)
DefTensor[Hb[],{M, t}, PrintAs->"H" ];
AutomaticRules[Hb,MakeRule[{PD[i][Hb[]], 0}, MetricOn->All, ContractMetrics->True]];
(* the scale factor *)
DefTensor[A[],{M, t}, PrintAs->"a" ];
AutomaticRules[A,MakeRule[{PD[i][A[]], 0}, MetricOn->All, ContractMetrics->True]];

(* background n *)
DefTensor[Nb[],{M, t}, PrintAs->"N_0" ];
AutomaticRules[Nb,MakeRule[{PD[i][Nb[]], 0}, MetricOn->All, ContractMetrics->True]];

(* the 3D metric *)
DefTensor[h[-i,-j],M, Symmetric[{-i,-j}],PrintAs->"h"];
(* define constant "de" to count order of the perturbations *)
DefConstantSymbol[de, PrintAs -> "δ"];

(*  define short-hand notations *)
ps = de Psi[];
ph = de Phi[];
bt = 0 de Bt[];

(* define lapse function and shift vector *)
(* lbf = A[] (1 + ps); *)
lbf = (1 + ps);

DefTensor[n[-i],M, PrintAs->"N" ];
n[i_] := A[]^2 PD[i][bt];

(* define components of the 3D metric *)

h[-i_,-j_] = A[]^2 (1 - 2 ph) g[-i,-j];
h[i_,-j_] = g[i,-j];
h[-i_,j_] = g[-i,j];
h[i_,j_] = A[]^(-2) (1 + 2 ph) g[i,j];

(* Gamma^i_{jk} *)
gam[i_, j_, k_] := (h[i,k20] / 2) (PD[j][h[k,-k20]] + PD[k][h[j,-k20]] - PD[-k20][h[j,k]]);

(* the 3D Ricci scalar *)
R3 = Simplification[SameDummies[ContractMetric[h[i,j] (PD[-k][gam[k,-i,-j]] - PD[-i][gam[k,-j,-k]])]
+  (ContractMetric[h[i,j] gam[k,-i,-j]] gam[k1,-k,-k1] - ContractMetric[h[i,j] gam[k,-i,-k1]] gam[k1,-j,-k] )
]];

(* define the extrinsic curvature K_{ij} *)
kij[i_, j_] :=  (1 / lbf) ((Hb[] - ParamD[t][ph]) h[i,j] - PD[i][n[j]] + ContractMetric[gam[k,i,j] n[-k]  ]);

(* compute K_{ij} K^{ij} *)
kijkij = NoScalar[Simplification[ContractMetric[kij[-k2,-k3] h[k2,k4] h[k3,k5] kij[-k4,-k5]] ]];

(* expand up to the 1st order in perturbation *)
kijkij = Simplify[Normal[Series[kijkij, {de, 0, 1}]] ];
CoefficientList[kijkij, de]

(* compute K = h^{ij}K_{ij} *)
kii = Simplification[ContractMetric[kij[-k2,-k3] h[k2,k3] ] ];
(* compute K^2 *)
K2 = NoScalar[Simplification[ReplaceDummies[kii] kii ]];

(* expand up to the 1st order in perturbation *)
K2 = Simplify[Normal[Series[K2, {de, 0, 1}]]];
CoefficientList[K2, de]

(* compute K_{ij} K^{ij} - K^2 *)
kijk2 = Simplification[kijkij - K2];
R4 = Simplification[R3 + kijk2];

(* expand up to the 1st order in perturbation *)
CoefficientList[Simplify[Normal[Series[R4, {de, 0, 1}]]], de]

(* G_{\mu\nu}  n^\mu n^\nu *)
Gnn = NoScalar[Simplification[Normal[Series[
 R3 / 2  - kijk2 / 2
, {de, 0, 1}]]]
];

(* We get the component 0-0 of the Einstein equation *)
G00 = CoefficientList[Gnn, de];

(* G_{\mu\nu} n^\mu P^\nu_\alpha *)
Gpn[k1_] := Simplification[ContractMetric[Normal[Series[
-PD[k1][kii] 
+ h[k2, k3] PD[-k2][kij[-k3, k1]]
-  h[k2, k3] gam[k4, -k2, -k3] kij[-k4, k1]
-  h[k2, k3] gam[k4, -k2, k1] kij[-k4, -k3]
, {de, 0, 1}]] 
]];

(* We get the component 0-i of the perturbed Einstein equation *)
Gnp = Simplification[CoefficientList[Gpn[-i], de]
/. ParamD[t][PD[-i][Zt[]]] -> PD[-i][ParamD[t][Zt[]]]
];