Higher-derivative terms in xPand

Quit
Clear["Global`*"]

Block[{Print},<<xAct`xPand`]//Quiet
$PrePrint=ScreenDollarIndices;
org[expr_]:=NoScalar@Collect[ContractMetric[expr],$PerturbationParameter,ToCanonical]
collect[expr_]:=NoScalar@Collect[expr,$PerturbationParameter,Identity]
(* Manifold *)
DefManifold[M,4,{μ,ν,α,β,λ,τ}]
(* Metric *)
DefMetric[-1,g[-μ,-ν],CD,{";","∇"},PrintAs->"g"]
SetSlicing[g,n,h,cd,{"|","(OverBar[D])"},"FLFlat"]
DefMetricFields[g,dg,h]
DefMatterFields[u,du,h]

$FirstOrderVectorPerturbations=False;$FirstOrderTensorPerturbations=True;$ConformalTime=False;
MyToxPand[expr_,gauge_,order_]:=ToxPand[expr,dg,u,du,h, gauge,order]

killsecondorder:=(bigexpr:(expr_[LI[x_],LI[y_],___]/;x>1)):>0
KillVectorsAndScalars=Flatten[{Outer[ϕh[LI[#1],LI[#2]]->0&,{1,2},{0,1,2}],
Outer[ψh[LI[#1],LI[#2]]->0&,{1,2},{0,1,2}]}]

(* Action *)
L0=Sqrt[-Detg[]] RicciScalarCD[];
EHaction=MyToxPand[L0,"SynchronousGauge",2]//org;
ExtractOrder[EHaction ,2]/.KillVectorsAndScalars/.killsecondorder//ToCanonical;
%//org;
Collect[%/(4 Sqrt[-(Detg[]*Scalar[ah[LI[0],LI[0]]]^6)])//org,{Eth[LI[1],LI[0],-α,-μ]*Eth[LI[1],LI[0],α,μ],Eth[LI[1],LI[0],α,μ]}]

L1=Sqrt[-Detg[]](RicciScalarCD[]RicciScalarCD[]); 
R2action=MyToxPand[L1,"SynchronousGauge",2]//org;
ExtractOrder[R2action ,2]/.KillVectorsAndScalars/.killsecondorder//ToCanonical;
%//org;
Collect[%/(4 Sqrt[-(Detg[]*Scalar[ah[LI[0],LI[0]]]^6)])//org,{Eth[LI[1],LI[0],-α,-μ]*Eth[LI[1],LI[0],α,μ],Eth[LI[1],LI[0],α,μ]}]

(* As you can see, this works fine for L0 and L1, however when considering 6-derivative terms ... *)

L2=Sqrt[-Detg[]] g[-μ,-ν]RicciScalarCD[]CD[μ][CD[ν][RicciScalarCD[]]];
RboxRaction=MyToxPand[L2,"SynchronousGauge",2]//org;
ExtractOrder[RboxRaction ,2]/.KillVectorsAndScalars/.killsecondorder//ToCanonical;
%//org;
Collect[%/(4 Sqrt[-(Detg[]*Scalar[ah[LI[0],LI[0]]]^6)])//org,{Eth[LI[1],LI[0],-α,-μ]*Eth[LI[1],LI[0],α,μ],Eth[LI[1],LI[0],α,μ]}]