<<xAct`xPert`

DefManifold[M,4,{a,b,c,d,e,f,h,i,j,k,l, a1, b1, c1, d1, e1, f1, h1, i1, j1,k1,l1,a2,b2,c2,d2,e2,f2,h2,i2,j2,k2,l2}]
DefMetric[-1,g[-a,-b],CD]

DefMetricPerturbation[g,gpert,ϵ]
PrintAs[gpert]^="h"

Unprotect[IndexForm]
IndexForm[LI[x_]]:=ColorString[ToString[x],RGBColor[0,0,1]]
Protect[IndexForm]

DefTensor[n[a],M]
DefTensorPerturbation[Pertn[LI[order]],n[a],M,PrintAs→"δn"]

DefTensor[v[a],M]
DefTensorPerturbation[Pertn[LI[order]],v[a],M,PrintAs→"δv"]

DefTensor[Q[a,b],M]
DefTensorPerturbation[Pertn[LI[order]],Q[a,b],M,PrintAs→"δQ"]

(*AutomaticRules[n,MakeRule[{n[a] n[-a],-1}]]*)
(*AutomaticRules[n,MakeRule[{g[-a,-b] n[a] n[b],-1}]]*)
(*AutomaticRules[v,MakeRule[{Sqrt[Scalar[g[-a,-b] v[a] v[b]]],1}]]*)

L1=Perturbed[Sqrt[Scalar[g[-a,-b]v[a]v[b]]],1]/.Perturbation[v[_],1]-> 0//ExpandPerturbation

L2=Perturbation[WeylCD[-c,-d,-e,-f]n[d]n[f]Q[c,e],2]/.{Perturbation[n[_],1]-> 0,Perturbation[n[_],2]-> 0, Perturbation[Q[__],1]-> 0, Perturbation[Q[__],2]-> 0}//ExpandPerturbation;

L=L1 L2//Expand//ContractMetric//ToCanonical

DefMetric[{3,1,0},Mink[-a,-b],PD,{",","∂"},FlatMetric->True,PrintAs-> "η"]

(*var1=VarD[gpert[LI[1],-a,-b],CD][L]/.delta[-LI[1],LI[1]]→1//Simplify//Expand//ContractMetric//ToCanonicalP*)

gtoη=Join[MakeRule[{g[-a,-b],Mink[-a,-b]},MetricOn-> None,ContractMetrics-> False],MakeRule[{g[a,b],Mink[a,b]},MetricOn-> None,ContractMetrics-> Flase]]

L;
SeparateMetric[]/@%;
ChangeCurvature/@%;
ChangeCovD[#,CD,PD]&/@%;
ChristoffelToGradMetric/@%;
%//.gtoη;
Expand/@%;
ContractMetric/@%;
ToCanonical[#,UseMetricOnVBundle-> None]&/@%;
Lflat=Simplify/@%

(*rule1=Join[MakeRule[{Sqrt[Scalar[Mink[-a,-b]*v[a]*v[b]]],1},MetricOn-> All,ContractMetrics-> True]]
rule2=MakeRule[{n[a] n[-a],-1}]
rule3=MakeRule[{Mink[-a,-b] n[a] n[b],-1}]*)

(*AutomaticRules[v,MakeRule[{Scalar[Sqrt[Mink[-a,-b] v[a] v[b]]],1}]]*)

(*Lflat1=Lflat/.{ϵ->1, Scalar[Sqrt[Mink[-a,-b] v[a] v[b]]]-> 1, 1/(Scalar[Sqrt[Mink[-a,-b] v[a] v[b]]])-> 1} *)

Lflat/.ϵ-> 1

%/.(1/(Sqrt[Scalar[Mink[-a_, -b_]*v[a_]*v[b_]]]))-> 1
(*(Sqrt[Scalar[Mink[-a, -b]*v[a]*v[b]])-> 1}*)

%/.(Sqrt[Scalar[Mink[-a_, -b_]*v[a_]*v[b_]]])-> 1

Lflat1=%

(*FirstOrderOnly1= h[LI[n_],__]:>0/;n>1;*)

(*DefConstantSymbol[ρ]*)

(*Lflat1/. ρ*gpert[LI[n_],__]:>0/;n>1*)

(*%%/.FirstOrderOnly1*)

(*DefConstantSymbol[ρ]
Lflat2=Coefficient[Series[%/. ρ*gpert[args,__]:>σ*gpert[args],{σ,0,2}],σ,2]*)

DefConstantSymbol[σ]
(*Coefficient[Series[Lflat1/. gpert[args__]:>σ*gpert[args],{σ,0,2}],σ,2]*)
(*Lflat2=Coefficient[Series[Lflat1/. gpert[args__]:>σ*gpert[args],{σ,0,2}],σ,2]*)
Lflat2=Coefficient[Series[Lflat1/. gpert[args__]:>σ*gpert[args],{σ,0,2}],σ,2]/.PD[_][PD[_][gpert[LI[2],__]]]-> 0

(*DefTensor[ϕ[f,h],M]*)

(*Lflat2*)

(*Lflat2/.gpert[LI[1],-a_,-b_]:>  ϕ[α_,β_]*)

(*({{h, ({{1,  ,  }, { , a, b}})}})//InputForm*)

(*({{h, ({{1,  ,  ,   ,   }, { , h, j, ,i, ,k}})}})//InputForm*)

(*({{h, ({{1,  ,  ,   }, { , b, k, ,i}})}})//InputForm*)

(*gpert[LI[1], -a, -b]*)

(*gpert[LI[1], -a, -b]/:VarD[gpert[LI[1],f,h],PD][gpert[LI[1], -a_, -b_], rest_]:= 1/2(Mink[-a,-f]Mink[-b,-h]+Mink[-b,-f]Mink[-a,-h]) rest;*)

(*gpert[LI[1], -a, -b]/:VarD[gpert[LI[1],l,n],PD][gpert[LI[1], -a_, -b_], rest_]:= 1/2(Mink[-a,-l]Mink[-b,-n]+Mink[-b,-n]Mink[-a,-l]) rest;*)

(*PD[-i][gpert[LI[1], -d, -k]]/:VarD[gpert[LI[1],f,h],PD][PD[-i_][gpert[LI[1], -d_, -k_]], rest_]:= p[-i]1/2(Mink[-d,-f]Mink[-k,-h]+Mink[-k,-f]Mink[-d,-h]) rest;*)
(*PD[-i][gpert[LI[1], -d, -k]]/:VarD[gpert[LI[1],l,n],PD][PD[-i_][gpert[LI[1], -d_, -k_]], rest_]:= q[-i]1/2(Mink[-d,-l]Mink[-k,-n]+Mink[-k,-l]Mink[-d,-n]) rest;*)

(*PD[-c][PD[-e][gpert[LI[1], -k, -j]]]/:VarD[gpert[LI[1],f,h],PD][PD[-c_][PD[-e_][gpert[LI[1], -k_, -j_]]], rest_]:= p[-c]p[-e]1/2(Mink[-k,-f]Mink[-j,-h]+Mink[-j,-f]Mink[-k,-h]) rest; 
PD[-c][PD[-e][gpert[LI[1], -k, -j]]]/:VarD[gpert[LI[1],l,n],PD][PD[-c_][PD[-e_][gpert[LI[1], -k_, -j_]]], rest_]:= p[-c]p[-e]1/2(Mink[-k,-l]Mink[-j,-n]+Mink[-j,-l]Mink[-k,-n]) rest; *)

DefTensor[p[a],M]

DefTensor[q[a],M]

VarD[gpert[LI[1],a1,b1],PD][gpert[LI[1], -a_, -b_], rest_]:= 1/2(Mink[-a,-a1]Mink[-b,-b1]+Mink[-b,-a1]Mink[-a,-b1]) rest;

VarD[gpert[LI[1],a2,b2],PD][gpert[LI[1], -a_, -b_], rest_]:= 1/2(Mink[-a,-a2]Mink[-b,-b2]+Mink[-b,-a2]Mink[-a,-b2]) rest;

VarD[gpert[LI[1],a1,b1],PD][PD[-i_][gpert[LI[1], -d_, -k_]], rest_]:= p[-i]1/2(Mink[-d,-a1]Mink[-k,-b1]+Mink[-k,-a1]Mink[-d,-b1]) rest

VarD[gpert[LI[1],a2,b2],PD][PD[-i_][gpert[LI[1], -d_, -k_]], rest_]:= q[-i]1/2(Mink[-d,-a2]Mink[-k,-b2]+Mink[-k,-a2]Mink[-d,-b2]) rest;

VarD[gpert[LI[1],a1,b1],PD][PD[-c_][PD[-e_][gpert[LI[1], -f_, -j_]]], rest_]:= p[-c]p[-e]1/2(Mink[-f,-a1]Mink[-j,-b1]+Mink[-j,-a1]Mink[-f,-b1]) rest;

VarD[gpert[LI[1],a2,b2],PD][PD[-c_][PD[-e_][gpert[LI[1], -f_, -j_]]], rest_]:= q[-c]q[-e]1/2(Mink[-f,-a2]Mink[-j,-b2]+Mink[-j,-a2]Mink[-k,-b2]) rest;

res1=VarD[gpert[LI[1],a1,b1],PD][Lflat2] (*first variation*)
res2=VarD[gpert[LI[1],a2,b2],PD][res1] (*second variation*)/.delta[-LI[1],LI[1]]-> 1//ExpandPerturbation//ContractMetric//ToCanonical

(*DefTensor[ψ[a,b,-d],M]*)

(*DefTensor[η[a,b],M]*)

(*Lflat2/.InvMink[a,b]-> η[a,b]*)

(*rule1=n[a_] n[-a_]:> -1*)
(*rule2= Mink[-a_,-b_] n[a_] n[b_]:> -1*)

(*Lflat2/.rule2*)

(*%/.rule1*)

(*Coefficient[Series[Lflat2/. gpert[LI[1]__]:>σ*gpert[LI[1]],{σ,0,2}],σ,2]*)
(*(({{h, ({{1,  ,  }, { , a, b}})}}) ({{v, ({{a}, { }})}}) ({{v, ({{b}, { }})}})) ({{h, ({{2,  ,  ,   ,   }, { , a, b, ,h, ,i}})}})//InputForm*)
(*(1/12) ({{iη, ({{a, b}, { ,  }})}}) ({{iη, ({{h, i}, { ,  }})}}) ({{η, ({{ ,  }, {c, d}})}}) ({{η, ({{ ,  }, {e, f}})}}) ({{n, ({{c}, { }})}}) ({{n, ({{d}, { }})}}) ({{Q, ({{e, f}, { ,  }})}}) (({{h, ({{1,  ,  }, { , a, b}})}}) ({{v, ({{a}, { }})}}) ({{v, ({{b}, { }})}})) ({{h, ({{2,  ,  ,   ,   }, { , a, b, ,h, ,i}})}})*)
(*%/.PD[_][PD[_][gpert[LI[2],__]]]-> 0*)

(*(Sqrt[(({{η, ({{ ,  }, {a, b}})}}) ({{v, ({{a}, { }})}}) ({{v, ({{b}, { }})}}))])//InputForm*)

(*AutomaticRules[v,MakeRule[{Sqrt[Scalar[Mink[-a,-b] v[a] v[b]]],1}]]*)

(*%/.Sqrt[Scalar[Mink[-a,-b]*v[a]*v[b]]]-> 1*)