xPert

High-order perturbation theory in General Relativity

David Brizuela, david.brizuela@ehu.eus
José M. Martín-García, jose@xact.es
Guillermo A. Mena Marugán, mena@iem.cfmac.csic.es

(c) 2005-2020, under GPL

http://www.xact.es/
http://groups.google.com/group/xact

October 2005 - May 2006: Version 0.5.0: First public version.
September 2007: Version 0.7.0: Notebook adapted for version 0.9.2 of xAct.
November 2007: Version 0.7.1: Notebook adapted for version 0.9.3 of xAct.
March 2008: Version 0.7.2: Notebook adapted for version 0.9.4 of xAct.
May-June 2008: Version 1.0.0. Substantial update for publication.

This package expands any finite-order perturbation of any of the GR-relevant curvature tensors in terms of metric perturbations.

The package xPert is distributed under the GNU GPL License, and runs on top of xTensor, a free package for fast manipulation of abstract tensor expressions in Mathematica.

DateList[]

xAct`xPert`$Version={"1.0.6",{2018,2,28}};

xAct`xPert`$xTensorVersionExpected={"1.1.4",{2020,2,16}};

1. Initialization

1.1. GPL

(* xPert, computer algebra for perturbation theory in General Relativity *)

(* Copyright (C) 2005-2020 David Brizuela, Jose M. Martin-Garcia and Guillermo A. Mena Marugan *)

(* This program is free software; you can redistribute it and/or
 modify it under the terms of the GNU General Public License as
 published by the Free Software Foundation; either version 2 of
 the License,or (at your option) any later version.

This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
 General Public License for more details.

You should have received a copy of the GNU General Public License
 along with this program; if not, write to the Free Software
 Foundation, Inc., 59 Temple Place-Suite 330, Boston, MA 02111-1307,
  USA. 
*)

1.2. Info package

(* :Title: xPert *)

(* :Author: David Brizuela, Jose M. Martin-Garcia and 
	Guillermo A. Mena Marugan *)

(* :Summary: Computer algebra for perturbation theory in General
    Relativity *)

(* :Brief Discussion:
   - Single-parameter metric perturbation theory
   - Inert-head Perturbation and explicit expansion formulas
     through ExpandPerturbation
   - Handles densities, derivatives and all index positions in
     tensors
   - Formula for arbitrary gauge change 
*)
  
(* :Context: xAct`xPert` *)

(* :Package Version: 1.0.6 *)

(* :Copyright: David Brizuela, Jose M. Martin-Garcia and 
	Guillermo A. Mena Marugan (2005-2020) *)

(* :History: see xPert.History *)

(* :Keywords: *)

(* :Source: xPert.nb *)

(* :Warning: ToCanonical can only act on expressions with the
    Perturbation head with option UseMetricOnVBundle->None *)

(* :Mathematica Version: 6.0 and later *)

(* :Limitations: 
	- Only one perturbation structure can be defined.
	- Only single-parameter metric perturbation theory *)

1.3. BeginPackage

Protect against multiple loading of the package:

With[{xAct`xPert`Private`xPertSymbols=DeleteCases[Join[Names["xAct`xPert`*"],Names["xAct`xPert`Private`*"]],"$Version"|"xAct`xPert`$Version"|"$xTensorVersionExpected"|"xAct`xPert`$xTensorVersionExpected"]},
Unprotect/@xAct`xPert`Private`xPertSymbols;
Clear/@xAct`xPert`Private`xPertSymbols;
]

Decide which is the last package being read:

If[Unevaluated[xAct`xCore`Private`$LastPackage]===xAct`xCore`Private`$LastPackage,xAct`xCore`Private`$LastPackage="xAct`xPert`"];

xAct`xCore`Private`$LastPackage

Explicit (not hidden) import of other packages. QUESTION: Do we really need ExpressionManipulation?

BeginPackage["xAct`xPert`",{"xAct`xTensor`","xAct`xPerm`","xAct`xCore`","xAct`ExpressionManipulation`"}]

Check version of xTensor. We simply compare dates:

If[Not@OrderedQ@Map[Last,{$xTensorVersionExpected,xAct`xTensor`$Version}],Throw@Message[General::versions,"xTensor",xAct`xTensor`$Version,$xTensorVersionExpected]]

Welcome message:

Print[xAct`xCore`Private`bars];
Print["Package xAct`xPert`  version ",$Version[[1]],", ",$Version[[2]]];
Print["CopyRight (C) 2005-2020, David Brizuela, Jose M. Martin-Garcia and Guillermo A. Mena Marugan, under the General Public License."];

We specify the context xAct`xPert` to avoid overriding the Disclaimer of the other packages. However we need to turn off the message General:shdw temporarily:

Off[General::shdw]
xAct`xPert`Disclaimer[]:=Print["These are points 11 and 12 of the General Public License:\n\nBECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW. EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM `AS IS´ WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU. SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING, REPAIR OR CORRECTION.\n\nIN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR REDISTRIBUTE THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGES."]
On[General::shdw]

If this is the last package show the GPL short disclaimer:

If[xAct`xCore`Private`$LastPackage==="xAct`xPert`",
Unset[xAct`xCore`Private`$LastPackage];
Print[xAct`xCore`Private`bars];
Print["These packages come with ABSOLUTELY NO WARRANTY; for details type Disclaimer[]. This is free software, and you are welcome to redistribute it under certain conditions. See the General Public License for details."];
Print[xAct`xCore`Private`bars]];

Note that symbols in the Global` context cannot be accessed while building the package:

$ContextPath

$Context

1.4. Set up

Screen all dollar indicess:

ReportSet[$PrePrint,ScreenDollarIndices];

Postfix format for derivatives:

ReportSet[$CovDFormat,"Postfix"];

In xPert we shall allow metric contractions everywhere:

ReportSetOption[ContractMetric,AllowUpperDerivatives->True];

ReportSetOption[MakeRule,MetricOn->All];

ReportSetOption[MakeRule,ContractMetrics->True];

We do not modify the the option MathLink of CanonicalPerm:

Options[CanonicalPerm]

$xpermQ

1.5. Usage messages

Define perturbation objects:

Perturbation::usage="Perturbation[expr] represents the perturbation of expr. Perturbation[expr, n] represents its n-th perturbation. Perturbation is an inert-head.";
Perturbed::usage="Perturbed[expr, n] returns the expansion of expr perturbed up to n-th order, using $PerturbationParameter.";
$PerturbationParameter::usage="$PerturbationParameter is a global variable containing the parameter used in the perturbative expansions.";
PerturbationParameter::usage="PerturbationParameter[metric] stores the perturbation parameter associated to the perturbation of the given metric.";
ExpandPerturbation::usage="ExpandPerturbation[expr] returns expr with all Perturbation expressions expanded in terms of metric perturbations. It uses fast pre-stored formulas for the general term of the perturbative expansion of the most common curvature tensors.";
DefTensorPerturbation::usage="DefTensorPerturbation[pert[LI[n], inds], tensor[inds], args] defines pert to be the perturbation of tensor, with the same index configuration. There is an internal call to DefTensor[pert[LI[n], inds], args].";
PerturbationOrder::usage="PerturbationOrder[expr] gives the perturbative order of expression, as given by the orders of its arguments. This function must be defined on those perturbations declared by the user.";
GaugeChange::usage="GaugeChange[pert, gen] returns the perturbative expression pert changed to a new gauge. The change of gauge is parametrized by the family of vector fields gen[LI[order], a]. Background objects are left untouched.";
GeneralPerturbation::usage="GeneralPerturbation[expr, n] returns ExpandPerturbation[Perturbation[expr, n]]. This is kept for backwards compatibility with older versions.";

Perturbation structure on already existing vbundle and metric:

DefMetricPerturbation::usage="DefMetricPerturbation[metric, pert, param] constructs a number of definitions for the perturbation of the given metric, which must exist already. A tensor pert[LI[order], -a, -b] is defined on the vbundle of the metric, and the parameter param will be defined and used to expand the perturbed expressions. DefMetricPerturbation[metric, pert] uses a default name for the parameter.";

Versions:

xAct`xPert`$Version::usage="$Version is a global variable giving the version of the package xPert in use.";
xAct`xPert`$xTensorVersionExpected::usage="$xTensorVersionExpected is a global variable giving the oldest possible version of the package xTensor which is required by the version of the package xPert in use.";

1.6. Begin private

Begin["xAct`xPert`Private`"]

Names["xAct`xPert`*"]

%//Length

There are 12 symbols in this package, but only 4 are nontrivial.

$ContextPath

$Context

2. Partitions

2.1. Sorted partitions

Generation of sorted partitions. The function group gives all partitions of length n-1 from a partition of length n. The function SortedPartitions simply computes all partitions starting from (1,1,...,1) applying group[ ] recursively until reaching length 0.

group[list_List]:=Table[Join[Take[list,i-1],{list[[i]]+list[[i+1]]},Drop[list,i+1]],{i,1,Length[list]-1}]

SortedPartitions[n_Integer?Positive]:=Flatten[FixedPointList[Union@Flatten[group/@#,1]&,{Table[1,{n}]}],1]

SortedPartitions[5]

Lengths 2^(i-1):

Table[{i,Length[SortedPartitions[i]]},{i,12}]

2.2. All partitions. Generalized Leibniz rule

addone[n_][{prev___Integer}]:=Table[{prev,i},{i,0,n-Total[{prev}]}];
onelevel[n_][list_List]:=Apply[Join,Map[addone[n],list]];
complete[n_][{prev___Integer}]:={prev,n-Total[{prev}]};
lastlevel[n_][list_List]:=Map[complete[n],list];

AllPartitions[0,n_Integer]:={};
AllPartitions[m_Integer,n_Integer]:=lastlevel[n][Nest[onelevel[n],{{}},m-1]];

AllPartitions[4,3]

LeibnizDistribute[operator_,n_,f_[args___],g_]:=Apply[g,Map[Multinomial@@# Inner[operator,{args},#,f]&,AllPartitions[Length[{args}],n]]];

Expand[(a+b+c+d+e+f)^20]//Length//AbsoluteTiming

AllPartitions[6,20]//Length//AbsoluteTiming

2.3. Folded partitions. Fa`a di Bruno rule

This are partitions {k1, ..., kn} of integer n obeying Sum[i*ki, {i,n}] == n .

FoldedPartitions[n_]:=Cases[Fold[step,fpart[{{},n},n],Reverse@Range[n-1]],{partition_,0}:>partition];

fpart[{{part___},0},d_]:={{{0,part},0}};
fpart[{{part___},m_},d_]:={{#,part},m-d #}&/@Range[0,m/d];
step[pairs_List,d_]:=Flatten[fpart[#,d]&/@pairs,1];

Folded partitions are currently being used by the GaugeChange formula.

FoldedPartitions[8]

Fa`a di Bruno formula:

faa[expr_][ki_,{0}]:=Scalar[expr]^ki;
faa[expr_][1,{i_}]:=Perturbation[expr,i]/i!;
faa[expr_][ki_,{i_}]:=(Scalar@Perturbation[expr,i]/i!)^ki;
FaadiBruno1[partition_,function_[expr_]]:=With[{n=partition.Range[Length[partition]],k=Plus@@partition},n!/Apply[Times,Factorial/@partition]Derivative[k][function][expr]Times@@MapIndexed[faa[expr],partition]];
FaadiBruno1[partition_,Power[expr_,exponent_]]:=With[{n=partition.Range[Length[partition]],k=Plus@@partition},n!/Apply[Times,Factorial/@partition]Derivative[k,0][Power][expr,exponent]Times@@MapIndexed[faa[expr],partition]];

FaadiBruno0[expr_,order_]:=Apply[Plus,FaadiBruno1[#,expr]&/@FoldedPartitions[order]]

A simpler and more general version based on the Dt builtin of Mathematica. Arguments of Power are wrapped with Scalar. Note however that this new version could be slower...

FaadiBruno[function_[args__],order_]:=Module[{vars=Table[Unique[x],{Length[{args}]}]},
Expand@Nest[Dt,function@@vars,order]/.Dt->Perturbation/.Power[x_Perturbation,n_]:>Scalar[x]^n/.g_?ScalarFunctionQ@@vars:>g[args]/.Inner[IndexRule,vars,{args},List]];

Tests:

FaadiBruno[function[x],7]-FaadiBruno0[function[x],7]//NoScalar

FaadiBruno[x^5,4]-FaadiBruno0[x^5,4]//NoScalar//Expand

3. Perturbation structure

Off[RuleDelayed::rhs]

3.1. DefTensorPerturbation

Given a tensor, we define its perturbation, with the same index structure and tensor properties:

makepattern[a_Symbol]:=a_Symbol;
makepattern[-a_Symbol]:=-a_Symbol;

Options[DefTensorPerturbation]=Options[DefTensor];
DefTensorPerturbation[pert_[LI[_],inds___],tensor_[inds___],args__]:=Module[{prot,pinds=makepattern/@{inds}},
DefTensor[pert[LI["order"],inds],args];
xAct`xTensor`Private`SymbolRelations[pert,Null,{tensor}];
prot=Unprotect[tensor,pert];
tensor/:Perturbation[tensor@@pinds,Optional[order_Integer]]:=pert[LI[order],inds];
pert[LI[0],indices___]:=tensor[indices];
pert/:PerturbationOrder[pert[LI[order_],indices___]]:=order;
pert/:Perturbation[pert[LI[order_],Sequence@@pinds],n_.]:=pert[LI[order+n],inds];
Protect[Evaluate[prot]];
];
SetNumberOfArguments[DefTensorPerturbation,{2,Infinity}];
Protect[DefTensorPerturbation];

3.2. DefMetricPerturbation

Jose: It is rather absurd that DefTensorPerturbation requires indices but DefMetricPerturbation does not.

?DefMetricPerturbation

Options[DefMetricPerturbation]={PrintAs->Identity};

DefMetricPerturbation[metric_,pert_Symbol,param_Symbol:Automatic,options___?OptionQ]:=Module[{vbundle,manifold,indices,covd,christoffel,riemann,ricci,ricciscalar,einstein,weyl,tfricci,kretschmann,pa},

{pa}={PrintAs}/.CheckOptions[options]/.Options[DefMetricPerturbation];

(* Check of metric *)
If[!MetricQ[metric],Throw@Message[DefMetricPerturbation::unknown,"metric",metric]];
(* Define parameter if it does not exist yet *)
If[!ParameterQ[param],
If[param===Automatic,
DefParameter[GiveSymbol[PerturbationParameter,metric],PrintAs:>GiveOutputString[PerturbationParameter,metric]],
DefParameter[param]
]
];

vbundle=VBundleOfMetric[metric];
manifold=BaseOfVBundle[vbundle];
indices=GetIndicesOfVBundle[vbundle,4];
DefTensorPerturbation[pert[LI["order"],-indices[[1]],-indices[[2]]],metric[-indices[[1]],-indices[[2]]],DependenciesOfTensor[metric],Symmetric[{2,3}],PrintAs->pa];
(* Store perturbation structure in global variables *)
covd=CovDOfMetric[metric];
christoffel=Christoffel[covd];
riemann=Riemann[covd];
ricci=Ricci[covd];
ricciscalar=RicciScalar[covd];
einstein=Einstein[covd];
weyl=Weyl[covd];
tfricci=TFRicci[covd];
kretschmann=Kretschmann[covd];
$PerturbationParameter=param;

(* Define perturbation definitions *)
Unprotect[Perturbation,GeneralPerturbation];
MakexTensions[DefMetricPerturbation,"Beginning",metric,pert,param];
With[{a=indices[[1]],b=indices[[2]],c=indices[[3]],d=indices[[4]]},
DefGenPertInvMetric[vbundle,metric[a,b],pert];
DefGenPertDet[vbundle,metric,pert];
DefGenPertInvMetric[vbundle,metric[a,b],pert];
DefGenPertChristoffel[vbundle,covd,pert,christoffel[a,-b,-c]];
DefGenPertRiemannFromChristoffel[vbundle,covd,christoffel,riemann[-a,-b,-c,d]];
DefGenPertRiemann[vbundle,covd,pert,riemann[-a,-b,-c,d]];
DefGenPertRicci[vbundle,covd,pert,riemann,ricci[-a,-b]];
DefGenPertRicciScalar[vbundle,metric,ricci,ricciscalar];
DefGenPertEinstein[vbundle,metric,ricci,ricciscalar,einstein[-a,-b]];
DefGenPertWeyl[weyl[-a,-b,-c,-d]];
DefGenPertTFRicci[tfricci[-a,-b]];
DefGenPertKretschmann[vbundle,metric,riemann,kretschmann[]];
];
MakexTensions[DefMetricPerturbation,"End",metric,pert,param];
Protect[Perturbation,GeneralPerturbation];

];
SetNumberOfArguments[DefMetricPerturbation,{2,3}];
Protect[DefMetricPerturbation];

Store the perturbation parameters of the various metrics:

PerturbationParameter[metric_?MetricQ]:=GiveSymbol[PerturbationParameter,metric];
SetNumberOfArguments[PerturbationParameter,1];
Protect[PerturbationParameter];

PrintAsCharacter[PerturbationParameter]="ε";

Metric perturbations:

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

IndexForm[LI[1]]

Additional object:

ThreePert[covd_,pert_][LI[order_],a_,b_,c_]:=1/2(covd[c][pert[LI[order],a,b]]+covd[b][pert[LI[order],a,c]]-covd[a][pert[LI[order],b,c]])

3.3. Perturbation

Main perturbation function, defined as inert-head. That means that the tensorial properties of Perturbation[expr] are exactly those of expr. The density definition must be given independently:

Options[DefInertHead]

(* Define linear inert head *)
DefInertHead[Perturbation,PrintAs->"△",DefInfo->False,ContractThrough->{delta},LinearQ->True];
(* Perturbation has two arguments, but the second defaults to 1 *)
Default[Perturbation,2]:=1;
Perturbation[expr_,0]:=expr;
Perturbation[expr_,1]:=Perturbation[expr];
(* Perturbation does not affect the density character *)
WeightOf[Perturbation[expr_,_.]]^:=WeightOf[expr];

It acts as a derivative and commutes with covariant partial derivatives or derivatives of scalars. As suggested by Teake, we block perturbations on derivatives of scalar densities:

Perturbation[_?ConstantQ,_.]:=0;
Perturbation[delta[a_,b_],_.]:=0;
Perturbation[PD[-a_][expr_],n_.]:=PD[-a][Perturbation[expr,n]];
Perturbation[cd_Symbol?CovDQ[-a_][expr_?ScalarQ],n_.]/;FreeQ[WeightOf@expr,WeightedWithBasis@cd]:=cd[-a][Perturbation[expr,n]];
Perturbation[ParamD[ps__][expr_],n_.]:=ParamD[ps][Perturbation[expr,n]];
Perturbation[OverDot[expr_],n_.]:=OverDot[Perturbation[expr,n]];
Perturbation[Scalar[expr_],n_.]:=Scalar[Perturbation[expr,n]];

It does not commute with covariant derivatives of general expressions:

Generalized Leibniz rule:

Perturbation[product_Times,n_.]:=LeibnizDistribute[Perturbation,n,product,Plus];

Extend to other algebra products:

Perturbation[product:_?ProductQ[__],n_.]:=LeibnizDistribute[Perturbation,n,product,Plus];

Perturbation/:Grade[Perturbation[expr_,n_.],prod_]:=Grade[expr,prod];

Perturbation of a scalar function:

Perturbation[expr:_?ScalarFunctionQ[__],order_.]:=FaadiBruno[expr,order];

Formatting (note the redundant check ?InertHeadQ to avoid problems with the new definition-ordering rules in Mathematica 6):

MakeBoxes[Perturbation?xAct`xTensor`Private`HeldInertHeadQ[expr_],StandardForm]:=xAct`xTensor`Private`interpretbox[Perturbation[expr],RowBox[{PrintAs[Perturbation],"[",MakeBoxes[expr,StandardForm],"]"}]];

MakeBoxes[Perturbation?xAct`xTensor`Private`HeldInertHeadQ[expr_,n_],StandardForm]:=xAct`xTensor`Private`interpretbox[Perturbation[expr,n],RowBox[{SuperscriptBox[PrintAs[Perturbation],ToString[n]],"[",MakeBoxes[expr,StandardForm],"]"}]];

Perturbations are joined together:

Perturbation[Perturbation[expr_,n_.],m_.]:=Perturbation[expr,n+m];

Finally:

SetNumberOfArguments[Perturbation,{1,2}];
Protect[Perturbation];

Examples:

Perturbation[a b c,10]

Perturbation[Sin[x],5]

Perturbation[Power[x,3],6]

Perturbation[Power[x,1/2],6]

General::canonpert="ToCanonical on Perturbation may incorrectly change index characters.";
xAct`xTensor`Private`Identify[Perturbation[expr_,n_.]]:=(If[(UseMetricOnVBundle/.xAct`xTensor`Private`$TCOptions)=!=None,Message[ToCanonical::canonpert]];xAct`xTensor`Private`addhead1A[Perturbation[expr,n],Perturbation[#,n]&,expr]);

3.4. Perturbed

This function gives the total perturbation up to a given order. It is constructed in terms of Perturbation, but could equally well be given in terms of GeneralPerturbation. Note the factor 1/n!

?Perturbed

Perturbed[expr_,order_Integer]:=Sum[$PerturbationParameter^n Perturbation[expr,n]/n!,{n,0,order}];
SetNumberOfArguments[Perturbed,2];
Protect[Perturbed];

3.5. ExpandPerturbation

This function knows nothing at the beginning. It must be trained little by little adding information to ExpandPerturbation1:

Options[ExpandPerturbation]:={SeparateMetric->True,OverDerivatives->True};
ExpandPerturbation[expr_,options___?OptionQ]:=expr/.expr1_Perturbation:>ExpandPerturbation1[expr1,options];
ExpandPerturbation[expr_,tensor_Symbol?xTensorQ,options___?OptionQ]:=expr/.expr1:Perturbation[tensor[___],n_.]:>ExpandPerturbation1[expr1,options];
ExpandPerturbation[expr_,list:{___Symbol?xTensorQ},options___?OptionQ]:=Fold[ExpandPerturbation[#1,#2,options]&,expr,list];
SetNumberOfArguments[ExpandPerturbation,{1,Infinity}];
Protect[ExpandPerturbation];

General definition, in use when none of the particular definitions work:

?ChangeCovD

ExpandPerturbation1[Perturbation[expr_,n_.],options___]:=Module[{sepmetric,od,tmp},
{sepmetric,od}={SeparateMetric,OverDerivatives}/.CheckOptions[options]/.Options[ExpandPerturbation];

(* Separate metrics *)
tmp=Perturbation[If[sepmetric,SeparateMetric[][expr],expr],n];
(* Derivatives *)
If[od,tmp=tmp/.expr1:HoldPattern[Perturbation[_Symbol?CovDQ[_][_]|LieD[_][_]|Bracket[_,_][_],_.]]:>ExpandPerturbationDer[expr1]];
(* Reexpand *)
If[tmp=!=Perturbation[expr,n],tmp=ExpandPerturbation[tmp,options]];

(* Return result *)
tmp
];
ExpandPerturbation1[expr_,options___]:=expr;

For backwards compatibility:

GeneralPerturbation[expr_,order_:1]:=ExpandPerturbation[Perturbation[expr,order]];

3.6. Perturbation of derivatives

ExpandPerturbationDer[Perturbation[cd_Symbol?CovDQ[-a_][expr_],n_.]]:=cd[-a][Perturbation[expr,n]]+With[{vbundles=VBundlesOfCovD[cd]},Plus@@Map[xAct`xTensor`Private`addChristoffel[expr,-a],xAct`xTensor`Private`selecton[Select[FindFreeIndices[expr],AIndexQ],vbundles]]/.expr1_ xAct`xTensor`Private`CHR[indices__]:>With[{chr=Christoffel[cd][indices]},Perturbation[expr1 chr,n]-chr Perturbation[expr1,n]]]+Module[{basis=WeightedWithBasis[cd],dummy=DummyIn@VBundleOfIndex[a],weight},weight=WeightOf[expr,basis];
If[weight===0,0,weight With[{chr=Christoffel[cd,xAct`xCoba`PDOfBasis[basis]][dummy,-a,-dummy]},-Perturbation[chr expr,n]+chr Perturbation[expr,n]]]];

ExpandPerturbationDer[Perturbation[LieD[v_][expr_],n_.]]:=Module[{lieD},
LeibnizDistribute[Perturbation,n,lieD[v,expr],Plus]/.lieD[V_,EXPR_]:>LieD[V][EXPR]
];

ExpandPerturbationDer[Perturbation[Bracket[v1_,v2_][a_Symbol],n_.],options___]:=Module[{bracket},
LeibnizDistribute[Perturbation,n,bracket[a][v1,v2],Plus]/.bracket[A_][V1_,V2_]:>Bracket[V1,V2][A]
];

3.7. PerturbationOrder

homorder[list_List]:=If[Length[Union[list]]===1,First[list],list];
PerturbationOrder[expr_Times]:=Apply[Plus,Map[PerturbationOrder,List@@expr]];
PerturbationOrder[expr_Plus]:=homorder@Map[PerturbationOrder,List@@expr];
PerturbationOrder[Power[expr_,n_]]:=n PerturbationOrder[expr];
PerturbationOrder[Perturbation[expr_,n_.]]:=n;
PerturbationOrder[_?FirstDerQ[expr_]]:=PerturbationOrder[expr];
PerturbationOrder[expr_]:=0;
SetNumberOfArguments[PerturbationOrder,1];
Protect[PerturbationOrder];

3.8. GaugeChange

Return perturbative expression in a different gauge. Background quantities are not changed. Nothing happens on inert heads.

GaugeChange[pert_Plus,gen_]:=GaugeChange[#,gen]&/@pert;
GaugeChange[pert_Times,gen_]:=GaugeChange[#,gen]&/@pert;
GaugeChange[func_Symbol?ScalarFunction[args__],gen_]:=Apply[func,GaugeChange[#,gen]&/@{args}];
GaugeChange[Scalar[pert_]^Optional[n_Integer],gen_]:=Scalar[xAct`xPert`GaugeChange[pert,gen]]^n;
GaugeChange[expr_,gen_]:=GaugeChange[expr,gen,PerturbationOrder[expr]];
GaugeChange[expr_,gen_,0]:=expr;
GaugeChange[expr_Perturbation,gen_,n_]:=Bruni[expr,gen,n];
GaugeChange[expr_,gen_,n_]:=With[{sexpr=SeparateMetric[][expr]},If[expr===sexpr,Bruni[expr,gen,n],GaugeChange[sexpr,gen]]];
Bruni[expr_,gen_,n_]:=expr+Expand@Sum[n!/(n-m)!
Apply[Plus,BruniTerm[Perturbation[expr,-m],gen,#]&/@FoldedPartitions[m]],{m,1,n}];
SetNumberOfArguments[GaugeChange,{2,3}];
Protect[GaugeChange];

Single term in Bruni's formula:

BruniTerm[expr_,gen_,partition_]:=With[{dummy=DummyIn@Last@SlotsOfTensor[gen],range=Range@Length@partition},1/Inner[Power,Factorial[range],partition,Times]/Times@@Factorial[partition]Fold[Nest[LieD[gen[LI[First@#2],dummy]],#1,Last@#2]&,expr,Reverse@Transpose[{range,partition}]]];

4. Perturbations of curvature tensors

4.1. Perturbation of the metric tensor. Not used now

DefPertMetric[metric_[-a_,-b_],pert_]:=(
Perturbation[metric[-a_Symbol,-b_Symbol],n_.]:=pert[LI[n],-a,-b];
Perturbation[pert[LI[n1_],-a_Symbol,-b_Symbol],n_.]:=pert[LI[n1+n],-a,-b]
);

4.2. Perturbation of the inverse metric tensor

In order to construct the GeneralPerturbation term we need to define

productPert[{p_},vbundle_,pert_][a_,b_]:=pert[LI[p],a,b]
productPert[partition_List,vbundle_,pert_][a_,b_]:=Module[{m=Length[partition],i,indices},
indices=Table[DummyIn[vbundle],{m-1}];pert[LI[partition[[1]]],a,indices[[1]]]Product[pert[LI[partition[[i+1]]],ChangeIndex[indices[[i]]],indices[[i+1]]],{i,1,m-2}]pert[LI[partition[[m]]],ChangeIndex[indices[[m-1]]],b]
];

DefGenPertInvMetric[vbundle_,metric_[a_,b_],pert_]:=(
ExpandPerturbation1[Perturbation[metric[a_Symbol,b_Symbol],order_.],options___]:=Plus@@(Map[(-1)^Length[#](Multinomial@@#)productPert[#,vbundle,pert][a,b]&,SortedPartitions[order]])
);

4.3. Perturbation of the determinant of the metric

Original code:

DefGenPertDet[vbundle_,metric_,pert_]:=With[{dim=DimOfVBundle[vbundle],metricepsilon=epsilon[metric],mdet=Determinant[metric][]},

ExpandPerturbation1[Perturbation[mdet,order_.]]:=With[{inds=Table[DummyIn[vbundle],{2dim}]},
With[{inds1=inds[[Range[dim]]],inds2=inds[[Range[dim+1,2dim]]]},ContractMetric[mdet/dim!SignDetOfMetric[metric]metricepsilon@@inds1 metricepsilon@@inds2 Perturbation[Times@@Apply[metric,Transpose[-{inds1,inds2}],{1}],order]]
]
];

ExpandPerturbation1[Perturbation[metricepsilon[superinds__?UpIndexQ],order_.]]:=metricepsilon[superinds]Sqrt[mdet]ExpandPerturbation[Perturbation[1/Sqrt[mdet],order]];
ExpandPerturbation1[Perturbation[metricepsilon[subinds__?DownIndexQ],order_.]]:=metricepsilon[subinds]/Sqrt[mdet]ExpandPerturbation[Perturbation[Sqrt[mdet],order]];

];

New code, contributed by Cyril Pitrou. Perturbation from the Log-Trace formula. Should work in any dimension, even abstract dimensions, but it does not use the identities when the order of perturbation is higher than the dimension:

productPertDet0[{p_},vbundle_,pert_]:=With[{ind=DummyIn[vbundle]},
pert[LI[p],-ind,ind]/p!
];

productPertDet0[partition_List,vbundle_,pert_]:=Module[{m=Length[partition],i,indices},
indices=Table[DummyIn[vbundle],{m}];
1/Times@@Map[#!&,partition]
 pert[LI[partition[[1]]],ChangeIndex[indices[[m]]],indices[[1]]]Product[pert[LI[partition[[i+1]]],ChangeIndex[indices[[i]]],indices[[i+1]]],{i,1,m-1}]
];

productPertDet1[partitionK_List,partitionL_List,vbundle_,pert_]:=With[{m=Length[partitionK]},
(-1)^m/m!/Times@@partitionK If[m≥2,Product[productPertDet0[Take[partitionL,{1+Plus@@Take[partitionK,i-1],Plus@@Take[partitionK,i]}],vbundle,pert],{i,2,m}],
1]
productPertDet0[Take[partitionL,{1,partitionK[[1]]}],vbundle,pert]
];

productPertDet2[partition_List,vbundle_,pert_]:=Plus@@Map[ productPertDet1[#,partition,vbundle,pert]&,SortedPartitions[Length[partition]]];

DefGenPertDet[vbundle_,metric_,pert_]:=With[{dim=DimOfVBundle[vbundle],metricepsilon=epsilon[metric],mdet=Determinant[metric][]},

ExpandPerturbation1[Perturbation[mdet,order_.]]:=With[{inds=Table[DummyIn[vbundle],{2dim}]},
With[{inds1=Take[inds,dim],inds2=Take[inds,-dim]},ContractMetric[mdet/dim!SignDetOfMetric[metric]metricepsilon@@inds1 metricepsilon@@inds2 Perturbation[Times@@Apply[metric,Transpose[-{inds1,inds2}],{1}],order]]
]
]/;IntegerQ[dim]&&order>dim;

ExpandPerturbation1[Perturbation[mdet,order_.]]:=
mdet order!Plus@@Map[(-1)^Length[#]productPertDet2[#,vbundle,pert]&,SortedPartitions[order]]/;!IntegerQ[dim]||order≤dim;

ExpandPerturbation1[Perturbation[metricepsilon[superinds__?UpIndexQ],order_.]]:=metricepsilon[superinds]Sqrt[mdet]ExpandPerturbation[Perturbation[1/Sqrt[mdet],order]];
ExpandPerturbation1[Perturbation[metricepsilon[subinds__?DownIndexQ],order_.]]:=metricepsilon[subinds]/Sqrt[mdet]ExpandPerturbation[Perturbation[Sqrt[mdet],order]];
];

4.4. Perturbation of the Christoffel tensor

General term:

productThreePert[{p_},vbundle_,covd_,pert_][a_,b_,c_]:=ThreePert[covd,pert][LI[p],a,b,c]
productThreePert[partition_List,vbundle_,covd_,pert_][a_,b_,c_]:=Module[{m=Length[partition],i,indices},
indices=Table[DummyIn[vbundle],{m-1}];pert[LI[partition[[1]]],a,indices[[1]]]Product[pert[LI[partition[[i+1]]],ChangeIndex[indices[[i]]],indices[[i+1]]],{i,1,m-2}]ThreePert[covd,pert][LI[partition[[m]]],ChangeIndex[indices[[m-1]]],b,c]]

DefGenPertChristoffel[vbundle_,covd_,pert_,christoffel_[a_,-b_,-c_]]:=(
ExpandPerturbation1[Perturbation[christoffel[a_Symbol,-b_Symbol,-c_Symbol],order_.]]:=Plus@@(Map[-(-1)^Length[#](Multinomial@@#)productThreePert[#,vbundle,covd,pert][a,-b,-c]&,SortedPartitions[order]])
);

4.5. Perturbation of the Riemann tensor

General term. There is first the expression in terms of the perturbations of the Christoffel tensor. We give it here, but do not include it in the package:

DefGenPertRiemannFromChristoffel[vbundle_,covd_,christoffel_,riemann_[-a_,-b_,-c_,d_]]:=(
GeneralPerturbation0[riemann[-a_Symbol,-b_Symbol,-c_Symbol,d_Symbol],order_Integer]:=Module[{dummy=DummyIn[vbundle]},2$RiemannSign Antisymmetrize[covd[-b][Perturbation[christoffel[d,-a,-c],order]]-Sum[Binomial[order,k]Perturbation[christoffel[d,-dummy,-a],k]Perturbation[christoffel[dummy,-b,-c],order-k],{k,1,order-1}],{-a,-b}]]);

Then we have an expression given directly in terms of the derivatives of the metric:

productRiemann1[{p_},vbundle_,covd_,pert_][a_,b_,c_,d_]:=covd[a][ThreePert[covd,pert][LI[p],d,c,b]]
productRiemann1[partition_List,vbundle_,covd_,pert_][a_,b_,c_,d_]:=Module[{m=Length[partition],i,indices},
indices=Table[DummyIn[vbundle],{m-1}];pert[LI[partition[[1]]],d,indices[[1]]]Product[pert[LI[partition[[i+1]]],ChangeIndex[indices[[i]]],indices[[i+1]]],{i,1,m-2}]covd[a][ThreePert[covd,pert][LI[partition[[m]]],ChangeIndex[indices[[m-1]]],c,b]]]

productRiemann2[partition_List,s_,vbundle_,covd_,pert_][a_,b_,c_,d_]:=Module[{m=Length[partition],i,indices},
indices=Table[DummyIn[vbundle],{s}];ThreePert[covd,pert][LI[partition[[s]]],indices[[s]],d,a]Product[pert[LI[partition[[i]]],ChangeIndex[indices[[i+1]]],indices[[i]]],{i,s-1,2,-1}]ThreePert[covd,pert][LI[partition[[1]]],ChangeIndex[indices[[2]]],b,c]]/;s==Length[partition]
productRiemann2[partition_List,s_,vbundle_,covd_,pert_][a_,b_,c_,d_]:=Module[{m=Length[partition],i,indices},
indices=Table[DummyIn[vbundle],{m}];pert[LI[partition[[m]]],d,indices[[m]]]Product[pert[LI[partition[[i]]],ChangeIndex[indices[[i+1]]],indices[[i]]],{i,m-1,s+1,-1}]ThreePert[covd,pert][LI[partition[[s]]],indices[[s]],ChangeIndex[indices[[s+1]]],a]Product[pert[LI[partition[[i]]],ChangeIndex[indices[[i+1]]],indices[[i]]],{i,s-1,2,-1}]ThreePert[covd,pert][LI[partition[[1]]],ChangeIndex[indices[[2]]],b,c]]/;2≤s<Length[partition]

genpertriemann[vbundle_,covd_,pert_,IndexList[-a_,-b_,-c_,d_],order_]:=2$RiemannSign Antisymmetrize[Plus@@(Map[(-1)^Length[#](Multinomial@@#)(productRiemann1[#,vbundle,covd,pert][-a,-b,-c,d]+Sum[productRiemann2[#,s,vbundle,covd,pert][-a,-b,-c,d],{s,2,Length[#]}])&,SortedPartitions[order]]),{-a,-b}];

DefGenPertRiemann[vbundle_,covd_,pert_,riemann_[-a_,-b_,-c_,d_]]:=(
ExpandPerturbation1[Perturbation[riemann[-a_Symbol,-b_Symbol,-c_Symbol,d_Symbol],order_.]]:=genpertriemann[vbundle,covd,pert,IndexList[-a,-b,-c,d],order];ExpandPerturbation1[Perturbation[riemann[-a_Symbol,-b_Symbol,-c_Symbol,-d_Symbol],order_.]]:=With[{dummy=DummyAs[c],metric=xAct`xTensor`Private`FirstMetricOfVBundle[VBundleOfIndex[d],True]},ExpandPerturbation@Perturbation[riemann[-a,-b,-c,dummy]metric[-dummy,-d],order]]
);

4.6. Perturbation of the Ricci tensor

We can get the GeneralPerturbation term directly from that of the Riemann tensor:

DefGenPertRicci[vbundle_,covd_,pert_,riemann_,ricci_[-a_,-b_]]:=(ExpandPerturbation1[Perturbation[ricci[-a_Symbol,-b_Symbol],order_.]]:=Module[{dummy=DummyIn[vbundle]},$RicciSign genpertriemann[vbundle,covd,pert,IndexList[-a,-dummy,-b,dummy],order]]);

4.7. Perturbation of the RicciScalar

The perturbation of the RicciScalar is constructed from that of Ricci and the inverse metric:

DefGenPertRicciScalar[vbundle_,metric_,ricci_,ricciscalar_]:=(ExpandPerturbation1[Perturbation[ricciscalar[],order_.]]:=Module[{dummy1=DummyIn[vbundle],dummy2=DummyIn[vbundle]},Sum[Binomial[order,k]ExpandPerturbation1@Perturbation[ricci[-dummy1,-dummy2],k]ExpandPerturbation1@Perturbation[metric[dummy1,dummy2],order-k],{k,0,order}]]);

4.8. Perturbation of the Einstein tensor

DefGenPertEinstein[vbundle_,metric_,ricci_,ricciscalar_,einstein_[-a_,-b_]]:=(ExpandPerturbation1[Perturbation[einstein[-a_Symbol,-b_Symbol],order_.]]:=ExpandPerturbation1@Perturbation[ricci[-a,-b],order]-1/2Sum[Binomial[order,k]ExpandPerturbation1@Perturbation[metric[-a,-b],k]ExpandPerturbation1@Perturbation[ricciscalar[],order-k],{k,0,order}]);

4.9. Perturbations of the Weyl tensor

DefGenPertWeyl[weyl_[-a_,-b_,-c_,-d_]]:=(ExpandPerturbation1[Perturbation[weyl[-a_Symbol,-b_Symbol,-c_Symbol,-d_Symbol],order_.]]:=ExpandPerturbation[Perturbation[WeylToRiemann[weyl[-a,-b,-c,-d]],order]]);

An expanded version of this formula was contributed by Cyril Pitrou (IAP, Paris) in Dec 2006:

DefGenPertWeyl[vbundle_,dimension_,metric_,ricci_,ricciscalar_,riemann_,weyl_[-a_,-b_,-c_,-d_]]:=(ExpandPerturbation1[Perturbation[weyl[-a_Symbol,-b_Symbol,-c_Symbol,-d_Symbol],order_.]]:= 
Module[{dummy=DummyIn[vbundle]},
Sum[Binomial[order,k]ExpandPerturbation1@Perturbation[metric[-d,-dummy],k]ExpandPerturbation1@Perturbation[riemann[-a,-b,-c,dummy],order-k],{k,0,order}]
-(Sum[Binomial[order,k]ExpandPerturbation1@Perturbation[metric[-a,-c],k]ExpandPerturbation1@Perturbation[ricci[-b,-d],order-k],{k,0,order}]
+Sum[Binomial[order,k]ExpandPerturbation1@Perturbation[metric[-b,-d],k]ExpandPerturbation1@Perturbation[ricci[-a,-c],order-k],{k,0,order}]
-Sum[Binomial[order,k]ExpandPerturbation1@Perturbation[metric[-a,-d],k]ExpandPerturbation1@Perturbation[ricci[-b,-c],order-k],{k,0,order}]
-Sum[Binomial[order,k]ExpandPerturbation1@Perturbation[metric[-b,-c],k]ExpandPerturbation1@Perturbation[ricci[-a,-d],order-k],{k,0,order}])/(dimension-2)
+(Sum[Binomial[order,k]ExpandPerturbation1@Perturbation[ricciscalar[],k](Sum[Binomial[order-k,m]ExpandPerturbation1@Perturbation[metric[-a,-c],m]ExpandPerturbation1@Perturbation[metric[-b,-d],order-k-m],{m,0,order-k}]),{k,0,order}]
-(Sum[Binomial[order,k]ExpandPerturbation1@Perturbation[ricciscalar[],k](Sum[Binomial[order-k,m]ExpandPerturbation1@Perturbation[metric[-a,-d],m]ExpandPerturbation1@Perturbation[metric[-b,-c],order-k-m],{m,0,order-k}]),{k,0,order}])
)/(dimension-1)/(dimension-2)
]);

4.10. Perturbations of the TFRicci tensor

DefGenPertTFRicci[tfricci_[-a_,-b_]]:=(ExpandPerturbation1[Perturbation[tfricci[-a_Symbol,-b_Symbol],order_.]]:=ExpandPerturbation[Perturbation[TFRicciToRicci[tfricci[-a,-b]],order]]);

4.11. Perturbations of the Kretchsmann scalar

This probably admits a simpler formula...

DefGenPertKretschmann[vbundle_,metric_,riemann_,kretschmann_[]]:=Module[{inds=GetIndicesOfVBundle[vbundle,8]},ExpandPerturbation1[Perturbation[kretschmann[],order_.]]:=ExpandPerturbation[Perturbation[metric[inds[[1]],inds[[5]]]metric[inds[[2]],inds[[6]]]metric[inds[[3]],inds[[7]]]metric[inds[[4]],inds[[8]]]riemann[-inds[[1]],-inds[[2]],-inds[[3]],-inds[[4]]]riemann[-inds[[5]],-inds[[6]],-inds[[7]],-inds[[8]]],order]]
];

5. End private and package

On[RuleDelayed::rhs]

End[]
EndPackage[]