<<xAct`xTensor`
<<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, PrintAs -> "g", FlatMetric-> True]; *)
 DefMetric[1, g[-i, -j], PD,PrintAs -> "g", FlatMetric-> True]; 
InputForm[PD[-i][g[i, j]]]

(* 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[Al[],{M, t}, PrintAs-> "α" ];
DefTensor[Bt[],{M, t}, PrintAs -> "β" ];
DefTensor[Zt[],{M, t}, PrintAs -> "ζ" ];
(* perturbations in scalar field *)
DefTensor[Pii[],{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]];
(* bar n *)
DefTensor[Nb[],{M, t}, PrintAs->"N_0" ];
AutomaticRules[Nb,MakeRule[{PD[i][Nb[]], 0}, MetricOn->All, ContractMetrics->True]];
(* the scalar field and its potential *)
DefTensor[phi[],{M, t}, PrintAs->"ϕ" ];
AutomaticRules[phi,MakeRule[{PD[i][phi[]], 0}, MetricOn->All, ContractMetrics->True]];

DefTensor[V[],{M, t}, PrintAs->"V" ];
DefTensor[Vp[],{M, t}, PrintAs->"V_ϕ" ];
DefTensor[Vpp[],{M, t}, PrintAs->"V_{ϕϕ}" ];
AutomaticRules[V,MakeRule[{PD[i][V[]], 0}, MetricOn->All, ContractMetrics->True]];
AutomaticRules[Vp,MakeRule[{PD[i][Vp[]], 0}, MetricOn->All, ContractMetrics->True]];
AutomaticRules[Vpp,MakeRule[{PD[i][Vpp[]], 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 *)
al = de Al[];
zt = de Zt[];
bt =  de Bt[];
pi = 0 de Pii[];

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

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

(* define components of the 3D metric  *)

h[-i_,-j_] = A[]^2 E^(2 zt) g[-i,-j];
h[i_,-j_] = g[i,-j];
h[-i_,j_] = g[-i,j];
h[i_,j_] = A[]^(-2) E^(- 2 zt) g[i,j];

(* define the connection *)
gamma[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 = Simplify[ContractMetric[h[i,j] (PD[-k][gamma[k,-i,-j]] - PD[-i][gamma[k,-j,-k]])]];
InputForm[R3 = Simplification[R3 + ContractMetric[h[i,j] (ContractMetric[gamma[k,-i,-j]] gamma[k1,-k,-k1] - ContractMetric[gamma[k,-i,-k1]] gamma[k1,-j,-k]) ]] ]

(* expand up to the 2nd order in perturbation *)
InputForm[R3se = Normal[Series[R3, {de,0,2}]]  ]
InputForm[R31 = Coefficient[R3se,de]]
InputForm[R32 = Coefficient[R3se,de,2]]

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

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

(* expand up to the 2nd order in perturbation *)
InputForm[kijkijor2 = Simplify[Normal[Series[kijkij, {de,0,2}]] ]]

InputForm[kijkij0 = Simplify[Coefficient[kijkijor2, de,0] ]]
InputForm[kijkij1 = Simplify[Coefficient[kijkijor2, de,1] ]]
InputForm[kijkij2 = Simplify[Coefficient[kijkijor2, de,2]]]

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

(* expand up to the 2nd order in perturbation *)
InputForm[k2or2 = Simplify[Normal[Series[K2, {de,0,2}]] ]]

InputForm[K20 = Simplify[Coefficient[k2or2, de,0] ]]
InputForm[K21 = Simplify[Coefficient[k2or2, de,1] ]]
InputForm[K22 = Simplify[Coefficient[k2or2, de,2] ]]

(* Ricci scalar in 4D *)
(* Gauss-Codazzi relation *)
InputForm[R4 = Simplify[R3 + kijkij - K2]]

(* expand up to the 2nd order in perturbation *)
InputForm[R4or2 = Simplify[Normal[Series[R4, {de,0,2}]] ]]

(* compute \sqrt{-g} *)
InputForm[sqrtg = A[]^3 Simplify[Normal[Series[lbf E^( 3 zt), {de, 0, 2}]]]] 

(* expand up to the 2nd order in perturbation *)
InputForm[sqrtgOr2 = Simplify[Normal[Series[sqrtg, {de,0,2}]] ]]

InputForm[lagG = Simplify[sqrtgOr2 R4]]

InputForm[lagSe = Simplify[Normal[Series[lagG, {de, 0, 2}]]]]

(* 
the kinetic term of the scalar field  X = - g^{\mu\nu} \partial_\mu\phi \partial_\nu \phi,
where g_{\mu\nu} is the ADM metric
*)
InputForm[Xf = Simplify[ 
(1 / (2 lbf^2)) (ParamD[t][(phi[] + pi)]^2 - ContractMetric[2 h[i, j] n[-i] PD[-j][pi] ParamD[t][(phi[] + pi)]] - ContractMetric[(lbf^2 h[i, j] - h[i, k1] h[j, k2] n[-k1] n[-k2]) PD[-i][pi] PD[-j][pi]])
]]

(* expand up to the 2nd order perturbation *)
InputForm[Xfse = Simplify[Normal[Series[Xf, {de, 0, 2}]]]] 

InputForm[Xf0 = Simplify[Coefficient[Xfse, de, 0]]]
InputForm[Xf1 = Simplify[Coefficient[Xfse, de, 1]]]
InputForm[Xf1 = Simplify[Coefficient[Xfse, de, 2]]]

(* the potential of the scalar field expanded up to the 2nd perturbation *)
InputForm[pot = Simplify[V[] + Vp[] pi + Vpp[] pi pi / 2]]

(* The Lagrangian of the scalar field *)
InputForm[lagphi = Simplify[sqrtgOr2 (Xfse - pot)]]
	 	 
(* The total Lagrangian of the second order perturbations *)
InputForm[l2 = Simplification[ScreenDollarIndices[SameDummies[Coefficient[Normal[Series[lagG /2 + lagphi, {de, 0, 2}]], de, 2]]] 
/. Nb[] -> A[]
/. V[] -> (3 Hb[]^2 - ParamD[t][phi[]]^2 / 2) / A[]^2
]]

(* integration by parts *)
InputForm[cf = Coefficient[l2, PD[i][Bt[]] PD[-i][Zt[]] ]]
InputForm[l2 = Simplify[l2 - cf (PD[i][Bt[]] PD[-i][Zt[]] + Zt[] PD[-i][PD[i][Bt[]]]  )  ]]

InputForm[cf = Coefficient[l2, PD[-j][PD[-i][Bt[]]] PD[j][PD[i][Bt[]]] ]]
InputForm[l2 = Simplify[l2 - cf (PD[-j][PD[-i][Bt[]]] PD[j][PD[i][Bt[]]] - PD[-i][PD[i][Bt[]]] PD[-j][PD[j][Bt[]]])  ]]

(* Vary the action with respect to \alpha  *)
InputForm[eq1 = Simplify[D[l2, Al[]] 
/. Scalar[ParamD[t][phi[]] ]^2 -> ParamD[t][phi[]]^2 
/. PD[-i][PD[i][Bt[] ]] -> p2bt 
]]

(* Vary the action with respect to \beta_{, i}^{, i} *)

InputForm[eq2 = Simplify[D[l2, PD[-i][PD[i][Bt[] ]] ] ]]

InputForm[sol = Solve[{eq1 == 0, eq2 == 0}, {p2bt, Al[]}]]

InputForm[sol = Solve[{eq1 == 0, eq2 == 0}, {p2bt, Al[]}]]

InputForm[p2btSol = Simplify[p2bt /. sol[[1]] ]]
InputForm[alSol = Simplify[Al[] /. sol[[1]] ]]

InputForm[l2 = Simplify[NoScalar[l2
/. Al[] -> alSol
/. PD[-i][PD[i][Bt[] ]] -> p2btSol
]]]

(* Integration by part *)
InputForm[cf = Simplify[Coefficient[l2, Zt[] PD[-i][PD[i][Zt[]]] ]]]
InputForm[l2 = Simplify[l2 - cf (Zt[] PD[-i][PD[i][Zt[]]] + PD[-i][Zt[]] PD[i][Zt[]]) ]]

(* Integration by part *)
InputForm[cf = Simplify[Coefficient[l2, Zt[] ParamD[t][Zt[]] ]]]
InputForm[l2 = Simplify[l2 - cf Zt[] ParamD[t][Zt[]] - ParamD[t][A[]^4 cf]  Zt[]^2 / (2  A[]^4)
/. ParamD[t][A[]] -> A[] Hb[]
]]

(* Integration by part *)
InputForm[cf = Simplify[Coefficient[l2, ParamD[t][Zt[]] PD[-i][PD[i][Zt[]]] ]]]
InputForm[l2 = Simplify[l2 - cf ParamD[t][Zt[]] PD[-i][PD[i][Zt[]]] + ParamD[t][A[]^4 cf] PD[-i][Zt[]] PD[i][Zt[]] / (2 A[]^4)
/. ParamD[t][A[]] -> A[] Hb[]
/. ParamD[t][Hb[]] -> - ParamD[t][phi[]]^2 / 2 + Hb[]^2
]]
Exit